Обращение матрицы через LDLT-разложение: итоги: nabbla1

nabbla1

Обращение матрицы через LDLT-разложение: итоги

Oct 17, 2019 20:53

Мы описали, как это делается:
раз
два
три.

Теперь подведём итоги по количеству необходимых операций. Поначалу кажется, что это удаление гланд армейским способом - мы выполняем операцию ну очень окольными путями. Что ж, сейчас сравним.

Немножко вспомним, как выполняется обычное обращение матрицы. Как ни странно, в этом методе симметричность матрицы никак не может нам помочь. По сути, мы берём "матрицу общего вида", напротив неё ставим единичную матрицу, и начинаем их преобразовывать, так что когда исходная матрица станет единичной, единичная станет обратной:

$\left(\begin{array}{cccc}a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}\end{array}\middle|\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)$

Первым делом надо верхнюю строчку поделить на a11:

$\left(\begin{array}{cccc}1&a'_{12}&a'_{13}&a'_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}\end{array}\middle|\begin{array}{cccc}1/a_{11}&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)$

По крайней мере, мы и здесь можем выполнять операции "на месте": после деления в левой части заведомо будет единица, её можно не хранить, и теперь мы здесь храним обратную величину (обзовём её b11). В то же время, другие элементы на строке в правой матрице пока нулевые.

Теперь нужно вычесть верхнюю строку, умноженную на a21, a31, a41, из 2-й, 3-й и 4-й строки соответственно. Увы, здесь это весьма хлопотно: мы перелопатим весь кусок (N-1)×(N-1):

$\left(\begin{array}{cccc}1&a'_{12}&a'_{13}&a'_{14}\\0&a_{22}-a_{21}a'_{12}&a_{23}-a_{21}a'_{13}&a_{24}-a_{21}a'_{14}\\0&a_{32}-a_{31}a'_{12}&a_{33}-a_{31}a'_{13}&a_{34}-a_{31}a'_{14}\\0&a_{42}-a_{41}a'_{12}&a_{43}-a_{41}a'_{13}&a_{44}-a_{41}a'_{14}\end{array}\middle|\begin{array}{cccc}1/a_{11}&0&0&0\\-a_{21}b_{11}&1&0&0\\-a_{31}b_{11}&0&1&0\\-a_{41}b_{11}&0&0&1\end{array}\right)$

Но принцип преобразования "на месте" вполне работает - теперь мы "забываем" первый столбец исходной матрицы, и потому с чистой совестью заполняем его столбцом преобразуемой единичной матрицы. Короче, когда мы смотрим на выкладки, то видим ровно N×N "нетривиальных" элементов (которые не равны единице или нулю).

Обозначим новые значения (чтобы не таскать за собой "старые" выкладки):
$\left(\begin{array}{cccc}1&a'_{12}&a'_{13}&a'_{14}\\0&a'_{22}&a'_{23}&a'_{24}\\0&a'_{32}&a'_{33}&a'_{34}\\0&a'_{42}&a'_{43}&a'_{44}\end{array}\middle|\begin{array}{cccc}b_{11}&0&0&0\\b_{21}&1&0&0\\b_{31}&0&1&0\\b_{41}&0&0&1\end{array}\right)$

Теперь делим вторую строку на a22:
$\left(\begin{array}{cccc}1&a'_{12}&a'_{13}&a'_{14}\\0&1&a'_{23}/a'_{22}&a'_{24}/a'_{22}\\0&a'_{32}&a'_{33}&a'_{34}\\0&a'_{42}&a'_{43}&a'_{44}\end{array}\middle|\begin{array}{cccc}b_{11}&0&0&0\\b_{21}/a'_{22}&1/a'_{22}&0&0\\b_{31}&0&1&0\\b_{41}&0&0&1\end{array}\right)$

Замечаем одну нехорошую вещь: количество требуемых операций не изменилось: мы по-прежнему должны поделить ВСЕ элементы строки - 2 в "исходной" матрице (a23 и a24) и 2 в преобразованной (b21 и b22, который изначально был единицей).

И затем мы должны вычесть вторую строку из ВСЕХ оставшихся, в том числе из первой, и сделать эту операцию по ВСЕМ столбцам: N-2 исходной матрицы и 2 из преобразованной.

Дальше можно не продолжать: становится понятно, что на каждом из N шагов количество операций одинаковое: сначала N делений (целая строка), затем N(N-1) умножений и столько же сложений.

Итого, имеем N2 делений, N2(N-1) умножений и столько же сложений.

В случае N=4, имеем 16 делений, 48 умножений и 48 сложений.
В случае N=6, имеем 36 делений, 180 умножений и 180 сложений.

А теперь просуммируем необходимые операции при трёхэтапном обращении матрицы через LDLT разложение.
Количество делений: (N-1)N/2 на первом этапе (само разложение) и ещё N на втором (обращение диагональной матрицы), итого N(N+1)/2.

Количество умножений: (N-1)N(N+4)/6 на первом этапе, (N-2)(N-1)N/6 на втором и (N-1)N(2N+5)/6 на третьем, итого (N-1)N(4N+7)/6.

Количество сложений на всех этапах, как это ни странно, ОДИНАКОВОЕ: (N-1)N(N+1)/6, итого (N-1)N(N+1)/2.

В случае N=4, имеем 10 делений, 46 умножений и 30 сложений.
В случае N=6, имеем 21 деление, 155 умножений и 105 сложений.

По асимптотике (при больших N) получится, что LDLT-метод требует вдвое меньше делений, на треть меньше умножений и вдвое меньше сложений.

Разумеется, такой метод использует практически вдвое меньше памяти: мы как хранили лишь N(N+1)/2 элементов, так и продолжаем их хранить на всём протяжении.

Так что метод определённо имеет преимущества перед вычислениями "в лоб", но не такие уж существенные.

Подозреваю, что главные преимущества - вовсе не в расходе памяти и вычислительных затратах, а скорее в устойчивости счёта и возможности "выловить блох" по ходу дела. А именно: уже на первом этапе, получая каждый следующий диагональный элемент, мы можем проверить: действительно ли он больше нуля. Если нет, значит матрица вовсе не является положительно определённой, и лучше нам "выругаться", пока не поздно!

Также мы можем практически бесплатно (затратив ещё N-1 умножений) посчитать определитель хоть исходной матрицы, хоть обратной, что бывает полезно.

Скорее всего, алгоритм можно ещё немножко "доработать напильником" - обратить диагональные элементы ещё в процессе LDLT-разложения, поскольку именно на них мы делим строку за строкой. Тем самым, сэкономим N делений. Возможно, есть и другие улучшения.

И ещё, для масштаба: обычное перемножение двух матриц N×N требует N3 умножений и N2(N-1) сложений, что даже больше, чем нам понадобилось для обращения. Почему-то человеку кажется, что нахождение обратной матрицы - операция СУЩЕСТВЕННО более сложная, но это только кажется :)

математика, программки, работа