Actualización densa de Cholesky en Python

¿Podría alguien dirigirme a 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 (que yo sepa) aún no se ha portado 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!

    Aquí hay un paquete de Python que hace actualizaciones de rango 1 y reduce las tasas 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)) 

    Esto debería hacer una actualización de rango 1 o una baja en las matrices de R y x con el signo ‘+’ o ‘-‘ correspondiente a la actualización o la baja. (Portado desde la base de datos de MATLAB 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 = xT 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 

    Este tipo está haciendo algo similar usando scikits y numpy / scipy.