msmbuilder / msmbuilder-legacy

Legacy release of MSMBuilder
http://msmbuilder.org
GNU General Public License v2.0
25 stars 28 forks source link

Force Dense Eigen for small matrices #406

Closed kyleabeauchamp closed 10 years ago

kyleabeauchamp commented 10 years ago

So with the addition of get_reversible_eigenvectors, we have a bug where small MSMs cause errors in CalculateImpliedTimescales:

ValueError: k must be between 1 and rank(A)-1
14:27:55 - BFGS likelihood maximization terminated after 7 function calls.  Initial and final log likelihoods: -36293.225243, -36293.218595.
14:27:55 - ARPACK cannot calculate 2 eigenvectors from a 2 x 2 matrix.
14:27:55 - Instead, calculating 0 eigenvectors.

The issue is that the get_eigenvectors explicitly converts small matrices to dense mode, while get_reversible_eigenvectors skips this check and always uses sparse eigen:

# only happens in get_eigenvectors()
    if n < dense_cutoff and scipy.sparse.issparse(t_matrix):
        t_matrix = t_matrix.toarray()
kyleabeauchamp commented 10 years ago

I think the best way to solve this is to just copy back in the dense checker from get_eigenvectors

kyleabeauchamp commented 10 years ago

There's actually some other minor issues in get_reversible_eigenvectors, such as this line:

elif k >= n - 1  and scipy.sparse.issparse(t_matrix):

We check for issparse(), but we don't have any cases for dense.

kyleabeauchamp commented 10 years ago

Also, ARPACK can't calculate 1 eigenvector from a 2x2 matrix:

In [15]: scipy.sparse.linalg.eigs(np.eye(2), k=1)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-04b219afc71c> in <module>()
----> 1 scipy.sparse.linalg.eigs(np.eye(2), k=1)

/home/kyleb/opt/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in eigs(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, OPpart)
   1269     params = _UnsymmetricArpackParams(n, k, A.dtype.char, matvec, mode,
   1270                                       M_matvec, Minv_matvec, sigma,
-> 1271                                       ncv, v0, maxiter, which, tol)
   1272 
   1273     while not params.converged:

/home/kyleb/opt/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in __init__(self, n, k, tp, matvec, mode, M_matvec, Minv_matvec, sigma, ncv, v0, maxiter, which, tol)
    683                              % ' '.join(_NEUPD_WHICH))
    684         if k >= n - 1:
--> 685             raise ValueError("k must be less than rank(A)-1, k=%d" % k)
    686 
    687         _ArpackParams.__init__(self, n, k, tp, mode, sigma,

ValueError: k must be less than rank(A)-1, k=1
kyleabeauchamp commented 10 years ago

See #407

kyleabeauchamp commented 10 years ago

I feel like some large fraction of my life has been spent working around ARPACK strangeness...