Overlaying Two Matrices So That Sum Of Squared Differences Is Minimized
I try to summarize the results of the discussion in the comments:


Using a heap instead of a list would be better suited for the proposed idea since all operations to be performed on the list are sorting, returning the smallest element and changing an element followed by resorting the list.  These are almost exactly the operations supported by heaps, which will deliver better performance.
The sums $\sum_{k,l} (L_{i+k,j+l}-S_{k,l})^2$ can be expressed as $A_{i,j}+B_{i,j}-2C_{i,j}$ where $A_{i,j}:=\sum_{k,l} L^2_{i+k,j+l}$, $B_{i,j}:=\sum_{k,l} S^2_{k,l}$ and $C_{i,j}:=\sum_{k,l}L_{i+k,j+l}S_{k,l}$.  Here $B_{i,j}$ does not depend on $i,j$ and $A_{i,j}$ can obtained from its predecessor $A_{i',j'}$ by only updating two rows or columns, if the order in which the $i,j$ are perused is chosen wisely.  Most importantly, $C_{i,j}$ describes a convolution, which can be calculated very fast by FFTs, convolution in the frequency domain corresponds to (pointwise) multiplication in the time domain!  But even without implementation of FFT, this algorithm scores better than the naive approach, as it only performs one operation (multiplication) instead of two (subtraction and squaring) in the large (size $HW$) sums.