3
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ You are benchmarking the format conversion to a csc_matrix together with the algorithm; is that intended? $\endgroup$ – Federico Poloni Jun 1 at 18:38
  • $\begingroup$ @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). $\endgroup$ – LAD Jun 1 at 18:52
  • $\begingroup$ @LAD Federico is talking about timing. You may want to move the conversion to csc_matrix outside of the timing procedure. $\endgroup$ – Anton Menshov Jun 1 at 19:38
  • $\begingroup$ Ah I got it. Thank you for your clarification! $\endgroup$ – LAD Jun 1 at 19:39
  • $\begingroup$ There is also the fact that you're creating a dense matrix that you then convert into the completely unsuitable CSC format. $\endgroup$ – Wolfgang Bangerth Jun 1 at 19:39
5
$\begingroup$

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
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1
export VECLIB_MAXIMUM_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)

This answer can also be useful.

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – LAD Jun 1 at 16:36
  • 1
    $\begingroup$ @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. $\endgroup$ – Anton Menshov Jun 1 at 16:42
  • $\begingroup$ 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}$$ $\endgroup$ – LAD Jun 1 at 16:46
  • $\begingroup$ 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. $\endgroup$ – Anton Menshov Jun 1 at 16:49
  • 1
    $\begingroup$ 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. $\endgroup$ – Anton Menshov Jun 1 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.