Re: [eigen] ordering of eigenvalues of EigenSolver

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2011/2/21 Susanne Suter <susanne.suter@xxxxxxxxx>:
> Hi
>
> Thanks Robert for pointing out my question, which was about sorting
> according to absolute value or signed value. Because in the current
> implementation the eigenvalues are ordered taking into account the
> sign (my example: -2.77687, -1.50582, -0.0370092,  0.848378 (last two
> values should be switched or if ascending order, the whole ordering
> should be reversed).

Indeed, both Eigen::SelfAdjointEigenSolver and LAPACK sorts in plain
ascending order, not in ascending order of absolute values.

>
> I'm actually using LAPACK SVD now, but I was considering to use eigen
> library for SVD and/or Eigenvalue decomposition. So, I'm verifying if
> eigen produces the results I expect.

For both SelfAdjointEigenSolver and (Jacobi)SVD, Eigen should have the
same sorting behavior as LAPACK.

>
> I think we should draw a differentiation between SVD and eigenvalue
> decomposition:
>
> For SVD, the solution is unique except for the sign and the singular
> values are always ordered descending. I.e., the singular values are
> ordered according to their absolute value.
>
> To my understanding, for eigenvalue decomposition the values don't
> necessarily have to be ordered.

Indeed: non-self-adjoint matrices may have non-real (complex)
eigenvalues. (Actually, in the case of normal matrices,
self-adjointness is equivalent to all eigenvalues being real). In this
case, there is no natural sorting.

> But I can only think of applications,
> which use them ordered and when ordering them, I would use it
> according to absolute value (for real value problems). I.e., use the x
> largest magnitude values.

Indeed, I understand your point. It is indeed very natural to sort
eigenvalues by absolute value. This is in particular a requirement if
one wants to see SelfAdjointEigenSolver as a special case of SVD.

However, doing this loses uniqueness of the sorting. Consider a
symmetry matrix, so its eigenvalues are +1 and -1 with certain
multiplicities. Since +1 and -1 have the same absolute value, sorting
by absolute value would be non-unique. A matrix with eigenvalue +1
with multiplicity 2 and eigenvalue -1 with multiplicity 3 could give
this eigenvalues vector,
+1
-1
-1
+1
-1

or if could equally well give that one:
-1
+1
+1
-1
-1

The advantage of sorting by plain value is that it guarantees this
unique result:
-1
-1
-1
+1
+1


>
> Furthermore, other eigenvalue decomposition libraries seem to order
> eigenvalues, too. E.g., MatLab (based on ARPACK)

Ah, I didn't know about Matlab.

> orders the eigenvalue
> according to absolute value and it seems that LAPACK SSYEV does the
> same but in reverse order.

No, if that link is accurate,
http://www.netlib.org/lapack/single/ssyev.f
then the eigenvalues are sorted "in ascending order" which, as far as
I understand, does not mean sorted by ascending absolute values.

Cheers,
Benoit

>
> So, is there a bug in the eigen library or are the values intended to
> be sorted like that (according to weight and direction/sign)?
>
> Thank you and best,
> Susanne
>
> On Mon, Feb 21, 2011 at 10:20 AM, Robert Bocquier
> <robert.bocquier@xxxxxxxxxxx> wrote:
>> Hi Benoit,
>>
>> I am interested in your answer to Susanne question, but I think you
>> didn't fully address it.
>> Her question wasn't about ascending versus descending ordering, but was
>> about the basis for the sort. Is it (and should it be) the eingenvalues
>> directly, or the absolute values of the eigenvalues ?
>>
>> Thx
>> Robert
>>
>> Le 18/02/2011 16:59, Benoit Jacob a écrit :
>>> Hi,
>>>
>>> Sorry for the belated answer, I hadn't actually checked carefully what
>>> both Eigen and others (LAPACK) were doing until today.
>>>
>>> As it turns out, we are doing exactly the same thing as LAPACK:
>>>
>>>  * for self-adjoint eigensolver, we sort eigenvalues in ascending
>>> order. For LAPACK, see:
>>>     http://www.netlib.org/lapack/single/ssyev.f
>>>
>>>  * for SVD, we sort singular values in descending order. For SVD, see:
>>>     http://www.netlib.org/lapack/single/sgesvd.f
>>>
>>> I agree that it's pretty weird to be using sometimes ascending and
>>> sometimes descending order. But since that's what both LAPACK and
>>> ourselves have been doing, we shouldn't change that now.
>>>
>>> Benoit
>>>
>>> 2011/2/8 Susanne Suter <susanne.suter@xxxxxxxxx>:
>>>> Hi
>>>>
>>>> Sorry, I accidently hit the "send button" too early. Here my complete
>>>> messge again.
>>>>
>>>> I'm testing the Eigen library in order to use it for eigenvalue
>>>> decomposition or SVD. In the end I need the eigenvectors (left
>>>> singular vectors) and the eigenvalues (singular values).
>>>>
>>>> I noticed that when using the EigenSolver classes, the ordering of the
>>>> eigenvalues is not as I expected. Normally, I would expect that I get
>>>> an ordering analogous to the singular values, i.e., an ordering with
>>>> the maximum absolute value first and then decreasing values ordered by
>>>> their absolute value (sometimes also called "largest magnitude
>>>> eigenvalues"). What I get now, looks a bit like an ordering
>>>> considering the sign (SelfAdjointEigenSolver). I had similar issues
>>>> with EigenSolver and ComplexEigenSolver (however, I found the ordering
>>>> not consistent in all methods).
>>>>
>>>> ### SelfAdjointEigenSolver: eigenvalues of A are:    -2.77687
>>>> -1.50582 -0.0370092   0.848378
>>>> ### ComplexEigenSolver: eigenvalues are: (-0.0370093,0) (0.848378,0)
>>>> (-1.50582,0) (-2.77687,0)
>>>> ### EigenSolver:  eigenvalues are : (-2.77687,0) (-0.0370092,0)
>>>> (0.848379,0)   (-1.50582,0)
>>>>
>>>> However, I would expect
>>>> ### eigenvalues of A are:   -2.77687   -1.50582 0.848378 -0.0370092
>>>>
>>>> Until now, I only was using eigenvalue decompositions, which used the
>>>> same ordering for eigenvalues as the SVD for the singular values (e.g.
>>>> Matlab). Is there any reason why the implementation of Eigen is
>>>> different? Or is there any option to change the ordering?
>>>>
>>>> I'm using the following test matrix:
>>>>
>>>> A =
>>>>
>>>>   -2.0000   -0.6714    0.8698    0.5792
>>>>   -0.6714   -1.1242   -0.0365   -0.5731
>>>>    0.8698   -0.0365   -0.4660   -0.8542
>>>>    0.5792   -0.5731   -0.8542    0.1188
>>>>
>>>> and I test it with the following code (using the Eigen3 version).
>>>>
>>>>        Eigen::SelfAdjointEigenSolver<Eigen::Matrix4f> es;
>>>>        es.compute(A);
>>>>        std::cout << "### SELFADJOINT EIGENVALUE DECOMPOSITION ### " << std::endl;
>>>>        std::cout << "** eigenvalues of A are: " <<
>>>> es.eigenvalues().transpose() << std::endl;
>>>>        std::cout << "** eigenvectors of A are: " << std::endl <<
>>>> es.eigenvectors() << std::endl;
>>>>
>>>> Thank you for any hints or advice on that topic.
>>>>
>>>> Best,
>>>> Susanne
>>>>
>>>>
>>>>
>>>
>>>
>>
>>
>>
>>
>
>
>