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

Question

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[
ArrayPad[ar, 1][[i, j]]*
  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!

Answer 1

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.

Answer 2

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}
   ]
  ][ArrayPad[ar, 1], Length[ar]]

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.

Answer 3

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_] := 
  -With[{ap = ArrayPad[ar, 1]}, 
     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}]
 ]

Nicely quadratic, as it should.