[Numpy-discussion] trace does not behave as advertised on arrays of rank > 2

[Harald Hanche-Olsen]

>> print Numeric.trace.__doc__
trace(a,offset=0, axis1=0, axis2=1) returns the sum along diagonals
    (defined by the last two dimenions) of the array.

For arrays of rank 2, trace does what you expect, but for arrays of
larger rank, it appears to simply sum along each of the two given
axes.  A simple experiment follows:

>>> B
array([[[       1,       10],
        [     100,     1000]],
       [[   10000,   100000],
        [ 1000000, 10000000]]])

>>> # What I thought trace(B) would be:
>>> B[0,0,0]+B[1,1,0], B[0,0,1]+B[1,1,1]
(1000001, 10000010)

>>> # But that is not what numpy thinks:
>>> Numeric.trace(B)
array([   10001, 10001000])

>>> # Instead, it must be computing it as follows:
>>> B[0,0,0]+B[1,0,0], B[0,1,1]+B[1,1,1]
(10001, 10001000)

That is, trace(B) is the vector C, given by C[i]=sum(B[j,i,i]: j=0,...).
A bit more experimentation reveals that trace ignores its fourth
argument, consistent with the above result:

>>> Numeric.trace(B,0,0,1)
array([   10001, 10001000])
>>> Numeric.trace(B,0,0,2)
array([   10001, 10001000])

Evidently, trace is going to need a rewrite.  It might perhaps also
benefit from further optional arguments in groups of three, e.g.,

  trace(A, p, 0, 3, q, 1, 2)[k,l,...] = A[i+p,j+q,j,i,k,l,...]

with summing over repeated indices (i, j) ala Einstein.

- Harald