I am trying to speed up my Mathematica code, so instead of testing and calculating the same thing again and again I want to save the result in a SparseArray. The goal is to generate a matrix which I can later use e.g. in a Matrix-Vector product. I am having trouble constructing the SparseArray - when trying to pass the indeces to a function I get a bunch of errors.
Rough description of the SparseArray $S$: The value at each position depends on the indexes of that position, $S$ looks somewhat like this:
SparseArray[{{jj_, kk_} /; Some test depending on jj and kk -> some value depending on jj and kk, {n, n}]
What works: My test is somewhat complicated, but it works fine (see below for the complete minimal example). I can generate a matrix like so:
SparseArray[{{jj_, kk_} /; Some test depending on jj and kk -> jj + kk, {n, n}]
e.g. for n=6 the MatrixForm looks like this:
$\begin{bmatrix}4&0&0&7&0&9 \\ 0&6&0&0&0&0\\ 0&0&8&0&0&0\\ 7&0&0&10&0&12\\ 0&0&0&0&12&0\\ 9&0&0&12&0&14\end{bmatrix}$
What does not work: Instead of just adding the two indexes $jj$ and $kk$ I want to pass them to a function getSval and use the value that this function returns, i.e.:
SparseArray[{{jj_, kk_} /; Some test depending on jj and kk -> getSval[degree,jj,kk], {n, n}]
getSval works fine when I call it outside the SparseArray definition, e.g. getSval[2, 4, 6] evaluates to $4\sqrt{2}\pi^{3/2}$. But using it in SparseArray throws a bunch of errors.
Here is the minimal working example:
(*Define IDX, essentially a list of indexes*)
Do[
IDX[n] =
Flatten[Table[
Table[{n - ii, ii - jj, jj}, {jj, 0, ii}], {ii, 0, n}], 1];
, {n, 0, 40}]
(*define the function getSval*)
getSval[degree_, j_, k_] := Block[{a, b, c, idx1, idx2},
idx1 := IDX[degree][[j]];
idx2 := IDX[degree][[k]];
a = 1/2 (idx1[[1]] + idx2[[1]]);
b = 1/2 (idx1[[2]] + idx2[[2]]);
c = 1/2 (idx1[[3]] + idx2[[3]]);
\[Pi]^(3/2) 2^(3/2)(a+b+c)
]
(*choose some setup-parameters for S*)
degree = 2;
length = Length[IDX[degree]]
(*try out the SparseArray function*)
MatrixForm[
SparseArray[{{jj_, kk_} /;
EvenQ[IDX[degree][[jj]][[1]] + IDX[degree][[kk]][[1]]] &&
EvenQ[IDX[degree][[jj]][[2]] + IDX[degree][[kk]][[2]]] &&
EvenQ[IDX[degree][[jj]][[3]] + IDX[degree][[kk]][[3]]]
->
degree + jj + kk}, {length, length}]] (*this works just fine*)
(*trying to use SparseArray with getSval*)
MatrixForm[
SparseArray[{{jj_, kk_} /;
EvenQ[IDX[degree][[jj]][[1]] + IDX[degree][[kk]][[1]]] &&
EvenQ[IDX[degree][[jj]][[2]] + IDX[degree][[kk]][[2]]] &&
EvenQ[IDX[degree][[jj]][[3]] + IDX[degree][[kk]][[3]]]
->
getSval[degree, jj, kk]}, {length, length}]] (*this crashes*)
I have tried several different things, e.g. including the calculation of getSval directly into the SparseArray, but so far nothing works. My impression is that $jj$ and $kk$ are handled differently after the arrow ->. For example
MatrixForm[
SparseArray[{{jj_, kk_} /;
EvenQ[IDX[degree][[jj]][[1]] + IDX[degree][[kk]][[1]]] &&
EvenQ[IDX[degree][[jj]][[2]] + IDX[degree][[kk]][[2]]] &&
EvenQ[IDX[degree][[jj]][[3]] + IDX[degree][[kk]][[3]]] ->
IDX[degree][[jj]][[1]] }, {length, length}]]
crashes with the complaint, that the value specified by the rule should not be a list, even though for example IDX[degree][[2]][[1]] evaluates to a number (in this case $1$).
What is the reason for this and is there some way to fix this?
Edit: I found the following work-around, but I am not happy with it - it certainly is not efficient:
- Step 1: Hold the evaluation when setting up the SparseArray:
S = SparseArray[{{jj_, kk_} /;
EvenQ[IDX[degree][[jj]][[1]] + IDX[degree][[kk]][[1]]] &&
EvenQ[IDX[degree][[jj]][[2]] + IDX[degree][[kk]][[2]]] &&
EvenQ[IDX[degree][[jj]][[3]] + IDX[degree][[kk]][[3]]] ->
Hold[getSval[degree, jj, kk]]}, {length, length}];
- Step 2: ReleaseHold of the Hold-expressions in $S$. Unfortunately
ReleaseHold[S]does not work (why?). However, I can parse $S$ as a list, ReleaseHold of that list and then parse the result back into a SparseArray:
SparseArray[ReleaseHold[Normal[S]]]
Obviously, this takes quite some time - is there a more efficient way?
