# Speeding up the evaluation of a product of 64 matrices, each 16x16, which takes 30 minutes

I need to calculate a very large product of matrices, (64 terms, each a 16x16 matrix). Specifically,

L = Simplify[
Normal[Series[
Dot @@ Table[MatrixExp[z a[i][R] Y[i]], {i, 1, 64}], {z, 0,
2}]]]; // AbsoluteTiming
Li = Simplify[
Normal[Series[
Dot @@ Table[
MatrixExp[-z a[(-1) (i - 64)][R] Y[(-1) (i - 64)]], {i, 0,
63}], {z, 0, 2}]]]; // AbsoluteTiming


where Y[m] are 16x16 matrices given as follows.

x = 8; nv = 8; Table[e[a, b] = Table[KroneckerDelta[a,c]KroneckerDelta[b, d], {c, 8 + nv}, {d, 8 + nv}], {a, 8 + nv}, {b, 8 + nv}]; Table[Y[i] = e[1, x + i] + e[x + i, 1], {i, nv}];Table[Y[i + nv] = e[2, x + i] + e[x + i, 2], {i, nv}]; Table[Y[i + 2 nv] = e[3, x + i] + e[x + i, 3], {i,nv}];Table[Y[i + 3 nv] = e[4, x + i] + e[x + i, 4], {i,nv}]; Table[Y[i + 4 nv] = e[5, x + i] + e[x + i, 5], {i,nv}]; Table[Y[i + 5 nv] = e[6, x + i] + e[x + i, 6], {i,nv}]; Table[Y[i + 6 nv] = e[7, x + i] + e[x + i, 7], {i,nv}]; Table[Y[i + 7 nv] = e[8, x + i] + e[x + i, 8], {i, nv}];


a[m][R] are just constants (unknowns). L and Li are inverse of each other, so their product is the 16x16 identity matrix.

The outputs of the first 2 commands are

{881.187, Null}
{1028.31, Null}


which is more than 30 minutes, taken together.

When the size of the product is smaller (42), with 13x13 matrices, the evaluation time is much shorter, but with the increased in matrix size as the number of terms, the time taken has increased by a lot.

Is there a way to significantly reduce the time for these operations? Possibly with the use of Activate/Inactive commands? I only have a vague idea of how Activate/Inactive works and I don't know how to apply those correctly to reduce the evaluating time.

Many thanks.

• What is a[i]? Please, always provide all relevant data. Commented Nov 22, 2019 at 9:36
• a[m][R] are just unknown variables. Commented Nov 22, 2019 at 9:40
• Then I am not surprised that this takes forever. What is the actual purpose that you try to achieve here? I mean: What are yo going to do with the result? Commented Nov 22, 2019 at 9:45
• L is a matrix encoding the 64 variables a[m] that are actually the scalars spanning a coset manifold of the form G/H, where the difference in dimensions between G and H is 64. L is a part of even more complicated tensors needed to calculate the scalar potential. Since it takes so long to calculate L, I now see that it's quite hopeless to calculate those more complicated tensors involving L. Commented Nov 22, 2019 at 9:55
• Unless I miss something, this looks like a matrix product that is entirely numeric, not symbolic. So use CUDA! Matrix multiplication is a great example for parallelization. Commented Nov 22, 2019 at 17:54

Since your interested in the series expansion of a dot product of matrix exponentials whose arguments are linear in the expansion variable, there is no need to actually perform the matrix exponential. Note that the derivative of such a matrix exponential is:

$$\frac{\partial e^{M x}}{\partial z}= M.e^{M x}$$

For example suppose the matrix is:

SeedRandom[1]
mat = RandomInteger[3, {3, 3}];


Check:

D[MatrixExp[mat x], x] == mat . MatrixExp[mat x] //FullSimplify


True

Let's define a function exp whose derivative satisfies the above relation. First, change a system option so that we can give a custom differentiation rule for exp:

old = OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"];
SetSystemOptions["DifferentiationOptions" -> "ExcludedFunctions" -> DeleteDuplicates @ Join[old, {exp}]];


Now, derivatives of exp don't evaluate:

D[exp[m x], x]


D[exp[m x], x]

So, we can define the desired derivative:

exp /: D[exp[m_ z_], z_] /; FreeQ[m, z] := m . exp[m z]
exp /: D[exp[m_ z_], {z_, n_Integer}] /; FreeQ[m, z] := Nest[D[#, z]&, exp[m z], n]


Check:

D[exp[m x], x]
D[exp[m x], {x, 2}]


m.exp[m x]

m.m.exp[m x]

Derivatives of dot products automatically work:

D[exp[m1 x] . exp[m2 x] . exp[m3 x], x]


m1.exp[m1 x].exp[m2 x].exp[m3 x] + exp[m1 x].m2.exp[m2 x].exp[m3 x] + exp[m1 x].exp[m2 x].m3.exp[m3 x]

It will be convenient to give exp a definition for a zero argument as well (since differentiating and setting the variable to 0 is how one computes a McLaurin series after all):

exp[0] = id;
id /: Dot[a___, id, b___] := Dot[a, b]


The matrix exponential of the zero matrix is the identity matrix, and the identity matrix drops out of dot expressions. Now, we are ready to consider finding the series of a dot product of matrix exponentials:

expr = exp[m1 x] . exp[m2 x] . exp[m3 x];
s = Sum[(D[expr, {x, n}] /. x->0) x^n/n!, {n, 0, 2}]


id + (m1 + m2 + m3) x + 1/2 x^2 (m1.m1 + 2 m1.m2 + 2 m1.m3 + m2.m2 + 2 m2.m3 + m3.m3)

Notice that the matrix exponential has completely disappeared. Let's check with some random matrices:

SeedRandom[1]
rules = {
id -> IdentityMatrix[3],
m1 -> RandomReal[1, {3, 3}],
m2 -> RandomReal[1, {3, 3}],
m3 -> RandomReal[1, {3, 3}]
};


r1 = Normal @ Series[
MatrixExp[m1 x] . MatrixExp[m2 x] . MatrixExp[m3 x] /. rules,
{x, 0, 2}
] //Chop


{{1. + 1.81073 x + 3.68536 x^2, 0.531297 x + 3.50467 x^2, 2.11866 x + 4.10715 x^2}, {0.739475 x + 1.83549 x^2, 1. + 0.795283 x + 2.21928 x^2, 1.75492 x + 2.37534 x^2}, {1.75799 x + 2.9816 x^2, 1.9764 x + 2.84889 x^2, 1. + 1.29941 x + 4.4811 x^2}}

Now for my version:

r2 = s /. rules


{{1 + 1.81073 x + 3.68536 x^2, 0.531297 x + 3.50467 x^2, 2.11866 x + 4.10715 x^2}, {0.739475 x + 1.83549 x^2, 1 + 0.795283 x + 2.21928 x^2, 1.75492 x + 2.37534 x^2}, {1.75799 x + 2.9816 x^2, 1.9764 x + 2.84889 x^2, 1 + 1.29941 x + 4.4811 x^2}}

Check:

r1 == r2


True

For a dot product of n terms, the second order series expansion is:

ser[matrices:{__?SquareMatrixQ}] := With[{id = IdentityMatrix[Length @ First @ matrices]},
id + Sum[m, {m, matrices}] x + Sum[m.m, {m, matrices}]x^2/2 + x^2 Total[Dot @@@ Subsets[matrices, {2}]]
]


Our previous result:

ser[{m1, m2, m3} /. rules]


{{1 + 1.81073 x + 3.68536 x^2, 0.531297 x + 3.50467 x^2, 2.11866 x + 4.10715 x^2}, {0.739475 x + 1.83549 x^2, 1 + 0.795283 x + 2.21928 x^2, 1.75492 x + 2.37534 x^2}, {1.75799 x + 2.9816 x^2, 1.9764 x + 2.84889 x^2, 1 + 1.29941 x + 4.4811 x^2}}

Using @Henrik's answer to find the Y matrices:

x = 8;
nv = 8;
esparse = SparseArray[
Flatten[Table[{a, b, a, b}, {a, 8 + nv}, {b, 8 + nv}], 1] -> 1.,
{1, 1, 1, 1} (8 + nv), 0.
];
Ysparse = Flatten[Table[esparse[[j, x + i]] + esparse[[x + i, j]], {j, 1, nv}, {i, 1, nv}], 1];


Then:

ser[Ysparse Array[a, 64]]; //AbsoluteTiming


{0.012494, Null}

should be what you're looking for.

• Thanks so much. This is indeed what I need ! There is just a tiny issue : I got this output (m1 + m2 + m3) z + Dot[] + 1/2 z^2 (m1.m1 + 2 m1.m2 + 2 m1.m3 + m2.m2 + 2 m2.m3 + m3.m3) when I ran the command for "expr", and instead of "id", I got "Dot[ ]", which later showed up in the output for the last command for "ser[Ysparse...]" as 1. all nonzero values. Commented Nov 23, 2019 at 8:58

Not directly the cause for the long computations, but the matrices e[a,b] and Y[i] can be built much faster with the following code:

esparse = SparseArray[
Flatten[Table[{a, b, a, b}, {a, 8 + nv}, {b, 8 + nv}], 1] -> 1.,
{1, 1, 1, 1} (8 + nv), 0.
];
Ysparse = Flatten[Table[esparse[[j, x + i]] + esparse[[x + i, j]], {j, 1, nv}, {i, 1, nv}], 1];


General good strategies for dealing with tensors:

1. If you want to store an array then store it in an array, not as DownValues of a symbol.

2. Try your best to store arrays in the most suitable array type: (i) packed array for numerical, machine type data (machine doubles, machine integers or machine complex doubles) or (ii) SparseArray for numerical or symbolic arrays with only few nonzero entries. The latter seems to be a good choice here.

1. When performing nonlinear tasks (such as MatrixExp), try to convert to machine precision as early as possible.
• Oh thanks, Carl. Indeed, it should be esparse. Commented Nov 22, 2019 at 17:31