Лапак: Проблемы факторизации матрицы Холецкого

Выпуск 1

Может ли кто-нибудь порекомендовать менее неуклюжий способ факторизации Холецкого в python? Особенно меня смущает последняя строчка.

SigmaSqrt = matrix(Sigma)
cvxopt.lapack.potrf(SigmaSqrt)
SigmaSqrt = matrix(np.tril(SigmaSqrt))

Выпуск 2

У меня проблема, что одна целая строка и столбец (например, все элементы в первой строке и все элементы в первом столбце) равны нулю, lapack терпит неудачу с ошибкой, что матрица не является положительно определенной. Каков наилучший способ справиться с этим?

В настоящее время я делаю это: (что кажется неудобным...)

try:
    SigmaSqrt = matrix(Sigma)
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
except ArithmeticError:
    SigmaSqrt = matrix(Sigma.ix[1:,1:])
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
    SigmaSqrt = sparse([[v0],[v0[1:].T, SigmaSqrt]])

python matrix lapack cvxopt
person ARF    schedule 16.04.2014    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Вы можете просто использовать numpy.linalg.cholesky. Также, если весь один столбец или вся одна строка являются нулями, матрица будет сингулярной, иметь хотя бы одно собственное значение, которое будет равно нулю, и, следовательно, не будет положительно определенной. Поскольку Холецкий определен только для матриц, которые являются «эрмитовыми (симметричными, если они действительны) и положительно определенными», это не сработает для него.

EDIT: "разобраться" с вашей проблемой зависит от того, чего вы хотите. Все, что вы делаете, чтобы заставить его работать, приведет к холески, которая не будет холески исходной матрицы. Если вы выполняете итеративный процесс, и его можно немного выдумать, если он уже симметричен, используйте numpy.linalg.eigvalsh, чтобы найти собственные значения матрицы. Пусть d будет самым отрицательным собственным значением. Затем установите A += (abs(d) + 1e-4) * np.indentity(len(A)). это сделает его положительно определенным.

EDIT: Это трюк, используемый в алгоритме Левенберга-Марквардта. Это ссылка на статью в Википедии о методе Ньютона, в которой он упоминается как статья о Левенберге. – Marquard не занимается этим. Кроме того, здесь на нем тоже есть бумага. По сути, это сдвинет все собственные значения на (abs(d) + 1e-4), что сделает их все положительными, что является достаточным условием положительной определенности матрицы.

person cdhagmann    schedule 16.04.2014
comment
Спасибо за помощь. Как называется процедура, которую вы описываете в своем редактировании, чтобы я мог читать, освежить свою линейную алгебру? - person ARF; 17.04.2014
comment
@ArikRaffaelFunke Добавлена ​​дополнительная информация. Это прием, используемый в алгоритме Левенберга-Марквардта. - person cdhagmann; 17.04.2014
comment
Кроме того, есть scipy.linalg.cho_factor и scipy.linalg.cho_solve, которые можно использовать в тандеме для решения Ax=b с помощью теста Холецкого. - person cdhagmann; 17.04.2014

arrow_upward
1
arrow_downward

Другим вариантом является модуль chompack: домашняя страница chompack chompack.cholesky Я сам использую его в сочетании с модуль cvxopt. Он отлично работает с (разреженными) матрицами от cvxopt.

person Angela    schedule 28.07.2014
comment
Вы случайно не знаете, есть ли разница в скорости для плотных матриц между chompack.cholesky и numpy.linalg.cholesky? Или вы используете chompack просто потому, что это удобно с матрицами cvxopt? - person ARF; 07.10.2014