Но особую индийскую магию мы пока что оставили без деталей реализации, продемонстрировав лишь сам принцип.
Как оказывается, в этом методе всё-таки нужна вспомогательная память, длиной в одну строку (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, то бишь менять местами строки и столбцы, так, чтобы поместить на диагональ большие элементы.
