numpy-discussion

[Numpy-discussion] tensor ops in dim > 3?

[Dr. D. G. <go...@sp...>]

I have a rank three n-dim tensor A.  For two of the axes I want to perform
v^t A v (a quadratic form with scalar output, where v is a vector). The
final output should be a vector.  I also need to compute the derivative of
this with respect to v.  This involves symmetrizing and matrix-vector
multiplication (2 sym(A)v using two axes of A only, which gives a vector)
with the final result being a matrix.

Is there something elegant I can do short of extending numpy with fortrash
routines? Perhaps something like SciPy's 3-d tensor ops but for n-dim?

Thanks in advance for any tips, - d

PS One more dumb question: I just installed the ScientificPython-2.4.1 rpm
on my reincarnated Mandrake linux machine running python2.2.  Do I need to
do something to configure it?  My scripts aren't finding things (e.g.
indexing.py).

Re: [Numpy-discussion] tensor ops in dim > 3?

[Travis O. <oli...@ee...>]

> I have a rank three n-dim tensor A.  For two of the axes I want to perform
> v^t A v (a quadratic form with scalar output, where v is a vector). The
> final output should be a vector.  I also need to compute the derivative of
> this with respect to v.  This involves symmetrizing and matrix-vector
> multiplication (2 sym(A)v using two axes of A only, which gives a vector)
> with the final result being a matrix.
>

I'm not exactly sure what you mean by a rank three n-dim tensor (in
Numeric the rank is the number of dimensions in the array).

But, I think you can accomplish what you desire in two different ways:

1) If A has N dimensions and you want to perform the reduction over axis
I'll label a1 and a2 (that is a1 is the axis for the v^t*A  sum while a2
is the axis for the A*v sum).  Then I think this should work (if a2 > a1).

ex1 = [Numeric.NewAxis]*(N-1)
ex1[a1] = slice(None)
ex2 = [Numeric.NewAxis]*N
ex2[a2] = slice(None)

Nar = Numeric.add.reduce
result = Nar(v[ex1]*Nar(A*v[ex2],axis=a2),axis=a1)

# I think you need a recent version of Numeric for the axis keyword
#   to be defined here.  Otherwise, just pass a2 and a1 as arguments
#   without the keyword.


2) Using dot (requires transposing the matrix) as dot only operates over
certain dimensions --- I would not try this.

-Travis Oliphant

Re: [Numpy-discussion] tensor ops in dim > 3?

[Konrad H. <hi...@cn...>]

"Dr. Dmitry Gokhman" <go...@sp...> writes:

> I have a rank three n-dim tensor A.  For two of the axes I want to perform
> v^t A v (a quadratic form with scalar output, where v is a vector). The
> final output should be a vector.  I also need to compute the derivative of
> this with respect to v.  This involves symmetrizing and matrix-vector
> multiplication (2 sym(A)v using two axes of A only, which gives a vector)
> with the final result being a matrix.

Whenever dealing with somewhat more complex operations of this time, I think
it's best to go back to the basic NumPy functionality rather then figuring
out if there happens to be a function that magically does it. In this case,
assuming that the first axis of A is the one that is not summed over:

   sum( sum(A*v[NewAxis, NewAxis, :], -1) * v[NewAxis, :], -1)

The idea is to align v with one of the dimensions of A, then multiply
elementwise and sum over the common axis. Note that the first (inner)
sum leaves a rank-2 array, so for the second multiplication v gets
extended to rank-2 only.

> PS One more dumb question: I just installed the ScientificPython-2.4.1 rpm
> on my reincarnated Mandrake linux machine running python2.2.  Do I need to
> do something to configure it?  My scripts aren't finding things (e.g.
> indexing.py).

If you took the binary RPM from my site, they might not work correctly
with Mandrake, as they were made for RedHat. The source RPM should
work with all RPM-based Linux distributions. There is nothing that needs
configuring with an RPM.

Also note that Scientific is a package, so the correct way to import the
indexing module is

   import Scientific.indexing

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hi...@cn...
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------