How to get the same eigenvectors of a degenerate eigenvalue as the matrix evolves slightly? I have a $n\times n$ symmetric normal complex matrix which always has $\frac 23 n$ degenerate eigenvalues. ($n$ is always in form of $3k \quad k=1,2,...)$
For example in $6\times 6$ case it has eigenvalues $\lambda_1,\lambda_2,\lambda_3,\lambda_3,\lambda_4,\lambda_4$.
I know a degenerate eigenvalue will create a $2D$ space of eigenvectors that can be associated to it. I need to measure the evolution of these eigenvectors as the matrix changes very slightly at each step (with a loop variable called $\omega$), which I have a computer code written for it. It calculates the following formula for eigenvectors and plots $\tau-\omega$.
$$f_{ij}(\alpha,\beta)=\sum_{l=1}^nF(\omega,\lambda_l)\frac {\langle i\alpha\lvert l\rangle \langle\bar l\lvert j\beta\rangle}{\langle l\lvert \bar l \rangle}$$
$$\tau_{ij}(\omega)=C(\omega)\sum_{\alpha \beta}\lvert f_{ij}(\alpha,\beta)\lvert^2$$
Where $\langle i\alpha\lvert$ and $\lvert j\beta\rangle$ are basis of an $nD$ space, $\langle \bar l\lvert$ is the notation for eigenvector $\lvert l\rangle$ transposed (not conjugated), $(i,j)$  are indexes from $1$ to $n$, so we would have $n^2$ values of $\tau$ changing as $\omega$ changes. (and $\tau_{ij}=\tau_{ji}$)
Problem is the computer code choosing one of infinite possible eigenvectors for degenerate eigenvalues in each iteration as $\omega$ changes so the final $\tau-\omega$ plot has noises.
First image is what the result is supposed to look like and the second one indicates the problem. These are not high-res outputs but  I think they do the job.


These noises vanish as I use only the non-degenerate eigenvalues and vectors or when I use a very bad approximation to make the matrix Hermitian. So the degeneracy problem is the case here.
Question: Is there any way around this problem to monitor the behavior of degenerate eigenvectors as the matrix changes slightly?
(I use LAPACK to calculate eigenvectors if it matters)
 A: Supposing the multiplicity $k$ of an eigenvalue $\lambda_i$ is fixed and known, you can represent the $k$-dimensional eigenspace as a unique line though $0$ in the exterior power $\wedge^k(V)$ of you $n$-dimensional space (which power has dimension $\binom nk$). Then you can represent the line by one of the two unit vectors it contains; hopefully you can control which one is chosen in a reasonable way. In this way you suppress the arbitrariness of the choice of a basis for the eigenspace.
For more details about how to do this, see Plücker embedding and Grassmann coordinates.
A: (Too long for a comment.)
Suppose you are able to trace the eigenvalues and the multiplicity of each eigenvalue remain unchanged. That is, suppose $A(t)$ is a normal matrix that varies continuously in time $t$, and $\lambda(t)$ is one of its eigenvalue that varies continuously in $t$ and whose multiplicity remain unchanged over some time window. Let $B$ and $C$ be two matrices, each with orthonormal columns, such that the columns of $B$ are the unit eigenvectors of $\lambda(t_1)$ and the columns of $C$ is a spanning set for the eigenspace for $\lambda(t_2)$ at some time point $t_2$ close to $t_1$. Then you are essentially asking how to find a set of unit vectors that spans the column space $C$ and is also close to the columns of $B$. In other words, you need to find a unitary matrix $Q$ that minimises the Frobenius norm of $B-CQ$.
If this is the case, the solution to the above classical least square problem is $Q=VU^\ast$, where $USV^\ast$ is a singular value decomposition of the matrix product $B^\ast C$. So, you may pick the columns of $B$ as the eigenvectors for $\lambda(t_1)$ and the columns of $CQ$ as the eigenvectors for $\lambda(t_2)$.
Things get nastier, however, if the eigenvalue paths cross each other.
