# Improving large array manipulation

I am doing an iterative program, which uses large arrays/vectors and where the output of one round is used as an input for the next steps but works slowly when dimensions are bigger. In brief, we create firstly the complex matrix s (dimension dos x dos) using SparseArray

    dos=101;
s = SparseArray[{{i_, i_} :> 1. + I, Band[{2, 1}] -> 2.,
Band[{1, 2}] -> 2., Band[{dos + 1, 1}] -> 2.0,
Band[{1, dos + 1}] -> 2.0}, {dos*dos, dos*dos}, 0.];


And initial condition

    input=ConstantArray[1,{dos,dos}];
input=Flatten[input];
u1=input;


And finally, we do the iterative steps following some conditions/rules related to which position in the array is used

    lr = Table[i, {i, 1 + dos, dos*dos - dos, dos}];
ll = Table[i, {i, dos + dos, dos*dos - dos, dos}];
lt = Table[i, {i, 1, dos, 1}];
lb = Table[i, {i, dos*dos - dos + 1, dos*dos, 1}];
lll = Sort[Join[lr, ll, lt, lb]]; (*rules*)

steps=50;
Do[f = SparseArray[{{i_} :>
If[MemberQ[lll, i],
0, -2.0*u1[[i - 1]] - 2.0*u1[[i + 1]] - 2.0*u1[[i + dos]] -
2.0*u1[[i - dos]] + (1.+I)*u1[[i]]]}, {dos*dos}];
u1 = LinearSolve[s, f];, {j, 1, steps, 1}]


This takes around 15 sec on my computer but with increasing dimensions, the time needed grows quickly. I noticed that applying LinearSolve is not slow, but creating the array f is where almost all the time is consumed. Is there a way to create f more efficiently? I tried to use Compile but not good results. Thanks in advance

• That code throws errors on my machine; an extra bracket at the end plus some lists of unequal length in the SparseArray. Have you enabled Evaluation>Debugger to see the timings of each call?
Apr 15 at 22:33
• Thanks @Adam, I realized I did some mistakes when I wrote this question. I am editing it now with no errors Apr 15 at 22:43
• Why are you Doing the LinearSolve 50 times? Presumably j should be referenced somewhere in the loop. I think your construction of f can be simplified by looking at small cases, and I think SparseArray isn't helping at all since you're manually specifying even 0 elements. The speedup of SparseArray comes from specifying only a few nonzero ones.
Apr 15 at 23:43
• @Adam I tried to minimize the question. In every step j of my original code, I am saving each u1 on a table in order to plot it lately. Here I put 50 steps, but it could be 5000 or higher. About SparseArray not helping at all since I am manually specifying elements, I didn't know that, but for the moment is my fastest way to create f. Apr 16 at 0:02
• I also suggest to call sinv = LinearSolve[s]; once. This will compute and store an LU-factorization of s. In the loop, you can simply use u1 = sinv[f]; unstead of u1 = LinearSolve[s, f];. On my machine, this gets down the solve time from 3 seconds to only 0.06 seconds for the factorization and 0.08 seconds all the solves in the loop. Apr 16 at 7:20

First OP's implementation with timing on my machine:

First@AbsoluteTiming[
u1 = input;
Do[f = SparseArray[{{i_} :>
If[MemberQ[lll, i],
0, -2.0*u1[[i - 1]] - 2.0*u1[[i + 1]] - 2.0*u1[[i + dos]] -
2.0*u1[[i - dos]] + (1. + I)*u1[[i]]]}, {dos*dos}];
u1 = LinearSolve[s, f];, {j, 1, steps, 1}];
result0 = u1;
]


17.1247

The computation of f from u1 is linear and as such can be represented by the following sparse matrix:

A = Times[
SparseArray[Partition[Complement[Range[dos dos], lll], 1] -> 1. + 0. I, {dos dos}, 0. I],
SparseArray[{
Band[{1, 1}] -> (1. + 1. I)
, Band[{1, 1 + 1}] -> -2.0 + 0. I
, Band[{1, 1 + dos}] -> -2.0 + 0. I
, Band[{1 + 1, 1}] -> -2.0 + 0. I
, Band[{1 + dos, 1}] -> -2.0 + 0. I
},
{dos dos, dos dos}, 0. I
]
]; // AbsoluteTiming // First


0.053559

So, using also one-time factorization of s, OP's implementation can be rewritten into

First@AbsoluteTiming[
sinv = LinearSolve[s];
u1 = input;
Do[u1 = sinv[A.u1];, {j, 1, steps, 1}];
result1 = u1;
]


0.141947

That is almost a 100-fold speedup.

This can also be formulated a bit shorter with Nest (which is not faster than Do):

First@AbsoluteTiming[
sinv = LinearSolve[s];
result2 = Nest[sinv[A.#] &, input, steps];
]


0.143743

Both lead to exactly the same result:

Max[Abs[result1 - result2]]


Comparison with OP's result:

Max[Abs[result1 - result0]]


541.03

Looks huge, but the relative error is much smaller:

Max[Abs[result1 - result0]]/Max[Abs[result0]]


1.64591*10^-14

What to learn from this: Pattern matching is all nice, but i maximizes rather programming speed, not execution speed.

• thanks @Henrik Schumacher for your answer. I have several questions, but maybe they can be answered knowing the following, what is idx? (it is inside Partition when you define A) Apr 18 at 4:31
• Uh, sorry. I had forgotton to copy its definition from my notebook. I edited my answer; maybe it is clearer now. (The sparse vector in which idx is used is supposed to nullify a couple of rows of the banded matrix to which it is multiplied.) Apr 18 at 8:00

The following (and my comment) overlooks the fact that u1 and f are dependent. Copy code, don't rewrite it -- also don't over edit once copied!

When I look at Normal@f for various values of dos, I find the "edge" elements of f (considered as a dos by dos 2d grid) are 0 and the "inner" elements of f are 1.0I-7.0. Here, then, is a quicker way to construct such arrays

f=Join[
0&/@Range@dos,
Join@@(Append[Prepend[1.0I-7.0&/@Range[dos-2],0],0]&/@Range[dos-2]),
0&/@Range@dos]
]


There are definitely faster/more efficient ways to construct f, but I think this may be a case of over-early-optimization. If you're interested in situations for arbitrary rules lll, it may simply not be possible to speed things up.

A good rule of thumb is to present the computationally hardest example you can think of in the question (i.e. a randomized set of rules, although this is probably not what you're interested in). Of course, simple cases are useful for understanding the problem at hand. And then there's the necessity of splitting up your problem into small, single-question-size parts.

I think there's more that can be elucidated with further code exchange. The crossed out text detracts from this post.

• thanks @Adam for your answer. It is true that inner elements are 1.0I-7.0 but only in the first step. In the following steps, those values are changing because f is not fixed, it is changing in every step Apr 18 at 4:42
• @DGtrauko my bad. Copied code wrong, but my responses are a bit disrespectful. Unfortunately I don't think I can improve on Henrik's answer!