```python
import sys

import numpy as np
import matplotlib.pyplot as plt

sys.path.append('../src')
import utils
```

# Iterative Methods for Optimization

A general optimization problem has the form

\begin{align}
\min_x & \quad f_0(x) \\
\text{subject to} & \quad f_i(x) \leq b_i, \quad i = 1, \dots, m,
\end{align}

with components
* $x \in \mathbf{R}^n$ - optimization variable,
* $f_0: \mathbf{R}^n \to \mathbf{R}$ - objective function,
* $f_i: \mathbf{R}^n \to \mathbf{R}$ - constraint functions, and
* $b_i \in \mathbf{R}$: constraint bounds.

Very few optimization problems have a closed-form solution, so most problems are solved using iterative methods.
This notebook introduces gradient descent, a fundamental iterative method for solving unconstrained optimization problems of the form

$$\min_x f(x)$$

with $x \in \mathbf{R}^n$ and $f: \mathbf{R}^n \to \mathbf{R}$, where $f$ is differentiable. Our examples use [convex functions](https://en.wikipedia.org/wiki/Convex_function), but the intuition gained from the exercises can also be helpful for nonconvex functions.

Iterative methods start with an initial guess $x^0$.
Each iteration, a new value $x^{k+1}$ is generated based on the previous value $x^k$, often using local information such as the gradient $\nabla f(x^k)$ or Hessian $\nabla^2 f(x^k)$.
The algorithm terminates when a stopping condition is satisfied, which can be specified in terms of the size of the gradient, $\| \nabla f(x^k) \| < \epsilon$, the difference between successive iterates, $\| x^{k+1} - x^k \| < \epsilon$ or $\| f(x^{k+1}) - f(x^k) \| < \epsilon$, or a maximum number of iterations, $k > k_{\max}$.

# Gradient Descent

[Gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) is built upon the fact that the gradient of a function always points in the direction of maximum ascent, so the negative gradient points in the direction of maximum descent.
Therefore, starting at the point $x^k$ and taking a step in the direction of the negative gradient, we should arrive at a point $x^{k+1}$ that has a smaller function value, i.e., $f(x^{k+1}) < f(x^k)$.
After multiple iterations, we hope to approximate the optimal solution $x^* = \text{argmin}_x f(x)$.

Here's the algorithm for gradient descent with initial guess $x^0$, step size $\alpha$, and stopping condition based on the size of the gradient:

> **input** $x^0$, $\alpha$, $\epsilon$  
> **set** $k = 0$  
> **while** $\| \nabla f(x^k) \| \geq \epsilon$ **do**  
> &nbsp;&nbsp;&nbsp;&nbsp;$x^{k+1} \leftarrow x^k - \alpha \nabla f(x^k)$  
> &nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$  
> **end while**

##### How do we choose the initial guess?

Depending upon the problem, the initial guess can influence the value of the solution, the number of iterations needed to converge to a solution, or whether or not the algorithm converges at all.
Determining a "good" initial guess is a difficult problem itself, especially when $x$ is high-dimensional, so heuristics are often useful:
* $x^0$ = a vector of ones or zeros
* $x^0$ = a random vector
* $x^0$ = the solution to a similar, but easier, problem

##### How do we choose the step size?

The step size influences the rate of convergence.
If the step size is too small, it may take many more iterations to converge.
However, if the step size is too large, we may overshoot the optimal point.
If you are really unlucky, you could choose a step size that results in the iterates cycling through the same points forever.

In some cases, we can guarantee convergence based on the step size. For example, if a function is convex,

$$f( \lambda x + (1 - \lambda) y ) \leq \lambda f(x) + (1 - \lambda) f(y) \quad \text{for all} \quad x, y \in \mathbf{R}^n \quad \text{and} \quad 0 \leq \lambda \leq 1,$$

and continuously differentiable with a $\beta$-Lipschitz gradient for some $\beta > 0$,

$$\| \nabla f(x) - \nabla f(y) \| \leq \beta \| x - y\| \quad \text{for all} \quad x, y \in \mathbf{R}^n,$$

then after $k$ iterations with a step size $\alpha \leq 1/\beta$, we will have

$$f(x^k) - f(x^*) \leq \frac{\| x^0 - x^* \|^2}{2\alpha k}.$$

This means that gradient descent will converge at a rate of $\mathcal{O}(1/k)$.

As an example, consider the linear least-squares problem, where $f(x) = \frac{1}{2}\|Ax - b\|^2$.
The objective is convex and differentiable, with $\nabla f(x) = A^T(Ax - b)$.
To find the Lipschitz constant for the gradient, note that

$$\| \nabla f(x) - \nabla f(y) \| = \| A^TA(x - y) \| \leq \| A^TA \| \| x - y\|,$$

which follows from the definition of the induced matrix norm

$$\|A\| = \max_{x \neq 0} \frac{\|Ax\|}{\|x\|}.$$

Therefore, using a step size $\alpha \leq 1/\beta$, where $\beta = \|A^TA\|$, we are guaranteed to converge at a linear rate.

Other approaches for adaptively choosing step sizes include:
* [line search methods](https://en.wikipedia.org/wiki/Line_search): optimizing the step size along a given descent direction
* [trust region methods](https://en.wikipedia.org/wiki/Trust_region): optimizing the descent direction for a given step size

##### How do we choose the stopping condition?

For an unconstrained problem where $f$ is convex and differentiable, a point $x$ is optimal if and only if $\nabla f(x) = 0$, so it's natural to terminate gradient descent once the gradient becomes "small enough", i.e., $\| \nabla f(x) \| < \epsilon$.
For constrained problems, stopping conditions can similarly be derived from the [KKT conditions](https://en.wikipedia.org/wiki/Karush-Kuhn-Tucker_conditions), which are necessary, and in convex cases sufficient, conditions for optimality.

Another option is to stop iterating once the difference between successive iterates becomes small, i.e., $\|x^{k+1} - x^k\| < \epsilon$ or $\|f(x^{k+1}) - f(x^k)\| < \epsilon$. In the case of gradient descent, where the difference between iterates is influenced by the size of the gradient, these stopping conditions turn out to be similar to stopping once the gradient becomes small.

When it comes to choosing $\epsilon$, we need to balance convergence with the needs of the problem. On the one hand, if $\epsilon$ is too large, the algorithm will terminate before converging to a good solution. On the other hand, setting $\epsilon$ too small will cause the algorithm to take an unnecessarily long time. For example, if our data was recorded up to three decimal points, there would be no need to iterate until $\|x^{k+1} - x^k\| <$ 1e-6 when $\epsilon =$ 1e-3 should be sufficient. In addition to the precision of our data, another consideration for choosing $\epsilon$ is the signal-to-noise ratio in the data.

## Exercise 1: Step Size

In the example below, experiment with the step size to answer the following questions:
* What happens if the step size is too small? Try 0.001.
* What happens if the step size is too big? Try 1.001.
* Can you find a step size that converges in one step?
* Can you find a step size that results in a cycle?
* What is the Lipschitz constant of the gradient? What happens if you set the step size to $1/\beta$ or $1/\beta \pm 0.1$?


```python
step = 0.1 # modify this variable
func = lambda x: x**2
grad = lambda x: 2*x
results = utils.gradient_descent(func, grad, step=step, x0=5)
utils.plot_1d(func, results)
```

## Exercise 2: Initial Guess

In the example below, experiment with the initial guess to answer the following questions:
* Trying different values of x0, can you find a relationship between the shape of the function, the initial guess, and the number of iterations?
* Trying different values of $c \in (-2, 2)$, can you find a relationship between the shape of the function, the initial guess, and the number of iterations?

##### Condition Number

For $c \in (-2, 2)$, the function $f(x,y) = x^2 + cxy + y^2$ describes an ellipsoid with a constant Hessian $\nabla^2 f(x, y)$.
The [condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices) of the matrix $\nabla^2 f(x, y)$, the ratio of its largest eigenvalue to its smallest eigenvalue, gives a measure of the ellipsoid's eccentricity:
* A condition number near one means that the ellipsoid is nearly spherical, and different initial conditions along a sublevel set will converge in a similar number of iterations. To see this, try values of $c$ near zero and initial guesses $x^0 = (\pm 1, \pm 1)$.
* A large condition number indicates that the sublevel sets are much wider in some directions than in others, so different initial conditions along a sublevel set could take much longer to converge than others. To see this, try values of $c$ near 2 and initial guesses $x^0 = (\pm 1, \pm 1)$. 


```python
x0 = [0, 1]; c = 1 # modify these variables
func = lambda x: x[0]**2 + c*x[0]*x[1] + x[1]**2
grad = lambda x: np.array([2*x[0] + c*x[1], 2*x[1] + c*x[0]])
hess = [[2, c], [c, 2]]
print(f'Condition number of the Hessian: {np.linalg.cond(hess):.2f}')
results = utils.gradient_descent(func, grad, x0=x0, step=0.1)
utils.plot_2d(func, results)
```

# Stochastic Gradient Descent

For many machine learning problems with large data sets, computing the gradient can be computationally expensive, so an approximation of the gradient is used instead.
This approach is called [stochastic gradient descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent).

Consider an unconstrained optimization problem of the form

$$\min_x \frac{1}{m} \sum_{i=1}^m f_i(x)$$

with $x \in \mathbf{R}^n$ and $f_i: \mathbf{R}^n \to \mathbf{R}$, where each $f_i$ is differentiable. In many cases, each function $f_i$ corresponds to the prediction error for a single training example.
Rather than computing the full gradient $\nabla f(x) = \frac{1}{m} \sum_{i=1}^m f_i(x)$ every iteration, we can approximate it using a sample of the constituent functions.
While true stochastic gradient descent only samples one training example per iteration, it's typically more effective to use multiple examples, called a mini-batch. 

Here's the algorithm for stochastic gradient descent with initial guess $x^0$, step size $\alpha^0$ (often called the [learning rate](https://en.wikipedia.org/wiki/Learning_rate) in the machine learning community), and mini-batch size $b$:

> **input** $x^0$, $\alpha^0$, $b$  
> **set** $k = 0$  
> **while** not converged **do**  
> &nbsp;&nbsp;&nbsp;&nbsp;Sample $b$ training examples $\mathcal{B} \subset [1, m]$  
> &nbsp;&nbsp;&nbsp;&nbsp;$x^{k+1} \leftarrow x^k - \frac{\alpha^k}{b} \sum_{i \in \mathcal{B}} \nabla f_i(x^k)$  
> &nbsp;&nbsp;&nbsp;&nbsp;$\alpha^{k+1} \leftarrow \frac{\alpha^k}{1 + dk}$  
> &nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$  
> **end while**

##### Why does the learning rate change each iteration?

One problem with approximating the gradient is that you get a noisy gradient.
Depending on the training examples used, the approximation can vary in both size and direction from the true gradient.
(For this reason, it no longer makes sense to use the size of the gradient as our stopping condition.)
One way to reduce this variance is to use a larger mini-batch, but using too large a batch size reduces the computational efficiency gained by approximating the gradient in the first place.

The same problems we had before with the step size for gradient descent are now magnified by the noisy gradient: too large a learning rate can cause the problem to diverge, while too small a learning rate can can make it slow to converge.
A common approach to this problem is to use an adaptive learning rate, defined by a learning rate schedule.
The main idea of this approach is to start with a larger learning rate, then slowly decrease it as the model nears a solution.
While this doesn't do anything to decrease the variance within the gradient estimates, it prevents the algorithm from oscillating around a solution once it's found one. Here we implement a learning rate schedule based on the previous iterate, the iterate number, and a decay parameter.

## Exercise 3: Mini-batch Size

The following two examples solve the linear least-squares problem

$$\min_x \frac{1}{m} \sum_{i=1}^m (a_i^T x - y_i)^2,$$

where $a_i \in \mathbf{R}^{784}$, $x \in \mathbf{R}^{784}$, $y_i \in \mathbf{R}$, and $m = 1000$.
Here each $a_i$ is a vector of data features, each $y_i$ is a data value, and $x$ is a vector of model parameters.
For illustrative purposes, we plot both the function value at each iteration, but it's worth noting that in practice, these values may not be available.
For large data sets, computing the function value can also be expensive, so estimates based on a subset of the observations or a test set can also be used, with the full function value being computed only periodically.

In the example below, experiment with the mini-batch size and maximum number of iterations to answer the following questions:
* Let max_iter = 1500. Using batch sizes 1 (true stochastic gradient descent), 10, 100, and 1000 (gradient descent), fill out the table below. Can you find a relationship between the batch size, minimum function value, and total computation time?

Batch Size | Function Value | Iterations | Total Time 
-|-|-|-
1 |  |  |
10 |  |  |
100 |  |  |
1000 |  |  |

* Let max_iter = 7500. Using batch sizes 1 (true stochastic gradient descent), 10, 100, and 1000 (gradient descent), fill out the table below. Can you find a relationship between the batch size, minimum function value, number of iterations, and total computation time?

Batch Size | Function Value | Iterations | Total Time 
-|-|-|-
1 |  |  | 
10 |  |  | 
100 |  |  | 
1000 |  |  | 

* Which of the batch sizes give a good compromise between the minimum function value and the total computation time?


```python
A = np.load('../data/train_images.npy')
y = np.load('../data/train_labels.npy')
```


```python
batch_size = 1; max_iter = 1500 # modify these variables
results = utils.stochastic_descent(A, y, batch_size=batch_size, max_iter=max_iter, learning_rate=1e-2)
utils.plot_sgd(results)
```

## Exercise 4: Decay Rate

In the example below, experiment with the decay parameter to answer the following questions:
* Using decay parameters None (constant learning rate), 1e-7, 5e-7, and 1e-6, fill out the table below. Can you find a relationship between the decay parameter, minimum function value, and total computation time?

Decay Parameter | Function Value | Iterations | Total Time 
-|-|-|-
None |  |  |
1e-7 |  |  |
5e-7 |  |  |
1e-6 |  |  |

* Which of the decay parameters give a good compromise between the minimum function value and the total computation time?


```python
decay = None # modify this variable
results = utils.stochastic_descent(A, y, decay=decay, batch_size=1, max_iter=7500, learning_rate=1e-2)
utils.plot_sgd(results)
```

# Further Reading

* Kelley, C. T. (1999). *Iterative methods for optimization*. Society for Industrial and Applied Mathematics.
* Boyd, S., Boyd, S. P., & Vandenberghe, L. (2004). *Convex optimization*. Cambridge University Press.
* Nocedal, J., & Wright, S. (2006). *Numerical optimization.* Springer Science & Business Media.
* Beck, A. (2014). *Introduction to nonlinear optimization: Theory, algorithms, and applications with MATLAB*. Society for Industrial and Applied Mathematics.

7fdc1900d09a8e2acdb8a4aff5f396cb58fe7185
ipynb
Jupyter Notebook
th-optimization-adsi/notebooks/Part1-IterativeMethods.ipynb
bastivkl/nh2020-curriculum
245a72af3f325495448cbf6c0c6baa2499d43d94
CC-BY-4.0
2020-06-27T19:04:11.000Z
2022-03-28T00:44:44.000Z
th-optimization-adsi/notebooks/Part1-IterativeMethods.ipynb
bastivkl/nh2020-curriculum
245a72af3f325495448cbf6c0c6baa2499d43d94
CC-BY-4.0
2020-07-23T02:11:40.000Z
2020-09-09T21:28:36.000Z
th-optimization-adsi/notebooks/Part1-IterativeMethods.ipynb
bastivkl/nh2020-curriculum
245a72af3f325495448cbf6c0c6baa2499d43d94
CC-BY-4.0
2020-07-15T03:37:49.000Z
2022-02-27T23:07:14.000Z
Qwen/Qwen-72B
1. YES
2. YES