Hi Jarrod,
any news with the 1.0.5? If you have same prerelease, I'd like to test it. Debian has just moved from python2.4 to python2.5 yesterday, so I'd like to test numpy in advance, I am sure there will be some issues to fix.
Ondrej
On Mon, Apr 14, 2008 at 12:06 PM, Ondrej Certik ondrej@certik.cz wrote:
Hi Jarrod,
any news with the 1.0.5? If you have same prerelease, I'd like to test it. Debian has just moved from python2.4 to python2.5 yesterday, so I'd like to test numpy in advance, I am sure there will be some issues to fix.
What is the plan with the release? There are some minor problems in the Debian package, some of which are fixed by the new release. I didn't fix those in Debian as I thought the new release is coming out. But if it's going to take let's say month or two, I'll fix the current Debian package now.
Thanks, Ondrej
2008/4/23 Ondrej Certik ondrej@certik.cz:
What is the plan with the release? There are some minor problems in the Debian package, some of which are fixed by the new release. I didn't fix those in Debian as I thought the new release is coming out. But if it's going to take let's say month or two, I'll fix the current Debian package now.
The question is: what blocks the release of 1.1? The following tickets deserve attention:
http://projects.scipy.org/scipy/numpy/ticket/750 http://projects.scipy.org/scipy/numpy/ticket/605 http://projects.scipy.org/scipy/numpy/ticket/551 http://projects.scipy.org/scipy/numpy/ticket/738 http://projects.scipy.org/scipy/numpy/ticket/743 http://projects.scipy.org/scipy/numpy/ticket/748 http://projects.scipy.org/scipy/numpy/ticket/751 http://projects.scipy.org/scipy/numpy/ticket/736
Further issues:
- Alan's matrix indexing: x[0][0] for matrices currently yield x[0]
Regards Stéfan
On Wed, Apr 23, 2008 at 6:21 AM, Stéfan van der Walt stefan@sun.ac.za wrote:
The question is: what blocks the release of 1.1? The following tickets deserve attention:
http://projects.scipy.org/scipy/numpy/ticket/750 http://projects.scipy.org/scipy/numpy/ticket/605 http://projects.scipy.org/scipy/numpy/ticket/551 http://projects.scipy.org/scipy/numpy/ticket/738 http://projects.scipy.org/scipy/numpy/ticket/743 http://projects.scipy.org/scipy/numpy/ticket/748 http://projects.scipy.org/scipy/numpy/ticket/751 http://projects.scipy.org/scipy/numpy/ticket/736
Further issues:
- Alan's matrix indexing: x[0][0] for matrices currently yield x[0]
I was still hoping to get 1.1.0 out ASAP. Unfortunately, I have been too busy for the last week to pay attention to the release blockers. The last time I looked the only issues I felt absolutely needed to be resolved before the release were: 1. Testing the Windows binary, which is now done. 2. Testing the Mac binary, which is now done. 3. Testing or removing fromregex, which is now done: http://projects.scipy.org/scipy/numpy/ticket/719
I haven't looked carefully at the above tickets, but at a cursory glance it didn't seem like there were any new regressions. Stefan could you confirm that? If so, I would like to branch 1.1.x and tag 1.1.0 on Friday night. Unless there is some concrete reason that we need to delay this release any longer, I will send out an email later today with an 'official' timeline for the release, which will probably be on Monday if I tag on Friday night.
Thanks,
Hi Jarrod
Of those tickets, the following are serious:
http://projects.scipy.org/scipy/numpy/ticket/605 (a patch is available?, David Huard) Fixing of histogram.
http://projects.scipy.org/scipy/numpy/ticket/551 (old regression, Chuck and Travis) Unpickled arrays don't work as expected, because of memory alignment issues.
http://projects.scipy.org/scipy/numpy/ticket/748 (old regression, Stefan) ifft pads incorrectly. I can fix this tonight, if someone would verify that this is, indeed, a bug.
http://projects.scipy.org/scipy/numpy/ticket/751 (old/new regression, Stefan, David C*) ALL section not recognised in site.cfg. In Python 2.6, the DEFAULT section no longer works. I replaced it with ALL, but apparently not everything works as expected. David has a way to reproduce the problem, so I hope he can have a look, otherwise I'll try and fix it tomorrow.
Regards Stéfan
2008/4/23 Jarrod Millman millman@berkeley.edu:
On Wed, Apr 23, 2008 at 6:21 AM, Stéfan van der Walt stefan@sun.ac.za wrote:
The question is: what blocks the release of 1.1? The following tickets deserve attention:
http://projects.scipy.org/scipy/numpy/ticket/750 http://projects.scipy.org/scipy/numpy/ticket/605 http://projects.scipy.org/scipy/numpy/ticket/551 http://projects.scipy.org/scipy/numpy/ticket/738 http://projects.scipy.org/scipy/numpy/ticket/743 http://projects.scipy.org/scipy/numpy/ticket/748 http://projects.scipy.org/scipy/numpy/ticket/751 http://projects.scipy.org/scipy/numpy/ticket/736
Further issues:
- Alan's matrix indexing: x[0][0] for matrices currently yield x[0]
I was still hoping to get 1.1.0 out ASAP. Unfortunately, I have been too busy for the last week to pay attention to the release blockers. The last time I looked the only issues I felt absolutely needed to be resolved before the release were:
- Testing the Windows binary, which is now done.
- Testing the Mac binary, which is now done.
- Testing or removing fromregex, which is now done:
http://projects.scipy.org/scipy/numpy/ticket/719
I haven't looked carefully at the above tickets, but at a cursory glance it didn't seem like there were any new regressions. Stefan could you confirm that? If so, I would like to branch 1.1.x and tag 1.1.0 on Friday night. Unless there is some concrete reason that we need to delay this release any longer, I will send out an email later today with an 'official' timeline for the release, which will probably be on Monday if I tag on Friday night.
Thanks,
-- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
2008/4/23, Stéfan van der Walt stefan@sun.ac.za:
Hi Jarrod
Of those tickets, the following are serious:
http://projects.scipy.org/scipy/numpy/ticket/605 (a patch is available?, David Huard) Fixing of histogram.
I haven't found a way to fix histogram reliably without breaking the current behavior. There is a patch attached to the ticket, if the decision is to break histogram.
David
Wed, 23 Apr 2008 16:20:41 -0400, David Huard wrote:
2008/4/23, Stéfan van der Walt stefan@sun.ac.za:
Of those tickets, the following are serious:
http://projects.scipy.org/scipy/numpy/ticket/605 (a patch is available?, David Huard) Fixing of histogram.
I haven't found a way to fix histogram reliably without breaking the current behavior. There is a patch attached to the ticket, if the decision is to break histogram.
Summary of the facts (again...):
a) histogram's docstring does not match its behavior wrt discarding data
b) given variable-width bins, histogram(..., normed=True) the results are wrong
c) it might make more sense to handle discarding data in some other way than what histogram does now
I think there are now a couple of choices what to do with this:
A) Change the semantics of histogram function. Old code using histogram will just simply break, maybe in mysterious ways.
B) Rename the bins parameter to bin_edges or something else, so that any old code using histogram immediately raises an exception that is easily understood.
C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.
D) Or, retain the old behavior (values below lowest bin ignored) and just fix the docstring and the normed=True bug? (I have a patch doing this.)
So which one (or something else) do we choose for 1.1.0?
I think a long term strategy needs to be adopted for histogram. Right now there is a great confusion in what the "bins" keyword does. Right now it is defined as the lower edge of each bin, meaning that the last bin is open ended and [inf,bin0> does not exist. While this may not be the right thing to fix in 1.1.0, I would really like to see it fixed somewhere down the line.
On Apr 24, 2008, at 10:28 AM, Pauli Virtanen wrote:
Wed, 23 Apr 2008 16:20:41 -0400, David Huard wrote:
I haven't found a way to fix histogram reliably without breaking the current behavior. There is a patch attached to the ticket, if the decision is to break histogram.
Summary of the facts (again...):
a) histogram's docstring does not match its behavior wrt discarding data
This is an easy fix and should definitively go into 1.1.0 :)
b) given variable-width bins, histogram(..., normed=True) the results are wrong
Also a quick fix that should be part of 1.1.0
c) it might make more sense to handle discarding data in some other way than what histogram does now
I would like to see this, but it does not have to happen in 1.1.0 :)
I think there are now a couple of choices what to do with this:
A) Change the semantics of histogram function. Old code using histogram will just simply break, maybe in mysterious ways
Not really a satisfactory approach. I really don't mind, even though it would break some code of mine. I would rather see a better function and have to do some code changes, than the current confusion. Other people will likely disagree.
B) Rename the bins parameter to bin_edges or something else, so that any old code using histogram immediately raises an exception that is easily understood.
Given this approach bin_edges would contain one more value than bins given that the right edge of the last bin has to be defined.
C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.
I think this is probably the best approach to accommodate everyone.
So which one (or something else) do we choose for 1.1.0?
-- Pauli Virtanen
Cheers Tommy
Pauli Virtanen wrote:
C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.
D) Or, retain the old behavior (values below lowest bin ignored) and just fix the docstring and the normed=True bug? (I have a patch doing this.)
So which one (or something else) do we choose for 1.1.0?
I like either C or D, but prefer C if it can be done before 1.1.0.
-Travis
The problem I see with C is that it will break compatibility with the other histogram functions, which also use bins.
So here is suggestion E:
The most common use case ( I think) is the following: h, b = histogram(r, number_of_bins, normed=True/False) for which the function behaves correctly.
Assuming we want the next version to : ignore values outside of range and accept and return the bin edges instead of the left edges, here could be the new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None, normed=False, new=False)
If new=False, return the histogram and the left edges, with a warning that in the next version, the edges will be returned. If new=True, return the histogram and the edges. If range is given explicitly , raise a warning saying that in the next version, the outliers will be ignored. To ignore outliers, use new=True. If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.
In 1.2, set new=True as the default, and in 2.3, remove new altogether.
David
2008/4/24 Travis E. Oliphant oliphant@enthought.com:
Pauli Virtanen wrote:
C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.
D) Or, retain the old behavior (values below lowest bin ignored) and just fix the docstring and the normed=True bug? (I have a patch doing this.)
So which one (or something else) do we choose for 1.1.0?
I like either C or D, but prefer C if it can be done before 1.1.0.
-Travis
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Thu, Apr 24, 2008 at 1:22 PM, David Huard wrote:
Assuming we want the next version to : ignore values outside of range and accept and return the bin edges instead of the left edges, here could be the new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None, normed=False, new=False)
If new=False, return the histogram and the left edges, with a warning that in the next version, the edges will be returned. If new=True, return the histogram and the edges. If range is given explicitly , raise a warning saying that in the next version, the outliers will be ignored. To ignore outliers, use new=True. If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.
In 1.2, set new=True as the default, and in 2.3, remove new altogether.
+1 That sounds fine to me assuming 2.3 is 1.3.
2008/4/24 Jarrod Millman millman@berkeley.edu:
On Thu, Apr 24, 2008 at 1:22 PM, David Huard wrote:
Assuming we want the next version to : ignore values outside of range
and
accept and return the bin edges instead of the left edges, here could be
the
new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None,
normed=False,
new=False)
If new=False, return the histogram and the left edges, with a warning
that
in the next version, the edges will be returned. If new=True, return the histogram and the edges. If range is given explicitly , raise a warning saying that in the next version, the outliers will be ignored. To ignore outliers, use
new=True.
If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.
In 1.2, set new=True as the default, and in 2.3, remove new altogether.
+1 That sounds fine to me assuming 2.3 is 1.3.
Indeed.
Done in r5085. I added a bunch of tests, but I'd appreciate if someone could double check before the release. This is not the time to introduce new bugs.
Hopefully this is the end of the histogram saga.
David
-- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/ _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
2008/4/25 David Huard david.huard@gmail.com:
2008/4/24 Jarrod Millman millman@berkeley.edu:
On Thu, Apr 24, 2008 at 1:22 PM, David Huard wrote:
Assuming we want the next version to : ignore values outside of range
and
accept and return the bin edges instead of the left edges, here could
be the
new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None,
normed=False,
new=False)
If new=False, return the histogram and the left edges, with a warning
that
in the next version, the edges will be returned. If new=True, return
the
histogram and the edges. If range is given explicitly , raise a warning saying that in the
next
version, the outliers will be ignored. To ignore outliers, use
new=True.
If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.
In 1.2, set new=True as the default, and in 2.3, remove new
altogether.
+1 That sounds fine to me assuming 2.3 is 1.3.
Indeed.
Done in r5085. I added a bunch of tests, but I'd appreciate if someone could double check before the release. This is not the time to introduce new bugs.
Hopefully this is the end of the histogram saga.
Well, it's not... there is still an issue... give me a couple of minutes to fix it.
David
-- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/ _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Fri, Apr 25, 2008 at 12:55 PM, David Huard david.huard@gmail.com wrote:
Done in r5085. I added a bunch of tests, but I'd appreciate if someone
could double check before the release. This is not the time to introduce new bugs.
Hopefully this is the end of the histogram saga.
Well, it's not... there is still an issue... give me a couple of minutes to fix it.
Hey David,
Excellent! I have one more request: Could you either send me a little blurb for the release notes or go ahead and edit it yourself: http://projects.scipy.org/scipy/numpy/milestone/1.1.0
Thanks,
On Fri, Apr 25, 2008 at 12:55 PM, Jarrod Millman millman@berkeley.edu wrote:
On Fri, Apr 25, 2008 at 12:55 PM, David Huard david.huard@gmail.com wrote:
Done in r5085. I added a bunch of tests, but I'd appreciate if someone
could double check before the release. This is not the time to introduce
new
bugs.
Hopefully this is the end of the histogram saga.
This one?
ERROR: Ticket #632 ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.5/site-packages/numpy/core/tests/test_regression.py", line 812, in check_hist_bins_as_list hist,edges = np.histogram([1,2,3,4],[1,2]) File "/usr/lib/python2.5/site-packages/numpy/lib/function_base.py", line 184, in histogram raise ValueError, 'Use new=True to pass bin edges explicitly.' ValueError: Use new=True to pass bin edges explicitly.
Chuck
Thanks Chuck,
I didn't know there were other tests for histogram outside of test_function_base.
The error is now raised only if bins are passed explicitly and normed=True.
David
2008/4/25 Charles R Harris charlesr.harris@gmail.com:
On Fri, Apr 25, 2008 at 12:55 PM, Jarrod Millman millman@berkeley.edu wrote:
On Fri, Apr 25, 2008 at 12:55 PM, David Huard david.huard@gmail.com wrote:
Done in r5085. I added a bunch of tests, but I'd appreciate if
someone
could double check before the release. This is not the time to
introduce new
bugs.
Hopefully this is the end of the histogram saga.
This one?
ERROR: Ticket #632
Traceback (most recent call last): File "/usr/lib/python2.5/site-packages/numpy/core/tests/test_regression.py", line 812, in check_hist_bins_as_list hist,edges = np.histogram([1,2,3,4],[1,2]) File "/usr/lib/python2.5/site-packages/numpy/lib/function_base.py", line 184, in histogram raise ValueError, 'Use new=True to pass bin edges explicitly.' ValueError: Use new=True to pass bin edges explicitly.
Chuck
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Wed, Apr 23, 2008 at 2:32 PM, Ondrej Certik ondrej@certik.cz wrote:
On Mon, Apr 14, 2008 at 12:06 PM, Ondrej Certik ondrej@certik.cz wrote:
Hi Jarrod,
any news with the 1.0.5? If you have same prerelease, I'd like to test it. Debian has just moved from python2.4 to python2.5 yesterday, so I'd like to test numpy in advance, I am sure there will be some issues to fix.
What is the plan with the release? There are some minor problems in the Debian package, some of which are fixed by the new release. I didn't fix those in Debian as I thought the new release is coming out. But if it's going to take let's say month or two, I'll fix the current Debian package now.
The problem is, that it was noticed that non-backwards compatible changes were committed into svn. Most noticeably a complete rewrite of the "ma" (Masked Array) package was done. As a consequence it was decided, they no one would be interested in "temporarily taking those changes out" "just to release 1.0.5".
Instead the next release will be called 1.1 -- which is (only; rather still) "mostly" compatible with 1.0.x What used to be referred to a the 1.1 version, that can break more stuff, to allow for a cleaner design, will now be 1.2
HTH, Sebastian Haase
On Wed, Apr 23, 2008 at 4:56 PM, Sebastian Haase haase@msg.ucsf.edu wrote:
On Wed, Apr 23, 2008 at 2:32 PM, Ondrej Certik ondrej@certik.cz wrote:
On Mon, Apr 14, 2008 at 12:06 PM, Ondrej Certik ondrej@certik.cz wrote:
Hi Jarrod,
any news with the 1.0.5? If you have same prerelease, I'd like to test it. Debian has just moved from python2.4 to python2.5 yesterday, so I'd like to test numpy in advance, I am sure there will be some issues to fix.
What is the plan with the release? There are some minor problems in the Debian package, some of which are fixed by the new release. I didn't fix those in Debian as I thought the new release is coming out. But if it's going to take let's say month or two, I'll fix the current Debian package now.
The problem is, that it was noticed that non-backwards compatible changes were committed into svn. Most noticeably a complete rewrite of the "ma" (Masked Array) package was done. As a consequence it was decided, they no one would be interested in "temporarily taking those changes out" "just to release 1.0.5".
Instead the next release will be called 1.1 -- which is (only; rather still) "mostly" compatible with 1.0.x What used to be referred to a the 1.1 version, that can break more stuff, to allow for a cleaner design, will now be 1.2
OK, so I think we are going to try to fix the package in Debian right now and not wait for the release.
I think this uncertainty about releases and about backward compatibility is not good. The same about uncertainty where f2py is going to end up. But you know the drill "talk is cheap, show me the code", so I know I could step up and help Jarrod with the release, but unfortunately, I don't have time for it. :)
Ondrej
On 23/04/2008, Alan Isaac aisaac@american.edu wrote:
On Wed, 23 Apr 2008, Sebastian Haase wrote:
What used to be referred to a the 1.1 version, that can break more stuff, to allow for a cleaner design, will now be 1.2
So ... fixing x[0][0] for matrices should wait until 1.2. Is that correct?
It seems to me that everyone agrees that the current situation is broken, but there is considerable disagreement on what the correct fix is. My inclination is to apply the minimal fix you suggest, with tests and documentation, as a stopgap, and let the different factions sort it out by discussion on the way to 1.2.
Objections?
Anne
On Wed, Apr 23, 2008 at 11:32 AM, Anne Archibald peridot.faceted@gmail.com wrote:
So ... fixing x[0][0] for matrices should wait until 1.2. Is that correct?
It seems to me that everyone agrees that the current situation is broken, but there is considerable disagreement on what the correct fix is. My inclination is to apply the minimal fix you suggest, with tests and documentation, as a stopgap, and let the different factions sort it out by discussion on the way to 1.2.
+1 That sound reasonable to me.
2008/4/23 Jarrod Millman millman@berkeley.edu:
On Wed, Apr 23, 2008 at 11:32 AM, Anne Archibald peridot.faceted@gmail.com wrote:
So ... fixing x[0][0] for matrices should wait until 1.2. Is that correct?
It seems to me that everyone agrees that the current situation is broken, but there is considerable disagreement on what the correct fix is. My inclination is to apply the minimal fix you suggest, with tests and documentation, as a stopgap, and let the different factions sort it out by discussion on the way to 1.2.
+1 That sound reasonable to me.
Done in r5072.
Cheers Stéfan
On Wed, 23 Apr 2008, Stéfan van der Walt wrote:
Done in r5072.
Much appreciated. I have updated URL:http://www.scipy.org/MatrixIndexing to reflect this change (and its provisional status). Alan
Alan Isaac wrote:
I have updated URL:http://www.scipy.org/MatrixIndexing to reflect this change (and its provisional status).
Thanks for writing this up -- it really clarifies what's being proposed. A few comments on that write up:
For matrix x, should x[0].A[0] == x.A[0][0]? That is, should the array attribute of a Row Vector or Column Vector be 1d?
Absolutely -- a vector is a 1-d array. The only difference here is that we can preserve the concept of a row vs. a column vector -- key to linear algebra, which is the whole point of the Matrix class. I'd be just as happy if a RowVector an ColumnVector were a more unique object -- a 1-d array with the operation overloaded, rather than a (1,n) or (n,1) matrix, but I understand that would be more work.
Deviations from the behavior of ndarrays. Specifically, iteration over a matrix will not yield 1d NumPy arrays.
Why should they? indexing a nested list gives a different result than indexing an ndarray, etc, etc, etc. The POINT of matrixes is that they are different! This proposal makes them a little more similar to ndarrays than the old behavior (always returning a matrix).
Inconsistency in the indexing style for producing row vectors and column vectors
not really:
row vector: M[i,:] column vector: M[:,i]
The fact that there is another way to get a row vector: M[i] is a bonus. I once proposed making M[i] raise an exception for the sake of enforcing consistency, but no one liked that!
Loss of a standard way to extract a row or a column as a submatrix (rather than a vector)
What's wrong with:
M[i:i+1, :] and M[:, i:i+1]
This is totally consistent with arrays (and other python sequences). Yes, it's a bit more awkward, but with the new vectors, it's also going to be needed far less.
No clear gain in functionality over the simpler proposal 2.
Well, I won't judge what's "clear" but I think this adds functionality. The concept of row and column vectors is a key part of linear algebra. Sure, you can model them as matrixes that happen to have only one row or one column, but that's what we had, but there seem to be real issues with use of that. Option (2), is "Let a Matrix be a container of 1d arrays.". Sounds simple, but do those arrays represent row or column vectors? Either would make sense. The proposal is for them to be row vectors, which is fine, but then if you want them to act like a row vectors, do you need to convert them back into matrices? Even if not, how do you get a column vector -- by indexing the current way, and getting a (n,1) matrix, so that now we have inconsistency between how row and column vectors are treated -- not a clean API.
And this is about a clean API for linear algebra -- otherwise, we'd just all use arrays (which many do!)
Aside from the fact that someone needs to write the code -- why don't people like the row/column vector idea? It just feels so natural to me:
A matrix is a 2-d array with some functionality overloaded to make it convenient to do linear algebra.
A vector is a 1-d array with some functionality overloaded to make it convenient to do linear algebra.
A scalar is the same in both contexts, so no need to change anything there.
And yes, maybe this will lead to tensors some day!
Another couple questions: -- would there be constructors for vectors?
np.RowVector((1,2,3,4,5)) np.ColumnVector((1,2,3,4,5,6))
and would: A_RowVector.T == A_ColumnVector
I would think so.
One more note:
All the example code I've seen during this discussion have been examples of iterating and indexing -- not one has actually been linear algebra -- perhaps we should see some of those before making any decisions.
-Chris
On Apr 23, 2008, at 3:55 PM, Christopher Barker wrote:
For matrix x, should x[0].A[0] == x.A[0][0]? That is, should the array attribute of a Row Vector or Column Vector be 1d?
Absolutely -- a vector is a 1-d array. The only difference here is that we can preserve the concept of a row vs. a column vector -- key to linear algebra, which is the whole point of the Matrix class. I'd be just as happy if a RowVector an ColumnVector were a more unique object -- a 1-d array with the operation overloaded, rather than a (1,n) or (n,1) matrix, but I understand that would be more work.
The way I envisioned this was that RowVector and ColumnVector would inherit from Matrix, so that you would inherit all of the Matrix functionality. They would just be special cases where ncol=0 or nrow=0, allowing single indexing.
If Matrix-Matrix multiplication is implemented properly, then the Row/ ColumnVector multiplication operators should automatically work.
Inconsistency in the indexing style for producing row vectors and column vectors
not really:
I agree with Alan here. M[i] returns a RowVector, but there is no corresponding single index way to retrieve a ColumnVector (unless we want to use __call__() to make M(j) return a ColumnVector . . . but this is highly unintuitive.)
Thinking out loud: we could support
M[i=#] return a RowVector M[j=#] return a ColumnVector M[#] equivalent to M[i=#]
The fact that there is another way to get a row vector: M[i] is a bonus.
Except that the behavior of M[i] is one of the driving issues of the conversation.
Loss of a standard way to extract a row or a column as a submatrix (rather than a vector)
If Row/ColumnVector inherit from Matrix, then this is not strictly true. There is no reason that these classes could not support [i,j] indexing, which means they would BE Matrices (inheritence) and they would BEHAVE like Matrices (except for V[i][j] . . . argh).
No clear gain in functionality over the simpler proposal 2.
Gains: (1) non-scalar extractions from linear algebra objects ARE and BEHAVE like linear algebra objects; (2) a clear path for dense and sparse matrices to have the same (or at least analogous) interfaces.
Aside from the fact that someone needs to write the code -- why don't people like the row/column vector idea? It just feels so natural to me:
Unless I'm misremembering, Alan is the only one who has expressed concerns and he is willing to concede to the design if others agree.
A matrix is a 2-d array with some functionality overloaded to make it convenient to do linear algebra.
A vector is a 1-d array with some functionality overloaded to make it convenient to do linear algebra.
Again, I would argue for Vectors inheriting from Matrix. I would make the case based on code reuse and elegance, although these might be overridden by some other concern.
And yes, maybe this will lead to tensors some day!
Don't see why not.
Another couple questions: -- would there be constructors for vectors?
np.RowVector((1,2,3,4,5)) np.ColumnVector((1,2,3,4,5,6))
Yes. They should be full-fledged classes, although with inheriting from Matrix, the implementation should be relatively small. What about
np.Matrix(<container of RowVector>) np.Matrix(<container of ColumnVector>)
?
and would: A_RowVector.T == A_ColumnVector
Yes.
All the example code I've seen during this discussion have been examples of iterating and indexing -- not one has actually been linear algebra -- perhaps we should see some of those before making any decisions.
Right. I have generally thought about this in the context of, say, a Krylov-space iterative method, and what that type of interface would lead to the most readable code.
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-0154 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
[CHOP]
The proposals thus far don't address two of the major issues I have with the matrix class:
1. The matrices and arrays should become more alike if possible and should share more of the same code base. From what I've seen, the people who write the code (for numpy) don't actually use matrices, they use the arrays. This is one of the main reasons matrices still have so many problems IMO; if one of the main developers were using them, they'd never have put up with it. This may be changing to some extent now, but the more we can make matrices and arrays share code and interface, the happier we'll be. None of the new rules strike me as helping in this arena and may in fact hurt. Again, IMO. 2. This doesn't support the concept of arrays of matrices/vectors, which is one thing I've always wanted out of the matrix class. For example it would be nice to be able to spell: 'A' is a 1D array of matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall, matrix multiply them together, yielding a 1D array of matrices (3D) overall. I have frequently run into various permutations of this problem. You can fake it with loops over dot, but the efficiency is very bad for small arrays, plus it's ugly.
I have lot's of vague ideas on this subject that have been floating around my head for the past couple of years. Basically, I think the way to go is to probably add some meta information to arrays that make them look a little bit more like tensors. Then, if one is careful (and things work out as I hope) matrices could be implemented as a thin layer on top of arrays. Or perhaps, if we are very lucky, done away with altogether.
Anyway, that's all very vague and hand-wavey; if I can spring free some time, I'll try to write something up in the not too distant future. I'm glad to see that any major changes are being put off for a while: I think we should be able to do better than just patch up the old matrix class.
-tim
+1
Nadav
-----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם Timothy Hochberg נשלח: ה 24-אפריל-08 04:16 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] numpy release
[CHOP]
The proposals thus far don't address two of the major issues I have with the matrix class:
1. The matrices and arrays should become more alike if possible and should share more of the same code base. From what I've seen, the people who write the code (for numpy) don't actually use matrices, they use the arrays. This is one of the main reasons matrices still have so many problems IMO; if one of the main developers were using them, they'd never have put up with it. This may be changing to some extent now, but the more we can make matrices and arrays share code and interface, the happier we'll be. None of the new rules strike me as helping in this arena and may in fact hurt. Again, IMO. 2. This doesn't support the concept of arrays of matrices/vectors, which is one thing I've always wanted out of the matrix class. For example it would be nice to be able to spell: 'A' is a 1D array of matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall, matrix multiply them together, yielding a 1D array of matrices (3D) overall. I have frequently run into various permutations of this problem. You can fake it with loops over dot, but the efficiency is very bad for small arrays, plus it's ugly.
I have lot's of vague ideas on this subject that have been floating around my head for the past couple of years. Basically, I think the way to go is to probably add some meta information to arrays that make them look a little bit more like tensors. Then, if one is careful (and things work out as I hope) matrices could be implemented as a thin layer on top of arrays. Or perhaps, if we are very lucky, done away with altogether.
Anyway, that's all very vague and hand-wavey; if I can spring free some time, I'll try to write something up in the not too distant future. I'm glad to see that any major changes are being put off for a while: I think we should be able to do better than just patch up the old matrix class.
-tim
2008/4/24 Nadav Horesh nadavh@visionsense.com:
+1
+1 to what? I'm glad that Tim chimed in -- his opinion is always valued -- but I didn't see any concrete suggestions on how to improve the situation. I'll gladly wait for his code, though :)
Regards Stéfan
Some, including me, stated that the matrix class suffer some deficiencies, so, after playing a bit with matrices, we always go back to the good old ndarray class. There were some arguments on this list what/if should be dome to correct it. I welcome any attempt to build a code that demonstrates an alternative.
Nadav.
-----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם St?fan van der Walt נשלח: ה 24-אפריל-08 13:13 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] numpy release
2008/4/24 Nadav Horesh nadavh@visionsense.com:
+1
+1 to what? I'm glad that Tim chimed in -- his opinion is always valued -- but I didn't see any concrete suggestions on how to improve the situation. I'll gladly wait for his code, though :)
Regards St?fan _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Hi, I would like to use the matrix functionality but I, like others, get by without it.
+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.
-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error. 0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.
Yes, I agree that talk is cheap but it avoids wasting time on code that will not be used and allows time to understand what people really want.
Bruce
Nadav Horesh wrote:
Some, including me, stated that the matrix class suffer some deficiencies, so, after playing a bit with matrices, we always go back to the good old ndarray class. There were some arguments on this list what/if should be dome to correct it. I welcome any attempt to build a code that demonstrates an alternative.
Nadav.
-----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם St?fan van der Walt נשלח: ה 24-אפריל-08 13:13 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] numpy release
2008/4/24 Nadav Horesh nadavh@visionsense.com:
+1
+1 to what? I'm glad that Tim chimed in -- his opinion is always valued -- but I didn't see any concrete suggestions on how to improve the situation. I'll gladly wait for his code, though :)
Regards St?fan _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Wed, 23 Apr 2008, Timothy Hochberg apparently wrote:
I think the way to go is to probably add some meta information
I have added this as Proposal 4 at URL:http://www.scipy.org/MatrixIndexing Forgive anything said the misrepresents your intent.
Cheers, Alan
On Thu, Apr 24, 2008 at 11:16 AM, Timothy Hochberg tim.hochberg@ieee.org wrote:
[CHOP]
The proposals thus far don't address two of the major issues I have with the matrix class:
The thing that seems missing to me is support for LAPACK's banded and packed (triangular) storage formats. I don't have an urgent need for it, but it seems like a serious front end for doing linear algebra should at least provide a convenient way to access all of what LAPACK offers.
Perhaps these other formats could/should be supported by different classes, though, the way it is done for sparse matrices. But the fact that more matrix types could exist some day is all the more reason to make sure the interface makes sense for more than just strided memory. My hunch is if the API is workable for both ndarray matrices and sparse matrices, then it should be ok for these other formats too.
--bb
On Wed, 23 Apr 2008, Bill Spotz apparently wrote:
Gains: (1) non-scalar extractions from linear algebra objects ARE and BEHAVE like linear algebra objects;
I believe this is not a gain, since there is a standard way to do this now, which would remain under the second proposal.
(2) a clear path for dense and sparse matrices to have the same (or at least analogous) interfaces.
This may prove true. I'm agnostic. What actually turns on this? I believe it is the following: if iterating over a sparse matrix yields sparse vectors that behave like a sparse matrix, then iterating over a matrix should return a vector that behaves like a matrix. But I have not seen any detailed discussion of this design issue. E.g., iterating over a sparse matrix could yield instead a sparse object that behaves like a 1d array.
I should add that if the first approach is taken, I agree that it is natural for these "vectors" to inherit from matrix.
With respect to the matrix object, my only "objection" has been that I'd like to hear someone state clearly what functionality is gained by following the more complex first proposal instead of the second proposal. The reason: the cost of complexity should be justified by a gain in functionality. It *may* be that the link between the behavior of matrices and that of sparse matrices is where the gain manifests, but that is not yet clear to me.
Unless I'm misremembering, Alan is the only one who has expressed concerns and he is willing to concede to the design if others agree.
I'm in no position to concede or not concede anything, since I am not writing code, but my core concerns will be met by either proposal.
Thanks! Alan
Alan G Isaac wrote:
the cost of complexity should be justified by a gain in functionality.
I don't think functionality is the right word here. the Matrix class(es) is all about clean, convenient API, i.e. style, not functionality -- we have all the functionality already, indeed we have it with plain old arrays, so I think that's really beside the point.
I just like the "feel" of Row and Column Vectors.
But this brings me back to:
We need examples of doing linear algebra operations, not just iterating operations. Frankly, if you have a matrix, and you want to iterate through all the rows, and you want to index into them as 1-d arrays, why they heck don't you just:
for r in MyMatrix.A: ....
The answer (I think), is that folks have simply been confused and annoyed by the fact that they need that extra ".A", which brings us back to style.
One extra functionally: you can index a ColumnVector with a scalar:
col = M[:,i] col[i]
You can't do that now.
This brings me to:
the people who write the code (for numpy) don't actually use matrices, they use the arrays
I don't really write numpy code, but I don't really use arrays either. Which, I suppose, means I should just shut up. But maybe I would if the API was better.
So we're back to: who does want matrices? is it only newbies coming from Matlab, where they were used to them?
For once, we do have someone that is writing the code that wants improvements, so I guess that's why this discussion has been revived.
Timothy Hochberg wrote:
- The matrices and arrays should become more alike if possible
I'm not sure I agree -- the more alike they are, the less point there is to having them.
should share more of the same code base.
Don't they share almost all the same code already? My understanding is that they are arrays with a couple operators overloaded -- how much more code could they share?
Nevertheless, I'm interested to see what Tim comes up with.
- This doesn't support the concept of arrays of matrices/vectors, which is one thing I've always wanted out of the matrix class. For example it would be nice to be able to spell: 'A' is a 1D array of matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall,
You could do this now with a 1-d object array, with a bunch of matrices in it.
matrix multiply them together, yielding a 1D array of matrices (3D) overall.
But this, or course, would fail. However, if we expand this scenario to A and B being 1-d arrays of other arrays that may not all be the same size (and thus you couldn't just use a 3-d array), and be able to run various ufuncs on them, that would be cool!
Bill Spotz wrote:
The fact that there is another way to get a row vector: M[i] is a bonus.
Except that the behavior of M[i] is one of the driving issues of the conversation.
well, yes, but a bit of a red herring, I think -- it's syntactic sugar. The issue, I think, is that even if you do:
r = M[:,i]
you can't then index that with a single index -- it's a matrix, NOT a 1-d array. Again, I proposed disallowing M[i] altogether, as it is ambiguous, but that really is gratuitous consistency.
A vector is a 1-d array with some functionality overloaded to make it convenient to do linear algebra.
Again, I would argue for Vectors inheriting from Matrix. I would make the case based on code reuse and elegance, although these might be overridden by some other concern.
I think that's a bit of an implementation detail -- we need to define behavior that we want, and decide how best to implement that. i.e.:
RowVector[i] -> scalar RowVector.A -> 1-d array (same for column)
How do we want it pretty-printed?
We need to focus on column vectors here -- It's easy to say: a row vector is a 1-d array -- that's easy and clean. It's column vectors that are tricky -- indeed they are a pain with regular old arrays.
What about
np.Matrix(<container of RowVector>) np.Matrix(<container of ColumnVector>)
I'd say you get a matrix from each of those, but they would be transposes of each-other -- another good reason for the Vector classes. though maybe those should be spelled:
np.row_stack and np.column_stack
I have generally thought about this in the context of, say, a Krylov-space iterative method, and what that type of interface would lead to the most readable code.
Can you whip up a small example, starting with the current implementation?
Bruce Southey wrote:
Hi, I would like to use the matrix functionality but I, like others, get by without it.
What does it lack that keeps you from using it? That's the key question?
+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.
This is all about being able to index a "vector" with a single index -- that's it, I think. Matlab handles this by special casing matrices that have a single dimension of one. It also has an easy way to spell a literal for column vector, though I suppose:
np.Matrix((1,2,3,4,5)).T
isn't so bad.
-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error.
But HOW can you know the orientation? a 1-d array has no orientation -- which is why the current version always returns a matrix when indexing.
0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.
I don't think we can have vector without ColumnVector, though I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.
By the way, maybe a column vector could be handy for doing array broadcasting, as an cleaner way to do:
x = np.arange(10) y = np.arange(20).reshape((-1,1))
z = FunctionOfXandY(x,y)
Final note:
If no one has an linear algebra examples, then we're wasting out time with even this cheap talk.
-Chris
On Thu, 24 Apr 2008, Christopher Barker apparently wrote:
for r in MyMatrix.A:
Of course: this possibility is mentioned explicitly in several of my posts. There are a few different answers to why asking users to do this annoying:
1. It is a deviation from the behavior of 2d arrays that does not seem to "buy" anything, so it pointlessly breaks a natural expectation. (But maybe some will show that it does buy a benefit.) So it breaks the rule: stick with the behavior of 2d arrays except for *, **, and nonscalar indexing, which is a pretty understandable rule. To this extent I agree with Tim that matrices and arrays should be as alike as it is feasible to make them, while still getting the convenience of *, **, and submatrices via nonscalar indexing. 2. Whenever you do not know ahead of time whether you will pass an array or a matrix, you must bear the (small but annoying) cost of working around your ignorance. 3. In my own experience, whenever I want to iterate over a matrix, I want the 1d arrays. Others may have different experiences, but they have not shared them. Unless we find users who actually want the vectors you want to provide them with, you seem to be promoting this approach based on fairly abstract design considerations. (Which is not necessarily bad of course.) 4. If users do want to iterate over rows and columns, it is more natural (i.e., symmetrical) to provide these as attributes of the matrix.
I do not mean for any of this to be determinative. Indeed, against this are two possible arguments.
1. For matrices and sparse matrices to nicely behave alike, your proposal will be necessary. (This case has not been made yet, however.) 2. I think Stefan might argue something like, to leave the "linear algebra world", you should have to be quite explicit, as x.A is. I find it explict enough if we use a matrix as an iterator---to me that means indeed I am leaving the matrices in favor of arrays.
Cheers, Alan Isaac
On Thu, 24 Apr 2008, Christopher Barker apparently wrote:
I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.
This bothers me.
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
Or are you still having separate classes and not just changing ``__getitem__``? OK, I see that is most likely what you meant.
Cheers, Alan
2008/4/25 Alan G Isaac aisaac@american.edu:
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
Yup. That's the current state of things.
Stéfan
2008/4/25 Alan G Isaac :
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
On Fri, 25 Apr 2008, Stefan van der Walt apparently wrote:
Yup. That's the current state of things.
I do not understand. The released "state of things" is that for matrix ``x`` we have that ``x[0]`` is a **matrix**.
I'm not working with SVN NumPy, but my understanding was that as of revision r5072 ``x[0]`` always returns a 1d array. The core requirement was that ``x[0][0]`` produce the first element of the matrix.
I do not have time to look at the revision right now, but if a matrix ``x`` we have that ``x[0]`` can return a scalar, that is very undesirable.
Cheers, Alan
2008/4/25 Alan G Isaac aisaac@american.edu:
2008/4/25 Alan G Isaac :
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
On Fri, 25 Apr 2008, Stefan van der Walt apparently wrote:
Yup. That's the current state of things.
I do not understand. The released "state of things" is that for matrix ``x`` we have that ``x[0]`` is a **matrix**.
I'm not working with SVN NumPy, but my understanding was that as of revision r5072 ``x[0]`` always returns a 1d array. The core requirement was that ``x[0][0]`` produce the first element of the matrix.
I do not have time to look at the revision right now, but if a matrix ``x`` we have that ``x[0]`` can return a scalar, that is very undesirable.
In current SVN:
In [3]: x = np.matrix(np.arange(9).reshape((3,3)))
In [5]: x Out[5]: matrix([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [6]: x[0] Out[6]: matrix([[0, 1, 2]])
In [7]: x[0,:] Out[7]: matrix([[0, 1, 2]])
In [8]: x[0][0] Out[8]: 0
In [9]: x[0,:][0] Out[9]: 0
In [10]: x[:,0] Out[10]: matrix([[0], [3], [6]])
In [11]: x[:,0][0] Out[11]: 0
In [12]: x[:2,:2][0] Out[12]: matrix([[0, 1]])
Regards Stéfan
2008/4/25 Alan G Isaac aisaac@american.edu:
I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.
The agreement was:
a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)
This breaks as little functionality as possible, at the cost of one (instead of two) API changes.
We should now discuss the proposals on the table, choose a good one, and implement all the API changes necessary for 1.2 or 2. It's a pity we have to change the API again, but the current situation is not tenable.
Cheers Stéfan
On 25/04/2008, Stéfan van der Walt stefan@sun.ac.za wrote:
2008/4/25 Alan G Isaac aisaac@american.edu:
I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.
The agreement was:
a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)
This breaks as little functionality as possible, at the cost of one (instead of two) API changes.
Hold on. There has definitely been some confusion here. This is not what I thought I was suggesting, or what Alan thought he was suggesting. I do not think special-casing matrices for which one dimension happens to be one is a good idea at all, even temporarily. This is the kind of thing that drives users crazy.
My suggested stopgap fix was to make x[0] return a 1D *array*; I feel that this will result in less special-casing. In fact I wasn't aware that anyone had proposed the fix you implemented. Can we change the stopgap solution?
We should now discuss the proposals on the table, choose a good one, and implement all the API changes necessary for 1.2 or 2. It's a pity we have to change the API again, but the current situation is not tenable.
Yes, well, it really looks unlikely we will be able to agree on what the correct solution is before 1.1, so I would like to have something non-broken for that release.
Anne
2008/4/25 Anne Archibald peridot.faceted@gmail.com:
The agreement was:
Clearly that was only the agreement in my mind. I apologise if I jumped the gun.
a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)
This breaks as little functionality as possible, at the cost of one (instead of two) API changes.
Hold on. There has definitely been some confusion here. This is not what I thought I was suggesting, or what Alan thought he was suggesting. I do not think special-casing matrices for which one dimension happens to be one is a good idea at all, even temporarily. This is the kind of thing that drives users crazy.
Apparently, not having x[0][0] == x[0,0] drives them even crazier :) I was trying to satisfy Alan's request without breaking the API in more than one place. The whole idea of matrices was that you always work with 2-dimensional arrays, even when you just extract a row. Unless you have a proper hierarchical container (as illustrated in the other patch I sent), returning a 1D array breaks at least some expectation. If that's what the matrix users want, though, then we should change it.
My suggested stopgap fix was to make x[0] return a 1D *array*; I feel that this will result in less special-casing. In fact I wasn't aware that anyone had proposed the fix you implemented. Can we change the stopgap solution?
Yes, of course. NumPy is a community project, after all.
We should now discuss the proposals on the table, choose a good one, and implement all the API changes necessary for 1.2 or 2. It's a pity we have to change the API again, but the current situation is not tenable.
Yes, well, it really looks unlikely we will be able to agree on what the correct solution is before 1.1, so I would like to have something non-broken for that release.
Unless we agree on something, and soon, that won't happen. Your workaround would break x[0] == x[0,:], so we're just swapping one set of broken functionality for another.
I'm starting to see Chris Barker's point; allowing x[0] is causing more problems than it is worth. On the other hand, how would you index into a vector (as in http://en.wikipedia.org/wiki/Vector_(spatial)) without it?
Regards Stéfan
2008/4/25 Stéfan van der Walt stefan@sun.ac.za:
I'm starting to see Chris Barker's point; allowing x[0] is causing more problems than it is worth. On the other hand, how would you index into a vector (as in http://en.wikipedia.org/wiki/Vector_(spatial)) without it?
To answer my own question: you wouldn't. Vectors won't exist, everything will be 2D, always. I could go with that.
Cheers Stéfan
2008/4/25 Alan G Isaac aisaac@american.edu:
I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.
On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:
The agreement was: a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)
1. This is **not** what I understood as the agreement (and I think the current solution is bad). I certainly did not mean to support the change that as implemented, and it is not clear to me that others did either.
2. These two goals are substantially in conflict for matrices, as we are seeing.
3. The important goal, (goal a., which everyone agrees on), has NOT been accomplished by the current change: x[0][0] raises a TypeError when x is a 1 by N matrix.
This breaks as little functionality as possible, at the cost of one (instead of two) API changes.
I cannot agree with this assessment.
the current situation is not tenable.
I agree with this! Better than the SVN behavior would be to raise an error in response to scalar indexing of matrices. I strongly hope that this behavior will not show up in a release!
Please return 1d arrays in response to scalar indexing as the provisional arrangement. This is certainly better than the new behavior.
Alan Isaac
On Fri, 25 Apr 2008, Jarrod Millman apparently wrote:
I would like, so I can't tell if there is anything so uncontroversial that it would make sense to change for the 1.1.0 release.
I think it is clear from reactions that the revision r5072 to matrix behavior should NOT go into any release.
If there is something that can be unanimously agreed on within the next 24 hours or so, I would be happy to have it included in 1.1.
I see 3 possibilities (in rank order of preference):
1. return 1d arrays in response to scalar indexing of matrices 2. raise an error on scalar indexing of matrices 3. simply revert the r5072 change.
Since there is universal accord that x[0][0] == x[0,0] is desirable, and we will get that with option 1, I think it is the right provisional behavior. I believe that option 1 provides a step in the direction of the ultimate behavior, no matter what that will be.
Cheers, Alan Isaac
On Fri, 25 Apr 2008, Stéfan van der Walt wrote:
workaround would break x[0] == x[0,:]
But there is not universal agreement that x[0] == x[0,:] is desirable. In contrast, there *is* universal agreement that x[0][0]==x[0,0] is desirable. Or so I've understood the discussion.
As you know, my view is that nonscalar indexing should return submatrices, while scalar indexing should return 1d arrays. Thus I do not see x[0] != x[0,:] as breakage: it is desirable. Indeed, without it, we will have an untenable situation, under all the different proposals. So it is a goner one way or another.
Cheers, Alan Isaac
On Fri, Apr 25, 2008 at 01:40:29PM -0400, Alan G Isaac wrote:
On Fri, 25 Apr 2008, Stéfan van der Walt wrote:
workaround would break x[0] == x[0,:]
But there is not universal agreement that x[0] == x[0,:] is desirable. In contrast, there *is* universal agreement that x[0][0]==x[0,0] is desirable. Or so I've understood the discussion.
I DON'T AGREE! Sorry for yelling, but I don't see where you get your universality. You want to break backward incompatibility to accomodated wrongful indexing, and break rightful indexing.
I don't know why people are indexing matrices with A[x][y], but they shouldn't. We have multidimensional objects, and a lot of effort has been put in that. Indexing with one dimension only is poor amongst other things for performance reasons. We also see that it breaks indexing semantics and fobrid doing clever indexing. That cannot be worked around, because by indexing with one dimension you are putting less information into the indexing routine than by indexing with two dimension.
Thus I do not see x[0] != x[0,:] as breakage: it is desirable.
You may think what you want. We may desagree on what is Right, but the breaks backward compatibility, and thus is a breakage and should be given a lot of thought.
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
I don't see where you get your universality.
Lists, tuples, arrays ... Where is the exception? Only matrices are the only 2d container I can think of that are broken this way.
Cheers, Alan Isaac
On Sat, Apr 26, 2008 at 11:15:14AM -0400, Alan G Isaac wrote:
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
I don't see where you get your universality.
Lists, tuples, arrays ... Where is the exception? Only matrices are the only 2d container I can think of that are broken this way.
List and Tuples are not 2D containers.
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
I don't see where you get your universality.
On Sat, Apr 26, 2008 at 11:15:14AM -0400, Alan G Isaac wrote:
Lists, tuples, arrays ... Where is the exception? Only matrices are the only 2d container I can think of that are broken this way.
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
List and Tuples are not 2D containers.
A list of lists is naturally conceived as 2d, which is why you can initialize a 2d matrix with a list of lists.
Try a Google search on "python two dimensional list" to confirm that I am not promoting an anomalous view.
So I do not buy the way you are trying to parse language here. (Although I'm open to a fuller argument.)
Cheers, Alan
On Sat, Apr 26, 2008 at 02:32:01PM -0400, Alan G Isaac wrote:
A list of lists is naturally conceived as 2d, which is why you can initialize a 2d matrix with a list of lists.
Try a Google search on "python two dimensional list" to confirm that I am not promoting an anomalous view.
So I do not buy the way you are trying to parse language here. (Although I'm open to a fuller argument.)
The object is not 2D. You can assemble several objects into something that looks 2D, but it is not 2D.
For exemple:
[[1, 2, 3], [1, 2]]
This is not 2D in an nd array sens. It how exposes a double successiv indexing, but doubly 1D, and no 2D.
For nD objects, I think people should rather use either nD indexing, or explicite 1D indexing of an nD object. By this is mean "A[1, ...]" rather than "A[1]". If the later breaks, I, for one, won't cry.
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
We may desagree on what is Right, but the breaks backward compatibility, and thus is a breakage and should be given a lot of thought.
I agree with this. I am a matrix user, and I have given it a lot of thought. I have been making this case for a *long* time. It has come to a head because of the announced desire to severely constrain API changes moving forward.
As I am best able to understand you abstract view, the right thing to do would be for matrices to raise an error in response to a scalar index. However, I deduce, you would just leave things alone to avoid backward incompatibility.
I weight the future more heavily. We are approaching a last chance to do things better, and we should seize it.
The right questions looking forward:
- what behavior allows the most generic code? - what behavior breaks fewest expectations? - what behavior is most useful?
Cheers, Alan Isaac
PS I cannot recall: do you use matrices?
Gael Varoquaux wrote:
On Fri, Apr 25, 2008 at 01:40:29PM -0400, Alan G Isaac wrote:
In contrast, there *is* universal agreement that x[0][0]==x[0,0] is desirable. Or so I've understood the discussion.
I don't know why people are indexing matrices with A[x][y], but they shouldn't.
I think there has been a misunderstanding here. I don't think anyone is suggesting that if a coder wants an element of a matrix, that s/he should write that:
element = M[i,j]
Rather, that one might want to extract a row from a Matrix, and then index THAT object with a scalar:
row = M[i]
now do something with each element in that row: something = row[j]
This is a rational, normal thing to do, and I think the desire for this is the core of this entire discussion, coming from the fact that in the current version of matrix both of M[i] and M[i,:] yield a matrix, which is 2-d, and cannot be indexed with a scalar. This is odd, as it is pretty natural to expect a single row (or column) to behave like a 1-d object.
Alan G Isaac wrote:
I believe that this conflicts with submatrix extraction.
If we want to keep submatrix extraction (as we should),
why? with 2-d arrays, you get a subarray with: A[i:j,:] and a 1-d array with A[i,:] -- why does that need to be different for Matrices?
- If ``v`` is a "RowVector" what behavior do we get?
Suppose the idea is that this will rescue ``x[0][0]==x[0,0]`` **and** ``x[0]=x[0,:]``. But then we must give up that ``x[0,:]`` is a submatrix.
Correct, but why should it be -- we can get a submatrix with slicing, as above. Indeed, I was about to post this comment the other day, as someone was concerned that there needs to be a distinction between a ColumnVector and a Matrix that happens to have a second dimension of one.
I think that the Vector proposal satisfies that:
if you have code like:
SubMatrix = M[i:j, k]
Then you will always get a SubMatrix, even if j == i+1. If you index like so:
Vector = M[i,k] you will always get a vector (a 1-d object)
Must ``v`` deviate from submatrix behavior in an important way? Yes: ``v[0][0]`` is an IndexError.
Correct, but if you're writing the code that generated that vector, you'd know it was a 1-d object.
Since submatrix extraction is fundamental, I think it is a *very* bad idea to give it up.
I agree -- but we're not giving it up, what we are doing is making a distinction between extracting a single row or column, and extracting a submatrix (that may or may not be a single row or column) -- just like that distinction is make for regular old ndarrays.
RowVector proposal we must give up ``x[0]==x[0,:]``. But this is just what we give up with the much simpler proposal that ``v`` be a 1d array.
no, we also give up that v acts like a row or column (particularly column) vector in computation (*, **)
We still need real linear algebra computation examples....
Alan G Isaac wrote:
I weight the future more heavily. We are approaching a last chance to do things better, and we should seize it.
Yes, but it seems that while a consensus will not be reached in time for 1.1, there is one that a change will probably occur in the next version, so there is a lot to be said for waiting until a proposal is settled on, then make the whole change at once.
The right questions looking forward:
- what behavior allows the most generic code? - what behavior breaks fewest expectations? - what behavior is most useful?
Yes, but I'm going to try to put down what I think are the key, very simple, questions:
1) Do we want to be able to extract a single row or column from a Matrix, and have it be indexable like a 1-d object?
2) Do we want to be able to do that without having to explicitly call the Array attribute: M.A[i]?
3) Do we want to have those 1-d objects act like Row or Column vectors in linear algebra operations (and broadcasting)?
4) Do we want to require slice notation to get a submatrix that happens to have a single row or column?
If we want (1) and (2), then we need to make a change. If we want (3) and (4), then something like the Row/ColumnVector proposal is needed.
I really think it's that simple.
-Chris
oops, typo!
Christopher Barker wrote:
Gael Varoquaux wrote:
I don't know why people are indexing matrices with A[x][y], but they shouldn't.
I think there has been a misunderstanding here. I don't think anyone is suggesting that if a coder wants an element of a matrix, that s/he should write that:
element = M[i,j]
I meant -- no one is suggesting:
element = M[i][j]
because, of course, M[i,j] is exactly what you should do do get an element.
Which, by the way, brings up another inconsistency in the current behavior:
M
matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
M[1,1]
5
M[1,:]
matrix([[4, 5, 6]])
so if you pull out a single element, you get a scalar(0-d), but if you pull out a row, you get a matrix(2-d). It makes me want vectors more...
-Chris
On Mon, 28 Apr 2008, Christopher Barker apparently wrote:
I'm going to try to put down what I think are the key, very simple, questions:
For useful reference, these are now here: URL:http://www.scipy.org/MatrixIndexing#guidelines
Cheers, Alan
On Tue, Apr 29, 2008 at 12:20 AM, Alan G Isaac aisaac@american.edu wrote:
On Mon, 28 Apr 2008, Christopher Barker apparently wrote:
I'm going to try to put down what I think are the key, very simple, questions:
For useful reference, these are now here: URL:http://www.scipy.org/MatrixIndexing#guidelines
Cheers, Alan
May I add that if I edit defmatrix.py to act like an array for scalar indexing, then the following works.
In [1]: a = matrix(eye(2))
In [2]: array([a,a]) Out[2]: array([[[ 1., 0.], [ 0., 1.]],
[[ 1., 0.], [ 0., 1.]]])
This generates an error with the current version of matrix and, frankly, I am not going to be bothered going all through the numpy c sources to special case matrices to fix that. Someone else can do it if they wish. There are recursive routines that expect the dimensions to decrease on each call.
Chuck
Hi Charles
2008/4/29 Charles R Harris charlesr.harris@gmail.com:
May I add that if I edit defmatrix.py to act like an array for scalar indexing, then the following works.
In [1]: a = matrix(eye(2))
In [2]: array([a,a]) Out[2]: array([[[ 1., 0.], [ 0., 1.]],
[[ 1., 0.], [ 0., 1.]]])
This generates an error with the current version of matrix and, frankly, I am not going to be bothered going all through the numpy c sources to special case matrices to fix that. Someone else can do it if they wish. There are recursive routines that expect the dimensions to decrease on each call.
I'd also like to see matrices become proper hierarchical containers -- the question is just how to do that. Thus far, I'm most convinced by the arguments for RowVectors/Columns, which leaves us with a sane model for doing linear algebra, while providing the enhancements you mentioned here and in comments to another ticket.
We were thinking of raising a warning on scalar indexing for 1.1, but given the above, would that be sensical?
Regards Stéfan
On Tue, 29 Apr 2008, Stéfan van der Walt apparently wrote:
We were thinking of raising a warning on scalar indexing for 1.1, but given the above, would that be sensical?
There seem to be three basic proposals for scalar indexing:
- raise an error - return a 1d array - return a new type of 1d container (apparently to be called a "vector" for some reason ...)
I do not think anyone is still arguing in favor of the keeping the current behavior (where scalar indexing returns a matrix)?
So the question seems to be: how best to warn users that the behavior of scalar indexing is very likely to change.
Cheers, Alan Isaac
On Tue, Apr 29, 2008 at 7:24 AM, Stéfan van der Walt stefan@sun.ac.za wrote:
Hi Charles
2008/4/29 Charles R Harris charlesr.harris@gmail.com:
May I add that if I edit defmatrix.py to act like an array for scalar indexing, then the following works.
In [1]: a = matrix(eye(2))
In [2]: array([a,a]) Out[2]: array([[[ 1., 0.], [ 0., 1.]],
[[ 1., 0.], [ 0., 1.]]])
This generates an error with the current version of matrix and, frankly,
I
am not going to be bothered going all through the numpy c sources to
special
case matrices to fix that. Someone else can do it if they wish. There
are
recursive routines that expect the dimensions to decrease on each call.
I'd also like to see matrices become proper hierarchical containers -- the question is just how to do that. Thus far, I'm most convinced by the arguments for RowVectors/Columns, which leaves us with a sane model for doing linear algebra, while providing the enhancements you mentioned here and in comments to another ticket.
We were thinking of raising a warning on scalar indexing for 1.1, but given the above, would that be sensical?
Hi Stefan,
The numpy c routines call PySequence_GetItem(s,i) as well as ask for the length (first index), so I think these should be left as they are for arrays in order to guarantee that matrices are compatible with all the normal array operations. This means either returning special row objects that index as expected, or returning 1D arrays. I don't think the '*' operator has these problems, but in any case that is a well documented feature of matrices.
Chuck
Charles R Harris wrote:
On Tue, Apr 29, 2008 at 7:24 AM, Stéfan van der Walt <stefan@sun.ac.za mailto:stefan@sun.ac.za> wrote:
Hi Charles 2008/4/29 Charles R Harris <charlesr.harris@gmail.com <mailto:charlesr.harris@gmail.com>>: > May I add that if I edit defmatrix.py to act like an array for scalar > indexing, then the following works. > > In [1]: a = matrix(eye(2)) > > In [2]: array([a,a]) > Out[2]: > array([[[ 1., 0.], > [ 0., 1.]], > > [[ 1., 0.], > [ 0., 1.]]]) > > This generates an error with the current version of matrix and, frankly, I > am not going to be bothered going all through the numpy c sources to special > case matrices to fix that. Someone else can do it if they wish. There are > recursive routines that expect the dimensions to decrease on each call. I'd also like to see matrices become proper hierarchical containers -- the question is just how to do that. Thus far, I'm most convinced by the arguments for RowVectors/Columns, which leaves us with a sane model for doing linear algebra, while providing the enhancements you mentioned here and in comments to another ticket. We were thinking of raising a warning on scalar indexing for 1.1, but given the above, would that be sensical?
Hi Stefan,
The numpy c routines call PySequence_GetItem(s,i) as well as ask for the length (first index), so I think these should be left as they are for arrays in order to guarantee that matrices are compatible with all the normal array operations. This means either returning special row objects that index as expected or returning 1D arrays. I don't think the '*' operator has these problems, but in any case that is a well documented feature of matrices.
Thanks for looking in to this Chuck,
I'm quite persuaded now that a[i] should return a 1-d object for arrays. In addition to the places Chuck identified, there are at least 2 other places where special code was written to work-around the expectation that item selection returns an object with reduced dimensionality (a special-cased .tolist for matrices and a special-cased getitem in the fancy indexing code).
As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.
The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there. I also don't think that whatever problem may actually be buried in the fact that type(a[i]) != type(a[i,:]) is worse than the problem that several pieces of NumPy code actually expect hierarchical container behavior of multi-dimensional sequences.
I don't think making the small change to have a[i] return 1-d arrays precludes us from that 1-d array being a bit more formalized in 1.2 as a RowVector should we choose that direction. It does, however, fix several real bugs right now.
So, I'm
+1 on making the following change:
def __getitem__(self, index): + if isscalar(index): + return self.__array__()[index]
-Travis
On Tue, Apr 29, 2008 at 10:49:50AM -0500, Travis E. Oliphant wrote:
As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.
The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there.
I think we have simply replaced one problem by another one. I think this is a bad route to go. I think the only proposition that makes sens so far is the RowVector/ColumnVector. Yes you return a 1D object, but one that knows it is a linear algebra object, and not an array. I _really_ don't like a[i] != a[i,:]. I also don't like loosing the information that you are doing linear algebra.
Gaël
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
I really don't like a[i] != a[i,:].
Tim H's proposal avoids that problem. What do you think of it?
I also don't like loosing the information that you are doing linear algebra.
Hmmm. Is it a given that asking for ``a[i]`` is "doing linear algebra"? Is ``a[i:i+1,:]`` so bad if you want a matrix? I think some more use cases, like Bill's, might prove clarifying. URL:http://www.scipy.org/ConjugateGradientExample
Alan
On Tue, Apr 29, 2008 at 12:41 PM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
I really don't like a[i] != a[i,:].
Tim H's proposal avoids that problem. What do you think of it?
Of course, now we will have to remove all that special casing code...
Chuck
On Tue, Apr 29, 2008 at 02:41:45PM -0400, Alan G Isaac wrote:
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
I really don't like a[i] != a[i,:].
Tim H's proposal avoids that problem. What do you think of it?
Sorry, I guess I am lost. Could you remind me which one it is? Given how well you are summing up the discussion on http://www.scipy.org/MatrixIndexing , I suspect it is already in the list of proposal.
I also don't like loosing the information that you are doing linear algebra.
Hmmm. Is it a given that asking for ``a[i]`` is "doing linear algebra"? Is ``a[i:i+1,:]`` so bad if you want a matrix?
Very, very bad. This is completly unreadable. Actually to the beginner, this looks like a 2*N matrix, not a 1*N matrix.
To me you are really replacing a problem by another. We should fix the problem, not replace it.
Gaël
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
Could you remind me which one it is?
Proposal #6 at URL:http://www.scipy.org/MatrixIndexing#list-of-proposals is my effort to summarize what Tim was proposing. It does not include the oriented 1d vectors that I think you need (if I have understood you).
Alan
PS For beginners, perhaps ``x[[i],:]`` would be better than ``x[i:i+1,:]``. It has other advantages too...
On Tue, Apr 29, 2008 at 03:32:17PM -0400, Alan G Isaac wrote:
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
Could you remind me which one it is?
Proposal #6 at URL:http://www.scipy.org/MatrixIndexing#list-of-proposals is my effort to summarize what Tim was proposing. It does not include the oriented 1d vectors that I think you need (if I have understood you).
I prefer proposal number 1.
PS For beginners, perhaps ``x[[i],:]`` would be better than ``x[i:i+1,:]``. It has other advantages too...
I think this will get beginners lost. From my experience, people take a long time to understand indexing with an iterable. It is actually a source of confusion, but I am not even remotely thinking of removing it, because it is soooo nice, when you know how to use it.
Gaël
Hi everyone,
Despite being a bit lost in the matrix debate, today I was working on something which might want to use what is being described. You can see an array only version at:
http://fable.svn.sourceforge.net/viewvc/fable/ImageD11/trunk/test/demo/latre...
The problem is in diffraction from crystals. We measure "scattering vectors" which are in "reciprocal space" and use these to find crystal structures in "real space". These spaces are a covariant/contravariat pair*. The purpose of the code in the script is to construct a lattice class which can work with vectors that are directly measured, or come from an FFT of those vectors, which averages lots of peaks into fewer peaks. [nb, for a single lattice this is a solved problem in many software packages].
I used a keyword argument space="real" or space="reciprocal" to indicate which space a vector is in. It might be a good case to think about putting in RowVector and ColumnVector and trying to deal with the consequences. It is not in this code yet, but the lattice symmetry should show up soon. I have not understood how to distinguish a metric tensor (flips RowVector to ColumnVector) from a reciprocal metric tensor (flips back) from a simple rotation matrix (applies a symmetry operator, like 2-fold rotation, so doesn't flip at all). I fear the limitation is in my maths?
Best,
Jon
* Some heretics have been known to scale reciprocal space by 2pi. See "Vectors and Tensors in Crystallography" by Donald Sands for an overview.
Gael Varoquaux wrote:
On Tue, Apr 29, 2008 at 10:49:50AM -0500, Travis E. Oliphant wrote:
As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.
The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there.
I think we have simply replaced one problem by another one. I think this is a bad route to go.
The problem is that I know exactly what the problem we are replacing, but I have no idea only vagueries about the problem we are getting. So, in the short term, it seems like the right thing to do, because we aren't preventing a change to RowVector/ColumnVector in the future.
-Travis
On 29/04/2008, Travis E. Oliphant oliphant@enthought.com wrote:
I'm quite persuaded now that a[i] should return a 1-d object for arrays. In addition to the places Chuck identified, there are at least 2 other places where special code was written to work-around the expectation that item selection returns an object with reduced dimensionality (a special-cased .tolist for matrices and a special-cased getitem in the fancy indexing code).
As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.
The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there. I also don't think that whatever problem may actually be buried in the fact that type(a[i]) != type(a[i,:]) is worse than the problem that several pieces of NumPy code actually expect hierarchical container behavior of multi-dimensional sequences.
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object? Does this not require many special-case bits of code as well? What about a[i,...]? That is what I would use to make a hierarchical bit of code, and I would be startled to find myself in an infinite loop waiting for the dimension to become one.
I don't think making the small change to have a[i] return 1-d arrays precludes us from that 1-d array being a bit more formalized in 1.2 as a RowVector should we choose that direction. It does, however, fix several real bugs right now.
+1
What's more, I think we need to have a matrix-cleanliness test suite that verifies that all basic numpy tools behave identically on matrices and preserve matrixiness. But that's a big job, and for 1.2.
Anne
On 29/04/2008, Gael Varoquaux gael.varoquaux@normalesup.org wrote:
On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?
It breaks A*B[i, :] where A and B are matrices.
Really? How?
In [26]: A = np.matrix([[1,0],[0,1],[1,1]])
In [28]: A*np.ones(2) Out[28]: matrix([[ 1., 1., 2.]])
In [29]: np.ones(3)*A Out[29]: matrix([[ 2., 2.]])
I guess you don't get dot products with *:
In [30]: np.ones(3).T*np.ones(3) Out[30]: array([ 1., 1., 1.])
Since the whole point of matrices it to avoid having to type "dot" I guess this should be convincing.
Anne
On Tue, Apr 29, 2008 at 11:16:15PM +0200, Anne Archibald wrote:
On 29/04/2008, Gael Varoquaux gael.varoquaux@normalesup.org wrote:
On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?
It breaks A*B[i, :] where A and B are matrices.
Really? How?
In [26]: A = np.matrix([[1,0],[0,1],[1,1]])
In [28]: A*np.ones(2) Out[28]: matrix([[ 1., 1., 2.]])
In [29]: np.ones(3)*A Out[29]: matrix([[ 2., 2.]])
Yes, sorry, I am shameful. I should have thought a bit more before posting.
There is no big problem with x[i,:] not being a 1-d object. The problem is for x[:, i]. However I would find it nice that, for linear algebra, x[i, :] == x[:, i].T
This is the kind of behavior I expect, and we won't be getting it with 1D arrays.
Cheers,
Gaël
Gael Varoquaux wrote:
However I would find it nice that, for linear algebra, x[i, :] == x[:, i].T
This is the kind of behavior I expect, and we won't be getting it with 1D arrays.
But you WOULD get it with 1-d row/column objects.
I'm going to try to summarize what I think is clear from the discussion:
1) There are significant problems with the fact that the current matrix object does not reduce rank with indexing.
2) people want a "row" or "column" extracted from a matrix to behave like a linear algebra object, for example: M[:,i] * M[i,:] is an matrix
As far as I can see, the ONLY way to satisfy both of these is with a row/column object of some sort.
You can get almost there by having M[i] return a 1-d array, while M[:,i] and M[i,:] continue to return Matrices, but this gives an asymmetry between extracting rows and columns, and, probably more important, breaks:
M[i] == M[i,:]
plus you don't get the benefit of being able to extract a row or column and index the result with a scalar.
-Chris
On Tue, Apr 29, 2008 at 03:03:51PM -0700, Christopher Barker wrote:
Gael Varoquaux wrote:
However I would find it nice that, for linear algebra, x[i, :] == x[:, i].T
This is the kind of behavior I expect, and we won't be getting it with 1D arrays.
But you WOULD get it with 1-d row/column objects.
Yes. I agree 100% with you that this is the best solution. I fully agree with your complete message.
Gaël
On Tue, Apr 29, 2008 at 2:07 PM, Gael Varoquaux < gael.varoquaux@normalesup.org> wrote:
On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?
It breaks A*B[i, :] where A and B are matrices.
Shouldn't that be B[i,:] * A? Or am I just confused?
In any event this wouldn't be a problem under row/column scheme or any of the other results that end up tagging the results with some metainformation. B[i] would be tagged as a row vector somehow and __mul__ would "do the right thing" even though it was 1D.
On Tue, Apr 29, 2008 at 3:07 PM, Gael Varoquaux < gael.varoquaux@normalesup.org> wrote:
On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?
It breaks A*B[i, :] where A and B are matrices.
A*B[:,i], surely.
So it's a notational problem, what is the easiest way to get row/col vectors out of a matrix. And it will be helpful if the notation is consistent and easy to understand. In any case, dot doesn't care, so the fault above comes from the matrix class itself.
Chuck
It breaks A*B[i, :] where A and B are matrices.
A*B[:,i], surely.
right. which is a really good argument for:
A * B.col[i]
(and .row, or course)
-Chris
On Tue, 29 Apr 2008, Christopher Barker apparently wrote:
a really good argument for: A * B.col[i]
Also see the syntax discussed in Proposal 5. (I am not expressing an opinion.) URL:http://www.scipy.org/MatrixIndexing#proposal-5
One possibility is to let the rows and cols methods take an argument with a default of None. The default is to return an iterator over the rows/columns. An iterator over subsets of rows/columns could be specified with a sequence. A single row/column could be specified with a scalar argument.
Cheers, Alan
I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object? Does this not require many special-case bits of code as well? What about a[i,...]? That is what I would use to make a hierarchical bit of code, and I would be startled to find myself in an infinite loop waiting for the dimension to become one.
The rationale is so you can write
x[i,:] * A * x[:,i]
and have it work correctly without the RowVector / ColumnVector objects.
Maybe it's just easier to create them and be done with it. But, nobody has yet, so I'm not sure.
-Travis
On Tue, 29 Apr 2008, "Travis E. Oliphant" apparently wrote:
x[i,:] * A * x[:,i]
This need is also fully met by Proposal 5. It is just a matter of syntax and sticking with matrices. URL:http://www.scipy.org/MatrixIndexing#proposal-5
E.g., x.rows(i) * A * x.cols(i)
Cheers, Alan Isaac
On Tue, 29 Apr 2008, Charles R Harris apparently wrote:
There are recursive routines that expect the dimensions to decrease on each call.
I have added your example as anomaly #2 at the very top of the discussion page: URL:http://www.scipy.org/MatrixIndexing
Cheers, Alan Isaac
Alan G Isaac wrote:
For useful reference, these are now here: URL:http://www.scipy.org/MatrixIndexing#guidelines
Thanks Alan, Here are a few comments on that text (I know, I really should go dig up my scipy password...)
- Do we want to be able to extract a single row or column from a
matrix, and have it be indexable like a 1-d object? (Answer:
certainly > yes for rows, but columns could be gotten from the transpose.
I don't get this -- transpose of what? if you get a 1-d ndarray from: M[i,:], that can't be transposed -- it's 1-d. You could do M.T[0] but then you've lost the concept of it being a column. If we stick with M[i,:] being a (n,1) matrix, then it doesn't index like a 1-d object, and it doesn't need transposing anyway.
This is most useful, allows the most generic code, and breaks fewest expectations.)
- Do we want to be able to do (1) without calling the Array
attribute: M.A[i]? (Yes: this is not less useful, allows the most generic code, and breaks fewest expectations.)
- Do we want to have those 1-d objects to have an orientation so
that they act like Row or Column vectors in linear algebra operations (and broadcasting)? (Unclear: the resulting functionality could be provided by rows and columns attributes that yield corresponding matrices.)
no, it couldn't -- we can already use M[:,i] (or M[:,i:i+1] to get a matrix that is essentially a column -- however that can't be indexed with a scalar (see (1). That's kind of the point of this discussion. Or did you mean that you get a column by: np.matrix(M[:,i]).T ? yow!
Maybe we need to separate the row and column questions -- particularly with Charles' discovery of how other parts of numpy expect rank reduction with indexing, I think we probably have a consensus that we need to be able to extract 1-d objects from matrixes with indexing -- at least rows.
That leaves the column question: -- do we need an object that represents a column, that otherwise acts like a 1-d array (returns a scalar with scalar indexing)?
I guess my key point is that the fact that seeing a Matrix as a collection of rows rather than columns is an accident of C memory arrangement -- if numpy (and python) had been written in Fortran, we'd have the opposite. I'd rather that not be built into the API more than required. We can only have one definition for M[i], it might as well be a row, but other than that, we can have symmetry:
M[i,:] is row M[:,i] is a column
and rows and columns behave the same way everywhere they "should": indexing, converting to ndarrays, ???, but not in linear algebra operations.
Then, or course there are implementation details...
-Chris
On Tue, 29 Apr 2008, Christopher Barker apparently wrote:
Here are a few comments on that text
In response, I have modified URL:http://www.scipy.org/MatrixIndexing#guidelines to try to clarify what was meant and perhaps to address your core concerns.
Alan
PS Perhaps I should note that I have come around to what I take to be Tim H's view that matrix indexing should should work just like ndarray indexing, except that when it produces a 2d result then a matrix is returned. (This means giving up the current behavior of ``x[0,:]``.) I also agree with those who want ``rows`` and ``cols`` attributes for iteration. A small difference: I would have these yield matrices. (I.e., no need for a new "vector" type.)
Alan G Isaac wrote:
Perhaps I should note that I have come around to what I take to be Tim H's view that matrix indexing should should work just like ndarray indexing, except that when it produces a 2d result then a matrix is returned. (This means giving up the current behavior of ``x[0,:]``.) I also agree with those who want ``rows`` and ``cols`` attributes for iteration.
+1 on all of that.
A small difference: I would have these yield matrices. (I.e., no need for a new "vector" type.)
So the only thing we "give up" (which we don't have now anyway) is something that acts like column and can be indexed like a 1-d object. I still think that would be nice, but maybe we're coming close to a consensus about the rest of it.
-Chris
On Tue, Apr 29, 2008 at 10:11 AM, Alan G Isaac aisaac@american.edu wrote:
indexing should should work just like ndarray indexing, except that when it produces a 2d result then a matrix is returned. (This means giving up the current behavior of ``x[0,:]``.)
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
1. Could you give an example of the circumstances of this use? 2. Would any or all of the following be just as good a result? a. 1d matrix b. row and column vectors (1d but "oriented" for linear algebra purposes) 3. If 2a. and 2b. are no good for you, a. would e.g. ``x[i:i+1,:]`` be too burdensome? b. would ``rows`` and ``columns`` attributes that yielded matrics be an offsetting convenience?
Thank you, Alan Isaac
PS It is nice to have another matrix user in the conversation!
On Tue, Apr 29, 2008 at 11:21 AM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
- Could you give an example of the circumstances of this
use?
In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix
list index idx = [0] x[idx,:] --> matrix?
scalar index idx = 0 x[idx,:] --> not matrix
- Would any or all of the following be just as good
a result? a. 1d matrix b. row and column vectors (1d but "oriented" for linear algebra purposes)
For me, having arrays and matrices is confusing enough. Adding row and column objects sounds even more confusing.
- If 2a. and 2b. are no good for you, a. would e.g. ``x[i:i+1,:]`` be too burdensome?
Yes.
b. would ``rows`` and ``columns`` attributes that yielded matrics be an offsetting convenience?
I like the idea, from earlier in the thread, of row and col being iterators for those time when you want to work on one column or row at a time.
Let me throw out a couple of more thoughts:
First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).
Second, most of the stuff that we have been discussing so far is primarily about notational convenience. However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors). As a concrete example, I at times need to work with stacks of small matrices. If I do the operations one by one, the overhead is prohibitive, however, most of that overhead can be avoided. For example, I rewrote some of the linalg routines to work on stacks of matrices and inverse is seven times faster for a 100x10x10 array (a stack of 100 10x10 matrices) when operating on a stack than when operating on the matrices one at a time. This is a result of sharing the setup overhead, the C routines that called are the same in either case.
On Tue, Apr 29, 2008 at 12:22:18PM -0700, Timothy Hochberg wrote:
First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).
Good summary. I support the 1D object with orientation, rather than the 2D object with special indexing. I would call the row_vector/column_vectors bras and kets rather than tensors, but that because I come from a quantum mechanics background.
Second, most of the stuff that we have been discussing so far is primarily about notational convenience. However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors). As a concrete example, I at times need to work with stacks of small matrices. If I do the operations one by one, the overhead is prohibitive, however, most of that overhead can be avoided. For example, I rewrote some of the linalg routines to work on stacks of matrices and inverse is seven times faster for a 100x10x10 array (a stack of 100 10x10 matrices) when operating on a stack than when operating on the matrices one at a time. This is a result of sharing the setup overhead, the C routines that called are the same in either case.
Good point. Do you have an idea to move away from this problem?
Gaël
Hi,
On Tue, Apr 29, 2008 at 2:22 PM, Timothy Hochberg tim.hochberg@ieee.org wrote:
Let me throw out a couple of more thoughts:
First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).
Actually I think that this simply stems from the different variants of linear (and more specifically matrix) algebra. My background is very heavy in the statistical favor of this (what the heck are tensors... http://mathworld.wolfram.com/Tensor.html). It is also clear (from that definition) that others are from more physics and engineering backgrounds. Hence all the back and forth because people have slightly different usage, terminology and expectations.
The ability to treat vectors as matrices would be sufficient for my needs because these are almost always used in the context of vector-matrix multiplication. There is no additional benefit from having row or column shapes or metadata because the row/column nature is usually predetermined and would be represented by the shape of the corresponding matrix. It really would be annoying to find that for an n by 1 vector/matrix, the product of X.T*X is a scalar or an error rather than an n by n matrix. Or requires additional code to get the desired effect.
Regards Bruce
On Tue, 29 Apr 2008, Bruce Southey apparently wrote:
There is no additional benefit from having row or column shapes or metadata because the row/column nature is usually predetermined and would be represented by the shape of the corresponding matrix.
The problem is that ``x[0]`` being 2d has produced a variety of anomalies, and the natural fix is for ``x[0]`` to be 1d.
Gael has argued strongly that she should be able to use the following notation: ``x[0,:]*A*x[:,0]``. But this will work only if ``x[0,:]`` is 2d or if it is 1d but has an "orientation".
So *if* you think ``x[0] == x[0,:]`` is desirable, *and* you want to satisfy Gael, *then* it seems you must introduce 1d "oriented" vectors.
I believe Travis is also suggesting that we travel that road, taking a first step as follows: for now let ``x[0]`` be a 1d array to quickly fix the anomalies, but let ``x[0,:]`` continue to be a matrix until the vector code is written, at which point ``x[0]`` and ``x[0,:]`` we be the same "row vector".
Or so I have understood things.
Cheers, Alan Isaac
The problem is that ``x[0]`` being 2d has produced a variety of anomalies, and the natural fix is for ``x[0]`` to be 1d.
Gael has argued strongly that she should be able to use the following notation: ``x[0,:]*A*x[:,0]``. But this will work only if ``x[0,:]`` is 2d or if it is 1d but has an "orientation".
So *if* you think ``x[0] == x[0,:]`` is desirable, *and* you want to satisfy Gael, *then* it seems you must introduce 1d "oriented" vectors.
I believe Travis is also suggesting that we travel that road, taking a first step as follows: for now let ``x[0]`` be a 1d array to quickly fix the anomalies, but let ``x[0,:]`` continue to be a matrix until the vector code is written, at which point ``x[0]`` and ``x[0,:]`` we be the same "row vector".
Or so I have understood things.
You've characterized my current thinking pretty well. I'm less concerned that x[0] != x[0,:] than I think Gael is.
-Travis
On Tue, Apr 29, 2008 at 2:43 PM, Travis E. Oliphant oliphant@enthought.com wrote:
The problem is that ``x[0]`` being 2d has produced a variety of anomalies, and the natural fix is for ``x[0]`` to be 1d.
Gael has argued strongly that she should be able to use the following notation: ``x[0,:]*A*x[:,0]``. But this will work only if ``x[0,:]`` is 2d or if it is 1d but has an "orientation".
So *if* you think ``x[0] == x[0,:]`` is desirable, *and* you want to satisfy Gael, *then* it seems you must introduce 1d "oriented" vectors.
I believe Travis is also suggesting that we travel that road, taking a first step as follows: for now let ``x[0]`` be a 1d array to quickly fix the anomalies, but let ``x[0,:]`` continue to be a matrix until the vector code is written, at which point ``x[0]`` and ``x[0,:]`` we be the same "row vector".
Or so I have understood things.
You've characterized my current thinking pretty well. I'm less concerned that x[0] != x[0,:] than I think Gael is.
I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.
What if the extant matrix class would continue to be available for awhile as "oldmatrix", which you could then import as "matrix" if you wished. Would that meet your needs? (I am simply assuming this would be feasible, without being sure how a lot of the special casing for matrices has been done.)
Cheers, Alan Isaac
On Tue, Apr 29, 2008 at 5:18 PM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.
What if the extant matrix class would continue to be available for awhile as "oldmatrix", which you could then import as "matrix" if you wished. Would that meet your needs? (I am simply assuming this would be feasible, without being sure how a lot of the special casing for matrices has been done.)
I never import matrix directly (is that what you are suggesting?). I usually create it with M.ones, M.rand, x * y, etc., where M is
import numpy.matlib as M
If a big change is made to matrix behavior could it be accessed in 1.X from the __future__ and, if successful, switched to the "present" in 2.X? I, obviously, have no idea what would be involved to make that happen.
In my use, changing x[0] is not a big deal, but changing x[0,:] is.
On Apr 29, 2008, at 6:01 PM, Keith Goodman wrote:
I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.
The only code that should break would be indexing the extracted row/ column with two indexes. Multiplication with the extracted row/column would still work as expected (and that sounded to me like what the application typically was . . . maybe I'm wrong).
** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-0154 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **
Bill Spotz wrote:
On Apr 29, 2008, at 6:01 PM, Keith Goodman wrote:
break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.
The only code that should break would be indexing the extracted row/ column with two indexes.
It will also break code that expects M[i,:] to be the same shape as M[[i],:] and M[i:i+1,:], though it sounded like that was only really expected during testing.
-Chris
On Tue, 29 Apr 2008, Alan G Isaac apparently wrote:
Gael has argued strongly that he should be able to use the following notation: ``x[0,:]*A*x[:,0]``.
By the way Gael, is x.rows(0) * A * x.cols(0) a good replacement in your view? I find it nicely explicit and, to meet another of your concerns, usable by newbies.
Cheers, Alan
On Tue, Apr 29, 2008 at 06:44:11PM -0400, Alan G Isaac wrote:
On Tue, 29 Apr 2008, Alan G Isaac apparently wrote:
Gael has argued strongly that he should be able to use the following notation: ``x[0,:]*A*x[:,0]``.
By the way Gael, is x.rows(0) * A * x.cols(0) a good replacement in your view? I find it nicely explicit and, to meet another of your concerns, usable by newbies.
I don't really this proposal. I find this harder to read and make the matrices different than the arrays for no good reason.
We have a good solution (the row/column), why not use it?
Hum, I guess the answer is "talk is cheap, show me the code". So as I have hands tied with other work, I should rather shut up :) and let people who are actually doing the work decide.
Gaël
Bruce Southey wrote:
The ability to treat vectors as matrices would be sufficient for my needs because these are almost always used in the context of vector-matrix multiplication. There is no additional benefit from having row or column shapes or metadata because the row/column nature is usually predetermined and would be represented by the shape of the corresponding matrix.
The benefit is that they can be indexed with a scalar and converted to a 1-d array with r.A, and no reshaping. Also that indexing a matrix reduces its rank, which is expected in a lot of places.
It really would be annoying to find that for an n by 1 vector/matrix, the product of X.T*X is a scalar or an error rather than an n by n matrix.
yes, it would, which is the whole point of the matrix object in the first place, and the point of the proposed row/column objects. They would provide 1-d object that behave as row an column vectors with * and **.
-Chris
Timothy Hochberg wrote:
However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors).
Tim, this is important, but also appears to be an orthogonal issue to me -- whatever we do with matrices, rows, columns, whatever, we still need to solve this problem, and any of the schemes being proposed strikes me as amenable to having a "array or matrices", and the LA functions that act on them. What am I missing?
-Chris
On Tue, Apr 29, 2008 at 2:43 PM, Christopher Barker Chris.Barker@noaa.gov wrote:
Timothy Hochberg wrote:
However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors).
Tim, this is important, but also appears to be an orthogonal issue to me -- whatever we do with matrices, rows, columns, whatever, we still need to solve this problem, and any of the schemes being proposed strikes me as amenable to having a "array or matrices", and the LA functions that act on them. What am I missing?
Ah, that's a very good question. I'll do my best to answer it, it's been a long time since I dug into this seriously, so I may make some missteps: holler if you see something odd. Basically, I believe that you need two things:
1. The ability to make matrix/row/column objects of arbitrary dimensionality, where the extra dimensions contain the matrices/rows/columns. Current matrix objects are always 2D. 2. A way to tell whether a given object represents a row, column or matrix. For instance, given a 2D object is it matrix or a stack of vectors. You might think that there's enough information implicit in the shape of the arrays in conjunction with the nature of the function called that this wouldn't be necessary, but it turns out not to be so. As I recall, I had to fudge things when I was trying to containerize linalg.solve, because there wasn't enough information available to figure out unambiguously the kinds of the arguments.
Any of the solutions that attach information to the array, such as your matrix/row/col approach, or the quasi-tensor approach that I (and someone else) have mentioned but never detailed, should work for this.
However, any of the solutions that involve trying to encode the information into the array/matrix shape will probably not work. At least not well. Fortunately, they seem to be losing traction.
On Tue, Apr 29, 2008 at 1:22 PM, Timothy Hochberg tim.hochberg@ieee.org wrote:
Let me throw out a couple of more thoughts:
First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).
Second, most of the stuff that we have been discussing so far is primarily about notational convenience. However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors). As a concrete example, I at times need to work with stacks of small matrices. If I do the operations one by one, the overhead is prohibitive, however, most of that overhead can be avoided. For example, I rewrote some of the linalg routines to work on stacks of matrices and inverse is seven times faster for a 100x10x10 array (a stack of 100 10x10 matrices) when operating on a stack than when operating on the matrices one at a time. This is a result of sharing the setup overhead, the C routines that called are the same in either case.
Some operations on stacks of small matrices are easy to get, for instance, +,-,*,/, and matrix multiply. The last is the interesting one. If A and B are stacks of matrices with the same number of dimensions with the matrices stored in the last two indices, then
sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)
is the matrix-wise multiplication of the two stacks. If B is replaced by a stack of 1D vectors, x, it is even simpler:
sum(A[...,:,:]*x[...,newaxis,:], axis=-1)
This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.
Chuck
2008/4/30 Charles R Harris charlesr.harris@gmail.com:
Some operations on stacks of small matrices are easy to get, for instance, +,-,*,/, and matrix multiply. The last is the interesting one. If A and B are stacks of matrices with the same number of dimensions with the matrices stored in the last two indices, then
sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)
is the matrix-wise multiplication of the two stacks. If B is replaced by a stack of 1D vectors, x, it is even simpler:
sum(A[...,:,:]*x[...,newaxis,:], axis=-1)
This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.
Yes and no. For the first operation, you have to create a temporary that is larger than either of the two input arrays. These invisible (potentially) gigantic temporaries are the sort of thing that puzzle users when as their problem size grows they suddenly find they hit a massive slowdown because it starts swapping to disk, and then a failure because the temporary can't be allocated. This is one reason we have dot() and tensordot() even though they can be expressed like this. (The other is of course that it lets us use optimized BLAS.)
This rather misses the point of Timothy Hochberg's suggestion (as I understood it), though: yes, you can write the basic operations in numpy, in a more or less efficient fashion. But it would be very valuable for arrays to have some kind of metadata that let them keep track of which dimensions represented simple array storage and which represented components of a linear algebra object. Such metadata could make it possible to use, say, dot() as if it were a binary ufunc taking two matrices. That is, you could feed it two arrays of matrices, which it would broadcast to the same shape if necessary, and then it would compute the elementwise matrix product.
The question I have is, what is the right mathematical model for describing these arrays-some-of-whose-dimensions-represent-linear-algebra-objects?
One idea is for each dimension to be flagged as one of "replication", "vector", or "covector". A column vector might then be a rank-1 vector array, a row vector might be a rank-1 covector array, a linear operator might be a rank-2 object with one covector and one vector dimension, a bilinear form might be a rank-2 object with two covector dimensions. Dimensions designed for holding repetitions would be flagged as such, so that (for example) an image might be an array of shape (N,M,3) of types ("replication","replication","vector"); then to apply a color-conversion matrix one would simply use dot() (or "*" I suppose). without too much concern for which index was which. The problem is, while this formalism sounds good to me, with a background in differential geometry, if you only ever work in spaces with a canonical metric, the distinction between vector and covector may seem peculiar and be unhelpful.
Implementing such a thing need not be too difficult: start with a new subclass of ndarray which keeps a tuple of dimension types. Come up with an adequate set of operations on them, and implement them in terms of numpy's functions, taking advantage of the extra information about each axis. A few operations that spring to mind:
* Addition: it doesn't make sense to add vectors and covectors; raise an exception. Otherwise addition is always elementwise anyway. (How hard should addition work to match up corresponding dimensions?) * Multiplication: elementwise across "repetition" axes, it combines vector axes with corresponding covector axes to get some kind of generalized matrix product. (How is "corresponding" defined?) * Division: mostly doesn't make sense unless you have an array of scalars (I suppose it could compute matrix inverses?) * Exponentiation: very limited (though I suppose matrix powers could be implemented if the shapes are right) * Change of basis: this one is tricky because not all dimensions need come from the same vector space * Broadcasting: the rules may become a bit byzantine... * Dimension metadata fiddling
Is this a useful abstraction? It seems like one might run into trouble when dealing with arrays whose dimensions represent vectors from unrelated spaces.
Anne
On Wed, Apr 30, 2008 at 8:16 PM, Anne Archibald peridot.faceted@gmail.com wrote:
2008/4/30 Charles R Harris charlesr.harris@gmail.com:
Some operations on stacks of small matrices are easy to get, for
instance,
+,-,*,/, and matrix multiply. The last is the interesting one. If A and B are stacks of matrices with the same number of dimensions with the
matrices
stored in the last two indices, then
sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)
is the matrix-wise multiplication of the two stacks. If B is replaced by
a
stack of 1D vectors, x, it is even simpler:
sum(A[...,:,:]*x[...,newaxis,:], axis=-1)
This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.
Yes and no. For the first operation, you have to create a temporary that is larger than either of the two input arrays. These invisible (potentially) gigantic temporaries are the sort of thing that puzzle users when as their problem size grows they suddenly find they hit a massive slowdown because it starts swapping to disk, and then a failure because the temporary can't be allocated. This is one reason we have dot() and tensordot() even though they can be expressed like this. (The other is of course that it lets us use optimized BLAS.)
This rather misses the point of Timothy Hochberg's suggestion (as I understood it), though: yes, you can write the basic operations in numpy, in a more or less efficient fashion. But it would be very valuable for arrays to have some kind of metadata that let them keep track of which dimensions represented simple array storage and which represented components of a linear algebra object. Such metadata could make it possible to use, say, dot() as if it were a binary ufunc taking two matrices. That is, you could feed it two arrays of matrices, which it would broadcast to the same shape if necessary, and then it would compute the elementwise matrix product.
The question I have is, what is the right mathematical model for describing these arrays-some-of-whose-dimensions-represent-linear-algebra-objects?
One idea is for each dimension to be flagged as one of "replication", "vector", or "covector". A column vector might then be a rank-1 vector array, a row vector might be a rank-1 covector array, a linear operator might be a rank-2 object with one covector and one vector dimension, a bilinear form might be a rank-2 object with two covector dimensions. Dimensions designed for holding repetitions would be flagged as such, so that (for example) an image might be an array of shape (N,M,3) of types ("replication","replication","vector"); then to apply a color-conversion matrix one would simply use dot() (or "*" I suppose). without too much concern for which index was which. The problem is, while this formalism sounds good to me, with a background in differential geometry, if you only ever work in spaces with a canonical metric, the distinction between vector and covector may seem peculiar and be unhelpful.
Implementing such a thing need not be too difficult: start with a new subclass of ndarray which keeps a tuple of dimension types. Come up with an adequate set of operations on them, and implement them in terms of numpy's functions, taking advantage of the extra information about each axis. A few operations that spring to mind:
- Addition: it doesn't make sense to add vectors and covectors; raise
an exception. Otherwise addition is always elementwise anyway. (How hard should addition work to match up corresponding dimensions?)
- Multiplication: elementwise across "repetition" axes, it combines
vector axes with corresponding covector axes to get some kind of generalized matrix product. (How is "corresponding" defined?)
- Division: mostly doesn't make sense unless you have an array of
scalars (I suppose it could compute matrix inverses?)
- Exponentiation: very limited (though I suppose matrix powers could
be implemented if the shapes are right)
- Change of basis: this one is tricky because not all dimensions need
come from the same vector space
- Broadcasting: the rules may become a bit byzantine...
- Dimension metadata fiddling
Is this a useful abstraction? It seems like one might run into trouble when dealing with arrays whose dimensions represent vectors from unrelated spaces.
Thanks Anne. That is exactly what I had in mind. Alas, every time I sit down to try to prototype some code, it collapses under its own weight. I'm becoming warmer to an extended version of the row/col/matrix idea just because its simpler to understand. It would be just as you describe above but would support only four cases:
1. scalar: all axes are 'replication' axes 2. vector: last axes is a 'vector' all other replication. 3. covector: last axes is a 'covector' all other replication 4. matrix: last two axes are covector and vector respectively. Others are replication.
Or something like that. It's basically the row/col/matrix formulation with stacking. I suspect in practice it gives us most of the power of the full proposal with less complexity (for the user anyway). Then again, if the implementation turns out not to be that bad, one could always implement some convenience types on top of it to make it usable for the masses with the fully general version available for the stout of heart.
Regards,
On Wed, Apr 30, 2008 at 9:16 PM, Anne Archibald peridot.faceted@gmail.com wrote:
2008/4/30 Charles R Harris charlesr.harris@gmail.com:
Some operations on stacks of small matrices are easy to get, for
instance,
+,-,*,/, and matrix multiply. The last is the interesting one. If A and
B
are stacks of matrices with the same number of dimensions with the
matrices
stored in the last two indices, then
sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)
is the matrix-wise multiplication of the two stacks. If B is replaced by
a
stack of 1D vectors, x, it is even simpler:
sum(A[...,:,:]*x[...,newaxis,:], axis=-1)
This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.
Yes and no. For the first operation, you have to create a temporary that is larger than either of the two input arrays. These invisible (potentially) gigantic temporaries are the sort of thing that puzzle users when as their problem size grows they suddenly find they hit a massive slowdown because it starts swapping to disk, and then a failure because the temporary can't be allocated. This is one reason we have dot() and tensordot() even though they can be expressed like this. (The other is of course that it lets us use optimized BLAS.)
But it is interesting that you can multiply stacks of matrices that way, is it not? I haven't seen it mentioned elsewhere.
Chuck
On Tue, Apr 29, 2008 at 11:46 AM, Keith Goodman kwgoodman@gmail.com wrote:
On Tue, Apr 29, 2008 at 11:21 AM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
- Could you give an example of the circumstances of this
use?
In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix
list index idx = [0] x[idx,:] --> matrix?
scalar index idx = 0 x[idx,:] --> not matrix
Here's another use case:
for i in xrange(x.shape[1]): xi = x[:,i] idx = M.where(M.isfinite(xi).A)[0] xi = xi[idx,:] # more code
Also the functions I have written work on matrices. If x[i,:] did not return a matrix, let's say it returned some other vector type object that you index with a scalar, then my functions will have to test whether the input is a full matrix or a vector since the indexing is different.
On Tue, Apr 29, 2008 at 12:50 PM, Keith Goodman kwgoodman@gmail.com wrote:
On Tue, Apr 29, 2008 at 11:46 AM, Keith Goodman kwgoodman@gmail.com wrote:
On Tue, Apr 29, 2008 at 11:21 AM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
- Could you give an example of the circumstances of this
use?
In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix
list index idx = [0] x[idx,:] --> matrix?
scalar index idx = 0 x[idx,:] --> not matrix
Here's another use case:
for i in xrange(x.shape[1]): xi = x[:,i] idx = M.where(M.isfinite(xi).A)[0] xi = xi[idx,:] # more code
Also the functions I have written work on matrices. If x[i,:] did not return a matrix, let's say it returned some other vector type object that you index with a scalar, then my functions will have to test whether the input is a full matrix or a vector since the indexing is different.
And here's a use case that currently doesn't work but would be nice (for me) if it did:
x[idx,i] = M.rand(2,1)
where x is a matrix, say 4x3, idx is an array with 2 elements, and i is a scalar
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
here's a use case that currently doesn't work but would be nice (for me) if it did: x[idx,i] = M.rand(2,1)
I think this is another example of the kind of thing we are trying to fix ...
Cheers, Alan Isaac
Keith Goodman wrote:
It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
I think this is where matrices should act like arrays:
a
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
a[:,[1]]
array([[1], [4], [7]])
a[:,1]
array([1, 4, 7])
in this case, it's not really different types, but totally different concepts -- indexing with a sequence is different than a scalar, just like indexing with a slice is different than a scalar.
There are many, many places in python code where we have to make distinctions between a single object and a sequence of objects that happens to be of length one.
-Chris
On 29/04/2008, Keith Goodman kwgoodman@gmail.com wrote:
On Tue, Apr 29, 2008 at 11:21 AM, Alan G Isaac aisaac@american.edu wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
- Could you give an example of the circumstances of this
use?
In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix
list index idx = [0] x[idx,:] --> matrix?
scalar index idx = 0 x[idx,:] --> not matrix
It is actually pretty unreasonable to hope that
A[0]
and
A[[1,2,3]] or A[[True,False,True]]
should return objects of the same rank.
Anne
On Tue, Apr 29, 2008 at 2:18 PM, Anne Archibald peridot.faceted@gmail.com wrote:
On 29/04/2008, Keith Goodman kwgoodman@gmail.com wrote:
In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:
array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix
list index idx = [0] x[idx,:] --> matrix?
scalar index idx = 0 x[idx,:] --> not matrix
It is actually pretty unreasonable to hope that
A[0]
and
A[[1,2,3]] or A[[True,False,True]]
should return objects of the same rank.
Why it unreasonable to hope that
x[0,:]
and
x[0, [1,2,3]]
or
x[0, [True,False,True]]
where x is a matrix, continue to return matrices?
On 30/04/2008, Keith Goodman kwgoodman@gmail.com wrote:
On Tue, Apr 29, 2008 at 2:18 PM, Anne Archibald peridot.faceted@gmail.com wrote:
It is actually pretty unreasonable to hope that
A[0]
and
A[[1,2,3]] or A[[True,False,True]]
should return objects of the same rank.
Why it unreasonable to hope that
x[0,:]
and
x[0, [1,2,3]]
or
x[0, [True,False,True]]
where x is a matrix, continue to return matrices?
Well, I will point out that your example is somewhat different from mine; nobody is arguing that your three examples should return objects of different ranks. There is some disagreement over whether
x[0,:]
should be a rank-1 or a rank-2 object, but x[0,[1,2,3]] and x[0, [True, False, True]] should all have the same rank as x[0,:]. Nobody's questioning that.
What I was pointing out is that x[0,0] should not have the same rank as x[0,:] or x[0,[0]]. In this case it's obvious; x[0,0] should be a scalar. But the same logic applies to x[0,:] versus x[[0,1],:] or even x[[0],:].
For arrays, if x has rank 2, x[0,:] has rank 1. if L is a list of indices, x[L,:] always has rank 2, no matter what is actually in L (even if it's one element). It is perfectly reasonable that they should yield objects of different ranks.
With matrices, we have the surprising result that x[0,:] is still a two-dimensional object; to get at its third element I have to specify x[0,:][0,2]. It seems to me that this is a peculiar hack designed to ensure that the result still has "*" defined in a matrix way.
Anne
I will answer this question too, because I feel like Keith.
On Tue, Apr 29, 2008 at 02:21:34PM -0400, Alan G Isaac wrote:
On Tue, 29 Apr 2008, Keith Goodman apparently wrote:
I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.
- Could you give an example of the circumstances of this
use?
x[i, :]*A*x[:, i]
I find this cleaner, because more explicit, than x[i]*A*x[:, i]. I'd be _very very_ upset if this broke. Breaking the explicit for the implicit is a _bad_ idea, and if you go down this path you go down the path of inconsistency.
- Would any or all of the following be just as good
a result? a. 1d matrix
Yes.
b. row and column vectors (1d but "oriented" for linear algebra purposes)
Yes.
- If 2a. and 2b. are no good for you, a. would e.g. ``x[i:i+1,:]`` be too burdensome?
No way. This is hard to read and the best way to make long calculations impossible to follow and debug.
b. would ``rows`` and ``columns`` attributes that yielded matrics be an offsetting convenience?
Hum, I don't love this because it introduces noise in the formulas, but at least it is explicit and readable.
I would really like us to go down the way of the row/vector columns. This is consistent at least.
Gaël
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
x[i, :]*A*x[:, i] I find this cleaner, because more explicit, than x[i]*A*x[:, i].
a. would e.g. ``x[i:i+1,:]``
No way.
So *if* ``x[0]`` is going to be 1d to fix the many anomalies that have surfaced, then you would still want ``x[0]==x[0,:]``, but to have those be 1d and have ``x[i, :]*A*x[:, i]`` work would require (oriented) "vector" objects, so you would want those.
Is that a fair summary of your view?
Cheers, Alan Isaac
On Tue, Apr 29, 2008 at 03:20:02PM -0400, Alan G Isaac wrote:
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
x[i, :]*A*x[:, i] I find this cleaner, because more explicit, than x[i]*A*x[:, i].
a. would e.g. ``x[i:i+1,:]``
No way.
So *if* ``x[0]`` is going to be 1d to fix the many anomalies that have surfaced, then you would still want ``x[0]==x[0,:]``, but to have those be 1d and have ``x[i, :]*A*x[:, i]`` work would require (oriented) "vector" objects, so you would want those.
Is that a fair summary of your view?
That an extremely good summary of my views. Better than I could have formulated it. Thank you.
Given this good short paragraph summing up the problem, as I see it, this is why I think the only correct option is the oriented 1D objects. The second best option is not to change anything in order not to break backward compatibility.
Gaël
On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:
I think the only correct option is the oriented 1D objects. The second best option is not to change anything in order not to break backward compatibility.
Maybe the extant ``matrix`` class should be available for a while as ``oldmatrix``?
Alan
2008/4/25 Alan G Isaac aisaac@american.edu:
- This is **not** what I understood as the agreement
(and I think the current solution is bad).
Reverted in r5084.
Cheers Stéfan
On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:
Reverted in r5084.
Thank you.
I think we have discovered that there is a basic conflict between two behaviors:
x[0] == x[0,:] vs. x[0][0] == x[0,0]
To my recollection, everyone has agree that the second behavior is desirable as a *basic expectation* about the behavior of 2d objects. This implies that *eventually* we will have ``x[0][0] == x[0,0]``. But then eventually we MUST eventually have x[0] != x[0,:].
Since a current behavior must disappear eventually, we should make it disappear as soon as possible: before the release. The question is how. I see two simple ways to move forward immediately:
1. scalar indexing of matrices returns a 1d array 2. scalar indexing of matrices raises a TypeError
If we raise a TypeError, we lose both behaviors. If we return a 1d array, we retain the behavior that has been seen as desirable (as a *basic expectation* about the behavior of 2d objects).
Stefan, do you agree, and if so, would you be willing to have a vote on these two options?
Alan
On Fri, Apr 25, 2008 at 12:02 PM, Alan G Isaac aisaac@american.edu wrote:
On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:
Reverted in r5084.
Thank you.
I think we have discovered that there is a basic conflict between two behaviors:
x[0] == x[0,:] vs. x[0][0] == x[0,0]
To my recollection, everyone has agree that the second behavior is desirable as a *basic expectation* about the behavior of 2d objects. This implies that *eventually* we will have ``x[0][0] == x[0,0]``. But then eventually we MUST eventually have x[0] != x[0,:].
Choices:
1) x[0] equiv x[0:1,:] -- current behavior 2) x[0] equiv x[0,:] -- array behavior, gives x[0][0] == x[0,0]
These are incompatible. I think 1) is more convenient in the matrix domain and should be kept, while 2) should be given up. What is the reason for wanting 2)? Maybe we could overload the () operator so that x(0) gives 2)? Or vice versa.
Chuck
On 25/04/2008, Charles R Harris charlesr.harris@gmail.com wrote:
On Fri, Apr 25, 2008 at 12:02 PM, Alan G Isaac aisaac@american.edu wrote:
I think we have discovered that there is a basic conflict between two behaviors:
x[0] == x[0,:] vs. x[0][0] == x[0,0]
To my recollection, everyone has agree that the second behavior is desirable as a *basic expectation* about the behavior of 2d objects. This implies that *eventually* we will have ``x[0][0] == x[0,0]``. But then eventually we MUST eventually have x[0] != x[0,:].
Choices:
- x[0] equiv x[0:1,:] -- current behavior
- x[0] equiv x[0,:] -- array behavior, gives x[0][0] == x[0,0]
These are incompatible. I think 1) is more convenient in the matrix domain and should be kept, while 2) should be given up. What is the reason for wanting 2)? Maybe we could overload the () operator so that x(0) gives 2)? Or vice versa.
We are currently operating in a fog of abstraction, trying to decide which operations are more "natural" or more in correspondence with our imagined idea of a user. We're not getting anywhere. Let's start writing up (perhaps on the matrix indexing wiki page) plausible use cases for scalar indexing, and how they would look under each proposal.
For example:
* Iterating over the rows of a matrix:
# array-return proposal; x is a scalar for row in M: if np.all(row==0): raise ValueError # exception-raising proposal for i in xrange(M.shape[0]): if np.all(M[i,:]==0): raise ValueError
Other applications?
Anne
Alan G Isaac wrote:
Please return 1d arrays in response to scalar indexing as the provisional arrangement.
+1 (for the provisional solution)
This has clarified it a bit for me. The current situation allows one to create "row vectors" and "column vectors" by generating matrices (in various ways) that happen to be shape (1,n) or (n,1). These objects behave as we all want in arithmetic operations (*, **). So what's missing:
***** The ability to index one of these "vector" objects by a single index, and get a scalar as a result. *****
If I'm not mistaken (and I probably am) -- that is the crux of this entire conversation.
Do we want that? I would say yes, it's a pretty natural and common thing to want.
Note that with the current version of matrix (numpy 1.1.0rc1), you can index a (n,1) matrix with a scalar, and get a (1,1) matrix, but you get an exception if you try to index a (1,n) matrix with a scalar:
rv
matrix([[3, 4, 5]])
rv[1]
Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/numpy/core/defmatrix.py", line 228, in __getitem__ out = N.ndarray.__getitem__(self, index) IndexError: index out of bounds
I know why that's happening, but I think it's a demonstration of why we might want some sort of "vector" object, for a cleaner API.
How do we achieve this?
One possibility is to special case (n,1) and (1,n) matrices, but I think there is consensus that that is a BAD IDEA (tm)
Which leaves us with needing a "vector" object -- which acts a lot like a 1-d array, except that it acts like a (n,1) or (1,n) matrix with * and **, and may get printed differently.
I think there are two ways to accomplish this: (1) RowVector or ColumnVector objects that are essentially 1-d ndarrays, with a few methods overridden.
(2) A Vector object that is a matrix of either (n,1) or (1,n) in shape, with a couple methods overridden.
I think the difference is an implementation detail. They should behave about the same either way.
In either case, we should probably have factory functions for these objects: RowVector() and ColumnVector()
By the way, I think this brings us closer to the goal of matrices "acting as much like arrays as possible" -- when you index into a 2-d array: A[:,i], you get a 1-d array. These "vectors" would act a lot like 1-d arrays, so that would be more similar than the current situation.
It also might be nice to have iterators over rows and columns, like:
for row in M.rows: ... and for col in M.columns: ...
We can get the rows as the default iterator, but it feels nice and symmetric to be able to iterate over columns too (though I have no idea how useful that would really be)
Stéfan van der Walt wrote:
The whole idea of matrices was that you always work with 2-dimensional arrays, even when you just extract a row. Unless you have a proper hierarchical container (as illustrated in the other patch I sent), returning a 1D array breaks at least some expectation.
I agree -- stopgap is OK (but having people do M.A[i] is a fine stopgap too), but if we really want this, we need some sort of "vector" object.
If that's what the matrix users want, though, then we should change it.
I still want to see some linear algebra examples -- that's what matrices are for. We can iterate over 2-d arrays and get 1-d arrays just fine now, if that's what you want.
Unless we agree on something, and soon, that won't happen. Your workaround would break x[0] == x[0,:], so we're just swapping one set of broken functionality for another.
Alan G Isaac wrote:
But then we MUST eventually have x[0] != x[0,:].
I don't think so -- I think a Vector object would allow both of:
M[i,j] == M[i][j] and M[i] == M[i,:]
however,
M[i] != M[:,i] -- M[i] can only mean one thing, and this doesn't hold for arrays, either.
and (if we do the iterators):
M.rows[i] == M[i,:] M.columns[i] == M[:,i]
You can get rows a bit more conveniently than columns, which is really just an accident of how C arranges memory, but why not? Again, to make this "pure" we'd disallow M[i], but practicality beats purity, after all.
how would you index into a vector (as in http://en.wikipedia.org/wiki/Vector_(spatial)) without it?
Right -- which is why a 1-d "vector" object has its uses.
Alan G Isaac wrote:
Since a current behavior must disappear eventually, we should make it disappear as soon as possible: before the release. The question is how. I see two simple ways to move forward immediately:
1. scalar indexing of matrices returns a 1d array 2. scalar indexing of matrices raises a TypeError
I don't agree -- what we have is what we have, not changing it will not break any existing code, so we should only make a change if it moves us closer to what we want in the future -- it does seem that most folks do want M[i] to return some sort of 1-d object, so I'm fine with (1). I'm also fine with (2), but I'm pretty sure that one was laughed out of the discussion a good while back!
-Chris
On Fri, 25 Apr 2008, Christopher Barker apparently wrote:
I think a Vector object would allow both of: M[i,j] == M[i][j] and M[i] == M[i,:]
The problem is that it would be a crime to give up the natural production of submatrices. The NATURAL RULE is: to get a submatrix, use nonscalar indices. We should *not* give up that x[0,:] is a sub*matrix* whose first element is x[0,0] and equivalently x[0][0]. *This* is why we must have x[0]!=x[0,:] if we want, as we do, that x[0][0]==x[0,0].
Note that the idea for attributes ``rows`` and ``columns`` is contained on the discussion page: URL:http://www.scipy.org/MatrixIndexing I claim that it is natural for these attributes to yield properly shaped *matrices*. Once we have these attributes, it is difficult to see what we gain by introducing the complication of a separate vector class.
I still see no real gain from a separate vector class (or row and column vector classes). Everything that has been proposed to do with these is achievable with the corresponding matrices with one exception: scalar indexing -> element so that x[0][0]==x[0,0]. But this outcome is achieved MUCH more simply by letting x[0] be a 1d array.
Anyway, I am glad you agree that letting ``x[0]`` be a 1d array is the proper provisional solution.
Cheers, Alan
On Fri, Apr 25, 2008 at 2:57 PM, Alan G Isaac aisaac@american.edu wrote:
On Fri, 25 Apr 2008, Christopher Barker apparently wrote:
I think a Vector object would allow both of: M[i,j] == M[i][j] and M[i] == M[i,:]
The problem is that it would be a crime to give up the natural production of submatrices. The NATURAL RULE is: to get a submatrix, use nonscalar indices. We should *not* give up that x[0,:] is a sub*matrix* whose first element is x[0,0] and equivalently x[0][0]. *This* is why we must have x[0]!=x[0,:] if we want, as we do, that x[0][0]==x[0,0].
I would give it up and let x(0) be a submatrix, the normal indexing is then exactly like an array. True, the result of x(0) can't be used as an lvalue (only getitem equivalence), but otherwise it should work fine.
Chuck
On Fri, 25 Apr 2008, Christopher Barker apparently wrote:
I think a Vector object would allow both of: M[i,j] == M[i][j] and M[i] == M[i,:]
On Fri, Apr 25, 2008 at 2:57 PM, Alan G Isaac aisaac@american.edu wrote:
The problem is that it would be a crime to give up the natural production of submatrices. The NATURAL RULE is: to get a submatrix, use nonscalar indices. We should not give up that x[0,:] is a sub*matrix* whose first element is x[0,0] and equivalently x[0][0]. This is why we must have x[0]!=x[0,:] if we want, as we do, that x[0][0]==x[0,0].
On Fri, 25 Apr 2008, Charles R Harris apparently wrote:
I would give it up and let x(0) be a submatrix, the normal indexing is then exactly like an array. True, the result of x(0) can't be used as an lvalue (only getitem equivalence), but otherwise it should work fine.
The "it" here is ambiguous. Which do you want to give up?
- x[0][0]==x[0,0] ? the argument has been that this is a fundamental expectations for 2d array-like objects, and I think that argument is right
- x[0]==x[0,:] ? I believe you mean we should give this up, and if so, I strongly agree. As you know. Enforcing this is proving untenable and costly.
Cheers, Alan Isaac
OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?
Those who have spoke up have generally favored letting x[0] return a 1d array, if I count correctly. And I think that is the right "provisional" behavior. (As well as ultimate.) But many have not spoken up.
So let's at least signal that this is an area of change. Right?
Thank you, Alan Isaac
On Fri, Apr 25, 2008 at 6:38 PM, Alan G Isaac aisaac@american.edu wrote:
OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?
Those who have spoke up have generally favored letting x[0] return a 1d array, if I count correctly. And I think that is the right "provisional" behavior. (As well as ultimate.) But many have not spoken up.
That would be the most compatible with arrays, but won't it break a lot of code?
Chuck
On Sat, Apr 26, 2008 at 10:27 AM, Charles R Harris charlesr.harris@gmail.com wrote:
On Fri, Apr 25, 2008 at 6:38 PM, Alan G Isaac aisaac@american.edu wrote:
OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?
Those who have spoke up have generally favored letting x[0] return a 1d array, if I count correctly. And I think that is the right "provisional" behavior. (As well as ultimate.) But many have not spoken up.
That would be the most compatible with arrays, but won't it break a lot of code?
Any change to the behavior of x[0] for matrices is going to break a lot of code. It actually seems like a good idea to me to make it an error for a while -- or maybe some kind of settable option -- to help people figure out where they need to fix their code. Better to get an error immediately.
--bb
Alan G Isaac wrote:
OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?
I think this is wise for the time being.
At this point, I'm leaning in the direction of the RowVector / ColumnVector approach (even if these are not really advertised and just used during indexing).
-Travis
On Fri, 25 Apr 2008, "Travis E. Oliphant" apparently wrote:
At this point, I'm leaning in the direction of the RowVector / ColumnVector approach (even if these are not really advertised and just used during indexing).
I believe that this conflicts with submatrix extraction. Details follow.
Suppose ``x`` is a matrix, and I say ``v=A[0]``. What is ``v``?
1. If ``v`` is a 1d array, its behavior is known. E.g., ``v[0]`` is the first element and ``v[0,0]`` is an IndexError. If we want to keep submatrix extraction (as we should), we give up ``x[0] == x[0,:]``.
2. If ``v`` is a matrix, its behavior can be learned by the experienced but breaks basic expectations (the x[0][0] problem).
3. If ``v`` is a "RowVector" what behavior do we get? Suppose the idea is that this will rescue ``x[0][0]==x[0,0]`` **and** ``x[0]=x[0,:]``. But then we must give up that ``x[0,:]`` is a submatrix. Must ``v`` deviate from submatrix behavior in an important way? Yes: ``v[0][0]`` is an IndexError. (Note that the same problem arises if we just override ``__getitem__`` for the matrix class to special case row and column matrices.)
Since submatrix extraction is fundamental, I think it is a *very* bad idea to give it up. If so, then under the RowVector proposal we must give up ``x[0]==x[0,:]``. But this is just what we give up with the much simpler proposal that ``v`` be a 1d array.
Cheers, Alan
On Sat, Apr 26, 2008 at 7:12 AM, Alan G Isaac aisaac@american.edu wrote:
On Fri, 25 Apr 2008, "Travis E. Oliphant" apparently wrote:
At this point, I'm leaning in the direction of the RowVector / ColumnVector approach (even if these are not really advertised and just used during indexing).
I believe that this conflicts with submatrix extraction. Details follow.
Suppose ``x`` is a matrix, and I say ``v=A[0]``. What is ``v``?
- If ``v`` is a 1d array, its behavior is known.
E.g., ``v[0]`` is the first element and ``v[0,0]`` is an IndexError. If we want to keep submatrix extraction (as we should), we give up ``x[0] == x[0,:]``.
- If ``v`` is a matrix, its behavior can be learned
by the experienced but breaks basic expectations (the x[0][0] problem).
- If ``v`` is a "RowVector" what behavior do we get?
Suppose the idea is that this will rescue ``x[0][0]==x[0,0]`` **and** ``x[0]=x[0,:]``. But then we must give up that ``x[0,:]`` is a submatrix. Must ``v`` deviate from submatrix behavior in an important way? Yes: ``v[0][0]`` is an IndexError. (Note that the same problem arises if we just override ``__getitem__`` for the matrix class to special case row and column matrices.)
Since submatrix extraction is fundamental, I think it is a *very* bad idea to give it up. If so, then under the RowVector proposal we must give up ``x[0]==x[0,:]``. But this is just what we give up with the much simpler proposal that ``v`` be a 1d array.
I may have missed something since I've been getting these emails all out of order for some reason and it's been hard to follow. Can you clarify what you mean by submatrix extraction? It sounds like you want to be able index into an MxN array and get out a 1xN or Mx1 matrix. If that's the case, wouldn't the natural way to spell that under the RowVector/ColumnVector approach (and most others that I can think of) be:
m2 = m1[:1]
or
m2 = m1[:,:1]
Here we're extracting the first row and column respectively as a matrix. This doesn't appear to have any of the conflicts with row/vector extraction and in general meshes much better with arrays than having some magic associated with appending an extra ':', which seems ill advised to me.
On Mon, 28 Apr 2008, Timothy Hochberg apparently wrote:
Can you clarify what you mean by submatrix extraction? It sounds like you want to be able index into an MxN array and get out a 1xN or Mx1 matrix. If that's the case, wouldn't the natural way to spell that under the RowVector/ColumnVector approach (and most others that I can think of) be: m2 = m1[:1] or m2 = m1[:,:1]
I have tried to capture your view in proposal #6 at URL:http://www.scipy.org/MatrixIndexing
Cheers, Alan
PS I would **love** to see things go this way! More consistency with ndarray behavior is better!! (Note: this conflicts with current behavior, where ``x[0]`` and ``x[0,:]`` return matrices.)
Alan G Isaac wrote:
OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?
On Fri, 25 Apr 2008, "Travis E. Oliphant" apparently wrote:
I think this is wise for the time being.
In response to Bill's concern, I chose a warning instead: URL:http://scipy.org/scipy/numpy/attachment/ticket/760
fwiw, Alan
On Fri, Apr 25, 2008 at 01:04:53PM -0400, Alan G Isaac wrote:
On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:
The agreement was: a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)
- This is **not** what I understood as the agreement
(and I think the current solution is bad). I certainly did not mean to support the change that as implemented, and it is not clear to me that others did either.
- These two goals are substantially in conflict for
matrices, as we are seeing.
- The important goal, (goal a., which everyone agrees on),
has NOT been accomplished by the current change: x[0][0] raises a TypeError when x is a 1 by N matrix.
I claim b is more important than a. IMHO, a is plain wrong: you should't be indexing x with x[0][0]. I am OK making a work, as long as it doesn't break b. In addition breaking be is backward incompatible changes for no good reason (replacing one missfeature with another).
I would like this issue addressed in next release, not this one. I have the feeling the discussion is not sane enough, right now.
Cheers,
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
I claim b is more important than a. IMHO, a is plain wrong: you should't be indexing x with x[0][0].
Why?? Would you say this about a 2d array? Why the difference?
The core argument has been that it is a **basic expectation** of the behavior of 2d array-like objects that you will be able to get, e.g., the first element with x[0][0]. (E.g., lists, tuples, arrays ... is there an exception?? It should not be broken!)
I teach with matrices and I thought you might too: if so, you surely have run into this expectation (which is **natural**). In fact I truly doubt anyone has not been puzzled on first encounter by this::
>>> x matrix([[1, 2], [3, 4]]) >>> x[0] matrix([[1, 2]]) >>> x[0][0] matrix([[1, 2]])
Another argument, which I agree with, has been that matrix behavior should match 2d array behavior except where there is a clear payoff from deviation.
What I do agree with, which I take to be implicit in your statement, is that if users of matrices should generally use multiple indexes. But that does not answer the question about what we should offer when x[0] is requested. (Or, more crucially, what you get when you iterate over a matrix.)
Cheers, Alan Isaac
PS If you have good arguments, please add them to URL:http://www.scipy.org/MatrixIndexing
On Sat, Apr 26, 2008 at 11:13:12AM -0400, Alan G Isaac wrote:
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
I claim b is more important than a. IMHO, a is plain wrong: you should't be indexing x with x[0][0].
Why??
Because a 2D object is not a list of list. It is more than that. The numpy array actually exposes this interface, but I think using it systematicaly when it is not required is abuse. Moreover you are creating temporary objects that are useless and harmful. They are harmful for performance: you have to create useless objects and call their __getitem__ functions.
Would you say this about a 2d array?
Of course. Even more. Suppose you have an nd array, with n=10, to index object using 10 different indexing, ie A[a][b]... you have to create 9 intermediate objects that get indexed. This is ridiculous and potentially very harmful for performance.
The core argument has been that it is a **basic expectation** of the behavior of 2d array-like objects that you will be able to get, e.g., the first element with x[0][0]. (E.g., lists, tuples, arrays ... is there an exception??
For me this is wrong. list and tuples are not 2D. Numpy arrays happen to offer this feature, but you should not use it do to multiple dimension indexing.
I teach with matrices and I thought you might too: if so, you surely have run into this expectation (which is **natural**).
Well, I don't find it that natural. I would naturally expect A[0][0] to be invalid. I have not heavily run into that expectation. People around me are used to indexing with several indexes. This probably comes from their fortran/Mathematica/Maple/Matlab background. I can see that someone coming from C would expect this multiple consecutive indexing. This is because C has no notion of multiple dimension objects. The C model is limited and not as powerful as modern array programming languages, let us just move along.
In fact I truly doubt anyone has not been puzzled on first encounter by this::
>>> x matrix([[1, 2], [3, 4]]) >>> x[0] matrix([[1, 2]]) >>> x[0][0] matrix([[1, 2]])
Well, actually, I can say that this behavior is surprising. When I teach numpy and somebody encounters this behavior (which doesn't happen often, because people use multiple dimension indexing), I have to point out that x[0] is a 2D object also.
I sympathize with your crusade to fix this inconvenience. This is why I support the RowVector/ColumnVector proposal. However your proposal breaks what I consider as normal and important for something I consider should be avoided.
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
For me this is wrong. list and tuples are not 2D. Numpy arrays happen to offer this feature, but you should not use it do to multiple dimension indexing.
But there is no proposal that people should index like this. The underlying issue is iteration. Currently::
>>> A matrix([[1, 2], [ 3, 4]]) >>> A = np.mat(x) >>> for row in A: ... for col in row: ... print col ... [[1 2]] [[3 4]]
So are you saying that one should not be able to iterate through to the elements of a matrix?
So many natural expectations are broken by the current behavior!
Cheers, Alan
On Sat, Apr 26, 2008 at 02:01:45PM -0400, Alan G Isaac wrote:
>>> A matrix([[1, 2], [ 3, 4]]) >>> A = np.mat(x) >>> for row in A: ... for col in row: ... print col ... [[1 2]] [[3 4]]
So are you saying that one should not be able to iterate through to the elements of a matrix?
I kinda expect this was the real issue.
I think the only way to get consistent behavior for such code is either
* to use a syntax similar to dictionnaries:
for row in A.rows(): for col in row.cols()
I actually think this is much better than the code you currently use,
* or implement row and column objects.
The problem in your code is that you do not distinguish iterating over rows and columns, and in linear algebra these are different objects. The two solutions I propose do this distinction. My favorite one is the first one (the rows and cols methods).
Cheers,
Gaël
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
for row in A.rows(): for col in row.cols()
Actually I am in favor of adding ``rows`` and ``cols`` attributes to allow such iteration. The only thing is, I say these should be matrices (i.e., 2d). I.e., this should provide a symmetric syntax (for rows and columns) for iteration across submatrices.
But once we have these, it is all the more reason not to object to letting iteration on the matrix yield 1d arrays. As others have observed: why break array behavior for no purpose?
Cheers, Alan Isaac
PS The possibility of adding these attributes is mentioned on the discussion page: URL:http://www.scipy.org/MatrixIndexing as are most of the issues that have been coming up again.
On Sat, Apr 26, 2008 at 02:42:04PM -0400, Alan G Isaac wrote:
As others have observed: why break array behavior for no purpose?
Because the behavior you are relying on is IMHO a coincidence.
Are other people on the mailing list relying on this? I strongly think it is wrong, but maybe I am the lone man.
Gaël
Gael Varoquaux wrote:
to use a syntax similar to dictionnaries:
for row in A.rows(): for col in row.cols()
I actually think this is much better than the code you currently use,
or implement row and column objects.
The problem in your code is that you do not distinguish iterating over rows and columns, and in linear algebra these are different objects. The two solutions I propose do this distinction. My favorite one is the first one (the rows and cols methods).
What would rows() and columns() return? I'd like them to return a vector object. If they return matrices, then you can't index them with a single value, if they return 1-d arrays, then you can't distinguish between rows and columns (which may be OK).
Travis E. Oliphant wrote:
I'm not personally persuaded by the iteration argument, because we can change iteration independently of mapping (__getitem__) access.
or just use:
for row in M.A: ...
How hard is that?
-Chris
Alan G Isaac wrote:
On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:
For me this is wrong. list and tuples are not 2D. Numpy arrays happen to offer this feature, but you should not use it do to multiple dimension indexing.
But there is no proposal that people should index like this. The underlying issue is iteration. Currently::
>>> A matrix([[1, 2], [ 3, 4]]) >>> A = np.mat(x) >>> for row in A: ... for col in row: ... print col ... [[1 2]] [[3 4]]
So are you saying that one should not be able to iterate through to the elements of a matrix?
Iteration can be handled separately with __iter__ method that returns the actual iterator object which may be different. Or, as others have proposed .rows and .cols iterators.
I'm not personally persuaded by the iteration argument, because we can change iteration independently of mapping (__getitem__) access.
-Travis
On Sat, 26 Apr 2008, "Travis E. Oliphant" apparently wrote:
I'm not personally persuaded by the iteration argument, because we can change iteration independently of mapping (__getitem__) access.
Would that not impose another deviation from array behavior? I just do not see the purpose.
I support .row and .cols iterators (returning submatrices), and have included the idea for some time on the discussion page URL:http://www.scipy.org/MatrixIndexing.
Cheers, Alan
I was hoping to get NumPy 1.1 tagged today, but it seems very unlikely at this point. Unfortunately, I haven't followed the matrix discussion as closely as I would like, so I can't tell if there is anything so uncontroversial that it would make sense to change for the 1.1.0 release. If there is something that can be unanimously agreed on within the next 24 hours or so, I would be happy to have it included in 1.1. If not, I would rather see us move on to working on 1.2 and releasing 1.1 ASAP.
If there is unanimous agreement on a minor fix, we will need document the changes on the release notes. So once an agreement is reached, could someone either send me some text to insert into the release notes or update the document themselves: http://projects.scipy.org/scipy/numpy/milestone/1.1.0
This also might be a good opportunity to start a new thread with a subject line like "1.2 matrices". The discussion about matrices in this thread makes it difficult for me to quickly see what still needs to done before the 1.1 release. And there may also be some members of the community who would be very interested in the matrices discussion, but aren't reading this thread.
Alan G Isaac wrote:
On Thu, 24 Apr 2008, Christopher Barker apparently wrote:
I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.
This bothers me.
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
exactly one of the things that really bugged me with MATLAB!
Or are you still having separate classes and not just changing ``__getitem__``? OK, I see that is most likely what you meant.
yes, it would still be a vector class, so a 1 by 2 matrix (or array) would still act like a 2-d object.
The only issue here is if you need a separate "row vector" and "column vector" class -- I think it's a matter of what whether these vectors are special versions of 1-d arrays or 2-d matrices. I kind of like them as special versions of 1-d arrays, as a vector is a 1-d object, but I see that the implementation is a bit cleaner if you just make them 2-d matrices with a different scalar indexing.
-Chris
On Fri, Apr 25, 2008 at 10:23 AM, Christopher Barker Chris.Barker@noaa.gov wrote:
Alan G Isaac wrote:
On Thu, 24 Apr 2008, Christopher Barker apparently wrote:
I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.
This bothers me.
So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.
exactly one of the things that really bugged me with MATLAB!
I thought MATLAB treated a single index as a flat index. So it yields a scalar no matter what (actually a 1x1 matrix no matter what, since it's all matrices in matlab).
--bb
Christopher Barker wrote: [cjop]
Bruce Southey wrote:
Hi, I would like to use the matrix functionality but I, like others, get by without it.
What does it lack that keeps you from using it? That's the key question?
I mainly do matrix-matrix and matrix-vector multiplications like in linear models and predictions. The main reason I have already learnt how to do to what I want using Numpy via mainly Numarray. So at the moment I am not interested in doing something else.
+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.
This is all about being able to index a "vector" with a single index -- that's it, I think. Matlab handles this by special casing matrices that have a single dimension of one. It also has an easy way to spell a literal for column vector, though I suppose:
np.Matrix((1,2,3,4,5)).T
isn't so bad.
-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error.
But HOW can you know the orientation? a 1-d array has no orientation --
I fully agree and which is why there is absolutely no sense in saying row or column vector/array! But vectors always have an orientation in the context of a matrix (as you also indicate below) or 2-d array so in that situation vectors are matrices. Numpy and other options have an implicit orientation because they usually always check that an operation is feasible first so errors are returned if the dimensions are not appropriate. So, as I understand it, Numpy's array do have an implicit orientation because Numpy's 1-d arrays act like 1 by n arrays when used with m by n arrays: import numpy a=numpy.ones((4)) b=numpy.ones((4,3)) numpy.dot(a,b) # gives array([ 4., 4., 4.]) numpy.dot(b,a) # gives error numpy.dot(b.T,a) #gives array([ 4., 4., 4.])
Thus, Numpy's 1-d arrays do have an orientation otherwise all three multiplications should return an error!
which is why the current version always returns a matrix when indexing.
As you point out, that is inconsistent with 1-d arrays having NO orientation. But this is fully consistent with the view that a vector is a special case of a matrix when we commonly refer to an n by 1 or 1 by n vector (as you also indicate below).
0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.
I don't think we can have vector without ColumnVector, though I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.
But didn't you say that '1-d array has no orientation'... :-) This why I am -2 on row and column vector because one should just be the transpose of the other.
By the way, maybe a column vector could be handy for doing array broadcasting, as an cleaner way to do:
x = np.arange(10) y = np.arange(20).reshape((-1,1))
z = FunctionOfXandY(x,y)
Final note:
If no one has an linear algebra examples, then we're wasting out time with even this cheap talk.
-Chris
Bruce
I think the use of the term 'vector' in this thread is becoming a bit confusing.
An M by N matrix is a vector. (I.e., it is an element of a vector space.)
Many people use the terms "row vector" and "column vector" to refer to special matrices. What is special is *not* that they are vectors (since all matrices are) but that they are a single row or a single column. Perhaps better terminology would have been "row matrix" and "column matrix". Using this terminology, we clearly expect a row matrix and a column matrix to have transposes. The question is, can their elements be indexed by scalars. The answer in all the proposals, I believe, is "yes".
I believe we are finding that this answer has some problems.
x is a 2 by 2 matrix x[0] is a 1 by 2 matrix x[0][0] is x[0,0] GOOD!! y = x[0] -> y is a 1 by 2 matrix y[0][0] is a TypeError BAD!!
So I am becoming more and more persuaded that for a matrix ``x``:
x[0] should be a 1d array (always!) x[0,:] should be a matrix
The rule:
use non-scalar indexing to extract submatrices
fwiw, Alan Isaac
On Fri, Apr 25, 2008 at 10:04 AM, Alan G Isaac aisaac@american.edu wrote:
I think the use of the term 'vector' in this thread is becoming a bit confusing.
An M by N matrix is a vector. (I.e., it is an element of a vector space.)
Sure, but the important thing is the multiplication. If it wasn't we could just flatten them all and be done with this discussion. Super bonus, no need for the '*' operator.
Chuck
2008/4/23 Christopher Barker Chris.Barker@noaa.gov:
Aside from the fact that someone needs to write the code -- why don't people like the row/column vector idea? It just feels so natural to me:
I wrote most of the code last week (see the previous thread on the mailing list for a patch). It currently only implements Vector (not RowVector and ColumnVector), but those two would be easy to add.
Cheers Stéfan