Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to do in place Cholesky factorization in Python

I'm trying to do an in place Cholesky factorization, but either that feature is not actually implemented in scipy or there is something I am not understanding. I am posting here in case it is the latter. Here is a simplified example of what I'm doing:

import numpy
import scipy.linalg

numpy.random.seed(0)
X = numpy.random.normal(size=(10,4))
V = numpy.dot(X.transpose(),X)
R = V.copy()
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

I think what should happen is for R to be overwritten with an upper triangular matrix. However, when I run this code, my V and R come out identical (cholesky is not overwriting R). Am I misunderstanding the purpose of overwrite_a, or making some other mistake? If this is just a bug in scipy, is there a workaround or another way I can do an in place Cholesky factorization in Python?

like image 999
jcrudy Avatar asked Feb 18 '23 19:02

jcrudy


1 Answers

Try it again, with

>>> scipy.__version__
'0.11.0'
>>> np.__version__
'1.6.2'

it works perfectly:

X = np.random.normal(size=(10,4))
V = np.dot(X.transpose(),X)
print V
R = V.copy()
print R
C = scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

Output is:

[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # R before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V after
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 3.35003677  1.52419728  0.20973811  0.93910339] # R after
 [ 0.          2.57332704 -1.35946252  0.08375069]
 [ 0.          0.          2.33703292 -0.99382158]
 [ 0.          0.          0.          2.52478036]]
like image 114
Jaime Avatar answered Feb 20 '23 11:02

Jaime