# Optimizing a series of basic matrix operations—scaling doesn't seem right with matrix size

I have a simple operation that I'd like to apply to an n×n matrix. My problem (well, one of them...) is that it runs in, say, 4 seconds for a 100×100 matrix, but it took five and a half minutes to run on a 300×300 matrix. It's a simple function I'm applying to each entry so I have no idea why this is happening, but I'm also pretty new to having to do anything efficiently.

The setup is this: I start with a random array of an equal number of 1s and -1s. I apply something like a Hamiltonian operator that looks at each element, does a calculation, and returns a number—namely, H[array] gives a matrix of the same dimensions of array.

All H does is look at the [[i,j]] entry of our input array ar, then multiplies that entry by all of its neighbors and takes the sum of those products. That sum is the [[i,j]] the entry in H[array].

H[ar_] := -ParallelTable[
Sum[ArrayPad[ar, 1][[i + m, j + n]],
{m, -1, 1}, {n, -1, 1}] - 1,
{i, 2, Length[ar] + 1}, {j, 2, Length[ar] + 1},
];


I use ArrayPad on the input to pad the matrix with 0s on all sides—this way, the entries on the edge (which have fewer than 8 neighbors) don't throw errors, and the "non-existent" entries contribute 0s to the sum.

Because I calculate the product of one value with 9 other values, I sped up the algorithm by just factoring out the first value and multiplying that by the sum instead of calculating the sum of 9 products. I subtract 1 from the final sum because this method includes a term comprising the product of the [[i,j]] entry of ar with itself, which is either 1×1=1 or -1×-1=1.

So, it's not particularly complicated—but I can't, for the life of me, figure out how this scales with size of this matrix. I would expect a ~9× increase in time when going from 100×100 to 300×300, but instead got an ~80× increase. Could anybody help me to figure out why this is and how to help optimize it?

Thank you!

• There's a superfluous comma at the end of your code. Copy 'n paste error? Dec 15, 2015 at 20:50
• Hello, Ben, and welcome to the Mathematica StackExchange site. Thanks for the good, properly formatted first question! Don't forget to upvote good answers (and other people's questions) using the triangle above the number next to the post, and use the checkmark to "accept" the answer to your question that you think best answers it. Dec 16, 2015 at 13:23

I believe that this is the best that can be done. Here we completely throw away your code and use only the structure of the problem, posing it as a convolution that can be attacked using the FFT.

HFFT[ar_] := ar*ListConvolve[-{{1, 1, 1}, {1, 0, 1}, {1, 1, 1}}, ar, {2, -2}, 0]


The FFT is (efficiently) parallelized by default and has mild scaling with the stencil size. A 9000×9000 matrix is no problem, taking under 5 seconds on my computer.

• +1 Very fast indeed: 40-200 times faster than my solution. Dec 15, 2015 at 22:58

Sjoerd's approach is the right way to go. It can be improved further, to secure an approximately 25-fold speedup (on my 4-core computer), by substituting Compile for ParallelTable:

Hc[ar_] := Compile[{{ap, _Integer, 2}, {len, _Integer, 0}},
-Table[
ap[[i, j]]*Total@Flatten[ap[[i - 1 ;; i + 1, j - 1 ;; j + 1]], 2] - 1,
{i, 2, len + 1}, {j, 2, len + 1}
]


For a 300×300 matrix, it runs in less than 0.1s on my computer. (Faster on yours: my computer is rather old.)

Note that this is not parallelized. Parallelizing it (not using ParallelTable, but rather the ability of compiled functions to evaluate over multiple elements in a list in parallel) might result in even greater improvements.

I don't know why your function scales so badly, but I know that it can be improved considerable. For instance, the padding needs to be done once only, whereas you do it every step. Additionally, a combination of Total and Span(;;) is much quicker than your Sum:

H[ar_] :=
ParallelTable[
ap[[i, j]]*Total[ap[[i - 1 ;; i + 1, j - 1 ;; j + 1]], 2] - 1,
{i, 2, Length[ar] + 1}, {j, 2, Length[ar] + 1}
]
]


This is about 100 times faster.

Scaling of timing with array size (size goes up as the square of the array width):

timingData =
Table[H[RandomChoice[{-1, 1}, i {100, 100}]];//AbsoluteTiming//First, {i, 10}]

(* {0.135327, 0.452207, 1.05277, 1.85567, 2.85895, 4.86328, 6.66683,
8.67205, 12.5623, 15.6128} *)

fit[x_] = Fit[timingData, {1, x, x^2}, x];

Show[
ListPlot[timingData],
Plot[fit[x], {x, 1, 10}]
] • The timings can be greatly improved by compilation in place of using ParallelTable. Dec 15, 2015 at 21:51
• I'm not surprised. Compile usually works wonders if you do simple arithmetical stuff. Dec 15, 2015 at 22:00