CODEWITHSUNDEEP
pythonnumpyscipy
jcrudy
jcrudy

Reputation: 4061

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?

Upvotes: 4

Views: 2895

Answers (3)

jcrudy
jcrudy

Reputation: 4061

For the sake of completeness, the following code based on Warren's comment solves the problem:

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('F')
print V.flags['C_CONTIGUOUS']
print R.flags['F_CONTIGUOUS']
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

All I've done is change the underlying storage format of R from C order to Fortran order, which causes overwrite_a to behave as expected.

If your matrix is real valued you can do this without making a copy by simply passing the transpose to cholesky. I don't think the transpose fix will work if your matrix is complex-valued, since the transpose of a complex hermitian matrix is not generally equal to its conjugate transpose. However, if you make sure your matrix is using fortran order (R.flags['F_CONTIGUOUS'] == True) in the first place then you shouldn't have a problem. But, see Setting the *default* data order (C vs. Fortran) in Numpy for difficulties with that strategy.

Upvotes: 3

Stefano M
Stefano M

Reputation: 4809

If you are brave enough you can avoid scipy and go for a low-level call to linalg.lapack_lite.dpotrf

import numpy as np

# prepare test data
A = np.random.normal(size=(10,10))
A = np.dot(A,A.T)
L = np.tril(A)

# actual in-place cholesky
assert L.dtype is np.dtype(np.float64)
assert L.flags['C_CONTIGUOUS']
n, m = L.shape
assert n == m
result = np.linalg.lapack_lite.dpotrf('U', n, L, n, 0)
assert result['info'] is 0

# check if L is the desired L cholesky factor
assert np.allclose(np.dot(L,L.T), A)
assert np.allclose(L, np.linalg.cholesky(A))

You have to understand lapack's DPOTRF, fortran calling conventions, memory layout. Don't forget upon exit to check result['info'] == 0. Nevertheless you see its just a line of code, and by throwing away all sanity checks and copying done by linalg.cholesky this could also be more efficient.

Upvotes: 2

Jaime
Jaime

Reputation: 67427

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]]

Upvotes: 4

Related Questions