的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如果你做到這一點。