[R] eigenvalues for a sparse matrix

[Douglas Bates]

Tamas Papp <tpapp at axelero.hu> writes:

> I have the following problem.  It has two parts.
>
> 1. I need to calculate the stationary probabilities of a Markov chain,
> eg if the transition matrix is P, I need x such that
>
> xP = x
>
> in other words, the left eigenvectors of P which have an eigenvalue of
> one.
>
> Currently I am using eigen(t(P)) and then pick out the vectors I need.
> However, this seems to be an overkill (I only need a single vector!)
> and takes a lot of time -- P is 1176 x 1176!  Is there a faster way?
>
>
> 2. In fact, P has a structure: it comes from the solution of a
> discrete dynamic optimzation problem.  There are exogenous (X) and
> endogenous (N) states, and I have a policy function X x N -> N, which
> gives the choice of the agent for any (x,n) in (X,N).  X has an
> exogenous transition matrix.  I use the following function to build
> the "global" transition matrix:
>
> globalTransition <- function(U, modelenv) {
>   G <- matrix(0, modelenv$Nn*modelenv$Xn, modelenv$Nn*modelenv$Xn)
>                                         # this matrix is very sparse
>   for (i in 1:modelenv$Nn) {
>     r <- (i - 1)*modelenv$Xn            # start at this row
>     for (j in 1:modelenv$Xn) {
>       G[r+j, 1:modelenv$Xn+modelenv$Xn*(U[j,i]-1)] <-
>         modelenv$Xtrans[j, 1:modelenv$Xn]
>     }
>   }
>   G
> }
>
> Nn and Xn are the number of engo- and exogenous states, U is the said
> policy function in a matrix form, and Xtrans is the transition matrix
> of exogenous states.
>
> The row (and column) indexes of G run like this:
>
> index Nn Xn
> 1     1  1
> 2     1  2
> 3     1  3
> ...
> Xn    1  Xn
> Xn+1  2  1
> ...
> Nn*Xn Nn Xn
>
> As you can see from the above function, each row contains just a few
> nonzero elements in columns Xn*(U[j,i]-1)+1, ...; this means that the
> agent chose the endogenous state n = U[j,i].

You should be able to generate the triplet form of a sparse matrix
from that representation then generate the compressed sparse
column-oriented form using facilities in the new Matrix package for
R-1.9.0 

I'm not sure about getting eigenvectors for specific eigenvalues from
that representation, although Simon Wood's suggestion of using inverse
iteration seems like it will be quite useful.  The point of using the
cscMatrix class would be that the iteration itself would be much
faster than the same iteration using the dense matrix representation.

Look for the classes "tripletMatrix" and "cscMatrix" in Matrix version
0.8-1 or later.

There is a method for crossprod(mm, x) where mm is a cscMatrix and x
is numeric.  Because mm is column oriented the greatest efficiency is
achieved by multiplying by t(mm) on the left or mm on the right.