# Building very large matrices - strange scaling behaviour

## Problem explanation

I work with a symmetric matrix $M$ that consists of four single matrices. I calculate three of them such that $M$ = $\left(\begin{matrix}a & b \\ b^T & c\end{matrix}\right)$ and transpose $b$. The calculation of each matrix element is constant(!) in time, the Matrix scales according to $M\propto N^6$ where $N$ denotes an input parameter to the function that generates the matrix. The dimensions of the matrices are $a$: $N\times N$, $b$: $N \times N^3$, $c:$ $N^3 \times N^3$.

## Expected as well as received scaling behaviour

My overall expectation is a scaling in time according to $N^6$ as this is the number of matrix elements that have to be calculated. I measured (expected):

• $N=13:$ 18.06 sec
• $N=14:$ 29.13 sec (28.17 sec)
• $N=20:$ 936.93 sec (239.46 sec)

As you can see for increasing $N$ the time noticeably deviates from the expectation. For small values of $N$ the approximation works quite well as can be seen with

exp[baseTime_, i_, n_] := baseTime*(i/6)^n;
baseTime = AbsoluteTiming[BuildMatrix][];
Table[{exp[baseTime, i, 2], exp[baseTime, i, 6], AbsoluteTiming[BuildMatrix[i]][]}, {i, 7, 12}]

(* Out: {{0.590009, 1.09306, 1.09899}, {0.770624, 2.43555, 2.45416}, {0.975321, 4.93756, 4.989}, {1.2041, 9.2909, 9.27697}, {1.45696, 16.4594,16.3822}, {1.7339, 27.7425, 27.5057}}*)


and a quite perfect scaling behaviour according to a $\propto N^6$ expectation (compared to a wrongly assumed $\propto N^2$ behaviour here).

## Code being used

The code I use is the following

t=1.;U=4.;
BuildMatrix[LatticeSize_] :=
Module[{m, n, o, i, j, k, matrix, aMatrix, bMatrix, bMatrixT, cMatrix},
aMatrix = ArrayReshape[Table[aElement[m, i, LatticeSize], {m, 0, LatticeSize - 1}, {i, 0, LatticeSize - 1}], {LatticeSize, LatticeSize}];
bMatrix = ArrayReshape[Table[bElement[i, j, k, p], {i, 0, LatticeSize - 1}, {j, 0, LatticeSize - 1}, {k, 0, LatticeSize - 1}, {p, 0, LatticeSize - 1}], {LatticeSize^3, LatticeSize}];

bMatrixT = Transpose[bMatrix];
cMatrix = ArrayReshape[Table[cElement[m, n, o, i, j, k, LatticeSize], {m, 0, LatticeSize - 1}, {n, 0, LatticeSize - 1}, {o, 0, LatticeSize - 1}, {i, 0, LatticeSize - 1}, {j, 0, LatticeSize - 1}, {k, 0, LatticeSize - 1}], {LatticeSize^3, LatticeSize^3}];

matrix = ArrayFlatten[{{aMatrix, bMatrixT}, {bMatrix, cMatrix}}];

matrix
];


where the methods defining the calculations are simple If-constructs of the form

aElement[m_, i_, LatticeSize_] := If[m == Mod[i - 1, LatticeSize], -t, 0.] + If[m == Mod[i + 1, LatticeSize], -t, 0.] + If[m == i, U/2, 0.]
bElement[i_, j_, k_, p_] := If[i == j == k == p, U/2, 0.];
cElement[m_, n_, o_, i_, j_, k_, LatticeSize_] := If[m == i && o == k && n == Mod[j + 1, LatticeSize] || n == j && o == k && m == Mod[i + 1, LatticeSize], -t, 0.] + If[m == i && o == k && n == Mod[j - 1, LatticeSize] || n == j && o == k && m == Mod[i - 1, LatticeSize], -t, 0.] + If[m == i && o == Mod[k + 1, LatticeSize] && n == j, t, 0.] + If[m == i && o == Mod[k - 1, LatticeSize] && n == j, t, 0.] + If[m == i && o == k && n == j, U/2, 0.];


where LatticeSize equals the formerly described input parameter $N$. Where is the bottleneck in my code? I am not able to see where I lose time in computation.

• This code does not run for me when I input BuildMatrix. Hard to give any ideas without working code. Very hard. – Daniel Lichtblau Dec 22 '16 at 20:10
• Are you sure you don't mean $N^6$? I get that the matrix has dimensions $(N+N^3)$ by $(N+N^3)$. – fred Dec 22 '16 at 20:29
• You have failed to define bElement and cElement. – Feyre Dec 22 '16 at 20:34
• @fred You are absolutely right. I noted $N^6$ in my handwritten sketches here according to exactly what you said but for some reason I didn't write it in my post. It's corrected now. Now the scaling is even worse as can be seen in the timing approximation for $N=20.$ – pbx Dec 22 '16 at 21:21
• @DanielLichtblau Depending on the simplicity of the …Element calculation methods I thought the problem has to be connected to my way of working with huge chunks of data (as can be seen in the given, not very extensive code). Usually, it's said here "focus on what is of interest". But I see it might help to also include the …Element calculations. See above now ;) – pbx Dec 22 '16 at 21:29

[Too long for a comment.]

It seems to scale more like O(n^2) and the timings I get reflect that.

--- edit #2 ---

That was not correct. Both speed and memory scale as O(n^6). I would surmise that eventually speed further degrades as size gets to the realm where RAM might be exceeded.

A related thing I should mention is that it is useful to make sure the result is a packed array. I changed the code slightly to this end.

t = 1.; U = 4.;
BuildMatrix[LatticeSize_] :=
Module[{m, n, o, i, j, k, matrix, aMatrix, bMatrix, bMatrixT,
cMatrix},
aMatrix =
DeveloperToPackedArray[
ArrayReshape[
Table[aElement[m, i, LatticeSize], {m, 0, LatticeSize - 1}, {i,
0, LatticeSize - 1}], {LatticeSize, LatticeSize}]];
bMatrix =
DeveloperToPackedArray[
ArrayReshape[
Table[bElement[i, j, k, p], {i, 0, LatticeSize - 1}, {j, 0,
LatticeSize - 1}, {k, 0, LatticeSize - 1}, {p, 0,
LatticeSize - 1}], {LatticeSize^3, LatticeSize}]];
bMatrixT = Transpose[bMatrix];
cMatrix =
ArrayReshape[
Table[cElement[m, n, o, i, j, k, LatticeSize], {m, 0,
LatticeSize - 1}, {n, 0, LatticeSize - 1}, {o, 0,
LatticeSize - 1}, {i, 0, LatticeSize - 1}, {j, 0,
LatticeSize - 1}, {k, 0, LatticeSize - 1}], {LatticeSize^3,
LatticeSize^3}];
matrix =
DeveloperToPackedArray[
ArrayFlatten[{{aMatrix, bMatrixT}, {bMatrix, cMatrix}}]];
matrix];


Still behaving at i=20 but I would not expect this to persist much beyond that level.

exp[baseTime_, i_] := baseTime*(i/6)^6;
Table[time = AbsoluteTiming[mat = BuildMatrix[i]][];
{{exp[baseTime, i], time}, {exp[baseSize, i], N@ByteCount[mat]}}, {i,
19, 20}]

(* Out= {{{317.340883389, 293.817071}, {3.97720005306*10^8,
3.78455224*10^8}}, {{431.702331962,
415.909048}, {5.41048010974*10^8, 5.14563352*10^8}}} *)


Again, creating this as a sparse array would almost certainly behave better in terms of speed and memory.

--- end edit #2 ---

Table[
time = AbsoluteTiming[mat = BuildMatrix[dim];];
dim = Dimensions[mat];
{time, dim, Times @@ dim, Tally[Flatten[mat]]}, {dim, 6, 12}]

(* Out= {{{0.314719, Null}, {222, 222},
49284, {{2., 234}, {-1., 876}, {0., 47742}, {1., 432}}}, {{0.754831,
Null}, {350, 350},
122500, {{2., 364}, {-1., 1386}, {0., 120064}, {1.,
686}}}, {{1.667838, Null}, {520, 520},
270400, {{2., 536}, {-1., 2064}, {0., 266776}, {1.,
1024}}}, {{3.393488, Null}, {738, 738},
544644, {{2., 756}, {-1., 2934}, {0., 539496}, {1.,
1458}}}, {{6.363055, Null}, {1010, 1010},
1020100, {{2., 1030}, {-1., 4020}, {0., 1013050}, {1.,
2000}}}, {{11.218442, Null}, {1342, 1342},
1800964, {{2., 1364}, {-1., 5346}, {0., 1791592}, {1.,
2662}}}, {{19.046299, Null}, {1740, 1740},
3027600, {{2., 1764}, {-1., 6936}, {0., 3015444}, {1., 3456}}}} *)


From the above one will notice that most entries are the same, 0. 0, so a possibility for improving speed and memory requirements would be to change the construction to produce a SparseArray.

--- edit ---

Here I repair the time estimator, which also works as a size estimator since the two in this case seem to scale linearly.

exp[baseTime_, i_] := baseTime*(i/6)^6;
Table[time = AbsoluteTiming[mat = BuildMatrix[i]][];
{{exp[baseTime, i, 6], time}, {exp[baseSize, i, 6],
N@ByteCount[mat]}}, {i, 5, 15}]
(* Out= {{{0.105396077139, 0.112266}, {132091.799554,
135352.}}, {{0.314711, 0.307111}, {394424.,
394424.}}, {{0.793583557077, 0.757408}, {994589.960048,
980152.}}, {{1.76825275171, 1.671213}, {2.21613265295*10^6,
2.163352*10^6}}, {{3.58475498438, 3.369494}, {4.492735875*10^6,
4.357304*10^6}}, {{6.7453489369, 6.317947}, {8.45387517147*10^6,
8.160952*10^6}}, {{11.949797108, 11.128567}, {1.49765555526*10^7,
1.4407864*10^7}}, {{20.141504, 18.785613}, {2.5243136*10^7,
2.4220952*10^7}}, {{32.5585109568, 30.261188}, {4.08052407625*10^7,
3.9072952*10^7}}, {{50.7893476529,
47.333028}, {6.36537574431*10^7, 6.0852664*10^7}}, {{76.8337402344,
71.695688}, {9.6294921875*10^7, 9.1936952*10^7}}} *)


--- end edit ---

• Maybe I got something completely wrong but please have a look at the code below which uses $N=6$ als calculation reference and on my machine produces: {{0.590009, 1.09306, 1.09899}, {0.770624, 2.43555, 2.45416}, {0.975321, 4.93756, 4.989}, {1.2041, 9.2909, 9.27697}, {1.45696, 16.4594, 16.3822}, {1.7339, 27.7425, 27.5057}} and thus a scaling that is much more like $N^6$ than $N^2$: exp[baseTime_, i_, n_] := baseTime*(i/6)^n; baseTime = AbsoluteTiming[BuildMatrix][]; Table[{exp[baseTime, i, 2], exp[baseTime, i, 6], AbsoluteTiming[BuildMatrix[i]][]}, {i, 7, 12}] – pbx Dec 23 '16 at 17:26
• It's a bit difficult to see all this code in a comment. Maybe better to edit the original question to show what it is you are measuring. – Daniel Lichtblau Dec 23 '16 at 17:36
• I included my thoughts into the initial question. Remark: I use the time for $N=6$ to make predictions according to assumed scaling behaviour. At least on my machine the assumption of a behaviour $\propto N^6$ fits much better. So I am not sure where your $\propto N^2$ scaling arises from. – pbx Dec 23 '16 at 17:47
• Umm, 6^n and n^6 are very different beasts. – Daniel Lichtblau Dec 23 '16 at 19:04
• I don't understand your comment to be honest. At first, your code definitely won't work this way (your exp awaits two arguments, you provide three), secondly, please have a more thorough look at what I am doing above, the n is simply a parameter for the estimated exponent in my exp function. And I compare scaling according to $N^2$ (n=2) and $N^6$ (n=6) with one another to show you your initial estimation with a quadratic scaling is wrong. The (i/6)^n` comes from the fact that I measure the baseTime in the $N=6$ case and has nothing to do with overall scaling. – pbx Dec 23 '16 at 22:55