3
$\begingroup$

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?

$\endgroup$
2
  • 2
    $\begingroup$ What happens if you use RuleDelayed (:>) instead of Rule (->)? $\endgroup$ – Carl Woll Jul 7 '20 at 17:27
  • $\begingroup$ Ahhh - yes, this solved the issue. I honestly wasn't aware of RuleDelayed, thank you :-) $\endgroup$ – Andrea Jul 7 '20 at 18:34
3
$\begingroup$

As Carl pointed out, RuleDelayed :> is your friend in this case and fixes the issue -

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}]]

gives the wanted result.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.