Актуализация на Dense Cholesky в Python

Може ли някой да ме насочи към библиотека/код, позволяващ ми да извършвам актуализации от нисък ранг на разлагане на Cholesky в python (numpy)? Matlab предлага тази функционалност като функция, наречена „cholupdate“. LINPACK също има тази функционалност, но тя (доколкото ми е известно) все още не е пренесена в LAPACK и следователно не е налична напр. scipy.

Разбрах, че scikits.sparse предлага подобна функция, базирана на CHOLMOD, но моите матрици са плътни.

Има ли наличен код за python с функционалност „cholupdate“, който е съвместим с numpy?

Благодаря!


python numpy matrix lapack factorization
person user1116403    schedule 26.12.2011    source източник

Отговори (3)


arrow_upward
5
arrow_downward

Ето пакет на Python, който класира актуализации и понижавания на 1-во място на факторите на Cholesky с помощта на Cython: 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))
person jcrudy    schedule 17.02.2013

arrow_upward
3
arrow_downward

Това трябва да направи актуализация от ранг-1 или понижаване на 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
person Andrew    schedule 23.04.2013

arrow_upward
1
arrow_downward

Този човек прави нещо подобно, използвайки scikits и numpy/scipy.

person Benjamin    schedule 03.01.2012