How to calculate the eigenvector of a column stochastic matrix in C++

I have a column stochastic matrix A(n-by-n real, non-negative matrix) and want to solve the following equation in C++: Ax=x

I am assuming that I need to find out an eigenvector, x, where eigenvalue wou开发者_如何学Gold have to be 1(right?) but I couldn't figure it out in C++. So far I have checked out some math-libs such as Seldon, CPPScaLapack, Eigen... Among them, Eigen seems a good option but I couldn't understand how to utilize any of them to solve the equation above.

Can you give me some suggestions/code snippets or ideas to solve the equation? Any help is highly appreciated.

Thanks.

Since the largest eigenvector of a stochastic matrix $M$ is unity, you can find this eigenvector by iteration (unless you're really bad at guessing initial values).

Start with some randomly chosen initial vector $v_1$ whose values (probabilities) sum to unity. Apply $M$ to $v_1$ to get $Mv_1$. Now renormalize this new vector $Mv_1$, that is, divide by the sum of its elements to get $v_2$. This is a new vector of probabilities and it will be closer to the desired eigenvector (unless your initial guess happened to be orthogonal to the eigenvector).

Repeat this process until $v_k$ approaches stability. That should be a vector with $Mv_k = v_k$, as desired.

Another method is to compute the kernel of your matrix minus the identity matrix. This may or may not be faster than using power iteration as explained by Carl, depending on the size of the matrix and the other eigenvalues. Power iteration is better when the matrix gets bigger and the second eigenvalue gets further away from one.

The idea is to rewrite Ax = x as Ax - x = 0. Then use that Ix = x, where I denotes the identity matrix. Thus, Ax - x = 0 is equivalent to Ax - Ix = 0 or (A-I) x = 0. So your eigenvector x lies in the kernel (or null space) of A-I.

This tutorial page explains how to compute the kernel of a matrix using Eigen. Some untested code:

MatrixXf M;
/* initialize M */
FullPivLU<MatrixXf> lu_decomp(M);
VectorXf x = lu_decomp.kernel().col(0);
/* x now contains the vector you want */

You may find that the kernel is empty. This means that either the matrix is not really stochastic, or that you need to adapt the threshold (see the page linked above).