📖 Algorithms and Complexity

📖 Algorithms and Complexity#

⏱ | words

Task1.1 Task1.2 Task1.3 Task1.4 Task1.5 References

Writing programs that work fast#

Example: evaluation of polynomials

Task: evaluate the polynomial of the form at a given $x$

\[y = a_1 + a_2 x + a_3 x^2 + \dots + a_k x^k\]

Explain the operation of this program.
Can you make this algorithm more efficient?

Example: next

Here is a better approach!

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

Why is this algorithm faster? What is the difference?

Practical Task 1.1: Evaluations the run time of the two polynomial algorithms.

Complete the coding assignment in the exercises repo in the Jupyter notebook 1_algorithms/task1.1polynomials_pre.ipynb

An algorithm is a method of solving a class of problems on a computer.

sequence of steps/commands for the computer to run

Relevant questions:

How much time does it take to run?
How much memory does it need?
What other resources may be limiting? (storage, communication, etc)

Smart algorithm is a lot more important that 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, heterogenous agents, overlapping generations, etc.)
Computationally heavy econometrics (Bayesian, MCMC, multi-dimensional fixed effects, etc.)
Structural estimation with the need to resolve the model solution many thousand times
Counterfactual analysis, sensitivity analysis and uncertainty quantification

Estimation of static and dynamic games is one of the areas of structural econometrics requiring quick computation $\implies$ smart algorithms

Algorithms with different complexity#

Complexity of an algorithms in the cost, measured in running time or in storage requirement, of using an 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)

# 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	=	1	1	1
4	=	1	0	0
7 AND 4	=	1	0	0 = 4


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


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

Bit shifts in Python

_images/20d91cc74c67709f1876c7fb179324af8472caac01cef99d7e87e2fd00bf0415.png

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

_images/8718dde4543a326f6100b99b26b31e2e07e0105dc44054d31a3a7e751134206a.png

Binary search in finite set#

Finding a discrete element between given boundaries

Example

Think of number between 1 and 100
How many guesses are needed to locate it if the only answers are “below” and “above”?
What is the optimal sequence of questions?

Explain the operation of the code below

Inputs: sorted list of numbers, and a value to find
Algorithm:
1. Find middle point  
1. If the sought value is below, reduce the list to the lower half
1. 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[7]=[84] in [ 6 15 21 34 35 61 80 84 91 92] is 7

_images/d3c6a2f4572ec400f2ea20f4f1a6206c3f8cc41f8f4352e643277407db35087a.png

_images/9488fd548e5c3f1b444ba8c6ca7e4e73fad3b46b947c349dbdadf36bf8676e1a.png

Rate of growth and big-O notation#

Very useful way to talk about the rate of growth $\leftrightarrow$ complexity of algorithms

Definition

\[f(n)=O\big(g(n)\big) \text{ as } n \to \infty \Leftrightarrow\]

\[\exists M>0 \text{ and } N \text{ such than } |f(n)| < M g(n) \text{ for all } n>N\]

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 order or 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. polinomial 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

Example

Examples of DAC algorithms:

Binary search
Quicksort and merge sort
Fast Fourier transform (FTT) algorithm
Karatsuba fast multiplication algorithm

Curse of dimensionality#

Example of a bad algorithm?

Definition

The term curse of dimensionality relates to the above exponential complexity of an algorithm.

Example

Many board games (checkers, chess, shogi, go) in their $n$-by-$n$ generalizations
Traveling salesman problem (TSP)
Many problems in economics are subject to curse of dimensionality 😢

Allocation of discrete good#

Maximize welfare $W(x_1,x_2,\dots,x_n)$ subject to $\sum_{i=1}^{n}x_i = A$ where $A$ is 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 maximize $W(\Lambda p_1,\Lambda p_2,\dots,\Lambda p_n)$ subject to above.

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

# 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]

_images/39369b27473a2f4807cbd5b2e50767c58a9a6b5e04402d3144b772c6dbebc05a.png

_images/10a1fe6f49720c0f9bea872bbd9812d3c4df4879cb31286a7222c3b10a256053.png

Hint

What to do with heavy to compute models?

Design of better solution algorithms
Analyze special classes of problems + rely on problem structure
Speed up the code (low level language, compilation to machine code)
Parallelize the computations
Bound the problem to maximize model usefulness while keeping it tractable
Wait for innovations in computing technology (quantum computing, etc.)

Classes of computational complexity in theoretical computer science

Thinking of all problems there are:

P can be solved in polynomial time
NP solution can checked in polynomial time, even if requires 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 Clay Mathematics Institute)

Recursion#

Definition

Recursive algorithm is an algorithm that calls itself in order to solve a problem

Surprisingly powerful technique in scientific programming!

Example

Fibonacci sequence defined as

\[x_k = x_{k-1} + x_{k-2}, k>2, \quad x_0 = 1, \quad x_1 = 1\]

Imagine a program that computes Fibonacci numbers using this definition, and calls itself in the process

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

Is this an efficient algorithm? Why or why not?

Towers of Hanoi problem#

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 place 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)

Practical Task 1.2

Code up the recursive solution using the algorithm above in the exercises repo in the Jupyter notebook 1_algorithms/task1.2_hanoi_pre.ipynb

Solution for 4 disks requires 13 steps:

Bisection method#

The first of two very important classic algorithms for equation solving

Solve equations of the form (focus on scalar case today)

\[f(x) = 0, \quad x \in [a,b] \subset \mathbb{R}, \; f(a)f(b)<0\]

The latter condition requires that the function $f(x)$ takes different signs at the endpoints $a$ and $b$

Algorithm is similar to binary search, but in continuous space

Input: function f(x)
       brackets [a,b] such that f(a)f(b)<0
       convergence tolerance epsilon
       maximum number of iterations max_iter

Algorithm:
  step 0: ensure all conditions are satisfied
  step 1: compute the sign of the function at (a+b)/2
  step 2: replace a with (a+b)/2 if f(a)f((a+b)/2)>0, otherwise replace b with (a+b)/2
  step 3: repeat steps 1-2 until |a-b|< epsilon, or max_iter number of iterations is reached
  step 4: return (a+b)/2

Converged in 22 steps

_images/0328d2b92421ab6115c22b3b3ff1a9fccf63bc50e56355725a2548cbeb849fe1.png

Practical Task 1.3: Implementing bisections method

Complete the coding assignment in the exercises repo in the Jupyter notebook 1_algorithms/task1.3_bisections_pre.ipynb

Newton-Raphson method#

The second of the two classic methods for solving an equation $f(x)=0$, gradient based

General form

\[f(x)=0\]

Equation solving
Finding maximum/minimum based on FOC, then $ f(x)=Q'(x) $

Derivation for Newton method using Taylor series expansion#

\[ f(x) = \sum_{k=0}^{\infty} \frac{f^{(k)}(x_0)}{k!} (x-x_0)^k \]

Take first two terms, assume $ f(x) $ is solution, and let $ x_0=x_i $ and $ x=x_{i+1} $

\[ 0 = f(x) = f(x_i) + f'(x_i) (x_{i+1}-x_i) \quad \Rightarrow \quad x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)} \]

The main ides of Newton-Raphson method is to iterate on the equation starting from some $x_0$

\[ x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}, \; i=1,2,\ldots \]

Applicable to the system of equations, in which case $x\in\mathbb{R}^n$ and $f: \mathbb{R}^n \to \mathbb{R}^n$

Input: function f(x)
       gradient function f'(x)
Algorithm:
1. Start with some good initial value
2. Update x using Newton step above
3. Iterate until convergence

_images/b3f89fad6c64a35ea7d1d7d1174904c778465250f620815cc8eb754981cee994.png

_images/b0a17936fa013cb781af885ae69ffacb5e4b2c37a1aee5384fe21043d276aedc.png

_images/86056dc4f11db934d9fd32cd4aef75c00f7494bd2fe052d968d2182a4ec3001a.png

_images/9855d8e6b07f3026bffe95c5c68f6b0e6abc13773253a2564d9fb33e22e71775.png

_images/c7573805d21852a492b4bc887f861e00a62c7be50fa0fc9832f83b5bd48c1176.png

_images/72063d0192a1944d7eee5acc82fa89b563a8ceab66994110892f1718a918b282.png

_images/614281950d9944f572c4d930dc0d308adf494e023d3ffbfc8f2db9f24d83d428.png

Converged in 7 steps

Practical Task 1.4: Implementing Newton-Raphson method

Complete the coding assignment in the exercises repo in the Jupyter notebook 1_algorithms/task1.4_newton_pre.ipynb

Practical Task 1.5: Multivariate Newton method [optional]

Complete the coding assignment in the exercises repo in the Jupyter notebook 1_algorithms/task1.5_multivariate_pre.ipynb

Measuring complexity of Newton and bisection methods#

What is the size of input $ n $?
Desired precision of the solution!
Thus, attention to the errors in the solution as algorithm proceeds
Rate of convergence is part of the computational complexity of the algorithms

Computational complexity

Calculating a root of a function f(x) with n-digit precision
Provided that a good initial approximation is known
Is $ O((logn)F(n)) $, where $ F(n) $ is the cost of
calculating $ f(x)/f'(x) $ with $ n $-digit precision

References and Additional Resources

📖 Wilf [2002] “Algorithms and Complexity” {download}`pdf of the book https://www2.math.upenn.edu/~wilf/AlgoComp.pdf ’
Complexity classes and P vs. NP
📺 Lecture on algorithm complexity by Erik Demaine, MIT Lecture recording, 50 min
Big-O cheat sheet https://www.bigocheatsheet.com
Bitwise operations post on Geeksforgeeks link
On computational complexity of Newton method link
“Improved convergence and complexity analysis of Newton’s method for solving equations” link
📺 Oscar Veliz videos on Newton method and its domains of attraction

📖 Algorithms and Complexity

Contents

📖 Algorithms and Complexity#

Writing programs that work fast#

Algorithms with different complexity#

Parity of a number#

Finding max/min of a list#

Binary search in finite set#

Rate of growth and big-O notation#

Classes of algorithm complexity#

How many operations as function of input size?#

Divide-and-conquer algorithms#

Curse of dimensionality#

Allocation of discrete good#

Recursion#

Towers of Hanoi problem#

Bisection method#

Newton-Raphson method#

Derivation for Newton method using Taylor series expansion#

Measuring complexity of Newton and bisection methods#