CODEWITHSUNDEEP
pythonsparse-matrixmatrix-factorization
Sgt. Pepper
Sgt. Pepper

Reputation: 177

sparse cholesky decomposition with interchanged rows and columns

I am using scikits.sparse.cholmod of python to get cholesky factorization of a symmetric matrix.

I compared the results from the cholesky() with matlab's chol(). The results have a disparity with some rows and columns interchanged. I am trying to iterate through the factorization to get the eigenvalues and this disparity seems to be problematic.

Here is my code:

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse import csc_matrix
from scikits.sparse.cholmod import cholesky

A = csr_matrix([[1,2,0,0], [0,0,3, 0], [4,0,5, 0], [0, 0, 1, 2]])
B = (A*A.T)
print "B: "
print B.todense()

for i in range(10):
    factor = cholesky(B.tocsc())
    l = factor.L()  #l is lower triangular
    B = (l.T*l)
    print l.todense()

And the lower triangular matrix for the first iteration is:

[[ 2.23606798  0.         0.          0.        ]
[ 0.          3.          0.          0.        ]
[ 0.          1.          2.          0.        ]
[ 1.78885438  5.          0.          3.57770876]]

And the matlab's lower triangular matrix is:

[2.2361        0         0         0
     0    3.0000         0         0
1.7889    5.0000    3.5777         0
     0    1.0000         0    2.0000]

The matlab result is plausible because it leads to correct eigenvalues. Am I doing something wrong with choosing the type of sparse matrix in python?

Upvotes: 4

Views: 1101

Answers (1)

ztik
ztik

Reputation: 3612

The cholesky algorithm is using a fill-reducing algorithm. Due to this it setups a permutation matrix P. So that LL'=PBP'.

You can refer to the factor documentation for more information.

If you print P you get:

>>> factor.P()
array([0, 1, 3, 2], dtype=int32)

Which is exactly the difference between the two matrices. The permutation of the last two rows and columns.

Upvotes: 1

Related Questions