О реализации "индийского" обращения матрицы
Когда мы рассматривали "классический" способ обращения симметричной положительно определённой матрицы через LDLT-разложение, обращение нижней унитреугольной матрицы и умножение трёх матриц, мы особенно подчёркивали, в каком порядке надо проводить вычисления, чтобы все они были выполнены "на месте", не требуя дополнительной памяти вообще!
Но особую индийскую магию мы пока что оставили без деталей реализации, продемонстрировав лишь сам принцип.
Как оказывается, в этом методе всё-таки нужна вспомогательная память, длиной в одну строку (N элементов), может хватить и N-1 элемента.
И ещё не исключено, что индийцы пошли тропой Рамануджана, заново переоткрыв в 2011 году метод, который давным-давно был реализован на Фортране, просто никто туда "не лез"... Впрочем, доказывая мой исходный тезис про Write-only :)
Речь идёт об этом алгоритме: https://people.sc.fsu.edu/~jburkardt/f_src/asa007/asa007.html
Там есть ссылки как на фортрановский код, так и на сишный, портабельный донельзя, в смысле что вместо выделения памяти в N элементов, он получает указатель на Workspace в аргументе, потому как alloc/free у нас может и не быть :) (и выделение на стеке N элементов, где N идёт аргументом функции, почему-то не прижилось)
Это, очевидно, не тот алгоритм, что мы обсуждали ранее: он использует LLT разложение (разложение Холецкого, оно же Cholesky или просто Chol), но явно похож!
Первым делом вызывается отдельная процедура для LLT-разложения, а потом из матрицы L мы получаем A-1 в один присест.
У индийцев приведён и вариант с LLT-разложением. Там мы снова пытаемся решить уравнение
относительно неизвестной матрицы X. (E - это единичная матрица)
Заменяем A на её разложение:
Помножим обе части слева на L-1:
Именно это уравнение мы и будем решать, причём находить целиком обратную матрицу L-1 нам не потребуется - только диагональные элементы.
Выпишем поэлементно:
%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}(l^{-1})_{11}&0&0&0%5C%5C(l^{-1})_{21}&(l^{-1})_{22}&0&0%5C%5C(l^{-1})_{31}&(l^{-1})_{32}&(l^{-1})_{33}&0%5C%5C(l^{-1})_{41}&(l^{-1})_{42}&(l^{-1})_{43}&(l^{-1})_{44}%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_{34}%5C%5Cx_{41}&x_{42}&x_{43}&x_{44}%5Cend{array}%5Cright)=%5Cleft(%5Cbegin{array}{cccc}1/l_{11}&0&0&0%5C%5C-&1/l_{22}&0&0%5C%5C-&-&1/l_{33}&0%5C%5C-&-&-&1/l_{44}%5Cend{array}%5Cright))
И ещё одно наблюдение, которое верно как для LLT-разложения (как сейчас), так и для LDLT: можно находить элементов обратной матрице не по столбцам (от последнего к первому, и снизу вверх в каждом), а по строкам: от нижней к верхней, справа налево внутри строки.
Давайте попробуем. Сначала выписываем уравнения, соотв. каждой строке, потом их решения.
Нижняя строка самая простая:
отсюда
(сначала извлекали квадратный корень, теперь вот назад в квадрат возводим!)
Теперь предпоследняя строка, справа налево:
откуда
(мы опять воспользовались симметрией обратной матрицы)
Перемещаемся на строку выше:
откуда
И ещё на строку выше:
откуда:
Как можно заметить, когда подходит черёд посчитать очередной элемент обратной матрицы, в правой части всё для этого уже есть (посчитано ранее). Поэтому и такой метод имеет право на существование.
Можно заметить и ещё одну вещь: когда мы вычисляем x14, он же x41, который должен занять место l41, мы должны сохранить где-то l41, ведь это значение нам понадобится ещё всю строку! То же самое потом происходит с l31, l21, и разве что l11 можно тут же заменить на x11, что и даёт нам необходимое хранилище в N-1 элементов.
Если теперь взглянуть на ASA007, то увидим что-то до боли знакомое!
А именно: в этом алгоритме мы сначала вызываем подпрограмму для LLT-разложения. Затем, начинаем двигаться по строкам, снизу вверх, справа налево.
Очередную строку мы первым делом заносим во временное хранилище, w (workplace). Затем переменная x хранит текущее вычисляемое значение. Заносим в неё либо ноль (когда считаем элементы вне диагонали), либо 1/ljj, когда считаем диагональный элемент j. После этого начинаем вычислять "скалярное произведение" между заблаговременно занесёнными в w элементами строки и ранее посчитанными элементами, причём имеет место "отражение от главной диагонали":
На этом рисунке мы показали элементы, которые нужны для вычисления x14 (самая правая линия), x13 (центральная, "переломленная"), x12 (левая, которая сразу же из вертикальной стала горизонтальной).
И наконец, посчитав "скалярное произведение" (числитель наших выражений), мы делим его на элемент главной диагонали.
Получается, что именно этот метод здесь и реализован. Утверждается, что он был описан здесь:
Michael Healy,
! Algorithm AS 7:
! Inversion of a Positive Semi-Definite Symmetric Matrix,
! Applied Statistics,
! Volume 17, Number 2, 1968, pages 198-199.
но в каком виде - не знаю, найти сей документ не удалось.
Сам фортрановский код, а также сишный код, как указано в комментариях к нему, был последний раз изменён в 2008 году, за 3 года до выхода индийской статьи. Исходный код для Fortran77, скорее всего был выпущен сильно раньше, где-нибудь так в 77 году :)
И ещё замечание. Мы всё время говорили, что исходная матрица должна быть положительно определённой. Это не совсем так: для методов, основанных на LDLT-разложении матрица может быть любой (просто симметричной), но только в случае положительной определённости метод обладает хорошей устойчивостью. В противном случае рекомендуется делать pivoting, то бишь менять местами строки и столбцы, так, чтобы поместить на диагональ большие элементы.
Но особую индийскую магию мы пока что оставили без деталей реализации, продемонстрировав лишь сам принцип.
Как оказывается, в этом методе всё-таки нужна вспомогательная память, длиной в одну строку (N элементов), может хватить и N-1 элемента.
И ещё не исключено, что индийцы пошли тропой Рамануджана, заново переоткрыв в 2011 году метод, который давным-давно был реализован на Фортране, просто никто туда "не лез"... Впрочем, доказывая мой исходный тезис про Write-only :)
Речь идёт об этом алгоритме: https://people.sc.fsu.edu/~jburkardt/f_src/asa007/asa007.html
Там есть ссылки как на фортрановский код, так и на сишный, портабельный донельзя, в смысле что вместо выделения памяти в N элементов, он получает указатель на Workspace в аргументе, потому как alloc/free у нас может и не быть :) (и выделение на стеке N элементов, где N идёт аргументом функции, почему-то не прижилось)
Это, очевидно, не тот алгоритм, что мы обсуждали ранее: он использует LLT разложение (разложение Холецкого, оно же Cholesky или просто Chol), но явно похож!
Первым делом вызывается отдельная процедура для LLT-разложения, а потом из матрицы L мы получаем A-1 в один присест.
У индийцев приведён и вариант с LLT-разложением. Там мы снова пытаемся решить уравнение
относительно неизвестной матрицы X. (E - это единичная матрица)
Заменяем A на её разложение:
Помножим обе части слева на L-1:
Именно это уравнение мы и будем решать, причём находить целиком обратную матрицу L-1 нам не потребуется - только диагональные элементы.
Выпишем поэлементно:
Диагональные элементы обратной матрицы найти легко: это просто обратные значения исходной. А все элементы обратной матрицы в нижнем треугольнике посчитаем неизвестными - поставим на их месте прочерки:
И ещё одно наблюдение, которое верно как для LLT-разложения (как сейчас), так и для LDLT: можно находить элементов обратной матрице не по столбцам (от последнего к первому, и снизу вверх в каждом), а по строкам: от нижней к верхней, справа налево внутри строки.
Давайте попробуем. Сначала выписываем уравнения, соотв. каждой строке, потом их решения.
Нижняя строка самая простая:
отсюда
(сначала извлекали квадратный корень, теперь вот назад в квадрат возводим!)
Теперь предпоследняя строка, справа налево:
откуда
(мы опять воспользовались симметрией обратной матрицы)
Перемещаемся на строку выше:
откуда
И ещё на строку выше:
откуда:
Как можно заметить, когда подходит черёд посчитать очередной элемент обратной матрицы, в правой части всё для этого уже есть (посчитано ранее). Поэтому и такой метод имеет право на существование.
Можно заметить и ещё одну вещь: когда мы вычисляем x14, он же x41, который должен занять место l41, мы должны сохранить где-то l41, ведь это значение нам понадобится ещё всю строку! То же самое потом происходит с l31, l21, и разве что l11 можно тут же заменить на x11, что и даёт нам необходимое хранилище в N-1 элементов.
Если теперь взглянуть на ASA007, то увидим что-то до боли знакомое!
А именно: в этом алгоритме мы сначала вызываем подпрограмму для LLT-разложения. Затем, начинаем двигаться по строкам, снизу вверх, справа налево.
Очередную строку мы первым делом заносим во временное хранилище, w (workplace). Затем переменная x хранит текущее вычисляемое значение. Заносим в неё либо ноль (когда считаем элементы вне диагонали), либо 1/ljj, когда считаем диагональный элемент j. После этого начинаем вычислять "скалярное произведение" между заблаговременно занесёнными в w элементами строки и ранее посчитанными элементами, причём имеет место "отражение от главной диагонали":
На этом рисунке мы показали элементы, которые нужны для вычисления x14 (самая правая линия), x13 (центральная, "переломленная"), x12 (левая, которая сразу же из вертикальной стала горизонтальной).
И наконец, посчитав "скалярное произведение" (числитель наших выражений), мы делим его на элемент главной диагонали.
Получается, что именно этот метод здесь и реализован. Утверждается, что он был описан здесь:
Michael Healy,
! Algorithm AS 7:
! Inversion of a Positive Semi-Definite Symmetric Matrix,
! Applied Statistics,
! Volume 17, Number 2, 1968, pages 198-199.
но в каком виде - не знаю, найти сей документ не удалось.
Сам фортрановский код, а также сишный код, как указано в комментариях к нему, был последний раз изменён в 2008 году, за 3 года до выхода индийской статьи. Исходный код для Fortran77, скорее всего был выпущен сильно раньше, где-нибудь так в 77 году :)
И ещё замечание. Мы всё время говорили, что исходная матрица должна быть положительно определённой. Это не совсем так: для методов, основанных на LDLT-разложении матрица может быть любой (просто симметричной), но только в случае положительной определённости метод обладает хорошей устойчивостью. В противном случае рекомендуется делать pivoting, то бишь менять местами строки и столбцы, так, чтобы поместить на диагональ большие элементы.