# Compiling matrix multiplication using a list of matrices

I'm trying to create compiled function which computes a list of Gamma matrices named GammaM and multiply them within the compiled code.

cf = Compile[{{n, _Integer}},
Module[{M, GammaM},
M = Table[0.0 + 0.0 I, {i, 1, n}, {j, 1, n}];
GammaM = Table[gamma[i, n], {i, 1, n}];
M = GammaM[[1]].GammaM[[2]];
M],
CompilationTarget -> "C", RuntimeOptions -> "Speed"];


However I am unable to get a compiled function; the error I get is:

Compile::cplist: GammaM[[1]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>

Compile::cset: "Variable M of type !({\"_Complex\", 2}) encountered in assignment of type {!(\"_Integer\", 0)}."

The function gamma[i,n] is defined recursively and returns a complex matrix:

s3 = {{1.0, 0}, {0.0, -1.0}};
s2 = {{0.0, -I}, {I, 0.0}};
s1 = {{0.0, 1.0}, {1.0, 0.0}};
gamma[mu_, N_] :=
If[mu <= N - 2,
-KroneckerProduct[gamma[mu, N - 2] , s3],
If[mu == N - 1, KroneckerProduct[IdentityMatrix[2^((N - 2)/2)], s2],
KroneckerProduct[IdentityMatrix[2^((N - 2)/2)], s1]]
]


Any idea what I may be doing wrong will be appreciated.

• I don't think Compile is the way to go here. Before you even get to the problems with different types, KroneckerProduct isn't compilable. – dr.blochwave Jul 28 '16 at 10:59
• @blochwave: Thanks for your comment, it indeed works if I define GammaM outside the Compile command. However if I replace GammaM = Table[gamma[i,n], {i,1,n}] by GammaM = Table[SparseArray[gamma[i,n]], {i,1,n}], it doesn't work and I get back the same errors. I want to make them sparse in order to speed up the computations (which indeed are sped up by defining them as sparse even without the compiled function). – MvP Jul 28 '16 at 12:04
• I think what @blochwave wanted to say is that you won't get a speedup from compiling. gamma cannot be compiled (at least not in the current form), and the matrix multiplication will not see any speedup. – sebhofer Jul 28 '16 at 12:12
• In some sense you are correct; however by doing this naively you will most likely loose performance as you generate calls to the Mathematica kernel. This could be worked around (e.g. by inlining the definitions of GammaM using With). But even then I expect no speedup. Where should it come from? The matrix multiplication will certainly be not any faster just because you wrap it in compile. – sebhofer Jul 28 '16 at 13:06
• Since the gamma matrices are very sparse, make sure everything is sparse. Don't use compilation. Don't declare M ahead of time the way you are now, with Table. And make s1/s2/s3 sparse, so that the kronecker products also stay sparse. Finally, Consider using memoization (e.g. gamma[mu_,N_] := gamma[mu,N] = If[....]) for the gamma matrices, as otherwise that recursive call could be nasty. That is, if you have the memory available for memoization. – Alex Meiburg Jul 28 '16 at 16:19