# Efficiently defining a SparseArray function

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"][];]
(* {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.

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.

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