Introduction

• Review of deep learning methods which provide insight into structured high-dimensional data.
• Deep learning uses layers of semi-affine input transformations to provide a predictive rule.
• Applying these layers of transformations leads to a set of attributes (or, features)
• To which probabilistic statistical methods can be applied.
• The best of both worlds can be achieved: scalable prediction rules fortified with uncertainty quantification, where sparse regularization finds the features.

Why DL?

• Deep learning provides a powerful pattern matching tool suitable for many AI applications.
• Image recognition and text analysis are probably two of the deep learning’s most successful.
• An image or a text as a high dimensional matrices and vectors, respectively.

What is DL?

• In general, a neural network can be described as follows. Let \(f_1 , \ldots , f_L\) be given univariate activation functions for each of the \(L\) layers.
• Activation functions are nonlinear transformations of weighted data.
• A semi-affine activation rule is then defined by \[ f_l^{W,b} = f_l \left ( \sum_{j=1}^{N_l} W_{lj} X_j + b_l \right ) = f_l ( W_l X_l + b_l )\,, \]

DL

\[ \hat{Y}(X) = F(X) = \left ( f_l^{W_1,b_1} \circ \ldots \circ f_L^{W_L,b_L} \right ) ( X)\,. \]

What is DL?

Similar to a classic basis decomposition, the deep approach uses univariate activation functions to decompose a high dimensional \(X\).

Let \(Z^{(l)}\) denote the \(l\)th layer, and so \(X = Z^{(0)}\). The final output \(Y\) can be numeric or categorical. The explicit structure of a deep prediction rule is then \[ \begin{aligned} \hat{Y} (X) & = W^{(L)} Z^{(L)} + b^{(L)} \\ Z^{(1)} & = f^{(1)} \left ( W^{(0)} X + b^{(0)} \right ) \\ Z^{(2)} & = f^{(2)} \left ( W^{(1)} Z^{(1)} + b^{(1)} \right ) \\ \ldots & \\ Z^{(L)} & = f^{(L)} \left ( W^{(L-1)} Z^{(L-1)} + b^{(L-1)} \right )\,. \end{aligned} \] Here \(W^{(l)}\) is a weight matrix and \(b^{(l)}\) are threshold or activation levels. Designing a good predictor depends crucially on the choice of univariate activation functions \(f^{(l)}\). The \(Z^{(l)}\) are hidden features which the algorithm will extract.

Deep approach employs hierarchical predictors comprising of a series of \(L\) nonlinear transformations applied to \(X\).

Deep Learning and Least Squares

• The goal is to find the optimal value of \(\theta\) that minimizes the expected loss function
• Given a training data set \(D = \{x_i,y_i\}_{i=1}^n\).
• The loss function is a measure of discrepancy between the true value of \(y\) and the predicted value \(f(x,\theta)\).

Model Estimation

Deep learning as well as a large number of statistical problems, can be expressed in the form \[ \min l(x) + \phi(x). \] In learning \(l(x)\) is the negative log-likelihood and \(\phi(x)\) is a penalty function that regularizes the estimate. From the Bayesian perspective, the solution to this problem may be interpreted as a maximum a posteriori \[p(y\mid x) \propto \exp\{-l(x)\}, ~ p(x) \propto \exp\{-\phi(x)\}.\]

The loss function

The loss function is usually defined as the negative log-likelihood function of the model. \[ l(\theta) = - \sum_{i=1}^n \log p(y_i | x_i, \theta), \] where \(p(y_i | x_i, \theta)\) is the conditional distribution of \(y_i\) given \(x_i\) and \(\theta\).

The loss function

Thus, in the case of regression, we have \[ y_i = f(x_i,\theta) + \epsilon, ~ \epsilon \sim N(0,\sigma^2), \] Thus, the loss function is \[ l(\theta) = - \sum_{i=1}^n \log p(y_i | x_i, \theta) = \sum_{i=1}^n (y_i - f(x_i, \theta))^2, \]

Optimisation

• Second order optimisation algorithms, such as BFGS used for traditional statistical models do not work well for deep learning models.
• The number of parameters a DL model has is large and estimating second order derivatives (Hessian or Fisher information matrix) becomes prohibitive from both computational and memory use standpoints.
  
• Instead, first order gradient descent methods are used for estimating parameters of a deep learning models.

Optimisation

The problem of parameter estimation (when likelihood belongs to the exponential family) is an optimisation problem

\[ \min_{\theta} l(\theta) := \dfrac{1}{n} \sum_{i=1}^n \log p(y_i, f(x_i, \theta)) \] where \(l\) is the negative log-likelihood of a sample, and \(\theta\) is the vector of parameters.

Optimisation: Gradient Descent

• Iterative algorithm
• Starts with an initial guess \(\theta^{0}\)
• Updates the parameter vector \(\theta\) at each iteration \(t\) as follows:

\[ \theta^{t+1} = \theta^t - \alpha_t \nabla l(\theta^t). \] - \(\alpha_t\) is the learning rate (step size). In DL model estimation it is a first class citizen parameter.

Gradient Descent for Linear Regression

\[ y_i = \theta_0 + \theta_1 x_i + \epsilon_i, \] or in matrix form \[ y = X \theta + \epsilon, \]

• \(\epsilon_i \sim N(0, \sigma^2)\)
• \(X\) is the design matrix

Regression

• Regression is simply a neural network which is wide and shallow.

data(iris)
y = iris$Petal.Length
x = cbind(rep(1,150),iris$Petal.Width)
# initialize theta
theta <- matrix(c(0, 0), nrow = 2, ncol = 1)
# learning rate
alpha <- 0.0001
# number of iterations
n_iter <- 1000
# gradient descent
for (i in 1:n_iter) {
    # compute gradient
    grad <- -2 * t(x) %*% (y - x %*% theta)
    # update theta
    theta <- theta - alpha * grad
}

Our gradient descent finds the following coefficients

Plot

Let’s plot the data and model estimated using gradient descent

plot(x[,2],y,pch=16, xlab="Petal.Width")
abline(theta[2],theta[1], lwd=3,col="red")

Compare

Let’s compare it to the standard estimation algorithm

m = lm(Petal.Length~Petal.Width, data=iris)

GD:

The values found by gradient descent are very close to the ones found by the standard OLS algorithm.

Logistic Regression

GLM with a logit link function \[ \log \left(\frac{p}{1-p}\right) = \theta_0 + \theta_1 x_1 + \ldots + \theta_p x_p, \]

The negative log-likelihood: \[ l(\theta) = - \sum_{i=1}^n \left[ y_i \log p_i + (1-y_i) \log (1-p_i) \right], \] where \(p_i = 1/\left(1 + \exp(-\theta_0 - \theta_1 x_{i1} - \ldots - \theta_p x_{ip})\right)\).

Derivative

The derivative of the negative log-likelihood function is \[ \nabla l(\theta) = - \sum_{i=1}^n \left[ y_i - p_i \right] \begin{pmatrix} 1 \\ x_{i1} \\ \vdots \\ x_{ip} \end{pmatrix}. \] In matrix notations, we have \[ \nabla l(\theta) = - X^T (y - p). \]

Code

Let’s implement gradient descent algorithm now.

y = ifelse(iris$Species=="setosa",1,0)
x = cbind(rep(1,150),iris$Sepal.Length)
lrgd  = function(x,y, alpha, n_iter) {
  theta <- matrix(c(0, 0), nrow = 2, ncol = 1)
  for (i in 1:n_iter) {
    # compute gradient
    p = 1/(1+exp(-x %*% theta))
    grad <- -t(x) %*% (y - p)
    # update theta
    theta <- theta - alpha * grad
  }
  return(theta)
}
theta = lrgd(x,y,0.005,20000)

The gradient descent parameters are

Plot

par(mar=c(4,4,0,0), bty='n')
p = 1/(1+exp(-x %*% theta))
plot(x[,2],y,pch=16, xlab="Sepal.Length")
lines(x[,2],p,type='p', pch=16,col="red")

Compare

Let’s compare it to the standard estimation algorithm

glm(y~x-1, family=binomial(link="logit"))

Call:  glm(formula = y ~ x - 1, family = binomial(link = "logit"))

Coefficients:
    x1      x2  
27.829  -5.176  

Degrees of Freedom: 150 Total (i.e. Null);  148 Residual
Null Deviance:      207.9 
Residual Deviance: 71.84    AIC: 75.84

Stochastic Gradient Descent

Stochastic gradient descent (SGD) is a variant of the gradient descent algorithm. \[ \nabla l(\theta) \approx \dfrac{1}{|B|} \sum_{i \in B} \nabla l(y_i, f(x_i, \theta)), \] where \(B \in \{1,2,\ldots,n\}\) is the batch samples from the data set. This method can be interpreted as gradient descent using noisy gradients, which are typically called mini-batch gradients with batch size \(|B|\).

Stochastic Gradient Descent

• In a small mini-batch regime, when \(|B| \ll n\) and typically \(|B| \in \{32,64,\ldots,1024\}\) it was shown that SGD converges faster than the standard gradient descent algorithm
• It does converge to minimizers of strongly convex functions
• It is robust to noise in the data
• It can avoid saddle-points, which is often an issue with deep learning log-likelihood functions.

Implementation

Now, we implement SGD for logistic regression and compare performance for different batch sizes

lrgd_minibatch  = function(x,y, alpha, n_iter, bs) {
  theta <- matrix(c(0, 0), nrow = 2, ncol = n_iter+1)
  n = length(y)
  for (i in 1:n_iter) {
    s = ((i-1)*bs+1)%%n
    e = min(s+bs-1,n)
    xl = x[s:e,]; yl = y[s:e]
    p = 1/(1+exp(-xl %*% theta[,i]))
    grad <- -t(xl) %*% (yl - p)
    # update theta
    theta[,i+1] <- theta[,i] - alpha * grad
  }
  return(theta)
}

Different Batch Sizes

Now run our SGD algorithm with different batch sizes.

set.seed(92) # kuzy
ind = sample(150)
y = ifelse(iris$Species=="setosa",1,0)[ind] # shuffle data
x = cbind(rep(1,150),iris$Sepal.Length)[ind,] # shuffle data
nit=200000
lr = 0.01
th1 = lrgd_minibatch(x,y,lr,nit,5)
th2 = lrgd_minibatch(x,y,lr,nit,15)
th3 = lrgd_minibatch(x,y,lr,nit,30)

Stochastic Gradient Descent

We run it with 2^{5} iterations and the learning rate of 0.01 and plot the values of \(\theta_1\) every 1000 iteration.

Important points

• We shuffle the data before using it.
• The larger the batch size, the smaller number of iterations are required for convergence
• The batch size does not change the number calculations required overall.
• Let’s look at the same plot, but scale the x-axis according to the amount of computations

Different Batch Sizes

plot(ind/1000,th1[1,ind], type='l', ylim=c(0,33), col=1, ylab=expression(theta[1]), xlab="Iteration")
abline(h=27.83, lty=2)
lines(ind/1000*3,th2[1,ind], type='l', col=2)
lines(ind/1000*6,th3[1,ind], type='l', col=3)
legend("bottomright", legend=c(5,15,30),col=1:3, lty=1, bty='n',title = "Batch Size")

Different Batch Sizes

There are several important considerations about choosing the batch size for SGD.

• The larger the batch size, the more memory is required to store the data.
• Parallelization is more efficient with larger batch sizes.
• Third, the larger the batch size, the less noise in the gradient.