Thread Subject: order of eigenvalues

Hi All,

I'm working on a problem that requires me to trace the
eigenvalues of the gain matrix of a linearization of a non-
linear ODE. However, the eig() function is causing some
problems.

Say we have the gain matrix in two points of time close
together : t_1 and t_2. For t_2 the gain matrix would be
slightly different from the gain matrix at t_1, eigenvalues
would also be slightly different. Unfortunately, eig() will
sometimes return the eigenvalues at two points in time
close together in a completely different order. For
instance, if the eigenvalues at t_1 were [1; -1+i;-1-i],
then, sometimes, the eigenvalues for t_2 are returned as [-
1.01+.99i; -1.01-.99i; -.99; ]. This makes a comparison
between the eigenvalues at two points in time very
difficult.

There is no natural order to the eigenvalues since they are
complex, so sorting doesn't work. Especially in higher
order systems this is futile. For now I've created a
workaround that sorts the eigenvalues at t_2 based on the
difference with the eigenvalues at t_1:
    
    [evOldmg, evNewmg] = meshgrid(evOld, evNew);
    the_diff = abs(evOldmg - evNewmg);
    % second return value of min function is location,
    % which can be used for permutation.
    [minims, perm] = min(the_diff, [], 1);

    % permute
    evperm = evNew(perm);

This, however, is a bit of a hack and still fails
sometimes. So does an equivalent method of sorting based on
the eigenvectors.

Does anyone know of a method of calculating the eigenvalues
that does not have the order problem? Or, perhaps a better
solution than my temporary workaround?

Regards,

Willem Geert

PS. I've sifted through news-group posts and google, but
could not find a solution to this problem. If I've missed a
previously posted solution, please point me to it.

"Willem Geert " <wgphaff@gmail.com> wrote in message
<fnv66a$cjt$1@fred.mathworks.com>...
> I'm working on a problem that requires me to trace the
> eigenvalues of the gain matrix of a linearization of a non-
> linear ODE. However, the eig() function is causing some
> problems.

As you have found out, eigenvalues have no
natural order that they come out in. They are
just the roots of a polynomial, and if you have
a pair of nearly similar eigensystems, the
corresponding eigenvalues may come out in
an arbitrarily different order.

Having said that, I've considered writing a
shuffle code for eigen-problems of this type.
The idea is that you have more information
available to you than the pure magnitude of
the eigenvalues. You also have the eigenvectors.

Assume that you have two sets of eigenvectors
and eigenvalues, related by the fact that they
came from systems which were close to each
other. Now shuffle the eigenvalues and vectors
from system 2 so that its eigenvectors are
maximally aligned with those of system 1.

This seems like it might be easy, except that
the tiny eigenvalues for some systems may
sometimes end up with virtually random
eigenvectors. I'd visualize some potentially
nasty problems to resolve.

John

 
> Assume that you have two sets of eigenvectors
> and eigenvalues, related by the fact that they
> came from systems which were close to each
> other. Now shuffle the eigenvalues and vectors
> from system 2 so that its eigenvectors are
> maximally aligned with those of system 1.
>
> This seems like it might be easy, except that
> the tiny eigenvalues for some systems may
> sometimes end up with virtually random
> eigenvectors. I'd visualize some potentially
> nasty problems to resolve.
>

Thanks for the reply!

I've already looked into the eigenvector-approach. It works
most of the time. Since the eigenvectors returned by matlab
always have length 1, you can figure out which eigenvectors
are close together by inspecting the norm of the
difference. Basically, I use

normOfDifference = sum(abs(evectorOld - evectorNew));

as the measure of difference. I then figure out which
eigenvectors are the best match*. The assumption underlying
this however is that:

sum(abs(evector_i - evector_j)) >>> sum(abs(evector_i_now -
evector_i_previous)).

This fails every once in a while though.

Using the angle between the differenct is also an option,
but not based on the cosine; round-off error is too large.

Willem Geert

* Compensating for the fact that eigenvectors are only
defined up unto a scalar multiple, sometimes eig() also
returns them as the negative of the previous eigenvectors.

"Willem Geert " <wgphaff@gmail.com> wrote in message
<fnv8sh$puv$1@fred.mathworks.com>...
>
> > Assume that you have two sets of eigenvectors
> > and eigenvalues, related by the fact that they
> > came from systems which were close to each
> > other. Now shuffle the eigenvalues and vectors
> > from system 2 so that its eigenvectors are
> > maximally aligned with those of system 1.
> >
> > This seems like it might be easy, except that
> > the tiny eigenvalues for some systems may
> > sometimes end up with virtually random
> > eigenvectors. I'd visualize some potentially
> > nasty problems to resolve.
> >
>
> Thanks for the reply!
>
> I've already looked into the eigenvector-approach. It works
> most of the time. Since the eigenvectors returned by matlab
> always have length 1, you can figure out which eigenvectors
> are close together by inspecting the norm of the
> difference. Basically, I use
>
> normOfDifference = sum(abs(evectorOld - evectorNew));
>
> as the measure of difference. I then figure out which
> eigenvectors are the best match*. The assumption underlying
> this however is that:
>
> sum(abs(evector_i - evector_j)) >>> sum(abs(evector_i_now -
> evector_i_previous)).
>
> This fails every once in a while though.
>
> Using the angle between the differenct is also an option,
> but not based on the cosine; round-off error is too large.
>
> Willem Geert
>
> * Compensating for the fact that eigenvectors are only
> defined up unto a scalar multiple, sometimes eig() also
> returns them as the negative of the previous eigenvectors.

This is why the norm of the difference is not
the right method. instead, use the dot product
of the two vectors, and take the absolute value.

Two vectors which are close, even if the sign is
wrong on one of them, will still have an absolute
dot product near 1.

John

"Willem Geert " <wgphaff@gmail.com> wrote in message <fnv8sh$puv
$1@fred.mathworks.com>...
> ...........
> Using the angle between the differenct is also an option,
> but not based on the cosine; round-off error is too large.
> ...........
----------
  Willem, if you are worried about accuracy in determining when vectors are
nearly parallel, there is a more accurate alternative as compared with finding
their scalar product and comparing its absolute value to 1, though it requires
more computation. The scalar product is the product of the norms of the
vectors multiplied by the cosine of the angle between them. As you
indicated, there is a loss of accuracy due to the infinite slope that is
approached by the arccosine function as its argument approaches 1 or -1.
Instead, we can directly compute the absolute value of the sine of the angle
between vectors in the following manner.

  The absolute value of the sine of the angle between two real n-dimensional
(n-element) column vectors v and w is given in matlab by

 norm(t(:))/sqrt(2)/norm(v)/norm(w)

where t is computed as

 t = v*w.'-w*v.';

Of course with unit-length eigenvectors, the division by norm(v) and norm(w)
can be omitted. Unfortunately, finding t and norm(t(:)) are order n^2
computations.

  When this sine quantity is nearly zero, the vectors are nearly parallel, even if
pointing in opposite directions. It should also be noted that this same
quantity approaches zero even when v and w are complex-valued and are
nearly parallel in the complex sense that one is nearly a complex scalar
multiple of the other, as can occur with eigenvectors of non-Hermitian
matrices.

  Note: Just for the record, you can get the actual angle between the real
vectors, lying in the range 0 to pi, using

 ang = atan2(norm(t(:))/sqrt(2),dot(v,w));

Roger Stafford