nabbla (nabbla1) wrote,
nabbla
nabbla1

Category:

Обращение симм. пол. опр. матрицы - последний этап

У нас была матрица A - симметричная N×N, положительно определённая. В памяти хранится только N(N+1)/2 её коэффициентов.

Чтобы обратить её, мы сначала сделали 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), поскольку циклы естественным образом организуются от максимального значения вниз до нуля.


В следующем посте немного подведём итоги и сравним с методом "в лоб".
Tags: tex, математика, программки, работа
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 

  • 0 comments