Python:将矩阵转换为正半定数 [英] Python: convert matrix to positive semi-definite
问题描述
我目前正在研究核方法，有时需要将一个非正半定矩阵(即相似矩阵)制成一个PSD矩阵. 我尝试过这种方法:
I'm currently working on kernel methods, and at some point I needed to make a non positive semi-definite matrix (i.e. similarity matrix) into one PSD matrix. I tried this approach:
def makePSD(mat):
#make symmetric
k = (mat+mat.T)/2
#make PSD
min_eig = np.min(np.real(linalg.eigvals(mat)))
e = np.max([0, -min_eig + 1e-4])
mat = k + e*np.eye(mat.shape[0]);
return mat
但是如果我使用以下函数测试结果矩阵，它将失败:
but it fails if I test the resulting matrix with the following function:
def isPSD(A, tol=1e-8):
E,V = linalg.eigh(A)
return np.all(E >= -tol)
I also tried the approach suggested in other related question (How can I calculate the nearest positive semi-definite matrix?), but the resulting matrix also failed to pass the isPSD test.
您对如何正确进行此类转换有何建议?
Do you have any suggestions on how to correctly make such transformation correctly?
推荐答案
我要说的第一件事是不要使用eigh来测试正定性，因为eigh假定输入为Hermitian.这可能就是为什么您认为您引用的答案不起作用的原因.
First thing I’d say is don’t use eigh for testing positive-definiteness, since eigh assumes the input is Hermitian. That’s probably why you think the answer you reference isn’t working.
我不喜欢这个答案，因为它有一个迭代(而且我不明白它的示例)，也没有 other答案，它不能保证为您提供最佳正定矩阵，即按照Frobenius范数(元素的平方和)最接近输入的矩阵. (我完全不知道您的问题中的代码应该做什么.)
I didn’t like that answer because it had an iteration (and, I couldn’t understand its example), nor the other answer there it doesn’t promise to give you the best positive-definite matrix, i.e., the one closest to the input in terms of the Frobenius norm (squared-sum of elements). (I have absolutely no idea what your code in your question is supposed to do.)
我确实喜欢Higham的 1988 论文的Matlab实现: https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd ，所以我将其移植到了Python:
I do like this Matlab implementation of Higham’s 1988 paper: https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd so I ported it to Python:
from numpy import linalg as la
def nearestPD(A):
"""Find the nearest positive-definite matrix to input
A Python/Numpy port of John D'Errico's `nearestSPD` MATLAB code [1], which
credits [2].
[1] https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd
[2] N.J. Higham, "Computing a nearest symmetric positive semidefinite
matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6
"""
B = (A + A.T) / 2
_, s, V = la.svd(B)
H = np.dot(V.T, np.dot(np.diag(s), V))
A2 = (B + H) / 2
A3 = (A2 + A2.T) / 2
if isPD(A3):
return A3
spacing = np.spacing(la.norm(A))
# The above is different from [1]. It appears that MATLAB's `chol` Cholesky
# decomposition will accept matrixes with exactly 0-eigenvalue, whereas
# Numpy's will not. So where [1] uses `eps(mineig)` (where `eps` is Matlab
# for `np.spacing`), we use the above definition. CAVEAT: our `spacing`
# will be much larger than [1]'s `eps(mineig)`, since `mineig` is usually on
# the order of 1e-16, and `eps(1e-16)` is on the order of 1e-34, whereas
# `spacing` will, for Gaussian random matrixes of small dimension, be on
# othe order of 1e-16. In practice, both ways converge, as the unit test
# below suggests.
I = np.eye(A.shape[0])
k = 1
while not isPD(A3):
mineig = np.min(np.real(la.eigvals(A3)))
A3 += I * (-mineig * k**2 + spacing)
k += 1
return A3
def isPD(B):
"""Returns true when input is positive-definite, via Cholesky"""
try:
_ = la.cholesky(B)
return True
except la.LinAlgError:
return False
if __name__ == '__main__':
import numpy as np
for i in xrange(10):
for j in xrange(2, 100):
A = np.random.randn(j, j)
B = nearestPD(A)
assert(isPD(B))
print('unit test passed!')
除了只查找最接近的正定矩阵外，上述库还包括isPD，它使用Cholesky分解来确定矩阵是否为正定矩阵.这样，您就不需要任何公差-任何需要正定的函数都可以在其上运行Cholesky，因此，这是确定正定性的绝对最佳方法.
In addition to just finding the nearest positive-definite matrix, the above library includes isPD which uses the Cholesky decomposition to determine whether a matrix is positive-definite. This way, you don’t need any tolerances—any function that wants a positive-definite will run Cholesky on it, so it’s the absolute best way to determine positive-definiteness.
最后还有一个基于蒙特卡洛的单元测试.如果将其放在posdef.py中并运行python posdef.py，它将运行一个单元测试，该测试在我的笔记本电脑上通过了约一秒钟.然后，您可以在代码中import posdef并调用posdef.nearestPD或posdef.isPD.
It also has a Monte Carlo-based unit test at the end. If you put this in posdef.py and run python posdef.py, it’ll run a unit-test that passes in ~a second on my laptop. Then in your code you can import posdef and call posdef.nearestPD or posdef.isPD.
如果您这样做的话，代码也位于要点中.
The code is also in a Gist if you do that.
这篇关于Python:将矩阵转换为正半定数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
