How can I calculate the nearest positive semi-definite matrix?

Question

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?

Answer 1

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.        ]]