# [R] Trace and inverse of big matrices

Martin Maechler maechler at stat.math.ethz.ch
Mon Aug 25 18:29:40 CEST 2014

```>>>>> 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.

To find their sparse, apply
nnzeros( M )
and possibly
nnzeros(zapsmall( 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á

```