Reply to topic  [ 1 post ] 
eigenvalues and vectors using SVdec 
Author Message
Yorick Master

Joined: Mon Nov 22, 2004 9:43 am
Posts: 354
Location: Livermore, CA, USA
Post 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;
}


Sat May 01, 2010 7:07 pm
Profile
Display posts from previous:  Sort by  
Reply to topic   [ 1 post ] 

Who is online

Users browsing this forum: No registered users and 1 guest


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group.
Designed by STSoftware for PTF.