eigenvalues and vectors using SVdec

Steve Pollaine pointed out that it is relatively easy to get the eigenvalues and eigenvectors of a symmetric matrix using the SVdec function. The singular values are the absolute values of the eigenvalues, and the singular vectors are the eigenvectors, except for the case of equal singular values corresponding to eigenvalues of opposite sign. Two tricks take care of these two problems: (1) By examining the columns of the U and V matrices returned by SVdec, you can figure out the sign of the eigenvalue. And (2), by adding a constant times the identity matrix, you can offset all the eigenvalues by that constant amount, separating any opposite-but-equal-magnitude pairs of eigenvalues.

Steve came up with a very short function only using trick (2): simply add a large enough constant that you ensure all eigenvalues are positive, and therefore equal the singular values -- subtract your constant and you are done. Unfortunately, the only obvious (to me or Steve) way to select a large enough constant is to call SVdec a first time, and use the biggest singular value as your constant for a second call.

So I offer the following compromise function, which uses both tricks, but is not mathematically certain to work (although I'm guessing failure is about as likely as SVdec itself failing). I may go ahead and put this into matrix.i in the distribution, if I can stand the inelegance...

| | |

| **Code:**
func SVeigen(a, &u) /* DOCUMENT s = SVeigen(a, u) * return eigenvalues S and eignvectors U of symmetric matrix A. * That is, * A(,+) * U(+,) = S(-,) * U * The eigenvalues in S are in decreasing order of absolute value. * The eigenvectors are orthonormal: * U(+,) * U(+,) = unit(dimsof(A)(2)) * If A is not symmetric, SVeigen returns garbage. * SEE ALSO: SVdec */ { local vt; n = dimsof(a)(2); u = array(0., n, n); u(1:numberof(u):n+1) = x = 0.4964726826642538*avg(abs(a)); /* this would work with x=0 except for case that A has * equal singular values corresponding to eigenvalues of opposite sign * adding essentially random x to all eigenvalues makes this very unlikely */ s = SVdec(a+u, u, vt); /* singular value = abs(eigenvalue) */ i = abs(u+transpose(vt))(sum,) > abs(u)(sum,); /* eigenvalue > 0 */ s = (i+i-1.0)*s - x; /* sort unnecessary when x=0 */ list = sort(-abs(s)); s = s(list); u = u(,list); return s; }
| |

| | |