Python: convierte la matriz a positivo semi-definido

Actualmente estoy trabajando en métodos de kernel, y en algún momento tuve que hacer una matriz semi-definida no positiva (es decir, una matriz de similitud) en una matriz de PSD. He intentado este enfoque:

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 

pero falla si pruebo la matriz resultante con la siguiente función:

 def isPSD(A, tol=1e-8): E,V = linalg.eigh(A) return np.all(E >= -tol) 

También probé el enfoque sugerido en otra pregunta relacionada ( ¿Cómo puedo calcular la matriz semidesfinida positiva más cercana? ), Pero la matriz resultante tampoco pasó la prueba de isPSD.

¿Tiene alguna sugerencia sobre cómo realizar correctamente dicha transformación correctamente?

Lo primero que diría es que no use eigh para probar el resultado positivo, ya que eigh asume que la entrada es Hermitian. Probablemente por eso crees que la respuesta a la que te refieres no funciona.

No me gustó esa respuesta porque tenía una iteración (y, no pude entender su ejemplo), ni la otra respuesta allí no promete darte la mejor matriz positiva-definida, es decir, la más cercana a la entrada en términos de la norma de Frobenius (sum cuadrada de elementos). (No tengo la menor idea de lo que su código en su pregunta debe hacer).

Me gusta esta implementación de Matlab del artículo de Higham de 1988 : https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd, así que lo porté a 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] NJ Higham, "Computing a nearest symmetric positive semidefinite matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6 """ B = (A + AT) / 2 _, s, V = la.svd(B) H = np.dot(VT, 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!') 

Además de encontrar la matriz positiva-definida más cercana, la biblioteca anterior incluye isPD que utiliza la descomposición de Cholesky para determinar si una matriz es positiva-definida. De esta manera, no necesita ninguna tolerancia, cualquier función que quiera una definición positiva se ejecutará en Cholesky, por lo que es la mejor manera de determinar la definición positiva.

También tiene una prueba unitaria basada en Monte Carlo al final. Si pones esto en posdef.py y ejecutas python posdef.py , se ejecutará una prueba de unidad que pasa en ~ un segundo en mi computadora portátil. Luego, en su código puede import posdef y llamar a posdef.nearestPD o posdef.isPD .

El código también está en un Gist si haces eso.