nabbla (nabbla1) wrote,
nabbla
nabbla1

Category:

LDLT-разложение

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

У нас есть симметричная положительно определённая матрица А:



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

Чтобы её обратить, можно воспользоваться довольно забавным трёхступенчатым методом, который начинается с LDLT-разложения:



Как видно, матрица L (от слова Lower) -нижняя треугольная матрица с единицами на гл. диагонали. D-диагональная матрица, причём в случае положительно определённой матрицы A, все её компоненты должны получиться ненулевыми. Данный метод - разновидность разложения Холецкого, но только в "обычном" разложении Холецкого мы записываем A = LLT, и приходится посчитать аж N (для матрицы N×N) квадратных корней. А если корни считать не хочется, то LDLT-разложение - то, что доктор прописал.

Две эти матрицы вполне помещаются в памяти, выделенной под матрицу A:



и это разложение может быть выполнено "на месте", не требуя временного хранилища.


Чтобы понять, как это сделать, помножим 3 матрицы и посмотрим, как их коэффициенты соотносятся с коэффициентами исходной матрицы:



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


Можно идти в разном порядке, но всегда это будет от левого верхнего угла к правому нижнему. Наверное, наиболее очевидный способ - идти по столбцам, слева направо, а в каждом столбце - сверху вниз.

Начинаем с первого столбца:




Каждый раз мы заменяем старое значение новым, по сути получив такую промежуточную матрицу:



Теперь обрабатываем второй столбец:




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


так что, посчитав его один раз, можно сэкономить 2 умножения.

Обрабатываем третий столбец:



Количество работы в каждом компоненте растёт, зато их количество падает. Здесь мы снова можем чуточку сэкономить, введя переменные


экономия снова составит 2 умножения.

Наконец, на 4-м столбце у нас все один, диагональный элемент:


На этом преобразование завершено.

Посчитаем необходимое количество операций для матрицы N×N.
Количество делений составляет

квадратичная зависимость, в нашем случае N=4 имеем 6 делений, для N=6 получается 15 делений.

Количество умножений составляет (если не вводить временных переменных)

здесь зависимость уже кубическая - ничего не поделаешь, но у нас, по счастью, размерности не очень большие. Для N=4 получаем 20 умножений, для N=6: 70 умножений.

Если ввести временные переменные, или поменять местами 2 цикла (вычитать из всех элементов на текущем столбце по одному множителю, вместо вычисления каждого из них "в один присест"), то получим значение чуть лучше:


Для N=4, получаем 16 умножений (да, как и обещали - сэкономили по 2 в двух строках), а для N=6: 50 умножений. Неплохо :)

И наконец, можно подсчитать количество сложений/вычитаний (мы разницы между ними не делаем):

для N=4 получаем 10 сложений, для N=6: 35 сложений.

Метод славится устойчивостью: нам не нужно предварительно менять местами строки матрицы, стараясь расположить самые большие элементы на главной диагонали. Считается, что здесь pivoting особенной пользы не приносит, и это очень большое подспорье.


Дальше расскажем, как обратить нижнюю треугольную матрицу с единицами на главной диагонали "на месте".
Tags: tex, математика, работа
Subscribe

  • Я создал монстра!

    Вот нормальная счастливая пара разъёмов ОНЦ-БС-1-10/14-Р12-2-В и ОНЦ-БС-1-10/14-В1-2-В: У розетки кроме основного выступа, отмечающего "верх",…

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

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

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

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

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 20 comments

  • Я создал монстра!

    Вот нормальная счастливая пара разъёмов ОНЦ-БС-1-10/14-Р12-2-В и ОНЦ-БС-1-10/14-В1-2-В: У розетки кроме основного выступа, отмечающего "верх",…

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

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

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

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