# [Numpy-discussion] Tensor contraction

Friedrich Romstedt friedrichromstedt@gmail....
Mon Jun 14 15:36:58 CDT 2010

2010/6/13 Alan Bromborsky <abrombo@verizon.net>:
> Friedrich Romstedt wrote:
>>> I am writing symbolic tensor package for general relativity.  In making
>>> symbolic tensors concrete
>>> I generate numpy arrays stuffed with sympy functions and symbols.
>>
>> That sound's interesting.

I'm not so sure if I can help you with the subject.  It seems general
relativity is a bit over my head, just a tiiiny bit ;-)

>>> The
>>> operations are tensor product
>>> (numpy.multiply.outer), permutation of indices (swapaxes),  partial and
>>> covariant (both vector operators that
>>> increase array dimensions by one) differentiation, and contraction.
>>
>> I would like to know more precisely what this differentiations do, and
>> how it comes that they add an index to the tensor.

Ok, this I understood now, nevertheless I'm completely lost with an
over-all understanding ...

>>> I think I need to do the contraction last
>>> to make sure everything comes out correctly.  Thus in many cases I would
>>> be performing multiple contractions
>>> on the tensor resulting from all the other operations.

This is up to you, as I see it now.

>> def contract(arr, *contractions):
>>     """*CONTRACTIONS is e.g.:
>>         (0, 1), (2, 3)
>>     meaning two contractions, one of 0 & 1, and one of 2 & 2,
>>     but also:
>>         (0, 1, 2),
>>     is allowed, meaning contract 0 & 1 & 2."""
>>
>>     # First, we check if we can contract using the *contractions* given ...
>>
>>     for contraction in contractions:
>>         # Extract the dimensions used.
>>         dimensions = numpy.asarray(arr.shape)[list(contraction)]
>>
>>         # Check if they are all the same.
>>         dimensionsdiff = dimensions - dimensions[0]
>>         if numpy.abs(dimensionsdiff).sum() != 0:
>>             raise ValueError('Contracted indices must be of same dimension.')
>>
>>     # So now, we can contract.
>>     #
>>     # First, pull the contracted dimensions all to the front ...
>>
>>     # The names of the indices.
>>     names = range(arr.ndim)
>>
>>     # Pull all of the contractions.
>>     names_pulled = []
>>     for contraction in contractions:
>>         names_pulled = names_pulled + list(contraction)
>>         # Remove the pulled index names from the pool:
>>         for used_index in contraction:
>>             # Some more sanity check
>>             if used_index not in names:
>>                 raise ValueError('Each index can only be used in one
>> contraction.')
>>             names.remove(used_index)
>>
>>     # Concatenate the pulled indices and the left-over indices.
>>     names_final = names_pulled + names
>>
>>     # Perform the swap.
>>     arr = arr.transpose(names_final)
>>
>>     # Perform the contractions ...
>>
>>     for contraction in contractions:
>>         # The respective indices are now, since we pulled them, the
>> frontmost indices:
>>         ncontraction = len(contraction)
>>         # The index array:
>>         # shape[0] = shape[1] = ... = shape[ncontraction - 1]
>>         I = numpy.arange(0, arr.shape[0])
>>         # Perform contraction:
>>         index = [I] * ncontraction
>>         arr = arr[tuple(index)].sum(axis=0)
>>
>>     # If I made no mistake, we should be done now.
>>     return arr

Btw, did this work for you?

I like Sebastian's stride tricks, but I think the poor-man's
.transpose() and contraction should also do it.  (and is maybe more

Btw, how are you planning to implement the differentiation?  Hmm ...
for the "normal" partial derivative ... it would be most elegant to
create the \partial_{l} operator as an sympy object, is that possible?
If not, one can think about a wrapper class, yielding the
differentiation upon multiplication with a normal sympy object.  If
one has such operators, one can put them in an array, say p for
"partial" and write:

partial_T = p[numpy.newaxis, numpy.newaxis, numpy.newaxis] * T

Or just write p as an instance of a class P with overloaded __mul__:

def __mul__(self, T):
return self.partials[tuple([numpy.newaxis] * T.ndim)] * T

Then partial differentiation is for all tensors of your
multidimensional world just:

partial_T = p * T

. ? For the Christoffel symbols, it's I think not so easy.  But since
I don't understand them, ...

Friedrich