如何计算最近的正半定矩阵? [英] How can I calculate the nearest positive semi-definite matrix?

问题描述

我从R来到Python，并尝试重现使用Python在R中做过的许多事情. R的Matrix库具有一个非常漂亮的函数，称为nearPD()，可找到与给定矩阵最接近的正半定(PSD)矩阵.尽管我可以编写一些东西，但是对于Python/Numpy来说是新手，如果已经有了一些东西，我对重新发明轮子并不感到兴奋.关于Python现有实现的任何提示?

I'm coming to Python from R and trying to reproduce a number of things that I'm used to doing in R using Python. The Matrix library for R has a very nifty function called nearPD() which finds the closest positive semi-definite (PSD) matrix to a given matrix. While I could code something up, being new to Python/Numpy I don't feel too excited about reinventing the wheel if something is already out there. Any tips on an existing implementation in Python?

推荐答案

我不认为有一个库可以返回您想要的矩阵，但是这里是近东正半定矩阵的只是为了好玩"编码Higham(2000)的算法

I don't think there is a library which returns the matrix you want, but here is a "just for fun" coding of neareast positive semi-definite matrix algorithm from Higham (2000)

import numpy as np,numpy.linalg

def _getAplus(A):
    eigval, eigvec = np.linalg.eig(A)
    Q = np.matrix(eigvec)
    xdiag = np.matrix(np.diag(np.maximum(eigval, 0)))
    return Q*xdiag*Q.T

def _getPs(A, W=None):
    W05 = np.matrix(W**.5)
    return  W05.I * _getAplus(W05 * A * W05) * W05.I

def _getPu(A, W=None):
    Aret = np.array(A.copy())
    Aret[W > 0] = np.array(W)[W > 0]
    return np.matrix(Aret)

def nearPD(A, nit=10):
    n = A.shape[0]
    W = np.identity(n) 
# W is the matrix used for the norm (assumed to be Identity matrix here)
# the algorithm should work for any diagonal W
    deltaS = 0
    Yk = A.copy()
    for k in range(nit):
        Rk = Yk - deltaS
        Xk = _getPs(Rk, W=W)
        deltaS = Xk - Rk
        Yk = _getPu(Xk, W=W)
    return Yk

在对本文示例进行测试时，它会返回正确答案

When tested on the example from the paper, it returns the correct answer

print nearPD(np.matrix([[2,-1,0,0],[-1,2,-1,0],[0,-1,2,-1],[0,0,-1,2]]),nit=10)
[[ 1.         -0.80842467  0.19157533  0.10677227]
 [-0.80842467  1.         -0.65626745  0.19157533]
 [ 0.19157533 -0.65626745  1.         -0.80842467]
 [ 0.10677227  0.19157533 -0.80842467  1.        ]]

这篇关于如何计算最近的正半定矩阵?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！