# Gradient descent

## Simple linear regression

In a previous notebook, we solved the problem of simple linear regression - finding a straight line that best fits a data set with two variables.
In that case, we were able to find the exact solution.
In this notebook, we'll use a common technique to approximate that solution.

Why would we want to approximate a solution when we can easily find an exact solution?
We don't - it's just that the technique we discuss here can also be used in situations where we can't find an exact solution or don't want to try, for whatever reason.
The approximation technique is called gradient descent.

### Preliminaries


```python
# numpy efficiently deals with numerical multi-dimensional arrays.
import numpy as np

# matplotlib is a plotting library, and pyplot is its easy-to-use module.
import matplotlib.pyplot as pl

# This just sets the default plot size to be bigger.
pl.rcParams['figure.figsize'] = (16.0, 8.0)
```

### Simple linear regression model
In simple linear regression, we have some data points $(x_i, y_i)$, and we decide that they belong to a straight line with a little bit of error involved.
Straight lines in two dimensions are of the form $y = mx + c$, and to fit a line to our data points we must find appropriate values for $m$ and $c$.
Numpy has a function called `polyfit` that finds such values for us.


```python
w = np.arange(1.0, 16.0, 1.0)
d = 5.0 * w + 10.0 + np.random.normal(0.0, 5.0, w.size)

m, c = np.polyfit(w, d, 1)
print("Best fit is m = %f and c = %f" % (m, c))

# Plot the best fit line.
pl.plot(w, d, 'k.', label='Original data')
pl.plot(w, m * w + c, 'b-', label='Best fit: $%0.1f x + %0.1f$' % (m,c))
pl.legend()
pl.show()
```

### Gradient descent
In gradient descent, we select a random guess of a parameter and iteratively improve that guess.
For instance, we might pick $1.0$ as our initial guess for $m$ and then create a `for` loop to iteratively improve the value of $m$.
The way we improve $m$ is to first take the partial derivative of our cost function with respect to $m$.

### Cost function
Recall that our cost function for simple linear regression is:

$$
Cost(m, c) = \sum_i (y_i - mx_i - c)^2
$$

### Calculate the partial derivatives
We calculate the partial derivative of $Cost$ with respect to $m$ while treating $c$ as a constant.
Note that the $x_i$ and $y_i$ values are all just constants.
We'll also calculate the partial derivative with respect to $c$ here.

$$
\begin{align}
Cost(m, c)                       &= \sum_i (y_i - mx_i - c)^2      \\[1cm]
\frac{\partial Cost}{\partial m} &= \sum 2(y_i - m x_i -c)(-x_i) \\
                                 &= -2 \sum x_i (y_i - m x_i -c) \\[0.5cm]
\frac{\partial Cost}{\partial c} & = \sum 2(y_i - m x_i -c)(-1) \\
                                 & = -2 \sum (y_i - m x_i -c)   \\
\end{align}
$$

### Code the partial derivatives
Once we've calculated the partial derivatives, we'll code them up in python.
Here we create two functions, each taking four parameters.
The first two parameters are arrays with our $x_i$ and $y_i$ data set values.
The second two are our current guesses for $m$ and $c$.


```python
def grad_m(x, y, m, c):
  return -2.0 * np.sum(x * (y - m * x - c))
```


```python
def grad_c(x, y, m , c):
  return -2.0 * np.sum(y - m * x - c)
```

### Iterate
Now we can run our gradient descent algorithm.
For $m$, we keep replacing its value with $m - \eta grad\_m(x, y, m, c)$ until it doesn't change.
For $c$, we keep replacing its value with $c - \eta grad\_c(x, y, m, c)$ until it doesn't change.

What is $\eta$? It is called the learning rate and we set it to a small value relative to the data points.

You can see on each iteration, $m$ and $c$ are getting closer to their true values.


```python
eta = 0.0001
m, c = 1.0, 1.0
delta = 0.0000001

mold, cold = m - 1.0, c - 1.0
i = 0
while abs(mold - m) > delta and abs(cold - c) > delta:
  mold, cold = m, c
  m = mold - eta * grad_m(w, d, mold, cold)
  c = cold - eta * grad_c(w, d, mold, cold)

  i = i + 1
  if i % 1000 == 0:
    print("m: %20.16f  c: %20.16f" % (m, c))
```

    m:   5.7149506777702488  c:   4.5518695120877757
    m:   5.5629397868829846  c:   6.1183980049248472
    m:   5.4852588608786546  c:   6.9189286986499470
    m:   5.4455621902673004  c:   7.3280175818245992
    m:   5.4252763139437015  c:   7.5370710455442040
    m:   5.4149097826482908  c:   7.6439019874021357
    m:   5.4096122559659570  c:   7.6984949605832087
    m:   5.4069051027099357  c:   7.7263931768575951
    m:   5.4055216875445105  c:   7.7406497821666767
    m:   5.4048147318094859  c:   7.7479352226744274
    m:   5.4044534617990276  c:   7.7516582438453998
    m:   5.4042688448364773  c:   7.7535607899748742


## Newton's method for square roots
Newton's method for square roots is a method for approximating the square root of a number $x$.
We begin with an initial guess $z_0$ of the square root - it doesn't have to be particularly good.
We then apply the following calculation repeatedly, to calculate $z_1$, then $z_2$, and so on:

$$
z_{i+1} = z_i - \frac{z_i^2 - x}{2z_i}
$$

### Coding the calculation
We can create a function that calculates the next value of $z$ based on the current value of $z$ as follows.


```python
def next_z(x, z):
  return z - (z**2 - x) / (2 * z)
```

### Calculating the square root of $x$
Suppose we want to calculate the square root of $x$.
We start with a random guess for the square root, $z_0$.
We then apply the `next_z` function repeatedly until the value of $z$ stops changing.
Let's create a function to do this.
We'll include the next_z function inside the `newtsqrt` function to make it one all-inclusive package.


```python
def newtsqrt(x):
  next_z = lambda x, z: z - (z**2 - x) / (2 * z)
  z = 2.0
  n = next_z(x, z)
  
  while z != n:
    z, n = n, next_z(x, n)
    print(z)
  
  return z

newtsqrt(11)
```

    3.75
    3.341666666666667
    3.316718620116376
    3.3166247916826186
    3.3166247903554





    3.3166247903554



### Comparison with the standard library
We can compare our square root method return value to the value calculated by Python's `math` standard library package. It has a `sqrt` function.


```python
import math
math.sqrt(11)
```




    3.3166247903554



### Being careful
Due to the complexities of floating point numbers, the `nextsqrt` function could get into an infinite loop.
For instance, calculating the square root of 10 gives an infinite loop.


```python
# Uncommenting will result in infinite loop.
# newtsqrt(10)
```

To counteract this problem, the condition of the loop is better written as:
```python
abs(z - n) > 0.001
```

## Gradient descent for square roots
Newton's method for square roots is efficient, but we can also use gradient descent to approximate the square root of a real number $x$.
Here, we use the following cost function.

$$
Cost(z \mid x) = (x - z^2)^2
$$

### Example value
Let's use it to calculate the square root of 20, i.e. $x = 20$.
Then the cost function is:

$$
Cost(z \mid x=20) = (20 - z^2)^2
$$

Our goal is to find the $z$ that minimises this.

### Plotting the cost function
Let's plot the cost function.
Given that we know the best $z$ will be between $4$ and $5$, we'll let $z$ range over 0 to 10.


```python
i = np.linspace(0.0, 10.0, 1000)
j = (20.0 - i**2)**2

pl.plot(i, j, 'k-', label='$(20-z^2)^2$')
pl.legend()
pl.show()
```

### The derivative
Looks like there's a low point at about $4.5$.
Let's take the derivative of the cost function with respect to $z$.

$$
\begin{align}
Cost(z) &= ( 20.0 - z^2 )^2 \\
\Rightarrow \frac{\partial Cost}{\partial z} &= 2(20.0 - z^2)(-2z) \\
                                              &= 4z(z^2 - 20) \\
                                              &= 4z^3 - 80z
\end{align}
$$

The derivative tells us what the slope of the tangent to the curve is at any point on the cost function.
What does that mean?
It means that if we pick a value of $z$, e.g. $8.0$, that the derivative tells us that a line going through the point $(8.0, (20.0 - (8.0)^2)^2)$ with the slope $4(8.0)^3 - 80(8.0)$ perfectly touches the graph above.
Let's plot that.

When you simply, the point $(8.0, (20.0 - (8.0)^2)^2)$ becomes $(8,1936)$.
The slope is $4(8.0)^3 - 80(8.0)$ which when simplified becomes $1408$.
So, the claim is that the line with slope $1408$ going through the point $(8,1936)$ touches the graph.
To calculate the equation of the line, we'll use $(y - y_1) = m(x - x_1)$:

$$
y - 1936 = 1408(x - 8) \\
\Rightarrow y = 1408x - 11264 + 1936 \\
\Rightarrow y = 1408x - 9328
$$

Let's plot that line and the cost function together.


```python
i = np.linspace(0.0, 10.0, 1000)
j = (20.0 - i**2)**2
k = 1408 * i - 9328

pl.plot(i, j, 'k-', label='$(20-z^2)^2$')
pl.plot(i, k, 'b-', label='$1408z - 9328$')
pl.legend()
pl.show()
```

### Why do we care about the slope?
It's a bit hard to see, but the blue line is perfectly touching the curve.
We care about this because the slope of the blue line tells us in which way to change $z$ in order to make the cost less.
If we increase $z$ the cost goes up.
If we decrease it the cost goes down.

### Gradient descent
Let's use the gradient descent algorithm to calculate the best $z$.
We'll start with the guess $z=8$, and then use the derivative to move $z$ ever so slightly in the direction that decreases the cost.
By ever so slightly, we mean $0.001$ times the slope:

$$
\begin{align}
z_{i+1} &= z_i - \eta \frac{\partial Cost}{\partial z} \\
        &= z_i - (0.001) (4 z_i^3 - 80 z_i)
\end{align}
$$

So, for our initial guess $z_0 = 8.0$ we get:

$$
\begin{align}
z_1 &= 8.0 - (0.001) (4 (8.0)^3 - 80 (8.0)) \\
    &= 8.0 - 1.408 = 6.592
\end{align}
$$

Let's code it up.


```python
def next_z(z, x, eta=0.001):
  return z - eta * (4.0 * z**3 - 80 * z)

def sqrt_grad_desc(x, z, verbose=False):
  while abs(z - next_z(z, x)) > 0.001:
    if verbose:
      print("Current: %14.8f\tNext: %14.8f" % (z, next_z(z, x)))
    z = next_z(z, x)
  return z

ans =sqrt_grad_desc(20.0, 8.0, True)
print("Square root:", ans, "\tSquared:", ans**2)
```

    Current:     8.00000000	Next:     6.59200000
    Current:     6.59200000	Next:     5.97355269
    Current:     5.97355269	Next:     5.59881186
    Current:     5.59881186	Next:     5.34469983
    Current:     5.34469983	Next:     5.16157297
    Current:     5.16157297	Next:     5.02444369
    Current:     5.02444369	Next:     4.91903017
    Current:     4.91903017	Next:     4.83645229
    Current:     4.83645229	Next:     4.77084541
    Current:     4.77084541	Next:     4.71815685
    Current:     4.71815685	Next:     4.67548576
    Current:     4.67548576	Next:     4.64069702
    Current:     4.64069702	Next:     4.61218330
    Current:     4.61218330	Next:     4.58871218
    Current:     4.58871218	Next:     4.56932433
    Current:     4.56932433	Next:     4.55326362
    Current:     4.55326362	Next:     4.53992784
    Current:     4.53992784	Next:     4.52883326
    Current:     4.52883326	Next:     4.51958845
    Current:     4.51958845	Next:     4.51187478
    Current:     4.51187478	Next:     4.50543157
    Current:     4.50543157	Next:     4.50004463
    Current:     4.50004463	Next:     4.49553736
    Current:     4.49553736	Next:     4.49176369
    Current:     4.49176369	Next:     4.48860255
    Current:     4.48860255	Next:     4.48595333
    Current:     4.48595333	Next:     4.48373229
    Current:     4.48373229	Next:     4.48186965
    Current:     4.48186965	Next:     4.48030717
    Current:     4.48030717	Next:     4.47899619
    Current:     4.47899619	Next:     4.47789603
    Square root: 4.477896027970839 	Squared: 20.05155283731702


### A question
Let's try some other initial guesses: 4.0, 1.0 and -1.0.
Can you explain the square root returned with -1.0?


```python
print("With initial guess %6.2f: %10.6f" % (4.0, sqrt_grad_desc(20.0, 4.0, False)))
print("With initial guess %6.2f: %10.6f" % (1.0, sqrt_grad_desc(20.0, 1.0, False)))
print("With initial guess %6.2f: %10.6f" % (-1.0, sqrt_grad_desc(20.0, -1.0, False)))
```

    With initial guess   4.00:   4.465951
    With initial guess   1.00:   4.465963
    With initial guess  -1.00:  -4.465963


### End

f0c2d20fc1537019ec6b3e812093d08cea42499e
ipynb
Jupyter Notebook
gradient-descent.ipynb
sean-meade/jupyter-teaching-notebooks
2f2a2bf6925fba479d1a73481122faebd201bdba
Unlicense
gradient-descent.ipynb
sean-meade/jupyter-teaching-notebooks
2f2a2bf6925fba479d1a73481122faebd201bdba
Unlicense
gradient-descent.ipynb
sean-meade/jupyter-teaching-notebooks
2f2a2bf6925fba479d1a73481122faebd201bdba
Unlicense
Qwen/Qwen-72B
1. YES
2. YES