1

I have the following piece of code where I compute the L matrix of a given square matrix using the command scipy.linalg.lu() and then I do the same thing again except then applied to the sparse form of the given matrix using scipy.sparse.linalg.splu(). This is the code:

import numpy as np
from scipy.sparse.linalg import splu
from scipy.sparse import csc_matrix
import scipy.linalg
A1 = csc_matrix([[1., 0, 0.], [5., 0, 2], [0, -1., 0]])
A2 = np.array([[1., 0, 0.], [5., 0, 2], [0, -1., 0]])
B = splu(A1)
P,L,U = scipy.linalg.lu(A2)
print(L);print(csr_matrix.todense(B.L))

Which returns the following:

[[ 1.   0.   0. ]
 [ 0.   1.   0. ]
 [ 0.2 -0.   1. ]]
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

As we can see, these matrices are not the same. Am I misunderstanding what both commands do or is something else going wrong? Any help is appreciated. Thanks!

Charlie Shuffler
  • 145
  • 7

1 Answers1

0

I think the answer here is that the "SuperLU" decomposition of a sparse matrix requires permutations of both the rows and columns (see the docs):

Pr * A * Pc = L * U

These are provided by the mapping of indices in the perm_r and perm_c attributes. So,

Pr = csc_matrix((3,3))
Pr[B.perm_r, np.arange(3)] = 1
Pc = csc_matrix((3,3))
Pc[np.arange(3), B.perm_c] = 1

(Pr.T @ B.U @ B.L @ Pc.T).A

gives, as required:

array([[ 1.,  0.,  0.],
       [ 5.,  0.,  2.],
       [ 0., -1.,  0.]])

The same as the non-sparse result which requires only permutation of the L-matrix, P @ L @ U.

xnx
  • 24,509
  • 11
  • 70
  • 109
  • So given my sparse matrix A1 and the sparse LU decomposition B = splu(A), could I somehow obtain the correct L matrix from B or are there other necessary steps I need to take for this? How would I have to treat B in order to still correctly "extract" the correct L matrix in sparse form? I hope I formulated my question clear enough, apologies. – Charlie Shuffler Dec 07 '18 at 09:25
  • I don't think this is possible in general: the optimization of the sparse matrix approach is such that the L, U matrices generated are not "compatible" with the dense representation. In a sense, the permutations are what keeps the representation sparse. It would help if you gave more information about *why* you think you need L explicitly. – xnx Dec 07 '18 at 09:36
  • I need the explicit sparse representation of L in order to solve a system Au = b where we decompose A using A = LU. If I don't get the proper representation of L I cannot get the right solution. – Charlie Shuffler Dec 07 '18 at 12:28
  • What’s wrong with [sparse.linalg.spsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html) for this? – xnx Dec 07 '18 at 12:54
  • There isn't anything wrong with it, it works. But as of right now I have to apply the LU decomposition on a non-sparse matrix, because splu doesn't return the right L matrix, which is very slow when the non-sparse matrix is large. I can gain significant speed if I could just make my non-sparse matrix sparse and then use splu and extract the correct L matrix. – Charlie Shuffler Dec 08 '18 at 09:43