I have a problem in hand that involves calculation of eigenvector of a large sparse matrix ($10^7 \times 10^7$). A first-order Markov chain is fit to the problem, and the state transition probability matrix has only 7 non-zero entries at each column. I used power method and Rayleigh quotient method to quickly calculate the steady state occupancy distribution, and finally MATLAB's eigs function with the option to return only the largest eigen pair.

Apart from which method is more efficient and all that, even the construction of the state transition probability matrix takes some time. Numerical analysis is not my field, but I was wondering if the the following facts can help me to improve the calculation speed:

  • I'm only interested in eigenvector corresponding to the largest eigenvalue (which is obviously 1)
  • the non-zero entries of the matrix can be calculated on the fly, I mean there is routine $f$ which $p_{ij} = f(i, j)$, $\text{P}=[p_{ij}]$ is the state transition probability matrix.

Any help/reference is much appreciated.

  • $\begingroup$ Do you have an idea of the order of magnitude of $1-\lambda_2$, where $\lambda_2$ is the second largest eigenvalue? Also, it may be possible to try to find the stationary distribution by actually solving the system $\pi P=P,\sum_i \pi_i=1$ using a sparse linear system routine instead of an eigenvalue routine. My concern with that is that you are solving $n+1$ equations in $n$ unknowns, which is a rank deficient problem, which can sometimes cause problems. The easy fixes (e.g. taking an SVD) destroy sparsity. Also, is your transition matrix symmetric or not? $\endgroup$ – Ian Sep 28 '16 at 0:43
  • $\begingroup$ @ByronSchmuland, thanks. I fixed that. $\endgroup$ – Ali Sep 28 '16 at 0:45
  • $\begingroup$ @Ian, no unfortunately I don't have any other side information, and the matrix is not symmetric. $\endgroup$ – Ali Sep 28 '16 at 0:46
  • $\begingroup$ You can solve $n$ equations in $n$ unknowns: $\sum_i \pi_i = 1$ and $n-1$ of the $n$ entries of $\pi P = \pi$; since the row sums are all $1$, the omitted equation is a consequence of these. $\endgroup$ – Robert Israel Sep 28 '16 at 0:59
  • $\begingroup$ @RobertIsrael, thanks for the comment. I actually tried this but didn't help me much as I have to construct P matrix again. As I mentioned this also takes quite some time, the routine $f$ in my post is not a simple on-line equation and has to be called $7 \times 10^7$ to construct the matrix :( $\endgroup$ – Ali Sep 28 '16 at 1:05

In this case the Arnoldi iteration method for calculating eigenvalues (i.e. the method used by eigs function) will probably be the fastest method to find stationary distribution.

You need to have an efficient method for forming a sparse-dense multiplication $Av$ many times for different $v$. If the matrix $A$ does not have any special structure, that can speed up forming $Av$, then you need to costruct this matrix explicitly. Calculating elements of $A$ on the fly will only slow down computations considerably, since all elements of $A$ will be recalculated many times.

  • $\begingroup$ As a remark, eigs can take a function handle representation of the matrix rather than the matrix itself. This might be faster than actually assembling A, though I'm a bit skeptical of that; if we naively assume that storing $7 \times 10^7$ entries of a sparse matrix takes $21 \times 10^7$ doubles, then that's about $1.7 \times 10^9$ bytes, which is large but fits in RAM. And I think the actual sparse matrix implementation is a bit more clever than that. By comparison, dynamically recomputing those entries using a function handle might very well be quite a bit slower. $\endgroup$ – Ian Sep 28 '16 at 2:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.