我从 R 转到 Python,并尝试使用 Python 重现我在 R 中习惯做的一些事情。 R 的 Matrix 库有一个非常漂亮的函数,称为 nearPD()
,它可以找到最接近给定矩阵的半正定 (PSD) 矩阵。虽然我可以编写一些代码,但作为 Python/Numpy 的新手,如果已经有了一些东西,我不会对重新发明轮子感到太兴奋。关于 Python 中现有实现的任何提示?
最佳答案
我不认为有一个库可以返回你想要的矩阵,但这里有一个来自 Higham (2000) 的近东半正定矩阵算法的“只是为了好玩”的编码
import numpy as np,numpy.linalg
def _getAplus(A):
eigval, eigvec = np.linalg.eig(A)
Q = np.matrix(eigvec)
xdiag = np.matrix(np.diag(np.maximum(eigval, 0)))
return Q*xdiag*Q.T
def _getPs(A, W=None):
W05 = np.matrix(W**.5)
return W05.I * _getAplus(W05 * A * W05) * W05.I
def _getPu(A, W=None):
Aret = np.array(A.copy())
Aret[W > 0] = np.array(W)[W > 0]
return np.matrix(Aret)
def nearPD(A, nit=10):
n = A.shape[0]
W = np.identity(n)
# W is the matrix used for the norm (assumed to be Identity matrix here)
# the algorithm should work for any diagonal W
deltaS = 0
Yk = A.copy()
for k in range(nit):
Rk = Yk - deltaS
Xk = _getPs(Rk, W=W)
deltaS = Xk - Rk
Yk = _getPu(Xk, W=W)
return Yk
当对论文中的例子进行测试时,它返回正确答案
print nearPD(np.matrix([[2,-1,0,0],[-1,2,-1,0],[0,-1,2,-1],[0,0,-1,2]]),nit=10)
[[ 1. -0.80842467 0.19157533 0.10677227]
[-0.80842467 1. -0.65626745 0.19157533]
[ 0.19157533 -0.65626745 1. -0.80842467]
[ 0.10677227 0.19157533 -0.80842467 1. ]]
关于python - 如何计算最近的半正定矩阵?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10939213/