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 假设输入是 Hermitian。这可能就是您认为您引用的 answer 不起作用的原因。

我不喜欢那个答案，因为它有一个迭代（而且，我无法理解它的例子），也不喜欢 other answer 那里它不承诺给你 best 正定矩阵，即根据 Frobenius 范数（元素的平方和）最接近输入的矩阵。（我完全不知道你问题中的代码应该做什么。）

我很喜欢 Higham 的 1988 论文的这个 Matlab 实现：https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd 所以我将它移植到 Python（编辑更新为 Python 3):

from numpy import linalg as la
import numpy as np


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 range(10):
        for j in range(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 和运行 python posdef.py 中，它会运行一个单元测试，在我的笔记本电脑上大约一秒钟就通过了。然后在您的代码中，您可以 import posdef 并调用 posdef.nearestPD 或 posdef.isPD.

如果您这样做，代码也在 Gist 中。

Answer 2

我知道这个帖子有点老了，但只想说@user1231818 link 提出的问题现在有了令人满意的答案，至少在我测试过的情况下是这样：

我将代码留在这里，但有关更多详细信息，请按照 link:

import numpy as np

def get_near_psd(A):
    C = (A + A.T)/2
    eigval, eigvec = np.linalg.eig(C)
    eigval[eigval < 0] = 0

    return eigvec.dot(np.diag(eigval)).dot(eigvec.T)

Answer 3

以防其他人遇到同样的问题。 @tjiagoM 和@Ahmed Fasih 上面解释的两种方法是等价的，应该都给出相同的 psd 矩阵，使 Frobenius 范数最小化。在论文的等式（2.1）和（2.2）中

'''
[2] N.J. Higham, "Computing a nearest symmetric positive semidefinite
matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6
'''

在 Ahmed 的回答中链接，描述了 tjiagoM 使用的计算方法，实际上被称为计算最接近的 psd 矩阵的首选方法。

Python: 将矩阵转换为半正定

Python: convert matrix to positive semi-definite

python

matrix

linear-algebra

svm