如何计算最近的正半定矩阵？

Question

我是从R来到Python并尝试使用Python重现我以前在R中做过的一些事情.用于R的矩阵库具有非常漂亮的函数,其被称为nearPD()找到给定矩阵的最接近的正半定(PSD)矩阵.虽然我可以编写一些代码,但是对于Python/Numpy来说是新手,如果有什么东西已经存在,我不会觉得重新发明轮子太兴奋了.关于Python中现有实现的任何提示？

Answer 1

我不认为有一个库可以返回你想要的矩阵,但这里是来自Higham(2000)的neareast正半定矩阵算法的"只是为了好玩"编码

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.        ]]

Answer 2

我会提交一个非迭代的方法.这与Rebonato和Jackel(1999)(第7-9页)稍作修改.迭代方法可能需要很长时间才能处理超过几百个变量的矩阵.

import numpy as np

def nearPSD(A,epsilon=0):
   n = A.shape[0]
   eigval, eigvec = np.linalg.eig(A)
   val = np.matrix(np.maximum(eigval,epsilon))
   vec = np.matrix(eigvec)
   T = 1/(np.multiply(vec,vec) * val.T)
   T = np.matrix(np.sqrt(np.diag(np.array(T).reshape((n)) )))
   B = T * vec * np.diag(np.array(np.sqrt(val)).reshape((n)))
   out = B*B.T
   return(out)

代码是从这个话题的讨论,修改这里周围R. nonPD/PSD矩阵

Answer 3

对于那些仍然在这里结束的人，您现在可以使用statsmodels.stats.correlation_tools.cov_nearest