4
$\begingroup$

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?

$\endgroup$

closed as off-topic by Daniel Lichtblau, Carl Woll, LCarvalho, J. M. is away Nov 9 '17 at 14:18

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Daniel Lichtblau, Carl Woll, LCarvalho, J. M. is away
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 3
    $\begingroup$ Using a sparse array is effective when the array is actually sparse. You are converting a dense array into a sparse array that has a sparsity of 0, so I'm not surprised that sparse methods are slower. $\endgroup$ – Carl Woll Nov 8 '17 at 23:28
  • 2
    $\begingroup$ If sparse matrix methods were faster than methods for dense matrices, when working with dense matrices, then the dense methods would disappear. Voting to close on grounds that there is no clear reason to expect behavior any different from what was observed. $\endgroup$ – Daniel Lichtblau Nov 9 '17 at 0:00
  • $\begingroup$ As already noted, trying to use sparse methods on a demonstrably dense matrix is a fool's errand. You had chosen an especially poor example to present your premise. $\endgroup$ – J. M. is away Nov 9 '17 at 0:55
  • $\begingroup$ I'm voting to close this question as off-topic because the issued raised is not really a problem; it is arises from the OP's not having the necessary understanding of the Mathematica objects he is working with. $\endgroup$ – m_goldberg Nov 9 '17 at 2:09
  • 1
    $\begingroup$ With only 2 kennels there is an upper bound how many dogs you can have running efficiently in parallell. $\endgroup$ – mathreadler Nov 9 '17 at 6:54
13
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ Thank you for pointing that out. Do u imply that linearsolve is not parallelized by saying "While UMFPACK is quite robust, its is rather slow, extremely memory hungry and it is not parallelized"? $\endgroup$ – Wilhelm Nov 8 '17 at 23:57
  • $\begingroup$ The reason I said that is that I set SetSystemOptions["ParallelOptions" -> "MKLThreadNumber" -> 1] and SetSystemOptions["ParallelOptions" -> "MKLThreadNumber" -> 16], respectively. However, there is little difference of computing time. $\endgroup$ – Wilhelm Nov 9 '17 at 0:05
  • 2
    $\begingroup$ Yes, that's what I say. Still, LinearSolve will perform magnitudes better for SparseArrays arising in FEM than for the according densified matrices. $\endgroup$ – Henrik Schumacher Nov 9 '17 at 0:12
  • $\begingroup$ I think you may assume that the default parameters "ParallelOptions" are already optimized for your system. $\endgroup$ – Henrik Schumacher Nov 9 '17 at 0:13
  • $\begingroup$ If this is the case, how the computing time is decreased from 2.155357 to 1.279845 when changing MKLThreadNumber from 1 to 2?@Henrik Schumacher $\endgroup$ – Wilhelm Nov 9 '17 at 16:46

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