Sparse matrix factorization of A = Bᵀ·B

Question

Let's say I have a sparse matrix A, and I want to compute a matrix B such that

B.T.dot(B) = A

Is there any function on scipy module that can do this? If not, is it possible in numpy?

Answer 1

This functionality is available in numpy and scipy for dense matrices and in scikits.sparse for sparse matrices, depending on your matrix. Is your matrix positive definite? Is it symmetric? If so, then you are essentially looking for a Cholesky decomposition.

>>> A = np.random.rand(100,100)  # Construct a dense matrix
>>> np.fill_diagonal(A,10) # Ensure the matrix is positive definite
>>> A = 0.5*(A+A.T) # Symmetrize the matrix
>>> B = np.linalg.cholesky(A) # Perform Cholesky decomposition 
>>> np.allclose(B.dot(B.T),A) # Verify
True

Normally B is a lower triangular matrix in numpy, but it is a good practice to be sure of this, as for an upper triangular matrix you need to change the multiplication order from B.dot(B.T) to B.T.dot(B). With the scipy version there is the keywork argument lower (which is False by default) that you can specify to always get the right order that you need.

>>> B = sp.linalg.cholesky(A,lower=False)
>>> np.allclose(B.T.dot(B),A)
True

Although this should not be a big deal. For sparse matrices, you can either transform your sparse matrix to dense matrix with todense() or toarray() (not recommended) or use the scikits.sparse module

>>> from scikits.sparse.cholmod import cholesky
>>> spA = csc_matrix(A)
>>> factor = cholesky(spA)
>>> spB = factor.L()
>>> np.allclose(spA.todense(),spB.dot(spB.T).todense()) # Just for verification
True

For semi-definite matrices, there is also a Cholesky type factorization available in the paper: Analysis of the Cholesky Decomposition of a Semi-definite Matrix, but not implemented in scipy or scikits, as far as I know.