![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
One thing keeps bugging me when I use numpy.matrix. All this is fine:: >>> x=N.mat('1 1;1 0') >>> x matrix([[1, 1], [1, 0]]) >>> x[1,:] matrix([[1, 0]]) But it seems to me that I should be able to extract a matrix row as an array. So this :: >>> x[1] matrix([[1, 0]]) feels wrong. (Similarly when iterating across rows.) Of course I realize that I can just :: >>> x.A[1] array([1, 0]) but since the above keeps feeling wrong I felt I should raise this as a possible design issue, better discussed early than latter. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/37c1dd52b545c64455f1de5977d16dff.jpg?s=120&d=mm&r=g)
Em Dom, 2007-03-25 às 13:07 -0400, Alan G Isaac escreveu:
I think the point here is that if you are using matrices, then all you "should" want are matrices, just like in MATLAB: >> A = [1 2; 3 4] A = 1 2 3 4 >> b = A(1, :) b = 1 2 >> size(b) ans = 1 2 >> b = A(:, 1) b = 1 3 >> size(b) ans = 2 1 >> b = 1 b = 1 >> size(b) ans = 1 1 You see, rows, columnes, and even numbers, are treated as matrices. Paulo -- Paulo José da Silva e Silva Professor Assistente do Dep. de Ciência da Computação (Assistant Professor of the Computer Science Dept.) Universidade de São Paulo - Brazil e-mail: rsilva@ime.usp.br Web: http://www.ime.usp.br/~rsilva Teoria é o que não entendemos o (Theory is something we don't) suficiente para chamar de prática. (understand well enough to call practice)
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Sun, 25 Mar 2007, Paulo Jose da Silva e Silva apparently wrote:
Yes, that is the idea behind this, which I am also accustomed to from GAUSS. But note again that the Matlab equivalent :: >>> x=N.mat('1 2;3 4') >>> x[0,:] matrix([[1, 2]]) does provide this behavior. The question I am raising is a design question and is I think really not addressed by the rule of thumb you offer. Specifically, that rule of thumb if it is indeed the justification of :: >>> x[1] matrix([[3, 4]]) finds itself in basic conflict with the idea that I ought to be able to iterate over the objects in an iterable container. I mean really, does this not "feel" wrong? :: >>> for item in x: print item.__repr__() ... matrix([[1, 2]]) matrix([[3, 4]]) Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
This may sound silly, but I really think seeing all those brackets is what makes it feel wrong. Matlab's output doesn't put it in your face that your 4 is really a matrix([[4]]), even though that's what it is to Matlab. But I don't see a good way to change that behavior. The other thing I find problematic about matrices is the inability to go higher than 2d. To me that means that it's impossible to go "pure matrix" in my code because I'll have to switch back to arrays any time I want more than 2d (or use a mixed solution like a list of matrices). Matlab allows allows >2D. --bb
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Colin J. Williams <cjw@sympatico.ca> wrote:
I'm not sure what you thought I meant, but all I meant by going "pure matrix" was having my Numpy code use the 'matrix' type exclusively instead of some mix of 'matrix' and the base 'ndarray' type. Things become messy when you mix and match them because you don't know any more if an expression like A[1] is going to give you a 1-D thing or a 2-D thing, and you can't be sure what A * B will do without always coercing A and B.
A list of matrices seems to be a logical structure.
Yes, and it's the only option if you want to make a list of matrices of different shapes, but I frequently have a need for things like a list of per-point transformation matrices. Each column from each of those matrices can be thought of as a vector. Sometimes its convenient to consider all the X basis vectors together, for instance, which is a simple and efficient M[:,:,0] slice if I have all the data in a 3-D array, but it's a slow list comprehension plus matrix constructor if I have the matrices in a list -- something like matrix([m[:,0] for m in M]) but that line is probably incorrect.
Numpy generally tries to treat all lists and tuples as array literals. That's not likely to change. --bb
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Bill Baxter wrote:
Yes, to my mind it's best to consider the multi-dimensional array and the matrix to be two distinct data types. In most cases, it's best that conversions between the two should be explicit.
Logically, this makes sense, where M is a list of matrices. My guess is that it would be a little faster to build one larger matrix and then slice it as needed.
That need no be a problem is there is clarity of thinking about the essential difference between the matrix data type (even if is is built as a sub-type of the array) and the multi-dimensional array.
--bb
Colin W.
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
I mean really, does this not "feel" wrong? ::
On Mon, 26 Mar 2007, Bill Baxter apparently wrote:
This may sound silly, but I really think seeing all those brackets is what makes it feel wrong.
I appreciate the agreement that it feels wrong, but I dispute the analysis of this symptom. What makes it "feel wrong", I contend, is that experience with Python make this a **surprising** behavior. And that is precisely why I suggest that this may point to a design issue. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
So you're saying this is what you'd find more pythonic?
Probably about half the bugs I get from mixing and matching matrix and array are things like row = A[i] ... z = row[2] Which works for an array but not for a matrix. I think Matlab makes it more bearable by having a single value index like X[i] be equivalent to X.flat()[i]. So X[2] is the same for row or col vec in Matlab. Now that I think about it, that's probably the main reason I feel more comfortable with array than matrix in Numpy. If I have a vector, I should only need one index to get the ith component. --bb
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Bill Baxter apparently wrote:
No; that is not possible, since a matrix is inherently 2d. I just want to get the constituent arrays when I iterate over the matrix object or use regular Python indexing, but a matrix when I use matrix/array indexing. That is :: >>> X[1] array([2,3]) >>> X[1,:] matrix([[3, 4]]) That behavior seems completely natural and unsurprising.
Exactly! That is the evidence of a "bad surprise" in the current behavior. Iterating over a Python iterable should provide access to the contained objects. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Oooops, they should match of course. :: >>> X[1] array([3,4]) >>> X[1,:] matrix([[3, 4]]) But again the point is: indexing for submatrices should produce matrices. Normal Python indexing should access the constituent arrays. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Alan G Isaac schrieb:
I think this is a tricky business. There's also the rule "indexing reduces the rank, slicing preserves it". Numpy-matrices form an exception to this rule, always being 2d, but the exception is consistently enforced. Now it seems you want to have an exception from the exception, correct? Isn't this why the .A1 method for numpy-matrices was introduced even after 1.0rc? -sven (matrix fan)
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sven Schreiber apparently wrote:
What I want: the best design. I did not claim that the current design is flawed, only to suspect it. Why I wrote: current behavior feels wrong -> suspect design flaw. What feels wrong: iterating over a container does not give access to the contained objects. This is not Pythonic. *Symptom* of the underlying problem: for matrix M, M[0] returns a matrix. Would the change I suggest mean that the behavior of the matrix class deviates less from the array class: yes. Does this mean the matrix class behavior would be less "consistent"? I have tried to explain why the answer is "no". hth, Alan Isaac
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Alan G Isaac schrieb:
What feels wrong: iterating over a container does not give access to the contained objects. This is not Pythonic.
If you iterate over the rows of the matrix, it feels natural to me to get the row vectors -- and as you know a 1d-array does not contain the information of "row-ness", so you need to get a 2d thing to properly return those "contained objects". So in my view the current behavior actually does exactly what you're calling for. (But I admit I'm a 2d fan so much so that I didn't even know that using a single index is possible with a matrix. I thought one always had to be explicit about both dimensions... So thanks for pointing this out. -- BTW, I'm aware that sticking to numpy-matrices exclusively is not possible. For example, eigenvalues are returned in a 1d-array even if you pass a matrix, and it's intended behavior.) cheers, sven
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sven Schreiber apparently wrote:
If you iterate over the rows of the matrix, it feels natural to me to get the row vectors
Natural in what way? Again, I am raising the question of what would be expected from someone familiar with Python. Abstractly, what do you expect to get when you iterate over a container? Seriously.
But I admit I'm a 2d fan so much so that I didn't even know that using a single index is possible with a matrix.
Exactly. When you want submatrices, you expect to index for them. EXACTLY!! Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/6194b135cba546afa82516de1537de49.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
If may chime in... I think Sven's argument in on the side saying, A "matrix" is an object that you expect a certain (mathematical !) behavior from. If some object behaves intuitively right -- that's ultimately pythonic ! The clash is, NOT to see a matrix "just as another container". Instead a matrix is a mathematical object , that has rows and columns. It is used in a field (lin-alg) where every vector is either a row or a column vector -- apparently that's big thing ;-) The whole reason to add a special matrix class to numpy in the first place, is to provide a better degree of convenience to lin-alg related applications. I would argue that it was just not consistently considered, that this should also come with "a column of a matrix is something else than a row -- (1,n) vs. (n,1) and not (n,) more notes/points: a) I have never heard about the m.A1 - what is it ? b) I don't think that if m[1] would return a (rank 2) matrix, that m[1].A could return a (rank 1) array ... c) I'm curious if there is a unique way to extend the matrix class into 3D or ND. -Sebastian
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sebastian Haase apparently wrote:
The problem is, as I am not the only one to point out, this particular behavior is NOT intuitively right.
The clash is, NOT to see a matrix "just as another container".
But be serious, no object is "just another container". Again, this just begs the question. The question is a design question. E.g., what is the principle of least surprise?
more notes/points: a) I have never heard about the m.A1 - what is it ?
It returns a 1d array holding the raveled matrix.
b) I don't think that if m[1] would return a (rank 2) matrix, that m[1].A could return a (rank 1) array ...
It does not, of course. (But both should, I believe.)
c) I'm curious if there is a unique way to extend the matrix class into 3D or ND.
Is that not what an array is for?? Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps things would be clearer if we thought of the constituent groups of data in a matrix as being themselves matrices. X[1] could represent the second row of a matrix. A row of a matrix is a row vector, a special case of a matrix. To get an array, I suggest that an explicit conversion X[1].A is a clearer way to handle things. Similarly, X[2, 3] is best returned as a value which is of a Python type. Colin W.
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps it would now help if you redefined the question. In an earlier posting, you appeared anxious that the matrix and the array behave in the same way. Since they are different animals, I see sameness of behaviour as being lower on the list of desirables than fitting the standard ideas of matrix algebra. Suppose that a is a row vector, b a column vector and A a conforming matrix then: a * A A * b and b.T * A are all acceptable operations. One would expect the iteration over A to return row vectors, represented by (1, n) matrices. Colin W.
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, "Colin J. Williams" apparently wrote:
One would expect the iteration over A to return row vectors, represented by (1, n) matrices.
This is again simple assertion. **Why** would "one" expect this? Some people clearly do not. One person commented that this unexpected behavior was a source of error in their code. Another person commented that they did not even guess that such a thing would be possible. Experience with Python should lead to the ability to anticipate the outcome. Apparently this is not the case. That suggests a design problem. What about **Python** would lead us to expect this behavior?? In *contrast*, everyone agrees that for a matrix M, we should get a matrix from M[0,:]. This is expected and desirable. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
Well, and what *should* they expect. I think it is expected because the for iterates over rows and rows of matrices are 1xn. Matrices and arrays, as has been stated, are different animals. Probably it would have been best to stick with arrays and I suspect that matrices appeared because of the dearth of Python operators, in particular to make matrix multiplication simpler. On the other hand, certain errors slip by when one is implementing matrix algebra with arrays, but they can be avoided by never using 1-d vectors. So all this mess results from catering to the matrix community. Matlab has the opposite problem, multidimensional arrays were tacked on later and they don't operate properly with everything. Chuck One person commented that this unexpected behavior was
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Charles R Harris apparently wrote:
Well, and what should they expect.
Since matrices are an iterable Python object, we *expect* to iterate over the contained objects. (Arrays.) I am not sure why this is not evident to all, but it is surely the sticking point in this discussion. A matrix is not a container of matrices. That it acts like one is surprsing. Surprises are bad unless they have a clear justification. Perhaps a clear justification exists, but it has not been offered in this discussion. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
I think that a clear justification has been offered several times on the list recently, though maybe not all in this thread. Matrices in numpy seem to exist as a construct to facilitate linear algebra. One such property is that row- and column-vector slices of matrices are (1xN) or (Nx1) matrices, because otherwise common linear algebra things -- like how you multiply a row-vector or a column vector by a matrix, and whether and when it needs to be transposed -- do not translate properly from "linear algebra notation" like in textbooks and papers. Once the matrix class is committed to providing row- and column- vector slices as other matrices, it makes no sense to have iteration over the matrix provide a different data type than slicing. Basically, my rule of thumb has been to *only* use matrices when I'm copying linear algebra operations out of textbooks and papers, and I want the notations to align. Doing other, non-linear-algebra operations with matrices -- like iterating over their elements -- is not really worth it. There's a meta question, of course: should things be changed to make it "worth it" to do "pythonic" tasks with matrices? Should there be special elementwise vs. matrix-operation operators? Should numpy work a lot more like matlab? That discussion is on-going in another thread, I think. Zach
![](https://secure.gravatar.com/avatar/4d021a1d1319f36ad861ebef0eb5ba44.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
It actually has been offered. You just don't accept it. Matrices are containers of matrices. If M is an (mxn) matrix then M[0] is a (1xn) matrix. Viewing this 1xn matrix as a 1-d array loses it's row-vectorness. This seems perfectly logical and acceptable to me. I'm waiting for a better explanation as to why this is not acceptable. Arguments that rest on what is and what isn't "Pythonic" are not very persuasive as this is very often in the eye of the beholder. -Travis
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Travis Oliphant wrote:
Again, I only raised a *question* about whether there might be a design problem here. My goal was only to have this explored, and I've tried to explain why. The evidence I offer: - it is surprising (at least to some) that iterating over a matrix yields matrices - I believe it is unexpected (prior to instruction) and that there is a natural more expected behavior - if that is right, deviation from the expected should have a good justification - this behavior has tripped up at least a couple people and I expect that to happen to many others over time (because I believe the behavior is unexpected) - when I desire to iterate over a matrix I always want the arrays. (Of course there is a way to have them; that's not the point). I'm interested to see a use case where the rows are desired as matrices As you indicate, none of this constitutes an "argument". And since nobody appears to agree, I should shut up. This will be my last post on this subject. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/4d021a1d1319f36ad861ebef0eb5ba44.jpg?s=120&d=mm&r=g)
Alan Isaac wrote:
Thanks for listing your points. I can see that this is an issue where reasonable people will disagree because there are multiple ways of looking at it. The idea that matrix selection would return matrices really comes from wanting to keep matrices more persistent through operations. M[0] could be seen either as a 1xn matrix or a n-length array. I agree that both concepts are possible. Seeing it as a 1xn matrix allows matrices to remain persistent more often. So, the arguments for the current approach and the arguments against it to me seem on the same level, so I don't see a reason to change the current behavior and see a lot of strong reasons not to change the behavior (we are in a 1.0 release and could not change anything until at least 1.1 anyway). With that said: One of my goals for the next year or two is to create a matrix class in C and incorporate CVXOPT matrices and it's coupled sparse matrix. We can re-visit this question in that process. I would like for there to be a sparse matrix implementation in NumPy (without the solver which will remain in SciPy). But, the sparse matrix and the matrix need to have the same behaviors and should be able to interoperate with each other. So, if you would like to help with that project all input is welcome. Best regards, -Travis
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
It seems to me that using shapes, (m,1) versus (1,n), to determine whether a vector is column- or row-oriented is a hack (or at least feels like one). If we know we have a vector, then we want to use a single index to obtain a scalar, and that extra "0," or ",0" shouldn't be necessary. It also seems to me that there is an object-oriented design solution to this, namely defining column_vector and row_vector classes that inherit from matrix, that accept a single index, but behave as matrices consistent with their type. I'm sure there are implications to this I haven't thought through yet, but at least it gives us logical indexing AND persistence of matrices through operations. On Mar 27, 2007, at 1:11 PM, Travis Oliphant wrote:
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
of the tensor product of a vector space with its dual, so one would expect either row or column vectors depending on the index position. But this discussion is not really relevant at this point. What happens is a convention and the numpy convention is out there. It isn't going to change now, it would break too much code. Chuck
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Charles R Harris wrote:
What happens is a convention
Certainly true.
It isn't going to change now, it would break too much code.
That is a different kind of argument. It might be true. May I see a use case where the desired return when iterating through a matrix is rows as matrices? That has never been what I wanted. Thank you, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/27/07, Alan Isaac <aisaac@american.edu> wrote:
Yeh, I was afraid you'd say that. :-) But maybe I've got a lot of points, and I don't feel like making a copy of the whole set. Or maybe it's not a linear transform, but instead xformedPt = someComplicatedNonLinearThing(pt) I do stuff like the above quite frequently in my code, although with arrays rather than matrices. :-) For instance in finite elements there's assembling the global stiffness matrix step where for each node (point) in your mesh you set some entries in the big matrix K. Something like for i,pt in enumerate(points): K[shape_fn_indices(i)] = stiffness_fn(pt) That's cartoon code, but I think you get the idea. I don't think there's any good way to make that into one vectorized expression. The indices of K that get set depend on the connectivity of your mesh.
Note that I am no longer disputing the convention, just trying to understand its usefulness.
--bb
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Bill Baxter apparently wrote:
xformedPt = someComplicatedNonLinearThing(pt)
I do stuff like the above quite frequently in my code, although with arrays rather than matrices.
Exactly: that was one other thing I found artificial. Surely the points will then be wanted as arrays. So my view is that we still do not have a use case for wanting matrices yielded when iterating across rows of a matrix. Alan Isaac PS In the "large number of points" case, the thing to do would be to extract modest sized submatrices.
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
It's pretty clear from my perspective: 1-D slices of matrices *must* be matrices. I gave an intuitive make-it-fit-with-known-notation explanation, and Charles gave a theoretically-grounded explanation. There was no objection to this from any quarter. Given that 1-D slices of matrices must be matrices for deep mathematical reasons as well as notational compatibility with the rest of linear algebra -- e.g. that m[0] yields a matrix if m is a matrix-- it almost certainly would violate the principle of least surprise for iteration over m (intuitively understood to be choosing m [0] then m[1] and so forth) to yield anything other than a matrix. This can't possibly be what you're asking for, right? You aren't suggesting that m[0] and list(iter(m))[0] should be different types? There are many clear and definite use cases for m[0] being a matrix, by the way, in addition to the theoretical arguments -- these aren't hard to come by. Whether or nor there are actual good use-cases for iterating over a matrix, if you believe that m[0] should be a matrix, it's madness to suggest that list(iter(m))[0] should be otherwise. My opinion? Don't iterate over matrices. Use matrices for linear algebra. There's no "iteration" construct in linear algebra. The behavior you find so puzzling is a necessary by-product of the fundamental behavior of the matrix class -- which has been explained and which you offered no resistance to. Zach
![](https://secure.gravatar.com/avatar/5dde29b54a3f1b76b2541d0a4a9b232c.jpg?s=120&d=mm&r=g)
Zachary Pincus wrote:
I don't think the OP was suggesting that. Rather, I think the suggestion was that for a mXn matrix, M: M[i].shape == (n,) M[i:i+1].shape == (1,n) that is, indexing (or iterating returns a vector, and slicing returns a matrix). This is, indeed exactly how numpy arrays behave! The problem with this is: numpy matrices were created specifically to support linear algebra calculations. For linear algebra, the distinction between row vectors and column vectors is critical. By definition, a row vector has shape: (1,n), and a column vector has shape (m,1). In this case, perhaps the OP is thinking that a shape (n,) array could be considered a row vector, but in that case: M[1,:].shape == (n,) M[:,1].shape == (m,) which of these is the row and which the column? This is why matrices index this way instead: M[1,:].shape == (1, n) M[:,1].shape == (m, 1) now we know exactly what is a row and what is a column. By the way, I think with the way numpy works: M[i] == M[i,:] by definition, so those couldn't yield different shaped results. Is that right? I think we got a bit sidetracked by the example given. If I have a bunch of points I want to store, I'm going to use an (m,2) *array*, not a matrix, then then A[i] will yield a (2,) array, which makes sense for (2-d) points. In fact, I do this a LOT. If I happen to need to do some linear algebra on that array of points, I'll convert to a matrix, do the linear algebra, then convert back to an a array (or just use the linear algebra functions on the array). I hope this helps -Chris
-- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
![](https://secure.gravatar.com/avatar/37ee74afc8985c99adac9eb37ab2b10e.jpg?s=120&d=mm&r=g)
I would like to hear your opinion on developing an explicit sparse/dense 2D matrix class with indexing similar to Matlab, and without significant differences between sparse and dense matrices from the user's perspective... I know that it's not one of Numpy/Scipy's goals to clone Matlab, but I think I represent a potentially large scientific Python user base, who would find Python matrices that feel a bit more like Matlab/Octave etc. extremely useful. I have a great respect for all the work in Numpy and Scipy, but at the same time I feel that numerical linear algebra in Python would benefit from a more dedicated matrix library, and I think that Numpy (or another single package) should provide that in a homogeneous way - without ties to how Numpy array works. - Joachim On 3/27/07, Christopher Barker <Chris.Barker@noaa.gov> wrote:
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
Hello all, I suspect my previous email did not contain the full chain of my reasoning, because I thought that some parts were basically obvious. My point was merely that given some pretty basic fundamental tenants of numpy, Alan's suggestions quickly lead to *serious issues* far worse than the original problem. Thanks Chris for making this more explicit; here I will do so further. Let A be an array of arbitrary shape, and M be m-by-n matrix. (1) A[i] equals by definition A[i,...] This is a pretty fundamental part of numpy, right? That property (and its relatives) along with the broadcasting behavior, are in large part what make numpy so very expressive compared to matlab, etc. I can write code in numpy to do multidimensional operations over arbitrary-shaped and arbitrary-dimensioned data that look *exactly* like the corresponding 1D operations. It's amazing, really, how numpy's "implicit indexing" (this property), broadcasting (sort of the inverse of this, really, where shape tuples are extended as necessary with implicit 1's, instead of slices being extended with :'s as necessary), and fancy indexing, allow complex operations to be written compactly, clearly, and in a fully-general manner. I thought that point (1) was obvious, as it's a pretty fundamental part of numpy (and Numeric and Numarray before that). It's in absolutely no way a "trivial convenience". Of course, broadcasting and ufuncs, etc, contribute a lot more to expressiveness than this property, but it's part of the family. (2) Thus, M[i] equals by definition M[i,:]. By the mathematical definition of matrices and for pragmatic reasons, M[i,:] is another matrix. Here the trouble starts. Should matrices, not being arbitrary in their dimension, even have the ability to be indexed as M[i]? Or should they be completely different beasts than arrays? Which violates the "principle of least surprise" more? Now, I think we all more or less agree that M[i,:] has to be a matrix, because otherwise the matrix class is pretty worthless for linear algebra, and doesn't operate according to the mathematical definition of a matrix. Indeed, this is one of the basic properties of the numpy matrix -- that, and that the multiplication operator is defined as matrix-multiplication. Remove that and there's almost no point in *having* a separate matrix class. Want a use-case? Do the singular value decomposition, and get the product of the second-smallest singular vectors as Vt[:,-2] * U [-2,:]. (Things vaguely akin to this come up a lot.) You can copy the operation from any linear algebra text if the row- and column-vectors are treated as above. Otherwise you have to wonder whether that's supposed to be an outer or inner product, and use the corresponding operator -- and at that point, why was I using matrix classes anyway? (3) If M[i] = M[i,:], and M[i,:] is a matrix, then the iteration behavior of matrices *has* to yield other matrices, unless Alan is willing to suggest that list(iter(m))[i] should have a different type than M[i] -- a clear absurdity. This is the point that I was trying to make previously, having mistakenly assumed that point (1) was clear to all, and that point (2) had been clearly established. So, if we trace the reasoning of Alan's suggestion, coupled with these above properties, we get just that: (a) Iteration over M should yield arrays. (b) Ergo M[i] should be an array. (c) Ergo M[i,:] should be an array. (d) Ergo the matrix class should be identical to the array class except for overloaded multiplication, and the only way to get a 'row' or 'column' vector from a matrix that operates as a proper row or column vector should is M[i:i+1,:]. Whicy is a lot of typing to get the 'i-th' vector from the matrix. I don't think these changes are worth it -- it basically discards the utility of the matrix class for linear algebra in order to make the matrix class more useful for general purpose data storage (a role already filled by the array type). Here is another set of changes which would make matrices fully consistent with their linear-algebra roots: (a) M[i,:] should be a matrix. (b) M[i] is an error. (c) Iteration over M is an error. That's kind of lame though, right? Because then matrices are completely inconsistent with their numpy roots. So we are left with the current behavior: (a) M[i,:] is a matrix. (b) M[i] is a matrix. (c) Iteration over M yields matrices. My view is that, fundamentally, they are linear algebra constructs. However, I also think it would be far more harmful to remove basic behavior common to the rest of numpy (e.g. that M[i] = M[i,:]), and behavior that comes along with that (iterability). Hence my suggestion: treat matrices like linear algebra constructs -- don't iterate over them (unless you know what you're doing). Don't index them like M[i] (unless you know what you're doing). Maybe it's "just a little bizarre" to suggest this, but I think the other options above are much more bizarre. Zach
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Hi Zach, The use case I requested was for iteration over a matrix where it is desirable that matrices are yielded. That is not what you offered. The context for this request is my own experience: whenever I have needed to iterate over matrices, I have always wanted the arrays. So I am simply interested in seeing an example of the opposite desire. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Gram-Schmidt orthogonalization. ortho = [mat[0] / sqrt(mat[0] * mat[0].T)] for rowv in mat[1:]: for base in ortho: rowv = rowv - base * float(rowv * base.T) rowv = rowv / sqrt(rowv * rowv.T) ortho.append(rowv) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Robert Kern apparently wrote:
Gram-Schmidt orthogonalization
I take it from context that you consider it desirable to end up with a list of matrices? I guess I would find it more natural to work with the arrays, but perhaps that starts just being taste. Thank you, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Honestly, I don't care. You asked about iteration, and I gave you an example where it was important that the iterates were matrices (such that .T worked correctly). Please stop moving the goal posts.
I guess I would find it more natural to work with the arrays, but perhaps that starts just being taste.
Personally, I find working with arrays from start to finish the most natural. If you really want iterates to be arrays, just use .A. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
No, I offered a thorough discussion of why the design of numpy, as I imperfectly understand it, make the trade-offs to achieve your desired goal unpalatable.
I'm glad that Robert provided such a nice example. Nevertheless, I think that even absent any example it should be pretty clear at this point that if you want iteration over matrices to return arrays, either: (1) you need to be willing to accept that list(iter(M))[i] != M[i] or (2) M[i], which equals M[i,:], should be an array. (or 3, that M[i] and iteration over M is not defined) Clearly, (1) is fundamentally ridiculous. Clearly use-cases (as well as the mathematical definition of matrices) abound to illustrate why (2) is no good. Option (3) helps nobody. So which is it? You have to choose one! It is my assertion that if you are iterating over matrices and hoping to get arrays, you are misusing the matrix class. Almost by definition if you expect M[i,:] or M[i] (or equivalently, the output of an iterator over M) to be an array, then M is being used outside of a linear algebra context -- and thus, in my view, misused. In those cases, the matrix should be cast to an array, which has the desired behavior. Robert nicely illustrated a linear-algebraic context for matrix iteration. Now, Bill offers up a different suggestion: indexing M yields neither a matrix nor an array, but a class that operates more or less like an array, except insofar as it interacts with other matrix objects, or other objects of similar classes. I'm interested in hearing more about this, what trade-offs or other compromises it might involve. Zach
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Zachary Pincus wrote:
M[i], which equals M[i,:]
Of course this equality must break. That was stated at the outset. As I said before, this may be impossible or undesirable. But, as I said before, it seems prima facie natural for M[i] to be ordinary Python indexing while M[i,:] is numpy indexing, each with its own natural interpretation. In any case, you may be spending too much energy on this. I am not a developer, and no developer has expressed a desire to make such a change. The current behavior is currently safe (if odd). Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
What I envisioned was that M[i,:] would return a row_vector and M [:,j] would return a column_vector, because this would be symmetric behavior. M[i], by convention, would behave the same as M[i,:]. But then I personally don't distinguish between "python indexing" and "numpy indexing". In both cases, __getitem__() (or __setitem__()) is called. For multiple indexes, the index object is a tuple. In any case, the behavior of "numpy indexing" as I have proposed it would return an object that inherits from matrix, thus would BE a matrix, although it would behave like an array. On Mar 29, 2007, at 3:24 PM, Alan G Isaac wrote:
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Bill Spotz <wfspotz@sandia.gov> wrote:
I'm thinking that a basic problem is not having row and column types distinct from matrices. Column types would represent elements of a vector space and row types elements of the dual, they would both be 1-dimensional. One could go full bore with tensors and index signatures, upper and lower, but I think plain old matrices with the two vector types would solve most ambiguities. Thus if v were a column vector, then v.T*v would be a scalar where the product in this case is shorthand for v.T(v). Likewise, v*v.Twould be a matrix (in tensor form the indices would be upper lower respectively, but ignore that). The default translation of a normal 1-D vector would preferably be a column vector in this case, opposite to current usage, but that is not really crucial. Note that nx1 and 1xn matrices are not quite the same thing as column and row vectors, as the latter would normally be considered linear maps and v.T*v would return a scalar matrix, quite a different thing from a scalar. Chuck
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Thu, 29 Mar 2007, Bill Spotz apparently wrote:
Can you please be explicit about the value added by that convention, as you see it? Thank you, Alan Isaac PS I assume your "vector" class would always accept a single index, for both row and column vectors? What is the return type of v[i]?
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
On Mar 29, 2007, at 6:48 PM, Alan G Isaac wrote:
I assume by "convention" you are referring to the convention that M [i] is M[i,:] is a row_vector. I see it as a design decision: Does M[i] mean something? No: Then raise an exception. If users want to iterate over M, this forces them to use the double index notation to specify what they are iterating over. I would find this acceptable. Yes: What should M[i] represent/return? Array representing row vector: Array representing column vector: I would argue against linear algebra objects returning raw arrays because their meaning/interface can be ambiguous, as this thread demonstrates row_vector: Reasonable choice, analogous to conventional mathematical notation column_vector: Reasonable choice, but less analogous to conventional mathematical notation Of all the choices, raising an exception and returning a row_vector are (in my opinion, anyway) the best. Unless someone has a better argument.
My best guess at a best design would be that it would return a scalar. And yes, vectors would always accept a single index. But because they are matrices, they would accept two indexes as well.
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/5b2449484c19f8e037c5d9c71e429508.jpg?s=120&d=mm&r=g)
"""matrix.py The discussion about matrix indexing has been interminible and for the most part pretty pointless IMO. However, it does point out one thing: the interaction between the matrix and array classes is still pretty klunky despite a fair amount of effort trying to make them interoperate. Here's one possible alternative approach to enabling matrix operations within the context of numpy arrays. I started with a few requirements: 1. Matrices should be proper subclasses of arrays in the sense that one should be able to use matrices wherever one uses arrays with no change in the result. 2. Indexing into matrices should produce row or column vectors where appropriate. 3. There should be some syntax for matrix multiplication. I know that there are other operations defined on the matrix class, but I don't consider them worth the headache of having two class hierarchies. This is in part inspired by some comments Bill Spotz, but after listening to his later comments, I don't think that this is really what he had in mind. So, please don't blame him for any deficiencies you find herein. The usual caveats apply: this is just a sketch, it's probably buggy, use at your own risk, etc. However, you should be able to save this mail as "matrix.py" and execute it to run all of the examples below. I enforce point #1 by not overriding any ndarray methods except for __array_finalize__ and __getitem__ and they are overridden in a way that shouldn't effect code unless it is specifically testing for types. With that summarily taken care of, let's move onto point #2: indexing. First let's create a matrix. Scratch that, let's create an array of matrices. One things that this approach does is let you create arrays of matrices and vectors in a natural way. The matrix below has a shape of (3,2,2), which correspons to an array of 3 2x2 matrices.
m = matrix([ [[1, 0],[0,1]], [[2, 0],[0,2]], [[3, 0],[0,3]] ])
If we index on the first axis, we are selecting a single matrix, so we would expect an ndmatrix back.
On the other hand, if we index along the second axis, we are indexing on the matrices columns and thus expect to get an array of row matrices.
Finally if we index on the last index, we get an array of column matrices.
Indexing of row and column matrices works similarly. It's not shown here and thus will probably turn out to be buggy. For matrix multiplication, I chose to abuse call. Perhaps someday we'll be able to convinve the powers that be to free up another operator, but for the time being this is better than "*" since that forces matrix and array into separate camps and arguably more readable than most of the new syntax proposals that I've seen. Here we multiply our set of three matrices 'm', with a single 2x2 matrix 'n':
Things behave similarly when multiplying a matrix with a row vector or a row vector with a column vector.
Note, however that you can't (for instance) multiply column vector with a row vector:
Finally, you'd like to be able to transpose a matrix or vector. Transposing a vector swaps the type from ndrow and ndcolumn while transposing a matrix transposes the last two axes of the array it is embedded in. Note that I use ".t" for this instead of ".T" to avoid interfering with the existing ".T". (Frankly, if it were up to me, I'd just deprecate .T).
That's about it. Enjoy (or not) as you will. """ from numpy import * _ULT = 1 # Utlimate AKA last _PEN = 2 # Penultimate AKA second to last # I'm sure there's a better way to do this and this is probably buggy... def _reduce(key, ndim): if isinstance(key, tuple): m = len(key) if m == 0: return 0 elif m == 1: key = key[0] else: if Ellipsis in key: key = (None,)*(ndim-m) + key if len(key) == ndim: return isinstance(key[-1], int) * _ULT + isinstance(key[-2], int) * _PEN if len(key) == ndim-1: return isinstance(key[-1], int) * _PEN if isinstance(key, int): if ndim == 1: return _ULT if ndim == 2: return _PEN return 0 class ndmatrix(ndarray): def __array_finalize__(self, obj): if obj.ndim < 2: raise ValueError("matrices must have dimension of at least 2") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT | _PEN: return value.view(ndarray) if rmask == _ULT: return value.view(ndcolumn) if rmask == _PEN: return value.view(ndrow) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other[...,newaxis,:], -1).view(ndcolumn) if isinstance(other, ndmatrix): return sum(self[...,newaxis] * other[...,newaxis,:], -2).view(ndmatrix) else: raise TypeError("Can only matrix multiply matrices by matrices or vectors") def transpose_vector(self): return self.swapaxes(-1,-2) t = property(transpose_vector) class ndcolumn(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): raise TypeError("Cannot matrix multiply columns with anything") def transpose_vector(self): return self.view(ndrow) t = property(transpose_vector) class ndrow(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other, -1).view(ndarray) if isinstance(othe, ndmatrix): raise notimplemented else: raise TypeError("Can only matrix multiply rows with matrices or columns") def transpose_vector(self): return self.view(ndcolumn) t = property(transpose_vector) def matrix(*args, **kwargs): m = array(*args, **kwargs) return m.view(ndmatrix) def column(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndcolumn) def row(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndrow) if __name__ == "__main__": import doctest, matrix doctest.testmod(matrix)
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/30/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
That should be allowed. (N,1)*(1,M) is just an (N,M) matrix with entries C[i,j] = A[i,0]*B[0,] I kind of like the idea of using call for multiply, though. If it doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot". --bb
![](https://secure.gravatar.com/avatar/5b2449484c19f8e037c5d9c71e429508.jpg?s=120&d=mm&r=g)
On 3/29/07, Bill Baxter <wbaxter@gmail.com> wrote:
I thought about that a little, and while I agree that it could be allowed, I'm not sure that it should be allowed. It's a trade off between a bit of what I would guess is little used functionality with some enhanced error checking (I would guess that usually row*column signals a mistake). However, I don't care much one way or the other; it's not hard to allow. I kind of like the idea of using call for multiply, though. If it
doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot".
We'll see how it goes down this time. I've proposed using call before, since I've thought the matrix situation was kind of silly for what seems like ages now, but it always sinks without a ripple. -- //=][=\\ tim.hochberg@ieee.org
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
It's really a sort of tensor product, so use outer(.,.). In my mind, row and column vectors are *not* matrices, they only have a single dimension. On the other hand (r)(c) is really the application of the dual vector r (a functional) to the vector c, i.e., r is a map from vectors into the reals (complex). However, I think overloading the multiplication in this case is reasonable. I kind of like the idea of using call for multiply, though. If it
doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot".
Hmmm, have to try it a bit to see how it looks. Overall, I like this approach. Chuck
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
Looks promising!
I'm also for allowing it; from the perspective of one writing code that "looks like" typical (if potentially technically incorrect) notation, being able to write direct analogs to a'b (to get a scalar) or ab' (to get an MxN matrix) and have everything "work" would be quite useful. Especially in various reconstruction tasks (e.g. using svd to approximate a given matrix with a lower-rank one), the "outer product" of a row and a column vector comes up often enough that I would be loathe to have it raise an error.
The call syntax is nice, if a bit opaque to someone looking at the code for the first time. On the other hand, it's explicit in a way that overloaded multiplication just isn't. I like it (for what little that's worth). Zach
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/30/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
It's useful for many things. You can use it in the computation of the least squares best fits between sets of points. The sum of p_i p_^t comes up in that context. (Which is also the covariance matrix of the points) The derivative of a unit vector (useful for mass spring dynamics calcs, among other things, I'm sure) is given by something like c * (I - x x^t). Householder reflections are useful for a lot of things such as implementing QR factorization. They're given by : H = I - 2 x x^t / x^t x http://www.cs.ut.ee/~toomas_l/linalg/lin2/node6.html I'm positive I've seen col * row come up in other contexts too, though I can't think of any other particulars right now. --bb
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
<snip>
Note, however that you can't (for instance) multiply column vector with a row vector:
Well, let's carry this to extremes. Strictly speaking, from the functional point of view, (r)(c) == (c)(r), the dual of the dual is the original vector space, for finite dimensional spaces anyway. However, custom reserves the latter form for the tensor product, and so that is probably a good thing to keep. However, to really make the transpose operate correctly it should also conjugate. This could just be a flag in the row vector, no changes in the data needed, but vdot would be called instead of dot when computing (c.t)(c). A similar thing could be done for matrices if vdot were extended to deal with matrices -- currently it only works for scalars and vectors. Now, products like (c.t)(M.t) would need to conjugate/transpose both, but this could be computed as ((M)(c)).t. So on and so forth. A lot of these types of operations are in BLAS, IIRC, so in some sense this is just a mapping to BLAS. Chuck
![](https://secure.gravatar.com/avatar/95198572b00e5fbcd97fb5315215bf7a.jpg?s=120&d=mm&r=g)
On 3/27/07, Alan G Isaac <aisaac@american.edu> wrote:
I'm personally not really partial to any side of this discussion, given how I don't use matrices at all (I'm content with a simple mental model of arrays, and use dot() as necessary). But it's probably worth mentioning, as a data point, that the python language has a well established precedent for 'iteration over a FOO where a FOO is yielded': In [1]: s='abc' In [2]: type(s) Out[2]: <type 'str'> In [3]: map(type,s) Out[3]: [<type 'str'>, <type 'str'>, <type 'str'>] I know strings aren't matrices, so take this as you will in terms of being plus, neutral or minus regarding your discussion. I'm just offering the data, not interpreting it :) Cheers, f
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Robert Kern apparently wrote:
Gram-Schmidt orthogonalization
Robert wrote:
It is really not my intent to be irritating or to move the posts. But I had two reasons to ask. First, your claim is false that this was important in your example. (To see this, replace M[1:] with M.A[1:].) So actually your code would work the same if iteration over matrices were producing arrays (after replacing M[0] with M[0,:] to match that case). Second, when I spoke of *desirability*, the output is relevant. Try nump.mat(ortho) to see what I mean. If iteration were to produce arrays, the outcome of implementing the algorithm (e.g., using ``dot``) would be I suggest more "desirable", in that numpy.mat(ortho) would work as hoped/expected. In this sense, the example perhaps speaks against your intent. You offer code that would work just fine if iteration yielded arrays. Apologies in advance if this again seems tangential to my original request. To me it does not.
If you really want iterates to be arrays, just use .A.
I have mentioned that several times. It is orthogonal to the design issue I raised. (I suggested that under the current behavior you do not usually end up with what you want when iterating over matrices.) But again, I recognize that I am alone in that view. I am just trying to understand why. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
You're right, it does, but only because the first place it is used, it gets operated with the matrix row-vector base before it does anything else, gets coerced into a matrix, and then gets assigned back to rowv. I don't tolerate that kind of implicitness in my code.
Second, when I spoke of *desirability*, the output is relevant. Try nump.mat(ortho) to see what I mean.
<shrug> If I wanted something other than a list of matrix row-vectors (which I might), then I'd post-process. Keeping the clarity of the main algorithm is also important.
Ditto for not using matrix objects at all.
People have been giving you reasons, over and over again. You are simply refusing to listen to them. You have a use case for arrays being the iterates. You are presuming that the only argument that can beat that is another use case for matrix objects being the iterates. This is not true; there are other principles at work. If we were to make matrix objects behave in different ways for different methods to satisfy what one person thinks is the most convenient way to use each method, we would have a patchwork object with no overall principle that people can use to understand what is going on. They would simply have to remember which behavior someone thought was the most convenient for each thing. Remember the lessons of rand() and randn(). -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Wed, 28 Mar 2007, Robert Kern wrote:
People have been giving you reasons, over and over again. You are simply refusing to listen to them.
Exploring whether the reasoning is adequate is not the same as refusing to listen. I do not presume my view is correct.
Put slightly differently: given the surprising passion of the attacks at the suggestion that perhaps iteration over a matrix might more consistently yield arrays, I presumed there must be *many* instances in which it was obviously desirable that such iteration should yield matrices. So I asked to see some. In the context of this discussion, I found the (lack of) responses very interesting. Even in your thoughtful response it proved irrelevant rather than important for iteration over matrices to yield matrices. I understand that some people claim that a general principle of consistency is involved. I have not been able to understand this particular design decision as a matter of consistency, and I have tried to say why. However I am just a user (and supporter) of numpy, and as indicated in other posts, I make no pretense of deep insight into the design decisions. In this case, I simply wanted to open a discussion of a design decision, not win an "argument". Anyway, I understand that I am being perceived as bull-headed here, so I'll let this go. Thanks for your attempt to help me see the virtues of the current design. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
It's because the property that A[i] == A[i,...] is much more important to most numpy users than the results of a particular (mis) use of the matrix class. This has been explained in many different contexts over many different email messages by many different people. You're not looking at the big picture, which has been fairly explicitly spelled out by myself and others. I'm not even sure why I'm replying at this point. Zach On Mar 28, 2007, at 3:25 PM, Alan Isaac wrote:
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Wed, Mar 28, 2007 at 05:25:57PM -0500, Alan Isaac wrote:
Matrices strike me as a bit of an anomaly. I would expect an N-dimensional container to contain (N-1)-dimensional objects. But matrices cannot have dimensions less than 2, (and in numpy, can't have dimensions higher than 2). So, to me these operations don't make much sense. Is a matrix then more of a collection of matrices than a container of matrices? I'm glad we have the ndarray concept in numpy where these concepts consistently make sense. Cheers Stéfan
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Wed, Mar 28, 2007 at 07:05:00PM -0500, Alan Isaac wrote:
Doesn't seem to be the way the matrix world works though: octave:2> x = zeros(3,3,5); octave:3> size(x) ans = 3 3 5 octave:4> size(x(:,:,1)) ans = 3 3 octave:5> size(x(:,1,1)) ans = 3 1 Cheers Stéfan
![](https://secure.gravatar.com/avatar/37ee74afc8985c99adac9eb37ab2b10e.jpg?s=120&d=mm&r=g)
On 3/27/07, Zachary Pincus <zpincus@stanford.edu> wrote:
I find it useful if M[i] indexes the matrix interpreted as a column-vector (vec(M), or M(:) in Matlab notation), e.g., you could pick out the diagonal as M[::n+1], or if N is an index-set of certain elements in M, then you could access those elements as M[N]. Then I would also say that iterating over a matrix should just return the elements of vec(M) one by one. In general, I think the Matlab notation works well: * M[I,J] where I and J and index-sets returns a len(I) x len(J) matrix * M[I] returns a len(I) vector (a column or row vector depending on orientation of I) The index-sets could be a single integer, a list, or a slice. But I can see there are good arguments for using the NumPy conventions as well...
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Hi Zachary, I think your response highlights very well the apparent design flaw. Here is your response to my request for a use case where iteration over a matrix should yield matrices: do not iterate! Does some part of you not find that just a little bizarre?! As for offering as a use case that with a matrix M you are allowed to use M[0] instead of M[0,:] when you want the first row as a matrix, I really cannot take that trivial convenience seriously as the basis of a fundamental design decision. However there is no sympathy for my view on this, so I am not debating it any more. Instead I have asked a simple question: what use case is this design decision supporting? I am interested in this so that I can see into the decision better. Nowith Chris proposes that "M[i] == M[i,:] by definition". If so, that is an design-based answer to my question. I agree that M[i,:] should return a matrix. But my understanding was different: I thought M[i] relied on standard Python idexing while M[i,:] was numpy indexing. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps our differences lies in two things: 1. the fact that the text books typically take the column vector as the default. For a Python version, based on C it makes more sense to treat the rows as vectors, as data is stored contiguously by row. 2. the convention has been proposed that the vector is more conveniently implemented as a matrix, where one dimension is one. The vector could be treated as a subclass of the matrix but this adds complexity with little clear benefit. PyMatrix has matrix methods isVector, isCVector and isRVector. I can see some merit in conforming to text book usage and would be glad to consider changes when I complete the port to numpy, in a few months. Colin W.
Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Sun, 25 Mar 2007, "Colin J. Williams" apparently wrote:
An array and a matrix are different animals. Conversion from one to the other should be spelled out.
But you are just begging the question here. The question is: when I iterate across matrix rows, why am I iterating across matrices and not arrays. This seems quite out of whack with general Python practice. You cannot just say "conversion should be explicit" because that assumes (incorrectly actually) that the rows are matrices. The "conversion should be explicit" argument actually cuts in the opposite direction of what you appear to believe. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Alan, Yes, this is where we appear to differ. I believe that vectors are best represented as matrices, with a shape of (1, n) or (m, 1). The choice of these determines whether we have a column or a row vectors. Thus any (m, n) matrix can be considered as either a collection of column vectors or a collection of row vectors. If the end result is required as an array or a list, this can be done explicitly with X[1].A or A[1].tolist(). Here, A is a property of the M (matrix) class.
Cheers, Alan Isaac
A long time ago, you proposed that PyMatrix should provide for matrix division in two way, as is done in MatLab. This was implemented, but PyMatrix has not yet been ported to numpy - perhaps this summer. Regards, Colin W.
![](https://secure.gravatar.com/avatar/37c1dd52b545c64455f1de5977d16dff.jpg?s=120&d=mm&r=g)
Em Dom, 2007-03-25 às 13:07 -0400, Alan G Isaac escreveu:
I think the point here is that if you are using matrices, then all you "should" want are matrices, just like in MATLAB: >> A = [1 2; 3 4] A = 1 2 3 4 >> b = A(1, :) b = 1 2 >> size(b) ans = 1 2 >> b = A(:, 1) b = 1 3 >> size(b) ans = 2 1 >> b = 1 b = 1 >> size(b) ans = 1 1 You see, rows, columnes, and even numbers, are treated as matrices. Paulo -- Paulo José da Silva e Silva Professor Assistente do Dep. de Ciência da Computação (Assistant Professor of the Computer Science Dept.) Universidade de São Paulo - Brazil e-mail: rsilva@ime.usp.br Web: http://www.ime.usp.br/~rsilva Teoria é o que não entendemos o (Theory is something we don't) suficiente para chamar de prática. (understand well enough to call practice)
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Sun, 25 Mar 2007, Paulo Jose da Silva e Silva apparently wrote:
Yes, that is the idea behind this, which I am also accustomed to from GAUSS. But note again that the Matlab equivalent :: >>> x=N.mat('1 2;3 4') >>> x[0,:] matrix([[1, 2]]) does provide this behavior. The question I am raising is a design question and is I think really not addressed by the rule of thumb you offer. Specifically, that rule of thumb if it is indeed the justification of :: >>> x[1] matrix([[3, 4]]) finds itself in basic conflict with the idea that I ought to be able to iterate over the objects in an iterable container. I mean really, does this not "feel" wrong? :: >>> for item in x: print item.__repr__() ... matrix([[1, 2]]) matrix([[3, 4]]) Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
This may sound silly, but I really think seeing all those brackets is what makes it feel wrong. Matlab's output doesn't put it in your face that your 4 is really a matrix([[4]]), even though that's what it is to Matlab. But I don't see a good way to change that behavior. The other thing I find problematic about matrices is the inability to go higher than 2d. To me that means that it's impossible to go "pure matrix" in my code because I'll have to switch back to arrays any time I want more than 2d (or use a mixed solution like a list of matrices). Matlab allows allows >2D. --bb
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Colin J. Williams <cjw@sympatico.ca> wrote:
I'm not sure what you thought I meant, but all I meant by going "pure matrix" was having my Numpy code use the 'matrix' type exclusively instead of some mix of 'matrix' and the base 'ndarray' type. Things become messy when you mix and match them because you don't know any more if an expression like A[1] is going to give you a 1-D thing or a 2-D thing, and you can't be sure what A * B will do without always coercing A and B.
A list of matrices seems to be a logical structure.
Yes, and it's the only option if you want to make a list of matrices of different shapes, but I frequently have a need for things like a list of per-point transformation matrices. Each column from each of those matrices can be thought of as a vector. Sometimes its convenient to consider all the X basis vectors together, for instance, which is a simple and efficient M[:,:,0] slice if I have all the data in a 3-D array, but it's a slow list comprehension plus matrix constructor if I have the matrices in a list -- something like matrix([m[:,0] for m in M]) but that line is probably incorrect.
Numpy generally tries to treat all lists and tuples as array literals. That's not likely to change. --bb
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Bill Baxter wrote:
Yes, to my mind it's best to consider the multi-dimensional array and the matrix to be two distinct data types. In most cases, it's best that conversions between the two should be explicit.
Logically, this makes sense, where M is a list of matrices. My guess is that it would be a little faster to build one larger matrix and then slice it as needed.
That need no be a problem is there is clarity of thinking about the essential difference between the matrix data type (even if is is built as a sub-type of the array) and the multi-dimensional array.
--bb
Colin W.
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
I mean really, does this not "feel" wrong? ::
On Mon, 26 Mar 2007, Bill Baxter apparently wrote:
This may sound silly, but I really think seeing all those brackets is what makes it feel wrong.
I appreciate the agreement that it feels wrong, but I dispute the analysis of this symptom. What makes it "feel wrong", I contend, is that experience with Python make this a **surprising** behavior. And that is precisely why I suggest that this may point to a design issue. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
So you're saying this is what you'd find more pythonic?
Probably about half the bugs I get from mixing and matching matrix and array are things like row = A[i] ... z = row[2] Which works for an array but not for a matrix. I think Matlab makes it more bearable by having a single value index like X[i] be equivalent to X.flat()[i]. So X[2] is the same for row or col vec in Matlab. Now that I think about it, that's probably the main reason I feel more comfortable with array than matrix in Numpy. If I have a vector, I should only need one index to get the ith component. --bb
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Bill Baxter apparently wrote:
No; that is not possible, since a matrix is inherently 2d. I just want to get the constituent arrays when I iterate over the matrix object or use regular Python indexing, but a matrix when I use matrix/array indexing. That is :: >>> X[1] array([2,3]) >>> X[1,:] matrix([[3, 4]]) That behavior seems completely natural and unsurprising.
Exactly! That is the evidence of a "bad surprise" in the current behavior. Iterating over a Python iterable should provide access to the contained objects. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Oooops, they should match of course. :: >>> X[1] array([3,4]) >>> X[1,:] matrix([[3, 4]]) But again the point is: indexing for submatrices should produce matrices. Normal Python indexing should access the constituent arrays. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Alan G Isaac schrieb:
I think this is a tricky business. There's also the rule "indexing reduces the rank, slicing preserves it". Numpy-matrices form an exception to this rule, always being 2d, but the exception is consistently enforced. Now it seems you want to have an exception from the exception, correct? Isn't this why the .A1 method for numpy-matrices was introduced even after 1.0rc? -sven (matrix fan)
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sven Schreiber apparently wrote:
What I want: the best design. I did not claim that the current design is flawed, only to suspect it. Why I wrote: current behavior feels wrong -> suspect design flaw. What feels wrong: iterating over a container does not give access to the contained objects. This is not Pythonic. *Symptom* of the underlying problem: for matrix M, M[0] returns a matrix. Would the change I suggest mean that the behavior of the matrix class deviates less from the array class: yes. Does this mean the matrix class behavior would be less "consistent"? I have tried to explain why the answer is "no". hth, Alan Isaac
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Alan G Isaac schrieb:
What feels wrong: iterating over a container does not give access to the contained objects. This is not Pythonic.
If you iterate over the rows of the matrix, it feels natural to me to get the row vectors -- and as you know a 1d-array does not contain the information of "row-ness", so you need to get a 2d thing to properly return those "contained objects". So in my view the current behavior actually does exactly what you're calling for. (But I admit I'm a 2d fan so much so that I didn't even know that using a single index is possible with a matrix. I thought one always had to be explicit about both dimensions... So thanks for pointing this out. -- BTW, I'm aware that sticking to numpy-matrices exclusively is not possible. For example, eigenvalues are returned in a 1d-array even if you pass a matrix, and it's intended behavior.) cheers, sven
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sven Schreiber apparently wrote:
If you iterate over the rows of the matrix, it feels natural to me to get the row vectors
Natural in what way? Again, I am raising the question of what would be expected from someone familiar with Python. Abstractly, what do you expect to get when you iterate over a container? Seriously.
But I admit I'm a 2d fan so much so that I didn't even know that using a single index is possible with a matrix.
Exactly. When you want submatrices, you expect to index for them. EXACTLY!! Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/6194b135cba546afa82516de1537de49.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
If may chime in... I think Sven's argument in on the side saying, A "matrix" is an object that you expect a certain (mathematical !) behavior from. If some object behaves intuitively right -- that's ultimately pythonic ! The clash is, NOT to see a matrix "just as another container". Instead a matrix is a mathematical object , that has rows and columns. It is used in a field (lin-alg) where every vector is either a row or a column vector -- apparently that's big thing ;-) The whole reason to add a special matrix class to numpy in the first place, is to provide a better degree of convenience to lin-alg related applications. I would argue that it was just not consistently considered, that this should also come with "a column of a matrix is something else than a row -- (1,n) vs. (n,1) and not (n,) more notes/points: a) I have never heard about the m.A1 - what is it ? b) I don't think that if m[1] would return a (rank 2) matrix, that m[1].A could return a (rank 1) array ... c) I'm curious if there is a unique way to extend the matrix class into 3D or ND. -Sebastian
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Sebastian Haase apparently wrote:
The problem is, as I am not the only one to point out, this particular behavior is NOT intuitively right.
The clash is, NOT to see a matrix "just as another container".
But be serious, no object is "just another container". Again, this just begs the question. The question is a design question. E.g., what is the principle of least surprise?
more notes/points: a) I have never heard about the m.A1 - what is it ?
It returns a 1d array holding the raveled matrix.
b) I don't think that if m[1] would return a (rank 2) matrix, that m[1].A could return a (rank 1) array ...
It does not, of course. (But both should, I believe.)
c) I'm curious if there is a unique way to extend the matrix class into 3D or ND.
Is that not what an array is for?? Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps things would be clearer if we thought of the constituent groups of data in a matrix as being themselves matrices. X[1] could represent the second row of a matrix. A row of a matrix is a row vector, a special case of a matrix. To get an array, I suggest that an explicit conversion X[1].A is a clearer way to handle things. Similarly, X[2, 3] is best returned as a value which is of a Python type. Colin W.
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps it would now help if you redefined the question. In an earlier posting, you appeared anxious that the matrix and the array behave in the same way. Since they are different animals, I see sameness of behaviour as being lower on the list of desirables than fitting the standard ideas of matrix algebra. Suppose that a is a row vector, b a column vector and A a conforming matrix then: a * A A * b and b.T * A are all acceptable operations. One would expect the iteration over A to return row vectors, represented by (1, n) matrices. Colin W.
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, "Colin J. Williams" apparently wrote:
One would expect the iteration over A to return row vectors, represented by (1, n) matrices.
This is again simple assertion. **Why** would "one" expect this? Some people clearly do not. One person commented that this unexpected behavior was a source of error in their code. Another person commented that they did not even guess that such a thing would be possible. Experience with Python should lead to the ability to anticipate the outcome. Apparently this is not the case. That suggests a design problem. What about **Python** would lead us to expect this behavior?? In *contrast*, everyone agrees that for a matrix M, we should get a matrix from M[0,:]. This is expected and desirable. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
Well, and what *should* they expect. I think it is expected because the for iterates over rows and rows of matrices are 1xn. Matrices and arrays, as has been stated, are different animals. Probably it would have been best to stick with arrays and I suspect that matrices appeared because of the dearth of Python operators, in particular to make matrix multiplication simpler. On the other hand, certain errors slip by when one is implementing matrix algebra with arrays, but they can be avoided by never using 1-d vectors. So all this mess results from catering to the matrix community. Matlab has the opposite problem, multidimensional arrays were tacked on later and they don't operate properly with everything. Chuck One person commented that this unexpected behavior was
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Charles R Harris apparently wrote:
Well, and what should they expect.
Since matrices are an iterable Python object, we *expect* to iterate over the contained objects. (Arrays.) I am not sure why this is not evident to all, but it is surely the sticking point in this discussion. A matrix is not a container of matrices. That it acts like one is surprsing. Surprises are bad unless they have a clear justification. Perhaps a clear justification exists, but it has not been offered in this discussion. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
I think that a clear justification has been offered several times on the list recently, though maybe not all in this thread. Matrices in numpy seem to exist as a construct to facilitate linear algebra. One such property is that row- and column-vector slices of matrices are (1xN) or (Nx1) matrices, because otherwise common linear algebra things -- like how you multiply a row-vector or a column vector by a matrix, and whether and when it needs to be transposed -- do not translate properly from "linear algebra notation" like in textbooks and papers. Once the matrix class is committed to providing row- and column- vector slices as other matrices, it makes no sense to have iteration over the matrix provide a different data type than slicing. Basically, my rule of thumb has been to *only* use matrices when I'm copying linear algebra operations out of textbooks and papers, and I want the notations to align. Doing other, non-linear-algebra operations with matrices -- like iterating over their elements -- is not really worth it. There's a meta question, of course: should things be changed to make it "worth it" to do "pythonic" tasks with matrices? Should there be special elementwise vs. matrix-operation operators? Should numpy work a lot more like matlab? That discussion is on-going in another thread, I think. Zach
![](https://secure.gravatar.com/avatar/4d021a1d1319f36ad861ebef0eb5ba44.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
It actually has been offered. You just don't accept it. Matrices are containers of matrices. If M is an (mxn) matrix then M[0] is a (1xn) matrix. Viewing this 1xn matrix as a 1-d array loses it's row-vectorness. This seems perfectly logical and acceptable to me. I'm waiting for a better explanation as to why this is not acceptable. Arguments that rest on what is and what isn't "Pythonic" are not very persuasive as this is very often in the eye of the beholder. -Travis
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Travis Oliphant wrote:
Again, I only raised a *question* about whether there might be a design problem here. My goal was only to have this explored, and I've tried to explain why. The evidence I offer: - it is surprising (at least to some) that iterating over a matrix yields matrices - I believe it is unexpected (prior to instruction) and that there is a natural more expected behavior - if that is right, deviation from the expected should have a good justification - this behavior has tripped up at least a couple people and I expect that to happen to many others over time (because I believe the behavior is unexpected) - when I desire to iterate over a matrix I always want the arrays. (Of course there is a way to have them; that's not the point). I'm interested to see a use case where the rows are desired as matrices As you indicate, none of this constitutes an "argument". And since nobody appears to agree, I should shut up. This will be my last post on this subject. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/4d021a1d1319f36ad861ebef0eb5ba44.jpg?s=120&d=mm&r=g)
Alan Isaac wrote:
Thanks for listing your points. I can see that this is an issue where reasonable people will disagree because there are multiple ways of looking at it. The idea that matrix selection would return matrices really comes from wanting to keep matrices more persistent through operations. M[0] could be seen either as a 1xn matrix or a n-length array. I agree that both concepts are possible. Seeing it as a 1xn matrix allows matrices to remain persistent more often. So, the arguments for the current approach and the arguments against it to me seem on the same level, so I don't see a reason to change the current behavior and see a lot of strong reasons not to change the behavior (we are in a 1.0 release and could not change anything until at least 1.1 anyway). With that said: One of my goals for the next year or two is to create a matrix class in C and incorporate CVXOPT matrices and it's coupled sparse matrix. We can re-visit this question in that process. I would like for there to be a sparse matrix implementation in NumPy (without the solver which will remain in SciPy). But, the sparse matrix and the matrix need to have the same behaviors and should be able to interoperate with each other. So, if you would like to help with that project all input is welcome. Best regards, -Travis
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
It seems to me that using shapes, (m,1) versus (1,n), to determine whether a vector is column- or row-oriented is a hack (or at least feels like one). If we know we have a vector, then we want to use a single index to obtain a scalar, and that extra "0," or ",0" shouldn't be necessary. It also seems to me that there is an object-oriented design solution to this, namely defining column_vector and row_vector classes that inherit from matrix, that accept a single index, but behave as matrices consistent with their type. I'm sure there are implications to this I haven't thought through yet, but at least it gives us logical indexing AND persistence of matrices through operations. On Mar 27, 2007, at 1:11 PM, Travis Oliphant wrote:
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/26/07, Alan G Isaac <aisaac@american.edu> wrote:
of the tensor product of a vector space with its dual, so one would expect either row or column vectors depending on the index position. But this discussion is not really relevant at this point. What happens is a convention and the numpy convention is out there. It isn't going to change now, it would break too much code. Chuck
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Mon, 26 Mar 2007, Charles R Harris wrote:
What happens is a convention
Certainly true.
It isn't going to change now, it would break too much code.
That is a different kind of argument. It might be true. May I see a use case where the desired return when iterating through a matrix is rows as matrices? That has never been what I wanted. Thank you, Alan Isaac
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Bill Baxter wrote:
This seems artificial to me for several reasons, but here is one reason: AllxformedPts = AllMyPoints * TransformMatrix Note that I am no longer disputing the convention, just trying to understand its usefulness. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/27/07, Alan Isaac <aisaac@american.edu> wrote:
Yeh, I was afraid you'd say that. :-) But maybe I've got a lot of points, and I don't feel like making a copy of the whole set. Or maybe it's not a linear transform, but instead xformedPt = someComplicatedNonLinearThing(pt) I do stuff like the above quite frequently in my code, although with arrays rather than matrices. :-) For instance in finite elements there's assembling the global stiffness matrix step where for each node (point) in your mesh you set some entries in the big matrix K. Something like for i,pt in enumerate(points): K[shape_fn_indices(i)] = stiffness_fn(pt) That's cartoon code, but I think you get the idea. I don't think there's any good way to make that into one vectorized expression. The indices of K that get set depend on the connectivity of your mesh.
Note that I am no longer disputing the convention, just trying to understand its usefulness.
--bb
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Bill Baxter apparently wrote:
xformedPt = someComplicatedNonLinearThing(pt)
I do stuff like the above quite frequently in my code, although with arrays rather than matrices.
Exactly: that was one other thing I found artificial. Surely the points will then be wanted as arrays. So my view is that we still do not have a use case for wanting matrices yielded when iterating across rows of a matrix. Alan Isaac PS In the "large number of points" case, the thing to do would be to extract modest sized submatrices.
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
It's pretty clear from my perspective: 1-D slices of matrices *must* be matrices. I gave an intuitive make-it-fit-with-known-notation explanation, and Charles gave a theoretically-grounded explanation. There was no objection to this from any quarter. Given that 1-D slices of matrices must be matrices for deep mathematical reasons as well as notational compatibility with the rest of linear algebra -- e.g. that m[0] yields a matrix if m is a matrix-- it almost certainly would violate the principle of least surprise for iteration over m (intuitively understood to be choosing m [0] then m[1] and so forth) to yield anything other than a matrix. This can't possibly be what you're asking for, right? You aren't suggesting that m[0] and list(iter(m))[0] should be different types? There are many clear and definite use cases for m[0] being a matrix, by the way, in addition to the theoretical arguments -- these aren't hard to come by. Whether or nor there are actual good use-cases for iterating over a matrix, if you believe that m[0] should be a matrix, it's madness to suggest that list(iter(m))[0] should be otherwise. My opinion? Don't iterate over matrices. Use matrices for linear algebra. There's no "iteration" construct in linear algebra. The behavior you find so puzzling is a necessary by-product of the fundamental behavior of the matrix class -- which has been explained and which you offered no resistance to. Zach
![](https://secure.gravatar.com/avatar/5dde29b54a3f1b76b2541d0a4a9b232c.jpg?s=120&d=mm&r=g)
Zachary Pincus wrote:
I don't think the OP was suggesting that. Rather, I think the suggestion was that for a mXn matrix, M: M[i].shape == (n,) M[i:i+1].shape == (1,n) that is, indexing (or iterating returns a vector, and slicing returns a matrix). This is, indeed exactly how numpy arrays behave! The problem with this is: numpy matrices were created specifically to support linear algebra calculations. For linear algebra, the distinction between row vectors and column vectors is critical. By definition, a row vector has shape: (1,n), and a column vector has shape (m,1). In this case, perhaps the OP is thinking that a shape (n,) array could be considered a row vector, but in that case: M[1,:].shape == (n,) M[:,1].shape == (m,) which of these is the row and which the column? This is why matrices index this way instead: M[1,:].shape == (1, n) M[:,1].shape == (m, 1) now we know exactly what is a row and what is a column. By the way, I think with the way numpy works: M[i] == M[i,:] by definition, so those couldn't yield different shaped results. Is that right? I think we got a bit sidetracked by the example given. If I have a bunch of points I want to store, I'm going to use an (m,2) *array*, not a matrix, then then A[i] will yield a (2,) array, which makes sense for (2-d) points. In fact, I do this a LOT. If I happen to need to do some linear algebra on that array of points, I'll convert to a matrix, do the linear algebra, then convert back to an a array (or just use the linear algebra functions on the array). I hope this helps -Chris
-- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
![](https://secure.gravatar.com/avatar/37ee74afc8985c99adac9eb37ab2b10e.jpg?s=120&d=mm&r=g)
I would like to hear your opinion on developing an explicit sparse/dense 2D matrix class with indexing similar to Matlab, and without significant differences between sparse and dense matrices from the user's perspective... I know that it's not one of Numpy/Scipy's goals to clone Matlab, but I think I represent a potentially large scientific Python user base, who would find Python matrices that feel a bit more like Matlab/Octave etc. extremely useful. I have a great respect for all the work in Numpy and Scipy, but at the same time I feel that numerical linear algebra in Python would benefit from a more dedicated matrix library, and I think that Numpy (or another single package) should provide that in a homogeneous way - without ties to how Numpy array works. - Joachim On 3/27/07, Christopher Barker <Chris.Barker@noaa.gov> wrote:
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
Hello all, I suspect my previous email did not contain the full chain of my reasoning, because I thought that some parts were basically obvious. My point was merely that given some pretty basic fundamental tenants of numpy, Alan's suggestions quickly lead to *serious issues* far worse than the original problem. Thanks Chris for making this more explicit; here I will do so further. Let A be an array of arbitrary shape, and M be m-by-n matrix. (1) A[i] equals by definition A[i,...] This is a pretty fundamental part of numpy, right? That property (and its relatives) along with the broadcasting behavior, are in large part what make numpy so very expressive compared to matlab, etc. I can write code in numpy to do multidimensional operations over arbitrary-shaped and arbitrary-dimensioned data that look *exactly* like the corresponding 1D operations. It's amazing, really, how numpy's "implicit indexing" (this property), broadcasting (sort of the inverse of this, really, where shape tuples are extended as necessary with implicit 1's, instead of slices being extended with :'s as necessary), and fancy indexing, allow complex operations to be written compactly, clearly, and in a fully-general manner. I thought that point (1) was obvious, as it's a pretty fundamental part of numpy (and Numeric and Numarray before that). It's in absolutely no way a "trivial convenience". Of course, broadcasting and ufuncs, etc, contribute a lot more to expressiveness than this property, but it's part of the family. (2) Thus, M[i] equals by definition M[i,:]. By the mathematical definition of matrices and for pragmatic reasons, M[i,:] is another matrix. Here the trouble starts. Should matrices, not being arbitrary in their dimension, even have the ability to be indexed as M[i]? Or should they be completely different beasts than arrays? Which violates the "principle of least surprise" more? Now, I think we all more or less agree that M[i,:] has to be a matrix, because otherwise the matrix class is pretty worthless for linear algebra, and doesn't operate according to the mathematical definition of a matrix. Indeed, this is one of the basic properties of the numpy matrix -- that, and that the multiplication operator is defined as matrix-multiplication. Remove that and there's almost no point in *having* a separate matrix class. Want a use-case? Do the singular value decomposition, and get the product of the second-smallest singular vectors as Vt[:,-2] * U [-2,:]. (Things vaguely akin to this come up a lot.) You can copy the operation from any linear algebra text if the row- and column-vectors are treated as above. Otherwise you have to wonder whether that's supposed to be an outer or inner product, and use the corresponding operator -- and at that point, why was I using matrix classes anyway? (3) If M[i] = M[i,:], and M[i,:] is a matrix, then the iteration behavior of matrices *has* to yield other matrices, unless Alan is willing to suggest that list(iter(m))[i] should have a different type than M[i] -- a clear absurdity. This is the point that I was trying to make previously, having mistakenly assumed that point (1) was clear to all, and that point (2) had been clearly established. So, if we trace the reasoning of Alan's suggestion, coupled with these above properties, we get just that: (a) Iteration over M should yield arrays. (b) Ergo M[i] should be an array. (c) Ergo M[i,:] should be an array. (d) Ergo the matrix class should be identical to the array class except for overloaded multiplication, and the only way to get a 'row' or 'column' vector from a matrix that operates as a proper row or column vector should is M[i:i+1,:]. Whicy is a lot of typing to get the 'i-th' vector from the matrix. I don't think these changes are worth it -- it basically discards the utility of the matrix class for linear algebra in order to make the matrix class more useful for general purpose data storage (a role already filled by the array type). Here is another set of changes which would make matrices fully consistent with their linear-algebra roots: (a) M[i,:] should be a matrix. (b) M[i] is an error. (c) Iteration over M is an error. That's kind of lame though, right? Because then matrices are completely inconsistent with their numpy roots. So we are left with the current behavior: (a) M[i,:] is a matrix. (b) M[i] is a matrix. (c) Iteration over M yields matrices. My view is that, fundamentally, they are linear algebra constructs. However, I also think it would be far more harmful to remove basic behavior common to the rest of numpy (e.g. that M[i] = M[i,:]), and behavior that comes along with that (iterability). Hence my suggestion: treat matrices like linear algebra constructs -- don't iterate over them (unless you know what you're doing). Don't index them like M[i] (unless you know what you're doing). Maybe it's "just a little bizarre" to suggest this, but I think the other options above are much more bizarre. Zach
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Hi Zach, The use case I requested was for iteration over a matrix where it is desirable that matrices are yielded. That is not what you offered. The context for this request is my own experience: whenever I have needed to iterate over matrices, I have always wanted the arrays. So I am simply interested in seeing an example of the opposite desire. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Gram-Schmidt orthogonalization. ortho = [mat[0] / sqrt(mat[0] * mat[0].T)] for rowv in mat[1:]: for base in ortho: rowv = rowv - base * float(rowv * base.T) rowv = rowv / sqrt(rowv * rowv.T) ortho.append(rowv) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Robert Kern apparently wrote:
Gram-Schmidt orthogonalization
I take it from context that you consider it desirable to end up with a list of matrices? I guess I would find it more natural to work with the arrays, but perhaps that starts just being taste. Thank you, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Honestly, I don't care. You asked about iteration, and I gave you an example where it was important that the iterates were matrices (such that .T worked correctly). Please stop moving the goal posts.
I guess I would find it more natural to work with the arrays, but perhaps that starts just being taste.
Personally, I find working with arrays from start to finish the most natural. If you really want iterates to be arrays, just use .A. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
No, I offered a thorough discussion of why the design of numpy, as I imperfectly understand it, make the trade-offs to achieve your desired goal unpalatable.
I'm glad that Robert provided such a nice example. Nevertheless, I think that even absent any example it should be pretty clear at this point that if you want iteration over matrices to return arrays, either: (1) you need to be willing to accept that list(iter(M))[i] != M[i] or (2) M[i], which equals M[i,:], should be an array. (or 3, that M[i] and iteration over M is not defined) Clearly, (1) is fundamentally ridiculous. Clearly use-cases (as well as the mathematical definition of matrices) abound to illustrate why (2) is no good. Option (3) helps nobody. So which is it? You have to choose one! It is my assertion that if you are iterating over matrices and hoping to get arrays, you are misusing the matrix class. Almost by definition if you expect M[i,:] or M[i] (or equivalently, the output of an iterator over M) to be an array, then M is being used outside of a linear algebra context -- and thus, in my view, misused. In those cases, the matrix should be cast to an array, which has the desired behavior. Robert nicely illustrated a linear-algebraic context for matrix iteration. Now, Bill offers up a different suggestion: indexing M yields neither a matrix nor an array, but a class that operates more or less like an array, except insofar as it interacts with other matrix objects, or other objects of similar classes. I'm interested in hearing more about this, what trade-offs or other compromises it might involve. Zach
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Zachary Pincus wrote:
M[i], which equals M[i,:]
Of course this equality must break. That was stated at the outset. As I said before, this may be impossible or undesirable. But, as I said before, it seems prima facie natural for M[i] to be ordinary Python indexing while M[i,:] is numpy indexing, each with its own natural interpretation. In any case, you may be spending too much energy on this. I am not a developer, and no developer has expressed a desire to make such a change. The current behavior is currently safe (if odd). Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
What I envisioned was that M[i,:] would return a row_vector and M [:,j] would return a column_vector, because this would be symmetric behavior. M[i], by convention, would behave the same as M[i,:]. But then I personally don't distinguish between "python indexing" and "numpy indexing". In both cases, __getitem__() (or __setitem__()) is called. For multiple indexes, the index object is a tuple. In any case, the behavior of "numpy indexing" as I have proposed it would return an object that inherits from matrix, thus would BE a matrix, although it would behave like an array. On Mar 29, 2007, at 3:24 PM, Alan G Isaac wrote:
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Bill Spotz <wfspotz@sandia.gov> wrote:
I'm thinking that a basic problem is not having row and column types distinct from matrices. Column types would represent elements of a vector space and row types elements of the dual, they would both be 1-dimensional. One could go full bore with tensors and index signatures, upper and lower, but I think plain old matrices with the two vector types would solve most ambiguities. Thus if v were a column vector, then v.T*v would be a scalar where the product in this case is shorthand for v.T(v). Likewise, v*v.Twould be a matrix (in tensor form the indices would be upper lower respectively, but ignore that). The default translation of a normal 1-D vector would preferably be a column vector in this case, opposite to current usage, but that is not really crucial. Note that nx1 and 1xn matrices are not quite the same thing as column and row vectors, as the latter would normally be considered linear maps and v.T*v would return a scalar matrix, quite a different thing from a scalar. Chuck
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Thu, 29 Mar 2007, Bill Spotz apparently wrote:
Can you please be explicit about the value added by that convention, as you see it? Thank you, Alan Isaac PS I assume your "vector" class would always accept a single index, for both row and column vectors? What is the return type of v[i]?
![](https://secure.gravatar.com/avatar/8a43a6806c8654a100a8cfb8e90fa73b.jpg?s=120&d=mm&r=g)
On Mar 29, 2007, at 6:48 PM, Alan G Isaac wrote:
I assume by "convention" you are referring to the convention that M [i] is M[i,:] is a row_vector. I see it as a design decision: Does M[i] mean something? No: Then raise an exception. If users want to iterate over M, this forces them to use the double index notation to specify what they are iterating over. I would find this acceptable. Yes: What should M[i] represent/return? Array representing row vector: Array representing column vector: I would argue against linear algebra objects returning raw arrays because their meaning/interface can be ambiguous, as this thread demonstrates row_vector: Reasonable choice, analogous to conventional mathematical notation column_vector: Reasonable choice, but less analogous to conventional mathematical notation Of all the choices, raising an exception and returning a row_vector are (in my opinion, anyway) the best. Unless someone has a better argument.
My best guess at a best design would be that it would return a scalar. And yes, vectors would always accept a single index. But because they are matrices, they would accept two indexes as well.
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
![](https://secure.gravatar.com/avatar/5b2449484c19f8e037c5d9c71e429508.jpg?s=120&d=mm&r=g)
"""matrix.py The discussion about matrix indexing has been interminible and for the most part pretty pointless IMO. However, it does point out one thing: the interaction between the matrix and array classes is still pretty klunky despite a fair amount of effort trying to make them interoperate. Here's one possible alternative approach to enabling matrix operations within the context of numpy arrays. I started with a few requirements: 1. Matrices should be proper subclasses of arrays in the sense that one should be able to use matrices wherever one uses arrays with no change in the result. 2. Indexing into matrices should produce row or column vectors where appropriate. 3. There should be some syntax for matrix multiplication. I know that there are other operations defined on the matrix class, but I don't consider them worth the headache of having two class hierarchies. This is in part inspired by some comments Bill Spotz, but after listening to his later comments, I don't think that this is really what he had in mind. So, please don't blame him for any deficiencies you find herein. The usual caveats apply: this is just a sketch, it's probably buggy, use at your own risk, etc. However, you should be able to save this mail as "matrix.py" and execute it to run all of the examples below. I enforce point #1 by not overriding any ndarray methods except for __array_finalize__ and __getitem__ and they are overridden in a way that shouldn't effect code unless it is specifically testing for types. With that summarily taken care of, let's move onto point #2: indexing. First let's create a matrix. Scratch that, let's create an array of matrices. One things that this approach does is let you create arrays of matrices and vectors in a natural way. The matrix below has a shape of (3,2,2), which correspons to an array of 3 2x2 matrices.
m = matrix([ [[1, 0],[0,1]], [[2, 0],[0,2]], [[3, 0],[0,3]] ])
If we index on the first axis, we are selecting a single matrix, so we would expect an ndmatrix back.
On the other hand, if we index along the second axis, we are indexing on the matrices columns and thus expect to get an array of row matrices.
Finally if we index on the last index, we get an array of column matrices.
Indexing of row and column matrices works similarly. It's not shown here and thus will probably turn out to be buggy. For matrix multiplication, I chose to abuse call. Perhaps someday we'll be able to convinve the powers that be to free up another operator, but for the time being this is better than "*" since that forces matrix and array into separate camps and arguably more readable than most of the new syntax proposals that I've seen. Here we multiply our set of three matrices 'm', with a single 2x2 matrix 'n':
Things behave similarly when multiplying a matrix with a row vector or a row vector with a column vector.
Note, however that you can't (for instance) multiply column vector with a row vector:
Finally, you'd like to be able to transpose a matrix or vector. Transposing a vector swaps the type from ndrow and ndcolumn while transposing a matrix transposes the last two axes of the array it is embedded in. Note that I use ".t" for this instead of ".T" to avoid interfering with the existing ".T". (Frankly, if it were up to me, I'd just deprecate .T).
That's about it. Enjoy (or not) as you will. """ from numpy import * _ULT = 1 # Utlimate AKA last _PEN = 2 # Penultimate AKA second to last # I'm sure there's a better way to do this and this is probably buggy... def _reduce(key, ndim): if isinstance(key, tuple): m = len(key) if m == 0: return 0 elif m == 1: key = key[0] else: if Ellipsis in key: key = (None,)*(ndim-m) + key if len(key) == ndim: return isinstance(key[-1], int) * _ULT + isinstance(key[-2], int) * _PEN if len(key) == ndim-1: return isinstance(key[-1], int) * _PEN if isinstance(key, int): if ndim == 1: return _ULT if ndim == 2: return _PEN return 0 class ndmatrix(ndarray): def __array_finalize__(self, obj): if obj.ndim < 2: raise ValueError("matrices must have dimension of at least 2") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT | _PEN: return value.view(ndarray) if rmask == _ULT: return value.view(ndcolumn) if rmask == _PEN: return value.view(ndrow) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other[...,newaxis,:], -1).view(ndcolumn) if isinstance(other, ndmatrix): return sum(self[...,newaxis] * other[...,newaxis,:], -2).view(ndmatrix) else: raise TypeError("Can only matrix multiply matrices by matrices or vectors") def transpose_vector(self): return self.swapaxes(-1,-2) t = property(transpose_vector) class ndcolumn(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): raise TypeError("Cannot matrix multiply columns with anything") def transpose_vector(self): return self.view(ndrow) t = property(transpose_vector) class ndrow(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other, -1).view(ndarray) if isinstance(othe, ndmatrix): raise notimplemented else: raise TypeError("Can only matrix multiply rows with matrices or columns") def transpose_vector(self): return self.view(ndcolumn) t = property(transpose_vector) def matrix(*args, **kwargs): m = array(*args, **kwargs) return m.view(ndmatrix) def column(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndcolumn) def row(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndrow) if __name__ == "__main__": import doctest, matrix doctest.testmod(matrix)
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/30/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
That should be allowed. (N,1)*(1,M) is just an (N,M) matrix with entries C[i,j] = A[i,0]*B[0,] I kind of like the idea of using call for multiply, though. If it doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot". --bb
![](https://secure.gravatar.com/avatar/5b2449484c19f8e037c5d9c71e429508.jpg?s=120&d=mm&r=g)
On 3/29/07, Bill Baxter <wbaxter@gmail.com> wrote:
I thought about that a little, and while I agree that it could be allowed, I'm not sure that it should be allowed. It's a trade off between a bit of what I would guess is little used functionality with some enhanced error checking (I would guess that usually row*column signals a mistake). However, I don't care much one way or the other; it's not hard to allow. I kind of like the idea of using call for multiply, though. If it
doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot".
We'll see how it goes down this time. I've proposed using call before, since I've thought the matrix situation was kind of silly for what seems like ages now, but it always sinks without a ripple. -- //=][=\\ tim.hochberg@ieee.org
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
It's really a sort of tensor product, so use outer(.,.). In my mind, row and column vectors are *not* matrices, they only have a single dimension. On the other hand (r)(c) is really the application of the dual vector r (a functional) to the vector c, i.e., r is a map from vectors into the reals (complex). However, I think overloading the multiplication in this case is reasonable. I kind of like the idea of using call for multiply, though. If it
doesn't turn out to have any major down sides it could be a good way to give ndarray a concise syntax for "dot".
Hmmm, have to try it a bit to see how it looks. Overall, I like this approach. Chuck
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
Looks promising!
I'm also for allowing it; from the perspective of one writing code that "looks like" typical (if potentially technically incorrect) notation, being able to write direct analogs to a'b (to get a scalar) or ab' (to get an MxN matrix) and have everything "work" would be quite useful. Especially in various reconstruction tasks (e.g. using svd to approximate a given matrix with a lower-rank one), the "outer product" of a row and a column vector comes up often enough that I would be loathe to have it raise an error.
The call syntax is nice, if a bit opaque to someone looking at the code for the first time. On the other hand, it's explicit in a way that overloaded multiplication just isn't. I like it (for what little that's worth). Zach
![](https://secure.gravatar.com/avatar/915f9e3213021d75d43bb4262167b067.jpg?s=120&d=mm&r=g)
On 3/30/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
It's useful for many things. You can use it in the computation of the least squares best fits between sets of points. The sum of p_i p_^t comes up in that context. (Which is also the covariance matrix of the points) The derivative of a unit vector (useful for mass spring dynamics calcs, among other things, I'm sure) is given by something like c * (I - x x^t). Householder reflections are useful for a lot of things such as implementing QR factorization. They're given by : H = I - 2 x x^t / x^t x http://www.cs.ut.ee/~toomas_l/linalg/lin2/node6.html I'm positive I've seen col * row come up in other contexts too, though I can't think of any other particulars right now. --bb
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 3/29/07, Timothy Hochberg <tim.hochberg@ieee.org> wrote:
<snip>
Note, however that you can't (for instance) multiply column vector with a row vector:
Well, let's carry this to extremes. Strictly speaking, from the functional point of view, (r)(c) == (c)(r), the dual of the dual is the original vector space, for finite dimensional spaces anyway. However, custom reserves the latter form for the tensor product, and so that is probably a good thing to keep. However, to really make the transpose operate correctly it should also conjugate. This could just be a flag in the row vector, no changes in the data needed, but vdot would be called instead of dot when computing (c.t)(c). A similar thing could be done for matrices if vdot were extended to deal with matrices -- currently it only works for scalars and vectors. Now, products like (c.t)(M.t) would need to conjugate/transpose both, but this could be computed as ((M)(c)).t. So on and so forth. A lot of these types of operations are in BLAS, IIRC, so in some sense this is just a mapping to BLAS. Chuck
![](https://secure.gravatar.com/avatar/95198572b00e5fbcd97fb5315215bf7a.jpg?s=120&d=mm&r=g)
On 3/27/07, Alan G Isaac <aisaac@american.edu> wrote:
I'm personally not really partial to any side of this discussion, given how I don't use matrices at all (I'm content with a simple mental model of arrays, and use dot() as necessary). But it's probably worth mentioning, as a data point, that the python language has a well established precedent for 'iteration over a FOO where a FOO is yielded': In [1]: s='abc' In [2]: type(s) Out[2]: <type 'str'> In [3]: map(type,s) Out[3]: [<type 'str'>, <type 'str'>, <type 'str'>] I know strings aren't matrices, so take this as you will in terms of being plus, neutral or minus regarding your discussion. I'm just offering the data, not interpreting it :) Cheers, f
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Tue, 27 Mar 2007, Robert Kern apparently wrote:
Gram-Schmidt orthogonalization
Robert wrote:
It is really not my intent to be irritating or to move the posts. But I had two reasons to ask. First, your claim is false that this was important in your example. (To see this, replace M[1:] with M.A[1:].) So actually your code would work the same if iteration over matrices were producing arrays (after replacing M[0] with M[0,:] to match that case). Second, when I spoke of *desirability*, the output is relevant. Try nump.mat(ortho) to see what I mean. If iteration were to produce arrays, the outcome of implementing the algorithm (e.g., using ``dot``) would be I suggest more "desirable", in that numpy.mat(ortho) would work as hoped/expected. In this sense, the example perhaps speaks against your intent. You offer code that would work just fine if iteration yielded arrays. Apologies in advance if this again seems tangential to my original request. To me it does not.
If you really want iterates to be arrays, just use .A.
I have mentioned that several times. It is orthogonal to the design issue I raised. (I suggested that under the current behavior you do not usually end up with what you want when iterating over matrices.) But again, I recognize that I am alone in that view. I am just trying to understand why. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
You're right, it does, but only because the first place it is used, it gets operated with the matrix row-vector base before it does anything else, gets coerced into a matrix, and then gets assigned back to rowv. I don't tolerate that kind of implicitness in my code.
Second, when I spoke of *desirability*, the output is relevant. Try nump.mat(ortho) to see what I mean.
<shrug> If I wanted something other than a list of matrix row-vectors (which I might), then I'd post-process. Keeping the clarity of the main algorithm is also important.
Ditto for not using matrix objects at all.
People have been giving you reasons, over and over again. You are simply refusing to listen to them. You have a use case for arrays being the iterates. You are presuming that the only argument that can beat that is another use case for matrix objects being the iterates. This is not true; there are other principles at work. If we were to make matrix objects behave in different ways for different methods to satisfy what one person thinks is the most convenient way to use each method, we would have a patchwork object with no overall principle that people can use to understand what is going on. They would simply have to remember which behavior someone thought was the most convenient for each thing. Remember the lessons of rand() and randn(). -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Wed, 28 Mar 2007, Robert Kern wrote:
People have been giving you reasons, over and over again. You are simply refusing to listen to them.
Exploring whether the reasoning is adequate is not the same as refusing to listen. I do not presume my view is correct.
Put slightly differently: given the surprising passion of the attacks at the suggestion that perhaps iteration over a matrix might more consistently yield arrays, I presumed there must be *many* instances in which it was obviously desirable that such iteration should yield matrices. So I asked to see some. In the context of this discussion, I found the (lack of) responses very interesting. Even in your thoughtful response it proved irrelevant rather than important for iteration over matrices to yield matrices. I understand that some people claim that a general principle of consistency is involved. I have not been able to understand this particular design decision as a matter of consistency, and I have tried to say why. However I am just a user (and supporter) of numpy, and as indicated in other posts, I make no pretense of deep insight into the design decisions. In this case, I simply wanted to open a discussion of a design decision, not win an "argument". Anyway, I understand that I am being perceived as bull-headed here, so I'll let this go. Thanks for your attempt to help me see the virtues of the current design. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/e815306c7885ec5d66dafb2f927c5faa.jpg?s=120&d=mm&r=g)
It's because the property that A[i] == A[i,...] is much more important to most numpy users than the results of a particular (mis) use of the matrix class. This has been explained in many different contexts over many different email messages by many different people. You're not looking at the big picture, which has been fairly explicitly spelled out by myself and others. I'm not even sure why I'm replying at this point. Zach On Mar 28, 2007, at 3:25 PM, Alan Isaac wrote:
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Wed, Mar 28, 2007 at 05:25:57PM -0500, Alan Isaac wrote:
Matrices strike me as a bit of an anomaly. I would expect an N-dimensional container to contain (N-1)-dimensional objects. But matrices cannot have dimensions less than 2, (and in numpy, can't have dimensions higher than 2). So, to me these operations don't make much sense. Is a matrix then more of a collection of matrices than a container of matrices? I'm glad we have the ndarray concept in numpy where these concepts consistently make sense. Cheers Stéfan
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Wed, Mar 28, 2007 at 07:05:00PM -0500, Alan Isaac wrote:
Doesn't seem to be the way the matrix world works though: octave:2> x = zeros(3,3,5); octave:3> size(x) ans = 3 3 5 octave:4> size(x(:,:,1)) ans = 3 3 octave:5> size(x(:,1,1)) ans = 3 1 Cheers Stéfan
![](https://secure.gravatar.com/avatar/37ee74afc8985c99adac9eb37ab2b10e.jpg?s=120&d=mm&r=g)
On 3/27/07, Zachary Pincus <zpincus@stanford.edu> wrote:
I find it useful if M[i] indexes the matrix interpreted as a column-vector (vec(M), or M(:) in Matlab notation), e.g., you could pick out the diagonal as M[::n+1], or if N is an index-set of certain elements in M, then you could access those elements as M[N]. Then I would also say that iterating over a matrix should just return the elements of vec(M) one by one. In general, I think the Matlab notation works well: * M[I,J] where I and J and index-sets returns a len(I) x len(J) matrix * M[I] returns a len(I) vector (a column or row vector depending on orientation of I) The index-sets could be a single integer, a list, or a slice. But I can see there are good arguments for using the NumPy conventions as well...
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
Hi Zachary, I think your response highlights very well the apparent design flaw. Here is your response to my request for a use case where iteration over a matrix should yield matrices: do not iterate! Does some part of you not find that just a little bizarre?! As for offering as a use case that with a matrix M you are allowed to use M[0] instead of M[0,:] when you want the first row as a matrix, I really cannot take that trivial convenience seriously as the basis of a fundamental design decision. However there is no sympathy for my view on this, so I am not debating it any more. Instead I have asked a simple question: what use case is this design decision supporting? I am interested in this so that I can see into the decision better. Nowith Chris proposes that "M[i] == M[i,:] by definition". If so, that is an design-based answer to my question. I agree that M[i,:] should return a matrix. But my understanding was different: I thought M[i] relied on standard Python idexing while M[i,:] was numpy indexing. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Perhaps our differences lies in two things: 1. the fact that the text books typically take the column vector as the default. For a Python version, based on C it makes more sense to treat the rows as vectors, as data is stored contiguously by row. 2. the convention has been proposed that the vector is more conveniently implemented as a matrix, where one dimension is one. The vector could be treated as a subclass of the matrix but this adds complexity with little clear benefit. PyMatrix has matrix methods isVector, isCVector and isRVector. I can see some merit in conforming to text book usage and would be glad to consider changes when I complete the port to numpy, in a few months. Colin W.
Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
On Sun, 25 Mar 2007, "Colin J. Williams" apparently wrote:
An array and a matrix are different animals. Conversion from one to the other should be spelled out.
But you are just begging the question here. The question is: when I iterate across matrix rows, why am I iterating across matrices and not arrays. This seems quite out of whack with general Python practice. You cannot just say "conversion should be explicit" because that assumes (incorrectly actually) that the rows are matrices. The "conversion should be explicit" argument actually cuts in the opposite direction of what you appear to believe. Cheers, Alan Isaac
![](https://secure.gravatar.com/avatar/b24e93182e89a519546baa7bafe054ed.jpg?s=120&d=mm&r=g)
Alan G Isaac wrote:
Alan, Yes, this is where we appear to differ. I believe that vectors are best represented as matrices, with a shape of (1, n) or (m, 1). The choice of these determines whether we have a column or a row vectors. Thus any (m, n) matrix can be considered as either a collection of column vectors or a collection of row vectors. If the end result is required as an array or a list, this can be done explicitly with X[1].A or A[1].tolist(). Here, A is a property of the M (matrix) class.
Cheers, Alan Isaac
A long time ago, you proposed that PyMatrix should provide for matrix division in two way, as is done in MatLab. This was implemented, but PyMatrix has not yet been ported to numpy - perhaps this summer. Regards, Colin W.
participants (17)
-
Alan G Isaac
-
Alan Isaac
-
Bill Baxter
-
Bill Spotz
-
Charles R Harris
-
Christopher Barker
-
Colin J. Williams
-
Fernando Perez
-
Joachim Dahl
-
Paulo Jose da Silva e Silva
-
Robert Kern
-
Sebastian Haase
-
Stefan van der Walt
-
Sven Schreiber
-
Timothy Hochberg
-
Travis Oliphant
-
Zachary Pincus