Python中的稠密Cholesky更新

Question

有没有人能告诉我，在Python（使用numpy）中，有什么库或者代码可以让我对Cholesky分解进行低秩更新？Matlab有一个叫做'cholupdate'的函数可以做到这一点。LINPACK也有类似的功能，但据我所知，它还没有被移植到LAPACK中，所以在像scipy这样的库里是找不到的。

我发现scikits.sparse提供了一个基于CHOLMOD的类似功能，但我的矩阵是稠密的。

有没有适合numpy的Python代码，可以实现'cholupdate'的功能？

谢谢！

Answer 1

这个人正在使用scikits和numpy/scipy做类似的事情。

Answer 2

这段代码应该可以对numpy数组R和x进行一次排名更新或降级，符号'+'表示更新，'-'表示降级。这个内容是从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

Answer 3

这里有一个Python包，它使用Cython来对Cholesky因子进行排名1的更新和降级操作。你可以在这个链接找到它：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))