# Performance of a Do-Loop of variable length

If I'm changing a basis depending on a (non quadratic) matrix $\mathbf{M}\in\mathbb Z^{d\times k}$ and some $\mu\in\mathbb N^k$, i.e. i have an array a of the same dimensions as mentioned in $\mu$, that is also given as an input parameter and I want to adress another array b at the same time after multiplication with $\mathbf{M}$. Hence i do need some Do-loop to run through all values $t = (t_1,\ldots,t_k),\ 1\leq t_i \leq \mu_i$.

Let's say d=3; k=2; mM = {{1, 0}, {1, 1}, {1, 0}}; and mu={512,1024}

Then my approach (that works) is something like

a = ConstantArray[0, mu]; (* Given as input for the function usually*)
b = ConstantArray[0, mM.mu]; (*to be computed/ "filled" with values *)
t1 = AbsoluteTiming[
Do[
b[[Sequence @@ ( mM.Table[Subscript[t, j], {j, 1, k}])]]
= a[[Sequence @@ (Table[Subscript[t, j], {j, 1, k}]) ]];
,Evaluate[Sequence @@ Table[{Subscript[t, j], 1, mu[[j]]}, {j, 1, k}]]];
]


Where the arrays in reality are of course not only zeros, usually they're not even sparse and i'm doing some computation not only the matching.

But even for this (i think) relatively small input the time needed for evaluation is very long (65 seconds on my Core2Duo 2.4 Ghz / 8GB Ram / Mathematica 8.1), e.g. compared to a multivariate Fourier transform on a (which i do before that in less than a second).

This one gets even worse, if $\mathbf{M}$ “spreads“ afurther, but that might also be due to Memory, perhaps.

In order to get that faster, I tried to extract the Table to be an Argument of an anonymous function (# & [...]), but that didn't gain any measurable time. Timing gets even worse, if i'm switching to my real application due to multiple access on the array a inside of Do.

What I am searching for now is a time optimization of that Do loop, whose Range is kind of variable with respect to the input $\mu$ and any ideas on faster access on a and b.

• without further reading into your problem as i can't test anything right now anyway, maybe you can adjust your problem to use NestWhile or NestWhileList. i'm pretty confident that they are significantly faster than Do. – mincos Nov 20 '12 at 10:04
• I tried to simplify the problem given here, in my application the right hand side is a Computation (certain sum of elements from a, but all “array index adresses” of course involve the actual one (written as the Sequence@@Table stuff. Thanks for the hint, I'll take a look at NestWhile. – Ronny Nov 20 '12 at 10:08
• I'm sorry, I tried adjusting the problem, bu i can't see how to do that - maybe i'm thinking too much in these Nested Do loops, but you're suggesting to replace Do kind of with NestWhile and i can't see how to do that. – Ronny Nov 20 '12 at 13:40

## 1 Answer

Setting this up to use ReplacePart seems much faster than replacing parts in a loop. The function f here creates a list of replacement rules and then performs the replacement. The uses of Thread are there to appropriately expand out the rules and lists.

Note that this is pretty memory intensive. You won't be able to go much larger than your sample problem without exhausting your memory at 8GB.

f[a_, b_, mM_] :=
Block[{t, mu, k, d, elem, index, rng, rep},
mu = Dimensions[a];
k = Length[mu];
d = Length[mM];
elem = mM.Table[t[j], {j, 1, k}];
rng = Sequence @@ Table[{t[j], 1, mu[[j]]}, {j, 1, k}];
index = Partition[Flatten[Table @@ {elem, rng}], d];
rep = Thread[index -> Flatten[a]];
ReplacePart[b, rep]
]

f[a, b, mM]; // AbsoluteTiming

(*{4.680007, Null}*)

• Uff, looks nice, but I don't understand Lines 3 and 4. Absolutely nothing. Can you explain them a little bit more? so {t1, t2} is my $\mu$ (because as input it doesn't matter whether you know a or $\mu$) but i can't recognize for example my “transformation rule” encoded in $\mathbf{M}$. – Ronny Nov 20 '12 at 16:36
• Well from Print[index] even for very small examples, i can't recognize the $\mathbf{M}$ or any useful structure to work with (and feel a little dumb if you mention “immediately”) - I'm sorry. The matrix should also be a variable of the function, that's why. I'll take a closer look on Partition an MapThread, though, just running your code with the numbers i posted as an example, it takes (on my humble machine, see above) your 4 seconds (while mine is at 18, maybe something messed up the timing this morning). – Ronny Nov 20 '12 at 16:50
• Ouh, we missunderstood each other, my implementation from above takes 18 seconds, while your code reproduces the 4 seconds exactly as you mentioned also on my machine. – Ronny Nov 20 '12 at 17:10
• Okay, I think I got your idea, but i think to implement that for other values of $k,d$ (or even have them variable too) might be a challenge. – Ronny Nov 20 '12 at 17:15
• I've added a more general solution. – Andy Ross Nov 20 '12 at 19:54