eigenvalues_sym and eigenvectors_sym currently return the values, and their associated vectors, in ascending order. How do I request an option be added to return the eigen decomposition as-is without sorting?
Looking at the c++ on the back-end, we don’t do any sorting for these functions, and just return the values calculated directly by Eigen, using:
Eigen::SelfAdjointEigenSolver<PlainMat> solver(m_eval,
Eigen::EigenvaluesOnly);
return solver.eigenvalues();
Then looking at Eigen’s documentation for that solver, the ordered eigenvalues appear to be correct. The full example from that link below:
Here is a random symmetric 5x5 matrix, A:
1.36 -0.816 0.521 1.43 -0.144
-0.816 -0.659 0.794 -0.173 -0.406
0.521 0.794 -0.541 0.461 0.179
1.43 -0.173 0.461 -1.43 0.822
-0.144 -0.406 0.179 0.822 -1.37
The eigenvalues of A are:
-2.65
-1.77
-0.745
0.227
2.29
The matrix of eigenvectors, V, is:
-0.326 -0.0984 0.347 -0.0109 0.874
-0.207 -0.642 0.228 0.662 -0.232
0.0495 0.629 -0.164 0.74 0.164
0.721 -0.397 -0.402 0.115 0.385
-0.573 -0.156 -0.799 -0.0256 0.0858
Consider the first eigenvalue, lambda = -2.65
If v is the corresponding eigenvector, then lambda * v =
0.865
0.55
-0.131
-1.91
1.52
... and A * v =
0.865
0.55
-0.131
-1.91
1.52
Finally, V * D * V^(-1) =
1.36 -0.816 0.521 1.43 -0.144
-0.816 -0.659 0.794 -0.173 -0.406
0.521 0.794 -0.541 0.461 0.179
1.43 -0.173 0.461 -1.43 0.822
-0.144 -0.406 0.179 0.822 -1.37
But if you’ve got an example where Stan’s decomposition returns different eigenvalues to another program, that would be great for debugging
“the ordered eigenvalues appear to be correct” – Yes, no doubt they are correct. However, I don’t want them sorted. I know that Eigen doesn’t offer this option. I had to hack the solver (21.4 KB). I have a particular algorithm in mind that needs the decomposition unsorted.
Oh I see what you mean sorry. For that you should open an issue in the Stan math library, and if you could attach that Eigen code there as well it would help speed things up.
I can add the functionality using the code you provided, but some of the other developers might have a different way of approaching it (e.g., whether it should be a pull to Eigen or Stan Math).
It looks like your only option in the meantime is to use the non-symmetric solver. But if you create that issue we can keep you posted on the progress on adding this.
I’m not sure what this means. There’s no natural or “unsorted” ordering to the eigenvalues. It’s just a set of numbers.
1 Like
By unsorted, I mean that the order of the eigenvalues should match where they appear on the original diagonal. So,
m1 <- matrix(rnorm(9,sd=.1), 3,3) + diag(3)
m1[2,] <- 0
m1[,2] <- 0
eigen(m1)
should result in eigenvalues c(X, 0, Y) instead of what R returns, c(X, Y, 0).
@andrjohns Thanks for your suggestions. I’ll follow up.
That still doesn’t answer the question because eigenvalues are not diagonal elements.
For example the matrix
\left[\begin{array}{cc}
\sin^{2}\theta & \frac{1}{2}\sin2\theta\\
\frac{1}{2}\sin2\theta & \cos^{2}\theta
\end{array}\right]
has eigenvalues 0 and 1 for any value of \theta but the order is different for \theta=0 and \theta=\frac{\pi}{2}. The matrix a continuous function of \theta so at what point should the eigenvalues suddenly change from [0,1] to [1,0]?
If the matrix is almost-diagonal and the eigenvalues are well-separated I think you could re-order the eigenvalues according to the diagonal. So like this:
vector eigenvalues_ordered(matrix x) {
vector[N] y = eigenvalues_sym(x);
// sort twice to invert the permutation
int idxs[N] = sort_indices_asc(sort_indices_asc(diagonal(x)));
return y[idxs];
}
But there’s probably some corner case where this produces the wrong result.
4 Likes
Hm, interesting example. The use case is for covariance matrices where zeros on the diagonal can correspond with zero eigenvalues.
The relevant c++ is on L527 of math\lib\eigen_3.3.7\Eigen\src\Eigenvalues\SelfAdjointEigenSolver.h:
// Sort eigenvalues and corresponding vectors.
// TODO make the sort optional ?
// TODO use a better sort algorithm !!
if (info == Success)
{
for (Index i = 0; i < n-1; ++i)
{
Index k;
diag.segment(i,n-i).minCoeff(&k);
if (k > 0)
{
std::swap(diag[i], diag[k+i]);
if(computeEigenvectors)
eivec.col(i).swap(eivec.col(k+i));
}
}
}
It looks they were planning on making the sort optional, but that it never happened
It might be easier for you to call Eigen’s RealSchur function. If your input matrix is symmetric, then the quasi-triangular matrix will be a diagonal matrix of real eigenvalues. They say that they chase zeros to the bottom right but I seem to remember the code not actually doing that.
OK, I won’t file a feature request until I try these other options. Thank you.