>>>>> Wagner Bonat <wbo...@gmail.com>
>>>>> on Mon, 25 Aug 2014 10:31:41 +0200 writes:
> I need to compute two equations related with trace and inverse of a around
> 30000 x 30000 density matrices. The equations are
> -trace( W_i %**% C) and -trace(W_i %**% C %*% W_j C)
[I assume that 2nd eq is missing a %*% ]
> I know W_i, W_j and inverse of C. These equations are related with Pearson
> estimating functions. I am trying to use R and package Matrix, but I
> couldn't compute the C matrix, using solve() or chol() and chol2inv(). I
do
> not know with is possible using solve() to solve a system of equation and
> after compute the trace. It is common to use solve function to compute
> something like C^{-1} W = solve(C, W), but my equation is a little bit
> different.
Not too much different, fortunately for you.
First, note that, mathematically,
tr(A %*% B) == tr(B %*% A)
whenever both matrix products are valid, e.g. when the matrices
are square of the same dimension.
Consequently, you typically can interchange the order of the
matrices in the product __ when inside the trace __
As you know D := C^{-1} and really need C = D^{-1}, let's
better use D notation, so you want
t1 <- -trace(W_i %**% D^{-1})
t2 <- -trace(W_i %**% D^{-1} %*% W_j %*% D^{-1})
so, if
CWi <- solve(D, W_i) {for 'i' and 'j' !}
t1 <- -trace(CWi)
t2 <- -trace(CWi %*% CWj)
Now, this alone will still not be good enough to get done
quickly, using Matrix:
The most important question really is if D ( = C^{-1} ) is
a *sparse* matrix and possibly the W_j are sparse as well.
In some (but not most) such cases, C will be sparse, too, and
the whole computations are done very efficiently using the
Matrix - underlying C libraries.
I'm interested to hear more about your matrices.
To find their sparse, apply
nnzeros( M )
and possibly
nnzeros(zapsmall( M ))
to your matrices M.
Also of interest here is
isSymmetric(D)
Martin Maechler,
ETH Zurich
and Maintainer of the 'Matrix' package
> Any help is welcome. I am thinking about to use RcppArmadillo, but I am
not sure that it is able to compute my equations.
> Thank you everyone.
> --
> Wagner Hugo Bonat
> LEG - Laboratório de Estatística e Geoinformação
> UFPR - Universidade Federal do Paraná
______________________________________________
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.
