19
$\begingroup$

I've got a fairly large SparseMatrix (~60,000 entries in a ~10,000-by-10,000 matrix) that depends on two parameters (im1 and im2). It's made from eight parts and is mostly banded, plus a top row.

Right now, I make it like this:

SetSystemOptions["SparseArrayOptions" -> "TreatRepeatedEntries" -> Total];

r1 = r2 = 1;
e = 0.1;
k1 = k2 = 50;
α12 = 0.5; α21 = 0.5;
m1 = m2 = 0.01;
n1max = 2 k1; n2max = 2 k2;

TM[im1_, im2_] := SparseArray[Join[
  Flatten[Table[
    {1 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} ->
    - r1 n1 - im1 - (r1 n1 (n1 + α12 n2)/k1) + m1 n1 
    - r2 n2 - im2 - (r2 n2 (α21 n1 + n2)/k2) + m2 n2 - e
  , {n1, 0, n1max}, {n2, 0, n2max}]],
  Flatten[Table[
    {2 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} -> r1 n1 + im1
  , {n1, 0, n1max - 1}, {n2, 0, n2max}]],
  Flatten[Table[
    {n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} -> r1 n1 (n1 + α12 n2)/k1 + m1 n1
  , {n1, 1, n1max}, {n2, 0, n2max}]],
  Flatten[Table[
    {1 + n1 + (n1max + 1) (1 + n2), 1 + n1 + (n1max + 1) n2} -> r2 n2 + im2
  , {n1, 0, n1max}, {n2, 0, n2max - 1}]],
  Flatten[Table[
    {1 + n1 + (n1max + 1) (-1 + n2), 1 + n1 + (n1max + 1) n2} -> r2 n2 (α21 n1 + n2)/k2 + m2 n2
  , {n1, 0, n1max}, {n2, 1, n2max}]],
  Table[
    {1 + (n1max + 1) n2, 1 + (n1max + 1) n2} -> -e
  , {n2, 1, n2max}],
  Flatten[Table[
    {1, 1 + n1 + (n1max + 1) n2} -> e
  , {n1, 1, n1max}, {n2, 0, n2max}]],
  Table[
    {1, 1 + (n1max + 1) n2} -> e
  , {n2, 1, n2max}]
  ], {(n1max + 1) (n2max + 1), (n1max + 1) (n2max + 1)}, 0.];

Sorry that it's a bit complicated but a more minimal example might not have the same performance tuning available.

I'll need to get a particular eigenvector of TM[im1, im2] for a number of {im1, im2} parameter values. Looking at timings, I see that setting up the matrix is 3X slower than finding the eigenvector:

RepeatedTiming[tm = TM[1., 2.];]
(* {0.32, Null} *)
RepeatedTiming[Eigenvectors[tm, -1, Method -> "Arnoldi"][[1]];]
(* {0.10, Null} *)

So, does anyone have ideas on how to speed up constructing TM?

I've had a look at @HenrikSchumacher's Q&A here but couldn't figure out if it was relevant to my problem or how to adapt it.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
22
$\begingroup$

The way the code is written, you can neither exploit packed arrays nor any vectorization. There are two major reasons:

  • Using Rule prevents using packed arrays since Rule enforces symbolic computations. It is more efficient to feed SparseArray with pat -> vals, where pat is a packed array of dimensions {nnz, 2} of integers and vals is a packed array of size {nnz} (nnz is the number of rules).

  • Using Tables; it is faster than filling arrays with Do or For loops but not as fast as using vectorized constructs such as Range, ConstantArray, or KroneckerProduct. When things become too complicated to be sorted out by KroneckerProduct, just compile the body of the Table into a CompiledFunction.

Just taking your code for specific values of im1, and im2. Let's see the timings.

SetSystemOptions["SparseArrayOptions" -> "TreatRepeatedEntries" -> Total];

r1 = r2 = 1.;
e = 0.1;
k1 = k2 = 50;
α12 = 0.5; α21 = 0.5;
m1 = m2 = 0.01;
n1max = 2 k1; n2max = 2 k2;
im1 = 10;
im2 = 10;
rules = Join[ 
     Flatten[Table[{1 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} -> -r1 n1 - im1 - (r1 n1 (n1 + α12 n2)/k1) + m1 n1 - r2 n2 - im2 - (r2 n2 (α21 n1 + n2)/k2) + m2 n2 - e, {n1, 0, n1max}, {n2, 0, n2max}]], 
     Flatten[Table[{2 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} -> r1 n1 + im1, {n1, 0, n1max - 1}, {n2, 0, n2max}]], 
     Flatten[Table[{n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2} -> (r1 n1 (n1 + α12 n2)/k1) + m1 n1, {n1, 1, n1max}, {n2, 0, n2max}]], 
     Flatten[Table[{1 + n1 + (n1max + 1) (1 + n2), 1 + n1 + (n1max + 1) n2} -> r2 n2 + im2, {n1, 0, n1max}, {n2, 0, n2max - 1}]], 
     Flatten[Table[{1 + n1 + (n1max + 1) (-1 + n2), 1 + n1 + (n1max + 1) n2} -> (r2 n2 (α21 n1 + n2)/k2) + m2 n2, {n1, 0, n1max}, {n2, 1, n2max}]],
     Table[{1 + (n1max + 1) n2, 1 + (n1max + 1) n2} -> -e, {n2, 1, n2max}], 
     Flatten[Table[{1, 1 + n1 + (n1max + 1) n2} -> e, {n1, 1, n1max}, {n2, 0, n2max}]], 
     Table[{1, 1 + (n1max + 1) n2} -> e, {n2, 1, n2max}]
     ]; // AbsoluteTiming // First
vals0 = Join[
     Flatten[Table[-r1 n1 - im1 - (r1 n1 (n1 + α12 n2)/k1) + m1 n1 - r2 n2 - im2 - (r2 n2 (α21 n1 + n2)/k2) + m2 n2 - e, {n1, 0, n1max}, {n2, 0, n2max}]],
     Flatten[Table[r1 n1 + im1, {n1, 0, n1max - 1}, {n2, 0, n2max}]], 
     Flatten[Table[(r1 n1 (n1 + α12 n2)/k1) + m1 n1, {n1, 1, n1max}, {n2, 0, n2max}]], 
     Flatten[Table[r2 n2 + im2, {n1, 0, n1max}, {n2, 0, n2max - 1}]], 
     Flatten[Table[(r2 n2 (α21 n1 + n2)/k2) + m2 n2, {n1, 0, n1max}, {n2, 1, n2max}]],
     Table[-e, {n2, 1, n2max}],
     Flatten[Table[e, {n1, 1, n1max}, {n2, 0, n2max}]],
     Table[e, {n2, 1, n2max} ]
     ]; // AbsoluteTiming // First
A0 = SparseArray[
     rules, {(n1max + 1) (n2max + 1), (n1max + 1) (n2max + 1)}, 0.]; //
   AbsoluteTiming // First

0.192017

0.101042

0.085255

Okay, here is

I optimized the creation of vals and pat only really crudely. There might be a good deal of more optimization available, in particular with more background info. (For example, you seem to do certain computations on a rectangular grid... Is there some geometry involved?)

First, a few compiled functions:

cf1 = Block[{n1, n2},
   With[{code = -r1 n1 - (r1 n1 (n1 + α12 n2)/k1) + m1 n1 - r2 n2 - (r2 n2 (α21 n1 + n2)/k2) + m2 n2 - e},
    Compile[{{n1, _Real}, {n2, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cf2 = Block[{n1, n2},
   With[{code = (r1 n1 (n1 + α12 n2)/k1) + m1 n1},
    Compile[{{n1, _Real}, {n2, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cf3 = Block[{n1, n2},
   With[{code = (r2 n2 (α21 n1 + n2)/k2) + m2 n2},
    Compile[{{n1, _Real}, {n2, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg1 = Block[{n1, n2},
   With[{code = {1 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2}},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg2 = Block[{n1, n2},
   With[{code = {2 + n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2}},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg3 = Block[{n1, n2},
   With[{code = {n1 + (n1max + 1) n2, 1 + n1 + (n1max + 1) n2}},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg4 = Block[{n1, n2},
   With[{code = {1 + n1 + (n1max + 1) (1 + n2), 1 + n1 + (n1max + 1) n2}},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg5 = Block[{n1, n2},
   With[{code = {1 + n1 + (n1max + 1) (-1 + n2), 1 + n1 + (n1max + 1) n2}},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

cg6 = Block[{n1, n2},
   With[{code = n1max},
    Compile[{{n1, _Integer}, {n2, _Integer, 1}},
     {0 n2 + 1, 1 + n1 + (code + 1) n2},
     CompilationTarget -> "C",
     RuntimeAttributes -> Listable,
     Parallelization -> True
     ]
    ]
   ];

Here the new timings:

n1range = Range[0., n1max];
n2range = Range[0., n2max];

pat = Join[
     Flatten[Transpose[cg1[n1range, n2range], {1, 3, 2}], 1],
     Flatten[Transpose[cg2[Most@n1range, n2range], {1, 3, 2}], 1],
     Flatten[Transpose[cg3[Rest@n1range, n2range], {1, 3, 2}], 1],
     Flatten[Transpose[cg4[n1range, Most@n2range], {1, 3, 2}], 1],
     Flatten[Transpose[cg5[n1range, Rest@n2range], {1, 3, 2}], 1],
     KroneckerProduct[(Rest@n2range), (n1max + 1) {1, 1}] + 1,
     Flatten[Transpose[cg6[Rest@n1range, n2range], {1, 3, 2}], 1],
     Transpose[{ConstantArray[1, n2max], Range[1 + (n1max + 1), 1 + (n1max + 1) n2max, (n1max + 1)]}]
     ]; // AbsoluteTiming // First


vals = Join[
     Flatten[cf1[n1range, n2range] - im1 - im2],
     Flatten[KroneckerProduct[r1 Most[n1range] + im1, ConstantArray[1., 1 + n2max]]],
     Flatten[cf2[Rest@n1range, n2range]],
     Flatten[ConstantArray[r2 Most[n2range] + im2, n1max + 1]],
     Flatten[cf3[n1range, Rest@n2range]],
     ConstantArray[-e, n2max],
     ConstantArray[e, n2max (n2max + 1)],
     ConstantArray[e, n2max]
     ]; // AbsoluteTiming // First
A = SparseArray[
     pat -> vals, {(n1max + 1) (n2max + 1), (n1max + 1) (n2max + 1)}, 
     0.]; // AbsoluteTiming // First
A == A0

0.002733

0.001045

0.006115

True

You see:

  • The actual assembly (the call to SparseArray) became 10 times faster.

  • The generation of vals becomes 100 times faster.

  • Overall, the code is about 30 times faster.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks! This is just the kind of thing I'd never come up with myself. $\endgroup$
    – Chris K
    Commented Jul 8, 2019 at 16:14
  • $\begingroup$ You're welcome. =) $\endgroup$
    – Henrik Schumacher
    Commented Jul 8, 2019 at 16:23
  • 3
    $\begingroup$ These all look like points that never get mentioned in the SparseArray Help. +1. $\endgroup$
    – berniethejet
    Commented Jul 9, 2019 at 14:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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