我从R来的是Python,试图重现我在Python中使用R所做的许多事情。R的Matrix库具有一个非常漂亮的函数nearPD(),该函数可找到与给定矩阵最接近的正半定(PSD)矩阵。虽然我可以编写一些东西,但是对于Python / Numpy来说是新手,如果已经有了一些东西,我对重新发明轮子并不感到兴奋。关于Python现有实现的任何提示?
nearPD()
我不认为有一个库可以返回您想要的矩阵,但是这里是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. ]]