# Calculating the trace of the product of several matrices

I want to calculate the trace of something like

$$\qquad\mathrm{Tr}(\Gamma_RG_d\Gamma_LG_d^\dagger)$$

In order to optimise my code, I found something like this Faster trace of product of two matrices, which greatly reduces calculation time for trace.

My question is - Is it possible to generalize this somehow to the above case? I tried looking at the Wikipedia with no luck

https://en.wikipedia.org/wiki/Trace_%28linear_algebra%29#Trace_of_a_product

Or is the most time efficient way just to say

$$\qquad M_1 = \Gamma_rG_d \qquad M_2=\Gamma_LG_d^\dagger$$

and then use the trick in the link and Wikipedia.

• That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations. Jan 20 '19 at 21:51

If one or more of these matrices have low rank, then one could try to employ low-rank factorizations. Of course, not all matrices have low rank. I just try to point our that low rank can be exploited nicely if present.

In order to demonstrate the effect, here is an implementation of Adaptive Cross Approximation (ACA) (see this paper by Mario Bebendorf for more details on the algorithm). Because this implementation uses LinearAlgebraBLASIAMAX which was introduced with version 11.2, this won't run on older versions of Mathematica.

Needs["LinearAlgebraBLAS"]

ClearAll[ACACompression];

Options[ACACompression] = {"MaxRank" -> 50, "Tolerance" -> 1. 10^-4,
"StartingIndex" -> 1,
"MaxPivotTrials" -> 100};

ACACompression::maxrank = "Warning: Computed factorization has maximal rank.";

ACACompression[A_?MatrixQ, opts : OptionsPattern[]] := ACACompression[i \[Function] A[[i]], j \[Function] A[[All, j]], opts];

ACACompression[row_, column_, OptionsPattern[]] :=
Module[{maxrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR, rank, σ, m, n, eps, maxiter}, eps = \$MachineEpsilon;
maxiter = OptionValue["MaxPivotTrials"];
maxrank = OptionValue["MaxRank"];
ϵ = OptionValue["Tolerance"];
k = 1;
ik = OptionValue["StartingIndex"];
vk = row[ik];
m = Length[vk];
v = ConstantArray[0., {maxrank, m}];
jk = IAMAX[vk];
v[[k]] = vk = vk/vk[[jk]];
uk = column[jk];
n = Length[uk];
u = ConstantArray[0., {maxrank, n}];
u[[k]] = uk;
test = uk.uk vk.vk;
norm2 = test;
maxrank = Min[m, n, maxrank];
While[test > ϵ^2 norm2 && k < maxrank,
k++;
ik = IAMAX[uk];
vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
jk = IAMAX[vk];
δ = vk[[jk]];
If[Abs[δ] < eps, w = uk;
δ = 0.;
iter = 0;
While[Abs[δ] < eps && iter < maxiter,
iter++;
w[[ik]] = 0.;
ik = IAMAX[w];
vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
jk = IAMAX[vk];
δ = vk[[jk]];
];
];
v[[k]] = vk = vk/δ;
uk = Subtract[column[jk], v[[;; k - 1, jk]].u[[;; k - 1]]];
u[[k]] = uk;
test = uk.uk vk.vk;
norm2 += test + Dot[u[[1 ;; k - 1]].uk, v[[1 ;; k - 1]].vk]
];
If[k == OptionValue["MaxRank"], Message[ACACompression::maxrank];];
{u[[1 ;; k]], v[[1 ;; k]]}]


Now let's consider the following setup. I assume that the matrix G has low rank; I ensure that by taking the squared distance matrix of a set of homogeneously distributed points. (I am well aware that $$G$$ in OP's example is probably thought of as a Gram matrix of an inner product and, as such, positive definite so that it may be not of low rank. Still, ΓR, ΓL may be of low rank.)

n = 2000;
ΓR = RandomReal[{-1, 1}, {n, n}];
ΓL = RandomReal[{-1, 1}, {n, n}];
G = DistanceMatrix[RandomReal[{-1, 1}, n]]^2;


Now, let's compute the trace in the straight-forward way and compare the result and its timing to an ACA followed by a careful evaluation of the trace. Notice that the order in which the Dot-operations are performed is crucial here.

a = Tr[ΓR.G.ΓL.Transpose[G]]; // RepeatedTiming // First
First@RepeatedTiming[
{u, v} = ACACompression[G, "Tolerance" -> 1. 10^-6];
b = Tr[(u.ΓR.Transpose[u]).(v.ΓL.Transpose[v])];
]

Abs[a - b]/Abs[a]


0.33

0.0033

4.71289*10^-16

Let's also check that the relative accuracy of ACA:

Max[Abs[Transpose[v].u - G]]/Max[Abs[G]]


3.33565*10^-16

Edit

Squared distance matrices in Euclidean space can actually factored explicitly:

FactoredSquaredDistanceMatrix[x_?VectorQ, y_?VectorQ] := {
{x^2, x, ConstantArray[1., {Length[x]}]},
{ConstantArray[1., {Length[y]}], -2. y, y^2}
}
FactoredSquaredDistanceMatrix[x_?MatrixQ, y_?MatrixQ] := {
Join[Partition[Total[x^2, {2}], 1], x,
ConstantArray[1., {Length[x], 1}], 2],
Join[ConstantArray[1., {Length[y], 1}], -2. y,
Partition[Total[y^2, {2}], 1], 2]
}


So the following can do the job without ever forming the squared distance matrix.

First@RepeatedTiming[
{u, v} = FactoredSquaredDistanceMatrix[x, x];
c = Tr[(u.\[CapitalGamma]R.Transpose[
u]).(v.\[CapitalGamma]L.Transpose[v])];
]
Abs[a - c]/Abs[a]

• In fact, LinearAlgebraBLASIAMAX[] is available in older versions, just undocumented. Jan 28 '20 at 11:38