Performance of a Do-Loop of variable length

Question

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.

Answer 1

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}*)