Python：将矩阵转换为正半定数

Python：将矩阵转换为正半定数

我要说的第一件事是不要

eigh

用于测试正定性，因为

eigh

假定输入为Hermitian。这可能就是为什么您认为引用的答案不起作用的原因。

我不喜欢那个答案，因为它有一个迭代（而且我不明白它的例子），也没有那里的另一个答案，它不能保证为您提供
最佳 的正定矩阵，即最接近正定矩阵的矩阵。以Frobenius范数（元素的平方和）表示的输入。（我完全不知道问题中的代码应该做什么。）

我喜欢Higham 1988年
论文的Matlab实现：https
:
//www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd，所以我将其移植到了Python：

from numpy import linalg as ladef nearestPD(A):    """Find the nearest positive-definite matrix to input    A Python/Numpy port of John D'Errico's `nearestSPD` MATLAB pre [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 A3def isPD(B):    """Returns true when input is positive-definite, via Cholesky"""    try:        _ = la.cholesky(B)        return True    except la.LinAlgError:        return Falseif __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，因此，这是确定正定的绝对最佳方法。

最后，它还具有基于Monte Carlo的单元测试。如果将其放入

posdef.py

并运行

pythonposdef.py

，它将在我的笔记本电脑上运行约一秒钟的单元测试。然后可以在代码中

importposdef

调用

posdef.nearestPD

或

posdef.isPD

。

如果您这样做的话，代码也要点。