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

The solution purely relies on the way you store your eigenvector matrix.

The best performance while sorting will be achieved if you can implement swap(evector1, evector2) so that it only rebinds the pointers and the real data is left unchanged.

This could be done using something like double* or probably something more complicated, depends on your matrix implementation.

If done this way, swap(...) wouldn't affect your sorting operation performance.

Answer 3

The idea of conglomerating your vector and matrix is probably the best way to do it in C++. I am thinking about how I would do it in R and seeing if that can be translated to C++. In R it's very easy, simply evec<-evec[,order(eval)]. Unfortunately, I don't know of any built in way to perform the order() operation in C++. Perhaps someone else does, in which case this could be done in a similar way.

Answer 4

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];
    }
};