Обращение симм. пол. опр. матрицы - последний этап
У нас была матрица A - симметричная N×N, положительно определённая. В памяти хранится только N(N+1)/2 её коэффициентов.
Чтобы обратить её, мы сначала сделали LDLT-разложение "на месте" (хотя можно и на новое место, если исходная матрица ещё пригодится).
Затем обратили матрицы L и D, также "на месте".
Кстати, матрицы L-1 с D-1 и сами по себе представляют интерес: если A-1, которую мы ищем - это ковариационная матрица, то L-1 совместно с D-1 позволяет генерировать вектор случайных величин, удовлетворяющий этой ковариационной матрице. Хотя LLT-разложение в этой задаче, наверное, чуть лучше подходит...
Наконец, нужно перемножить матрицы, чтобы найти обратную матрицу:
^TD^{-1}L^{-1})
Посмотрим, как можно произвести это умножение "за один этап", и притом по-прежнему "на месте".
Выкинем "-1" (обращение) - неудобно его за собой таскать. Как бы рассмотрим отдельную задачу: как посчитать выражение
для нижней унитреугольной матрицы L и диагональной матрицы D.
Выпишем такое умножение в явном виде для матриц 4×4. (почему мы выбрали 4×4 - потому что это минимальный размер матрицы, где начинает проглядываться весь алгоритм целиком, все 3 вложенных цикла. До этого кажется, что циклов всего 2). Слева - мы уже помножили диагональную матрицу на нижнюю треугольную, справа - итоговый результат:
%5Cleft(%5Cbegin{array}{cccc}d_1&0&0&0%5C%5Cd_2l_{21}&d_2&0&0%5C%5Cd_3l_{31}&d_3l_{32}&d_3&0%5C%5Cd_4l_{41}&d_4l_{42}&d_4l_{43}&d_4%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{cccc}d_1+d_2l_{21}^2+d_3l_{31}^2+d_4l_{41}^2&d_2l_{21}+d_3l_{31}l_{32}+d_4l_{41}l_{42}&d_3l_{31}+d_4l_{41}l_{43}&d_4l_{41}%5C%5Cd_2l_{21}+d_3l_{31}l_{32}+d_4l_{41}l_{42}&d_2+d_3l_{32}^2+d_4l_{42}^2&d_3l_{32}+d_4l_{42}l_{43}&d_4l_{42}%5C%5Cd_3l_{31}+d_4l_{41}l_{43}&d_3l_{32}+d_4l_{42}l_{43}&d_3+d_4l_{43}^2&d_4l_{43}%5C%5Cd_4l_{41}&d_4l_{42}&d_4l_{43}&d_4%5Cend{array}%5Cright))
Это выражение уже посложнее... Заметим, что изначально мы представляли матрицу A как LDLT, т.е нижнетреугольная умножалась на верхнетреугольную. Здесь, из-за транспонирования, они поменялись местами. Как и прежде, результирующая матрица - симметричная, иначе и быть не могло. Но теперь самые простые выражения оказываются в правом нижнем углу, а к левому верхнему усложняются.
Снова отбросим верхний треугольник, чтобы не мешался - его мы всё равно считать не будем (незачем). Кроме того, мы отметили зелёным, значения, которые стоят на этом месте изначально:
)
В предыдущих частях мы всё время начинали от простого к сложному: сначала элемент вообще не требуется менять, потом нужно вычесть один множитель, потом два - и так далее. Здесь тоже поначалу кажется, что начать стоит с правого нижнего угла (так, d4 уже стоит на своём месте!). Но так не получится - куда бы мы ни двинулись из правого нижнего угла - начнём перезаписывать значения, которые требуются для других вычислений!
Как ни странно, здесь надо, как и прежде, начать с левого верхнего угла, хоть выражение для него и самое длинное.
Но если присмотреться, можно заметить: двигаясь в естественном порядке, слева направо, сверху вниз (или как в ГОСТах просят нумеровать элементы на принципиальных схемах: сверху вниз, слева направо), мы каждый раз будем перезаписывать элемент, который в дальнейшем нам не понадобится!
Так что остаётся лишь "формализовать", как у нас меняются элементы - и алгоритм готов.
Проще всего обрабатываются диагональные элементы:
От текущего элемента "откладываем" две линии. Первая - вертикальная, сразу под текущим элементом, и до низу. По ней берутся множители, которые возводятся в квадрат. Затем помножаем результат на элемент со второй линии - диагональной - и всё это складываем, ещё и добавив исходный элемент.
Всего для вычисления диагонального элемента под номером i (нумерация от 1 до N) уходит 2(N-i) умножений и (N-i) сложений. Всего на диагональные элементы наберётся N(N-1) умножений и N(N-1)/2 сложений.
Все остальные элементы формируются чуть сложнее: нужно помножать по 3 компоненты, из разных мест:
Один множитель - с диагонали, начиная с текущей строки и донизу. Второй - на текущем столбце, начиная с "самого же себя". Третий - на столбце, номер которого совпадает с номером текущей строки, причём первый раз вместо диагонального элемента берём единицу, а потом уже спускаемся вниз.
На обработку элемента на строке j требуется 2(N-j)+1 умножений и (N-j) сложений. Всего их на строке j содержится (j-1) штук. Итого, на обработку всех недиагональных элементов требуется ещё (N-1)N(2N-1)/6 умножений и (N-2)(N-1)N/6 сложений.
Полное количество операций (и по диагонали, и вне её): (N-1)N(2N+5)/6 умножений и (N-1)N(N+1)/6 сложений.
В случае N=4, получаем 26 умножений и 10 сложений. Для N=6, получится уже 85 умножений и 35 сложений.
Ещё из хорошего: все алгоритмы могут проходить элементы в одном и том же порядке. Если вместо разложения LDLT взять разложение UDUT, получится всё то же самое, но теперь придётся идти не от первого элемента "слева направо сверху вниз", а наоборот - всегда из правого нижнего угла "справо налево снизу вверх". Это нам больше по душе при реализации на QuatCore (ℍ-core), поскольку циклы естественным образом организуются от максимального значения вниз до нуля.
В следующем посте немного подведём итоги и сравним с методом "в лоб".
Чтобы обратить её, мы сначала сделали LDLT-разложение "на месте" (хотя можно и на новое место, если исходная матрица ещё пригодится).
Затем обратили матрицы L и D, также "на месте".
Кстати, матрицы L-1 с D-1 и сами по себе представляют интерес: если A-1, которую мы ищем - это ковариационная матрица, то L-1 совместно с D-1 позволяет генерировать вектор случайных величин, удовлетворяющий этой ковариационной матрице. Хотя LLT-разложение в этой задаче, наверное, чуть лучше подходит...
Наконец, нужно перемножить матрицы, чтобы найти обратную матрицу:
Посмотрим, как можно произвести это умножение "за один этап", и притом по-прежнему "на месте".
Выкинем "-1" (обращение) - неудобно его за собой таскать. Как бы рассмотрим отдельную задачу: как посчитать выражение
для нижней унитреугольной матрицы L и диагональной матрицы D.
Выпишем такое умножение в явном виде для матриц 4×4. (почему мы выбрали 4×4 - потому что это минимальный размер матрицы, где начинает проглядываться весь алгоритм целиком, все 3 вложенных цикла. До этого кажется, что циклов всего 2). Слева - мы уже помножили диагональную матрицу на нижнюю треугольную, справа - итоговый результат:
Это выражение уже посложнее... Заметим, что изначально мы представляли матрицу A как LDLT, т.е нижнетреугольная умножалась на верхнетреугольную. Здесь, из-за транспонирования, они поменялись местами. Как и прежде, результирующая матрица - симметричная, иначе и быть не могло. Но теперь самые простые выражения оказываются в правом нижнем углу, а к левому верхнему усложняются.
Снова отбросим верхний треугольник, чтобы не мешался - его мы всё равно считать не будем (незачем). Кроме того, мы отметили зелёным, значения, которые стоят на этом месте изначально:
В предыдущих частях мы всё время начинали от простого к сложному: сначала элемент вообще не требуется менять, потом нужно вычесть один множитель, потом два - и так далее. Здесь тоже поначалу кажется, что начать стоит с правого нижнего угла (так, d4 уже стоит на своём месте!). Но так не получится - куда бы мы ни двинулись из правого нижнего угла - начнём перезаписывать значения, которые требуются для других вычислений!
Как ни странно, здесь надо, как и прежде, начать с левого верхнего угла, хоть выражение для него и самое длинное.
Но если присмотреться, можно заметить: двигаясь в естественном порядке, слева направо, сверху вниз (или как в ГОСТах просят нумеровать элементы на принципиальных схемах: сверху вниз, слева направо), мы каждый раз будем перезаписывать элемент, который в дальнейшем нам не понадобится!
Так что остаётся лишь "формализовать", как у нас меняются элементы - и алгоритм готов.
Проще всего обрабатываются диагональные элементы:
От текущего элемента "откладываем" две линии. Первая - вертикальная, сразу под текущим элементом, и до низу. По ней берутся множители, которые возводятся в квадрат. Затем помножаем результат на элемент со второй линии - диагональной - и всё это складываем, ещё и добавив исходный элемент.
Всего для вычисления диагонального элемента под номером i (нумерация от 1 до N) уходит 2(N-i) умножений и (N-i) сложений. Всего на диагональные элементы наберётся N(N-1) умножений и N(N-1)/2 сложений.
Все остальные элементы формируются чуть сложнее: нужно помножать по 3 компоненты, из разных мест:
Один множитель - с диагонали, начиная с текущей строки и донизу. Второй - на текущем столбце, начиная с "самого же себя". Третий - на столбце, номер которого совпадает с номером текущей строки, причём первый раз вместо диагонального элемента берём единицу, а потом уже спускаемся вниз.
На обработку элемента на строке j требуется 2(N-j)+1 умножений и (N-j) сложений. Всего их на строке j содержится (j-1) штук. Итого, на обработку всех недиагональных элементов требуется ещё (N-1)N(2N-1)/6 умножений и (N-2)(N-1)N/6 сложений.
Полное количество операций (и по диагонали, и вне её): (N-1)N(2N+5)/6 умножений и (N-1)N(N+1)/6 сложений.
В случае N=4, получаем 26 умножений и 10 сложений. Для N=6, получится уже 85 умножений и 35 сложений.
Ещё из хорошего: все алгоритмы могут проходить элементы в одном и том же порядке. Если вместо разложения LDLT взять разложение UDUT, получится всё то же самое, но теперь придётся идти не от первого элемента "слева направо сверху вниз", а наоборот - всегда из правого нижнего угла "справо налево снизу вверх". Это нам больше по душе при реализации на QuatCore (ℍ-core), поскольку циклы естественным образом организуются от максимального значения вниз до нуля.
В следующем посте немного подведём итоги и сравним с методом "в лоб".