[Numpy-discussion] missing array type

Christopher Barker Chris.Barker at noaa.gov
Tue Feb 28 09:28:07 EST 2006

Travis Oliphant wrote:
> [Soap Box]
> I've been annoyed for several years that the idea of linear operators is 
> constrained in most libraries to 2 dimensions.  There are many times I 
> want to find an inverse of an operator that is most naturally expressed 
> with 6 dimensions.

Yes, yes, yes!

"numpy is not matlab"

One of the things I love most about numpy is that it is an n-d array 
package, NOT a matrix package. I also love broadcasting. Similar to 
Travis, I was recently helping out a friend using Matlab for a graduate 
structural mechanics course. The machinations required to shoehorn the 
natural tensor math into 2-d matrices was pretty ugly indeed. I'd much 
rather see numpy encourage the use of higher dimension arrays and 
broadcasting over traditional 2-d matrix solutions.


> I have to myself play games with indexing to give 
> the computer a matrix it can understand.  Why is that?

One of the reasons is that we want to use other people already optimized 
code (i.e. LAPACK). They only work with the 2-d data structures. I 
suppose we could do the translation to the LAPACK data structures under 
the hood, but that would take some work.

However, this makes me wonder....

I'm unclear on the details, but from what I understand of the post that 
started this thread, one use repmat is used in order to turn some 
operations into standard linear algebra operations, and that's done for 
performance purposes. The repmat matrix would therefore need to be in a 
form usable by LAPACK and friends, and thus would need to be dense 
anyway ... a zero-stride array would not work, so maybe the potential 
advantages of the compact storage wouldn't really be realized (until we 
write out own LAPACK)

This also brings me to...

Sasha wrote:
> Desired:
>>>> x = zeros(5)
>>>> x.strides=0
>>>> x += 1
>>>> x
> array([1, 1, 1, 1, 1])
>>>> x += arange(5)
>>>> x
> array([1, 2, 3, 4, 5])

So what the heck is a zero-strided array? My understanding was that the 
whole point was the what looked like multiple values, were really a 
single, shared, value. In this case, is shouldn't be possible to 
in-place add more than one value. I wouldn't say that what Sasha 
presented as "desired" is desired.. an in=place operation shouldn't 
fundamentally change the nature of the array. That array should ALWAYS 
remain single-valued.

So what should the result of x += arange(5) be? I say it should raise an 

Maybe zero-stride arrays are only really useful read-only?

This is a complicated can of worms.....


Christopher Barker, Ph.D.
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov

More information about the NumPy-Discussion mailing list