📖 Algorithms and complexity

Everything later in this course — solving Bellman equations, nesting a solver inside a likelihood, finding all equilibria of a game — is limited by how fast the inner loop runs. This class is about what makes an algorithm fast, and how to tell before you write it.

Writing programs that work fast¶

def calc_polynomial_faster(qs=[0,], x=0.0):
  '''Evaluates the polynomial given by coefficients qs at given x.
  First coefficient qs[0] is a constant, last coefficient is for highest power.
  Faster than before!
  '''
  res, xpw = qs[0], x  # init result and power of x
  for i in range(1,len(qs)):  # start with second coefficient
    res += xpw * qs[i]
    xpw *= x
  return res

An algorithm is a method of solving a class of problems on a computer — a sequence of steps/commands for the computer to run.

Relevant questions:

1. How much time does it take to run?

2. How much memory does it need?

3. What other resources may be limiting? (storage, communication, etc.)

A smart algorithm is a lot more important than a fast computer

Professor Martin Grötschel, Konrad-Zuse-Zentrum für Informationstechnik Berlin, expert in optimization

“a benchmark production planning model solved using linear programming would have taken 82 years to solve in 1988, using the computers and the linear programming algorithms of the day. Fifteen years later — in 2003 — this same model could be solved in roughly 1 minute, an improvement by a factor of roughly 43 million. Of this, a factor of roughly 1,000 was due to increased processor speed, whereas a factor of roughly 43,000 was due to improvements in algorithms!”

Algorithms are behind any computation done in economics:

• Macro simulation models (growth, heterogeneous agents, overlapping generations, etc.)

• Computationally heavy econometrics (Bayesian, MCMC, multi-dimensional fixed effects, etc.)

• Structural estimation with the need to re-solve the model many thousands of times

• Counterfactual analysis, sensitivity analysis and uncertainty quantification

Structural estimation of dynamic models is one of the areas of econometrics requiring quick computation $\implies$ smart algorithms.

Algorithms with different complexity¶

Complexity of an algorithm is the cost, measured in running time or in storage requirement, of using the algorithm to solve one of the problems in the relevant class.

Let’s look at some particular algorithms.

Parity of a number¶

Check whether an integer is odd or even.

Algorithm:
Convert the number to binary
Check whether the last digit is 0 (number is even) or 1 (number is odd)

def parity(n, verbose=False):
  '''Returns 1 if passed integer number is odd
  '''
  if not isinstance(n, int): raise TypeError('Only integers in parity()')
  if verbose: print('n = ', format(n, "b"))  # print binary form of the number
  return n & 1  # bitwise and operation returns the value of last bit

# check parity of various numbers
for n in [2,4,7,32,543,671,780]:
  print('n = {0:5d} ({0:08b}), parity={1:d}'.format(n,parity(n)))

n =     2 (00000010), parity=0
n =     4 (00000100), parity=0
n =     7 (00000111), parity=1
n =    32 (00100000), parity=0
n =   543 (1000011111), parity=1
n =   671 (1010011111), parity=1
n =   780 (1100001100), parity=0

Some details on bitwise operations

Bitwise operations in Python

• bitwise AND &

• bitwise OR |

• bitwise XOR ^

• bitwise NOT ~ (including sign bit!)

• right shift >>

• left shift << (without overflow!)

Bitwise AND, OR and XOR

7	=	0	1	1	1
4	=	0	1	0	0
7 AND 4	=	0	1	0	0 = 4

7	=	0	1	1	1
4	=	0	1	0	0
7 OR 4	=	0	1	1	1 = 7

7	=	0	1	1	1
4	=	0	1	0	0
7 XOR 4	=	0	0	1	1 = 3

Bit shifts in Python

Source

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [9, 6]

N = 50
kk = lambda i: 10**(i+1)+i  # step formula
n,x,std = [0]*N,[0]*N,[0]*N # initialize data lists
for i in range(N):
  k = kk(i)  # input value for testing
  n[i] = k.bit_length() # size of problem = bits in number
  t = %timeit -n5000 -r100 -o -q parity(k)
  x[i] = t.average
  std[i] = t.stdev

plt.errorbar(n,x,std)
plt.xlabel('number of bits in the input argument', fontsize=14)
plt.ylabel('run time, sec', fontsize=14)
plt.title("Run times for parity check as function of number length in bits",fontsize=14)
plt.show()

Finding max/min of a list¶

Find max or min in an unsorted list of values.

Algorithm:
cycle through the list once saving the current extremum value

def maximum_from_list(vars):
  '''Returns the maximum from a list of values
  '''
  m=float('-inf')  # init with the worst value
  for v in vars:
    if v > m: m = v
  return m

Source

import numpy as np
N = 50
kk = lambda i: 2*i  # step formula
n,x,std = [0]*N,[0]*N,[0]*N # initialize data lists
for i in range(N):
  n[i] = kk(i) # size of the array
  vv = np.random.uniform(low=0.0, high=100.0, size=n[i])
  t = %timeit -n1000 -r100 -o -q maximum_from_list(vv)
  x[i] = t.average
  std[i] = t.stdev

plt.errorbar(n,x,std)
plt.xlabel('number of elements in the list', fontsize=14)
plt.ylabel('run time, sec', fontsize=14)
plt.title("Run times for max finder as function of the array length",fontsize=14)
plt.show()

Binary search in a finite set¶

Finding a discrete element between given boundaries.

def binary_search(grid=[0,1],val=0):
  '''Returns the index of val on the sorted grid
  '''
  i1,i2 = 0,len(grid)-1
  if val==grid[i1]: return i1
  if val==grid[i2]: return i2
  j=(i1+i2)//2
  while grid[j]!=val:
    if val>grid[j]:
      i1=j
    else:
      i2=j
    j=(i1+i2)//2  # divide in half
  return j

Inputs: sorted list of numbers, and a value to find
Algorithm:
1. Find middle point
2. If the sought value is below, reduce the list to the lower half
3. If the sought value is above, reduce the list to the upper half

import numpy as np
N = 10
# random sorted sequence of integers up to 100
x = np.random.choice(100,size=N,replace=False)
x = np.sort(x)
# random choice of one number/index
k0 = np.random.choice(N,size=1)
k1 = binary_search(grid=x,val=x[k0])
print(f'Index of x{k0}={x[k0]} in {x} is {k1}')

Index of x[9]=[99] in [ 5  6 10 26 32 38 58 67 98 99] is 9

Source

N = 50  # number of points
kk = lambda i: 100+(i+1)*500  # step formula
# precompute the sorted sequence of integers of max length
vv = np.random.choice(10*kk(N),size=kk(N),replace=False)
vv = np.sort(vv)

n,x,std = [0]*N,[0]*N,[0]*N   # initialize lists
for i in range(N):
  n[i] = kk(i)  # number of list elements
  # randomize the choice in each run to smooth out simulation error
  t = %timeit -n10 -r100 -o -q binary_search(grid=vv[:n[i]],val=vv[np.random.choice(n[i],size=1)])
  x[i] = t.average
  std[i] = t.stdev

plt.errorbar(n,x,std)
plt.xlabel('number of elements in the list', fontsize=14)
plt.ylabel('run time, sec', fontsize=14)
plt.title("Run times for binary search as function of the array length",fontsize=14)
plt.show()

plt.errorbar(n,x,std)
plt.xscale('log')
plt.xlabel('log(number of elements in the list)', fontsize=14)
plt.ylabel('run time, sec', fontsize=14)
plt.title("Run times for binary search as function of the LOG of array length",fontsize=14)
plt.show()

Rate of growth and big-O notation¶

A very useful way to talk about the rate of growth $\leftrightarrow$ complexity of algorithms.

In words, $f(x) = O\big(g(x)\big)$ simply means that as $x$ increases, $f(x)$ certainly does not grow at a faster rate than $g(x)$.

In measuring solution time we may distinguish performance in

• best (easiest to solve) case

• average case

• worst case ($\leftarrow$ the focus of the theory!)

Constants and lower terms are ignored because we are only interested in the order of growth.

Classes of algorithm complexity¶

• $O(1)$ constant time

• $O(\log_{2}(n))$ logarithmic time

• $O(n)$ linear time

• $O(n \log_{2}(n))$ quasi-linear time

• $O(n^{k}), k>1$ quadratic, cubic, etc. polynomial time $\uparrow$ tractable

• $O(2^{n})$ exponential time $\downarrow$ curse of dimensionality

• $O(n!)$ factorial time

How many operations as function of input size?¶

• Parity: just need to check the lowest bit, does not depend on input size $\Rightarrow O(1)$

• Maximum element: need to loop through elements once $\Rightarrow O(n)$

• Binary search: divide the problem in 2 each step $\Rightarrow O(\log(n))$

Divide-and-conquer algorithms¶

Divide-and-conquer structure is what typically marks an excellent algorithm.

Curse of dimensionality¶

An example of a bad algorithm?

Allocation of a discrete good¶

Maximize welfare $W(x_1,x_2,\dots,x_n)$ subject to $\sum_{i=1}^{n}x_i = A$ where $A$ is a discrete good that is only divisible in steps of $\Lambda$.

Let $M=A/\Lambda \in \mathbb{N}$. Let $p_i \in \{0,1,\dots,M\}$ such that $\sum_{i=1}^{n}p_i = M$.

Then the problem is equivalent to maximizing $W(\Lambda p_1,\Lambda p_2,\dots,\Lambda p_n)$ subject to the above.

$(p_1,p_2,\dots,p_n)$ is a composition of the number $M$ into $n$ parts.

def compositions(N,m):
    '''Iterable on compositions of N with m parts
    Returns the generator (to be used in for loops)
    '''
    cmp=[0,]*m
    cmp[m-1]=N  # initial composition is all to the last
    yield cmp
    while cmp[0]!=N:
        i=m-1
        while cmp[i]==0: i-=1  # find lowest non-zero digit
        cmp[i-1] = cmp[i-1]+1  # increment next digit
        cmp[m-1] = cmp[i]-1    # the rest to the lowest
        if i!=m-1: cmp[i] = 0  # maintain cost sum
        yield cmp

# example of compositions generation
for c in compositions(5,3) : print(c)

[0, 0, 5]
[0, 1, 4]
[0, 2, 3]
[0, 3, 2]
[0, 4, 1]
[0, 5, 0]
[1, 0, 4]
[1, 1, 3]
[1, 2, 2]
[1, 3, 1]
[1, 4, 0]
[2, 0, 3]
[2, 1, 2]
[2, 2, 1]
[2, 3, 0]
[3, 0, 2]
[3, 1, 1]
[3, 2, 0]
[4, 0, 1]
[4, 1, 0]
[5, 0, 0]

Source

N = 10  # number of points
kk = lambda i: 2+i  # step formula
M = 20  # quantity of indivisible good in units of lambda

n,x,std = [0]*N,[0]*N,[0]*N   # initialize lists
for i in range(N):
    n[i] = kk(i)  # number of list elements
    t = %timeit -n2 -r10 -o -q for c in compositions(M,n[i]) : pass
    x[i] = t.average
    std[i] = t.stdev

plt.errorbar(n,x,std)
plt.xlabel('Number of elements in compositions',fontsize=14)
plt.ylabel('run time, sec',fontsize=14)
plt.title('Run time as function of number of compositions',fontsize=14)
plt.show()

plt.errorbar(n,x,std)
plt.yscale('log')
plt.xlabel('Number of elements in compositions',fontsize=14)
plt.ylabel('log(run time)',fontsize=14)
plt.title('Curse of dimensionality in composition generation',fontsize=14)
plt.show()

Classes of computational complexity in theoretical computer science

Thinking of all problems there are:

• P can be solved in polynomial time

• NP solution can be checked in polynomial time, even if it requires an exponential solution algorithm

• NP-hard as complex as any NP problem (including all exponential and combinatorial problems)

• NP-complete both NP and NP-hard (tied via reductions)

NP stands for non-deterministic polynomial time $\leftrightarrow$ ‘magic’ guess algorithm.

P vs. NP

Unresolved question of whether P = NP or P $\ne$ NP ($1 mln. prize by the Clay Mathematics Institute)

Recursion¶

A surprisingly powerful technique in scientific programming — and the structure of every dynamic programming solver we write later in the course.

def fibonacci(n):
    if n == 0:
        return 1
    elif n == 1:
        return 1
    else:
        return fibonacci(n - 1) + fibonacci(n - 2)

for i in range(10):
    print(fibonacci(i),end=' ')

1 1 2 3 5 8 13 21 34 55 

Towers of Hanoi problem¶

A classic puzzle: given a board with three pegs, move a stack of disks of different size from the left-most peg to the right-most peg, moving one disk at a time and following the rule that no larger disk can be placed on top of a smaller one.

The problem can be solved nicely by breaking it into small parts using the following algorithm:

def move(from,to):
  move one disk from --> to

def move_via(from,via,to):
  move(from,via)
  move(via,to)

def main_algorithm(n,source,aux,target):
  '''
  Inputs: number of disks n
        source peg
        auxiliary peg
        target peg
  '''
  if n==0:
    do nothing, return
  if n==1:
    move(source,target)
  if n>0:
    main_algorithm(n-1,source,target,aux)
    move(source,target)
    main_algorithm(n-1,aux,source,target)

The solution for $n$ disks requires $2^n-1$ moves, so 15 for the four disks below — the illustration stops at configuration 13, two moves short of the goal.