>>>>> peter dalgaard <pda...@gmail.com>
>>>>> on Sun, 19 Oct 2014 21:26:39 +0200 writes:
>> On 19 Oct 2014, at 19:00 , Spencer Graves
<spencer.gra...@structuremonitoring.com> wrote:
>>
>> On 10/19/2014 8:42 AM, peter dalgaard wrote:
>>>> On 19 Oct 2014, at 16:43 , Wagner Bonat <wbo...@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 <maech...@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 <><
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
