Actualización de Cholesky denso en Python

Question

¿Alguien podría indicarme una biblioteca/código que me permita realizar actualizaciones de bajo rango en una descomposición de Cholesky en python (numpy)? Matlab ofrece esta funcionalidad como una función llamada 'cholupdate'. LINPACK también tiene esta funcionalidad, pero (hasta donde sé) todavía no se ha transferido a LAPACK y, por lo tanto, no está disponible en, por ejemplo, scipy.

Descubrí que scikits.sparse ofrece una función similar basada en CHOLMOD, pero mis matrices son densas.

¿Hay algún código disponible para python con la funcionalidad 'cholupdate' que sea compatible con numpy?

Gracias!

Answer 1

This guy está haciendo algo similar con scikits y numpy/scipy.

Answer 2

Aquí es un paquete Python que hace Posición 1 actualizaciones y downdates sobre los factores de Cholesky utilizando Cython: https://github.com/jcrudy/choldate

Ejemplo:

from choldate import cholupdate, choldowndate 
import numpy 

#Create a random positive definite matrix, V 
numpy.random.seed(1) 
X = numpy.random.normal(size=(100,10)) 
V = numpy.dot(X.transpose(),X) 

#Calculate the upper Cholesky factor, R 
R = numpy.linalg.cholesky(V).transpose() 

#Create a random update vector, u 
u = numpy.random.normal(size=R.shape[0]) 

#Calculate the updated positive definite matrix, V1, and its Cholesky factor, R1 
V1 = V + numpy.outer(u,u) 
R1 = numpy.linalg.cholesky(V1).transpose() 

#The following is equivalent to the above 
R1_ = R.copy() 
cholupdate(R1_,u.copy()) 
assert(numpy.all((R1 - R1_)**2 < 1e-16)) 

#And downdating is the inverse of updating 
R_ = R1.copy() 
choldowndate(R_,u.copy()) 
assert(numpy.all((R - R_)**2 < 1e-16)) 

Answer 3

Esto debería hacer una actualización de rango 1 o downdate en matrices numpy R yx con el signo '+' o '-' correspondiente a la actualización o al descenso. (Portado de MATLAB cholupdate en la página de Wikipedia: http://en.wikipedia.org/wiki/Cholesky_decomposition):

def cholupdate(R,x,sign): 
    import numpy as np 
     p = np.size(x) 
     x = x.T 
     for k in range(p): 
     if sign == '+': 
      r = np.sqrt(R[k,k]**2 + x[k]**2) 
     elif sign == '-': 
      r = np.sqrt(R[k,k]**2 - x[k]**2) 
     c = r/R[k,k] 
     s = x[k]/R[k,k] 
     R[k,k] = r 
     if sign == '+': 
      R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c 
     elif sign == '-': 
      R[k,k+1:p] = (R[k,k+1:p] - s*x[k+1:p])/c 
     x[k+1:p]= c*x[k+1:p] - s*R[k, k+1:p] 
     return R