[Numpy-discussion] Tensor contraction

[Friedrich Romstedt]

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