Lol after debugging and fixing this and then reading #10715 it's obvious that @sky88088 diagnosed the problem very accurately in his initial post.

You may also want to completely remove the sign flip in the Laplacian

For speed and simplicity, you may want just to call a dense solver if n_nodes < 5 * n_components no matter what the other conditions are.

The conditions are really complicated to follow. Recalling and trusting my old self (probably I should not do that) #10715 (comment), we could simplify with something like:

if not sparse.isparse(laplacian) and eigen_solver == 'amg':
    # warns that we switch to arpack
    eigen_solver = 'arpack'

if eigen_solver == 'arpack':
    # solve the problem
    try:
        ...
    except RuntimeError:
        eigen_solver = 'lobcpg'

# from that point eigen_solver will be either amg or lobcpg
# check if we have enough n_nodes
if n_nodes < 5 * n_components:
    # use eigh
else:
    # check that we should precondition
    if eigen_solver == 'amg':
        # preconditionne and get M
    # call lobcpg with M if available otherwise go for default

Anyway, we should probably merge the regression test and later refactor the code to be readable.

@lobpcg dense solver meaning arpack?

@glemaitre your code has the same issue I'm fixing here. There is no return anywhere, so "# from that point eigen_solver will be either amg or lobcpg" is not true?

Ok I also really don't understand the sign flip tbh.

@lobpcg dense solver meaning arpack?

dense solver meaning scipy.linalg.eigh or numpy.linalg.eigh (unsure about the difference), also turning a sparse Laplacian into dense if needed. The point being that if n_nodes < 5 * n_components it it likely cheaper to compute all eigenvectors and drop some.

@lobpcg? Shouldn't the fallback for arpack just be the standard eigsh, not lobpcg? The benefits of shift-invert mode are not within my expertise, can you comment on the benefits of using that instead of using the standard eigsh?

Ok I also really don't understand the sign flip tbh.

laplacian *= -1 in several places makes no sense whatsoever to me.
There used to be a bug in lobpcg, now fixed, for small matrices requiring this, but not any more.

@lobpcg did you see the comment above the arpack magic? That explains what's happening. I don't think it's related to lobpcg.

There's a bunch of other issues here, but I think this PR fixes one obvious bug and adds one clear regression test and we should merge it. Then we can triage the other bugs.

@glemaitre would love to increase coverage but calling lobpcg actually fails, see #13393 (comment)

@lobpcg? Shouldn't the fallback for arpack just be the standard eigsh, not lobpcg? The benefits of shift-invert mode are not within my expertise, can you comment on the benefits of using that instead of using the standard eigsh?

For "small" matrices, scipy.sparse.linalg.eigsh is typically the best in all respects. For large matrices:

1. scipy.sparse.linalg.eigsh is the most expensive, scaling cubically, but never fails.
2. arpack shift and invert typically scales quadritically, and may fail and miss eigenvalues
3. lobpcg without amg scales linearly (with the number of nnz in the space matrix), but the convergence may be slow. The latest version [MRG] multiple stability updates in lobpcg scipy/scipy#10621 fixes stability issues
4. lobpcg with amg should scale linearly and converge very fast in most cases, but needs Fix for spectral clustering error when using 'amg' solver #13707 to be merged and also may need [MRG] multiple stability updates in lobpcg scipy/scipy#10621

(4 should be with amg I think ;)

@jnothman and binds stronger than or so it's parsed as

    if (eigen_solver == 'arpack' or (eigen_solver != 'lobpcg' and
       ((not sparse.isspmatrix(laplacian)) or n_nodes < 5 * n_components)):

(4 should be with amg I think ;)

thanks, fixed

@lobpcg did you see the comment above the arpack magic? That explains what's happening. I don't think it's related to lobpcg.

Do you mean

            laplacian *= -1
            v0 = random_state.uniform(-1, 1, laplacian.shape[0])
            _, diffusion_map = eigsh(
                laplacian, k=n_components, sigma=1.0, which='LM',
                tol=eigen_tol, v0=v0)

It's just a bad way to call arpack shift-and-invert in this case. A good way is just to change sigma, replacing the above with something like

            v0 = random_state.uniform(-1, 1, laplacian.shape[0])
            _, diffusion_map = eigsh(
                laplacian, k=n_components, sigma=-1e-5, which='LM',
                tol=eigen_tol, v0=v0)

to find the eigenvalues closest to zero. You want sigma to be a bit negative, to avoid LU factorization failures, so that arpack never fails. Above I use a safe choice sigma=-1e-5 that should work in single precision as well (if arpack supports it). If you make sigma even more negative, it should slow down the convergence a bit.

I'd still like to see this merged before we go into rewriting the logic.

addes whatsnew, opened #14713