Question

The first major optimization is to see if we can reduce the algorithmic complexity of our solution, say from \(\mathcal{O}(n^2)\) to \(\mathcal{O}(n \log(n))\). Unless you are going to invent an entirely new algorithm (possible but uncommon), this involves research into whether the data structures used are optimal, or whetther there is a way to reformulate the problem to take advantage of better algorithms. If your inital solution is by “brute force”, there is sometimes room for huge performacne gains here.

Taking a course in data structures and algorithms is a very worthwhile investment of your time if you are developing novel statsitical algorithms - perhaps Bloom filters, locality sensitive hashing, priority queues, Barnes-Hut partitionaing, dynamic programming or minimal spanning trees can be used to solve your problem - in which case you can expect to see dramatic improvements over naive brute-force implementations.

Suppose you were given two lists xs and ys and asked to find the unique elements in common between them.

xs = np.random.randint(0, 1000, 10000)
ys = np.random.randint(0, 1000, 10000)

# This is easy to solve using a nested loop

def common1(xs, ys):
    """Using lists."""
    zs = []
    for x in xs:
        for y in ys:
            if x==y and x not in zs:
                zs.append(x)
    return zs

%timeit -n1 -r1 common1(xs, ys)

1 loops, best of 1: 14.7 s per loop

# However, it is much more efficient to use the set data structure

def common2(xs, ys):
    return list(set(xs) & set(ys))

%timeit -n1 -r1 common2(xs, ys)

1 loops, best of 1: 2.82 ms per loop

assert(sorted(common1(xs, ys)) == sorted(common2(xs, ys)))

We have seen many such examples in the course - for example, numerical quadrature versus Monte Carlo integration, gradient desceent versus conjugate gradient descent, random walk Metropolis versus Hamiltonian Monte Carlo.