[Numpy-discussion] Tensor contraction
David Goldsmith
d.l.goldsmith at gmail.com
Sun Jun 13 15:59:51 EDT 2010
Is this not what
core.numeric.tensordot<http://docs.scipy.org/numpy/docs/numpy.core.numeric.tensordot/>does?
DG
On Sun, Jun 13, 2010 at 12:37 PM, Friedrich Romstedt <
friedrichromstedt at gmail.com> 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
>
--
Mathematician: noun, someone who disavows certainty when their uncertainty
set is non-empty, even if that set has measure zero.
Hope: noun, that delusive spirit which escaped Pandora's jar and, with her
lies, prevents mankind from committing a general suicide. (As interpreted
by Robert Graves)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20100613/d153c4d2/attachment.html>
More information about the NumPy-Discussion
mailing list