Python：行列を正の半定理に変換する

Question

私は現在、カーネルの方法を研究しています。ある時点で、非正準確定行列（すなわち類似行列）を1つの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の実装のようです（私はあなたの質問であなたのコードが行うことになっているものを絶対にないアイデアを持っていません。）

：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が含まれています。このようにすれば、公差を必要としない関数はコレスキーを実行するので、正の確定性を決める最善の方法です。

また、最後にモンテカルロベースの単体テストがあります。これをposdef.pyに入れてpython posdef.pyを実行すると、ラップトップで1秒後に通過する単体テストが実行されます。その後、あなたのコードでimport posdefとposdef.nearestPDまたはposdef.isPDと呼ぶことができます。

コードはGistにもあります。