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