[Numpy-discussion] Tensor contraction

[Alan Bromborsky]

Friedrich Romstedt wrote:
> 2010/6/13 Alan Bromborsky <abrombo at verizon.net>:
>   
>> 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.
>
>   
>> 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.
>
>   
>> 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.
>>     
>
> Hm, ok, so I guess I shall give my 1 cent now.
>
> Ok.
>
>     # First attempt (FYI, failed):
>
>     # The general procedure is, to extract a multi-dimensional diagonal array.
>     # The sum \sum_{ij = 0}^{M} \sum_{kl = 0}^{N} is actually the sum over a
>     # 2D array with indices I \equiv i \equiv j and K \equiv k \equiv
> l.  Meaning:
>     # \sum_{(I, K) = (0, 0)}^{(M, N)}.
>     # Thus, if we extract the indices with 2D arrays [[0], [1], ...,
> [N - 1]] for I and
>     # [[0, 1, ..., M - 1]] on the other side for K, then numpy's broadcasting
>     # mechanism will broadcast them to the same shape, yielding (N, M) arrays.
>     # Then finally we sum over this X last dimensions when there were X
>     # contractions, and we're done.
>
>     # Hmmm, when looking closer at the problem, it seems that this isn't
>     # adequate.  Because we would have to insert open slices, but cannot
>     # create them outside of the [] operator ...
>
>     # So now follows second attemt:
>
> 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
>
> Ok, it didn't get much shorter than Pauli's solution, so you decide ...
>
>   
>> One question to
>> ask would be considering that I am stuffing
>> the arrays with symbolic objects and all the operations on the objects
>> would be done using the sympy modules,
>> would using numpy operations to perform the contractions really save any
>> time over just doing the contraction in
>> python code with a numpy array.
>>     
>
> I don't know anything about sympy.  I think there's some typo around:
> I guess you mean creating some /sympy/ array and doing the operations
> using that instead of using a numpy array having sympy dtype=object
> content?
>
> Friedrich
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>   
No I stuff a numpy array with sympy objects.  Except for contraction and 
differentiation numpy does everything I need to do with the tensors. 
Namely -

addition
subtraction
multiplication by a scalar (sympy object)
tensor product (multiply.outer)
index permutation (swapaxes)

Plus -

contraction
differentiation

There are two types of vector derivatives (think gradient) partial and 
covariant.  In index notation the partial derivative of an array T_{ijk} 
is the array

\partial_{l}T_{ijk} = \frac{\partial T_{ijk}}{\partial x^{l}}

The covariant derivative is

D_{l}T_{ijk} = \partial_{l}T_{ijk} 
+T_{uvw}\Gamma_{i}^{u}_{l}\Gamma_{j}^{v}_{l}\Gamma_{k}^{w}_{l}

Where the \Gamma_{i}^{j}_{k} are the Christoffle symbols with 
appropriately raised or lowered indices.