LDLT-разложение: особая индийская магия!
Мы изложили "стандартный" подход к обращению симметричной положительно определённой матрицы через LDLT-разложение с последующим обращением матриц L и D, и получением итогового результата как результат умножения 3 матриц.
В итоге количество умножений получилось порядка 2/3N3, что конечно получше, чем N3 в алгоритме обращения матрицы "общего вида", но не настолько лучше, как хотелось бы. Казалось, что когда у тебя вдвое меньше "уникальных" элементов как на входе, так и на выходе, то и количество операций должно упасть почти вдвое. Но всё-таки "окольные пути" взяли своё.
Но не так давно (в 2011 году с правками в 2013) два индийских математика опубликовали статью, в которой заявили, что их метод позволяет обратить такую матрицу всего за 1/2N3 умножений, что на 1/6 лучше существующих методов. Что ещё интереснее, они проверяли работу "с фиксированной точкой" и получили, что такой метод даёт чуть меньше ошибок вычислений (так я на эту статью и вышел - искал подсказки, как это всё получше выполнить без FPU). Я, конечно, не удивлюсь, если всё это уже было у Гаусса, но пока не слышал о таком :)
Сейчас попытаюсь объяснить, что же там происходит...
Первый этап всё тот же самый, сделать LDLT-разложение:
Дальше, мы хотим решить уравнение
относительно неизвестной матрицы X. (E - это единичная матрица)
Подставляем наше разложение:
Умножаем обе части слева на L-1:
И ещё умножим обе части слева на D-1:
Это и есть то выражение, из которого мы хотим найти матрицу X, то бишь, обратную матрицу!
Кажется, что в этом нет никакого смысла: нам по-прежнему надо обращать матрицу L, а когда мы её обратим - то почему бы не посчитать X напрямую!?
В этом и состоит великая индийская хитрость. Изобразим это уравнение поэлементно:
%5Cleft(%5Cbegin{array}{cccc}x_{11}&x_{12}&x_{13}&x_{14}%5C%5Cx_{21}&x_{22}&x_{23}&x_{24}%5C%5Cx_{31}&x_{32}&x_{33}&x_{43}%5C%5Cx_{41}&x_{42}&x_{43}&x_{44}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{cccc}1/d_1&0&0&0%5C%5C0&1/d_2&0&0%5C%5C0&0&1/d_3&0%5C%5C0&0&0&1/d_4%5Cend{array}%5Cright)%5Cleft(%5Cbegin{array}{cccc}1&0&0&0%5C%5C(l^{-1})_{21}&1&0&0%5C%5C(l^{-1})_{31}&(l^{-1})_{32}&1&0%5C%5C(l^{-1})_{41}&(l^{-1})_{42}&(l^{-1})_{43}&1%5Cend{array}%5Cright))
Помножим матрицы в правой части, причём на нижнем треугольнике поставим прочерки - как оказывается, он нам не нужен!
%5Cleft(%5Cbegin{array}{cccc}x_{11}&x_{12}&x_{13}&x_{14}%5C%5Cx_{21}&x_{22}&x_{23}&x_{24}%5C%5Cx_{31}&x_{32}&x_{33}&x_{43}%5C%5Cx_{41}&x_{42}&x_{43}&x_{44}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{cccc}1/d_1&0&0&0%5C%5C-&1/d_2&0&0%5C%5C-&-&1/d_3&0%5C%5C-&-&-&1/d_4%5Cend{array}%5Cright))
Начинаем решать эту штуку, начиная с правого столбца:
%5Cleft(%5Cbegin{array}{c}x_{14}%5C%5Cx_{24}%5C%5Cx_{34}%5C%5Cx_{44}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}0%5C%5C0%5C%5C0%5C%5C1/d_4%5Cend{array}%5Cright))
Решаем это прямой подстановкой, снизу вверх:
Пока всё было понятно. Но теперь берёмся за следующий столбец, в котором одного, самого нижнего значения мы не знаем:
%5Cleft(%5Cbegin{array}{c}x_{13}%5C%5Cx_{23}%5C%5Cx_{33}%5C%5Cx_{43}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}0%5C%5C0%5C%5C1/d_3%5C%5C-%5Cend{array}%5Cright))
Как будто бы значение x43 найти не удастся, а без него "посыпется" и всё остальное... Да вот только мы же знаем, что наша матрица симметричная, поэтому x43 = x34, а его мы уже нашли в прошлый раз! А зная его, можно выразить и все остальные:
"Для закрепления пройденного", выпишем следующий по счёту столбец, второй:
%5Cleft(%5Cbegin{array}{c}x_{12}%5C%5Cx_{22}%5C%5Cx_{32}%5C%5Cx_{42}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}0%5C%5C1/d_2%5C%5C-%5C%5C-%5Cend{array}%5Cright))
Теперь мы "не знаем" уже двух переменных: x42 и x32. А в действительности нашли уже обе: x42=x24 и x32=x23. Выразим оставшиеся две:
И наконец, берёмся за первый столбец:
%5Cleft(%5Cbegin{array}{c}x_{11}%5C%5Cx_{21}%5C%5Cx_{31}%5C%5Cx_{41}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{c}1/d_1%5C%5C-%5C%5C-%5C%5C-%5Cend{array}%5Cright))
Три нижних значения нам уже известны, это x12, x13 и x14, найденные ранее. Осталось выразить самую последнюю неизвестную:
На этом работа завершена.
Итак, вместо двух подстановок (как здесь), удалось обойтись только одной, благодаря тому, что мы не испугались наполовину не заполненной матрицы в правой части, зная, что в своё время сможем заполнить эти пробелы!
В итоге, количество умножений на этом шаге составляет порядка N3/3, что прибавляясь к нашим исходным (на этапе LDLT - разложения) N3/6, и даёт заявленное значение N3/2.
Шикарный метод, как по мне!
В итоге количество умножений получилось порядка 2/3N3, что конечно получше, чем N3 в алгоритме обращения матрицы "общего вида", но не настолько лучше, как хотелось бы. Казалось, что когда у тебя вдвое меньше "уникальных" элементов как на входе, так и на выходе, то и количество операций должно упасть почти вдвое. Но всё-таки "окольные пути" взяли своё.
Но не так давно (в 2011 году с правками в 2013) два индийских математика опубликовали статью, в которой заявили, что их метод позволяет обратить такую матрицу всего за 1/2N3 умножений, что на 1/6 лучше существующих методов. Что ещё интереснее, они проверяли работу "с фиксированной точкой" и получили, что такой метод даёт чуть меньше ошибок вычислений (так я на эту статью и вышел - искал подсказки, как это всё получше выполнить без FPU). Я, конечно, не удивлюсь, если всё это уже было у Гаусса, но пока не слышал о таком :)
Сейчас попытаюсь объяснить, что же там происходит...
Первый этап всё тот же самый, сделать LDLT-разложение:
Дальше, мы хотим решить уравнение
относительно неизвестной матрицы X. (E - это единичная матрица)
Подставляем наше разложение:
Умножаем обе части слева на L-1:
И ещё умножим обе части слева на D-1:
Это и есть то выражение, из которого мы хотим найти матрицу X, то бишь, обратную матрицу!
Кажется, что в этом нет никакого смысла: нам по-прежнему надо обращать матрицу L, а когда мы её обратим - то почему бы не посчитать X напрямую!?
В этом и состоит великая индийская хитрость. Изобразим это уравнение поэлементно:
Помножим матрицы в правой части, причём на нижнем треугольнике поставим прочерки - как оказывается, он нам не нужен!
Начинаем решать эту штуку, начиная с правого столбца:
Решаем это прямой подстановкой, снизу вверх:
Пока всё было понятно. Но теперь берёмся за следующий столбец, в котором одного, самого нижнего значения мы не знаем:
Как будто бы значение x43 найти не удастся, а без него "посыпется" и всё остальное... Да вот только мы же знаем, что наша матрица симметричная, поэтому x43 = x34, а его мы уже нашли в прошлый раз! А зная его, можно выразить и все остальные:
"Для закрепления пройденного", выпишем следующий по счёту столбец, второй:
Теперь мы "не знаем" уже двух переменных: x42 и x32. А в действительности нашли уже обе: x42=x24 и x32=x23. Выразим оставшиеся две:
И наконец, берёмся за первый столбец:
Три нижних значения нам уже известны, это x12, x13 и x14, найденные ранее. Осталось выразить самую последнюю неизвестную:
На этом работа завершена.
Итак, вместо двух подстановок (как здесь), удалось обойтись только одной, благодаря тому, что мы не испугались наполовину не заполненной матрицы в правой части, зная, что в своё время сможем заполнить эти пробелы!
В итоге, количество умножений на этом шаге составляет порядка N3/3, что прибавляясь к нашим исходным (на этапе LDLT - разложения) N3/6, и даёт заявленное значение N3/2.
Шикарный метод, как по мне!