R help - Apr 2004 - eigenvalues for a sparse matrix

Hi,

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

Thanks,

Tamas

-- 
Tam??s K. Papp
E-mail: tpapp at axelero.hu
Please try to send only (latin-2) plain text, not HTML or other garbage.

If you know the eigenvalues then you can use inverse iteration to get the
eigenvectors you need: I think that this only requires you to
be able to multiply a vector by your matrix, so that you can take
advantage of the sparse structure, but I'm away from home and can't
check
references. However Golub and van Loan (1996) will probably tell you what
to do.

Simon

_____________________________________________________________________
> Simon Wood simon at stats.gla.ac.uk        www.stats.gla.ac.uk/~simon/
>>  Department of Statistics, University of Glasgow, Glasgow, G12 8QQ
>>>   Direct telephone: (0)141 330 4530          Fax: (0)141 330 4814

>
> 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].
>
> Thanks,
>
> Tamas
>
> --
> Tam??s K. Papp
> E-mail: tpapp at axelero.hu
> Please try to send only (latin-2) plain text, not HTML or other garbage.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>