Re: eigs compatibility

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

Rik-4
On 07/03/2013 12:08 PM, [hidden email] wrote:

> Message: 7
> Date: Wed, 03 Jul 2013 15:08:04 -0400
> From: "John W. Eaton" <[hidden email]>
> To: Jordi Guti?rrez Hermoso <[hidden email]>
> Cc: David Bateman <[hidden email]>, octave maintainers mailing list
> <[hidden email]>
> Subject: Re: eigs fails for scalars and 2x2 matrices
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
> On 06/26/2013 09:59 AM, Jordi Guti?rrez Hermoso wrote:
>> > On 26 June 2013 04:28, John W. Eaton<[hidden email]>  wrote:
>>> >> eigs fails for scalars and 2x2 matrices, apparently intentionally
>>> >> given the error message.  Is this a limitation of ARPACK?  Should we
>>> >> work around this problem by calling eig for these cases?
>> >
>> > This was Bateman's response:
>> >
>> >      https://mailman.cae.wisc.edu/pipermail/octave-maintainers/2011-July/024272.html
> OK, to make the discussion easier to follow, I'm reproducing that
> message here and I'm adding David to the Cc: line since he wrote
> the eigs function in Octave.
>
>    On 07/24/2011 09:10 PM, Jordi Guti?rrez Hermoso wrote:
>    > While porting Matlab code, I noticed that someone had used eigs for a
>    > 2x2 matrix, which produces an error to use eig instead. I patched the
>    > code to do that, but I inadverently introduced a bug because I forgot
>    > that eigs sorts its eigenvalues but eig does not.
>    >
>    > Should eigs be patched to call eig for small matrices and just sort
>    > them? Also, for the case of when you need all eigenvalues, should eigs
>    > just silently call eig instead of warning about it?
>    >
>    > - Jordi G. H.
>    >
>
>    The sorting rules aren't that obvious as the order of the values
>    returned might be largest eignevalue first, smallest first or even
>    sorted relative to a distance from a particular eigenvalue. So if matlab
>    does this it might be better to relax the error in eigs-base.cc to a
>    warning and letter ARPACK do the work.
>
>    David
>
> It does seem that eigs in Matlab works for small matrices.
>
> I'm not sure that it will work to change the error to a warning and
> ask ARPACK to do the job because ARPACK doesn't seem to be able to
> return all eigenvalues, and in Matlab, eigs (rand (2)) returns 2
> eigenvalues.  It also seems that it only starts warning when the size
> of the matrix is greater than 8.
>
> I agree that getting all the sorting options right will take some
> effort.  It seems like the easiest way to deal with this would be to
> write a wrapper eigs.m function that checks the size of the input
> matrix (or matrices) and the number of requested eigenvalues and then
> calls the current eigs function if possible, otherwise issues a
> warning (if appropriate) and calls eig + sort as directed by the
> options.
>
> Does that seem like a reasonable solution?  If so, I'll take a shot at
> doing it.  If not, do you see a better way?
Seems like a low value compatibility item.  'eigs' is supposed to be for
determining a few eigenvectors of a very large, sparse matrix.  Asking it
for the eigenvalues of a 2x2 matrix is attempting to use a sledge hammer to
swat a fly.  I would be more inclined to just check the dimensions of the
input and if it is less than 10 print a warning suggesting that the user
consider 'eig' for their problem.  However, if someone is paying for this
compatibility fix then your method seems reasonable.

--Rik

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

John W. Eaton
Administrator
On 07/03/2013 05:17 PM, Rik wrote:

 > However, if someone is paying for this
> compatibility fix then your method seems reasonable.

I would probably not be attempting to fix it otherwise.

jwe


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

Daniel Sebald
In reply to this post by Rik-4
On 07/03/2013 04:17 PM, Rik wrote:

> On 07/03/2013 12:08 PM, [hidden email] wrote:
>> Message: 7
>> Date: Wed, 03 Jul 2013 15:08:04 -0400
>> From: "John W. Eaton"<[hidden email]>
>> To: Jordi Guti?rrez Hermoso<[hidden email]>
>> Cc: David Bateman<[hidden email]>, octave maintainers mailing list
>> <[hidden email]>
>> Subject: Re: eigs fails for scalars and 2x2 matrices
>> Message-ID:<[hidden email]>
>> Content-Type: text/plain; charset=UTF-8; format=flowed
>>
>> On 06/26/2013 09:59 AM, Jordi Guti?rrez Hermoso wrote:
>>>> On 26 June 2013 04:28, John W. Eaton<[hidden email]>   wrote:
>>>>>> eigs fails for scalars and 2x2 matrices, apparently intentionally
>>>>>> given the error message.  Is this a limitation of ARPACK?  Should we
>>>>>> work around this problem by calling eig for these cases?
>>>>
>>>> This was Bateman's response:
>>>>
>>>>       https://mailman.cae.wisc.edu/pipermail/octave-maintainers/2011-July/024272.html
>> OK, to make the discussion easier to follow, I'm reproducing that
>> message here and I'm adding David to the Cc: line since he wrote
>> the eigs function in Octave.
>>
>>     On 07/24/2011 09:10 PM, Jordi Guti?rrez Hermoso wrote:
>>     >  While porting Matlab code, I noticed that someone had used eigs for a
>>     >  2x2 matrix, which produces an error to use eig instead. I patched the
>>     >  code to do that, but I inadverently introduced a bug because I forgot
>>     >  that eigs sorts its eigenvalues but eig does not.
>>     >
>>     >  Should eigs be patched to call eig for small matrices and just sort
>>     >  them? Also, for the case of when you need all eigenvalues, should eigs
>>     >  just silently call eig instead of warning about it?
>>     >
>>     >  - Jordi G. H.
>>     >
>>
>>     The sorting rules aren't that obvious as the order of the values
>>     returned might be largest eignevalue first, smallest first or even
>>     sorted relative to a distance from a particular eigenvalue. So if matlab
>>     does this it might be better to relax the error in eigs-base.cc to a
>>     warning and letter ARPACK do the work.
>>
>>     David
>>
>> It does seem that eigs in Matlab works for small matrices.
>>
>> I'm not sure that it will work to change the error to a warning and
>> ask ARPACK to do the job because ARPACK doesn't seem to be able to
>> return all eigenvalues, and in Matlab, eigs (rand (2)) returns 2
>> eigenvalues.  It also seems that it only starts warning when the size
>> of the matrix is greater than 8.
>>
>> I agree that getting all the sorting options right will take some
>> effort.  It seems like the easiest way to deal with this would be to
>> write a wrapper eigs.m function that checks the size of the input
>> matrix (or matrices) and the number of requested eigenvalues and then
>> calls the current eigs function if possible, otherwise issues a
>> warning (if appropriate) and calls eig + sort as directed by the
>> options.
>>
>> Does that seem like a reasonable solution?  If so, I'll take a shot at
>> doing it.  If not, do you see a better way?
> Seems like a low value compatibility item.  'eigs' is supposed to be for
> determining a few eigenvectors of a very large, sparse matrix.  Asking it
> for the eigenvalues of a 2x2 matrix is attempting to use a sledge hammer to
> swat a fly.  I would be more inclined to just check the dimensions of the
> input and if it is less than 10 print a warning suggesting that the user
> consider 'eig' for their problem.  However, if someone is paying for this
> compatibility fix then your method seems reasonable.

I guess a sparse matrix of 2x2 isn't really all that sparse.  Anyway, if
it is important to have compatibility or just handle all dimensions, and
it is only 2x2 and 1x1 that must be dealt with, rather than require a
different math package be used to offer those results, it might be just
as easy to write out the formulas for the solution to the eigenvalues.
If one were to study the math, it probably breaks down into some simple
inversion formulas which aren't so challenging until getting above 3x3.
  Why isn't there support in the sparse math package?

Dan
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

John W. Eaton
Administrator
On 07/04/2013 03:03 PM, Daniel J Sebald wrote:

> I guess a sparse matrix of 2x2 isn't really all that sparse. Anyway, if
> it is important to have compatibility or just handle all dimensions, and
> it is only 2x2 and 1x1 that must be dealt with, rather than require a
> different math package be used to offer those results, it might be just
> as easy to write out the formulas for the solution to the eigenvalues.
> If one were to study the math, it probably breaks down into some simple
> inversion formulas which aren't so challenging until getting above 3x3.
> Why isn't there support in the sparse math package?

I have no clue why ARPACK has the limitation of providing at most N-2
eigenvalues, but that's the way it works.  I'm not going to try to fix
that problem.

Apparently Matlab's eigs function automatically handles the cases of
computing N or N-1 eigenvalues.  I don't know how to do that other than
by calling eig and sorting the results as needed.

I think Matlab does not issue a warning until N >= 8, then it says that
you should use "eig (full (A))" or something similar, but still goes
ahead and does the calculation anyway.  Would someone like to confirm
that behavior?

I don't think it makes sense to try to write out the analyical solution
for eigenvalues for N == 6, or 5, ...

I'm also not trying to defend this behavior.  I just want to make it
work because a paying customer wants it to work.

jwe
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

Daniel Sebald
On 07/04/2013 02:21 PM, John W. Eaton wrote:

> On 07/04/2013 03:03 PM, Daniel J Sebald wrote:
>
>> I guess a sparse matrix of 2x2 isn't really all that sparse. Anyway, if
>> it is important to have compatibility or just handle all dimensions, and
>> it is only 2x2 and 1x1 that must be dealt with, rather than require a
>> different math package be used to offer those results, it might be just
>> as easy to write out the formulas for the solution to the eigenvalues.
>> If one were to study the math, it probably breaks down into some simple
>> inversion formulas which aren't so challenging until getting above 3x3.
>> Why isn't there support in the sparse math package?
>
> I have no clue why ARPACK has the limitation of providing at most N-2
> eigenvalues, but that's the way it works. I'm not going to try to fix
> that problem.

I'm not familiar with the ARPACK library, but it looks like the number
of subdimensions is selectable as 'k'.  I'm going to guess that these
eigenvalue techniques of ARPACK are fast at finding the eigenvalues of
higher strength but then become slower with weaker eigenvalues.  So, if
one were to ask for all the eigenvalues, the routine my get bogged down,
in a sense.  Not sure.

But here

http://www.caam.rice.edu/software/ARPACK/UG/node13.html#SECTION00570000000000000000

the developers indicate that BLAS and LAPACK are requirements and if
there is something in those libraries to substitute or utilize, do so.
Given that, I guess falling back on those libraries for small matrix
eigenvalues is fine.


> Apparently Matlab's eigs function automatically handles the cases of
> computing N or N-1 eigenvalues. I don't know how to do that other than
> by calling eig and sorting the results as needed.

Oh, so it isn't just the 2x2 case.  It might be if someone has a 20x20
matrix and asks for 19 eigenvalues/vectors?  I would say try a timing
experiment comparing:

20x20 non-sparse
20x20 sparse requesting 18 eigenvalues/vectors

If those two are close then, i.e., the ARPACK really slows down as the
number of subdimensions approaches the dimension of the matrix, then the
ARPACK developers probably figured it is best to have the users of the
software fall back on the LAPACK or BLAS full matrix equivalent...as you
are proposing.  I've a feeling that there really isn't a problem here,
it just that the original Octave programmer didn't realize ARPACK
suggests falling back on LAPACK or BLAS.


> I think Matlab does not issue a warning until N >= 8, then it says that
> you should use "eig (full (A))" or something similar, but still goes
> ahead and does the calculation anyway. Would someone like to confirm
> that behavior?
>
> I don't think it makes sense to try to write out the analyical solution
> for eigenvalues for N == 6, or 5, ...

One wouldn't want to write out analytical solutions for N equal 5 or 6,
just 2 or in a pinch 3.  Anything above that is best handled with
rotation, pivoting, whatever.

Also, I don't know why Matlab would issue a warning about N > 7.  That
doesn't seem like a such a big and CPU consuming problem.  I would guess
somewhere above 20 is when things become noticeably slow for modern
processors.  Maybe 15, but certainly not 7.

Dan
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

eem2314


On Thu, Jul 4, 2013 at 12:54 PM, Daniel J Sebald <[hidden email]> wrote:
On 07/04/2013 02:21 PM, John W. Eaton wrote:
On 07/04/2013 03:03 PM, Daniel J Sebald wrote:

I have no clue why ARPACK has the limitation of providing at most N-2
eigenvalues, but that's the way it works. I'm not going to try to fix
that problem.

they use a shifting strategy to get eigenvalues clustered nearest a given shift, so
if you want all the eigenvalues for a large sparse matrix you can call arpack twice,
once for k/2 eigenvalues at the low end of the spectrum, and again for the k/2 at
the high end. "Large" is greater than 300th order even on my klunky laptop; otherwise
it makes more sense to just call the lapack routines if the user wants a large percentage
of the eigenvalues.


--
Ed Meyer
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: eigs compatibility

Daniel Sebald
On 07/04/2013 05:52 PM, Ed Meyer wrote:

>
>
> On Thu, Jul 4, 2013 at 12:54 PM, Daniel J Sebald <[hidden email]
> <mailto:[hidden email]>> wrote:
>
>     On 07/04/2013 02:21 PM, John W. Eaton wrote:
>
>         On 07/04/2013 03:03 PM, Daniel J Sebald wrote:
>
>         I have no clue why ARPACK has the limitation of providing at
>         most N-2
>         eigenvalues, but that's the way it works. I'm not going to try
>         to fix
>         that problem.
>
>
> they use a shifting strategy to get eigenvalues clustered nearest a
> given shift, so
> if you want all the eigenvalues for a large sparse matrix you can call
> arpack twice,
> once for k/2 eigenvalues at the low end of the spectrum, and again for
> the k/2 at
> the high end.

Creative idea.  I might try a time comparison.


>  "Large" is greater than 300th order even on my klunky
> laptop; otherwise
> it makes more sense to just call the lapack routines if the user wants a
> large percentage
> of the eigenvalues.

Thanks Ed.

Dan
Loading...