Mfrhtyj

I want to calculate the trace of something like

$qquadmathrm{Tr}(Gamma_RG_dGamma_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.

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 Approximiation (ACA) (see this paper by Mario Bebendorf for more details on the algorithm). Because this implementation uses LinearAlgebra`BLAS`IAMAX which was introduced with version 11.2, this won't run on older versions of Mathematica.

Needs["LinearAlgebra`BLAS`"]



ClearAll[ACACompression];



Options[ACACompression] = {

   "MaxRank" -> 50,

   "Tolerance" -> 1. 10^-4,

   "StartingIndex" -> 1

   };



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

  eps = $MachineEpsilon;

  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++;

     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 straigth-forward way and compare the result and its timing to an ACA followed by a careful evaluation of the trace. Notice that the order of the Dot-operations 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]

0.311

0.0031

8.14907*10^-10

By the way, let's check the accuracy of the ACA is actually far better:

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

1.33227*10^-15