Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Best way to do contractions (arbitrary Tables with a Sum)?

38 views
Skip to first unread message

Erik Max Francis

unread,
Dec 27, 2009, 2:28:03 AM12/27/09
to
I'm trying to do arbitrary contractions with tensors, which basically
amounts to taking an (arbitrarily) large multi-dimensional array,
iterating over the uncontracted indices, and then summing over the two
(and only two) indices to be contracted. If I were dealing with a
specific case, I'd use Table with Sum:

Table[
Sum[
a[[i1]][[i2]]...[[j]]...[[j]]...[[im]]],
{j, n}],
{i1, n}, {i2, n}, ... {im, n}]

That is, iterating over the indices i1, i2, through im (all taking on
values 1 through n) and summing over two of the indices (as j). I'm
trying to figure out the most elegant way to do this in Mathematica and
I'm only coming up with ugly solutions which are basically arbitrary
reimplementations of Table-like functionality.

I figure there's probably some more elegant way to approach this.
Anyone have any ideas?

--
Erik Max Francis && m...@alcyone.com && http://www.alcyone.com/max/
San Jose, CA, USA && 37 18 N 121 57 W && AIM/Y!M/Skype erikmaxfrancis
Without love, benevolence becomes egotism.
-- Dr. Martin Luther King, Jr.

Andrzej Kozlowski

unread,
Dec 28, 2009, 4:56:24 AM12/28/09
to
Read the documentation for Transpose and Inner. Essentially you can
perform arbitrary contractions by first transposing the pair of indices
you wish to contract into the first and last position, then applying
Inner and then Transposing again back to the original position.

Andrzej Kozlowski

Norbert P.

unread,
Dec 28, 2009, 4:57:07 AM12/28/09
to
On Dec 26, 11:28 pm, Erik Max Francis <m...@alcyone.com> wrote:
> I'm trying to do arbitrary contractions with tensors, which basically
> amounts to taking an (arbitrarily) large multi-dimensional array,
> iterating over the uncontracted indices, and then summing over the two
> (and only two) indices to be contracted. If I were dealing with a
> specific case, I'd use Table with Sum:
>
> Table[
> Sum[
> a[[i1]][[i2]]...[[j]]...[[j]]...[[im]]],
> {j, n}],
> {i1, n}, {i2, n}, ... {im, n}]
>
> That is, iterating over the indices i1, i2, through im (all taking on
> values 1 through n) and summing over two of the indices (as j). I'm
> trying to figure out the most elegant way to do this in Mathematica and
> I'm only coming up with ugly solutions which are basically arbitrary
> reimplementations of Table-like functionality.
>
> I figure there's probably some more elegant way to approach this.
> Anyone have any ideas?
>
> --
> Erik Max Francis && m...@alcyone.com &&http://www.alcyone.com/max/

> San Jose, CA, USA && 37 18 N 121 57 W && AIM/Y!M/Skype erikmaxfrancis
> Without love, benevolence becomes egotism.
> -- Dr. Martin Luther King, Jr.

Hi,

not sure if this is the best way, but you don't need to use Table,
use

Sum[a[[All, ..., All, j, All, ..., All, j, All, ...]], {j,n}]

But it is still a Table-like functionality...

Norbert

David Park

unread,
Dec 28, 2009, 4:58:05 AM12/28/09
to
The Tensorial package has facilities for converting indexed tensor
expressions to Mathematica arrays. In any case, the following are some tips
that we list in the "Tensors and Mathematica Arrays" tutorial.

The Prime Rule for Products of 'Tensor' Arrays in Mathematica:
S.T dots the lowest level of S with the highest level of T,
or equivalently
S.T dots the last level of S with the first level of T.

The Mathematica Transpose[T,{n1,n2,n3,...}] moves levels {1,2,3,...} to
levels {n1,n2,n3,...}. We will always want to move the contracted level to
the first or last level when doing Dot products and to the first two levels
when doing single array contractions.

If R, S, T,... are Mathematica tensor arrays, then their direct product is
given by Outer[Times,R,S,T,...]. This will produce a single Mathematica
array. The levels are in the same order as the levels in the successive
arrays.
When expanding tensors to be put in Outer it is best to keep the indicies in
strictly ascending sort order with the slots, or at least within each
tensor.

The basic Mathematica command for contraction of the top two levels in a
single array T is Tr[T,Plus,2]. We will have to use Transpose on T to put
the contraction slots in the first two levels. We will have to repeat the
operation if we want to do multiple contractions.

So essentially you have to Transpose your array to get the contracted
indices into the first two levels of the array and then use Tr.


David Park
djm...@comcast.net
http://home.comcast.net/~djmpark/


From: Erik Max Francis [mailto:m...@alcyone.com]

I'm trying to do arbitrary contractions with tensors, which basically
amounts to taking an (arbitrarily) large multi-dimensional array,
iterating over the uncontracted indices, and then summing over the two
(and only two) indices to be contracted. If I were dealing with a
specific case, I'd use Table with Sum:

Table[
Sum[
a[[i1]][[i2]]...[[j]]...[[j]]...[[im]]],
{j, n}],
{i1, n}, {i2, n}, ... {im, n}]

That is, iterating over the indices i1, i2, through im (all taking on
values 1 through n) and summing over two of the indices (as j). I'm
trying to figure out the most elegant way to do this in Mathematica and
I'm only coming up with ugly solutions which are basically arbitrary
reimplementations of Table-like functionality.

I figure there's probably some more elegant way to approach this.
Anyone have any ideas?

--
Erik Max Francis && m...@alcyone.com && http://www.alcyone.com/max/

Leonid Shifrin

unread,
Dec 28, 2009, 4:59:11 AM12/28/09
to
Hi Erik,

The following function will presumably do what you want and should be
reasonably efficient:


Clear[contract];
contract[tensor_?ArrayQ, firstIndex_Integer, secondIndex_Integer] :=
Module[{dims = Dimensions[tensor], levs, mpos, npos, len},
len = Length[dims];
levs = Range[len];
Map[Tr,
Transpose[tensor,
ReplacePart[
levs - UnitStep[levs - firstIndex] -
UnitStep[levs - secondIndex], {firstIndex -> len - 1,
secondIndex -> len}]], {len - 2}] /;
0 < firstIndex <= len && 0 < secondIndex <= len &&
firstIndex != secondIndex &&
dims[[firstIndex]] == dims[[secondIndex]]];

The main idea is this: the built-in Transpose has an optional second
argument which is a list and allows one to reorder levels in your structure
(tensor). This is equivalent to reordering the indices, so for example by
calling

Transpose[tens,{1,3,2,4}],

you go from T<i,j,k,l> to T<i,k,j,l>.

Then, what is done is the following: make you contraction indices the two
last ones by "removing" them from the tensor index list and adding to the
end, while shifting other indices without any further reordering. Using
Transpose in the above manner with a properly constructed permutation of
levels allows us to implement this strategy. Then Map the Trace (Tr)
function on the appropriate level (which is the number of tensor indices
minus two) of resulting structure (transposed tensor). What you get is the
contracted tensor.

I also inserted the error-checking conditions so that contraction will only
happen when the dimensions of the contracted indices are the same (to
clarify the code a bit - I used a rather rarely used but actually very
useful construct when you share local variables between the body and the
condition, by placing condition inside the body. This is helpful when the
conditions are most easily expressed by using some results which are only
available after some part of the code in the body of the function has been
executed. The advantage is that, if the condition gets violated, the
rule-based dispatch goes on trying next available rules as if no execution
in the body took place, which means that you can reuse some parts of the
code in the body of your definitions in rule-based dispatch mechanism).

A simple test:

In[1]:=
(test = Array[a[#]*b[#2]*c[#3]*d[#4]&,{3,2,3,2}])//Short[#,10]&

Out[1]//Short= {{{{a[1] b[1] c[1] d[1],a[1] b[1] c[1] d[2]},{a[1] b[1] c[2]
d[1],a[1] b[1] c[2] d[2]},{a[1] b[1] c[3] d[1],a[1] b[1] c[3] d[2]}},{{a[1]
b[2] c[1] d[1],a[1] b[2] c[1] d[2]},{a[1] b[2] c[2] d[1],a[1] b[2] c[2]
d[2]},{a[1] b[2] c[3] d[1],a[1] b[2] c[3] d[2]}}},<<1>>,{{{a[3] b[1] c[1]
d[1],a[3] b[1] c[1] d[2]},{a[3] b[1] c[2] d[1],a[3] b[1] c[2] d[2]},{a[3]
b[1] c[3] d[1],a[3] b[1] c[3] d[2]}},{{a[3] b[2] c[1] d[1],a[3] b[2] c[1]
d[2]},{<<1>>,<<1>>},{<<1>>,<<1>>}}}}


In[2]:=
Dimensions[test]

Out[2]= {3,2,3,2}

In[3]:=
contract[test,2,4]//Simplify

Out[3]= {{a[1] c[1] (b[1] d[1]+b[2] d[2]),a[1] c[2] (b[1] d[1]+b[2]
d[2]),a[1] c[3] (b[1] d[1]+b[2] d[2])},{a[2] c[1] (b[1] d[1]+b[2] d[2]),a[2]
c[2] (b[1] d[1]+b[2] d[2]),a[2] c[3] (b[1] d[1]+b[2] d[2])},{a[3] c[1] (b[1]
d[1]+b[2] d[2]),a[3] c[2] (b[1] d[1]+b[2] d[2]),a[3] c[3] (b[1] d[1]+b[2]
d[2])}}


In[4]:=
contract[test,1,3]//Simplify

Out[4]=
{{b[1] (a[1] c[1]+a[2] c[2]+a[3] c[3]) d[1],b[1] (a[1] c[1]+a[2] c[2]+a[3]
c[3]) d[2]},{b[2] (a[1] c[1]+a[2] c[2]+a[3] c[3]) d[1],b[2] (a[1] c[1]+a[2]
c[2]+a[3] c[3]) d[2]}}


I expect <contract> to be reasonably fast since it avoids explcit looping
and uses efficient operaqtions such as Transpose, Map and Tr, on an entire
structure at once. I did not do any timing comparisons or benchmarks
however.

Hope this helps.

Regards,
Leonid

Erik Max Francis

unread,
Dec 29, 2009, 1:22:36 AM12/29/09
to
Leonid Shifrin wrote:
> I expect <contract> to be reasonably fast since it avoids explcit looping
> and uses efficient operaqtions such as Transpose, Map and Tr, on an entire
> structure at once. I did not do any timing comparisons or benchmarks
> however.

Thanks to Leonid and everyone else who responded. I was indeed unaware
of the extra functionality available in Transpose and Tr, which is
exactly what I needed. (I figured there was something lurking
somewhere, but wasn't sure where.) I use the approach of transposing
the two relevant indices to the beginning of the array, then using
Tr[..., Plus, 2]. Works just fine.

For raising and lowering indices, I'll just do a tensor (outer) product
with the metric (or its inverse), do the contraction trick, and then
transpose the results back to the original place. For covariant
derivatives, I'll, uh, figure something out :-). Then I should be ready
do add some documentation and a few examples and post it to see what
people think of my first shot at a package.

Weird way to spend a vacation, I know ... :-).

--
Erik Max Francis && m...@alcyone.com && http://www.alcyone.com/max/
San Jose, CA, USA && 37 18 N 121 57 W && AIM/Y!M/Skype erikmaxfrancis

I have not yet begun to right!
-- John Paul Jones

0 new messages