# Homework 05: Conjugate Gradient Descent

In this homework, we will implement the conjugate graident descent algorithm. While you should nearly always use an optimization routine from a library for practical data analyiss, this exercise is useful because it will make concepts from multivariatble calculus and linear algebra covered in the lectrures concrete for you. Also, it brings you up the learning curve for the implementaiton of more complex algorithms than the ones you have been exposed to so far.

Note: The exercise assumes that we can calculate the gradient and Hessian of the fucntion we are trying to minimize. This can be computationally expensive or not even possible for soeme functions. Approximate methods can then be used; we do not go into such complexities here.

**Conjugate gradient descent**

We want to implement the line search method 

$$ 
x_{k+1} = x_k + \alpha_k p_k
$$

where $\alpha_k$ is the step size and $p_k$ is the search direction. 

In particular, we want the search directions $p_k$ to be *conjugate*, as this will allow us to find the minimum in $n$ steps for $x \in \mathbb{R}^n$ if $f(x)$ is a quadratic function.

The following exercises will unpack this:

- What quadratic functions are
- What conjugate vectors are
- How to find conjugate vectors by Gram-Schmidt process
- How to find the step size $\alpha_k$

and finally wrap them all into a conjugate gradient algorithm.

**Quadratic function surfaces**

Recall that our objective is to minimize a scalar valued function which maps $\mathbb{R}^n \mapsto \mathbb{R}$, for example, a log likelihoood function (for MLE) or unnormalized posterior distribution (for MAP). Geometrically, we are tring to find the value of the lowest point of some surface. The conjugate gradient algorihtm assumes that the surface can be approximated by the quadratic expression (say, by using a Taylor series expansion about $x$)

$$
f(x) = \frac{1}{2}x^TAx - b^Tx + c
$$

and that 

$$
\nabla f = Ax - b = 0
$$

at the minimum (if A is positive definite). Note that $A$ is a matrix, $b$ is a vector, and $c$ is a scalar.  Also, note that the matrix $A$ is the Hessian of the quadratic function.For simplicity, we'll work in $\mathbb{R}^2$ so we can visualize the surface, so that $x$ is a 2-vector.

Note: A **form** is a polynomial function where every term has the same degree - for example, $x^2 + 2xy + y^2$ is a quadratic form, whcih can be rewritten as 
$$
\begin{pmatrix}
x & y
\end{pmatrix}
\begin{pmatrix}
    1 & 1\\
    1 & 1
\end{pmatrix}
\begin{pmatrix}
x \\
y
\end{pmatrix}
$$

That is, $x^TAx$ is a quadratic form.


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg as la
```

**Exercise 1 (20 points)** 

We will work with function $f_1$

$$
f1(x) = \frac{1}{2} x^T \pmatrix{1 & 0 \\ 0 & 1}x
$$

and function $f_2$

$$
f2(x) = \frac{1}{2} x^T \pmatrix{1 & 0 \\ 0 & 3}x
$$

- Plot the labeled contours of the quadratic functions
- Use a `streamplot` to show the gradient vector field of the above quadratic functions.


```python
def f1(x, y):
    '''function f1'''

    return 0.5 * (x**2 + y**2)

def f2(x, y):
    '''fucntion f2'''
    
    return 0.5 * (x**2 + 3 * y**2)

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
plt.contour(X, Y, f1(X, Y), 10)
plt.contour(X, Y, f2(X, Y), 10)
plt.axis('square')
pass
```


```python
def u1(x, y):
    '''df1/px'''
    
    return x

def v1(x, y):
    '''df1/dy'''
    
    return y
    
def u2(x, y):
    '''df2/dx'''
    
    return x
    
def v2(x, y):
    '''df2/dy'''
    
    return 3 * y

plt.streamplot(X, Y, u1(X, Y), v1(X, Y), density = 1)
plt.streamplot(X, Y, u2(X, Y), v2(X, Y), density = 1)
pass
```

**Gram-Schmidt**

The way to numerically find conjugate vectors is to use the Gram-Schmidt process. Here, instead of the usual projection 

$$
\text{proj}_u(v) = \frac{u \cdot v}{u \cdot u} \, u
$$

we use the generalized projection
$$
\text{proj}_u(v) = \frac{uA^Tv}{uA^Tu} \, u
$$

**Exercise 2 (30 points)**

The vectors $u$ and $v$ are orthogonal i.e. $u^Tv = 0$ and conjugate with respect to $A$ if $u^TAv = 0$. The geometric intuition for conjugate vecrors $u$ and $v$ is that $u$ and $v$ would be orthogonal if we stretched the contour plots so that it became isotropic (same in all directions, just like when A=𝟙). Write a Gram-Schmidt function to find orthogonal and conjuate vectors with the following signature

```python
def gram_schmidt(U, inner):
    """Return an orthogonal matrix.
    
    U is a matrix of (column) vectors.
    inner is a function that calculates the inner product.
    
    Returns an orthogonal matrix of the same shape as U.
    """
```

Use this function and the appropiate inner product to plot

- An orhtogonal set of basis vectors for $f_1$
- A conjugate set of basic vectors for $f_2$ 

where the first basis vector is to parallel to $\pmatrix{1 \\ 1}$.


```python
def inner_f1(x, y):
    
    A = np.array([[1,0],[0,1]])
    return (x.T @ A.T @ y).reshape(-1)

def inner_f2(x, y):
    
    A = np.array([[1,0],[0,3]])
    return (x.T @ A.T @ y).reshape(-1)

def gram_schmidt(U, inner):
    """Return an orthogonal matrix.

    U is a matrix of (column) vectors.
    inner is a function that calculates the inner product.

    Returns an orthogonal matrix of the same shape as U.
    """
    
    V = np.zeros((U.shape[0], U.shape[1]), dtype = 'float64')
    V[:,0] = U[:,0]
    
    for i in range(1, U.shape[1]):
        V[:,i] = U[:,i]
        
        for j in range(0, i):
            V[:,i] -= (inner(V[:,i], V[:,j]) / inner(V[:,j], V[:,j])) * V[:,j]
    
    return V
```


```python
U = np.array([[1,0],[1,1]])
basis_f1 = gram_schmidt(U, inner_f1)
basis_f2 = gram_schmidt(U, inner_f2)
print(basis_f1)
print(basis_f2)

plt.figure(figsize=(4,4))
plt.quiver(0, 0, basis_f1[0,0],basis_f1[1,0],color='red',angles='xy',scale_units='xy',scale=2);
plt.quiver(0, 0, basis_f1[0,1],basis_f1[1,1],color='red',angles='xy',scale_units='xy',scale=2);
plt.quiver(0, 0, basis_f2[0,0],basis_f2[1,0],color='blue',angles='xy',scale_units='xy',scale=2);
plt.quiver(0, 0, basis_f2[0,1],basis_f2[1,1],color='blue',angles='xy',scale_units='xy',scale=2);
```

**Exercise 3 (20 points)**

We now need to find the "step size" $\alpha$ to take in the direction of the search vector $p$. We can get a quadratic approximation to a general nonliner function $f$ by taking the Taylor series in the driection of $p$

$$
f(x + \alpha p) = f(x) + \alpha [f'(x)]^T p + \frac{\alpha^2}{2} p^T f''(x) p
$$

Find the derivative with respect to $\alpha$ and use this to find the optimal value for $\alpha$ with respect to the quadratic approcimaiton. 

* Write the derivation of the optimal step size $\alpha$ using LaTeX in a Markdown cell.
* Write a funciton that returns $\alpha$ for a quadratic funciton with the following signature

```python
def step(x, p, A, b):
    """Returns the optimal step size to take in line search on a quadratic.
    
    A and b are the coefficients of the quadartic expression 
    $$
    f(x) = \frac{1}{2}x^TAx - b^Tx + c
    $$
    p is the search direction
    x is the current location
    """
```

$$
\begin{align}
f(x+\alpha p)&=f(x)+\alpha[f'(x)]^Tp+\frac{\alpha^2}{2}p^Tf''(x)p\\
\frac{df(x+\alpha p)}{d\alpha}&=[f'(x)]^Tp+\alpha p^Tf''(x)p\\
\end{align}
$$
Solve for $\frac{df(x+\alpha p)}{d\alpha}=0$, we get
$$
\begin{align}
\alpha&=-\frac{[f'(x)]^Tp}{p^Tf''(x)p}\\
&=-\frac{\nabla f^T p}{p^T A p}\\
&=\frac{(b - Ax)^T p}{p^T A p}
\end{align}
$$





```python
def step(x, p, A, b):
    """Returns the optimal step size to take in line search on a quadratic.

    A and b are the coefficients of the quadartic expression 
    $$
    f(x) = \frac{1}{2}x^TAx - b^Tx + c
    $$
    p is the search direction
    x is the current location
    """
    
    alpha = ((b - A @ x).T @ p) / (p.T @ A @ p)
    
    return alpha
```

**Line search**

We now know how to find a search direction $p_k$ - this is a vector that is conjugate to the previous search direction. The first search direction is usually set to be the gradient. Next we need to find out how far along $p_k$ we need to travel, i.e., we need to find $\alpha_k$. First we take a Taylor expansion in the direction of $p$

$$
f(x + \alpha p) = f(x) + \alpha [f'(x)]^T p + \frac{\alpha^2}{2} p^T f''(x) p
$$

followed by finding the derivative with respect to $\alpha$

$$
\frac{d}{d\alpha} f(x + \alpha p) = [f'(x)]^T p + \alpha p^T f''(x) p
$$

Solvign for $\frac{d}{d\alpha} f(x + \alpha p) = 0$, we get

$$
\alpha = - \frac{[f'(x)]^T p}{p^T f''(x) p} \\
= - \frac{\nabla f^T p}{p^T A p} \\
= \frac{(b - Ax)^T p}{p^T A p}
$$



**Exercise 4 (30 points)**

Implement the conjugate grdient descent algorithm with the following signature

```python
def cg(x, A, b, c, max_iter=100, tol=1e-3):
    """Conjugate gradient descent on a quadratic function surface.
    
    x is the starting position
    A, b and c are the coefficients of the quadartic expression 
    $$
    f(x) = \frac{1}{2}x^TAx - b^Tx + c
    $$
    max_iter is the maximum number of iterations to take
    tol is the tolerance (stop if the length of the gradient is smaller than tol)

    Returns the number of steps taken and the list of all positions visited.
    """
```

- Use cg to find the minimum of the funciton $f_2$ from Exercise 1, starting from $\pmatrix{6 \\ 7}$.

- Plot the contour of the funciton f and the trajectory taken from the inital starting poitn $x$ to the final position, inlcuding all the intermediate steps. 

We are not particularly concerned about efficiency here, so don't worry about JIT/AOT/C++ level optimization. 


```python
def cg(x, A, b, c, max_iter=100, tol=1e-3):
    """Conjugate gradient descent on a quadratic function surface.

    x is the starting position
    A, b and c are the coefficients of the quadartic expression 
    $$
    f(x) = \frac{1}{2}x^TAx - b^Tx + c
    $$
    max_iter is the maximum number of iterations to take
    tol is the tolerance (stop if the length of the gradient is smaller than tol)

    Returns the number of steps taken and the list of all positions visited.
    """
    
    def f(x):
        '''return f(x)'''
        
        return 0.5 * x.T @ A @ x - b.T @ x + c 
    
    orbit = np.zeros((max_iter + 1, len(x)))
    p = np.zeros((max_iter + 1, len(x)))
    r = np.zeros((max_iter + 1, len(x)))
    
    orbit[0] = x.ravel()
    r[0] = b - A @ orbit[0]
    if la.norm(r[0]) < tol:
        return f(orbit[i])
    
    p[0] = r[0]
    for i in range(max_iter):
        
        alpha = step(orbit[i], p[i], A, b)
        orbit[i+1] = (orbit[i] + alpha * p[i]).ravel()
        r[i+1] = r[i] - alpha * (A @ p[i])
        
        if la.norm(r[i+1]) < tol:
            break
        
        beta = (r[i+1].T @ r[i+1]) / (r[i].T @ r[i])
        p[i+1] = r[i+1] + beta * p[i]
    
    return orbit[:i+2,], f(orbit[i+1])
```


```python
x0 = np.array([6, 7])
A = np.array([[1, 0], [0, 3]])
b = np.array([0,0])
c = 0.0

xs, minf = cg(x0, A, b, c)
print(xs)
print(minf)
```

    [[ 6.00000000e+00  7.00000000e+00]
     [ 3.89403974e+00 -3.70860927e-01]
     [ 4.44089210e-16  2.22044605e-16]]
    1.7256332301709633e-31



```python
x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
plt.contour(X, Y, f2(X, Y), 10)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
pass
```


```python

```

6f5b4fef3458b1d60813d288f43eb27675cc5cee
ipynb
Jupyter Notebook
homework/Homework08_Boyao.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
homework/Homework08_Boyao.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
homework/Homework08_Boyao.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
Qwen/Qwen-72B
1. YES
2. YES