I think I understand what you are trying to do. I have no idea of the relative performance of our two machines so maybe you can benchmark it yourself.

from PIL import Image
import numpy as np

# Load images, convert to RGB, then to numpy arrays and ravel into long, flat things
a=np.array(Image.open('a.png').convert('RGB')).ravel()
b=np.array(Image.open('b.png').convert('RGB')).ravel()

# Calculate the sum of the absolute differences divided by number of elements
MAE = np.sum(np.abs(np.subtract(a,b,dtype=np.float))) / a.shape[0]

The only "tricky" thing in there is the forcing of the result type of np.subtract() to a float which ensures I can store negative numbers. It may be worth trying with dtype=np.int16 on your hardware to see if that is faster.

A fast way to benchmark it is as follows. Start ipython and then type in the following:

from PIL import Image
import numpy as np

a=np.array(Image.open('a.png').convert('RGB')).ravel()
b=np.array(Image.open('b.png').convert('RGB')).ravel()

Now you can time my code with:

%timeit np.sum(np.abs(np.subtract(a,b,dtype=np.float))) / a.shape[0]
6.72 µs ± 21.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Or, you can try an int16 version like this:

%timeit np.sum(np.abs(np.subtract(a,b,dtype=np.int16))) / a.shape[0]
6.43 µs ± 30.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

If you want to time your code, paste in your function then use:

%timeit compare_images_pil(img1, img2)

You can simply use NumPy functions to get the desired result:

np.sum(np.abs(np.diff(time)))

This works according to your desired formula even though np.diff calculates the difference time[i+1] - time[i] (instead of time[i] - time[i+1]) because you're using the absolute value.

Because this uses NumPy functions on a NumPy array it's probably much faster than any Python comprehension and/or functions.

I can offer an O(n log n) solution for a start: Let fi be the i-th number of the result. We have:

When walking through the array from left to right and maintain a binary search tree of the elements a0 to ai-1, we can solve all parts of the formula in O(log n):

• Keep subtree sizes to count the elements larger than/smaller than a given one
• Keep cumulative subtree sums to answer the sum queries for elements larger than/smaller than a given one

We can replace the augmented search tree with some simpler data structures if we want to avoid the implementation cost:

• Sort the array beforehand. Assign every number its rank in the sorted order
• Keep a binary indexed tree of 0/1 values to calculate the number of elements smaller than a given value
• Keep another binary indexed tree of the array values to calculate the sums of elements smaller than a given value

TBH I don't think this can be solved in O(n) in the general case. At the very least you would need to sort the numbers at some point. But maybe the numbers are bounded or you have some other restriction, so you might be able to implement the sum and count operations in O(1).

An implementation:

# binary-indexed tree, allows point updates and prefix sum queries
class Fenwick:
  def __init__(self, n):
    self.tree = [0]*(n+1)
    self.n = n
  def update_point(self, i, val):  # O(log n)
    i += 1
    while i <= self.n:
      self.tree[i] += val
      i += i & -i
  def read_prefix(self, i):        # O(log n)
    i += 1
    sum = 0
    while i > 0:
      sum += self.tree[i]
      i -= i & -i
    return sum

def solve(a):
  rank = { v : i for i, v in enumerate(sorted(a)) }
  res = []
  counts, sums = Fenwick(len(a)), Fenwick(len(a))
  total_sum = 0
  for i, x in enumerate(a):
    r = rank[x]
    num_smaller = counts.read_prefix(r)
    sum_smaller = sums.read_prefix(r)
    res.append(total_sum - 2*sum_smaller + x * (2*num_smaller - i))
    counts.update_point(r, 1)
    sums.update_point(r, x)
    total_sum += x
  return res

print(solve([3,5,6,7,1]))  # [0, 2, 4, 7, 17]
print(solve([2,0,1]))      # [0, 2, 2]

Is there a faster way to do the calculations? The lists are read and parsed from text files and an increased preprocessing time would be no problem (such as creating numpy arrays). The lists contain floating point numbers and are of equal size which is known beforehand.

PyPy is very good at optimizing list accesses, so you should probably stick to using lists.

One thing that will help PyPy optimize things is to make sure your lists always have only one type of objects. I.e. if you read strings from a file, don't put them in a list, then parse them into floats in-place. Rather, create the list with floats, for example by parsing each string as soon as it is read. Likewise, never try to preallocate a list, especially with [None,]*N, or PyPy will not be able to guess that all the elements have the same type.

Second, iterate the list as few times as possible. Your np_sum function walks both arrays three times (subtract, abs, sum) unless PyPy notices and can optimize it. Both 1. and 2. walk the list once, so they are faster.

First of all, I want to say this is really well presented. It just looks neat, with docstrings and comments and such. If you really want to perfect it, check the Python conventions in PEP 8 and linked documents. For example,

An inline comment is a comment on the same line as a statement. Inline comments should be separated by at least two spaces from the statement. They should start with a # and a single space.

and

The docstring is a phrase ending in a period. It prescribes the function or method's effect as a command ("Do this", "Return that"), not as a description; e.g. don't write "Returns the pathname ...".

I also very much appreciate the case analysis that has gone into the code to define those four cases, and so to reduce the naively quadratic function to a linear one. (Refering back to comments, that case analysis would not be at all amiss in a comment just after the docstring explaining the internal logic of the function)

Because, per the code and your case analysis, max1 and min1 always contain the same data, you could make the code shorter and reduce its space requirements by only having one of those. (And likewise for max2 and min2) However, see the next section.

In terms of the code itself, I suggest getting rid of those arrays. If they are only filled so that you can reduce them to a single maximum or minimum value, you might as well just track the rolling maximum or minimum as you go. This definitely saves space and probably saves time.

I would name the rolling reductions max_sum and max_difference rather than max1 and max2.

In terms of testing examples, do look for pathological cases. One case that always deserves to be checked, although it's not always obvious what the behaviour should be, is the case of an empty list.

Finally, it is generally worth using Python 3 in preference to Python 2 for new code, and especially practice code.

If you want the absolute element-wise difference between both matrices, you can easily subtract them with NumPy and use numpy.absolute on the resulting matrix.

import numpy as np

X = [[12,7,3],
[4 ,5,6],
[7 ,8,9]]

Y = [[5,8,1],
[6,7,3],
[4,5,9]]

result = np.absolute(np.array(X) - np.array(Y))

Outputs:

[[7 1 2]
 [2 2 3]
 [3 3 0]]

Alternatively (although unnecessary), if you were required to do so in native Python you could zip the dimensions together in a nested list comprehension.

result = [[abs(a-b) for a, b in zip(xrow, yrow)]
          for xrow, yrow in zip(X,Y)]

Outputs:

[[7, 1, 2], [2, 2, 3], [3, 3, 0]]

If m is a square matrix, this will do the job:

import numpy as np
np.sum((m-m.T)**2/(m+m.T))/2

Here is a function that covers the case in which there is 0 in the denominator:

def find_s(m):
    d=(m+m.T)
    off_diag_indices=np.triu_indices(len(d),1)
    if 0 in d[off_diag_indices]:
        return 'NA'
    else:
        numerator=(m-m.T)**2
        denominator=m+m.T
        return np.sum(numerator[off_diag_indices]/denominator[off_diag_indices])

The reason that I used off_diag_indices is because we actually do allow to have 0 on the diagonal of m+m.T, because we never sum the elements on the diagonal.