23

我从 R 开始使用 Python，并尝试使用 Python 重现我习惯于在 R 中做的许多事情。R 的 Matrix 库有一个非常漂亮的函数，称为nearPD()找到最接近给定矩阵的正半正定 (PSD) 矩阵。虽然我可以编写一些代码，但作为 Python/Numpy 的新手，如果已经存在某些东西，我不会对重新发明轮子感到太兴奋。关于 Python 中现有实现的任何提示？

4

4 回答 4

34

我不认为有一个库可以返回你想要的矩阵，但这里是 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.        ]]
于 2012-06-07T21:40:03.380 回答
17

我会提交一个非迭代的方法。这是从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 中的非 PD/PSD 矩阵的讨论中修改的。

于 2013-08-30T22:03:00.547 回答
3

这可能是对考虑相关矩阵和协方差矩阵的 DomPazz 答案的愚蠢扩展。如果您要处理大量矩阵，它也会提前终止。

def near_psd(x, epsilon=0):
    '''
    Calculates the nearest postive semi-definite matrix for a correlation/covariance matrix

    Parameters
    ----------
    x : array_like
      Covariance/correlation matrix
    epsilon : float
      Eigenvalue limit (usually set to zero to ensure positive definiteness)

    Returns
    -------
    near_cov : array_like
      closest positive definite covariance/correlation matrix

    Notes
    -----
    Document source
    http://www.quarchome.org/correlationmatrix.pdf

    '''

    if min(np.linalg.eigvals(x)) > epsilon:
        return x

    # Removing scaling factor of covariance matrix
    n = x.shape[0]
    var_list = np.array([np.sqrt(x[i,i]) for i in xrange(n)])
    y = np.array([[x[i, j]/(var_list[i]*var_list[j]) for i in xrange(n)] for j in xrange(n)])

    # getting the nearest correlation matrix
    eigval, eigvec = np.linalg.eig(y)
    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)))
    near_corr = B*B.T    

    # returning the scaling factors
    near_cov = np.array([[near_corr[i, j]*(var_list[i]*var_list[j]) for i in xrange(n)] for j in xrange(n)])
    return near_cov
于 2016-07-25T12:50:54.197 回答
0

我知道这个线程很旧，但是这里提供的解决方案对我的协方差矩阵并不满意：转换后的矩阵看起来总是与原始矩阵大不相同（至少对于我测试的情况）。因此，基于此答案中提供的解决方案，我将在这里留下一个非常简单的答案：

import numpy as np

def get_near_psd(A):
    C = (A + A.T)/2
    eigval, eigvec = np.linalg.eig(C)
    eigval[eigval < 0] = 0

    return eigvec.dot(np.diag(eigval)).dot(eigvec.T)

这个想法很简单：我计算对称矩阵，然后进行特征分解以获得特征值和特征向量。我将所有负特征值归零并构造回矩阵，该矩阵现在将是半正定的。

为了完整起见，我留下一个简单的代码来检查矩阵是否使用 numpy 是半正定的（基本上检查所有特征值是否都是非负的）：

def is_pos_semidef(x):
    return np.all(np.linalg.eigvals(x) >= 0)
于 2020-07-28T09:39:53.170 回答