That's not really the definition of numerical stability. You're thinking of linear system stability, which is a whole other topic. Numerical stability is how resilient a computation is to computational error. These computations are taking place using floating point values, the issues arise in the error terms for floats. The author's algorithm is reproduced below:

    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()

Immediately concerning is the use of a dot product and a sum, which will lose a lot of information when the values being added are of different orders of magnitude. (For example `M + epsilon = M` in many floating point computations).

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.

Me: https://twitter.com/tmarkhor/status/842347955390685184

Why do you say it will scale far beyond? Mat mul is N^3 as is eigenvalue solving.

It's actually the second highest eigenvalue. The highest eigenvalue is always 1 for stochastic matrices.

Power method is not matrix-matrix multiplication (which is not N^3, BTW [1]), but rather matrix-vector multiplication. So the power method is N^2*k where k is the number of iterations required to reach precision (usually polylogarithmic).

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 :)

Oh right, of course, because you're iterating the distribution. Duh.

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.