Using dense matrix is faster than sparse matrix in LinearSolve [closed]

Question

I'm doing FEM with Mathematica and have to deal with LinearSolve of a large matrix inside a Newton loop constantly. Hence I need to speed up the computation efficiency of LinearSolve. I'm aware of this post.

Similarly, I test with the following code. First I tried with dense matrix:

Clear["Global`*"];
dim = 5000;
a = RandomReal[{1, 2}, {dim, dim}];
b = RandomReal[{1}, {dim}];
Table[SetSystemOptions["ParallelOptions" -> "MKLThreadNumber" -> i];
  Print["Case=", i];
  Print[SystemOptions["ParallelOptions" -> "MKLThreadNumber"]];
  t = AbsoluteTime[];
  LinearSolve[a, b];
  time2 = AbsoluteTime[] - t;
  Print["t(", i, ")=", time2];
  Print["******"], {i, 4}];

The output is:

Case=1

{ParallelOptions->{MKLThreadNumber->1}}

t(1)=2.155357

******

Case=2

{ParallelOptions->{MKLThreadNumber->2}}

t(2)=1.279845

******

Case=3

{ParallelOptions->{MKLThreadNumber->3}}

t(3)=1.525075

******

Case=4

{ParallelOptions->{MKLThreadNumber->4}}

t(4)=1.246443

******

From the above test, it seems that the LinearSolve is automatically parallelized, which is good. My laptop has only 2 kennels hence no much improvement for 2,3, and 4 threads. Then I run it again with sparse matrix:

Clear["Global`*"];
dim = 5000;
a = RandomReal[{1, 2}, {dim, dim}];
b = RandomReal[{1}, {dim}];
sa = SparseArray[a];
sb = SparseArray[b];
Table[SetSystemOptions["ParallelOptions" -> "MKLThreadNumber" -> i];
  Print["Case=", i];
  Print[SystemOptions["ParallelOptions" -> "MKLThreadNumber"]];
  t = AbsoluteTime[];
  LinearSolve[sa, sb];
  time2 = AbsoluteTime[] - t;
  Print["t(", i, ")=", time2];
  Print["******"], {i, 4}];

However, it seems that LinearSolve needs much more time to deal with sparse matrix then with dense matrix.

Case=1

{ParallelOptions->{MKLThreadNumber->1}}

t(1)=11.346750

******

Case=2

{ParallelOptions->{MKLThreadNumber->2}}

t(2)=5.508803

******

Case=3

{ParallelOptions->{MKLThreadNumber->3}}

t(3)=5.551608

******

Case=4

{ParallelOptions->{MKLThreadNumber->4}}

t(4)=5.518430

******

This seems a little bit weird to me since I have experienced the improved efficiency by converting a dense matrix into sparse matrix in Matlab. Can anyone explain this? Is there any other way to speed up LinearSolve?

Answer 1

You cannot benefit from SparseArray if your matrix is not sparse! The data structure storing a SparseArray needs a certain amount of overhead which only pays off if your matrix is really, really sparse (density below a permille or so). Moreover, whenever you find yourself converting a large matrix to a sparse matrix you have actually done something wrong. SparseArrays have to be constructed in a sparse way by special constructor functions (the constructor functions are called SparseArray in Mathematica and sparse in MATLAB).

But admittedly, a point that is bugging me is that Mathematica uses UMFPACK for factorizations of SparseArrays (as backend of LinearSolve). While UMFPACK is quite robust, it is rather slow, extremely memory hungry, and not parallelized. However, we have the alternative to use the mostly undocumented option Method->"Pardiso" in LinearSolve which uses the Intel MKL Pardiso solver; and Pardiso is (partially) parallelized (back substitution is hardly parallelizable).