[Numpy-discussion] Syntax Improvement for Array Transpose

[Andras Deak]

Dear Ilhan,

Thanks for writing these up.
I feel that from a usability standpoint most people would support #3
(.H/.mH), especially considering Marten's very good argument about @.
Having to wrap your transposed matrices in function calls half defeats
the purpose of being able to write stacked matrix operations elegantly
within the ndarray class. The question is of course whether it's
feasible from a project management/API design stand point (just to
state the obvious).
Regarding #1 (1d transpose): I just want to make it clear as someone
who switched from MATLAB to python (and couldn't be happier) that we
should treat MATLAB's behaviour as more of a cautionary tale rather
than design ethos. I paused for exactly 5 seconds the first time I ran
into the no-op of 1d transposes, and then I thought "yes, this makes
sense", and that was it. To put it differently, I think it's more
about MATLAB injecting false assumptions into users than about numpy
behaving surprisingly. (On a side note, MATLAB's quirks are one of the
reasons that the Spyder IDE, designed to be a MATLAB replacement, has
very weird quirks that regularly trip up python users.)
Regards,

András

On Wed, Jun 26, 2019 at 9:04 AM Ilhan Polat <ilhanpolat at gmail.com> wrote:
>
> Maybe a bit of a grouping would help, because I am also losing track here. Let's see if I could manage to get something sensible because, just like Marten mentioned, I am confusing myself even when I am thinking about this
>
> 1- Transpose operation on 1D arrays:
>     This is a well-known confusion point for anyone that arrives at NumPy usage from, say matlab background or any linear algebra based user. Andras mentioned already that this is a subset of NumPy users so we have to be careful about the user assumptions. 1D arrays are computational constructs and mathematically they don't exist and this is the basis that matlab enforced since day 1. Any numerical object is an at least 2D array including scalars hence transposition flips the dimensions even for a col vector or row vector. That doesn't mean we cannot change it or we need to follow matlab but this is kind of what anybody kinda sorta wouda expect. For some historical reason, on numpy side transposition on 1D arrays did nothing since they have single dimensions. Hence you have to create a 2D vector for transpose from the get go to match the linear algebra intuition. Points that has been discussed so far are about whether we should go further and even intercept this behavior such that 1D transpose gives errors or warnings as opposed to the current behavior of silent no-op. as far as I can tell, we have a consensus that this behavior is here to stay for the foreseeable future.
>
> 2- Using transpose to reshape the (complex) array or flip its dimensions
>     This is a usage that has been mentioned above that I don't know much about. I usually go the "reshape() et al." way for this but apparently folks use it to flip dimensions and they don't want the automatically conjugation which is exactly the opposite of a linear algebra oriented user is used to have as an adjoint operator. Therefore points that have been discussed about are whether to inject conjugation into .T behavior of complex arrays or not. If not can we have an extra .H or something that specifically does .conj().T together (or .T.conj() order doesn't matter). The main feel (that I got so far) is that we shouldn't touch the current way and hopefully bring in another attribute.
>
> 3- Having a shorthand notation such as .H or .mH etc.
>     If the previous assertion is true then the issue becomes what should be the new name of the attribute and how can it have the nice properties of a transpose such as returning a view etc. However this has been proposed and rejected before e.g., GH-8882 and GH-13797. There is a catch here though, because if the alternative is .conj().T then it doesn't matter whether it copies or not because .conj().T doesn't return a view either and therefore the user receives a new array anyways. Therefore no benefits lost. Since the idea is to have a shorthand notation, it seems to me that this point is artificial in that sense and not necessarily a valid argument for rejection. But from the reluctance of Ralf I feel like there is a historical wear-out on this subject.
>
> 4- transpose of 3+D arrays
>     I think we missed the bus on this one for changing the default behavior now and there are glimpses of confirmation of this above in the previous mails. I would suggest discussing this separately.
>
> So if you are not already worn out and not feeling sour about it, I would like to propose the discussion of item 3 opened once again. Because the need is real and we don't need to get choked on the implementation details right away.
>
> Disclaimer: I do applied math so I have a natural bias towards the linalg-y way of doing things. And sorry about that if I did that above, sometimes typing quickly loses the intention.
>
>
> Best,
> ilhan
>
>
> On Wed, Jun 26, 2019 at 4:39 AM Ralf Gommers <ralf.gommers at gmail.com> wrote:
>>
>>
>>
>> On Wed, Jun 26, 2019 at 3:56 AM Marten van Kerkwijk <m.h.vankerkwijk at gmail.com> wrote:
>>>
>>> Hi Ralf,
>>>
>>> On Tue, Jun 25, 2019 at 6:31 PM Ralf Gommers <ralf.gommers at gmail.com> wrote:
>>>>
>>>>
>>>>
>>>> On Tue, Jun 25, 2019 at 11:02 PM Marten van Kerkwijk <m.h.vankerkwijk at gmail.com> wrote:
>>>>>
>>>>>
>>>>> For the names, my suggestion of lower-casing the M in the initial one, i.e., `.mT` and `.mH`, so far seemed most supported (and I think we should discuss *assuming* those would eventually involve not copying data; let's not worry about implementation details).
>>>>
>>>>
>>>> For the record, this is not an implementation detail. It was the consensus before that `H` is a bad idea unless it returns a view just like `T`: https://github.com/numpy/numpy/issues/8882
>>>
>>>
>>> Is there more than an issue in which Nathaniel rejecting it mentioning some previous consensus?
>>
>>
>> Yes, this has been discussed in lots of detail before, also on this list (as Nathaniel mentioned in the issue). I spent 10 minutes to try and find it but that wasn't enough. I do think it's not necessarily my responsibility though to dig up all the history here - that should be on the proposers of a new feature ....
>>
>>> I was part of the discussion of the complex conjugate dtype, but do not recall any consensus beyond a "wish to keep properties simple". Certainly the "property does not do any calculation" rule seems arbitrary; the only strict rule I would apply myself is that the computation should not be able to fail (computationally, out-of-memory does not count; that's like large integer overflow). So,  I'd definitely agree with you if we were discussion a property `.I` for matrix inverse (and indeed have said so in related issues). But for .H, not so much. Certainly whoever wrote np.Matrix didn't seem to feel bound by it.
>>>
>>> Note that for *matrix* transpose (as opposed to general axis reordering with .tranpose()), I see far less use for what is returned being a writable view. Indeed, for conjugate transpose, I would rather never allow writing back even if it we had the conjugate dtype since one would surely get it wrong (likely, `.conj()` would also return a read-only view, at least by default; perhaps one should even go as far as only allowing `a.view(conjugate-dtype)` as they way to get a writable view).
>>>
>>>>
>>>> So, specific items to confirm:
>>>>
>>>> 1) Is this a worthy addition? (certainly, their existence would reduce confusion about `.T`... so far, my sense is tentative yes)
>>>>
>>>> 2) Are `.mT` and `.mH` indeed the consensus? [1]
>>>
>>>
>>> > I think `H` would be good to revisit *if* it can be made to return a view. I think a tweak on `T` for >2-D input does not meet the bar for inclusion.
>>>
>>> Well, I guess it is obvious I disagree: I think this more than meets the bar for inclusion. To me, this certainly is a much bigger deal that something like oindex or vindex (which I do like).
>>
>>
>> Honestly, I don't really want to be arguing against this (or even be forced to spend time following along here). My main problem with this proposal right now is that we've had this discussion multiple times, and it was rejected with solid arguments after taking up a lot of time. Restarting that discussion from scratch without considering the history feels wrong. It's like a democracy voting on becoming a dictatorship repeatedly: you can have a "no" vote several times, but if you rerun the vote often enough at some point you'll get a "yes", and then it's a done deal.
>>
>> I think this requires a serious write-up, as either a NEP or a GitHub issue with a good set of cross-links and addressing all previous arguments.
>>
>>>
>>> Indeed, it would seem to me that if a visually more convenient way to do (stacks of) matrix multiplication for numpy is good enough to warrant changing the python syntax, then surely having a visually more convenient standard way to do matrix transpose should not be considered off-limits for ndarray; how often do you see a series matrix manipulations that does not involve both multiplication and transpose?
>>>
>>> It certainly doesn't seem to me much of an argument that someone previously decided to use .T for a shortcut for the computer scientist idea of transpose to not allow the mathematical/physical-scientist one - one I would argue is guaranteed to be used much more.
>>>
>>> The latter of course just repeats what many others have written above, but since given that you call it a "tweak", perhaps it is worth backing up. For astropy, a quick grep gives:
>>>
>>> - 28 uses of the matrix_transpose function I wrote because numpy doesn't have even a simple function for that and the people who wrote the original code used the Matrix class which had the proper .T (but doesn't extend to multiple dimensions; we might still be using it otherwise).
>>
>>
>> A utility function in scipy.linalg would be a more low-API-impact approach to addressing this.
>>
>>>
>>> - 11 uses of .T,  all of which seem to be on 2-D arrays and are certainly used as if they were matrix transpose (most are for fitting). Certainly, all of these are bugs lying in waiting if the arrays ever get to be >2-D.
>>
>>
>> Most linalg is 2-D, that's why numpy.matrix and scipy.sparse matrices are 2-D only. If it's a real worry for those 11 cases, you could just add some comments or tests that prevent introducing bugs.
>>
>> More importantly, your assumption that >2-D arrays are "stacks of matrices" and that other usage is for "computer scientists" is arguably incorrect. There are many use cases for 3-D and higher-dimensional arrays that are not just "vectorized matrix math". As a physicist, I've done lots of work with 3-D and 4-D grids for everything from quantum physics to engineering problems in semiconductor equipment. NumPy is great for that, and I've never needed >=3-D linalg for any of it (and transposing is useful). So please don't claim the physicial-scientist view for this:)
>>
>> Cheers,
>> Ralf
>>
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at python.org
>> https://mail.python.org/mailman/listinfo/numpy-discussion
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion