Python: convert matrix to positive semi-definite -

i'm working on kernel methods, , @ point needed make non positive semi-definite matrix (i.e. similarity matrix) 1 psd matrix. tried approach:

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 

but fails if test resulting matrix following function:

def ispsd(a, tol=1e-8):     e,v = linalg.eigh(a)     return np.all(e >= -tol) 

i tried approach suggested in other related question (how can calculate nearest positive semi-definite matrix?), resulting matrix failed pass ispsd test.

do have suggestions on how correctly make such transformation correctly?

first thing i’d don’t use eigh testing positive-definiteness, since eigh assumes input hermitian. that’s why think answer reference isn’t working.

i didn’t answer because had iteration (and, couldn’t understand example), nor other answer there doesn’t promise give best positive-definite matrix, i.e., 1 closest input in terms of frobenius norm (squared-sum of elements). (i have absolutely no idea code in question supposed do.)

i matlab implementation of higham’s 1988 paper: https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd ported python:

from numpy import linalg la  def nearestpd(a):     """find nearest positive-definite matrix input      python/numpy port of john d'errico's `nearestspd` matlab code [1],     credits [2].      [1] https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd      [2] n.j. higham, "computing 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))     # above different [1]. appears matlab's `chol` cholesky     # decomposition accept matrixes 0-eigenvalue, whereas     # numpy's not. [1] uses `eps(mineig)` (where `eps` matlab     # `np.spacing`), use above definition. caveat: our `spacing`     # larger [1]'s `eps(mineig)`, since `mineig` on     # order of 1e-16, , `eps(1e-16)` on order of 1e-34, whereas     # `spacing` will, gaussian random matrixes of small dimension, on     # othe order of 1e-16. in practice, both ways converge, unit test     # below suggests.     = np.eye(a.shape[0])     k = 1     while not ispd(a3):         mineig = np.min(np.real(la.eigvals(a3)))         a3 += * (-mineig * k**2 + spacing)         k += 1      return a3  def ispd(b):     """returns true when input positive-definite, via cholesky"""     try:         _ = la.cholesky(b)         return true     except la.linalgerror:         return false  if __name__ == '__main__':     import numpy np     in xrange(10):         j in xrange(2, 100):             = np.random.randn(j, j)             b = nearestpd(a)             assert(ispd(b))     print('unit test passed!') 

in addition finding nearest positive-definite matrix, above library includes ispd uses cholesky decomposition determine whether matrix positive-definite. way, don’t need tolerances—any function wants positive-definite run cholesky on it, it’s absolute best way determine positive-definiteness.

it has monte carlo-based unit test @ end. if put in posdef.py , run python posdef.py, it’ll run unit-test passes in ~a second on laptop. in code can import posdef , call posdef.nearestpd or posdef.ispd.

the code in gist if that.