[R] Trace and inverse of big matrices

[Wagner Bonat]

Thank you Martin for your advices.

Let me give more information about my matrices. I am working on a
space-time model, it is a very simple model based on second-moment
assumption. My model is

E(Y) = g(X%*%beta) -> g is a suitable link function
V(Y) = C = V^{1/2} Omega V^{1/2} -> It is my variance-covariance matrix


It is my general formulation, for this specific case I am trying to use
similar structure used in Bayesian inference, based on Gaussian Markov
Random Fields, in general these models are specified by their inverse or
precision matrix, and all these matrices are sparse. Specifically  I am
using a CAR (Conditional autoregressive structure) for spatial effects and
RW1( first order random walk) for time effects. Then, my
variance-covariance matrix perhaps I should call dispersion matrix is the
form,

C^{-1} = V^{-1/2} Omega^{-1} V^{-1/2} where V = diag(E[Y]^power) is a
diagonal matrix with the variance function. In this case power variance
function, because I am fitting a Tweedie model. I assume that the
Omega^{-1} can be written as a linear combination of known matrices,

Omega^{-1} = tau0*I + tau1*R_t + tau2*R_s + tau3*R_ts

where R_t, R_s and R_ts are suitable matrices to describe the dependence
structure of the time, space and interaction (space*time), all matrices are
sparse. I use these functions to build these matrices (see below). The full
matrix is given by

Omega^{-1} = tau0* Diagonal(n.obs, 1) + tau1*kronecker(Diagonal(n.time,1),
R_s) + tau2*kronecker(R_t, Diagonal(n.space,1) + tau3*kronecker(R_t, R_s)

This is my space-time model with interaction term. To make inference I am
using the quasi-score function for regression (similar GEE approach) and
Pearson estimating fnuction for covariance parameters, in this case power,
tau0, tau1,tau2 and tau3. The Pearson estimating functions are

phi(p, tau0,tau1,tau2,tau3)_p = t(res)%*%W_p%*%res - tr(W_p%*%C) -> It is
for the power parameter
phi(p, tau0,tau1,tau2,tau3) = t(res)%*%W_tau0%*%res - tr(W_tau0%*%C) -> It
is for the tau0 parameter
phi(p, tau0,tau1,tau2,tau3) = t(res)%*%W_tau1%*%res - tr(W_tau1%*%C) -> It
is for the tau1 parameter

Similar for tau2 and tau3. The W matrices are weights matrices, defined by
the negative the derivative of C matrix with respect every parameter, in
this case very easy for example for tau0 is

W_tau0 = V^{-1/2} I V{-1/2}
W_tau1 = V^{-1/2} kronecker(Diagonal(n.time,1), R_s) V{-1/2} and similar
for tau2, tau3.

The weights matrices are simple and sparse, my problem is to compute the
trace tr(W%*%C) where I need to solve of the C^{-1} or as you called D
matrix.

Thank you for your time !


## Space matrix
mat.espacial <- function(list.viz,ncol,nrow){
vizi <- Matrix(0,nrow=nrow,ncol=ncol, sparse = TRUE)
diagonal <- c()
for(i in 1:ncol){
  diagonal[i] <- length(list.viz[[i]])}
diag(vizi) <- diagonal
for(i in 1:ncol){
vizi[i,c(list.viz[[i]])] <- -1}
return(vizi)}

## Time matrix
mat.temporal <- function(nrow,ncol){
R <- Matrix(0,nrow=nrow,ncol=ncol, sparse = TRUE)
## Restrições de borda
R[1,c(1,2)] <- c(1,-1)
R[ncol,c(ncol-1,ncol)] <- c(-1,1)
## Corpo da matriz
n <- ncol-1
for(i in 2:n){
R[i,c(i-1,i,i+1)] <- c(-1,2,-1)}
R <- as(R, "symmetricMatrix")
return(R)}







2014-08-25 18:29 GMT+02:00 Martin Maechler <maechler at stat.math.ethz.ch>:

> >>>>> 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á
>



-- 
Wagner Hugo Bonat
LEG - Laboratório de Estatística e Geoinformação
UFPR - Universidade Federal do Paraná

	[[alternative HTML version deleted]]