I have $n$-dimensional matrices $\mathrm{\hat{H}}(\vec{k})$ depending on vector parameter $\vec{k}$.
Now, eigenvalue routines return eigenvalues in no particular order (they are usually sorted), but I want to trace eigenvalues $E_i$ as smooth functions of $\vec{k}$. Because eigenvalues are not returned in any particular order, just tracing $E_i$ for some particular index $i\in\{1,..,n\}$ will return set of lines which are not smooth, as shown on the picture bellow
My idea to trace continuous lines was to use eigenvectors. Namely, for two close points $\vec{k}$ and $\vec{k}+d\vec{k}$ eigenvectors should be approximately orthonormal so that $v_i(\vec{k})\cdot v_j(\vec{k}+d\vec{k})\sim \delta_{p_i p_j}$ where $p_i, p_j\in \pi(\{1,...,n\})$, and $\pi$ is some permutation. Then I would use given permutation to reorder the eigenvalues and thus trace smooth lines.
I other words, I would trace continuity of eigenvectors.
However, I run into some problems with numerical routines. At a given small subset of points I use, few eigenvectors at nearby points are not almost orthonormal. My first suspicion was that those eigenvectors correspond to a degenerate eigenvalue, but that is not always true.
This also holds true if I reduce $d\vec{k}$ to be really small.
Is such thing allowed to happen. Or, is it possible to guarantee that numerical routines return continuous eigenvectors? Routine I use is numpy.linalg.eigh which is an interface for zheevd from LAPACK.
(Physicists amongst you will recognize that I am talking about the band structure)
numpy.linalg.svd
to generate your eigenvectors. At least in Matlab, the underlying routine forsvd
always returns the eigenvalues and eigenvectors in decreasing order. $\endgroup$