Обращение нижней унитреугольной матрицы
Нам хочется обратить симметричную положительно определённую матрицу А.
Первым делом мы сделали LDLT-разложение:
Здесь D-диагональная матрица, а L-нижняя треугольная матрица с единицами на главной диагонали, что иногда называют термином унитреугольная матрица.
Выразим обратную матрицу:
^TD^{-1}L^{-1})
Обратить диагональную матрицу D легко, мы просто заменяем каждое число обратной величиной:
)
С нижней унитреугольной матрицей L самую чуточку посложнее.
Попробуем обратить её "вручную" по обычным правилам: справа записываем единичную матрицу и преобразуем их одновременнно:
)
Верхняя строка уже в полном порядке. Далее вычитаем верхнюю строку из всех остальных, помножая на первые элементы строк. Как видно, здесь деление уже не требуется! Получаем:
)
Снова наблюдаем замечательную "синхроничность": когда появляются ненулевые значения в правой части (в преобразуемой единичной матрицы), в левой части (в исходной матрице) появляются нули, т.е эти значения "выходят из игры". Так что мы снова можем орудовать "на месте". В данном случае мы попросту сменили знак всех элементов первого столбца. Переобзовём их как m21, m31 и т.д. - элементы обратной матрицы, пока ещё не окончательные:
)
Теперь и вторая строка в полном порядке. Вычитаем её из оставшихся, помножив на l32 и l42, соответственно:
)
Всё идёт по плану: ещё один столбец исходной матрицы "освободился", на его место записываем значения из "правой" матрицы (а также перезаписались значения m21 и m31:
)
Остаётся преобразовать последнюю строку, а именно, вычесть из неё третью, помноженную на l43:
)
Вот, собственно, и всё. Как и обещалось, такая операция не требует дополнительной памяти для работы, всё происходит на том же самом месте. Но заметим, что мы можем в качестве аргументов для алгоритма (будь то LDLT-разложение, будь обращение матрицы L) задать указатель на исходную матрицу и указатель на результирующую (с уже выделенной памятью), с комментарием: они могут указывать на одно и то же место, и результат по-прежнему будет корректным!
Оценим вычислительную сложность этого метода. На первом шаге (на первом столбце) мы просто поменяли знак для N-1 чисел, будем считать это N-1 "сложениями".
На втором шаге нам понадобилось N-2 умножений и столько же сложений, чтобы обновить элементы первого столбца, и ещё N-2 "сложений" для второго столбца.
На третьем шаге: 2(N-3) умножений и столько же сложений, чтобы обновить элементы 1-го и 2-го столбца, и ещё N-3 "сложений" для третьего.
И всего итераций N-1, в нашем случае как раз три.
Полное количество умножений составляет
+1%5Ccdot(N-2)+2%5Ccdot(N-3)+...+(N-1)%5Ccdot{1}=%5Cfrac{(N-2)(N-1)N}{6})
Да, немного странно видеть здесь кубическую зависимость - кажется, что вычислений очень мало, и для N=4 так оно и выходит - всего 4 умножения :) Для N=6 чуть похуже: 20 умножений!
Полное количество сложений (в том числе обращения знака) составляет
+(N-2)+...+1+%5Cfrac{(N-2)(N-1)N}{6}=%5Cfrac{(N-1)N(N+1)}{6})
Для N=4 получится 10 сложений, а для N=6: 35 сложений.
Остаётся последний этап: помножить обращённую матрицу L-1 на D-1 и на себя же транспонированную, и тоже очень желательно "на месте", раз уж начали!
Первым делом мы сделали LDLT-разложение:
Здесь D-диагональная матрица, а L-нижняя треугольная матрица с единицами на главной диагонали, что иногда называют термином унитреугольная матрица.
Выразим обратную матрицу:
Обратить диагональную матрицу D легко, мы просто заменяем каждое число обратной величиной:
С нижней унитреугольной матрицей L самую чуточку посложнее.
Попробуем обратить её "вручную" по обычным правилам: справа записываем единичную матрицу и преобразуем их одновременнно:
Верхняя строка уже в полном порядке. Далее вычитаем верхнюю строку из всех остальных, помножая на первые элементы строк. Как видно, здесь деление уже не требуется! Получаем:
Снова наблюдаем замечательную "синхроничность": когда появляются ненулевые значения в правой части (в преобразуемой единичной матрицы), в левой части (в исходной матрице) появляются нули, т.е эти значения "выходят из игры". Так что мы снова можем орудовать "на месте". В данном случае мы попросту сменили знак всех элементов первого столбца. Переобзовём их как m21, m31 и т.д. - элементы обратной матрицы, пока ещё не окончательные:
Теперь и вторая строка в полном порядке. Вычитаем её из оставшихся, помножив на l32 и l42, соответственно:
Всё идёт по плану: ещё один столбец исходной матрицы "освободился", на его место записываем значения из "правой" матрицы (а также перезаписались значения m21 и m31:
Остаётся преобразовать последнюю строку, а именно, вычесть из неё третью, помноженную на l43:
Вот, собственно, и всё. Как и обещалось, такая операция не требует дополнительной памяти для работы, всё происходит на том же самом месте. Но заметим, что мы можем в качестве аргументов для алгоритма (будь то LDLT-разложение, будь обращение матрицы L) задать указатель на исходную матрицу и указатель на результирующую (с уже выделенной памятью), с комментарием: они могут указывать на одно и то же место, и результат по-прежнему будет корректным!
Оценим вычислительную сложность этого метода. На первом шаге (на первом столбце) мы просто поменяли знак для N-1 чисел, будем считать это N-1 "сложениями".
На втором шаге нам понадобилось N-2 умножений и столько же сложений, чтобы обновить элементы первого столбца, и ещё N-2 "сложений" для второго столбца.
На третьем шаге: 2(N-3) умножений и столько же сложений, чтобы обновить элементы 1-го и 2-го столбца, и ещё N-3 "сложений" для третьего.
И всего итераций N-1, в нашем случае как раз три.
Полное количество умножений составляет
Да, немного странно видеть здесь кубическую зависимость - кажется, что вычислений очень мало, и для N=4 так оно и выходит - всего 4 умножения :) Для N=6 чуть похуже: 20 умножений!
Полное количество сложений (в том числе обращения знака) составляет
Для N=4 получится 10 сложений, а для N=6: 35 сложений.
Остаётся последний этап: помножить обращённую матрицу L-1 на D-1 и на себя же транспонированную, и тоже очень желательно "на месте", раз уж начали!