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

nabbla1

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

Oct 16, 2019 18:08

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

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

$A=\left(\begin{array}{cccc}a_{11}&a_{21}&a_{31}&a_{41}\\a_{21}&a_{22}&a_{32}&a_{42}\\a_{31}&a_{32}&a_{33}&a_{43}\\a_{41}&a_{42}&a_{43}&a_{44}\end{array}\right)=\left(\begin{array}{cccc}a_{11}&-&-&-\\a_{21}&a_{22}&-&-\\a_{31}&a_{32}&a_{33}&-\\a_{41}&a_{42}&a_{43}&a_{44}\end{array}\right)$

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

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

$A=LDL^T=\left(\begin{array}{cccc}1&0&0&0\\l_{21}&1&0&0\\l_{31}&l_{32}&1&0\\l_{41}&l_{42}&l_{43}&1\end{array}\right)\left(\begin{array}{cccc}d_1&0&0&0\\0&d_2&0&0\\0&0&d_3&0\\0&0&0&d_4\end{array}\right)\left(\begin{array}{cccc}1&l_{21}&l_{31}&l_{41}\\0&1&l_{32}&l_{42}\\0&0&1&l_{43}\\0&0&0&1\end{array}\right)$

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

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

$A^*=\left(\begin{array}{cccc}d_1&-&-&-\\l_{21}&d_2&-&-\\l_{31}&l_{32}&d_3&-\\l_{41}&l_{42}&l_{43}&d_4\end{array}\right)$

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

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

$LDL^T=\left(\begin{array}{cccc}d_1&d_1l_{21}&d_1l_{31}&d_1l_{41}\\d_1l_{21}&d_1l_{21}^2+d_2&d_1l_{21}l_{31}+d_2l_{32}&d_1l_{21}l_{41}+d_2l_{42}\\d_1l_{31}&d_1l_{21}l_{31}+d_2l_{32}&d_1l_{31}^2+d_2l_{32}^2+d_3&d_1l_{31}l_{41}+d_2l_{32}l_{42}+d_3l_{43}\\d_1l_{41}&d_1l_{21}l_{41}+d_2l_{42}&d_1l_{31}l_{41}+d_2l_{32}l_{42}+d_3l_{43}&d_1l_{41}^2+d_2l_{42}^2+d_3l_{43}^2+d_4\end{array}\right)$

разумеется, и здесь мы получили симметричную матрицу, так что уберём ненужный верхний треугольник:
$A=LDL^T=\left(\begin{array}{cccc}d_1&-&-&-\\d_1l_{21}&d_1l_{21}^2+d_2&-&-\\d_1l_{31}&d_1l_{21}l_{31}+d_2l_{32}&d_1l_{31}^2+d_2l_{32}^2+d_3&-\\d_1l_{41}&d_1l_{21}l_{41}+d_2l_{42}&d_1l_{31}l_{41}+d_2l_{32}l_{42}+d_3l_{43}&d_1l_{41}^2+d_2l_{42}^2+d_3l_{43}^2+d_4\end{array}\right)$

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

Начинаем с первого столбца:
$d_1=a_{11},$

$l_{21}=a_{21}/d_1,$

$l_{31}=a_{31}/d_1,$

$l_{41}=a_{41}/d_1,$

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

$\left(\begin{array}{cccc}d_1&-&-&-\\l_{21}&a_{22}&-&-\\l_{31}&a_{32}&a_{33}&-\\l_{41}&a_{42}&a_{43}&a_{44}\end{array}\right)$

Теперь обрабатываем второй столбец:
$d_2=a_{22}-d_1l_{21}^2,$

$l_{32}=\frac{a_{32}-d_1l_{21}l_{31}}{d_2},$

$l_{42}=\frac{a_{42}-d_1l_{21}l_{41}}{d_2}.$

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

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

Обрабатываем третий столбец:
$d_3=a_{33}-d_1l_{31}^2-d_2l_{32}^2,$

$l_{43}=\frac{a_{43}-d_1l_{31}l_{41}-d_2l_{32}l_{42}}{d_3}$

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

$t_2=d_2l_{32},$

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

Наконец, на 4-м столбце у нас все один, диагональный элемент:
$d_4=a_{44}-d_1l_{41}^2-d_2l_{42}^2-d_3l_{43}^2$

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

Посчитаем необходимое количество операций для матрицы N×N.
Количество делений составляет
$(N-1)+(N-2)+...+0=\frac{(N-1)N}{2}$

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

Количество умножений составляет (если не вводить временных переменных)
$2(0\cdot{N}+1\cdot(N-1)+2\cdot(N-2)+...+N\cdot{0})=\frac{(N-1)N(N+1)}{3}$

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

Если ввести временные переменные, или поменять местами 2 цикла (вычитать из всех элементов на текущем столбце по одному множителю, вместо вычисления каждого из них "в один присест"), то получим значение чуть лучше:
$1\cdot(1+N-1)+2\cdot(1+N-2)+...+(N-1)\cdot{2})=\frac{(N-1)N(N+4)}{6}$

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

И наконец, можно подсчитать количество сложений/вычитаний (мы разницы между ними не делаем):
$1\cdot(N-1)+2\cdot(N-2)+...+(N-1)\cdot{1}=\frac{(N-1)N(N+1)}{6}$

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

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

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

tex, математика, работа