sparse diagonalization

[ExpHP]

Just recording notes here

[ExpHP]

Turns out ARPACK + shift-invert can give non-eigenvectors as results:

https://gist.github.com/ExpHP/3faeedfaa3b75236eb8036f48dbc2d56

To my best understanding, this matrix has no true negative non-acoustic modes.

Here I look at abs(np.vdot(normalize(m @ v), normalize(v))) to estimate the quality of an eigenvector. (for a true eigenvector it should be 1, except for acoustics where it is ill-defined and m @ v - val * v is more reliable). I compare the modes produced by which='SA' with and without shift-invert mode.

Shift-invert mode sometimes finds a negative mode (once it finds one around -6 cm, another time around -9cm). These modes fail both tests miserably.

Modes found by non-shift-invert all satisfy both tests (except for one acoustic mode which fails the normalize test because its eigenvalue is so impressively close to zero)

[ExpHP]

Next question: when there is a legitimate negative mode, does shift-invert give a good eigenvector?

[ExpHP]

Structure with two legitimate negative modes at -3.08101606647

$ python3
Python 3.6.5 (default, May 11 2018, 04:00:52) 
[GCC 8.1.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> d=json.load(open('whyyyyyyy-019-a/gamma-dynmat.json'))
>>> import scipy.sparse.linalg as spla
>>> from scipy.sparse import bsr_matrix
>>> m=bsr_matrix((np.array(d['complex-blocks'])[:,0], d['col'], d['row-ptr']), tuple(3*x for x in d['dim']))
>>> def f():
...   try:
...     return spla.eigsh(m, k=8, OPinv=OPinv, sigma=0, maxiter=30, which='SA')
...   except scipy.sparse.linalg.eigen.arpack.ArpackNoConvergence as e:
...     return (e.eigenvalues, e.eigenvectors)                         
... 
>>> OPinv = spla.LinearOperator(m.shape, matvec=spla.eigen.arpack.get_OPinv_matvec(m, M=None, sigma=0, tol=0), dtype=m.dtype)
/home/lampam/asd/clone/scipy/scipy/sparse/linalg/dsolve/linsolve.py:295: SparseEfficiencyWarning: splu requires CSC matrix format
  warn('splu requires CSC matrix format', SparseEfficiencyWarning)
>>> evals,evecs = f(); len(evals)
31 iterations
5
>>> 15.6333043006705 * 33.3564095198152 * np.abs(evals) ** 0.5 * np.signum(evals)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: module 'numpy' has no attribute 'signum'
>>> 15.6333043006705 * 33.3564095198152 * np.abs(evals) ** 0.5 * np.sign(evals)  
array([-3.08201022e+00, -3.08101607e+00, -8.17189576e-03, -8.17189139e-03,
       -3.52617685e-05])
>>> def normalize(a): return a / np.sqrt(np.sum(np.square(a)))
... 
>>> np.vdot(normalize(m @ evecs[:,0]), normalize(evecs[:,0]))
0.013210816815992654
>>> np.vdot(normalize(m @ evecs[:,1]), normalize(evecs[:,1]))
-0.9999995816404131
>>> np.vdot(normalize(m @ evecs[:,2]), normalize(evecs[:,2]))
-7.86575381979795e-05
>>>

Shift-invert gets one of them right on the nose, but the other one has a slightly wrong eigenvalue and a totally garbage eigenvector.

[ExpHP]

triage: not an actual issue