Jacobi Method for finding eigenvalues of symmetric matrix.
Stephen Montgomery-Smith
I have been looking at the Jacobi Method for finding eigenvalues of a
symmetric matrix. Using google, I found references that said it
converges with quadratic order when the eigenvalues are distinct.
What is known about the speed of convergence when the eigenvalues are
not necessarily distinct?
Thanks, Stephen
P.S. I only learned this method the day before yesterday. It really is
a very cool algorithm!
user923005
wrote:
> (I posted this on sci.math.research, but it didn't cross-post properly!)
>
> I have been looking at the Jacobi Method for finding eigenvalues of a
> symmetric matrix. Using google, I found references that said it
> converges with quadratic order when the eigenvalues are distinct.
That's odd. I thought that the convergence was linear. Is it a
parallel solver?
> What is known about the speed of convergence when the eigenvalues are
> not necessarily distinct?
>
> Thanks, Stephen
>
> P.S. I only learned this method the day before yesterday. It really is
> a very cool algorithm!
Aside:
http://mymathlib.webtrellis.net/matrices/eigen/symmetric.html
Stephen Montgomery-Smith
> On Dec 26, 5:18 pm, Stephen Montgomery-Smith <step...@missouri.edu>
> wrote:
>> (I posted this on sci.math.research, but it didn't cross-post properly!)
>>
>> I have been looking at the Jacobi Method for finding eigenvalues of a
>> symmetric matrix. Using google, I found references that said it
>> converges with quadratic order when the eigenvalues are distinct.
>
> That's odd. I thought that the convergence was linear. Is it a
> parallel solver?
I got a post on sci.math.research pointing me to the Wikipedia entry,
which said it was quadratic even if the eigenvalues are not distinct.
Hans Mittelmann
That is incorrect.
http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
The basic method requires separation of eigenvalues. For further
results, see http://math.fullerton.edu/mathews/n2003/jacobimethod/JacobiMethodBib/Links/JacobiMethodBib_lnk_3.html
Stephen Montgomery-Smith
Any chance of getting a brief summary? Is there an example of a matrix
with non-distinct eigenvalues which converges with non-quadratic order?
Is there an enhanced method that does give quadratic order?
Thanks, Stephen
(Incidentally, maybe user923005 was momentarily confusing this Jacobi
method with the other one that solves sparse linear systems of equations.)
Peter Spellucci
In article <S9udnfn4odnKCOja...@centurytel.net>,
the jacobi eigenvalue method _is_ quadratically convergent, see the original
work of
schoenhage and wilkinson and also the (somewaht abbreviated) discussion
in parlett: the symmetric eigenvalue problem, reprint, page 197-199.
however, although it is indeed an easily described and compactly to code
method you should not forget that it always works on the full A
(or in the tricky implemetation of Rutishauser: the upper triangle)
since zeroes once produced are later destroyed again. hence its usefulness is
restricted to small matrices.
hth
peter
Stephen Montgomery-Smith
So it is quadratically convergent even when the matrix has repeated
eigenvalues? (I can sort of see why it would be true if the matrix has
distinct eigenvalues.)
Also, my application is for 3 by 3 matrices only. Is there a quicker
method for 3 by 3 matrices?
Peter Spellucci
In article <dbWdnfx7fP8oIOba...@centurytel.net>,
Stephen Montgomery-Smith <ste...@math.missouri.edu> writes:
>Peter Spellucci wrote:
(snip)
>> in parlett: the symmetric eigenvalue problem, reprint, page 197-199.
>> however, although it is indeed an easily described and compactly to code
>> method you should not forget that it always works on the full A
>> (or in the tricky implemetation of Rutishauser: the upper triangle)
>> since zeroes once produced are later destroyed again. hence its usefulness is
>> restricted to small matrices.
>> hth
>> peter
>>
>
>So it is quadratically convergent even when the matrix has repeated
>eigenvalues? (I can sort of see why it would be true if the matrix has
>distinct eigenvalues.)
yes, it is quadratically convergent also for multiple eigenvalues.
to make it practically useful it is essential to follow the tricks of
rutishauser to compute the cos and sin values (without computing the angle itself)
and to use the modified transformation rule (type: newrow=oldrow + correction etc)
in order to cope with accumulation of roundoff.
but in you case, there is really a much faster way:
step1: annihilate a_13 by one Givens rotation. (from left and right).
step2: A is now tridiagonal.
compute one eigenvalue using the bisection scheme
(gaussian elimantion on A-mu*I, counting the negative values =
number of eigenvalues <mu, replace zero pivot by eps.
take the gerschgorin discs for inclusion of all eigenvalues, then bisect
until one eigenvalue is found.
step3: compute one eigenvector by applying "inverse iteration" to A-lambda*I
with (0,0,1) as initial vector.
step4: apply a Householdertransformation which takes the eigenvector to (1,0,0)
as similarity transform to A (this is one step of transformation to Schur Normalform)
step: now A is block diagonal. on the lower right yuo have a 2 by 2 minor
from which the eigenvalues and eigenvectors are found directly.
step6: backtransform the eigenvectors.
here every step has operation count order 3 multiplications and 3 additions.
no more than six small loops and a little preparation
hth
peter
Stephen Montgomery-Smith
> In article <dbWdnfx7fP8oIOba...@centurytel.net>,
> Stephen Montgomery-Smith <ste...@math.missouri.edu> writes:
>> Peter Spellucci wrote:
>
> (snip)
>
>>> in parlett: the symmetric eigenvalue problem, reprint, page 197-199.
>>> however, although it is indeed an easily described and compactly to code
>>> method you should not forget that it always works on the full A
>>> (or in the tricky implemetation of Rutishauser: the upper triangle)
>>> since zeroes once produced are later destroyed again. hence its usefulness is
>>> restricted to small matrices.
>>> hth
>>> peter
>>>
>> So it is quadratically convergent even when the matrix has repeated
>> eigenvalues? (I can sort of see why it would be true if the matrix has
>> distinct eigenvalues.)
>
>
> yes, it is quadratically convergent also for multiple eigenvalues.
> to make it practically useful it is essential to follow the tricks of
> rutishauser to compute the cos and sin values (without computing the angle itself)
> and to use the modified transformation rule (type: newrow=oldrow + correction etc)
> in order to cope with accumulation of roundoff.
> but in you case, there is really a much faster way:
Thank you very much. I have a couple of questions.
But first let me say these things. It would take me quite a while to
understand the method you provided to me, since I lack expertize, but it
definitely looks interesting. Let me also say that the reason I wanted
to write my own 3 by 3 eigenvalue solver was because (a) using
BLAS/LAPACK routines was awkward because it is a little tricky to mix C
with FORTRAN and (b) the stuff at http://www.gnu.org/software/gsl comes
with a GNU license, and I really want to release my final code under a
BSD license. Super-duper speed wasn't essentially, because
diagonalizing the matrix was about 0.1% of the total calculation, so a
lot of my questions are motivated by curiosity rather than need.
My first question to you is this - does BLAS/LAPACK do it the way you
suggest?
My second question is getting back to the Jacobi method. I have to
admit that I found the explanations in terms of angles to be rather
awkward, because what you are really doing is diagonalizing each 2 by 2
diagonal submatrix. So it would seem to me that the tricks of
Ritishauer (presuming that I am identifying them correctly) really boil
down to just going through the steps of finding eigenvalues and
eigenvectors of a 2 by 2 symmetric matrix.
Anyway, this is the code I came up with. I wrote this on Christmas day
while my kids were playing with their new toys. I appreciate that it
probably isn't super optimal, but I was quite proud of it, and maybe it
presents a slightly different viewpoint than the usual. Note I try
quite hard to reduce floating point errors associated with diagonalizing
the 2 by 2 submatrix.
#include <math.h>
/*
* Find the eigenvalues of a symmetric matrix A using Jacobi Method.
* The eigenvalues are returned in the diagonal entries of A.
*/
void diagonalize_sym(int n, double A[n][n]) {
int i, j, k, done;
double discr, aii, aij, ajj, aik, ajk, v1, v2, norm;
do {
done = 1;
for (i=0;i<n;i++) for (j=0;j<i;j++) {
/* Diagonalize the diagonal 2 by 2 sub-matrix along the ith and jth axes. */
aii = A[i][i];
aij = A[i][j];
ajj = A[j][j];
if (fabs(aij)>1e-100) {
discr = hypot(aii-ajj,2*aij);
/* Compute eigenvalues. */
if (aii+ajj>0) {
A[i][i] = (aii+ajj+discr)/2;
A[j][j] = (aii*ajj-aij*aij)/A[i][i];
} else {
A[j][j] = (aii+ajj-discr)/2;
A[i][i] = (aii*ajj-aij*aij)/A[j][j];
}
A[i][j] = A[j][i] = 0;
/* Compute normalized eigenvector corresponding to first eigenvalue. */
if (aii>ajj) {
v1 = (aii-ajj+discr)/2;
v2 = aij;
} else {
v1 = aij;
v2 = (ajj-aii+discr)/2;
}
norm = hypot(v1,v2);
v1 /= norm;
v2 /= norm;
/* Apply change of basis to the rest of the matrix. */
for (k=0;k<n;k++) if (k!=i && k!=j) {
aik = A[i][k];
ajk = A[j][k];
A[i][k] = A[k][i] = v1*aik+v2*ajk;
A[j][k] = A[k][j] = -v2*aik+v1*ajk;
}
done = 0;
}
}
} while (!done);
}
Gordon Sande
<ste...@math.missouri.edu> said:
As a latecomer to this discussion I wonder why if all you want are the
three eigenvalues of a symmetric 3x3 matrix without several auxiliary
concerns why you don't just expand out the characteristic polynomial
(it is a cubic!) and solve it? You are still down in the sizes with
algebraic closed forms available. You seem to have wondered of the path
of simplicity with concerns that apply to much larger problems.
It is like figuring out which are the best passenger buses when all
you are trying to do is walk next door without even the complication
of having to cross the street. It is fun but what is the point!
Stephen Montgomery-Smith
Have you ever seen the formula for the cubic written out in full? It is
a mess.
Also, my personal experience (but I admit that it is very small
experience) is that the cubic formula is not very numerically stable.
(Or maybe that was the quartic formula.) So if I was going to solve the
cubic, I would actually use Newton-Raphson. And then you have to worry
about finding the other two roots by some kind of deflation. So in the
end, while solving the cubic looks attractive, the other methods end up
looking much nicer!
Peter Spellucci
In article <VLOdnUPdw8Kb9uPa...@centurytel.net>,
Stephen Montgomery-Smith <ste...@math.missouri.edu> writes:
>Peter Spellucci wrote:
>> In article <dbWdnfx7fP8oIOba...@centurytel.net>,
>> Stephen Montgomery-Smith <ste...@math.missouri.edu> writes:
>>> Peter Spellucci wrote:
>>
>> (snip)
>suggest?
>
no. but Lapack also contains these techniques. the first step - reduction
to tridiagonal form , is in lapack, then you also could call for ql for
tridiagonal matrices again in lapack, but I think in this special case
the technique is simpler than you want to write it yourself.
>
>My second question is getting back to the Jacobi method. I have to
>admit that I found the explanations in terms of angles to be rather
>awkward, because what you are really doing is diagonalizing each 2 by 2
>diagonal submatrix.
not by using the eigenvaectors of this 2 by 2 matrix !
this seems o.k. formally but contains just the points Rutishauser avoids in order to
keep the computation stable. for your convenience I include here a matlab translation
of rutishausers code, you can see there that c (=cos(phi)) and s(=sin(phi))
are computed without computing phi at all and that the formula for the
transformation is a little bit rewritten to make it more robust aginst
roundoff accumulation.
function [v,d,history,historyend,numrot]=jacobi(a_in,itermax)
% [v,d,history,historyend,numrot]=jacobi(a_in,itermax)
% computes the eigenvalues d and
% eigenvectors v of the real symmetric matrix a_in,
% using rutishausers modfications of the classical
% jacobi rotation method with treshold pivoting.
% history(1:historyend) is a column vector of the length of
% total sweeps used containing the sum of squares of
% strict upper diagonal elements of a. a is not
% touched but copied locally
% the upper triangle is used only
% itermax is the maximum number of total sweeps allowed
% numrot is the number of rotations applied in total
% check arguments
siz=size(a_in);
if siz(1) ~= siz(2)
error('jacobi : matrix must be square ' );
end
if norm(a_in-a_in',inf) ~= 0
error('jacobi ; matrix must be symmetric ');
end
if ~isreal(a_in)
error(' jacobi : valid for real matrices only');
end
n=siz(1);
v=eye(n);
a=a_in;
history=zeros(itermax,1);
d=diag(a);
bw=d;
zw=zeros(n,1);
iter=0;
numrot=0;
while iter < itermax
iter=iter+1;
history(iter)=sqrt(sum(sum(triu(a,1).^2)));
historyend=iter;
tresh=history(iter)/(4*n);
if tresh ==0
return;
end
for p=1:n
for q=p+1:n
gapq=10*abs(a(p,q));
termp=gapq+abs(d(p));
termq=gapq+abs(d(q));
if iter>4 & termp==abs(d(p)) & termq==abs(d(q))
% annihilate tiny elements
a(p,q)=0;
else
if abs(a(p,q)) >= tresh
%apply rotation
h=d(q)-d(p);
term=abs(h)+gapq;
if term == abs(h)
t=a(p,q)/h;
else
theta=0.5*h/a(p,q);
t=1/(abs(theta)+sqrt(1+theta^2));
if theta < 0
t=-t;
end
end
c=1/sqrt(1+t^2);
s=t*c;
tau=s/(1+c);
h=t*a(p,q);
zw(p)=zw(p)-h; %accumulate corrections to diagonal elements
zw(q)=zw(q)+h;
d(p)=d(p)-h;
d(q)=d(q)+h;
a(p,q)=0;
%rotate, use information from the upper triangle of a only
%for a pipelined cpu it may be better to work
%on full rows and columns instead
for j=1:p-1
g=a(j,p);
h=a(j,q);
a(j,p)=g-s*(h+g*tau);
a(j,q)=h+s*(g-h*tau);
end
for j=p+1:q-1
g=a(p,j);
h=a(j,q);
a(p,j)=g-s*(h+g*tau);
a(j,q)=h+s*(g-h*tau);
end
for j=q+1:n
g=a(p,j);
h=a(q,j);
a(p,j)=g-s*(h+g*tau);
a(q,j)=h+s*(g-h*tau);
end
% accumulate information in eigenvectormatrix
for j=1:n
g=v(j,p);
h=v(j,q);
v(j,p)=g-s*(h+g*tau);
v(j,q)=h+s*(g-h*tau);
end
numrot=numrot+1;
end
end %if
end % for q
end % for p
bw=bw+zw;
d=bw;
zw=zeros(n,1);
end %while