[R] Trace and inverse of big matrices

[Martin Maechler]

>>>>> Wagner Bonat <wbonat at 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á