LLT-Факторизация

[yog]

Задача:
C = A + B
A, B, C - матрицы симметричные, положительно определенные
матрица A представляется в виде A=L_AL_A^T
матрица B представляется в любом удобном для решения виде, предположительно также B=L_B L_B^T
Известно: A, B и их разложения
Вопрос: Каким образом можно факторизовать матрицу С в виде С=L_C L_C^T не проводя повторную процедуру факторизации.
Или формально: factor( A+B ) = function( factor( A ), factor( B ) ).
Найти функцию или метод function(). Достаточно в виде алгоритма, можно приближенного.

P.S. Подобная задача возникает при решении численными методами некоторых задач в механике твёрдого тела.

[незваный гость]

Первая идея -- увы, Вы ищете то, чего не существует. Пример, который я приведу, возможно тривиален, и не отражает всей сложности задачи. Возьмем одномерную систему, $A = a^2$, $B = b^2$. Тогда разложения $L_A =a$, $L_B = b$ очевидны. Но вот для $L_C = \sqrt{a^2+b^2}$ надо снова извлекать корень. Я понимаю, что это аналогия, и что она может быть не верна, и тем не менее...

Вы можете попытаться сконструировать итеративный алгоритм для специальных случаев, например, когда $||A|| \gg||B||$. Это может пройти... но и это не очевидно.

[Yuri Gendelman]

незваный гость писал(а):

Первая идея -- увы, Вы ищете то, чего не существует.

У меня тоже есть подозрение, что если даже такой алгоритм существует, то он не проще, чем прямое разложение матрицы C. В конце концов алгоритм Холецкого не такой уж сложный.

Единственный опубликованный алгоритм такого рода, который я знаю: "Корректировка обратной матрицы после изменения одного элемента в прямой матрице". То есть, если известны обратная матрица $A = M^{-1}$ и изменение одного элемента исходной матрицы $d = M^{new}_{i,j} - M_{i,j}$, то вычислить $B = (M^{new})^{-1}$ можно быстрее, чем непосредственно обращать новую матрицу. Недостаток метода - потеря точности.
Ссылки, если кому интересно (Algol-60 !):
- Herndon J.R., CACM, 1961, 4, algorithm 51
- М.И.Агеев, В.П.Алик, Ю.И.Марков, Библиотека алгоритмов 51б-100б, М., "Советское Радио", 1976 (алгоритм 51б).
Обобщение алгоритма - изменение строки или столбца, а еще "общЕе" - прибавление матрицы ранга 1, описывается кажется в книге Фаддеевых.

[yog]

Yuri Gendelman писал(а):

незваный гость писал(а):

Первая идея -- увы, Вы ищете то, чего не существует.

У меня тоже есть подозрение, что если даже такой алгоритм существует, то он не проще, чем прямое разложение матрицы C. В конце концов алгоритм Холецкого не такой уж сложный.

В том то и дело, что Холецкого применять немного накладно.
Матрица размерностью 10^6. Правда сильно разреженная. А вторая матрица B и вовсе почти полностью из нулей состоит. Мне кажется, необходимо копать в сторону приближенных вычислений.

Большое спасибо за подсказки и ссылки!

[Yuri Gendelman]

yog писал(а):

Матрица размерностью 10^6. Правда сильно разреженная. А вторая матрица B и вовсе почти полностью из нулей состоит. Мне кажется, необходимо копать в сторону .

Если приближенными Вы называете в итерационные методы, то да. Вообще большие разреженные матрицы типа используемых в МКЭ - это особый раздел вычислительной математики. Стандартные методы для них не очень эффективны. В программах, реализующих МКЭ, матрицы иногда стараются привести к ленточным. Для этого применяются алгоритмы "уменьшения ширины" перестановкой строк и столбцов. В более сложных алгоритмах используется упорядочение дающее минимальное число ненулевых элементов треугольника Холецкого.

Ссылки:
1) А.Джордж, Дж.Лю. Численное решение больших разреженных систем уравнений. - М.: Мир. 1984
2) Fundamentals of matrix computations, by Watkins D. - Wiley. 2002

[yog]

По вышеуказанным ссылкам наткнулся на одну полезную формулу Shermann - Morrison:
(A + $\vec u \otimes \vec v$)^{-1} = A^{-1} - \frac{(A^{-1}\vec u)\otimes(\vec v A^{-1})}{1+\lambda}
где \lambda = \vec v A^{-1}\vec u
в моей задаче матрицу B можно представить в простом виде:
B = $\vec b \otimes \vec b$
где вектор $\vec b$ размерностью 10^6 имеет порядка 10^3 ненулевых элементов
вводим новый вектор $\vec k = A^{-1}\vec b$
и формула приобретает простой и удобный для расчётов (в сравнении с алгоритмом Холецкого) вид:
y = A_{new}^{-1}\vec f = A^{-1}\vec f - \frac{\vec k(\vec k, \vec f)}{1+(\vec b, \vec f)}
y - решение системы A$\vec y = \vec f$

Возможно, это решение не самое оптимальное...

[Yuri Gendelman]

yog писал(а):

...наткнулся на одну полезную формулу Shermann - Morrison:
(A + $\vec u \otimes \vec v$)^{-1} = A^{-1} - \frac{(A^{-1}\vec u)\otimes(\vec v A^{-1})}{1+\lambda}

Да, этот алгоритм описан в книге Фаддеевых. Работает он быстро, вряд ли Вы найдете что-нибудь лучше. Недостаток метода - потеря точности, если алгоритм многократно применяется к одной и той же исходной матрице. Тогда желательно время от времени пересчитывать обратную матрицу стандартным алгоритмом.