The nice thing about the power method is its conceptual simplicity. In cases like this, it's quite hard to screw it up. And it will scale far beyond those numpy functions (not that this is needed for this example.)
Also, did anyone mention PageRank yet?
m /= m.sum(axis=0)[numpy.newaxis,:]
u = numpy.ones(len(items))
for i in xrange(100):
u = numpy.dot(m, u)
u /= u.sum()
Also, while scalability for this size computation is way overkill, it is precisely a problem where the numerical stability problem gets even worse. Imagine if I have data on some kind of power law, and compute left-to-right `M+epsilon_0+...+epsilon_n`. No matter now large `n` is, for sufficiently different order of magnitude M and epsilon all the information is lost -- `epsilon_0+...+epsilon_n+M` could be an entirely different number. Highly recommend checking out the LAPACK stability guide here http://www.netlib.org/lapack/lug/node72.html for more.
It's actually the second highest eigenvalue. The highest eigenvalue is always 1 for stochastic matrices.
All this being said, scalability is _obviously_ a non-issue when talking about a matrix of programming languages. All methods are constant time.
[1]: https://en.wikipedia.org/wiki/Matrix_multiplication_algorith... Interesting tidbit: nobody can even prove it's not N^2 :)
And yea, mat mul is not N^3 theoretically, but most implementations are. I've heard that some (mkl maybe) are 2.8, but haven't had someone point code to me. My personal attempts at implementing Strassen were slower than a tuned N^3 implementation, at least for matrices that fit into memory.