Calculating the trace of the product of several matrices












2












$begingroup$


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.










share|improve this question









New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.







$endgroup$








  • 1




    $begingroup$
    That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
    $endgroup$
    – Henrik Schumacher
    6 hours ago
















2












$begingroup$


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.










share|improve this question









New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.







$endgroup$








  • 1




    $begingroup$
    That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
    $endgroup$
    – Henrik Schumacher
    6 hours ago














2












2








2





$begingroup$


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.










share|improve this question









New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.







$endgroup$




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.







matrix performance-tuning






share|improve this question









New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.











share|improve this question









New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









share|improve this question




share|improve this question








edited 4 hours ago









m_goldberg

84.9k872196




84.9k872196






New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









asked 6 hours ago









ThomasThomas

111




111




New contributor




Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.





New contributor





Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.






Thomas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.








  • 1




    $begingroup$
    That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
    $endgroup$
    – Henrik Schumacher
    6 hours ago














  • 1




    $begingroup$
    That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
    $endgroup$
    – Henrik Schumacher
    6 hours ago








1




1




$begingroup$
That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
$endgroup$
– Henrik Schumacher
6 hours ago




$begingroup$
That depends on the matrices. If one or more of these matrices have low rank, then one could try to employ low-rank factorizations.
$endgroup$
– Henrik Schumacher
6 hours ago










1 Answer
1






active

oldest

votes


















3












$begingroup$

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, cutrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR,rank, σ, m, n, eps},
eps = $MachineEpsilon;
maxrank = OptionValue["MaxRank"];
cutrank = OptionValue["CutRank"];
ϵ = OptionValue["Tolerance"];
k = 1;
ik = OptionValue["StartingIndex"];
vk = row[ik];
m = Length[vk];
v = ConstantArray[0., {maxrank, m}];
jk = LinearAlgebra`BLAS`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 = LinearAlgebra`BLAS`IAMAX[uk];
vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
jk = LinearAlgebra`BLAS`IAMAX[vk];
δ = vk[[jk]];
If[Abs[δ] < eps,
w = uk;
δ = 0.;
iter = 0;
While[Abs[δ] < eps,
iter++;
w[[ik]] = 0.;
ik = LinearAlgebra`BLAS`IAMAX[w];
vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
jk = LinearAlgebra`BLAS`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.



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







share|improve this answer









$endgroup$













    Your Answer





    StackExchange.ifUsing("editor", function () {
    return StackExchange.using("mathjaxEditing", function () {
    StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
    StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
    });
    });
    }, "mathjax-editing");

    StackExchange.ready(function() {
    var channelOptions = {
    tags: "".split(" "),
    id: "387"
    };
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function() {
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled) {
    StackExchange.using("snippets", function() {
    createEditor();
    });
    }
    else {
    createEditor();
    }
    });

    function createEditor() {
    StackExchange.prepareEditor({
    heartbeatType: 'answer',
    autoActivateHeartbeat: false,
    convertImagesToLinks: false,
    noModals: true,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: null,
    bindNavPrevention: true,
    postfix: "",
    imageUploader: {
    brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
    contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
    allowUrls: true
    },
    onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    });


    }
    });






    Thomas is a new contributor. Be nice, and check out our Code of Conduct.










    draft saved

    draft discarded


















    StackExchange.ready(
    function () {
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f189908%2fcalculating-the-trace-of-the-product-of-several-matrices%23new-answer', 'question_page');
    }
    );

    Post as a guest















    Required, but never shown

























    1 Answer
    1






    active

    oldest

    votes








    1 Answer
    1






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes









    3












    $begingroup$

    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, cutrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR,rank, σ, m, n, eps},
    eps = $MachineEpsilon;
    maxrank = OptionValue["MaxRank"];
    cutrank = OptionValue["CutRank"];
    ϵ = OptionValue["Tolerance"];
    k = 1;
    ik = OptionValue["StartingIndex"];
    vk = row[ik];
    m = Length[vk];
    v = ConstantArray[0., {maxrank, m}];
    jk = LinearAlgebra`BLAS`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 = LinearAlgebra`BLAS`IAMAX[uk];
    vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
    jk = LinearAlgebra`BLAS`IAMAX[vk];
    δ = vk[[jk]];
    If[Abs[δ] < eps,
    w = uk;
    δ = 0.;
    iter = 0;
    While[Abs[δ] < eps,
    iter++;
    w[[ik]] = 0.;
    ik = LinearAlgebra`BLAS`IAMAX[w];
    vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
    jk = LinearAlgebra`BLAS`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.



    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







    share|improve this answer









    $endgroup$


















      3












      $begingroup$

      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, cutrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR,rank, σ, m, n, eps},
      eps = $MachineEpsilon;
      maxrank = OptionValue["MaxRank"];
      cutrank = OptionValue["CutRank"];
      ϵ = OptionValue["Tolerance"];
      k = 1;
      ik = OptionValue["StartingIndex"];
      vk = row[ik];
      m = Length[vk];
      v = ConstantArray[0., {maxrank, m}];
      jk = LinearAlgebra`BLAS`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 = LinearAlgebra`BLAS`IAMAX[uk];
      vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
      jk = LinearAlgebra`BLAS`IAMAX[vk];
      δ = vk[[jk]];
      If[Abs[δ] < eps,
      w = uk;
      δ = 0.;
      iter = 0;
      While[Abs[δ] < eps,
      iter++;
      w[[ik]] = 0.;
      ik = LinearAlgebra`BLAS`IAMAX[w];
      vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
      jk = LinearAlgebra`BLAS`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.



      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







      share|improve this answer









      $endgroup$
















        3












        3








        3





        $begingroup$

        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, cutrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR,rank, σ, m, n, eps},
        eps = $MachineEpsilon;
        maxrank = OptionValue["MaxRank"];
        cutrank = OptionValue["CutRank"];
        ϵ = OptionValue["Tolerance"];
        k = 1;
        ik = OptionValue["StartingIndex"];
        vk = row[ik];
        m = Length[vk];
        v = ConstantArray[0., {maxrank, m}];
        jk = LinearAlgebra`BLAS`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 = LinearAlgebra`BLAS`IAMAX[uk];
        vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
        jk = LinearAlgebra`BLAS`IAMAX[vk];
        δ = vk[[jk]];
        If[Abs[δ] < eps,
        w = uk;
        δ = 0.;
        iter = 0;
        While[Abs[δ] < eps,
        iter++;
        w[[ik]] = 0.;
        ik = LinearAlgebra`BLAS`IAMAX[w];
        vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
        jk = LinearAlgebra`BLAS`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.



        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







        share|improve this answer









        $endgroup$



        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, cutrank, ϵ, u, v, k, ik, jk, uk, vk, test, norm2, δ, w, iter, uQ, uR, vQ, vR,rank, σ, m, n, eps},
        eps = $MachineEpsilon;
        maxrank = OptionValue["MaxRank"];
        cutrank = OptionValue["CutRank"];
        ϵ = OptionValue["Tolerance"];
        k = 1;
        ik = OptionValue["StartingIndex"];
        vk = row[ik];
        m = Length[vk];
        v = ConstantArray[0., {maxrank, m}];
        jk = LinearAlgebra`BLAS`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 = LinearAlgebra`BLAS`IAMAX[uk];
        vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
        jk = LinearAlgebra`BLAS`IAMAX[vk];
        δ = vk[[jk]];
        If[Abs[δ] < eps,
        w = uk;
        δ = 0.;
        iter = 0;
        While[Abs[δ] < eps,
        iter++;
        w[[ik]] = 0.;
        ik = LinearAlgebra`BLAS`IAMAX[w];
        vk = Subtract[row[ik], u[[;; k - 1, ik]].v[[;; k - 1]]];
        jk = LinearAlgebra`BLAS`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.



        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








        share|improve this answer












        share|improve this answer



        share|improve this answer










        answered 4 hours ago









        Henrik SchumacherHenrik Schumacher

        50.9k469145




        50.9k469145






















            Thomas is a new contributor. Be nice, and check out our Code of Conduct.










            draft saved

            draft discarded


















            Thomas is a new contributor. Be nice, and check out our Code of Conduct.













            Thomas is a new contributor. Be nice, and check out our Code of Conduct.












            Thomas is a new contributor. Be nice, and check out our Code of Conduct.
















            Thanks for contributing an answer to Mathematica Stack Exchange!


            • Please be sure to answer the question. Provide details and share your research!

            But avoid



            • Asking for help, clarification, or responding to other answers.

            • Making statements based on opinion; back them up with references or personal experience.


            Use MathJax to format equations. MathJax reference.


            To learn more, see our tips on writing great answers.




            draft saved


            draft discarded














            StackExchange.ready(
            function () {
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f189908%2fcalculating-the-trace-of-the-product-of-several-matrices%23new-answer', 'question_page');
            }
            );

            Post as a guest















            Required, but never shown





















































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown

































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown







            Popular posts from this blog

            How to make a Squid Proxy server?

            Is this a new Fibonacci Identity?

            Touch on Surface Book