most optimal way to perform double summation

Question

For example below, how could I optimally calculate Gxx assuming that there are easily hundred thousand operations performed? If I just do nested for loops to perform summation, it takes a very long time. Is there a way to use the sum function or something else to make this process much more efficient?:

For reference, my current psuedo code does this:

rrange=range(r-(rows-1)/2,r+(rows-1)/2)
crange=range(c-(cols-1)/2,c+(cols-1)/2)
Gxx=0
for rval,cval in product(rrange,crange):
        #sum( for x in range())
        Gxx+=(someval(rval,cval)-someval2)^2

Answer 1

For calculating a single element of G_xx(r,c) element, you can't optimize anything. That's logical: after all you don't know anything about the structure of x so you will have to read all x_j,i elements in the range.

However, things change if you need to calculate the entire matrix. In that case, you can reuse work from calculating the previous element.

First some basic mathematics. Since x-bar(r,c) doesn't depend on j or i, it's a constant in the process. Now we know that:

(a-b)^2 = a^2+b^2-2*a*b

With b if you substitute it back to x-bar a constant. So if you apply this to a sommation, you can state that:

---                ---
\                  \
/     (x_ji-b)^2 = /   x_ji^2-2*b*x_ji+b^2   = S-2*s*b-D*b^2
---                ---
i,j                i,j

With S the sum of squares of the x-elements and s the sum of the x elements.

Now if you look to a matrix, the first iteration, you use a certain domain of the matrix:

x  x  x  x  x  x  x  x  x
  /-------\
x |x  x  x| x  x  x  x  x
  |       |
x |x  x  x| x  x  x  x  x
  |       |
x |x  x  x| x  x  x  x  x
  \-------/
x  x  x  x  x  x  x  x  x

The next iteration, you only move the scope only one element further so:

x  x  x  x  x  x  x  x  x
     /-------\
x  x |x^ x^ x| x  x  x  x
     |       |
x  x |x^ x^ x| x  x  x  x
     |       |
x  x |x^ x^ x| x  x  x  x
     \-------/
x  x  x  x  x  x  x  x  x

The xs denoted with ^ are in other words reused. So what one can do is use some kind of sliding window.

The first time you thus calculate first the sum and the sum of squares of the elements (you store them). Next each time you move the "cursor" you subtract again the line column that is no longer in the scope, and add the column that appears in scope. A basic algorithm is thus:

for r in range (rmin,rmax) :
    sum = 0
    sumsq = 0
    jmin = r-(rows-1)/2
    jmax = r+(rows-1)/2

    #calculate sum and sum of square of the first
    imin = cmin-(cols-1)/2
    imax = cmin+(cols-1)/2
    for j,i in product(range(jmin,jmax),range(imin,imax)) :
        xji = x(j,i) #cache xji
        sum += xji
        sumsq += xji * xji
    d = (jmax-jmin)*(imax-imin)
    #now we can calculate the first element of the row
    xb = xbar(r,cmin)
    Gxx(r,cmin) = sumsq-2*sum*xb+d*xb*xb

    #now iterate over all elements (except the first)
    for c in range(cmin+1,cmax) :
        isub = c-1-(cols-1)/2 #column to remove, (previous column = -1)
        iadd = c+(cols-1)/2 #column to add
        for j in range(jmin,jmax) :
            xji = x(j,isub)
            sum -= xji
            sumsq -= xji*xji
            xji = x(j,iadd)
            sum += xji
            sumsq += xji*xji
        #Now the sums and the sum of squares are updated
        xb = xbar(r,c)
        Gxx(r,c) = sumsq-2*sum*xb+d*xb*xb

I think there will be some work to adapt the algorithm, but it should be doable. Furthermore please check first on a small instance whether it works correct. Small rounding errors are possible.

This will not give much difference if the cols and rows are small, but it case these are large, it can result in a huge boost.