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 > MessageID: <[hidden email]> > ContentType: text/plain; charset=UTF8; 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/octavemaintainers/2011July/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 eigsbase.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? 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 
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 
In reply to this post by Rik4
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 >> MessageID:<[hidden email]> >> ContentType: text/plain; charset=UTF8; 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/octavemaintainers/2011July/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 eigsbase.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 
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 N2 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 N1 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 
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 N2 > 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 N1 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 nonsparse 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 
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: 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 
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 N2 > 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 
Free forum by Nabble  Edit this page 