2
$\begingroup$

According to the Perron-Frobenius theorem, a real matrix with only positive entries (or one with non-negative entries with a property called irreducibility) will have a unique eigenvector that contains only positive entries. Its corresponding eigenvalue will be real and positive, and will be the eigenvalue with greatest magnitude.

I have a situation where I'm interested in such an eigenvector. I'm currently using numpy to find all the eigenvalues, then taking the eigenvector corresponding to the one with largest magnitude. The trouble is that for my problem, when the size of the matrix gets large, the results start to go crazy, e.g. the eigenvector found that way might not have all positive entries. I guess this is due to rounding errors.

Because of this, I'm wondering if there's an algorithm that can give better results by making use of the facts that $(i)$ the matrix has non-negative entries and is irreducible, and $(ii)$ we're only looking for the eigenvector whose entries are positive. Since there are algorithms that can make use of other matrix properties (e.g. symmetry), it seems reasonable to think this might be possible.

While writing this question it occurred to me that just iterating $\nu_{t+1} = \frac{A\nu_t}{|A\nu_t|}$ will work (starting with an initial $\nu_0$ with positive entries), but imagine with a large matrix the convergence will be very slow, so I guess I'm looking for a more efficient algorithm than this. (I'll try it though!)

Of course, if the algorithm is easy to implement and/or has been implemented in a form that can easily be called from Python, that's a huge bonus.

Incidentally, in case it makes any kind of difference, my problem is this one. I'm finding that as I increase the matrix size (finding the eigenvector using Numpy as described above) it looks like it's converging, but then suddenly starts to jump all over the place. This instability gets worse the smaller the value of $\lambda$.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

The algorithm you describe that computes $x^{(k+1)}=\frac{Ax^{(k)}}{\|Ax^{(k)}\|}$ is of course what is called the the power method. It will converge in your case if you have a non-degenerate largest eigenvalue. Furthermore, if you start with an initial guess $x^{(0)}$ that has only positive entries, you are guaranteed that all future iterates are also strictly positive and, moreover, that round-off should not be an issue since in every matrix-vector product you only ever add up positive terms. In other words, this method must work if your problem has the properties you describe.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Won't the convergence be quite slow, though? My matrices are potentially huge, since I'm trying to approximate the limit of infinite matrix size. (But then, maybe it's not slow, I haven't tried yet. I'll do it this afternoon if I get a chance.) $\endgroup$
    – N. Virgo
    Feb 10, 2014 at 4:47
  • $\begingroup$ The convergence rate depends on the ratio $\lambda_1/\lambda_2$, where $\lambda_1$ is the largest eigenvalue and $\lambda_2$ the next largest, provided that the eigenspace corresponding to $\lambda_1$ is 1-dimensional. For some families of matrices, the ratio may stay the same even as the matrix size increases. $\endgroup$
    – Daniel Shapero
    Feb 10, 2014 at 5:31
  • $\begingroup$ It turns out that it isn't all that slow for my problem - but I found a way to make it much faster anyway, which I've posted as an answer. $\endgroup$
    – N. Virgo
    Feb 10, 2014 at 11:50
1
$\begingroup$

I tried the power method, as described in my question and in Wolfgang Bangerth's answer. It turns out that the convergence isn't unacceptably slow for my application, though it does depend on the parameters. So in a way this was kind of a dumb question.

However, I noticed that there's a way to exponentially increase the speed of this algorithm, by doing repeated squaring of the matrix. That is, let $B_0=A$ and iterate $B_{n+1}=\frac{B_n^2}{\|B_n^2\|}$, where $\|\cdot\|$ is whatever matrix norm you feel like. (I just summed all the elements, since they're all positive anyway.) This very rapidly converges to a matrix whose columns are each proportional to the leading eigenvector. This is because $B_n \nu_0 \propto A^{2^n}\nu_0$, so iterating this algorithm for $n$ steps is the same as doing power iteration for $2^n$ steps.

Although multiplying large matrices can be slow (especially in numpy, unfortunately), I'm finding that this tends to converge pretty nicely after around 10 to 15 iterations, so on the whole it's pretty fast.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Alternatively, you could accelerate the method by using inverse power iteration. In your case, it amounts to solving the system A - I (or A - 0.999 I) at each iteration. $\endgroup$
    – Juan M. Bello-Rivas
    Feb 10, 2014 at 14:53
  • 1
    $\begingroup$ Are you multiplying your matrices explicitly? I can not imaging a scenario where it is faster that matrix-vector multiplications. The algorithm you describe might converge in fewer iterations, but each iteration involves to calculate $B_n \times B_n$, and it should be much more expensive than $B_n v$. Maybe I am missing something... $\endgroup$
    – sebas
    Feb 10, 2014 at 15:28
  • $\begingroup$ @sebas the choice is between 15 matrix multiplications, or $2^{15}$ matrix-vector ones. I guess for large matrices the vector version will win out, but for smallish ones this is way faster. $\endgroup$
    – N. Virgo
    Feb 10, 2014 at 15:54
  • $\begingroup$ @Nathaniel It depends on the size of the matrix. I was just wondering how large are the matrices for your problem. I know that for large enough matrices just 2 or 3 matrix-matrix multiplications is completely prohibitive. Moreover if your matrices are sparse, because the matrix-matrix product will increase the fill-in of the resulting matrix. Just a comment I thought can be useful. $\endgroup$
    – sebas
    Feb 10, 2014 at 16:02
  • $\begingroup$ @sebas honestly I think you're right, my idea isn't worthwhile once the matrices are above 500x500 or so. But for exploring the parameter space with a size of around 100 it's quite handy. $\endgroup$
    – N. Virgo
    Feb 10, 2014 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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