Re: Best way to do contractions (arbitrary Tables with a
- To: mathgroup at smc.vnet.net
- Subject: [mg105973] Re: [mg105937] Best way to do contractions (arbitrary Tables with a
- From: Leonid Shifrin <lshifr at gmail.com>
- Date: Mon, 28 Dec 2009 04:58:36 -0500 (EST)
- References: <200912270724.CAA21087@smc.vnet.net>
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
On Sat, Dec 26, 2009 at 11:24 PM, Erik Max Francis <max at 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 && max at 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.
>
>
- References:
- Best way to do contractions (arbitrary Tables with a Sum)?
- From: Erik Max Francis <max@alcyone.com>
- Best way to do contractions (arbitrary Tables with a Sum)?