Проблема с разложением Холецкого полуопределенных матриц заключается в том, что 1) он не уникален 2) алгоритм Краута не работает.
Существование разложения обычно доказывается неконструктивно с помощью ограничивающего аргумента (если M_n-> M и M_n = U ^ t_n U_n, тогда || U_n || = || M_n || ^ 1/2 где ||. || = норма Гильберта-Шмидта, а U_n ограниченная последовательность. Извлечение подпоследовательностичтобы найти предел U, удовлетворяющий U ^ t U = M и U треугольной формы.)
Я обнаружил, что в интересующих меня случаях было бы достаточно умножить диагональные элементы на 1 + эпсилон с эпсилономмаленький (возьмите несколько тысяч раз машинного эпсилона), чтобы получить вполне приемлемое разложение.
Действительно, если М положительно полуопределен, то для каждого эпсилона> 0 М + эпсилон I определенно положителен.
Поскольку схема сходится, когда эпсилон стремится к нулю, вы можете рассмотреть возможность вычисления разложения для нескольких эпсилонов и выполнить экстраполяцию Ричардсона.
Что касаетсяВ определенном случае вы могли бы реализовать алгоритм Crout самостоятельно (в Numeric Recipes есть пример кода), но я бы настоятельно рекомендовал не писать его самостоятельно и рекомендовать вместо этого использовать LAPACK.
Это может включать в себя оплату вашего босса за Intel MKL, если он обеспокоен потенциально плохой реализацией LAPACK.Большую часть времени я слышал такую речь, обоснование было «но мы не можем контролировать код, мы хотим написать его самостоятельно, чтобы мы могли отладить его в случае возникновения проблемы».Тупой аргумент.LAPACK 40 лет и тщательно протестирован.
Требовать не использовать LAPACK так же глупо, как требование не использовать стандартную библиотеку для синусов, косинусов и логарифмов.