Numerical computation of Perron-Frobenius eigenvector

Question

I would like to efficiently and (to the extent possible) reliably find the Perron-Frobenius eigenvector of non-negative matrices. These are not stochastic matrices, they are typically dense, and their entries differ by many orders of magnitude. They are probably usually irreducible, though if you think of the very small entries as zero they will not be. They are not symmetric. They will be as large as I can make them with reasonable computation time on a laptop.

I'm not much of a numerics person and up to now I've tended to treat eigenvector routines as a "back box". However, in this case none of the off-the-shelf solutions I've tried work very well. Numpy's linalg.eig is very slow for matrices bigger than $2000\times 2000$ or so, presumably because it always computes every eigenvector, not just the top one. Moreover, the eigenvalue with the largest real part often does not correspond to an eigenvector with all non-negative (or even real) entries, which messes everything up. On the other hand, scipy.sparse.linalg.eigs is great when it works (despite my matrices being dense), but sometimes I just can't get it to converge no matter what I do.

I tried implementing power iteration by myself, but it's very slow and also suffers from convergence issues - it leaves me feeling that there ought really to be a better way.

Ideally I would like an algorithm with the following properties:

• faster than power iteration and with better convergence properties
• able to give me some kind of approximate solution in the case where it can't converge
• guaranteed not to give me any eigenvector other than the Perron-Frobenius one, i.e. not to give me something with complex elements or a mixture of positive and negative ones
• able to compute the left and right Perron-Frobenius eigenvectors simultaneously, if that would be more efficient than simply running it twice
• either available off-the-shelf (in Python or some other language), or explained in a form such that I, as a non-numerics person, can understand how to implement it without worrying too much about the theory behind it

Does such a thing exist? References to literature giving an overview of how to solve this kind of problem would also be welcome, as I'm sure I'm not the first person to need to do this.

Answer 1

The algorithm used by eigs is called Arnoldi iteration (in shift-and-invert mode). You could think of it as a sophisticated way of doing an inverse power iteration.

The inverse power iteration method works as follows: if $A$ is your matrix and you know that the Perron root (the spectral radius of $A$) is $\lambda$, then a power iteration on the matrix $B = (A - \lambda I)^{-1}$, where $I$ is the identity matrix, will converge to the Perron-Frobenius eigenvector of $A$, let's call it $v$, much faster than a simple power iteration on $A$. This is so because $v$ is an eigenvector of the largest eigenvalue of $B$, which will be close to infinity by construction and hence much larger than the next largest eigenvalue.

You still have to obtain an estimate of $\lambda$ but you could use the Gershgorin circle theorem for this purposeor use other techniques (such as the bounds given in [Duan2013]).

Therefore, my suggestion is to use something like:

v = np.ones(A.shape[0])
eigenvalues, eigenvectors = linalg.eigs(A, k=2, sigma=shift, which='LM', v0=v)

where the variable shift equals your estimate for $\lambda$. Then check that the resulting eigenvalues are sufficiently far apart and pick the eigenvector corresponding to the Perron root.

Note that if the matrix is not irreducible, you will get more than one eigenvector corresponding to the largest eigenvalue. Also note that even though $B$ involves a matrix inverse, it is rarely ever necessary (or convenient) to actually compute the inverse.

[Duan2013] Duan, X., & Zhou, B. (2013). Sharp bounds on the spectral radius of a nonnegative matrix. Linear Algebra and Its Applications, 439(10), 2961–2970. http://doi.org/10.1016/j.laa.2013.08.026