Application: Optimization and Approximation

This chapter includes some examples of using our Python knowledge to solve some challenging problems.

Assume that we have a "nice" function f(x) that we want to find the output of for a specific x, but for some reason we can't actually compute f(x). If we know the output of f for some input a, and if we can continously take the derivative of f with respect to x, then we can use the Taylor series expansion of f(x) to approximate that value.

Let's approximate a nice function f at point x using the Taylor series:

$$ f(x) = f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \frac{f'''(a)}{3!}(x-a)^3 + ... $$

where f is a simple, and continuous function that we can continuously differentiate to get the next terms. Normally, the Taylor series goes on forever, but generally we only compute a fixed number of terms because that is a good enough approximation for us.

To simplify things, let us assume that f(x) is a polynomial, so that the differentiation is easy. Polynomials are made of terms that are powers of $x$.

A term is just something like $a x^{n}$, where $a$ is a coefficient and $n$ represents the power of $x$. Its main use is twofold: A term is easily derivable, and a term needs to be able to return a value for a substition for $x$. .derive_many(num) is just a shorthand for applying .derive() in a loop.

class Term:
    '''Represents something like coef * x^{power}'''
    def __init__(self, coef, power):
        self.coef = coef
        self.power = power
    def derive(self):
        if self.power == 0:  # We just have a constant
            return Term(0, 0) # represents zero
        else:
            return Term(self.coef * self.power, self.power - 1)
    def derive_many(self, i):
        if i == 0:
            return self
        next_term = self.derive()
        return next_term.derive_many(i-1)
    def __repr__(self):
        return f'{self.coef}x^({self.power})'
    def subst(self, x):
        return self.coef * (x ** self.power)

Here are some examples of using this term class:

term1 = Term(3, 5) 

term1

3x^(5)

term1.derive()

15x^(4)

term1

3x^(5)

term1.subst(1)

3

term1.subst(2)

96

term1 = Term(3, 5)
print(term1)
print(term1.derive())
print(term1.derive().derive())
print(term1.derive().derive().derive())
print(term1.derive().derive().derive().derive())
print(term1.derive().derive().derive().derive().derive())
print(term1.derive().derive().derive().derive().derive().derive())
print(term1.derive().derive().derive().derive().derive().derive().derive())

3x^(5)
15x^(4)
60x^(3)
180x^(2)
360x^(1)
360x^(0)
0x^(0)
0x^(0)

As we see, the derivative of a constant polynomial is zero.

term1.derive_many(8)

0x^(0)

Given Term, a Polynomial is just a sum of different terms. We technically need to have a unique Term for each power, but this definition is enough for now:

class Polynomial:
    '''Represents a*x^{n} + b*x^{n-1} + ... + '''
    def __init__(self, term_list):
        self.term_list = term_list
    def __repr__(self):
        if not self.term_list:
            return '0'
        return ' + '.join(f'[{term.__repr__()}]' for term in self.term_list)
    def derive(self):
        new_terms = [term.derive() for term in self.term_list]
        return Polynomial([term for term in new_terms 
                           if not (term.coef == 0 and term.power == 0)])
    def derive_many(self, i):
        if i == 0:
            return self
        next_poly = self.derive()
        return next_poly.derive_many(i-1)
    def subst(self, x):
        return sum([term.subst(x) for term in self.term_list])

Writing .derive() and .subst() are easy, because the main task is just calling .derive() and .subst() of the Terms inside.

Let's try it out in action:

poly1 = Polynomial([Term(1, 5), Term(2, 4), Term(3, 3)])  # x^5 + 2x^4 + 3x^3
print(poly1)
print(poly1.derive())

print(poly1.derive_many(3))
print(poly1.subst(2))

[1x^(5)] + [2x^(4)] + [3x^(3)]
[5x^(4)] + [8x^(3)] + [9x^(2)]
[60x^(2)] + [48x^(1)] + [18x^(0)]
88

Now that our polynomial is ready, we can code the rest of the Taylor approximation code as follows:

def factorial(x):
    return 1 if x == 0 else x * factorial(x-1)

def approx(poly, x, a, m):
    """Approximates poly(x) by using fixed point a, with up to m terms."""
    total = poly.subst(a)  # f(a)
    for i in range(1, m+1):
        # Compute the next term
        f_derived = poly.derive_many(i)
        value = (f_derived.subst(a) * (x - a) ** i) / factorial(i)
        total += value
    return total
        
        

terms = [
    Term(1, 3),  # x^3
    Term(2, 5),  # 2x^5
    Term(4, 4)   # 4x^4
]
poly = Polynomial(terms)
a = 1
m = 5

x = 5

print('Actual:', poly.subst(x))

for i in range(6):
    print(f'Approximate ({i+1} terms):', approx(poly, x, a, i))

Actual: 8875
Approximate (1 terms): 7
Approximate (2 terms): 123.0
Approximate (3 terms): 875.0
Approximate (4 terms): 3243.0
Approximate (5 terms): 6827.0
Approximate (6 terms): 8875.0

Looks like it works well. Let's also draw graphs showing the approximation performance of the Taylor series with increasing number of used terms:

import numpy as np
xs = np.linspace(-4, 4, 50)

yhats_m1 = [approx(poly, x, 1, 1) for x in xs] # 1 term
yhats_m2 = [approx(poly, x, 1, 2) for x in xs] # 2 terms
yhats_m3 = [approx(poly, x, 1, 3) for x in xs] # 3 terms
yhats_m5 = [approx(poly, x, 1, 5) for x in xs] # 5 terms

import matplotlib.pyplot as plt

poly

[1x^(3)] + [2x^(5)] + [4x^(4)]

plt.plot(xs, [poly.subst(x) for x in xs], label='actual')
plt.plot(xs, yhats_m1, label='approximate (m=1)')
plt.plot(xs, yhats_m2, label='approximate (m=2)')
plt.plot(xs, yhats_m3, label='approximate (m=3)')
plt.plot(xs, yhats_m5, label='approximate (m=5)')

plt.legend()

<matplotlib.legend.Legend at 0x7fbb002c92b0>

It's a bit hard to see, but the blue and purple lines are on top of each other, meaning that the 5-term approximator fits this polynomial perfectly.

Finding a root of a polynomial

A similar application is using iterative refinement to find a root to a polynomial expression that is set to zero. This time, following the explanation in the textbook, the core idea that we need to code is as follows: We start off with a random initial guess for a root ($x_{0}$), and then in a loop we continuously obtain better estimates of a root through this equation:

$$ x_{i+1} = x_{i} - \frac{f(x_{i})}{f'(x_{i})} $$

Let's show an example polynomial that has some "nice" roots:

poly = Polynomial([
    Term(1, 2),
    Term(4, 1),
    Term(-5, 0)
])

xs = np.linspace(-10, 10, 100)
ys = [poly.subst(x) for x in xs]
plt.plot(xs, ys)
plt.hlines(y=0, xmin=min(xs), xmax=max(xs), color='black')
plt.vlines(x=0, ymin=min(ys), ymax=max(ys), color='black')

<matplotlib.collections.LineCollection at 0x7fbb20b5e820>

The parabola above has roots at -5 and 1. The function below is an implementation of the "find-a-root" idea that we talked about.

One detail is that since we're updating our root guess every loop, we would like to stop if we're "stuck" in some value. However, since we're working with floating point numbers, checking for equality is error-prone: instead, we should be checking if our next guess is smaller than our current guess by some fixed margin (set to 0.001 below):

def find_a_root(poly, initial_guess, num_iters):
    x = initial_guess
    delta = 0.001
    for i in range(1, num_iters + 1):
        # Get the next derivative
        p = poly.derive()
        result = poly.subst(x)
        new_result = p.subst(x)
        new_x = x - result / new_result
        if abs(new_x - x) < delta:
            # We can say that we've converged
            return new_x
        else:
            # Keep on updating x
            x = new_x
    return x
        

find_a_root(poly, initial_guess=10, num_iters=10)

1.000000000000038

The algorithm is pretty cool because even if our initial guess is very off, the updating mechanism quickly converges to some root of the equation:

find_a_root(poly, initial_guess=-123134583945824721131, num_iters=100)

-5.000000000001934