3
$\begingroup$

Consider the following function that generates a sparse diagonal matrix:

W[n_] := W[n] = Exp[(2 Pi I)/n]//N;

Clear[X];

X[n_, L_, power_ : 1][i_, chain_] := 
  X[n, L, power][i, chain] =
    DiagonalMatrix@
     SparseArray@
       Table[
         W[n]^(-IntegerDigits[a, n, 2 L][[chain*L - i + 1]] power),
         {a, 0, n^(2 L) - 1}
       ]

Due to the f[x_] := f[x] = SlowFunction definition, I expect this code to be much faster on a second run. Indeed, if I evaluate the following on my laptop several times

X[3, 7][1, 1] // Timing

I get $12$ seconds on the first run then around $2.8$ seconds on any evaluation after that. Clearly the memoization trick seems to have made this faster, but it is still much slower than it should be. For example if I run

a = X[3, 7][1, 1]
a // Timing

I get $10^{-6}$ seconds. Running X[3, 7][1, 1] the second time should also be this fast, since the matrix is already computed and saved. But it seems that instead it is still doing some computation.

Why does this happen, and how could I avoid it so I can take full advantage of memoization to speed up my repeated calculations?

$\endgroup$
2
  • 4
    $\begingroup$ I think the problem is with the missing optional argument in the memoized version. $\endgroup$
    – mikado
    Feb 11, 2021 at 6:57
  • $\begingroup$ I just tried to make the "power' argument not be optional and it worked! Is it possible to keep the argument as optional and still make it work? If you make your comment into an answer (and maybe elaborate on why the optional variable causes this) I will be happy to accept it. $\endgroup$
    – Heidar
    Feb 11, 2021 at 7:06

1 Answer 1

6
$\begingroup$

One way is to use multiple dispatch instead of optional arguments:

X[n_, L_, power_][i_, chain_] := X[n, L, power][i, chain] = 
  DiagonalMatrix@SparseArray@
    Table[W[n]^(-IntegerDigits[a, n, 2 L][[chain*L-i+1]]*power), {a, 0, n^(2L)-1}]
X[n_, L_][i_, chain_] := X[n, L, 1][i, chain]

Here the memoization always happens in the full-argument form. Missing arguments are filled in by the multiple dispatch and forwarded to the full-argument form via delayed assignment.

$\endgroup$
4
  • $\begingroup$ What's the purpose of the X[n, L, power][i, chain] = in your first line? It seems to work the same if you remove it, i.e., if you define the recursion directly without the intermediate expression: X[n_, L_, power_][i_, chain_] := DiagonalMatrix@SparseArray@Table[W[n]^(-IntegerDigits[a, n, 2 L][[chain*L - i + 1]]*power), {a, 0, n^(2 L) - 1}] $\endgroup$
    – theorist
    Mar 26, 2022 at 7:34
  • 1
    $\begingroup$ @theorist the question was about memoization. Yes it works without, but it will be much slower. $\endgroup$
    – Roman
    Mar 26, 2022 at 9:55
  • $\begingroup$ You're referring to the time to rerun the recursion, not the time to run the recursion initially, right? On a fresh kernel, X[3, 7, 1][1, 1] takes ~10 s to run on my computer irrespective of whether the recursion is defined with or without memoization. Rerunning X[3, 7, 1][1, 1] on the same kernel is of course much faster with memoization, because the rule assigning it a value has been cached. I've found memoization only speeds up a recursion on a fresh kernel if you are calling more than one prior value. E.g., f[n]=f[n-1]^2 won't be faster, but f[n]=f[n-1]^2+f[n-2] will be. $\endgroup$
    – theorist
    Mar 27, 2022 at 8:50
  • $\begingroup$ ...And this makes sense to me since, in computing a recursion that only calls one prior value, there's no speed-up with memoization because MMA only needs that one prior value, which is going to be stored anyways. It's when you need both, say, the last value f(n-1) and the one before it f(n-2) that memoization helps, since without memoization you would need to recompute the latter. What's curious is that memorization helps when I need to call a prior value twice, e.g., f(n) = f(n-1)^2 + f(n-1). That behavior I can't rationalize. [Maybe I should assemble this as a separate question...] $\endgroup$
    – theorist
    Mar 27, 2022 at 8: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.