# Why is 'scipy.sparse.linalg.spilu' less efficient than 'scipy.linalg.lu' for sparse matrix?

I have a matrix B which is sparse and try to utilize a function scipy.sparse.linalg.spilu specialized for sparse matrix to factorize B. Could you please explain why this function is significantly less efficient than the function scipy.linalg.lu for general matrix? Thank you so much!

import numpy as np
import scipy.linalg as la
import scipy.sparse.linalg as spla
import time
from scipy import sparse
from scipy.sparse import csc_matrix
A = np.random.randint(100, size=(10000, 10000))
B = np.triu(A, -100)

start = time.time()
(P, L, U) = la.lu(B)
end = time.time()
print('Time to decompose B with lu =', end - start)

start = time.time()
mtx = spla.spilu(csc_matrix(B))
end = time.time()
print('Time to decompose B with spilu =', end - start)


The computation time is

Time to decompose B with lu = 4.7765138149261475
Time to decompose B with spilu = 14.165712594985962

• You are benchmarking the format conversion to a csc_matrix together with the algorithm; is that intended? Jun 1, 2020 at 18:38
• @FedericoPoloni If I do not use csc_matrix, there is a warning SparseEfficiencyWarning: splu requires CSC matrix format warn('splu requires CSC matrix format', SparseEfficiencyWarning). Jun 1, 2020 at 18:52
• @LAD Federico is talking about timing. You may want to move the conversion to csc_matrix outside of the timing procedure. Jun 1, 2020 at 19:38
• Ah I got it. Thank you for your clarification! Jun 1, 2020 at 19:39
• There is also the fact that you're creating a dense matrix that you then convert into the completely unsuitable CSC format. Jun 1, 2020 at 19:39

This particular effect is highly likely to come from parallelization.

In many setups, numpy will use multiple threads for invoked BLAS/LAPACK calls. In the default setting on my laptop (Mac OS, native Apple python):

('Time to decompose B with lu =', 9.530492067337036)
('Time to decompose B with spilu =', 20.418880939483643)


and the Activity Monitor shows multiple threads invoked during the lu call, and only a single one during spilu call.

After explicitly specifying the number of threads (notice, on Mac you have to do this as well) with an "overkill set of commands":

export MKL_NUM_THREADS=1


The timings changed for the vanilla dense LU:

('Time to decompose B with lu =', 25.678237199783325)
('Time to decompose B with spilu =', 21.03290104866028)


General comment:

scipy.sparse.linalg.spilu corresponds to sparse incomplete LU decomposition, which is usually used as a preconditioner. Consider using devoted sparse direct solvers in case it fits your needs better.

• Do you mean by "devoted sparse direct solvers" that I should figure out the special structure of my matrix and look for a solver dedicated for this structure? Jun 1, 2020 at 16:36
• @LAD not necesserily. I mean to use, say, scipy.sparse.linalg.spsolve. or scipy.sparse.linalg.splu if you need the factorization. As spilu is not intended to be used out of the preconditioner realm in most cases. I added a link to the post, which you might find helpful. Jun 1, 2020 at 16:42
• Do you have any suggestion on the module to factorize the matrix whose many rows are already of upper triangular form? $$\begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} & x_{5} \\ 0 & x_{22} & x_{23} & x_{24} & x_{25} \\ 0 & 0 & x_{23} & x_{34} & x_{35} \\ 0 & 0 & 0 & 0 & x_{45} \\ x_{51} & x_{52} & x_{53} & x_{54} & x_{55} \\ x_{61} & x_{62} & x_{63} & x_{64} & x_{65} \end{bmatrix}$$ Jun 1, 2020 at 16:46
• I've seen your question. I would first try to use standard sparse direct solvers. It might be good enough Because if you want to squeeze all juices, you might have to do a lot of manual work yourself. I would also look into doing direct interaction/configuring of PARDISO solver which I am the most familiar. Not that I know right away how to deal with your particular case, though. Jun 1, 2020 at 16:49
• Also, I would warn you that a triangular (or close to) matrix is not sparse enough to be efficient. You need much more sparsity to take advantage of it using sparse linear algebra. Otherwise, dense operations might be a better fit. Jun 1, 2020 at 16:50