¿Cómo puedo calcular la matriz semi-definida positiva más cercana?

Vengo a Python desde R y trato de reproducir una serie de cosas que estoy acostumbrado a hacer en R usando Python. La biblioteca de Matrix para R tiene una función muy ingeniosa llamada nearPD() que encuentra la matriz semi-definida positiva (PSD) más cercana a una matriz dada. Si bien podría codificar algo, ser nuevo en Python / Numpy no me entusiasma mucho reinventar la rueda si ya hay algo ahí afuera. ¿Algún consejo sobre una implementación existente en Python?

No creo que exista una biblioteca que devuelva la matriz que desea, pero aquí hay una encoding “solo por diversión” del algoritmo de matriz semi-definido positivo de Neareast de 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*QT 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 

Cuando se prueba en el ejemplo del papel, devuelve la respuesta correcta

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

Yo sometería un enfoque no iterativo. Esto se modifica ligeramente de Rebonato y Jackel (1999) (página 7-9). Los enfoques iterativos pueden tardar mucho tiempo en procesar matrices de más de unos pocos cientos de variables.

 import numpy as np def nearPSD(A,epsilon=0): n = A.shape[0] eigval, eigvec = np.linalg.eig(A) val = np.matrix(np.maximum(eigval,epsilon)) vec = np.matrix(eigvec) T = 1/(np.multiply(vec,vec) * val.T) T = np.matrix(np.sqrt(np.diag(np.array(T).reshape((n)) ))) B = T * vec * np.diag(np.array(np.sqrt(val)).reshape((n))) out = B*BT return(out) 

El código se modifica a partir de una discusión de este tema aquí en torno a matrices no PD / PSD en R.

Esta es quizás una extensión tonta de la respuesta de DomPazz para considerar las matrices de correlación y covarianza. También tiene una terminación temprana si está tratando con un gran número de matrices.

 def near_psd(x, epsilon=0): ''' Calculates the nearest postive semi-definite matrix for a correlation/covariance matrix Parameters ---------- x : array_like Covariance/correlation matrix epsilon : float Eigenvalue limit (usually set to zero to ensure positive definiteness) Returns ------- near_cov : array_like closest positive definite covariance/correlation matrix Notes ----- Document source http://www.quarchome.org/correlationmatrix.pdf ''' if min(np.linalg.eigvals(x)) > epsilon: return x # Removing scaling factor of covariance matrix n = x.shape[0] var_list = np.array([np.sqrt(x[i,i]) for i in xrange(n)]) y = np.array([[x[i, j]/(var_list[i]*var_list[j]) for i in xrange(n)] for j in xrange(n)]) # getting the nearest correlation matrix eigval, eigvec = np.linalg.eig(y) val = np.matrix(np.maximum(eigval, epsilon)) vec = np.matrix(eigvec) T = 1/(np.multiply(vec, vec) * val.T) T = np.matrix(np.sqrt(np.diag(np.array(T).reshape((n)) ))) B = T * vec * np.diag(np.array(np.sqrt(val)).reshape((n))) near_corr = B*BT # returning the scaling factors near_cov = np.array([[near_corr[i, j]*(var_list[i]*var_list[j]) for i in xrange(n)] for j in xrange(n)]) return near_cov