Developers of Scientific Python

Re: Ball Tree code updated (ticket 1048)

Jake VanderPlas wrote:
> Hello,
> I have had comments from a few people over the last two months on the
> Ball Tree code that I submitted (ticket 1048).  I cleaned up the code
> a bit and posted the changes on the tracker.  Any other comments would
> be appreciated!
>    -Jake
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev <at> scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
> 

Hi! Your code works great, with impressive speed improvements in high 
dimensions. I could just find a minor thing: leaf_size in the docstring 
should read leafsize.

If this turns out to be too specific for scipy, would you mind that it 
get's included in scikit.learn [1] (BSD-license)?

Thanks,

fabian

[1] http://scikit-learn.sourceforge.net

headers

David Cournapeau | 4 Feb 07:27

Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

Hi,

I have played a bit with adding support for computing a few eigenvalues 
of full symmetric matrices. I would like some comments on the current API:

import numpy as np
from scipy.linalg import eigs
x = np.random.randn(10, 10)
x = np.dot(x.T, x)
# Retrieve the 3 biggest eigenvalues
eigs(x, 3)[0]
# Retrieve the 3 smallest eigenvalues
eigs(x, -3)[0]
# Retrieve the 3 biggest eigenvalues
eigs(x, [0, 3], mode="index")[0]
# Retrieve the 2nd and 3rd biggest
eigs(x, [1, 3], mode="index")[0]
# Retrieve all the eigenvalues in the range [1.5, 3.5[
eigs(x, [1.5, 3.5], mode="range")[0]

One thing which does not feel right is that that the range in the 
"index" mode is exactly inverted compared to the output (i.e. if you ask 
for the range [0, 3], you get the last three items from what you would 
get if you asked for the full range [0, 10]), but this is because I kept 
compatibility with Octave (always showing from biggest to smallest). It 
would be easy to always do the contrary (from smallest to biggest) - 
besides consistency, it has the advantages that it is the actual order 
you get back from the underlying LAPACK function.

Also, I needed to modify the f2py file for the related LAPACK functions,

(Continue reading)

headers

Charles R Harris | 4 Feb 07:37

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

On Wed, Feb 3, 2010 at 11:27 PM, David Cournapeau <david <at> silveregg.co.jp> wrote:

Hi,

I have played a bit with adding support for computing a few eigenvalues
of full symmetric matrices. I would like some comments on the current API:

import numpy as np
from scipy.linalg import eigs
x = np.random.randn(10, 10)
x = np.dot(x.T, x)
# Retrieve the 3 biggest eigenvalues
eigs(x, 3)[0]
# Retrieve the 3 smallest eigenvalues
eigs(x, -3)[0]
# Retrieve the 3 biggest eigenvalues
eigs(x, [0, 3], mode="index")[0]
# Retrieve the 2nd and 3rd biggest
eigs(x, [1, 3], mode="index")[0]
# Retrieve all the eigenvalues in the range [1.5, 3.5[
eigs(x, [1.5, 3.5], mode="range")[0]

One thing which does not feel right is that that the range in the
"index" mode is exactly inverted compared to the output (i.e. if you ask

Why not use separate keywords for index and range. Changing the meaning of an argument using another keyword is just weird. I've seen it used elsewhere, but still... At some point you might even want the largest three in a range. So something like

eigs(x, index=3)
eigs(x, range=[1,4])

Why the index on the end? Is the return a list?

for the range [0, 3], you get the last three items from what you would
get if you asked for the full range [0, 10]), but this is because I kept
compatibility with Octave (always showing from biggest to smallest). It

Why follow Octave, isn't Octave like Matlab? Matlab functions are a mess.

Chuck

_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

headers

David Cournapeau | 4 Feb 07:49

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

Charles R Harris wrote:
> 
> 
> On Wed, Feb 3, 2010 at 11:27 PM, David Cournapeau <david <at> silveregg.co.jp 
> <mailto:david <at> silveregg.co.jp>> wrote:
> 
>     Hi,
> 
>     I have played a bit with adding support for computing a few eigenvalues
>     of full symmetric matrices. I would like some comments on the
>     current API:
> 
>     import numpy as np
>     from scipy.linalg import eigs
>     x = np.random.randn(10, 10)
>     x = np.dot(x.T, x)
>     # Retrieve the 3 biggest eigenvalues
>     eigs(x, 3)[0]
>     # Retrieve the 3 smallest eigenvalues
>     eigs(x, -3)[0]
>     # Retrieve the 3 biggest eigenvalues
>     eigs(x, [0, 3], mode="index")[0]
>     # Retrieve the 2nd and 3rd biggest
>     eigs(x, [1, 3], mode="index")[0]
>     # Retrieve all the eigenvalues in the range [1.5, 3.5[
>     eigs(x, [1.5, 3.5], mode="range")[0]
> 
>     One thing which does not feel right is that that the range in the
>     "index" mode is exactly inverted compared to the output (i.e. if you ask
> 
> 
> Why not use separate keywords for index and range.

Then I am not sure how to handle the "give me the k biggest/smallest" 
case, which is the most common I think (that's the only one I care 
personally :) ).

Maybe this just warrants several functions ? But then there is the issue 
of naming (supporting non-symmetric matrices would be nice, but requires 
a totally different implementation, as LAPACK does not support it AFAIK 
- the easiest way would be to use ARPACK ATM).

> Changing the meaning 
> of an argument using another keyword is just weird. 

Agreed - this is just the best I could came up with one function to 
support all cases while keeping the common case simple. But I feel a bit 
like I outsmarted myself here.

> Why the index on the end? Is the return a list?

Yes, it also returns eigenvectors

> 
>     for the range [0, 3], you get the last three items from what you would
>     get if you asked for the full range [0, 10]), but this is because I kept
>     compatibility with Octave (always showing from biggest to smallest). It
> 
> 
> Why follow Octave, isn't Octave like Matlab? Matlab functions are a mess.

Actually, there is another reason I forgot to mention: that's how eigen 
in scipy.sparse does as well (in decreasing order). I think that the 
ARPACK wrappers are not that well thought out, though.

cheers,

David

headers

Charles R Harris | 4 Feb 08:06

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

On Wed, Feb 3, 2010 at 11:49 PM, David Cournapeau <david <at> silveregg.co.jp> wrote:

Charles R Harris wrote:
>
>
> On Wed, Feb 3, 2010 at 11:27 PM, David Cournapeau <david <at> silveregg.co.jp

> <mailto:david <at> silveregg.co.jp>> wrote:
>
> Hi,
>
> I have played a bit with adding support for computing a few eigenvalues
> of full symmetric matrices. I would like some comments on the
> current API:
>
> import numpy as np
> from scipy.linalg import eigs
> x = np.random.randn(10, 10)
> x = np.dot(x.T, x)
> # Retrieve the 3 biggest eigenvalues
> eigs(x, 3)[0]
> # Retrieve the 3 smallest eigenvalues
> eigs(x, -3)[0]
> # Retrieve the 3 biggest eigenvalues
> eigs(x, [0, 3], mode="index")[0]
> # Retrieve the 2nd and 3rd biggest
> eigs(x, [1, 3], mode="index")[0]
> # Retrieve all the eigenvalues in the range [1.5, 3.5[
> eigs(x, [1.5, 3.5], mode="range")[0]
>
> One thing which does not feel right is that that the range in the
> "index" mode is exactly inverted compared to the output (i.e. if you ask
>
>
> Why not use separate keywords for index and range.

Then I am not sure how to handle the "give me the k biggest/smallest"
case, which is the most common I think (that's the only one I care
personally :) ).

That would be eigs(x, index=3) or eigs(x, index=-3) respectively, the default value of both range and index would be None, which could possibly return all eigenvalues. I'm not sure that index is the best word but I can't think of a better at the moment.

<snip>

Chuck

_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

headers

Dag Sverre Seljebotn | 4 Feb 09:12

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

David Cournapeau wrote:
> Charles R Harris wrote:
>>
>> On Wed, Feb 3, 2010 at 11:27 PM, David Cournapeau <david <at> silveregg.co.jp 
>> <mailto:david <at> silveregg.co.jp>> wrote:
>>
>>     Hi,
>>
>>     I have played a bit with adding support for computing a few eigenvalues
>>     of full symmetric matrices. I would like some comments on the
>>     current API:
>>
>>     import numpy as np
>>     from scipy.linalg import eigs
>>     x = np.random.randn(10, 10)
>>     x = np.dot(x.T, x)
>>     # Retrieve the 3 biggest eigenvalues
>>     eigs(x, 3)[0]
>>     # Retrieve the 3 smallest eigenvalues
>>     eigs(x, -3)[0]
>>     # Retrieve the 3 biggest eigenvalues
>>     eigs(x, [0, 3], mode="index")[0]
>>     # Retrieve the 2nd and 3rd biggest
>>     eigs(x, [1, 3], mode="index")[0]
>>     # Retrieve all the eigenvalues in the range [1.5, 3.5[
>>     eigs(x, [1.5, 3.5], mode="range")[0]
>>
>>     One thing which does not feel right is that that the range in the
>>     "index" mode is exactly inverted compared to the output (i.e. if you ask
>>
>>
>> Why not use separate keywords for index and range.
> 
> Then I am not sure how to handle the "give me the k biggest/smallest" 
> case, which is the most common I think (that's the only one I care 
> personally :) ).
> 
> Maybe this just warrants several functions ? But then there is the issue 
> of naming (supporting non-symmetric matrices would be nice, but requires 
> a totally different implementation, as LAPACK does not support it AFAIK 
> - the easiest way would be to use ARPACK ATM).
> 
>> Changing the meaning 
>> of an argument using another keyword is just weird. 
> 
> Agreed - this is just the best I could came up with one function to 
> support all cases while keeping the common case simple. But I feel a bit 
> like I outsmarted myself here.
> 
>> Why the index on the end? Is the return a list?
> 
> Yes, it also returns eigenvectors
> 
>>     for the range [0, 3], you get the last three items from what you would
>>     get if you asked for the full range [0, 10]), but this is because I kept
>>     compatibility with Octave (always showing from biggest to smallest). It
>>
>>
>> Why follow Octave, isn't Octave like Matlab? Matlab functions are a mess.
> 
> Actually, there is another reason I forgot to mention: that's how eigen 
> in scipy.sparse does as well (in decreasing order). I think that the 
> ARPACK wrappers are not that well thought out, though.

I feel there's also a certain convention of listing the largest 
eigenvalues first; "let \lambda_i be the eigenvalues...assume \lambda_2 
< \lambda_1 ...".

I know I'd assume that the largest eigenvalue came out first, and I've 
barely used either MATLAB or Octave.

(Also with SVD one would typically put the largest singular values 
first, and so on.)

Just a data-point, I don't really care either way...

--

-- 
Dag Sverre

headers

Sturla Molden | 4 Feb 11:40

Re: compilation with fort77


>>
>>
>>
>> the problem is that fort77 cannot deal with those variable sized
>> input arguments:
>>
>>      SUBROUTINE mvnun(d, n, lower, upper, means, covar, maxpts,
>>     &                   abseps, releps, value, inform)
>> ...
>>      integer n, d, infin(d), maxpts, inform, tmpinf
>>      double precision lower(d), upper(d), releps, abseps,
>>     &                 error, value, stdev(d), rho(d*(d-1)/2),
>>     &                 covar(d,d),
>>     &                 nlower(d), nupper(d), means(d,n), tmpval
>>      integer i, j
>>
>> Could some Fortran expert please help me to make it fort77  
>> compatible?

The code posted here (I haven't looked in SVN) looks like valid  
Fortran 77 to me. Dummy arguments can have variable size. f2c might  
not accept it because variable-size arrays are not supported in C89,  
and the Fortran to C conversion is rather primitive.

Sturla

>
> This file is using some F90/F95 - I am not sure what the policy is on
> this point, but don't think we want to guarantee that we will never  
> use
> fortran > F77.
>
> Can't you get gfortran running on your platform ? Cross-compiling
> gfortran  for your platform should not be difficult
>
> David
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev <at> scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev

headers

Pauli Virtanen | 4 Feb 12:30

Re: compilation with fort77

Thu, 04 Feb 2010 11:40:52 +0100, Sturla Molden wrote:
>>> the problem is that fort77 cannot deal with those variable sized input
>>> arguments:
>>>
>>>      SUBROUTINE mvnun(d, n, lower, upper, means, covar, maxpts,
>>>     &                   abseps, releps, value, inform)
>>> ...
>>>      integer n, d, infin(d), maxpts, inform, tmpinf double precision
>>>      lower(d), upper(d), releps, abseps,
>>>     &                 error, value, stdev(d), rho(d*(d-1)/2), &       
>>>              covar(d,d),
>>>     &                 nlower(d), nupper(d), means(d,n), tmpval
>>>      integer i, j
>>>
>>> Could some Fortran expert please help me to make it fort77 compatible?
> 
> The code posted here (I haven't looked in SVN) looks like valid Fortran
> 77 to me. Dummy arguments can have variable size. f2c might not accept
> it because variable-size arrays are not supported in C89, and the
> Fortran to C conversion is rather primitive.

Note that infin, stdev, rho, nlower, and nupper are not dummy variables, 
so it's not valid F77.

--

-- 
Pauli Virtanen

headers

David Warde-Farley | 4 Feb 20:31

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

On 4-Feb-10, at 3:12 AM, Dag Sverre Seljebotn wrote:

> I feel there's also a certain convention of listing the largest
> eigenvalues first; "let \lambda_i be the eigenvalues...assume  
> \lambda_2
> < \lambda_1 ...".

+1, though I've seen them listed from smallest to largest as well.

David

headers

Rob Clewley | 4 Feb 21:01

Re: Comments on API for Matlab's eigs equivalent (computing a few eigenvalues only)

On Thu, Feb 4, 2010 at 2:31 PM, David Warde-Farley <dwf <at> cs.toronto.edu> wrote:
> On 4-Feb-10, at 3:12 AM, Dag Sverre Seljebotn wrote:
>
>> I feel there's also a certain convention of listing the largest
>> eigenvalues first; "let \lambda_i be the eigenvalues...assume
>> \lambda_2
>> < \lambda_1 ...".
>
> +1, though I've seen them listed from smallest to largest as well.
>
> David

I'm +1 too, but are these going to be ordered by norm of the
eigenvalue if they are complex, or by the real part only? I suggest by
norm, and the docstring needs to make the choice clear. This
functionality has to be consistent for all uses of eig()!

-Rob

Mon	Tue	Wed	Thu	Fri	Sat	Sun

1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28

88	February 2012
196	January 2012
65	December 2011
44	November 2011
85	October 2011
122	September 2011
95	August 2011
33	July 2011
67	June 2011
82	May 2011
61	April 2011
62	March 2011
128	February 2011
99	January 2011
52	December 2010
123	November 2010
62	October 2010
92	September 2010
85	August 2010
335	July 2010
613	June 2010
207	May 2010
194	April 2010
101	March 2010
202	February 2010
113	January 2010
201	December 2009
178	November 2009
301	October 2009
166	September 2009
328	August 2009
162	July 2009
215	June 2009
183	May 2009
177	April 2009
251	March 2009
523	February 2009
227	January 2009
306	December 2008
417	November 2008
174	October 2008
116	September 2008
158	August 2008
130	July 2008
213	June 2008
160	May 2008
333	April 2008
170	March 2008
145	February 2008
200	January 2008
152	December 2007
42	November 2007
116	October 2007
124	September 2007
126	August 2007
192	July 2007
171	June 2007
180	May 2007
197	April 2007
130	March 2007
66	February 2007
116	January 2007
181	December 2006
51	November 2006
151	October 2006
64	September 2006
45	August 2006
71	July 2006
85	June 2006
45	May 2006
194	April 2006
230	March 2006
159	February 2006
606	January 2006
358	December 2005
477	November 2005
536	October 2005
118	September 2005
51	August 2005
45	July 2005
38	June 2005
51	May 2005
26	April 2005
52	March 2005
55	February 2005
66	January 2005
47	December 2004
254	November 2004
143	October 2004
94	September 2004
21	August 2004
12	July 2004
20	June 2004
28	May 2004
50	April 2004
43	March 2004
9	February 2004
16	January 2004
22	December 2003
42	November 2003
77	October 2003
46	September 2003
22	August 2003
14	July 2003
10	June 2003
48	May 2003
23	April 2003
99	March 2003
39	February 2003
34	January 2003
56	December 2002
64	November 2002
89	October 2002
95	September 2002
76	August 2002
16	July 2002
24	June 2002
31	May 2002
1	January 2000