4 ms·

I appreciate that the author wanted to implement their own eigen vector/value method, but really they should use: numpy.linalg.eig(x)[0] numpy.linalg.e

• mathgenius 10y ago
  Yes, but... The matrix has all non-negative entries, and the author is after the highest eigenvalue/vector so I think this means stability is just not an issue. The only possible issue is time to convergence. 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?
  • vmarkovtsev 10y ago
    Me: https://twitter.com/tmarkhor/status/842347955390685184 https://twitter.com/tmarkhor/status/842347955390685184
  • howeman 10y ago
    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.
    • j2kun 10y ago
      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_algorithm#Sub-cubic_algorithms https://en.wikipedia.org/wiki/Matrix_multiplication_algorith... Interesting tidbit: nobody can even prove it's not N^2 :)
      • howeman 10y ago
        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.
  • jlrubin 10y ago
    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 http://www.netlib.org/lapack/lug/node72.html for more.
• emcq 10y ago
  They are just describing power iteration which is a standard technique for finding the primary eigenvector. It is commonly used for web scale recommendation problems which likely explains it's usage given the author's prior background. Also underlies pagerank. AFAIK it is quite stable but converges slowly. From my experience within 50 iterations you will have converged to a stable result.