Optimization¶

In The Mean and Slopes, we used a simple but slow way to find the slope that best predicted one vector of values from another vector of values.

First we go back to find that slope.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Make plots look a little bit more fancy
plt.style.use('fivethirtyeight')
# Print to 2 decimal places, show tiny values as 0
np.set_printoptions(precision=2, suppress=True)
import pandas as pd

Download the ckd_clean.csv file to the same directory as this notebook, if you are running on your own computer.

We fetch and process the data. See mean and slopes for a slower description of this processing.

# Load the data file
ckd = pd.read_csv('ckd_clean.csv')
pcv = np.array(ckd['Packed Cell Volume'])
hgb = np.array(ckd['Hemoglobin'])

Our criterion is the sum of squared error:

def sos_error(slope):
    fitted = hgb * slope  # 'hgb' comes from the top level
    error = pcv - fitted     # 'pcv' comes from the top level
    return np.sum(error ** 2)

We found the best slope by trying a very large number of slopes, and recording, for each slope, the sum of squared error. We chose the slope from the slopes that we tried, that gave us the lowest sum of squared error.

# Slopes to try
some_slopes = np.arange(2, 4, 0.01)
n_slopes = len(some_slopes)
# Try all these slopes, calculate and record sum of squared error
sos_errors = np.zeros(n_slopes)
for i in np.arange(n_slopes):
    slope = some_slopes[i]
    sos_errors[i] = sos_error(slope)
# The slope minimizing the sum of squared error
best_slope = some_slopes[np.argmin(sos_errors)]
best_slope

3.0499999999999776

At the end, of the mean and slopes notebook, we saw that a function in Scipy called minimize can do this work for us, relatively quickly.

from scipy.optimize import minimize
minimize(sos_error, 3)

      fun: 3619.6157878183294
 hess_inv: array([[0.]])
      jac: array([0.])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 10
      nit: 2
     njev: 5
   status: 2
  success: False
        x: array([3.05])

What we are doing, with our slow dumb technique, and with the minimize function, is something called mathematical optimization. We use optimization when we have some function that takes one or more parameters. We want to chose, or optimize the parameters to give us some desired output from the function.

In our case our function is the sum of squared error, sos_error. The parameter is the slope. We are trying to find the value for the parameter that minimizes the result of calling the function sos_error.

One way of doing this minimization, is the slow dumb way. We just try a huge number of values for the parameter (the slope), and chose the value that gives us the lowest output value (the sum of squared error).

This is such a common problem, that there has been an enormous amount of theoretical and practical work on building algorithms to make process of searching for the minimum value more efficient.

This notebook is to give you an idea of how you might do this, and therefore, the kind of things that minimize can do, to search quickly for the best parameter.

Let’s look again at the shape of the curve relating the slope to the sum of squared error:

plt.plot(some_slopes, sos_errors)
plt.xlabel('Candidate slopes')
plt.ylabel('Sum of squared error')
plt.title('SSE as a function of slope')

Text(0.5, 1.0, 'SSE as a function of slope')

This is the function we are trying minimize. Specifically, we are trying to optimize the function that gives the SSE as a function of the slope parameter.

We want to avoid trying every possible value for the slope.

To do this, we are going to start with one value for the slope, say 100, then see if there is a good way to chose the next value to try.

Looking at the graph, we see that, when the slope is far away from the minimum, the sum of squared error (on the y axis) changes very quickly as the slope changes. That is, the function has a steep gradient.

Maybe we could check what the gradient is, at our starting value of 100, by calculating the sum of squares (y) value, and then calculating the sum of squares (y) value when we increase the slope by a tiny amount. This is the change in y for a very small change in x. We divide the change in y by the change in x, to get the gradient:

def sos_error_gradient(x, dx=0.0001):
    # Gradient of the sos_error at this value of x
    # sos_error at this x value.
    sos_0 = sos_error(x)
    # sos_error a tiny bit to the right on the x axis.
    sos_1 = sos_error(x + dx)
    # gradient is y difference divided by x difference.
    return (sos_1 - sos_0) / dx

# The y value of the function.
sos_error(4)

31657.96

# The gradient of the function at this point.
sos_error_gradient(4)

58876.17044587387

A large positive gradient means the x value (slope) that we tried is still far to the right of the minimum. This might encourage us to try an x value that is well to the left. We could call this a large step in x, and therefore a large step size.

Let’s try another value:

# The y value of the function.
sos_error(2)

37529.64

# The gradient of the function at this point.
sos_error_gradient(2)

-64741.669554132386

A large negative gradient means the x value (slope) that we tried is still far to the left of the minimum. This might encourage us to try an x value that is well to the right.

As the gradients get small, we want to take smaller steps, so we don’t miss the minimum.

The general idea then, is to chose our step sizes in proportion to the gradient of the function.

This is the optimization technique known as gradient descent.

Here it is in action. We try new x-axis values by making big jumps when the gradient is steep, and small jumps when the gradient is shallow. For each new value of x, we check the gradient, and choose our new jump size.

Remember each value on the x-axis is another value we want to try, in order to find the best-fitting slope for our original problem.

next_x = 4 # We start the search at x=4
gamma = 0.00001 # Step size multiplier
precision = 0.00001 # Desired precision of result
max_iters = 1000 # Maximum number of iterations

for i in np.arange(max_iters):
    # Go to the next x value
    current_x = next_x
    # Estimate the gradient
    gradient = sos_error_gradient(current_x)
    # Use gradient to choose the next x value to try.
    # This takes negative steps when the gradient is positive
    # and positive steps when the gradient is negative.
    next_x = current_x - gamma * gradient
    step = next_x - current_x
    print('x: {:0.5f}; step {:0.5f}; gradient {:0.2f}'.format(
        current_x, step, gradient))
    # When the step size is equal to or smaller than the desired
    # precision, we are near enough.
    if abs(step) <= precision:
        # Break out of the loop.
        break

print("Minimum at", next_x)

x: 4.00000; step -0.58876; gradient 58876.17
x: 3.41124; step -0.22485; gradient 22485.45
x: 3.18638; step -0.08587; gradient 8587.43
x: 3.10051; step -0.03280; gradient 3279.63
x: 3.06771; step -0.01253; gradient 1252.53
x: 3.05519; step -0.00478; gradient 478.35
x: 3.05040; step -0.00183; gradient 182.69
x: 3.04858; step -0.00070; gradient 69.77
x: 3.04788; step -0.00027; gradient 26.65
x: 3.04761; step -0.00010; gradient 10.18
x: 3.04751; step -0.00004; gradient 3.89
x: 3.04747; step -0.00001; gradient 1.48
x: 3.04746; step -0.00001; gradient 0.57
Minimum at 3.0474521484512476

As you can see, by doing this, we save ourselves from trying a very large number of other potential x values, and we focus in on the minimum very quickly.

This is just one method among many for optimizing our search for the minimum of a function. Now you know what kind of thing it is doing, we will just let miminize do its job for us:

# Use float to show us the final result in higher precision
result = minimize(sos_error, 100)
float(result.x)

3.0474986373317794

Coding for Data - 2020 edition

Optimization¶