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)

我也尝试了其他相关问题中建议的方法(如何计算最近的正半正定矩阵？),但得到的矩阵也未能通过isPSD测试.

对于如何正确地进行正确的转换,您有什么建议吗？

Answer 1

我要说的第一件事是不要eigh用于测试正定性,因为eigh假设输入是Hermitian.这可能就是为什么你认为你提到的答案不起作用的原因.

我不喜欢那个答案,因为它有一个迭代(并且,我无法理解它的例子),也没有其他答案,它不承诺给你最好的正定矩阵,即最接近的那个根据Frobenius范数(元素的平方和)的输入.(我完全不知道你的问题代码应该做什么.)

我喜欢Higham 1988年论文的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使用Cholesky分解来确定矩阵是否为正定矩阵.这样,你不需要任何容差 - 任何想要肯定的函数都会运行Cholesky,因此它是确定正定性的绝对最佳方法.

它最后还有一个基于蒙特卡罗的单元测试.如果你把它放入posdef.py并运行python posdef.py,它将运行一个单元测试,在我的笔记本电脑上传递〜秒.然后在你的代码中你可以import posdef打电话posdef.nearestPD或posdef.isPD.

如果你这样做,代码也在Gist中.