```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, fixed
from matplotlib import cm
plt.rcParams['font.size'] = 14
plt.rcParams['axes.spines.right'] = False
plt.rcParams['ytick.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['xtick.top'] = False

# Get n grayscale colors
def getGrayColors(n):
    colors = [];
    for i in range(1, n+1):
        colors.append(i/(n+1.)*np.ones(3))
    return colors
```

### Example case 

Suppose we have observed spike responses $y_i$ to some one-dimensional stimulus $x_i$, $i \in 1,\dots, N$. Furthermore, suppose that the cell exhibited some thresholding mechanism: hardly any spikes were observed for stimulus values below a certain threshold, whereas the response grew roughly linearly above the threshold. We would like to describe the cell's stimulus response mapping using a simple model, and one alternative is a soft-rectifying function:
\begin{equation}
 \hat{y_i} = f(z_i) = 
 \begin{cases}
  \exp(z_i),& \text{if  } z_i \leq 0\\
  z_i + 1,& \text{if  } z_i > 0\\
 \end{cases}
\end{equation}
where $\hat{y_i}$ is the predicted output, $z_i=w_0+w_1x_i$, $w_0$ the threshold, and $w_1$ the slope. However, it is often simpler to work with vector notations, in such cases, we define $z_i$ as a dot product between a parameter vector and a stimulus vector $\mathbf{w}^T\mathbf{x}_i$, where $\mathbf{w}=[w_0,\; w_1]$ and $\mathbf{x}_i=[1; x_i]$.

Below, we generate simulated data from the proposed model, and we illustrate what the proposed model looks like for various choices of $w_0$ and $w_1$.


```python
# The soft-rectifying function (f), also known as the inverse of the link function for generalized linear models
def invLinkFun(X, w):
    z = np.dot(X, w)
    mu = np.zeros(z.shape)
    mu[z<0] = np.exp(z[z<0])
    mu[z>=0] = z[z>=0] + 1
    return mu

# Parameters
n = 200;           # Number of data points
wTrue = [-1, 0.5]  # w0 and w1 for the simulated data

# Generate example data
x = np.random.rand(n)*30 - 15
x = np.sort(x)
X = np.vstack([-np.ones(n), x]).T
mu = invLinkFun(X, wTrue)
y = np.random.poisson(mu)

# Callback plotting function for interactive plots
def plot_fun(x, y, w0, w1):
    fig = plt.figure(figsize=(7.5, 5))
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(x, y, 'k.', label='data')
    wTmp = np.array([w0, w1])
    muTmp = invLinkFun(X, wTmp)
    ax.plot(x, muTmp, '-', color=0.5*np.ones(3))
    ax.set_ylabel('Spike count')
    ax.set_xlabel('x')
    ax.set_xticks([-10, 0, 10]);
    ax.set_ylim([-1, 11]);

interact(plot_fun, x=fixed(x), y=fixed(y), w0=(-5., 5.), w1=(-5., 5.));
```


    interactive(children=(FloatSlider(value=0.0, description='w0', max=5.0, min=-5.0), FloatSlider(value=0.0, desc…


### Objective functions

A manual search for good parameter values can be quite tricky and cumbersome, especially if there are more than just two model parameters. An automatic way would be a lot simpler, but for that we need some measure of how "good" a specific selection of model parameters are. One natural way of quantifying the "goodness" of parameters is through an objective function which scores each parameter combination. It might seen tricky to come up with such a function, but one common approach is to define it to be a likelihood function. The idea being that the likelihood function tells you how likely you are to observe the current data given some model parameters. For example, if we assume that the observed spike counts come from a Poisson distribution whose mean is given by our soft-rectifying function, then the probability $p$ of observing a specific data point ($x_i, y_i$) is given by:
\begin{equation}
 p(y_i|x_i, w_0, w_1) = \frac{f(z_i)^{y_i} \exp(-f(z_i))}{y_i!}.
\end{equation}
The likelihood function $l$ is then defined as the probability of observing all data points, where independence is essentially always assumed to simplify the expression into a product: 
\begin{equation}
 l(w_0, w_1) = \prod_i^N \frac{f(z_i)^{y_i} \exp(-f(z_i))}{y_i!}.
\end{equation}
The likelihood function is, however, often cumbersome to work with, and people tend to use the log of the likelihood function instead. Logarithms turns products into sums, and the log-likelihood function $ll$ thus has the form:
\begin{equation}
 ll(w_0, w_1) = \sum_i^N y_i\log(f(z_i)) -f(z_i) -\log(y_i!).
\end{equation}
Finally, we note that the last term can be ignored (it does not depend on our model parameters), whereupon we are left with:
\begin{equation}
 ll(w_0, w_1) = \sum_i^N y_i\log(f(z_i)) -f(z_i),
\end{equation}
which also corresponds to the objective function for **Poisson regression**.

Equipped with a function that quantifies the "goodness" of parameter choices, we can now select parameter values that maximizes the log-likelihood function.


```python
# Poisson regression, log-likelihood
logLikFun = lambda X, w, y: np.sum(y*np.log(invLinkFun(X, w))-invLinkFun(X, w))
```

##### Exampel 1, find w0 given w1


```python
# Example parameters
w0Ex = np.array([wTrue[0]-4, wTrue[0], wTrue[0]+4])
w0Vals = np.linspace(wTrue[0]-5, wTrue[0]+5, 101)

# Log-likelihood function over the whole interval
ll = np.empty_like(w0Vals)
for i, w0 in enumerate(w0Vals):
    wTmp = np.array([w0, wTrue[1]])
    ll[i] = logLikFun(X, wTmp, y)

# Log-likelihood values for our example parameters
llEx = np.empty_like(w0Ex)
for i, w0 in enumerate(w0Ex):
    wTmp = np.array([w0, wTrue[1]])
    llEx[i] = logLikFun(X, wTmp, y)

# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(15, 5))
# Log-likelihood
ax = fig.add_subplot(1, 2, 1)
ax.plot(w0Vals, ll, 'k-')
for i in range(w0Ex.size):
    ax.plot(w0Ex[i], llEx[i], 'o', ms=10, color=grayColors[i])
ax.set_xlabel('$w_0$')
ax.set_ylabel('Log-likelihood')
# Data and example response function mappings
ax = fig.add_subplot(1, 2, 2)
ax.plot(x, y, 'k.', label='data')
for i in range(w0Ex.size):
    wTmp = np.array([w0Ex[i], wTrue[1]])
    pTmp = invLinkFun(X, wTmp)
    ax.plot(x, pTmp, '-', color=grayColors[i])
ax.set_ylabel('Spike count')
ax.set_xlabel('x')
ax.set_xticks([-10, 0, 10]);
```

##### Example 2, find both w0 and w1


```python
# Example parameters
wEx = np.array([[wTrue[0]-0.9, wTrue[0], wTrue[0]+0.9], [wTrue[1]+0.3, wTrue[1], wTrue[1]-0.3]])

# Get w0 and w1 combinations over a grid
nGrid = 51
W0, W1 = np.meshgrid(np.linspace(wTrue[0]-1, wTrue[0]+1, nGrid), np.linspace(wTrue[1]-0.5, wTrue[1]+0.5, nGrid))

# Get the log-likelihood for each parameter combination
llVals = np.zeros([nGrid, nGrid])
for i in range(nGrid):
    for j in range(nGrid):
        wTmp = np.array([W0[i, j], W1[i, j]])
        llVals[i, j] = logLikFun(X, wTmp, y)

# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(15, 5))
# Log-lieklihood and example parameters
ax = plt.subplot(1, 2, 1)
contourHandle = ax.contourf(W0, W1, llVals, 50, cmap=cm.coolwarm)
for i in range(wEx.shape[1]):
    ax.plot(wEx[0, i], wEx[1, i], 'o', ms=10, color=grayColors[i])
ax.set_xlabel('$w_0$')
ax.set_ylabel('$w_1$');
cBarHandle = plt.colorbar(contourHandle)
cBarHandle.set_label('Log-likelihood')
# Data and example response function mappings
ax = fig.add_subplot(1, 2, 2)
ax.plot(x, y, 'k.', label='data')
for i in range(wEx.shape[1]):
    wTmp = wEx[:, i]
    pTmp = invLinkFun(X, wTmp)
    ax.plot(x, pTmp, '-', color=grayColors[i])
ax.set_ylabel('Spike count')
ax.set_xlabel('x')
ax.set_xticks([-10, 0, 10]);
```

### Optimization

At first sight, it might look like the objective function only slightly simplifies our problem of finding good parameter values: finding the maxima is easy when you see the objective function, but visualizing the function for more than three parameters is extremely difficult. Luckily, however, there is a whole optimization field devoted to finding the maxima or minima of functions using methods that don't require you to visualize the objective function. However, before continuing, we will modify our objective functions once more by multiplying it with -1. This is just to turn the maximization problem into a minimization problem, which is usually the default notation. Finding the best parameters values thus corresponds to minimizing the negative log-likelihood $nll$: 

\begin{equation}
 nll(w_0, w_1) = \sum_i^N f(z_i) - y_i\log(f(z_i)).
\end{equation}


```python
# Poisson regression, negative log-likelihood
negLogLikFun = lambda X, w, y: np.sum(invLinkFun(X, w) - y*np.log(invLinkFun(X, w)))
```

### Gradient

Several commonly used optimization methods utilize gradients. The gradient of a function is by definition a vector consisting of the partial derivatives with respect to each parameter. For our example from above, the gradient thus becomes:
\begin{equation}
 \nabla nll = \frac{\partial nll}{\partial \mathbf{w}} = \left[ \frac{\partial nll}{\partial w_0}, \frac{\partial nll}{\partial w_1} \right], 
\end{equation}
where $\nabla$ denotes the linear nabla operator. The fact that $\nabla$ is linear means that we can obtain the gradient of the negative log-likelihood function by just summing up the gradients of each term in the sum:
\begin{align}
 \nabla nll &= \sum_i^N \frac{\partial f(z_i) - y_i\log(f(z_i))}{\partial \mathbf{w}}, \\
 \frac{\partial f(z_i) - y_i\log(f(z_i))}{\partial \mathbf{w}} &= 
 \begin{cases}
  (\exp(z_i) - y_i)\mathbf{x}_i,& \text{if  } z_i \leq 0\\
  \left( 1-\frac{y_i}{z_i+1} \right) \mathbf{x}_i,& \text{if  } z_i > 0\\
 \end{cases}.
\end{align}
We can thus compute the gradient for any parameter values of our choice, but the real benefit of this possibility only becomes apparent when considering what the gradient really denotes. The partial derivative of any parameter tells us how quickly the function increases when that parameter is changed slightly, and consequently the gradient, which is the vector containing all partial derivatives, indicates in which direction in parameter space that the function value increases the most.


```python
# Poisson regression, gradient of the negative log-likelihood
def negLogLikDerFun(X, w, y):
    z = np.dot(X, w)
    der = np.dot(np.exp(z[z<0])-y[z<0], X[z<0, :])
    der += np.dot(1-y[z>=0]/(z[z>=0]+1), X[z>=0, :])
    return der
```

##### Example 3, the gradient of the negative log-likelihood function


```python
# Use a lower resolution for gradient evaluations
nGrid = 21
W0q, W1q = np.meshgrid(np.linspace(wTrue[0]-1, wTrue[0]+1, nGrid), np.linspace(wTrue[1]-0.5, wTrue[1]+0.5, nGrid))

# Get the gradient for each combination
nllVals = np.zeros([nGrid, nGrid])
W0der = np.zeros([nGrid, nGrid])
W1der = np.zeros([nGrid, nGrid])
for i in range(nGrid):
    for j in range(nGrid):
        wTmp = np.array([W0q[i, j], W1q[i, j]])
        derTmp =  negLogLikDerFun(X, wTmp, y)
        W0der[i, j] = derTmp[0]
        W1der[i, j] = derTmp[1]

# Normalize to unit length
W0derNorm = W0der / np.sqrt(W0der**2. + W1der**2.)
W1derNorm = W1der / np.sqrt(W0der**2. + W1der**2.)
        
# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(25, 7))
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
# Negative log-likelihood
contourHandle = ax.contourf(W0, W1, -llVals, 50, cmap=cm.coolwarm)
# Gradient
ax.quiver(W0q, W1q, W0derNorm, W1derNorm, angles='xy', scale_units='xy', scale=30, width=2e-3)
# Axes settings
ax.set_xlabel('$w_0$')
ax.set_ylabel('$w_1$');
cBarHandle = plt.colorbar(contourHandle)
cBarHandle.set_label('negative log-likelihood')
```

### Gradient descent

A stupidly simple idea for finding the minima of the negative log-likelihood function would then be to start from some parameters $\mathbf{w}_\mathrm{init}$ and then just iteratively move in the opposite direction to the gradient. That is, we evaluate the gradient for our current parameter guess and then we update our guess by subtracting the gradient times a constant $\eta$. More specifically, at each step $k$ of our iterative algorithm, we update our parameters as:
\begin{equation}
 \mathbf{w}_{k+1} = \mathbf{w}_k - \eta \nabla nll(\mathbf{w}_k).
\end{equation}
For sufficiently small $\eta$, this will guarantee that we backtrack along the gradients all the way to the minima of the negative log-likelihood function, as illustrated below.


```python
# Gradient descent parameters
eta = 2e-4
nIterations = 25
wInit = np.array([-1.8, 0.85])
wUpdates = np.zeros([nIterations, 2])
wUpdates[0, :] = wInit
nllVals = np.zeros(nIterations)
nllVals[0] = negLogLikFun(X, wUpdates[0, :], y)
   
# Gradient descent loop
for i in range(1, nIterations):
    wUpdates[i, :] = wUpdates[i-1, :] - eta*negLogLikDerFun(X, wUpdates[i-1, :], y)
    nllVals[i] = negLogLikFun(X, wUpdates[i, :], y)
    
# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(25, 7))
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
# Negative log-likelihood
contourHandle = ax.contourf(W0, W1, -llVals, 50, cmap=cm.coolwarm)
# Gradient
ax.quiver(W0q, W1q, W0derNorm, W1derNorm, angles='xy', scale_units='xy', scale=30, width=2e-3)
# Gradient descent progress
ax.plot(wUpdates[:, 0], wUpdates[:, 1], '-o', ms=6, color='white', lw=3)
# Axes settings
ax.set_xlabel('$w_0$')
ax.set_ylabel('$w_1$');
cBarHandle = plt.colorbar(contourHandle)
cBarHandle.set_label('negative log-likelihood')
```

### Backtracking line search

Finding a suitable $\eta$ value for the gradient descent approach might not be trivial: a too small value of $\eta$ will hamper progress, whereas a too large value will cause the algorithm to diverge. It would therefore be very convenient to have an adaptive method that automatically adjusted $\eta$ to take large but not too large steps on each iteration of gradient descent. One common method for convex problems is to use something called backtracking line search. The idea being that we start of with a large $\eta$ value and then successively decrease it by a factor $\beta$. This continues until the observed decrease of the objective function is within a factor $\alpha$ of what would be expected for a linear approximation of the objective function. This is illustrated below for a one-dimensional case (find $w_0$ given $w_1$) by showing the real negative log-likelihood and the linear approximation around our current guess of $w_0$. 


```python
w0Init = w0Vals.min()
nllInit = negLogLikFun(X, np.array([w0Init, wTrue[1]]), y)
w0InitDer = negLogLikDerFun(X, np.array([w0Init, wTrue[1]]), y)[0]
offset = nllInit - w0Init*w0InitDer;

# Calculate real and expected decreases
realNllDecrease =  -ll - nllInit
expectedNllDecrease = offset  + w0InitDer*w0Vals - nllInit

# Plotting
fig = plt.figure(figsize=(15, 5))
# Log-likelihood
ax = fig.add_subplot(1, 2, 1)
ax.plot(w0Vals, -ll, 'k-', label='Real nll')
ax.plot(w0Vals, offset  + w0InitDer*w0Vals, 'k:', label='linear approx.')
ax.plot(w0Init, nllInit, 'o', ms=10, color='k')
ax.set_xlabel('$w_0$')
ax.set_ylabel('Negative log-likelihood')
ax.legend()
# Log-likelihood decrease
ax = fig.add_subplot(1, 2, 2)
ax.plot(w0Vals[1:], realNllDecrease[1:] / expectedNllDecrease[1:], 'k-')
ax.set_xlabel('$w_0$')
ax.set_ylabel('nll / linear approx.');
```

The linear approximation will always lie below the real negative log-likelihood for convex objective functions, but it provides a good approximation close to $w_0$. The idea of backtracking line search is thus to take the largest step possible while still ensuring that the linear approximation is good enough. Goodness is measured by how close the real objective function is to the linear approximation, and the acceptance criteria is set by $\alpha$. $\alpha=0.5$ thus means that we want the largest step size for which the decrease observed for the real objective function is at least half of what would be predicted by the linear approximation. However, we don't want to spend time solving this line search problem exactly, rather, we start of with a high value of $\eta$ and decrease it successively by a factor $\beta$ until the goodness criteria is met. That is, we decrease $\eta$ to $\eta\beta$ until
\begin{equation}
 nll(\mathbf{w}_{k}-\eta \nabla nll(\mathbf{w}_k)) > nll(\mathbf{w}_k) - \alpha \eta \nabla nll(\mathbf{w}_k)^T\nabla nll(\mathbf{w}_k),
\end{equation}
where $\eta \nabla nll(\mathbf{w}_k)^T\nabla nll(\mathbf{w}_k)$ is the predicted decrease for our linear approximation. Standard values for $\alpha$ and $\beta$ are 0.5 and 0.8, respectively. The illustration below use these standard values to show how gradient descent progresses for our example problem of finding both $w_0$ and $w_1$ using a backtracking line search implementation.


```python
# Gradient descent parameters
eta = 0.1
alpha = 0.5
beta = 0.8
wUpdates = np.zeros([nIterations, 2])
wUpdates[0, :] = wInit
nllVals = np.zeros(nIterations)
nllVals[0] = negLogLikFun(X, wUpdates[0, :], y)

# Gradient descent with backtracking
for i in range(1, nIterations):
    wUpdates[i, :] = wUpdates[i-1, :]
    delta = negLogLikDerFun(X, wUpdates[i, :], y)
    tooLarge = True
    # Backtracking loop
    while tooLarge:
        wUpdates[i, :] -= eta*delta
        nllTmp = negLogLikFun(X, wUpdates[i, :], y)
        if nllTmp > nllVals[i-1] - alpha*eta*np.dot(delta, delta):
            wUpdates[i, :] += eta*delta
            eta *= beta
        else:
            tooLarge = False
            nllVals[i] = nllTmp
    
# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(25, 7))
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
# Negative log-likelihood
contourHandle = ax.contourf(W0, W1, -llVals, 50, cmap=cm.coolwarm)
# Gradient
ax.quiver(W0q, W1q, W0derNorm, W1derNorm, angles='xy', scale_units='xy', scale=30, width=2e-3)
# Gradient descent progress
ax.plot(wUpdates[:, 0], wUpdates[:, 1], '-o', ms=6, color='white', lw=3)
# Axes settings
ax.set_xlabel('$w_0$')
ax.set_ylabel('$w_1$');
cBarHandle = plt.colorbar(contourHandle)
cBarHandle.set_label('negative log-likelihood')
```

### Acceleration

Gradient descent can be made to converge faster using a method called Nesterov accelerated gradient, named after Yurii Nesterov that introduced it. The basic idea behind Nesterov's accelerated gradient method is to first make a move in the direction of the previous parameter change, and then evaluate the gradient at this intermediate position $v$ in order to let the standard gradient descent rule make a correction step. That is, at each iteration, we first determine intermediate parameter values $v$, and then take a step in the direction of the negative gradient evaluated at $v$:
\begin{align}
 \mathbf{v}_{k+1} &= \mathbf{w}_k + \frac{k-2}{k+1}(\mathbf{w}_k - \mathbf{w}_{k-1}), \\
 \mathbf{w}_{k+1} &= \mathbf{v}_{k+1} - \eta \nabla nll(\mathbf{v}_{k+1}). 
\end{align}
We thus carry with us some momentum of previous weight updates, and the coefficient $\frac{k-2}{k+1}$ tells us to use more momentum later on (it converged to 1 for large $k$). However, the momentum term gives rise to oscillations in the objective function value, and the Nesterov's accelerated gradient is therefore not a strict descent method. Acceleration can nonetheless be combined with backtracking line search, as for the standard gradient descent method, but now from the intermediate parameter values $v$.


```python
# Gradient descent parameters
eta = 0.1
wUpdates = np.zeros([nIterations, 2])
wUpdates[0:2, :] = [-1.8, 0.85]
nllVals = np.zeros(nIterations)
nllVals[0:2] = negLogLikFun(X, wUpdates[0, :], y)

# Nesterov accelerated gradient with backtracking
for i in range(2, nIterations):
    v = wUpdates[i-1, :] + (i-2.)/(i+1)*(wUpdates[i-1, :]-wUpdates[i-2, :])
    delta = negLogLikDerFun(X, v, y)
    nllTmp = negLogLikFun(X, v , y)
    tooLarge = True
    # Backtracking loop
    while tooLarge:
        wUpdates[i, :] = v - eta*delta
        nllNew = negLogLikFun(X, wUpdates[i, :] , y)
        if nllNew > nllTmp - eta*np.dot(delta, delta)/2.:
            wUpdates[i, :] = wUpdates[i-1, :]
            eta *= beta
        else:
            tooLarge = False
            nllVals[i] = nllNew

# Plotting
grayColors = getGrayColors(3)
fig = plt.figure(figsize=(25, 7))
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
# Negative log-likelihood
contourHandle = ax.contourf(W0, W1, -llVals, 50, cmap=cm.coolwarm)
# Gradient
ax.quiver(W0q, W1q, W0derNorm, W1derNorm, angles='xy', scale_units='xy', scale=30, width=2e-3)
# Gradient descent progress
ax.plot(wUpdates[:, 0], wUpdates[:, 1], '-o', ms=6, color='white', lw=3)
# Axes settings
ax.set_xlabel('$w_0$')
ax.set_ylabel('$w_1$');
cBarHandle = plt.colorbar(contourHandle)
cBarHandle.set_label('negative log-likelihood')
```

81f4111e56b98e109647c2eeca7e74f75a85a818
ipynb
Jupyter Notebook
ConvexOptimization/Part 1, basic ideas and gradient descent.ipynb
ala-laurila-lab/jupyter-notebooks
c7fac1ee74af8e61832dad8536b223a205e79bf2
MIT
ConvexOptimization/Part 1, basic ideas and gradient descent.ipynb
ala-laurila-lab/jupyter-notebooks
c7fac1ee74af8e61832dad8536b223a205e79bf2
MIT
ConvexOptimization/Part 1, basic ideas and gradient descent.ipynb
ala-laurila-lab/jupyter-notebooks
c7fac1ee74af8e61832dad8536b223a205e79bf2
MIT
2019-02-21T17:03:39.000Z
2019-02-21T17:03:39.000Z
Qwen/Qwen-72B
1. YES
2. YES