Fast construction and reshaping of large tensors

Question

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?

Answer 1

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.