5
$\begingroup$

I'm trying to numerically solve some equations using splitting operator method. The solver I construct iteratively constructs a matrix and feeds it to LinearSolve. After profiling with the solver, I found half of the time is spent on simple linear algebra operations before LinearSolve, and I would like to speedup these process if possible.

Consider this simple example of matrix operation

Needs["CompiledFunctionTools`"]
Needs["Experimental`"]

lth = 200;
mtx = RandomReal[{0, 1}, {lth, lth}];
ls = RandomReal[{0, 1}, {lth}];

Et = Function[{t}, Sin[(π t)/20] Sin[2 t]];

and I'm trying to construct a matrix IdentityMatrix[lth] + I (DiagonalMatrix[ls] + Et[t]*mtx), repeatly with different t:

Table[
       IdentityMatrix[lth] + I (DiagonalMatrix[ls] + Et[t]*mtx);, {t, 0., 
        20, 0.01}]; // AbsoluteTiming
(* {1.33086, Null} *)

As a comparison, LinearSolve is about 2X faster for solving a matrix equation with the same size

Table[
   LinearSolve[mtx, ls];, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.883315, Null} *)

So my goal is to speed up the matrix operation, especially considering my matrix operation is much simpler than what is required to solve a dense matrix equation in LinearSolve.

We can try to copmile the functions:

Etc = Compile[{{t, _Real}}, Et[t], 
   CompilationOptions -> {"InlineCompiledFunctions" -> True, 
     "InlineExternalDefinitions" -> True}];
Htc = Compile[{{t, _Real}}, 
  IdentityMatrix[lth] + I (DiagonalMatrix[ls] + Et[t]*mtx), 
  CompilationOptions -> {"InlineCompiledFunctions" -> True, 
    "InlineExternalDefinitions" -> True}];

and compare these different ways:

Table[
   IdentityMatrix[lth] + I (DiagonalMatrix[ls] + Et[t]*mtx);, {t, 0., 
    20, 0.01}]; // AbsoluteTiming
Table[IdentityMatrix[lth] + I (DiagonalMatrix[ls] + Etc[t]*mtx);, {t, 
    0., 20, 0.01}]; // AbsoluteTiming
Table[Htc[t];, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {1.33086, Null} *)
(* {1.29705, Null} *)
(* {2.26304, Null} *)

We have two observations:

  1. We gain only very little speedup in by compiling the function Et.
  2. We slow down almost 70% by compiling the whole matrix operations.

For the first point, I'm guessing this may due to the automatic compilation of Table. For the second point, we may look at the compiled code of the function Htc

"(*omitted*)
1   T(I2)0 = MainEvaluate[ Hold[IdentityMatrix][ I0]]
2   T(R2)2 = MainEvaluate[ Hold[DiagonalMatrix][ T(R1)1]]
3   R1 = R0
4   R5 = I1
5   R3 = Reciprocal[ R5]
6   R5 = R2 * R1 * R3
7   R3 = Sin[ R5]
8   R5 = I2
9   R5 = R5 * R1
10  R6 = Sin[ R5]
11  R3 = R3 * R6
12  T(R2)4 = R3 * T(R2)3
13  T(R2)2 = T(R2)2 + T(R2)4
14  T(C2)4 = C0 * T(R2)2
15  T(C2)2 = CoerceTensor[ I3, T(I2)0]]
16  T(C2)2 = T(C2)2 + T(C2)4
17  Return
"

We can see that it has two invokes to the MainEvaluate to compute the IdentityMatrix and DiagonalMatrix. So this slowdown may come from the repeatedly invoking of the MainEvaluator, although in the comment here Oleksandr R. pointed out that IdentityMatrix is called using opcode 47 which should give very small overhead.

We can try to remove the MainEvaluate completely by expanding the matrix using Silvia's construction (this construction is a workaround for a problem of OptimizedExpression):

Htc2 = 
  Compile[{{t, _Real}}, 
   Evaluate[
    OptimizeExpression[
      Hold[IdentityMatrix[lth] + 
          I (DiagonalMatrix[ls] + Et[t]*mtx)] /. Et -> EtTemp // 
       ReleaseHold] /. 
     EtTemp[x_] :> With[{val = Et[x]}, val /; True]]];

Table[Htc2[t];, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {5.30743, Null} *)

However, we see an almost 2X slowdown for some reasons.

We can also try to simplify the matrix operations manually by removing IdentityMatrix, and it gains almost 2X speedup

Table[ 
   DiagonalMatrix[1. + I ls] + I Etc[t]*mtx;, {t, 0., 20, 
    0.01}]; // AbsoluteTiming
(* {0.728859, Null} *)

But I'm sill wondering:

Is it possible to further speedup the matrix operations?


version: 10.4 on OS X 10.11.4.

Note: This is a follow-up question to the question here.

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.