2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Histograms May 21 '15 at 10:23
  • $\begingroup$ 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. $\endgroup$ – user12876 May 21 '15 at 12:29
  • $\begingroup$ 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. $\endgroup$ – user12876 May 21 '15 at 13:14
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.