Fast construction and reshaping of large tensors

I'm trying to construct a large matrix which is derived from some higher rank tensor (the rank of interest to me changes case by case, so it needs to be a general method). Currently, the process of building the tensor is taking a huge amount of time. Schematically, the thing I'm trying to build (I think) looks like this.

latsize = 20;
d = 3;
t = Table[Symbol["t" <> ToString[i]], {i, d}];
s = Table[Symbol["s" <> ToString[i]], {i, d}];
dotprod = (t - s).(t - s);

Print["Time to do the whole thing = ", First[AbsoluteTiming[
ulimit = ConstantArray[latsize, 2 d];
llimit = ConstantArray[1, 2 d];
tablelimits = Transpose[{Join[t, s], llimit, ulimit}];
Print["Time to construct tensor = ", First[AbsoluteTiming[
tensor = Table[dotprod, ##] & @@ tablelimits;
]]];
Print["Time to construct matrix = ", First[AbsoluteTiming[
matrix = ArrayReshape[tensor, {latsize^d, latsize^d}];
]]];
]]];
Print["Total number of elements = ", Length@Flatten@tensor];

The process quickly becomes lengthy as d and latsize are increased. Any suggestions for how to speed this up?

• Have you looked at the Tensor functions: reference.wolfram.com/language/guide/SymbolicTensors.html ? I'm not sure what general tensor you're trying to build but multiple calls to KroneckerProduct ,TensorProduct and Outer might speed things up. – Histograms May 21 '15 at 10:23
• Hi Histograms, I'm only concerned with building a large array with numerical entries. Right now the slowest part by far is the line: tensor = Table[dotprod, ##] & @@ tablelimits; So in the absence of a better way of doing the full calculation, really my question boils down to finding a faster way of constructing this Table with numerical entries. – user12876 May 21 '15 at 12:29
• Could the post "What is the most efficient way to construct large block matrices in Mathematica?" be useful? So far I'm still struggling with it. – user12876 May 21 '15 at 13:14

The following approach is faster, but the sizes of your output matrices will quickly overwhelm Mathematica's memory. At any rate:

constructMatrix[l_,d_]:=With[{one=Table[(t-s)^2,{t,l},{s,l}]},
ArrayReshape[
Transpose[
Outer[Plus, Sequence @@ Table[one, d]],
Flatten @ Transpose[{Range[d], Range[d]+d}]
],
{l^d,l^d}
]
]

Compared to your solution:

latsize=20;
d=3;
t=Table[Symbol["t"<>ToString[i]],{i,d}];
s=Table[Symbol["s"<>ToString[i]],{i,d}];
dotprod=(t-s).(t-s);

Print["Time to do the whole thing = ",First[AbsoluteTiming[ulimit=ConstantArray[latsize,2 d];
llimit=ConstantArray[1,2 d];
tablelimits=Transpose[{Join[t,s],llimit,ulimit}];
Print["Time to construct tensor = ",First[AbsoluteTiming[tensor=Table[dotprod,##]&@@tablelimits;]]];
Print["Time to construct matrix = ",First[AbsoluteTiming[matrix=ArrayReshape[tensor,{latsize^d,latsize^d}];]]];]]];
Print["Total number of elements = ",Length@Flatten@tensor];

Time to construct tensor = 6.27127

Time to construct matrix = 0.336778

Time to do the whole thing = 6.6083

Total number of elements = 64000000

m = constructMatrix[20, 3]; //AbsoluteTiming

m === matrix

{0.690541, Null}

True

So, about an order of magnitude faster. For larger sizes:

r1 = constructMatrix[30, 3]; //AbsoluteTiming
r2 = constructMatrix[10, 4]; //AbsoluteTiming

{11.7492, Null}

{61.3809, Null}

Going much larger will either take too long, or kill my kernel due to a memory exception.