Em Dom, 2007-03-25 às 13:07 -0400, Alan G Isaac escreveu:

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.

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)

On 3/26/07, Alan G Isaac wrote:

...

Em Dom, 2007-03-25 às 13:07 -0400, Alan G Isaac escreveu:

...

>>> x[1] matrix([[1, 0]]) feels wrong. (Similarly when iterating across rows.)

On Sun, 25 Mar 2007, Paulo Jose da Silva e Silva apparently wrote:

...

I think the point here is that if you are using matrices, then all you "should" want are matrices, just like in MATLAB: >> b = A(1, :) b = 1 2

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]])

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

Alan G Isaac wrote:

...

Em Dom, 2007-03-25 Ã s 13:07 -0400, Alan G Isaac escreveu:

...

>>> x[1] matrix([[1, 0]]) feels wrong. (Similarly when iterating across rows.)

On Sun, 25 Mar 2007, Paulo Jose da Silva e Silva apparently wrote:

...

I think the point here is that if you are using matrices, then all you "should" want are matrices, just like in MATLAB: >> b = A(1, :) b = 1 2

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

Perhaps this would be clearer with: >>> for rowVector in x: print item.__repr__() ... matrix([[1, 2]]) matrix([[3, 4]]) Colin W.

On 3/26/07, Colin J. Williams wrote:

Bill Baxter 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 "pure matrix" seems to me an area of exploration, does it have any application in numerical computation at this time?

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.

PyMatrix deals with lists in building a larger matrix from sub-matrices.

Suppose that we have matrices A (3, 4), B (3, 6), C (4, 2) and D (4, 8).

Then E= M([[A, B], [C, D]]) gives E (7, 10).

Numpy generally tries to treat all lists and tuples as array literals. That's not likely to change. --bb

On 3/26/07, Alan G Isaac wrote:

...

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]])

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

...

On 3/26/07, Alan G Isaac wrote:

...

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]])

On Mon, 26 Mar 2007, Bill Baxter apparently wrote:

So you're saying this is what you'd find more pythonic?

...

...

...

X[1] matrix([2,3]) X[:,1] matrix([[3, 4]]) Just trying to make it clear what you're proposing.

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.

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.

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

Bill Baxter wrote:

On 3/26/07, Colin J. Williams wrote:

...

Bill Baxter 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 "pure matrix" seems to me an area of exploration, does it have any application in numerical computation at this time?

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. It was a term I had not come across before but I assumed that you were referring to something like this link - beyond my comprehension.

http://72.14.203.104/search?q=cache:Yu9gbUQEfWkJ:math.ca/Events/winter05/abs/pdf/ma-df.pdf+pure+matrix&hl=en&ct=clnk&cd=4&gl=ca&lr=lang_en 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.

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.

...

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.

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.

...

PyMatrix deals with lists in building a larger matrix from sub-matrices.

Suppose that we have matrices A (3, 4), B (3, 6), C (4, 2) and D (4, 8).

Then E= M([[A, B], [C, D]]) gives E (7, 10).

Numpy generally tries to treat all lists and tuples as array literals. That's not likely to change.

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.

Alan G Isaac schrieb:

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.

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)

On Mon, 26 Mar 2007, "Colin J. Williams" apparently wrote:

Perhaps things would be clearer if we thought of the constituent groups of data in a matrix as being themselves matrices.

This "thinking of" is what you have suggested before. You need to explain why it is not begging the question. Cheers, Alan Isaac

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

On 3/26/07, Alan G Isaac wrote:

...

Alan G Isaac schrieb:

...

What feels wrong: iterating over a container does not give access to the contained objects. This is not Pythonic.

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!!

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

Alan G Isaac wrote:

On Mon, 26 Mar 2007, "Colin J. Williams" apparently wrote:

...

Perhaps things would be clearer if we thought of the constituent groups of data in a matrix as being themselves matrices.

This "thinking of" is what you have suggested before. You need to explain why it is not begging the question.

Cheers, Alan Isaac

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.

On 3/26/07, Alan G Isaac wrote:

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.

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

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

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

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.

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

On 3/26/07, Alan G Isaac wrote:

...

...

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.

...

On 3/26/07, Alan G Isaac wrote:

...

This is again simple assertion. **Why** would "one" expect this? Some people clearly do not.

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.

Matrices don't contain arrays. Indeed, in linear algebra they are elements

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

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

On 3/27/07, Alan Isaac wrote:

...

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.

On Tue, 27 Mar 2007, Bill Baxter wrote:

AllMyPoints = mat(rand(100,2)) # 100 two-d points for pt in AllMyPoints: xformedPt = pt * TransformMatrix # do something to transformed point

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

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.

Zachary Pincus wrote:

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.

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

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

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- 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

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 wrote:

Zachary Pincus wrote:

...

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.

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

...

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

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- 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 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

Alan Isaac wrote:

On Mon, 26 Mar 2007, Travis Oliphant 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.

Again, I only raised a *question* about whether there might be a design problem here.

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

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

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

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

Alan G Isaac wrote:

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.

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

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.

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.

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.

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

Cheers, Alan Isaac

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

On 3/27/07, Zachary Pincus wrote:

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.

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...

On Tue, 27 Mar 2007, Fernando Perez wrote:

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'

Yes. And I recall being bitten by this long ago. Now it feels natural enough, I confess. But then, on reflection, it is an oddity. ;-) Cheers, Alan Isaac

Alan G Isaac wrote:

...

...

On Tue, 27 Mar 2007, Robert Kern apparently wrote:

...

Gram-Schmidt orthogonalization

...

Alan G Isaac wrote:

...

I take it from context that you consider it desirable to end up with a list of matrices?

Robert 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.

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).

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.

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.

Ditto for not using matrix objects at all.

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.

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

On Wed, Mar 28, 2007 at 05:25:57PM -0500, Alan Isaac wrote:

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.

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

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.

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.

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

On Wed, 28 Mar 2007, Stefan van der Walt wrote:

Matrices strike me as a bit of an anomaly. I would expect an N-dimensional container to contain (N-1)-dimensional objects.

Yes indeed. Cheers, Alan Isaac

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:

On Tue, 27 Mar 2007, Zachary Pincus apparently wrote:

...

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.

So would you agree that in this case M[0] and M[0,:] would have to be different things? Or would you then want even such numpy indexing not to produce matrices??

Cheers, Alan Isaac

** 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 **

On Thu, 29 Mar 2007, Bill Spotz apparently wrote:

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,:].

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]?

"""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.

...

...

m[0] ndmatrix([[1, 0], [0, 1]])

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.

...

...

m[:,0] ndrow([[1, 0], [2, 0], [3, 0]])

Finally if we index on the last index, we get an array of column matrices.

...

...

m[...,0] ndcolumn([[1, 0], [2, 0], [3, 0]])

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':

...

...

n = matrix([[0, 1], [3, 0]]) (m)(n) ndmatrix([[[0, 1], [3, 0]], <BLANKLINE> [[0, 2], [6, 0]], <BLANKLINE> [[0, 3], [9, 0]]])

Things behave similarly when multiplying a matrix with a row vector or a row vector with a column vector.

...

...

c = column([3, 5]) (m)(c) ndcolumn([[ 3, 5], [ 6, 10], [ 9, 15]]) (n)(c) ndcolumn([5, 9])

...

...

r = row([[2, 1], [3, 4]]) (r)(c) array([11, 29]) (r[0])(c) 11

Note, however that you can't (for instance) multiply column vector with a row vector:

...

...

(c)(r) Traceback (most recent call last): ... TypeError: Cannot matrix multiply columns with anything

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).

...

...

r ndrow([[2, 1], [3, 4]]) r.t ndcolumn([[2, 1], [3, 4]]) n.t ndmatrix([[0, 3], [1, 0]])

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)

On 3/29/07, Bill Baxter wrote:

On 3/30/07, Timothy Hochberg wrote:

...

Note, however that you can't (for instance) multiply column vector with a row vector:

...

...

...

(c)(r) Traceback (most recent call last): ... TypeError: Cannot matrix multiply columns with anything

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 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

On 3/29/07, Timothy Hochberg wrote:

On 3/29/07, Bill Baxter wrote:

...

On 3/30/07, Timothy Hochberg wrote:

...

Note, however that you can't (for instance) multiply column vector with a row vector:

...

...

...

(c)(r) Traceback (most recent call last): ... TypeError: Cannot matrix multiply columns with anything

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 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.

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

On 3/30/07, Timothy Hochberg wrote:

On 3/29/07, Bill Baxter wrote:

...

On 3/30/07, Timothy Hochberg wrote:

...

Note, however that you can't (for instance) multiply column vector with a row vector:

...

...

...

(c)(r) Traceback (most recent call last): ... TypeError: Cannot matrix multiply columns with anything

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 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.

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

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.

Add to that quasi-Newton updates for quadratic models in optimization, like the BFGS. Paulo