Предобусловливание СЛАУ

[i_m_mortal]

Здравствуйте.
Помогите, пожалуйста, разобраться со следующей проблемой.
Решаю систему уравнений $Ax = b$, где А - разреженная матрица. Для решения применяю итерационный метод Гаусс-Зейделя:
$A = L + D + R$,
где L - нижнетреугольная, D - диагональная, R - верхнетреугольная матрицы.
Далее
$(L+D)x_{n+1} + Rx_n = b$,
$x_{n+1} = D^{-1}(b - Rx_n - Lx_{n+1})$.

Хочу использовать предобусловливание исходной системы ($Ax = b$) для ускорения счета:
$P^{-1}Ax = P^{-1}b$,
где P - предобусловливатель.
Применяю метод неполной факторизации LU
$P = ILU(0)$,
получая
$(LU)^{-1}Ax = (LU)^{-1}b$.
Но дальше я зашел в тупик. Совсем не представляю как от последнего уравнения перейти к итерационному методу.
Отсюда, у меня возникло несколько вопросов:
Как перейти от последнего уравнения к следующему:
$x_{n+1} = f(x_n, x_{n+1})$?
И действительно ли необходимо находить обратные матрицы к L и U, ведь тогда будет проще решить систему прямым методом LU-факторизации?

Спасибо огромное за внимание!
С уважением, Артем.

[Kallikanzarid]

1) Заменяете исходную задачу на предобусловленную - и вперед!
2) Обращать $L$ и $U$ вам не придется, у вас же написано: $x_{n+1} = D^{-1}(b - Rx_n - Lx_{n+1})$. Чтобы посчитать $(Lx_{n+1})^{i+1}$, вам достаточно знать $(Lx_{n+1})^1, \ldots, (Lx_{n+1})^i$. Тупо делаете цикл по массиву, в котором храните компоненты $x$ - простейший итерационный метод в этой галактике!
3) Прямые методы очень медленные, и нет гарантии, что при LU-факторизации у вас получатся разреженные матрицы.

[i_m_mortal]

Здравсвуйте.
Спасибо за ответ!

Простите, Kallikanzarid, у меня вышло повторение обозначения L. Я имел в виду следующее:
$Ax = b$
$A = E + D + F$,
где E - нижнетреугольная, D - диагональная, F - верхнетреугольная матрицы.
Тогда
$x_{n+1} = D^{-1}(b - Fx_n - Ex_{n+1})$.
А после предобусловливания получаем:
$(LU)^{-1}Ax = (LU)^{-1}b$.
Т.е. простой переход $A = E + D + F$ не проясняет ситуацию:
$(LU)^{-1}(E + D + F)x = (LU)^{-1}b$.

Из литературы я вычитал, что предобусловливание бывает левым, правым и разнесенным. Вот что выходит при попытке применить каждый из них:
1. Левое предобусловливание
$(LU)^{-1}Ax = (LU)^{-1}b$
Приведем уравнение к виду
$x_{n+1} = Gx_n + f$,
где $G = M^{-1}N = I - M^{-1}A$.
$(I-G)x = f$ или $M^{-1}Ax = M^{-1}b$.
Если предположить, что $M^{-1} = (LU)^{-1}$, то обратной подстановкой прийдем к
$x_{n+1} = (I - (LU)^{-1}A)x_n + (LU)^{-1}b$.

2. Правое предобусловливание
$A(LU)^{-1}y = b$,
$x = (LU)^{-1}y$.

3. Разнесенное предобусловливание
$L^{-1}AU^{-1}y = L^{-1}b$,
$x = U^{-1}y$.

Т.е. во всех случаях получается необходимо находить обратные матрицы к L и U? Ладно, когда операция обращения применяется к матрице D (из первого абзаца), в уравнении просто будет присутствовать $1/D_{ii}$. А ведь матрицы L и U непросто обратить.

Спасибо огромное за внимание!
С уважением, Артем.

[Kallikanzarid]

Я занимался такими задачами, так что не волнуйтесь, я все правильно понял :) Неполная LU-факторизация ищет приближенное разложение исходной матрицы, обычно на ее вычисление уходит $O(n)$ операций, и полученные матрицы - разреженные. Но если бы вы использовали прямой метод, то на вычисление точных $L$ и $U$ у вас бы ушло $O(n^3)$ операций, и из разреженность не была бы гарантирована.

Непосредственно обращать матрицы $L$ и $U$ вам не нужно, просто на каждой итерации вычисляйте $(LU)^{-1}x = U^{-1}L^{-1}x$, в силу специального вида матриц на это уйдет всего лишь $O(n)$ операций.

[i_m_mortal]

Добрый день, Kallikanzarid.
Огромное спасибо за ответ!

Т.е. вы предлагаете перейти от $(LU)^{-1}Ax=(LU)^{-1}b$ к $Cx=f$, заменив и вычислив предварительно $C = (LU)^{-1}A = L^{-1}U^{-1}A$ и $f = (LU)^{-1}b = L^{-1}U^{-1}b$? Я попробую сравнить быстродействие такого подхода и обычным непредобусловленным Гаусс-Зейделем (GS).

А что вы скажете про это?
Нашел в теме о применении методов GS и ILU в качестве сглаживателей в многосеточном методе: $M(x^{k+1} - x^k) = b - Ax^k$, где M - предобусловливатель.
GS: $x_i^{k+1} = (b_i - \sum_{j<i}a_{ij}x_j^{k+1} - \sum_{j>i}a_{ij}x_j^{k})/a_{ii}$ или
$(D_A + L_A)D_A^{-1}(D_A + U_A)(x^{k+1}-x^k) = b - Ax^k$.
ILU: $x_i^{k+1} = (b_i - \sum_{j<i}a_{ij}x_j^{k+1} - \sum_{j>i}a_{ij}x_j^{k})/d_{ii}$ или
$(D + L_A)D^{-1}(D + U_A)(x^{k+1}-x^k) = b - Ax^k$,
где $d_{ii} = a_{ii} - \sum_{j<i}\frac{a_{ij}a{ji}}{d{jj}}$.
Получается, что здесь не надо явно находить $L^{-1}U^{-1}$, но является ли это реализацией предобусловленного GS с помощью ILU?

Спасибо огромное за внимание!
С уважением, Артем.

[Kallikanzarid]

i_m_mortal в сообщении #445239 писал(а):

Т.е. вы предлагаете перейти от $(LU)^{-1}Ax=(LU)^{-1}b$ к $Cx=f$, заменив и вычислив предварительно $C = (LU)^{-1}A = L^{-1}U^{-1}A$ и $f = (LU)^{-1}b = L^{-1}U^{-1}b$? Я попробую сравнить быстродействие такого подхода и обычным непредобусловленным Гаусс-Зейделем (GS).

Не совсем: обычно не считают матрицу $C = (LU)^{-1}A$, а вместо этого на каждой итерации считают непосредственно $U^{-1}(L^{-1}(Ax))$. ИМХО в интернете должны быть примеры реализации или псевдокода предобусловленных методов, посмотрите их для примера.