[R] Trace of product of matrices

[Martin Maechler]

>>>>> peter dalgaard <pdalgd at gmail.com>
>>>>>     on Sun, 19 Oct 2014 21:26:39 +0200 writes:

    >> On 19 Oct 2014, at 19:00 , Spencer Graves <spencer.graves at structuremonitoring.com> wrote:
    >> 
    >> On 10/19/2014 8:42 AM, peter dalgaard wrote:
    >>>> On 19 Oct 2014, at 16:43 , Wagner Bonat <wbonat at gmail.com> wrote:
    >>>> 
    >>>> Dear,
    >>>> 
    >>>> I have to compute the trace of a product between four matrices. For
    >>>> example, I know the matrices Wi, Wj and C, I need to compute this
    >>>> 
    >>>> -trace(Wi%*%C^-1%*%Wj%*%C^-1)
    >>>> 
    >>>> 
    >>>> I would like to avoid compute the complete matrix and after take the
    >>>> diagonal, something like
    >>>> 
    >>>> sum(diag( solve(Wi,C)%*% solve(Wj,C)))
    >>> <this can't be right: it is C that is the invertible matrix>
    >>> 
    >>>> Any idea is welcome.
    >>>> 
    >>> The usual "trick" is that the trace of a matrix product is the inner product in matrix space, which is just the sum of the  elementwise products
    >>> 
    >>> tr(AB) = tr(BA) = sum_i sum_j a_ij b_ij.
    >>> 
    >>> In R, this becomes simply sum(A*B) -- notice that the ordinary product is used, not %*%. So presumably, you are looking for
    >>> 
    >>> sum(solve(C, Wi) * solve(C, Wj))
    >> 
    >> missing a transpose?

    > Yep... 

    > tr(AB) = tr(BA) = sum_i sum_j a_ij b_ji

    > which is of course sum(A*t(B)) or vice versa. Thanks.

Thank you, Peter, and Spencer.

For a few years now, I have had in my TODO file for the Matrix
package:

** TODO tr(A %*% B) {and even  tr(A %*% B %*% C) ...} are also needed
  frequently in some computations {conditional normal distr. ...}.
  Since this can be done faster than by
    sum(diag(A %*% B))  even for traditional matrices, e.g.
    	       sum(A * t(B)) or {even faster for "full" mat}
	       crossprod(as.vector(A), as.vector(B))
  and even more so for, e.g.  <sparse> %*% <dense>
  {used in Soeren's 'gR' computations},
  we should also provide a generic and methods.

** TODO diag(A %*% B) might look like a "generalization" of tr(A %*% B),
  but as the above tricks show, is not really.
  Still, it's well worth to provide  diag.prod(A, B):

  Well, if A %*% B is square,   diag(A %*% B)  ===  colSums(t(A) * B)
  and we should probably teach people about that !

-----------

Are there good suggestions for a sensible function name for
these potential matrix utility function?

  trprod()
  traceprod()

  diagprod()
?

--
Martin <Maechler at stat.math.ethz.ch>  http://stat.ethz.ch/people/maechler
Seminar für Statistik, ETH Zürich  HG G 16      Rämistrasse 101
CH-8092 Zurich, SWITZERLAND
phone: +41-44-632-3408       fax: ...-1228      <><