您好, 欢迎来到 !    登录 | 注册 | | 设为首页 | 收藏本站

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

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

我要说的第一件事是不要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,因此,这是确定正定的绝对最佳方法

最后,它还具有基于Monte Carlo的单元测试。如果将其放入posdef.py并运行python posdef.py,它将在我的笔记本电脑上运行约一秒钟的单元测试。然后可以在代码import posdef调用posdef.nearestPDposdef.isPD

如果您这样做的话,代码要点

python 2022/1/1 18:37:17 有433人围观

撰写回答


你尚未登录,登录后可以

和开发者交流问题的细节

关注并接收问题和回答的更新提醒

参与内容的编辑和改进,让解决方法与时俱进

请先登录

推荐问题


联系我
置顶