[SciPy-dev] Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

[Charles R Harris]

On Thu, Feb 4, 2010 at 8:08 PM, <josef.pktd at gmail.com> wrote:

> On Thu, Feb 4, 2010 at 9:54 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> >
> > On Thu, Feb 4, 2010 at 7:24 PM, <josef.pktd at gmail.com> wrote:
> >>
> >> On Thu, Feb 4, 2010 at 8:27 PM, David Cournapeau <cournape at gmail.com>
> >> wrote:
> >> > On Fri, Feb 5, 2010 at 9:43 AM, Warren Weckesser
> >> > <warren.weckesser at enthought.com> wrote:
> >> >
> >> >> David C,
> >> >>
> >> >> You said this was for symmetric matrices, but do you envision later
> >> >> allowing nonsymmetric matrices?
> >> >
> >> > Yes. I have only implemented symmetric/Hermitian because that's the
> >> > only solver that handles getting only a few eigenvalues in LAPACK.
> >> > Matlab own eigs function uses ARPACK for both symmetric/unsymmetric
> >> > cases, which is not good according to one of the Lapack developer
> >> > (http://mail.scipy.org/pipermail/scipy-dev/2007-March/006778.html).
> >> > But we could use ARPACK for the general case in eigs.
> >> >
> >> >>
> >> >> If not, then perhaps the name of the function should be 'eigsh',
> >> >> following the precedent set by numpy.linalg and scipy.linalg.
> >> >
> >> > I wonder whether those eigh functions are a good idea: I fear that
> >> > most people will always use eig - maybe one could use the underlying
> >> > eig*h solver in eig* if the matrix is detected as being symmetric ? I
> >> > am not really knowledgeable about those issues, though. For example, I
> >> > don't know whether the symmetric aspect should be checked exactly, or
> >> > if it is better to use a symmetric solver  even if there are very
> >> > small errors in say A'-A.
> >> >
> >> >>
> >> >> In particular, the choice of ordering by magnitude or by real part is
> >> >> convenient.
> >> >
> >> > It seems that one would need to implement the non-symmetric
> >> > capabilities to sort this out. I fear that those options are solver
> >> > specific, though - maybe the solution is to have two levels of API,
> >> > one low-level and one high level.
> >> >
> >> > cheers,
> >> >
> >> > David
> >> > _______________________________________________
> >> > SciPy-Dev mailing list
> >> > SciPy-Dev at scipy.org
> >> > http://mail.scipy.org/mailman/listinfo/scipy-dev
> >> >
> >>
> >> the current version of scipy.linalg.eigh sorts from smallest to largest
> >>
> >> >>> scipy.linalg.eigh(x)[0]
> >> array([  4.04457427e-03,   1.84073286e-01,   6.74875960e-01,
> >>         3.23328824e+00,   4.00741304e+00,   6.98333680e+00,
> >>         1.45314842e+01,   1.49260377e+01,   2.47166702e+01,
> >>         2.98955755e+01])
> >> >>> scipy.linalg.eigh(x, eigvals=(0,3))[0]
> >> array([ 0.00404457,  0.18407329,  0.67487596,  3.23328824])
> >> >>> scipy.linalg.eigh(x, eigvals=(len(x)-3, len(x)-1))[0]
> >> array([ 14.92603769,  24.71667021,  29.89557547])
> >>
> >>
> >> >>> x = np.random.randn(10, 10) + 1j*np.random.randn(10, 10)
> >> >>> x = np.dot(x.T, x)
> >> >>> scipy.linalg.eigh(x, eigvals=(0,3))[0]
> >> array([-36.82039662, -18.20362967, -12.21716065,  -9.79752274])
> >> >>> scipy.linalg.eigh(x, eigvals=(len(x)-3, len(x)-1))[0]
> >> array([ 14.44815917,  28.8394026 ,  42.45471058])
> >>
> >
> > Yeah. It's my fault ;) It didn't use to sort at all.
>
> I thought that's how Tiziano added them, I didn't see a sort when I
> briefly browsed the source.
> But I'm looking at lapack only from the outside.
>
>
I was thinking of the numpy version. But now that I think about it, it was a
test that needed the eigenvalues sorted because the routine didn't sort
them. I seem to be developing one of those trick memories...

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20100204/79b81c78/attachment.html>