From: eric <er...@en...> - 2002-03-07 07:26:28
|
While we're discussing several big issues, it would be worthwhile to step back and get an idea of what people feel is still missing for numeric computation in the Python language and also what is missing from Numeric itself. I'm not talking about libraries here, but issues with notation and array functionality that you run into in day-to-day programming. Things are pretty dang good right now, but there are some areas (array indexing, matrix multiply, etc.) that some people see as sub-optimal or better implemented in other langauges. A list will provide a "big picture" of where we want to go (maybe we are there...) and also help us pick our battles on language changes, etc. So, what mathematical expressions are commonly used and yet difficult to write in Python? I don't mean integrals divergence, etc. I mean things like matrix multiply and transpose. Here is a beginning to the list: 1. Matrix Multiply -- should we ask for ~*? 2. Transpose -- In a perfect world, we'd have an operator for this. 3. complex conjugate -- An operator for this would also be welcomed. 4. Others?? These three have all been discussed on this list or on the SciPy list in the last month, so they are obvious. I don't think there is a solution for 2 and 3 besides using the current function or method calls (but they are still on my list). As I mentioned in my last post, 1 might be fixable. As far as core Numeric functionality: 1. Array indexing with arrays. 2. .M attributes -- an alternative to (1) in language changes. And I'll add a third, that I'd like. 3. tensor notation indexing as in the Blitz++ array library http://www.oonumerics.org/blitz/manual/blitz03.html#l75 NewAxis and Ellipses allow for the same functionality, but the tensor notation is much easier to read. This requires yet more indexing trickery though... Please limit this thread to language changes or Numeric enhancements. Desired changes to the current behavior or interface of Numeric should be saved for a different discussion. thanks, eric -- Eric Jones <eric at enthought.com> Enthought, Inc. [www.enthought.com and www.scipy.org] (512) 536-1057 |
From: Andrew P. L. <bs...@al...> - 2002-03-07 10:30:51
|
On Thu, 7 Mar 2002, eric wrote: > 1. Matrix Multiply -- should we ask for ~*? > 2. Transpose -- In a perfect world, we'd have an operator for this. > 3. complex conjugate -- An operator for this would also be welcome I, personally, don't find the arguments particularly compelling for extra operators for numeric stuff. While extra operators may make code more "math-like", "MATLAB-like" or "Fortran-like", it won't help with efficiency. If I have code to compute A*x+B*y+C, I'm going to have to call out the A*x+Z and Z=B*y+C primitives as functions anyway. No set of binary operators will work out that optimization. Requesting domain specific operators actually scares me. The main problem is that it is *impossible* to remove them if your choices later turn out to be confusing or wrong. If operators must be added, I would rather see a generic operator mechanism in place in Python. Choice 1 would be a fixed set of operators getting allocated (~* ~+ ~- etc.) which the core language *does not use*. Then any domain can override with their special meaning without collapsing the base language under the weight of domain specific extensions. Now the specific domains can make their changes and only break their own users rather than the Python community at large. Choice 2 would be for a way for Python to actually adjust the interpretation semantics and introduce new operators from inside code. This is significantly trickier and more troublesome, but has the potential of being a much more generally useful solution (far beyond the realm of numerics). Furthermore, it allows people to make code look like whatever they choose. -a |
From: Perry G. <pe...@st...> - 2002-03-08 01:28:09
|
Eric makes a good point about stepping back and thinking about these issues in a broader context. Along these lines I'd like to make a proposal and see what people think. I think Konrad made a very good point about matrix vs array representation. If we made it illegal to combine them in expressions without explicit conversions, we could prevent much confusion about what kind of operations would be performed. An attempt to use one kind of in place of an other would trigger an exception and thus users would always know when that was a problem. Implementing this behavior in numarray would be simple, as would having both share the same implementation for common operations (without any extra performance penalty). That still leaves the question of how to do the conversions, i.e., one of the following options matrix(a) * b # matrix multiply of array (a) with matrix (b) a.M * b a.M() * b likewise: a * array(b) # element-wise multiply of array (a) with matrix (b) a * b.A a * b.A() I strongly prefer the first (functional) form. Rick White has also convinced me that this alone isn't sufficient. There are numerous occasions where people would like to use matrix multiply, even in a predominately "array" context, enough so that this would justify a special operator for matrix multiplication. If the Numeric community is united on this, I think Guido would be receptive. We might suggest a particular operator symbol or pair (triple) but leave him some room to choose alternatives he feels are better for Python (he could well come up with a better one). It would be nice if it were a single character (such as @) but I'd be happy with many of the other operator suggestions (~*, (*), etc.) Note that this does not imply we don't need a seperate matrix object. I think it is clear that simply providing a matrix multiply operator is not going to answer all their needs. As to the other related issues that Eric raises, in particular: operators for transpose and complex conjugate, I guess I don't see these as so important. Both of these are unary operators, and as such either of the following options does not seem to be notationally much worse (whereas using binary functions in place of binary operators is much less readable) transpose(x) conjugate(x) x.transpose() x.conjugate() x.T() x.C() x.T x.C (Personally, I prefer the first two) Perry |
From: Konrad H. <hi...@cn...> - 2002-03-08 08:05:25
|
"Perry Greenfield" <pe...@st...> writes: > That still leaves the question of how to do the conversions, i.e., one > of the following options ... > I strongly prefer the first (functional) form. Me too. I wouldn't call it "functional" though, it's exactly the way object constructors are written. > Rick White has also convinced me that this alone isn't sufficient. > There are numerous occasions where people would like to use matrix > multiply, even in a predominately "array" context, enough so that this > would justify a special operator for matrix multiplication. If the Could you summarize those reasons please? I know that there are applications of matrix multiplication in array processing, but in my experience they are rare enough that writing dot(a, b) is not a major distraction. Maybe we need to take another step back as well: Python is a general-purpose language, with several specialized subcommunities such as ours, some of them even more important in size. Most likely they are having similar discussions. Perhaps the database guys are discussing why they need two more special operators for searching and concatenating databases. I don't think such requests are reasonable. It is tempting to think that it doesn't matter, if you don't need that operator, you just don't use it. But a big advantage of Python is readability. If we get our (well, *yours*, I don't want it ;-) matrix multiply operator, a month later someone will decide that it's just great for his database application, and the database community will have to get used to it as well. > Numeric community is united on this, I think Guido would be receptive. > We might suggest a particular operator symbol or pair (triple) but Actually I feel quite safe: there might be a majority for another operator, but I don't expect we'd ever agree on a symbol :-) Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hi...@cn... Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais ------------------------------------------------------------------------------- |
From: Pearu P. <pe...@ce...> - 2002-03-08 10:26:24
|
Hi, On Fri, 8 Mar 2002, Konrad Hinsen wrote: <removed relevant comments> > > Numeric community is united on this, I think Guido would be receptive. > > We might suggest a particular operator symbol or pair (triple) but > > Actually I feel quite safe: there might be a majority for another > operator, but I don't expect we'd ever agree on a symbol :-) Thanks Konrad for this excellent point. Seeing all these proposals about solving our matrix-array issues makes also me feel safer with the current situation. My general point is that a _good_ solution is simple, if a solution is not simple, then it is probably a bad solution. I find separating array and matrix instances (in a sense of raising exception when doing <matrix> <op> <array>) not a very simple solution: New concepts are introduced that actually do not solve the simplicity problem of representing matrix operations. As I see it, they only introduce restrictions and the main assumption behind the rationale is that "users are dumb and they don't know what is best for them". This is how I interpret the raised exception as behind the scenes matrix and array are the same (in the sense of data representation). Let me remind at least to myself that the discussion started from a proposal that aimed to write matrix multiplication of arrays Numeric.dot(a,b) in somewhat simpler form. Current solution is to use Matrix class that being just a wrapper of arrays, redefines __mul__ method; and after one has defined a = Matrix.Matrix(a) the matrix multiplication of arrays looks simple: a * b Travis's proposal was to reduce the first step and have it inside an expression in a short form: a.M * b (there have been two implementation approaches proposed for this: (i) a.M returns a Matrix instance, (ii) a.M returns the same array with a temporarily set bit saying that the following operation is somehow special). To me, this looks like a safe solution. Though it is a hack, at least it is simple and understandable in anywhere where it is used (having a * b where b can be either matrix or array, it is not predictable from just looking the code what the result will be -- not very pythonic indeed). The main objection to this proposal seems to be that it deviates from a good pythonic style (ie don't mess with attributes in this way). I'd say that if python does not provide a good solution to our problem, then we are entitled to deviate from a general style. After all, in doing numerics the efficiency issue has a rather high weight. And a generally good style of Python cannot always support that. I guess I am missing here a goal to get Numeric or numarray joined to Python. With this perspective the only efficient solution seems to be introducing a new operator (or new operators). Few candidates have been proposed: a ~* b - BTW, to me this looks like dot(conjugate(a),b). a (*) b - note the in situ version of it: a (*)= b a (**) b - looks ugly enough? ;-) Actually, why not a [*] b, a {*} b for direct products of matrices (BTW (*) seems more appropriate here). So, my ideal preference would be: a .* b - element-wise multiplication of arrays, 2nd pref.: a * b a * b - matrix multiplication of arrays, 2nd preference: a [*] b a (*) b - direct matrix multiplication (also know as tensor product) of arrays a~ - conjugate of arrays a` - transpose of arrays This looks great but requires many new features to Python (new operators, the concept of right-hand unary operator). I don't think that Python should introduce these new operators just because of Numeric community. It is fine if they get used in other fields as well that suffer from the lack of operators. About unary operations: transpose and conjugate. BTW, in complex linear algebra their composition is equally frequent operator. Let me propose the following solution: to have a ** T for Numeric.transpose(a) a ** H for Numeric.transpose(Numeric.conjugate(a)) define T = TransposeOp() H = TransposeOp(conjugate=1) where class TransposeOp: def __init__(self, conjugate=0): self.conjugate = conjugate def __rpow__(self,arr): if self.conjugate: return Numeric.transpose(Numeric.conjugate(a)) return Numeric.transpose(arr) Looks Pythonic to me;-) Regards, Pearu |
From: Andrew P. L. <bs...@al...> - 2002-03-08 23:35:19
|
On Fri, 8 Mar 2002, Pearu Peterson wrote: > This is how I interpret the raised exception as behind the scenes matrix > and array are the same (in the sense of data representation). Matricies can have many different storage formats. Sparse, Banded, Dense, Triangular, etc. Arrays and Matricies are the same behind the scenes *for the moment*. For examples of matrix storage formats, check the BLAST Technical Forum documentation at netlib. Unfortunately, www.netlib.org appears to be down right now. Locking the assumption that Matrix and Array are "the same behind the scenes" into the main Python specification is not a good idea. -a |
From: Rick W. <rl...@st...> - 2002-03-08 13:21:53
|
On Fri, 8 Mar 2002, Konrad Hinsen wrote: > "Perry Greenfield" <pe...@st...> writes: > > > Rick White has also convinced me that this alone isn't sufficient. > > There are numerous occasions where people would like to use matrix > > multiply, even in a predominately "array" context, enough so that this > > would justify a special operator for matrix multiplication. If the > > Could you summarize those reasons please? I know that there are > applications of matrix multiplication in array processing, but in my > experience they are rare enough that writing dot(a, b) is not a major > distraction. A couple of quick examples: I do lots of image processing (e.g. deconvolution) using arrays. It is often helpful to take the outer product of two 1-D vectors; e.g. if there is a separable function f(x,y) = g(x)*h(y), you can compute separate g & h vectors and then combine them with outer product (a special case of matrix multiply) to get the desired 2-D image. Another example: when I'm working with either 2-D images or 1-D vectors, it is helpful to be able to compute projections using a set of basis vectors (e.g. for singular value decomposition, eigenvectors, etc.) This is most easily expressed using matrix multiplies - but most uses of the data still treat them as simple arrays instead of matrices. Being able to group these operations together is helpful both for readability of the code and for efficiency of execution. Having said that, I think I actually agree with Konrad that these sorts of operations are rare enough (in the data processing context) that it is no great burden to write them using function calls instead of operators. If we could agree on a matrix-multiply operator, that would be nice -- but if we can't, I can live with that too. For my purposes, I certainly don't see the need to add special operations to do things like transpose. Those should be limited to a separate matrix class as Konrad proposes and should be available as function calls for arrays. Rick ------------------------------------------------------------------ Richard L. White rl...@st... http://sundog.stsci.edu/rick/ Space Telescope Science Institute Baltimore, MD |
From: Paul F D. <pa...@pf...> - 2002-03-08 16:24:41
|
To sum up my previous postings: I prefer the "constructor" notation, not the funny-attribute notation. However, I believe an efficient matrix class can be done in Python on top of a numeric array object. Shadow classing just doesn't seem to be an overhead problem in most cases. The coding for MA is very complicated because of its semantics are so different; for Matrix it would be much less complicated. When I designed Basis (1984) I was faced with the operator issue. The vast majority of the operations were going to be elementwise but some users would also want matrix multiply and the solution of linear systems. (Aside: a/b meaning the solution x of bx = a, is best calculated not as (b**-1) * a but by solving the system without forming the inverse.) My choice was: matrix multiply and divide were *!, /!. This was successful in two senses: the users found it easy to remember, and you could implement it in the tokenizer or just in the grammar. I chose to make it a token so as to forbid internal spaces, but for Python it could be done without touching the tokenizer, and it doesn't use up any new symbols. Symbols are precious; when I designed Basis I had a keyboard map and I would cross out the keys I had "used up". If I were Guido I would be very reluctant to give up anything valuable like @ for our purposes. One should not have any illusions: putting such operators on the array class is just expediency, a way to give the arrays a bit of a dual life. But a real matrix facility would have an abstract base class, be restricted to <= 2 dimensions, have realizations including symmetric, Hermitian, sparse, tridiagonal, yada yada yada. Another aside: There are mathematical operations whose output type depends not just on the input type but some of these other considerations. This led me to believe that the Numpy approach is essentially correct, that the type of the elements be variable rather than having separate classes for each type. |
From: Konrad H. <hi...@cn...> - 2002-03-08 17:55:08
|
> The Python core has long had at least 2 examples of operators which > act as object constructors: 'j' which performs complex() and 'L' which > performs long() (you can't get much more `pythonic' than a built-in > type). Those are suffixes for constants, not operators. If they were operators, you could apply them to variables - which you can't. More importantly, the L suffix wouldn't even work as an operator, as the preceding number might extend the range of integers before it has a chance of being converted to a long integer. > I would venture to say that the numeric community is pretty high up > there in importance if not size, given the early appearance of the > complex number type and strong math capacity not to mention GvR's The complex type was introduced for the benefit of NumPy (I remember it all too well, as I did the initial implementation), but after a long discussion on the Python list, with many expressing disapprovement because of its special-need status. I'd say it shows the limits of what one can get accepted. Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hi...@cn... Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais ------------------------------------------------------------------------------- |