Suppose we are interested in the eigenvalues and eigenvectors of a hermitian matrix h(t) that depends on a parameter t. My matrix is large and sparse and hence needs to be treated numerically.
A naive approach is to evaluate the matrix h(t_k) at discretized parameter values t_k. Is it possible to sort the eigenvectors and eigenvalues according to the "character of the eigenvector"?
Let me illustrate what I mean by "character of the eigenvector" with the following simple example (i denotes the imaginary unit).
h(t) = {{1, it}, {-it, 1}}
The eigenvalues are 1-t and 1+t with the corresponding eigenvectors {-i, 1} and {i, 1}. Hence sorting according to the "eigenvector character" the eigenvalues should cross at t = 0. However most eigensolvers sort them by increasing eigenvalues exchanging the eigenvectors from negative to positive t (see code and output plot).
import numpy as np
import scipy.sparse.linalg as sla
import matplotlib.pyplot as plt
def h(t):
# parametrized hermitian matrix
return np.array([[1, t*1j], [-t*1j, 1]])
def eigenvalues(t):
# convert to tuple for np.vectorize to work
return tuple(sla.eigsh(h(t), k=2, return_eigenvectors=False))
eigenvalues = np.vectorize(eigenvalues)
t = np.linspace(-1, 1, num=200)
ev0, ev1 = eigenvalues(t)
plt.plot(t, ev0, 'r')
plt.plot(t, ev1, 'g')
plt.xlabel('t')
plt.ylabel('eigenvalues')
plt.show()
Idea Most eigensolvers iteratively approximate the eigenvectors and eigenvalues. By feeding the eigensystem of the matrix h(t_k) to the solver as an initial guess to diagonalize h(t_{k+1}) one might obtain a result ordered by the "character of the eigenvector".
Is it possible to achieve this with scipy or more generally with python? Preferentially the heavy diagonalization should be delegated to a dedicated compiled library (e.g. Lapack in scipy). Is there a suitable Lapack routine maybe already wrapped in scipy?
Is there an alternative method that achieves the same? How can it be implemented in python?
