的Python：转换矩阵半正定

Question

我目前正在核方法，并在某些时候我需要做出非半正定矩阵（即相似矩阵）转换成一个PSD矩阵。 我尝试这样的做法：

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 

，但如果我用下面的函数测试结果矩阵失败：

def isPSD(A, tol=1e-8): 
    E,V = linalg.eigh(A) 
    return np.all(E >= -tol) 

我也想尽了办法在其他相关的问题（How can I calculate the nearest positive semi-definite matrix?）建议，但结果矩阵也未能通过isPSD测试。

你有关于如何正确正确地做出这样的转变有什么建议？

Answer 1

我要说的第一件事就是不要用eigh来测试正定性，因为eigh假设输入是厄米特。这可能就是你认为你参考的answer不起作用的原因。

我不喜欢那个答案，因为它有一个迭代（并且，我不明白它的例子），也不是other answer有它不答应给你最好的正定矩阵，即，根据Frobenius范数（元素的平方和）最接近输入的那个。 （我绝对不知道你在你的问题的代码是应该做的。）

我喜欢这个matlab实现海厄姆的1988年纸：https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd所以我把它移植到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，所以它是确定正定性的绝对最佳方式。

它还具有在端部具有基于蒙特卡罗的单元测试。如果你把它放在posdef.py并运行python posdef.py，它会运行一个单元测试，在我的笔记本电脑上传递一秒钟。然后在你的代码，你可以import posdef，并呼吁posdef.nearestPD或posdef.isPD。

的代码也是一个Gist如果你做到这一点。