Look into LOBPCG (sparse diagonalizer method)

LOBPCG is a sparse matrix diagonalization method that produces the minimal eigenvalue. This could possibly be a nice alternative to ARPACK for detecting negative modes? It's worth looking into, anyways.

I see conflicting statements on wikipedia about the suitability of this method for our problem: (emphasis mine)

Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) is a matrix-free method for finding the largest (or smallest) eigenvalues and the corresponding eigenvectors of a symmetric positive definite generalized eigenvalue problem

A simple work-around is to negate the function, substituting -DT(D X) for DT(D X) and thus reversing the order of the eigenvalues, since LOBPCG does not care if the matrix of the eigenvalue problem is positive definite or not.

(I think the first quote might just mean that B (the extra matrix in the generalized eigenvalue problem) is positive definite, not A...)

Some choice quotes (note: ARPACK uses something similar to the Lanczos method, with "implicit restarting" to reduce memory requirements):

In contrast to the Lanczos method, LOBPCG rarely exhibits asymptotic superlinear convergence in practice

LOBPCG allows warm starts and computes an approximation to the eigenvector on every iteration. It has no numerical stability issues similar to those of the Lanczos method.

LOBPCG is reasonably easy to implement, so many implementations have appeared.

(with that last quote: All it does is minimize the function (x^T A x) / (x^T x), with some special considerations for numerical stability)

~~Har har har. Nope.~~

I took gamma-dynmat-01.npz from 1261-3-a (known to have a buckling mode) and tried asking scipy.sparse.linalg.lobpcg to take care of it. I tried asking for 1 mode, and later asking for 10. With a random initial state, it always produces eigenvalues around 10.9. If you start with a vector full of ones, then it almost immediately returns with eigenvalue 10^-16... unless you set tol=1e-10, niter=100, in which case it once again finds the modes around 10.9. (why did it move from a smaller value to a larger value!?)

~~There aren't many options I could potentially change to fix this; I can try supplying a preconditioner, but I think that only improves performance and not reliability?~~

I am a doofus

I also have an experimental implementation that simply uses rsp2's conjugate gradient on the rayleigh quotient. If this is ran on a random state, it quickly finds the acoustic modes and then gets stuck for thousands of iterations as the value decreases by tiny amounts around 1e-8. Eventually the step size grows rather alarmingly to values like 1e28. I haven't tried letting it run all the way to completion yet but I doubt the end result will be any good...

(it's possible my gradient has the wrong scale factor? I haven't checked it numerically)

I haven't tried letting it run all the way to completion yet but I doubt the end result will be any good...

...the hell? It worked:

[5566.134s][TRACE]  i:  12239  v:  -0.000  dv: -5.9e-18  g: 4.2e-24  max: 1.8e-25  a:  1.2e6 cos: +0.99 +0.99 +0.99
[5566.443s][TRACE]  i:  12240  v:  -0.000  dv: +1.0e-17  g: 4.2e-24  max: 1.6e-25  a:  1.3e6 cos: +0.99 +0.99 +0.99
[5566.754s][TRACE]  i:  12241  v:  -0.000  dv: -2.4e-18  g: 4.1e-24  max: 1.5e-25  a:  1.3e6 cos: +0.99 +0.99 +0.99
[5567.064s][TRACE]  i:  12242  v:  -0.000  dv: -2.4e-18  g: 4.9e-24  max: 1.8e-25  a:  1.5e6 cos: +1.00 +0.99 +0.99
[5567.376s][TRACE]  i:  12243  v:  -0.000  dv: +1.0e-17  g: 5.0e-24  max: 2.0e-25  a:  1.7e6 cos: +1.00 +0.99 +0.99
[5567.700s][TRACE]  i:  12244  v:  -0.000  dv: -1.2e-17  g: 4.8e-24  max: 1.7e-25  a:  1.7e6 cos: +1.00 +0.99 +0.99
[5568.013s][TRACE]  i:  12245  v:  -0.000  dv: +2.5e-17  g: 4.9e-24  max: 1.9e-25  a:  1.8e6 cos: +1.00 +0.99 +0.99
[5568.014s][ INFO] ACGSD Finished.
[5568.014s][ INFO] Iterations: 12245
[5568.014s][ INFO]      Value: -0.00002049885474252192
[5568.014s][ INFO]  Grad Norm: 4.932045034235915e-24
[5568.014s][ INFO]   Grad Max: 1.875005300249886e-25
[5568.022s][TRACE] Done diagonalizing dynamical matrix
[5568.027s][TRACE] ============================
[5568.027s][TRACE] Finished diagonalization
[5569.035s][TRACE] Examining mode 1 (-2.3609940) (ddot = 0.999987)...
[5569.035s][TRACE] Computing eigensystem info
[5569.036s][TRACE] computing EvAcousticness
[5569.036s][TRACE] computing EvPolarization
[5569.037s][TRACE] computing EvLayerAcousticness
[5569.037s][TRACE] not computing UnfoldProbs due to missing requirement RequestToUnfoldBands
[5569.037s][TRACE] computing EvRamanTensors
[5569.038s][ INFO] # (C)  Frequency(cm-1)    Acoust. RamnA RamnB  Layer  [ X  ,  Y  ,  Z  ]
[5569.038s][ INFO]   (‼) -2.3609939680669703   1e-09 0     0       1e-05 [0.00, 0.00, 1.00]
[5569.409s][TRACE] Chasing 1 bad eigenvectors...

That is indeed around the frequency of the negative modes in the system. Notice how the value of the objective function (v column) is very small; this is likely due to the eigenvector getting a large magnitude, which appears in the denominator of the rayleigh quotient. I'm not sure whether that slowed it down or what. (no dummy, that's because it's eigenvalue vs frequency)

Almost certainly, it probably ran far, far more steps than it needed to, because I gave it a stop condition which was overkill and then some:

  eigensolver:
    experimental-lobpcg:
      cg:
        stop-condition:
          value-delta:
            rel-greater-than: 0
            steps-ago: 50

Killed it after 100 iterations from a random state and it produced an expectation value of 0.4.

Probably best to start with something near the acoustic modes, with a small random perturbation...

Uhhhhhhhhhhhhhh

I don't think I can facepalm hard enough. An eigenvalue of 10.9 apparently corresponds to a frequency of around 1700 wavenumber. Turns out I was asking scipy to find the LARGEST eigenvalues. Yep. (in my defense, the documentation doesn't mention the default value of largest, and is written in a way that suggests it defaults to False...)

That said, now that I'm correctly asking for the smallest eigenvectors, the results of scipy.sparse.linalg.eigen.lobpcg look quite a bit more promising on the same dynamical matrix that was tested before. A fairly long run with n=3, maxiter=100000, tol=1e-13 takes about 700 iterations to find a rayleigh quotient of about the magnitude of the true most negative mode, takes about 1500 iterations to converge onto the eigenvalue up to 3 digits, and... well, never really reaches 1e-13 tolerance, so these are poor stop conditions.

More research is still necessary to find good parameters.

tol - Tiny tolerances waste a lot of time for little benefit. I don't think we need the eigenvalues and eigenvectors to be very precise; as long as the Rayleigh quotient is sufficiently negative, it's probably a good direction to relax along. BUT, that does not mean that we can use a large tolerance like 1e-3 or even the default of n*sqrt(eps); large tolerances generally cause it to stop on a positive eigenvalue. Basically we need an idea of what is the largest we can push tol without it potentially stopping too soon?
maxiter - I'm not sure I like this as a stopping condition either, since it may need to scale with system size. This will require some methodological benchmarking.
n - it's not yet clear to me whether increasing or decreasing the number of requested kets improves convergence time for the most negative eigenvalue. (I think there might be a sweet spot at n=2; n=2 seems to reach landmarks in orders of magnitude fewer iterations than n=1 while only taking twice as long each iteration, while n larger than 2 don't seem to significantly affect iteration count)
starting input - random starting points seem to work pretty good. How does it compare to starting with acoustic modes?

It'd be nice if we could stop on sufficiently negative eigenvalues as a stop condition, but that wouldn't work too well for the case where there's no negative modes!

ExpHP / rsp2

Look into LOBPCG (sparse diagonalizer method) #107