Pythonで密集したCholeskyの更新

Question

誰かが私にライブラリ/コードを教えて、コレスキー分解で低ランクの更新を実行できるようにすることができますか？ Matlabはこの機能を「cholupdate」という機能として提供しています。 LINPACKにもこの機能がありますが、LAPACKにまだ移植されていないため、（私の知る限り）LAPACKには移植されていません。サイフィ。

私はscikits.sparseがCHOLMODに基づいて同様の機能を提供することを知りましたが、私の行列は密集しています。

numupと互換性のある 'cholupdate'の機能を備えたPython用のコードはありますか？

ありがとうございます！

Answer 1

This guyは、scikitsとnumpy/scipyを使用して同様のことを行っています。ここで

Answer 2

はCythonを使用してコレスキー因子に関するランク1回の更新とdowndatesを行いPythonパッケージです： https://github.com/jcrudy/choldate

例：

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

これはnumpyのアレイ上のランク1更新またはdowndateを行う必要がありますRおよびxは、更新またはダウンデートに対応する記号「+」または「 - 」で示されます。 （WikipediaのページのMATLAB cholupdateから移植されました：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