[Numpy-discussion] .T Transpose shortcut for arrays again

[Sasha]

On 7/6/06, Bill Baxter <wbaxter at gmail.com> wrote:
> ...
> Yep, like Tim said.  The usage is say a N sets of basis vectors.  Each set
> of basis vectors is a matrix.

This brings up a feature that I really miss from numpy: an ability to do

array([f(x) for x in a])

without python overhead.  APL-like languages have a notion of "adverb"
- a higher level operator that maps a function to a function. Numpy
has some adverbs implemented as attributes to ufuncs: for example
add.reduce is the same as +/ in K and add.accumulate is the same as +\
('/' and '\' are 'over' and 'scan' adverbs in K). However, there is no
way to do f/ or f\ where f is an arbitrary dyadic function.

The equivalent of array([f(x) for x in a]) is spelled f'(a) in K (' is
an adverb 'each'). The transpose operator (+) is swaps the first two
axes, so in order to apply to the array of matrices, one would have to
do +:'a (: in +: disambiguates + as a unary operator).

I don't know of a good way to introduce adverbs in numpy, nor can I
think of a good way to do list comprehensions, but array friendly
versions of map, filter and reduce may be a good addition.  These
higher order functions may take an optional axes argument to deal with
the higher rank arrays and may be optimized to recognize ufuncs so
that map(f, a) could call f(a) and reduce(f, a) could do f.reduce(a)
when f is a ufunc.

[snip]
> Either way swapaxes(-2,-1) is likely more likely to be what you want than
> .transpose().
>

Agree, but swapaxes(0, 1) is a close runner-up which is also known as
zip in python.

> Well, I would be really happy for .T to return an (N,1) column vector if
> handed an (N,) 1-d array.  But I'm pretty sure that would raise more furuor
> among the readers of the list than leaving it 1-d.
>

Would you be even happier if .T would return a matrix? I hope not
because my .M objection will apply.  Maybe we can compromize by
implementing a.T so that it raises ValueError unless rank(a) == 2 or
at least unless rank(a) <= 2?

> I have serious reservations about a function called t().  x,y,z, and t are
> probably all in the top 10 variable names in scientific computing.
>

What about T()?

>
> > K (an APL-like language) overloads
> > unary '+' to do swapaxes(0,1) for rank>=2 and nothing for lower rank.
>
>
> Hmm.  That's kind of interesting, it seems like an abuse of notation to me.
> And precedence might be an issue too.  The precedence of unary + isn't as
> high as attribute access.

It is high enough AFAICT - higher than any binary operator.

>  Anyway, as far as the meaning of + in K, I'm
> guessing K's arrays are in Fortran order, so (0,1) axes vary the fastest.

No, K has 1d arrays only, but they can be nested.  Matrices are arrays
of arrays and tensors are arrays of arrays of arrays ..., but you are
right (0,1)  swap is faster than (-2,-1) swap and this motivated the
choice for the primitive.

> I couldn't find any documentation for the K language from a quick search,
> though.

Kx Systems, the company behind K has replaced K with Q and pulled old
manuals from the web.  Q is close enough to K: see
http://kx.com/q/d/k.txt for a terse summary.

[snip]
> > Why would anyone do that if b was a matrix?
> Maybe because, like you, they think "that a.T is fairly cryptic".
>
If they are like me, they will not use numpy.matrix to begin with :-).

>
> > > But probably a better solution
> > > would be to have matrix versions of these in the library as an optional
> > > module to import so people could, say, import them as M and use
> M.ones(2,2).
> > >
> >
> > This is the solution used by ma, which is another argument for it.
>
> Yeh, I'm starting to think that's better than slapping an M attribute on
> arrays, too.  Is it hard to write a module like that?
>

Writing matrixutils with

def zeros(shape, dtype=float):
      return asmatrix(zeros(shape, dtype))

is trivial, but matrixutils.zeros will have two python function calls
overhead.  This may be a case for making zeros a class method of
ndarray that can be written in a way that will make inherited
matrix.zeros do the right thing with no overhead.

[snip]
> * +A implies addition.
No, it does not.  Unary '+' is a noop.  Does * imply multiplication or
** imply pow in f(*args, **kwds) to you?

> The general rule with operator overloading is that
> the overload should have the same general meaning as the original operator.
Unary '+' has no preset meaning in plain python. It can be interpreted
as transpose if you think of scalars as 1x1 matrices.

> So overloading * for matrix multiplication makes sense.

It depends on what you consider part of "general meaning".  If the
commutativity property is part of it then overloading * for matrix
multiplication doesn't make sense. If the "general meaning" of unary +
includes x = +x invariant, then you are right, but I am willing to
relax that to x = ++x invariant when x is a non-symmetric matrix.

>  ... New users looking at something like A + +B are pretty
> certain to be confused because they think they know what + means, but
> they're wrong.

In my experience new users don't realize that unary + is defined for
arrays.  Use of unary + with non-literal numbers is exotic enough that
new users seeing "something like A + +B" will not assume that they
know what it means.

[snip]
> * +A has different precedence than the usual transpose operator.  (But I
> can't think of a case where that would make a difference now.)
>
Maybe you can't because it doesn't? :-)

> I would be willing to accept a .T that just threw an exception if ndim were
> > 2.

Aha!  Let's start with an error unless ndim != 2.  It is always easier
to add good features than to remove bad ones.