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?
A <- matrix(1:4, 2)
sum(diag(A %*% A))
[1] 29
sum(A*A)
[1] 30
sum(diag(A %*% t(A)))
[1] 30
and sum(A*t(A))
[1] 29
This is the answer you want, I think. [Sorry: I hit send before I
thought through what I was doing. Suppose A is a row vector and B a
column vector. Then the elementwise product A*B returns an error in R
but A%*%B is a scalar while B%*%A is not, and sum(A*t(B)) = sum(t(A)*B)
= trace(A%*%B).
When A and B are not square but the product is, then A*B = the
elementwise product of A and B is not defined, but A*t(B) is, and so is
t(A)*B. When A
See "Trace of a product" in the Wikipedia article on "Trace
(linear algebra)"
(https://en.wikipedia.org/wiki/Trace_%28linear_algebra%29#Trace_of_a_product).
Spencer
Thanks
--
Wagner Hugo Bonat
LEG - Laboratório de Estatística e Geoinformação
UFPR - Universidade Federal do Paraná
[[alternative HTML version deleted]]
______________________________________________
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.
______________________________________________
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.
