Finding a unique definition of eigenspace for degenerate eigenvalues I am diagonalizing an Hermitian matrix numerically, and this results in a set of eigenvalues $\varepsilon_1, ..., \varepsilon_n$ and a set of eigenvectors $v_1, ..., v_n$.
If the eigenvalues $\varepsilon_i = \varepsilon_{i+1} = \varepsilon_{i+2}$ are degenerate this results in an eigenspace, spanned by $v_i, v_{i+1}, v_{i+2}$.
The Problem is, that unlike the eigenvalues, $v_i, v_{i+1}, v_{i+2}$ are not uniquely defined and they differ between different Lapack and ScaLapack implementations, which makes debugging very hard.
Is there an established algorithm to transform the eigenspace, so that it always results in a
presentation that is:

*

*still orthonormal

*unique regardless of the original state

I have been toying with different ideas, but I would prefer to use something established, that then is also numerically stable etc.
 A: Imagine that if we were given any vector space $V$, we had a heuristic by which to "order" the non-zero vectors of $V$. Now suppose (like the problem asks for) that we wanted a procedure so that, no matter what spanning set $S_1$ of $V$ we input, we will get the same (ordered) orthonormal basis $B_2$ of $V$. Let us call this process "finding the canonical basis" of $V$. We could devise a procedure to find the canonical basis of $V$ just by doing the following:

*

*Start with the "first" vector $v_1 \in V$ (where "first" is defined by the order we give to the vectors in $V$), and normalize it. This vector will $\hat{v}_1$ will be the first element in our canonical basis.


*Find the canonical basis of $V_2 = \{v ~|~ v \in V \text{ and } \langle v, v_1\rangle = 0\}$ (in other words, the subspace of $V$ orthogonal to $v_1$). Given $B_1$, finding a spanning set for $V_2$ is easy: we just subtract off the projection onto $v_1$ of every element.
Since in every iteration $V_2$ has dimension one less than $V$, we only need to iterate $\dim V$ times.
It remains to determine a procedure by which to find the "first" vector in some arbitrary subspace $V$, given any set of vectors $\{v_1, v_2, \cdots, v_n\}$ which span $V$. Fortunately, the RREF gives us a way of doing so! Let us construct the matrix $M$ where $v_1, \cdots, v_n$ are the rows. Then we can define the "first" vector in a subspace to be the first row of $\text{rref}(M)$. Since $\text{rref}(M)$ is the unique rref of any matrix whose row space is $V$, we can be guaranteed that "first" vector by our definition will be the same for any $\{v_1, v_2, \cdots, v_n\}$ that span $V$.
Some (rough) pseudocode for the algorithm:

Given: $S_1$, a set of vectors spanning $V$
Output: Canonical Basis of $V$
$S \gets S_1$
$B \gets [~]$
while there exists non-zero $v_i \in S$:
$\quad$Remove all zero vectors from $S$
$\quad$$M \gets $ the matrix whose rows are the elements of $S$
$\quad$$r_1 \gets $ first row of $\text{rref}(M)$
$\quad$Add $r_1 / \|r_1\|$ to $B$
$\quad$$S \gets$ all non-zero rows of $\text{rref}(M)$
$\quad$for each $v_i \in S$:
$\qquad$$v_i \gets v_i - \frac{\langle v_i, ~r_1 \rangle}{\|r_1\|^2} r_1$
return $B$

I am sure optimizations can be made (especially if the input is always an orthonormal basis of the subspace), but this is a working start!
A: I came up with a similar solution. First I find the vector in the space which is "closest" to some target vector using a least squares and the I perform a gram-schmidt on it, which gives me an orthonormal basis. Here I use a set of sinus-functions as target vectors:
def t(n):
    tar = np.sin(np.arange(old_eig_vecs.shape[0])*(n+1) * old_eig_vecs.shape[0]/2.0)
    return tar

def unify_eigvec(old):
    new = np.zeros(old.shape)
    for i in range(3):
        target_vec = t(i)
        res = np.linalg.lstsq(old, target_vec, rcond=None)
        new[:,i] = old@res[0]
    return gram_schmidt_columns(new)

def gram_schmidt_columns(X):
    Q, _ = np.linalg.qr(X)
    return Q
```

