nabbla (nabbla1) wrote,
nabbla
nabbla1

Category:

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

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

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


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



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



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

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



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

Обозначим новые значения (чтобы не таскать за собой "старые" выкладки):


Теперь делим вторую строку на a22:


Замечаем одну нехорошую вещь: количество требуемых операций не изменилось: мы по-прежнему должны поделить ВСЕ элементы строки - 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) сложений, что даже больше, чем нам понадобилось для обращения. Почему-то человеку кажется, что нахождение обратной матрицы - операция СУЩЕСТВЕННО более сложная, но это только кажется :)
Tags: математика, программки, работа
Subscribe

  • Нахождение двух самых отдалённых точек

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

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 7 comments