Sorting eigenvectors by their eigenvalues (associated sorting)

Question

I have an unsorted vector of eigenvalues and a related matrix of eigenvectors. I'd like to sort the columns of the matrix with respect to the sorted set of eigenvalues. (e.g., if eigenvalue[3] moves to eigenvalue[2], I want column 3 of the eigenvector matrix to move over to column 2.)

I know I can sort the eigenvalues in O(N log N) via std::sort. Without rolling my own sorting algorithm, how do I make sure the matrix's columns (the associated eigenvectors) follow along with their eigenvalues as the latter are sorted?

Answer 1

Typically just create a structure something like this:

struct eigen { 
    int value;
    double *vector;

    bool operator<(eigen const &other) const { 
        return value < other.value;
    }
};

Alternatively, just put the eigenvalue/eigenvector into an std::pair -- though I'd prefer eigen.value and eigen.vector over something.first and something.second.

Answer 2

I've done this a number of times in different situations. Rather than sorting the array, just create a new array that has the sorted indices in it.

For example, you have a length n array (vector) evals, and a 2d nxn array evects. Create a new array index that has contains the values [0, n-1].

Then rather than accessing evals as evals[i], you access it as evals[index[i]] and instead of evects[i][j], you access it evects[index[i]][j].

Now you write your sort routine to sort the index array rather than the evals array, so instead of index looking like {0, 1, 2, ... , n-1}, the value in the index array will be in increasing order of the values in the evals array.

So after sorting, if you do this:

for (int i=0;i<n;++i)
{
  cout << evals[index[i]] << endl;
}

you'll get a sorted list of evals.

this way you can sort anything that's associated with that evals array without actually moving memory around. This is important when n gets large, you don't want to be moving around the columns of the evects matrix.

basically the i'th smallest eval will be located at index[i] and that corresponds to the index[i]th evect.

Edited to add. Here's a sort function that I've written to work with std::sort to do what I just said:

template <class DataType, class IndexType>
class SortIndicesInc
{
protected:
    DataType* mData;
public:
    SortIndicesInc(DataType* Data) : mData(Data) {}
    Bool operator()(const IndexType& i, const IndexType& j) const
    {
        return mData[i]<mData[j];
    }
};