# Regression

Regression is the process of approximating noisy data with functions.  The same numerical methods can also be used to approximate a complicated function with a simpler function.  We'll begin by looking at some limited cases in terms of linear algebra.

## Polynomial regression

Long ago (in the Linear Algebra notebook), we solved an over-determined linear system to compute a polynomial approximation of a function.
Sometimes we make approximations because we are interested in the polynomial coefficients.
That is usually only for very low order (like linear or quadratic fits).
Inferring higher coefficients is ill-conditioned (as we saw with Vandermonde matrices) and probably not meaningful.
We now know that Chebyshev bases are good for representing high-degree polynomials, but if the points are arbitrarily spaced, how many do we need for the Chebyshev approximation to be well-conditioned?


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

def cosspace(a, b, n=50):
    return (a + b)/2 + (b - a)/2 * (
        np.cos(np.linspace(-np.pi, 0, n)))

def vander_chebyshev(x, n=None):
    if n is None:
        n = len(x)
    T = np.ones((len(x), n))
    if n > 1:
        T[:,1] = x
    for k in range(1,n-1):
        T[:,k+1] = 2 * x * T[:,k] - T[:,k-1]
    return T

def runge1(x):
    return 1 / (1 + 10*x**2)
```


```python
def chebyshev_regress_eval(x, xx, n):
    V = vander_chebyshev(x, n)
    Q, R = np.linalg.qr(V)
    return vander_chebyshev(xx, n) @ np.linalg.solve(R, Q.T)

xx = np.linspace(-1, 1, 100)
x = np.linspace(-1, 1, 80)
plt.plot(x, runge1(x), '.k')
plt.plot(xx, chebyshev_regress_eval(x, xx, 40) @ runge1(x), label='p(x)')
plt.plot(xx, runge1(xx), '-.', label='runge1(x)')
plt.legend(loc='upper left');
```

* What is the degree $k$ of the polynomial that is used?
* What distribution of points is used for the data?
* Would this have artifacts if we used a degree $k$ polynomial to interpolate noise-free samples at $k+1$ points using the same distribution?


```python
ns = np.geomspace(5, 1000, dtype=int)
conds = [np.linalg.cond(vander_chebyshev(np.linspace(-1, 1, n),
                                         int(n**(.5))))
         for n in ns]
plt.semilogy(ns, conds)
plt.xlabel('$n$')
plt.ylabel('cond');
```

* If we have $n$ data points, can we use a polynomial of degree $k = \lfloor n/5 \rfloor$?
* What about $k = n^{3/4}$?
* What expression $k = ?(n)$ appears to be sufficient?

## Noisy data

Regression really comes into its own when the data is noisy.


```python
def runge1_noisy(x, sigma):
    return runge1(x) + np.random.randn(*x.shape)*sigma

x = np.linspace(-1, 1, 500)
y = runge1_noisy(x, 0.2)
plt.plot(x, y, '.')
plt.plot(x, runge1(x), 'k', label='runge1(x)')
plt.plot(x, chebyshev_regress_eval(x, x, 7) @ y, label='regress(x)')
plt.legend(loc='upper left');
```

### Probability distributions

In order to interpret real data, we need a model for the noise.  We have chosen the most common and computationally convenient choice when creating the synthetic data above.
The function `randn` draws samples from the "standard normal" or Gaussian distribution.

$$ p(t) = \frac{1}{\sqrt{2\pi}} e^{-t^2/2} . $$


```python
def stdnormal(t):
    return np.exp(-t**2/2) / np.sqrt(2*np.pi)

n = 1000
w = np.random.randn(n)
plt.hist(w, bins=40, range=(-3,3), density=True)
t = np.linspace(-3, 3)
plt.plot(t, stdnormal(t))
plt.xlabel('$t$')
plt.ylabel('$p(t)$');
w.mean()
```

### Regression on noisy data

We can just go run our regression algorithm on the noisy data.


```python
x = np.linspace(-1, 1, 500)
y = runge1_noisy(x, sigma=0.1)
```


```python
plt.plot(x, y, '.')
plt.plot(xx, chebyshev_regress_eval(x, xx, 9) @ y, label='p(x)')
plt.plot(xx, runge1(xx), label='runge1(x)')
plt.legend(loc='upper left');
```

## What problem are we solving?

### Why do we call it a linear model?

We are currently working with algorithms that express the regression as a linear function of the model parameters.  That is, we search for coefficients $c = [c_0, c_1, \dotsc]^T$ such that

$$ V(x) c \approx y $$

where the left hand side is linear in $c$.  In different notation, we are searching for a predictive model

$$ f(x_i, c) \approx y_i \text{ for all $(x_i, y_i)$} $$

that is linear in $c$.

### Assumptions

1. The independent variables $x$ are error-free
1. The prediction (or "response") $f(x,c)$ is linear in $c$
1. The noise in the measurements $y$ is independent (uncorrelated)
1. The noise in the measurements $y$ has constant variance

There are reasons why all of these assumptions may be undesirable in practice, thus leading to more complicated methods.

### Loss functions

The error in a single prediction $f(x_i,c)$ of an observation $(x_i, y_i)$ is often measured as
$$ \frac 1 2 \big( f(x_i, c) - y_i \big)^2, $$
which turns out to have a statistical interpretation when the noise is normally distributed.
It is natural to define the error over the entire data set as
\begin{align} L(c; x, y) &= \sum_i \frac 1 2 \big( f(x_i, c) - y_i \big)^2 \\
&= \frac 1 2 \lVert f(x, c) - y \rVert^2
\end{align}
where I've used the notation $f(x,c)$ to mean the vector resulting from gathering all of the outputs $f(x_i, c)$.
The function $L$ is called the "loss function" and is the key to relaxing the above assumptions.

### Partial derivatives

Let's step back from optimization and consider how to differentiate a function of several variables.  Let $f(\boldsymbol x)$ be a function of a vector $\boldsymbol x$.  For example,

$$ f(\boldsymbol x) = x_0^2 + \sin(x_1) e^{3x_2} . $$

We can evaluate the **partial derivative** by differentiating with respect to each component $x_i$ separately (holding the others constant), and collect the result in a vector,

\begin{align}
\frac{\partial f}{\partial \boldsymbol x} &= \begin{bmatrix} \frac{\partial f}{\partial x_0} & \frac{\partial f}{\partial x_1} & \frac{\partial f}{\partial x_2} \end{bmatrix} \\
&= \begin{bmatrix} 2 x_0 & \cos(x_1) e^{3 x_2} & 3 \sin(x_1) e^{3 x_2} \end{bmatrix}.
\end{align}

Now let's consider a vector-valued function $\boldsymbol f(\boldsymbol x)$, e.g.,

$$ \boldsymbol f(\boldsymbol x) = \begin{bmatrix} x_0^2 + \sin(x_1) e^{3x_2} \\ x_0 x_1^2 / x_3 \end{bmatrix} . $$

and write the derivative as a matrix,

\begin{align}
\frac{\partial \boldsymbol f}{\partial \boldsymbol x} &=
\begin{bmatrix} \frac{\partial f_0}{\partial x_0} & \frac{\partial f_0}{\partial x_1} & \frac{\partial f_0}{\partial x_2} \\
\frac{\partial f_1}{\partial x_0} & \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} \\
\end{bmatrix} \\
&= \begin{bmatrix} 2 x_0 & \cos(x_1) e^{3 x_2} & 3 \sin(x_1) e^{3 x_2} \\
x_1^2 / x_3 & 2 x_0 x_1 / x_3 & -x_0 x_1^2 / x_3^2
\end{bmatrix}.
\end{align}

#### Geometry of derivatives



* Handy resource on partial derivatives for matrices and vectors: https://explained.ai/matrix-calculus/index.html#sec3

#### Derivative of dot product

Let $f(\boldsymbol x) = \boldsymbol y^T \boldsymbol x = \sum_i y_i x_i$ and compute the derivative

$$ \frac{\partial f}{\partial \boldsymbol x} = \begin{bmatrix} y_0 & y_1 & \dotsb \end{bmatrix} = \boldsymbol y^T . $$

Note that $\boldsymbol y^T \boldsymbol x = \boldsymbol x^T \boldsymbol y$ and we have the product rule,

$$ \frac{\partial \lVert \boldsymbol x \rVert^2}{\partial \boldsymbol x} = \frac{\partial \boldsymbol x^T \boldsymbol x}{\partial \boldsymbol x} = 2 \boldsymbol x^T . $$

Also,
$$ \frac{\partial \lVert \boldsymbol x - \boldsymbol y \rVert^2}{\partial \boldsymbol x} = \frac{\partial (\boldsymbol x - \boldsymbol y)^T (\boldsymbol x - \boldsymbol y)}{\partial \boldsymbol x} = 2 (\boldsymbol x - \boldsymbol y)^T .$$

### Optimization
Given data $(x,y)$ and loss function $L(c; x,y)$, we wish to find the coefficients $c$ that minimize the loss, thus yielding the "best predictor" (in a sense that can be made statistically precise).  I.e.,
$$ \bar c = \arg\min_c L(c; x,y) . $$

It is usually desirable to design models such that the loss function is differentiable with respect to the coefficients $c$, because this allows the use of more efficient optimization methods.  Recall that our forward model is given in terms of the Vandermonde matrix,

$$ f(x, c) = V(x) c $$

and thus

$$ \frac{\partial f}{\partial c} = V(x) . $$

We can now differentiate our loss function
$$ L(c; x, y) = \frac 1 2 \lVert f(x, c) - y \rVert^2 = \frac 1 2 \sum_i (f(x_i,c) - y_i)^2 $$
term-by-term as
\begin{align} \nabla_c L(c; x,y) = \frac{\partial L(c; x,y)}{\partial c} &= \sum_i \big( f(x_i, c) - y_i \big) \frac{\partial f(x_i, c)}{\partial c} \\
&= \sum_i \big( f(x_i, c) - y_i \big) V(x_i)
\end{align}
where $V(x_i)$ is the $i$th row of $V(x)$.
Alternatively, we can take a more linear algebraic approach to write the same expression is
\begin{align} \nabla_c L(c; x,y) &= \big( f(x,c) - y \big)^T V(x) \\
&= \big(V(x) c - y \big)^T V(x) \\
&= V(x)^T \big( V(x) c - y \big) .
\end{align}
A necessary condition for the loss function to be minimized is that $\nabla_c L(c; x,y) = 0$.

* Is the condition sufficient for general $f(x, c)$?
* Is the condition sufficient for the linear model $f(x,c) = V(x) c$?
* Have we seen this sort of equation before?

##### Uniqueness
We can read the expression $ V(x)^T \big( V(x) c - y \big) = 0$ as saying that the residual $V(x) c - y$ is orthogonal to the range of $V(x)$.  Suppose $c$ satisfies this equation and $c' \ne c$ is some other value of the coefficients.  Then
\begin{align}
V(x)^T \big( V(x) c' - y \big) &= V(x)^T V(x) c' - V(x)^T y \\
&= V(x)^T V(x) c' - V(x)^T V(x) c \\
&= V(x)^T V(x) (c' - c) \ne 0
\end{align}
whenever $V(x)^T V(x)$ is nonsingular, which happens any time $V(x)$ has full column rank.

##### Gradient descent
Instead of solving the least squares problem using linear algebra (QR factorization), we could solve it using gradient descent.  That is, on each iteration, we'll take a step in the direction of the negative gradient.


```python
def grad_descent(loss, grad, c0, gamma=1e-3, tol=1e-5):
    """Minimize loss(c) via gradient descent with initial guess c0
    using learning rate gamma.  Declares convergence when gradient
    is less than tol or after 500 steps.
    """
    c = c0.copy()
    chist = [c.copy()]
    lhist = [loss(c)]
    for it in range(500):
        g = grad(c)
        c -= gamma * g
        chist.append(c.copy())
        lhist.append(loss(c))
        if np.linalg.norm(g) < tol:
            break
    return c, np.array(chist), np.array(lhist)

class quadratic_loss:
    """Test problem to give example of gradient descent."""
    def __init__(self):
        self.A = np.array([[1, 1], [1, 4]])
    def loss(self, c):
        return .5 * c @ self.A @ c
    def grad(self, c):
        return self.A @ c

def test(gamma):
    q = quadratic_loss()
    c, chist, lhist = grad_descent(q.loss, q.grad, .9*np.ones(2), gamma=gamma)
    plt.semilogy(lhist)
    plt.ylabel('loss')
    plt.xlabel('cumulative iterations')

    plt.figure()
    l = np.linspace(-1, 1)
    x, y = np.meshgrid(l, l)
    z = [q.loss(np.array([x[i,j], y[i,j]])) for i in range(50) for j in range(50)]
    plt.contour(x, y, np.reshape(z, x.shape))
    plt.plot(chist[:,0], chist[:,1], 'o-')
    plt.title('gamma={}'.format(gamma))
    
test(.45)
```


```python
class chebyshev_regress:
    def __init__(self, x, y, n):
        self.V = vander_chebyshev(x, n)
        self.y = y
        self.n = n
        
    def init(self):
        return np.zeros(self.n)

    def loss(self, c):
        r = self.V @ c - y
        return 0.5 * (r @ r)
    
    def grad(self, c):
        r = self.V @ c - y
        return self.V.T @ r
    
reg = chebyshev_regress(x, y, 6)
c, _, lhist = grad_descent(reg.loss, reg.grad, reg.init(), gamma=2e-3)
plt.semilogx(np.arange(1, 1+len(lhist)), lhist)
plt.ylabel('loss')
plt.xlabel('cumulative iterations')
plt.ylim(bottom=0);
```


```python
np.linalg.cond(reg.V.T @ reg.V)
```




    8.108087894474187



* How does changing the "learning rate" `gamma` affect convergence?
* How does changing the number of basis functions affect convergence rate?

#### Did we find the same solution?
We intend to solve the same problem as we previously solved using QR, so we hope to find the same solution.


```python
plt.plot(xx, chebyshev_regress_eval(x, xx, 6) @ y, '-k', label='QR')
plt.plot(xx, vander_chebyshev(xx, 6) @ c, '--', label='grad_descent')
plt.plot(x, y, '.')
plt.legend();
```

#### Observations

* This algorithm is kind of finnicky:
 * It takes many iterations to converge
 * The "learning rate" needs to be empirically tuned
 * When not tuned well, the algorithm can diverge
* It could be made more robust using a line search
* The $QR$ algorithm is a more efficient and robust way to solve these problems

## Nonlinear regression

Instead of the linear model
$$ f(x,c) = V(x) c = c_0 + c_1 \underbrace{x}_{T_1(x)} + c_2 T_2(x) + \dotsb $$
let's consider a rational model with only three parameters
$$ f(x,c) = \frac{1}{c_0 + c_1 x + c_2 x^2} = (c_0 + c_1 x + c_2 x^2)^{-1} . $$
We'll use the same loss function
$$ L(c; x,y) = \frac 1 2 \lVert f(x,c) - y \rVert^2 . $$

We will also need the gradient
$$ \nabla_c L(c; x,y) = \big( f(x,c) - y \big)^T \nabla_c f(x,c) $$
where
\begin{align}
\frac{\partial f(x,c)}{\partial c_0} &= -(c_0 + c_1 x + c_2 x^2)^{-2} = - f(x,c)^2 \\
\frac{\partial f(x,c)}{\partial c_1} &= -(c_0 + c_1 x + c_2 x^2)^{-2} x = - f(x,c)^2 x \\
\frac{\partial f(x,c)}{\partial c_2} &= -(c_0 + c_1 x + c_2 x^2)^{-2} x^2 = - f(x,c)^2 x^2 .
\end{align}


```python
class rational_regress:
    def __init__(self, x, y):
        self.x = x
        self.y = y
        self.n = 3
        
    def init(self):
        return np.ones(self.n)
    
    def f(self, c):
        x = self.x
        return 1 / (c[0] + c[1]*x + c[2]*x**2)
        
    def df(self, c):
        x = self.x
        f2 = self.f(c)**2
        return np.array([-f2, -f2*x, -f2*x**2]).T

    def loss(self, c):
        r = self.f(c) - self.y
        return 0.5 * (r @ r)
    
    def grad(self, c):
        r = self.f(c) - self.y
        return r @ self.df(c)

reg = rational_regress(x, y)
c, _, lhist = grad_descent(reg.loss, reg.grad, reg.init(), gamma=2e-2)
plt.semilogx(np.arange(1, 1+len(lhist)), lhist)
plt.ylabel('loss')
plt.xlabel('cumulative iterations')
plt.ylim(bottom=0)
plt.title('rational c={}'.format(c));
```


```python
x = np.linspace(-1, 1, 500)
y = runge1_noisy(x, sigma=0.1)

c0 = np.array([1, 0, 1.])
c, _, lhist = grad_descent(reg.loss, reg.grad, c0, gamma=2e-2)
plt.semilogx(np.arange(1, 1+len(lhist)), lhist)
plt.ylabel('loss')
plt.xlabel('cumulative iterations')
plt.ylim(bottom=0)
plt.title('rational c={}'.format(c));
```


```python
plt.plot(x, y, '.')
plt.plot(xx, runge1(xx), label='runge1(x)')
plt.plot(xx, chebyshev_regress_eval(x, xx, 40) @ y, label='Polynomial')
plt.plot(xx, 1/(c[0] + c[1]*xx + c[2]*xx**2), '--k', label='Rational')
plt.legend();
```

#### Observations

* There can be local minima or points that look an awful lot like local minima.
* Convergence is sensitive to learning rate.
* It takes a lot of iterations to converge.
* A well-parametrized model (such as the rational model above) can accurately reconstruct the mean even with very noisy data.

### Outlook

* The optimization problem can be solved using a Newton method.  It can be onerous to implement the needed derivatives.
* The [Gauss-Newton method](https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm) (see activity) is often more practical than Newton while being faster than gradient descent, though it lacks robustness.
* The [Levenberg-Marquardt method](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) provides a sort of middle-ground between Gauss-Newton and gradient descent.
* Many globalization techniques are used for models that possess many local minima.
* One pervasive approach is stochastic gradient descent, where small batches (e.g., 1 or 10 or 20) are selected randomly from the corpus of observations (500 in our current example), and a step of gradient descent is applied to that reduced set of observations.  This helps to escape saddle points and weak local minima.
* Among expressive models $f(x,c)$, some may converge much more easily than others.  Having a good optimization algorithm is essential for nonlinear regression with complicated models, especially those with many parameters $c$.
* Classification is a very similar problem to regression, but the observations $y$ are discrete, thus

 * models $f(x,c)$ must have discrete output
 * the least squares loss function is not appropriate.
* [Why momentum really works](https://distill.pub/2017/momentum/)

95beca73965f65f472957065640a055460c8a7a8
ipynb
Jupyter Notebook
Regression.ipynb
cu-numcomp/numcomp-class
8d9f5944774082ad4d8e8462e3afae7a87cddf4a
BSD-2-Clause
2020-01-14T17:00:27.000Z
2020-07-03T14:17:07.000Z
Regression.ipynb
cu-numcomp/numcomp-class
8d9f5944774082ad4d8e8462e3afae7a87cddf4a
BSD-2-Clause
2020-02-28T08:39:35.000Z
2020-03-17T20:42:56.000Z
Regression.ipynb
cu-numcomp/numcomp-class
8d9f5944774082ad4d8e8462e3afae7a87cddf4a
BSD-2-Clause
2020-01-16T19:52:55.000Z
2021-02-09T18:30:56.000Z
Qwen/Qwen-72B
1. YES
2. YES