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 :
> >>>>> Wagner Bonat
> >>>>> on Mon, 25 Aug 2014 10:31:41 +0200 writes:
>
> > I need to compute two equations related with trace and inverse of a
> around
> > 3 x 3 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' !}
>
