Лапак: проблемы факторизации матрицы Холецкого
Выпуск 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]])
2 ответа
Вы могли бы просто использовать numpy.linalg.cholesky
, Также, если все один столбец или все строки являются нулями, матрица будет сингулярной, будет иметь по крайней мере собственное значение, равное нулю, и, следовательно, не будет положительно определенным. Поскольку Холецкий определен только для матриц, которые являются "эрмитовыми (симметричными, если они действительными) и положительно определенными", он не будет работать для него.
РЕДАКТИРОВАТЬ: "справиться" с вашей проблемой зависит от того, что вы хотите. Все, что вы сделаете, чтобы заставить его работать, приведет к холески, которые не будут холесками оригинальной матрицы. If you are doing an iterative process and it can be fudge a little, if it is already symmetric then use numpy.linalg.eigvalsh
to find the eigenvalues of the matrix. Let d be the most negative eigenvalue. Затем установите A += (abs(d) + 1e-4) * np.indentity(len(A))
, this will make it positive definite.
EDIT: It is a trick used in the Levenberg–Marquardt algorithm. This links to a Wikipedia article on Newton's Method that mentions it as the article on Levenberg–Marquard is doesn't go into this.Also, here is a paper on it as well. Basically this will shift all the eigenvalues by (abs(d) + 1e-4)
which will thereby make them all positive, which is a sufficient condition to the matrix being positive definite.
Другой вариант - модуль chompack: домашняя страница chompack chompack.cholesky
Я использую это сам в сочетании с модулем cvxopt. Он отлично работает с (разреженными) матрицами из cvxopt.