Outline

• Univariate Optimization
  • Optimization in R
  • Gradient decent
  • Newton’s method
• Multivariate Optimization
  • Gradient decent
  • Newton’s method
  • Coordinate decent

Optimization

• In mathematics, “optimization” or “mathematical programming” refers to the selection of a best element (with regard to some criterion) from some set of available alternatives.

• A typical optimization problem usually consists of maximizing or minimizing a real function (objective function) by systematically choosing input values under certain constraint.

• “Convex programming” studies the case when the objective function is convex (minimization) or concave (maximization) and the constraint set is convex.

Convex function

Convex function properties:

• Second derivative always greater than 0
• Tangent lines always underestimate the function
• \(\forall t \in (0,1)\), \(tf(a) + (1 - t)f(b) > f(ta + (1 - t)b)\)

• Convex optimization always leads to a global minimizer
• Non-convex optimization usually leads to a local minimizer

Optimization functions in R

• optimize(): One dimensional optimization, no gradient or Hessian
• optim(): General purpose optimization, five possible methods, gradient optional
• constrOptim(): Minimization of a function subject to linear inequality constraints, gradient optional
• nlm(): Non-linear minimization, can optionally include the gradient and hessian of the function as attributes of the objective function
• nlminb(): Minimization using PORT routines, can optionally include the gradient and Hessian of the objective function as additional arguments

Univariate optimization example

Suppose our objective function is \[f(x) = e^x + x^4\] What is the minimum value of \(f(x)\), what is the correponding \(x\)?

f <- function(x) exp(x) + x^4
curve(f, from = -1, to = 1)

Use R optimize() function to perform the optimization

f <- function(x) exp(x) + x^4
xmin <- optimize(f, interval = c(-10, 10))
xmin

## $minimum
## [1] -0.5282468
## 
## $objective
## [1] 0.6675038

curve(f, from = -1, to = 1)
points(xmin$minimum, xmin$objective, col="red", pch=19)

Gradient decent

• Consider an unconstrained minimization of \(f\), we want to find \(x^*\) such that \[f(x^*) = \min_{x\in \mathbb{R}} f(x)\]
• Another way to formulate the problem is to find \(x^*\) such that \[x^* = \arg_{x\in \mathbb{R}} \min f(x)\]

Gradient decent: choose initial \(x^{(0)} \in \mathbb{R}\), repeat: \[x^{(k)} = x^{(k - 1)} - t\times f'(x^{(k - 1)}), k = 1,2,3,\ldots\]

• Until \(|x^{(k)} - x^{(k - 1)}| < \varepsilon\), (e.g. \(\varepsilon = 10^{-6}\)).
• \(f'(x^{(k - 1)})\) is the gradient (derivative of \(f\)) evaluated at \(x^{(k - 1)}\).
• \(t\) is a step size for gradient decent procedure.

Interpretation

• At each iteration \(k\), consider Taylor expansion at \(x^{(k - 1)}\): \[f(y) = f(x^{(k - 1)}) + f'(x^{(k - 1)}) \times (y - x^{(k - 1)}) + \frac{1}{2} f''(x^{(k - 1)}) (y - x^{(k - 1)})^2 + \ldots\]
• Quadratic approximation:
  • Replacing \(f''(x^{(k - 1)})\) with \(\frac{1}{t}\)
  • Ingore higher order terms

\[f(y) \approx f(x^{(k - 1)}) + f'(x^{(k - 1)}) \times (y - x^{(k - 1)}) + \frac{1}{2t} (y - x^{(k - 1)})^2\]

• Minimizing the quadratic approximation. \[x^{(k)} = \arg_y \min f(y)\]

  • Set \(f'(x^{(k)}) = 0 \Leftrightarrow f'(x^{(k - 1)}) + \frac{1}{t} (x^{(k)} - x^{(k - 1)}) = 0 \Leftrightarrow x^{(k)} = x^{(k - 1)} - t\times f'(x^{(k - 1)})\)

  • This is exactly the same as the gradient decent procedure.

Visual interpretation

• iterate 1: \(x^{(1)} = x^{(0)} - t\times f'(x^{(0)})\)
  • quadratic approximation at \(x_0\),
  • minimizer of the quadratic approximation is \(x_1\)
• iterate 2: \(x^{(2)} = x^{(1)} - t\times f'(x^{(1)})\)
  • quadratic approximation at \(x_1\),
  • minimizer of the quadratic approximation is \(x_2\)
• …

Gradient decent on our motivating example

f <- function(x) exp(x) + x^4 ## original function
g <- function(x) exp(x) + 4 * x^3 ## gradient function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 0.8; x_old <- 0; t <- 0.1; k <- 0; error <- 1e-6
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error){ ## there is a scale issue. 
  ## E.g., if x = 1e-7, x_old = 1e-6, relative change is large, absolute change is small
  ## We can monitor abs(x - x_old) / (abs(x) + abs(x_old))
  k <- k + 1
  x_old <- x
  x <- x_old - t*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
}

trace_grad <- trace
print(trace)

##  [1]  0.80000000  0.37264591  0.20678990  0.08028037 -0.02828567 -0.12548768
##  [7] -0.21290391 -0.28986708 -0.35496119 -0.40719158 -0.44673749 -0.47504573
## [13] -0.49435026 -0.50702281 -0.51511484 -0.52018513 -0.52332289 -0.52524940
## [19] -0.52642641 -0.52714331 -0.52757914 -0.52784381 -0.52800441 -0.52810183
## [25] -0.52816091 -0.52819673 -0.52821844 -0.52823161 -0.52823959 -0.52824443
## [31] -0.52824736 -0.52824914 -0.52825021 -0.52825087

A non-convex example

• If the objective function has multiple local minimums.
• Gradient will still work, but may fall into a local minimum (might be global minimum).
• kmeans algorithm is a non-convex minimization problem.

Discussion, will gradient decent algorithm always converge?

No, see the example

• Minimize \(f(x) = x^2\)
• Initial point \(x_0 = 1\)
• stepsize \(t = 1.1\)

f <- function(x) x^2 ## original function
g <- function(x) 2*x ## gradient function
curve(f, from = -10, to = 10, lwd=2) ## visualize the objective function

x <- x0 <- 1; x_old <- 0; t <- 1.1; k <- 0; error <- 1e-6
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error & k < 30){ ## there is a scale issue. 
  ## E.g., if x = 1e-7, x_old = 1e-6, relative change is large, absolute change is small
  ## We can monitor abs(x - x_old) / (abs(x) + abs(x_old))
  k <- k + 1
  x_old <- x
  x <- x_old - t*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
  segments(x0=x, y0 = f(x), x1 = x_old, y1 = f(x_old))
}

print(trace)

##  [1]    1.000000   -1.200000    1.440000   -1.728000    2.073600   -2.488320
##  [7]    2.985984   -3.583181    4.299817   -5.159780    6.191736   -7.430084
## [13]    8.916100  -10.699321   12.839185  -15.407022   18.488426  -22.186111
## [19]   26.623333  -31.948000   38.337600  -46.005120   55.206144  -66.247373
## [25]   79.496847  -95.396217  114.475460 -137.370552  164.844662 -197.813595
## [31]  237.376314

Change to another stepsize

• Minimize \(f(x) = x^2\)
• Initial point \(x_0 = 1\)
• stepsize \(t = 0.2\)

f <- function(x) x^2 ## original function
g <- function(x) 2*x ## gradient function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 1; x_old <- 0; t <- 0.2; k <- 0; error <- 1e-6
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error & k < 30){
  k <- k + 1
  x_old <- x
  x <- x_old - t*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
  segments(x0=x, y0 = f(x), x1 = x_old, y1 = f(x_old))
}

print(trace)

##  [1] 1.000000e+00 6.000000e-01 3.600000e-01 2.160000e-01 1.296000e-01
##  [6] 7.776000e-02 4.665600e-02 2.799360e-02 1.679616e-02 1.007770e-02
## [11] 6.046618e-03 3.627971e-03 2.176782e-03 1.306069e-03 7.836416e-04
## [16] 4.701850e-04 2.821110e-04 1.692666e-04 1.015600e-04 6.093597e-05
## [21] 3.656158e-05 2.193695e-05 1.316217e-05 7.897302e-06 4.738381e-06
## [26] 2.843029e-06 1.705817e-06 1.023490e-06

Change to another stepsize

• Minimize \(f(x) = x^2\)
• Initial point \(x_0 = 1\)
• stepsize \(t = 0.01\)

f <- function(x) x^2 ## original function
g <- function(x) 2*x ## gradient function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 1; x_old <- 0; t <- 0.01; k <- 0; error <- 1e-6
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error & k < 30){
  k <- k + 1
  x_old <- x
  x <- x_old - t*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
  segments(x0=x, y0 = f(x), x1 = x_old, y1 = f(x_old))
}

print(trace)

##  [1] 1.0000000 0.9800000 0.9604000 0.9411920 0.9223682 0.9039208 0.8858424
##  [8] 0.8681255 0.8507630 0.8337478 0.8170728 0.8007314 0.7847167 0.7690224
## [15] 0.7536419 0.7385691 0.7237977 0.7093218 0.6951353 0.6812326 0.6676080
## [22] 0.6542558 0.6411707 0.6283473 0.6157803 0.6034647 0.5913954 0.5795675
## [29] 0.5679762 0.5566167 0.5454843

How to select stepsize \(t\), backtracking line search

How to select stepsize \(t\), backtracking line search

• Linear approximation (first order Taylor expansion) at the current point \(x\). \[f(x + t\Delta x) = f(x) + tf'(x)\Delta x\]

• Decrease the slope of the approximation by \(\alpha\)

\[f(x) + t\alpha f'(x)\Delta x\]

• A sufficient condition for a proposed \(t\) to guarantee decent (decrease objective function): \[f(x + t\Delta x) < f(x) + t\alpha f'(x)\Delta x\]

• Plugin \(\Delta x = -f'(x)\) \[f(x - tf'(x)) < f(x) - t\alpha (f'(x))^2\]

Backtracking line search Algorithm

• Fix parameter \(0 < \beta < 1\) and \(0 < \alpha \le 1/2\)
• At each gradient decent iteration, start with \(t=1\), while \[f(x - t\times f'(x)) > f(x) - \alpha t \times (f'(x))^2\] Update \(t = \beta t\)
• When the above interation stops, use \(t\) as the step size

Note: In each iteration, you may always want to initialize \(t = 1\) before backtracking line search.

Implementation of back-tracking for gradient decent on the motivating example

• Minimize \(f(x) = \exp(x) + x^4\)
• gradient \(g(x) = \exp(x) + 4 x^3\)
• Initial point \(x_0 = 1\)
• initial stepsize \(t_0 = 1\), \(\alpha = 1/3\), \(\beta = 1/2\)
  • At each step, start with \(t = t_0\),
  • if \(f(x - tg(x)) > f(x) - \alpha t (g(x))^2\), set \(t = \beta t\)
  • otherwise, use \(t\) as step size

f <- function(x) exp(x) + x^4 ## original function
g <- function(x) exp(x) + 4 * x^3 ## gradient function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 1; x_old <- 0; t0 <- 1; k <- 0; error <- 1e-6; beta = 0.8; alpha <- 0.4
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error & k < 30){
  k <- k + 1
  x_old <- x
  
  ## backtracking
  t <- t0
  while(f(x_old - t*g(x_old)) > f(x_old) - alpha * t * g(x_old)^2){
    t <- t * beta
  }
  
  x <- x_old - t*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
  segments(x0=x, y0 = f(x), x1 = x_old, y1 = f(x_old))
}

print(trace)

## [1]  1.00000000  0.09828748 -0.61024238 -0.53353118 -0.52803658 -0.52825877
## [7] -0.52825165 -0.52825188

Exercise (will be on HW), solve for logistic regression

library(ElemStatLearn)
data2 <- prostate[,c("svi", "lcavol")]
str(data2)

## 'data.frame':    97 obs. of  2 variables:
##  $ svi   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ lcavol: num  -0.58 -0.994 -0.511 -1.204 0.751 ...

• Y: svi, binary 0/1.
• X: lcavol, continuous.
• logistic regression

\[\log \frac{E(Y|x)}{1 - E(Y|x)} = \beta_0 + x\beta_1\]

Exercise (will be on HW), solve for logistic regression

glm_binomial_logit <- glm(svi ~ lcavol, data = prostate,  family = binomial())
summary(glm_binomial_logit)

## 
## Call:
## glm(formula = svi ~ lcavol, family = binomial(), data = prostate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.79924  -0.48354  -0.21025  -0.04274   2.32135  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0296     1.0429  -4.823 1.42e-06 ***
## lcavol        1.9798     0.4543   4.358 1.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.35  on 96  degrees of freedom
## Residual deviance:  64.14  on 95  degrees of freedom
## AIC: 68.14
## 
## Number of Fisher Scoring iterations: 6

• We set \(\beta_0 = -5.0296\).
• Try to optimize \(\beta_1\), the coefficient for lcavol using gradient decent method
• compare with glm result.
• compare with the result using optimize() function.

Exercise (will be on HW), some hints

• density function: \[f(y) = p(x)^y(1 - p(x))^{1 - y}\]

• likelihood function: \[L(\beta_0, \beta_1) = \prod_{i = 1}^n p(x_i)^{y_i}(1 - p(x_i))^{1 - y_i}\]

• logistic regression, logit link \[\log \frac{p(x)}{1 - p(x)} = \beta_0 + x\beta_1\]

• log likelihood function:

\(\begin{aligned} l(\beta_0, \beta_1) &= \log\prod_{i = 1}^n p(x_i)^{y_i}(1 - p(x_i))^{(1 - y_i)} \\ & = \sum_{i=1}^n y_i \log p(x_i) + (1 - y_i) \log (1 - p(x_i)) \\ & = \sum_{i=1}^n \log (1 - p(x_i)) + \sum_{i=1}^n y_i \log \frac{p(x_i)}{1 - p(x_i)} \\ & = \sum_{i=1}^n -\log(1 + \exp(\beta_0 + x_i \beta_1)) + \sum_{i=1}^n y_i (\beta_0 + x_i\beta_1)\\ \end{aligned}\)

Newton’s method

• For unconstrained, smooth univariate convex optimization \[\min f(x)\] where \(f\) is convex, twice differentable, \(x \in \mathbb{R}\), \(f \in \mathbb{R}\). For gradient descent, start with intial value \(x^{(0)}\) and repeat the following \((k = 1,2,3,\ldots)\) until converge \[x^{(k)} = x^{(k - 1)} - t_k f'(x^{(k - 1)})\]

• For Newton’s method, start with intial value \(x^{(0)}\) and repeat the following \((k = 1,2,3,\ldots)\) until converge \[x^{(k)} = x^{(k - 1)} - (f''(x^{(k - 1)}))^{-1} f'(x^{(k - 1)})\] Where \(f''(x^{(k - 1)})\) is the second derivative of \(f\) at \(x^{(k - 1)}\). It is also refered as Hessian matrix for higher dimension (e.g. \(x \in \mathbb{R}^p\))

Newton’s method interpretation

• For gradient decent step at \(x\), we minimize the quadratic approximation \[f(y) \approx f(x) + f'(x)(y - x) + \frac{1}{2t}(y - x)^2\] over \(y\), which yield the update \(x^{(k)} = x^{(k - 1)} - tf'(x^{(k-1)})\)

• Newton’s method uses a better quadratic approximation: \[f(y) \approx f(x) + f'(x)(y - x) + \frac{1}{2}f''(x)(y - x)^2\] minimizing over \(y\) yield \(x^{(k)} = x^{(k - 1)} - (f''(x^{(k-1)}))^{-1}f'(x^{(k-1)})\)

Newton’s method on our motivating example

f <- function(x) exp(x) + x^4 ## original function
g <- function(x) exp(x) + 4 * x^3 ## gradient function
h <- function(x) exp(x) + 12 * x^2 ## Hessian function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 0.8; x_old <- 0; t <- 0.1; k <- 0; error <- 1e-6
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error){
  k <- k + 1
  x_old <- x
  x <- x_old - 1/h(x_old)*g(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
}
lines(trace, f(trace), lty = 2)

trace_newton <- trace
print(trace)

## [1]  0.8000000  0.3685707 -0.1665553 -0.8686493 -0.6362012 -0.5432407 -0.5285880
## [8] -0.5282520 -0.5282519

Compare Gradient decent with Newton’s method

f <- function(x) exp(x) + x^4 ## original function
par(mfrow=c(1,2))
title_grad <- paste("gradient decent, nstep =", length(trace_grad))
curve(f, from = -1, to = 1, lwd=2, main = title_grad) 
points(trace_grad, f(trace_grad), col=seq_along(trace_grad), pch=19)
lines(trace_grad, f(trace_grad), lty = 2)

title_newton <- paste("Newton's method, nstep =", length(trace_newton))
curve(f, from = -1, to = 1, lwd=2, main = title_newton) 
points(trace_newton, f(trace_newton), col=seq_along(trace_newton), pch=19)
lines(trace_newton, f(trace_newton), lty = 2)

Alternative interpretation about Newton’s method

• Aternative interpretation of Newton’s step: we seek a direction \(v\) so that \(f'(x + v) = 0\).
• By Taylor expansion: \[0 = f'(x + v) \approx f'(x) + f''(x)v\]
• Solve for \(v\), we have \(v = -(f''(x))^{-1} f'(x)\)

Backtracking line search for Newton’s method

• Usually pure Newton’s method works very well (do not necessary need to use backtracking)
• We have seen pure Newton’s method, which need not to decent at all iterations.
• In practice, we can use Newton’s method with backtracking search, which repeats \[x^{(k)} = x^{(k - 1)} - t(f''(x^{(k-1)}))^{-1}f'(x^{(k-1)})\]

• Note that the pure Newton’s method uses \(t = 1\)

• Algorithm:
  • Fix parameter \(0 < \beta < 1\) and \(0 < \alpha \le 1/2\)
  • At each Newton’s iteration \(x\), start with \(t=1\), while \[f(x + t \Delta x) > f(x) + \alpha t \times f'(x) \Delta x\] Update \(t = \beta t\)
  • In other word (set \(\Delta x = - (f''(x))^{-1} f'(x)\)), while \[f(x - t (f''(x))^{-1} f'(x)) > f(x) - \alpha t (f''(x))^{-1} (f'(x))^2\] Update \(t = \beta t\)

Implementation of back-tracking for Newton’s method

• Minimize \(f(x) = \exp(x) + x^4\)
• gradient \(g(x) = \exp(x) + 4 x^3\)
• Hessian \(h(x) = \exp(x) + 12 x^2\)
• Initial point \(x_0 = 1\)
• initial stepsize \(t_0 = 1\), \(\alpha = 1/3\), \(\beta = 1/2\)
  • At each step, start with \(t = t_0\),
  • if \(f(x - t g(x) / h(x) ) > f(x) - \alpha t (g(x))^2/h(x)\), set \(t = \beta t\)
  • otherwise, use \(t\) as step size

f <- function(x) exp(x) + x^4 ## original function
g <- function(x) exp(x) + 4 * x^3 ## gradient function
h <- function(x) exp(x) + 12 * x^2 ## gradient function
curve(f, from = -1, to = 1, lwd=2) ## visualize the objective function

x <- x0 <- 1; x_old <- 0; t0 <- 1; k <- 0; error <- 1e-6; beta = 0.8; alpha <- 0.4
trace <- x
points(x0, f(x0), col=1, pch=1)

while(abs(x - x_old) > error & k < 30){
  k <- k + 1
  x_old <- x
  
  ## backtracking
  t <- t0
  while(f(x_old - t*g(x_old)/h(x_old)) > f(x_old) - alpha * t * g(x_old)^2/h(x_old)){
    t <- t * beta
  }
  
  x <- x_old - t*g(x_old)/h(x_old)
  trace <- c(trace, x) ## collecting results
  points(x, f(x), col=k, pch=19)
  segments(x0=x, y0 = f(x), x1 = x_old, y1 = f(x_old))
}

print(trace)

## [1]  1.00000000  0.54354171  0.09465801 -0.63631402 -0.54326938 -0.52858926
## [7] -0.52825205 -0.52825187

Exercise (will be on HW), solve for logistic regression

library(ElemStatLearn)
glm_binomial_logit <- glm(svi ~ lcavol, data = prostate,  family = binomial())
summary(glm_binomial_logit)

## 
## Call:
## glm(formula = svi ~ lcavol, family = binomial(), data = prostate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.79924  -0.48354  -0.21025  -0.04274   2.32135  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0296     1.0429  -4.823 1.42e-06 ***
## lcavol        1.9798     0.4543   4.358 1.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.35  on 96  degrees of freedom
## Residual deviance:  64.14  on 95  degrees of freedom
## AIC: 68.14
## 
## Number of Fisher Scoring iterations: 6

• We set \(\beta_0 = -5.0296\).
• Try to optimize \(\beta_1\), the coefficient for lcavol using Newton’s method
• compare with gradient method

Comparison between gradient decent and Newton’s method

Multivariate Case

• What if \(\beta\) is not a scalar but instead a vector (i.e. \(\beta \in \mathbb{R}^p\))?

• Will gradient decent still work?
  • Yes

• How to calculate the gradient (derivative) for multivariate case?

• Will Newton’s method still work?
  • Yes

• How to calculate the Hessian matrix for multivariate case?

Gradient in multivariate case

• \(\beta = (\beta_1, \beta_2, \ldots, \beta_p)^\top \in \mathbb{R}^p\) is a \(p\)-dimensional column vector
• \(f(\beta) \in \mathbb{R}\) is a function which map \(\mathbb{R}^p \rightarrow \mathbb{R}\)
• Suppose \(f(\beta)\) is differentiable with respect to \(\beta\).

Then \[\nabla_\beta f(\beta) = \frac{\partial f(\beta)}{\partial \beta} = (\frac{\partial f(\beta)}{\partial \beta_1}, \frac{\partial f(\beta)}{\partial \beta_2}, \ldots, \frac{\partial f(\beta)}{\partial \beta_p})^\top \in \mathbb{R}^p\]

• Example \[f(x,y) = 4x^2 + y^2 + 2xy - x - y\]
  • \(\frac{\partial f}{\partial x} = 8x + 2y - 1\)
  • \(\frac{\partial f}{\partial y} = 2y +2x - 1\)
  • \(\nabla f = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})^\top = \binom{8x + 2y - 1}{2y +2x - 1} \in \mathbb{R}^2\)

Visualization of the example function

f <- function(x,y){
  4 * x^2 + y^2 + 2 * x * y - x - y 
}


xx <- yy <- seq(-20, 20, len=1000)
zz <- outer(xx, yy, f)
# Contour plots
contour(xx,yy,zz, xlim = c(-10,10), ylim = c(-20,20), nlevels = 10, levels = seq(1,100,length.out = 10), main="Contour Plot")

Hessian in multivariate case

• \(\beta = (\beta_1, \beta_2, \ldots, \beta_p)^\top \in \mathbb{R}^p\) is a \(p\)-dimensional column vector
• \(f(\beta) \in \mathbb{R}\) is a function which map \(\mathbb{R}^p \rightarrow \mathbb{R}\)
• Suppose \(f(\beta)\) is twice differentiable with respect to \(\beta\).

Then \[\nabla_\beta f(\beta) = \frac{\partial f(\beta)}{\partial \beta} = (\frac{\partial f(\beta)}{\partial \beta_1}, \frac{\partial f(\beta)}{\partial \beta_2}, \ldots, \frac{\partial f(\beta)}{\partial \beta_p})^\top \in \mathbb{R}^p\]

\[\Delta_\beta f(\beta) = \nabla^2_\beta f(\beta) = \nabla_\beta\frac{\partial f(\beta)}{\partial \beta} = (\nabla_\beta \frac{\partial f(\beta)}{\partial \beta_1}, \nabla_\beta \frac{\partial f(\beta)}{\partial \beta_2}, \ldots, \nabla_\beta \frac{\partial f(\beta)}{\partial \beta_p})^\top \]

\[\Delta_\beta f(\beta) = \begin{pmatrix} \frac{\partial^2 f(\beta)}{\partial \beta_1^2} & \frac{\partial^2 f(\beta)}{\partial \beta_1 \partial \beta_2} & \ldots & \frac{\partial^2 f(\beta)}{\partial \beta_1 \partial\beta_p}\\ \frac{\partial^2 f(\beta)}{\partial \beta_2 \partial \beta_1 } & \frac{\partial^2 f(\beta)}{\partial \beta_2^2} & \ldots & \frac{\partial^2 f(\beta)}{\partial \beta_2 \partial \beta_p} \\ \ldots &\ldots &\ldots &\ldots\\ \frac{\partial^2 f(\beta)}{\partial \beta_p \partial \beta_1} & \frac{\partial^2 f(\beta)}{\partial \beta_p \partial \beta_2} & \ldots & \frac{\partial^2 f(\beta)}{\partial\beta_p^2} \end{pmatrix} \]

Hessian in multivariate case example

\[f(x,y) = 4x^2 + y^2 + 2xy - x - y\]

• Gradient

  • \(\frac{\partial f}{\partial x} = 8x + 2y - 1\)
  • \(\frac{\partial f}{\partial y} = 2y +2x - 1\)

\[\nabla f = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})^\top = \binom{8x + 2y - 1}{2y +2x - 1} \in \mathbb{R}^2\]

• Hessian
  • \(\frac{\partial^2 f}{\partial x^2} = 8\)
  • \(\frac{\partial^2 f}{\partial y^2} = 2\)
  • \(\frac{\partial^2 f}{\partial x \partial y} = 2\)

\[\Delta f = \begin{pmatrix} 8 & 2 \\ 2 & 2 \end{pmatrix} \in \mathbb{R}^{2\times 2}\]

Gradient decent for multivariate case

Gradient decent: choose initial \(\beta^{(0)} \in \mathbb{R}^p\), repeat: \[\beta^{(k)} = \beta^{(k - 1)} - t\times \nabla_\beta f(\beta^{(k - 1)}), k = 1,2,3,\ldots\]

• Until \(\frac{\|\beta^{(k)} - \beta^{(k - 1)}\|_2}{\|\beta^{(k)} + \beta^{(k - 1)}\|_2} < \varepsilon\), (e.g. \(\varepsilon = 10^{-6}\)).
• \(\nabla_\beta f(\beta^{(k - 1)})\) is the gradient (derivative of \(f\)) evaluated at \(\beta^{(k - 1)}\).
• \(t\) is a step size for gradient decent procedure.

Gradient decent on our motivating example

f <- function(x,y){
  4 * x^2 + y^2 + 2 * x * y - x - y 
}

g <- function(x,y){
  c(8*x + 2*y - 1, 2*y +2*x - 1)
}

xx <- yy <- seq(-20, 20, len=1000)
zz <- outer(xx, yy, f)
# Contour plots
contour(xx,yy,zz, xlim = c(-10,10), ylim = c(-20,20), nlevels = 10, levels = seq(1,100,length.out = 10), main="Contour Plot")


x <- x0 <- 8; x_old <- 0; 
y <- y0 <- -10; y_old <- 0; 
curXY <- c(x,y)
preXY <- c(x_old, y_old)
t <- 0.1; k <- 0; error <- 1e-6
trace <- list(curXY)
points(x0,y0, col=1, pch=1)

l2n <- function(avec){
  sqrt(sum(avec^2))
}
diffChange <- function(avec, bvec){
  deltaVec <- avec - bvec
  sumVec <- avec + bvec
  l2n(deltaVec)/l2n(sumVec)
}

while(diffChange(curXY, preXY) > error){
  k <- k + 1
  preXY <- curXY
  curXY <- preXY - t*g(preXY[1], preXY[2])
  trace <- c(trace, list(curXY)) ## collecting results
  points(curXY[1], curXY[2], col=k, pch=19)
  segments(curXY[1], curXY[2], preXY[1], preXY[2])
}

print(k)

## [1] 97

How to select step size \(t\), backtracking line search

How to select step size \(t\), backtracking line search

• Linear approximation (first order Taylor expansion) at the current point \(x\). \[f(x + t\Delta x) = f(x) + t\nabla f(x)^\top\Delta x\]

• Decrease the slope of the approximation by \(\alpha\)

\[f(x) + t\alpha \nabla f(x)^\top\Delta x\]

• A sufficient condition for a proposed \(t\) to guarantee decent (decrease objective function): \[f(x + t\Delta x) < f(x) + t\alpha \nabla f(x)^\top\Delta x\]

• Plugin \(\Delta x = -f'(x)\) \[f(x - t\nabla f(x)) < f(x) - \alpha t \times \|(\nabla f(x)\|_2^2\]

• Algorithm:
  • Fix parameter \(0 < \gamma < 1\) and \(0 < \alpha \le 1/2\)
  • At each gradient decent iteration, start with \(t=1\), while \[f(\beta - t\times \nabla f(\beta)) > f(\beta) - \alpha t \times \|(\nabla f(\beta)\|_2^2\] Update \(t = \gamma t\)

Implementation of back-tracking for gradient decent

• Minimize \(f(x,y) = 4x^2 + y^2 + 2xy - x - y\)
• gradient \(\nabla f = \binom{8x + 2y - 1}{2y +2x - 1}\)
• Initial point \(x_0 = 8\); \(y_0 = -10\)
• initial step size \(t_0 = 1\), \(\alpha = 1/3\), \(\gamma = 1/2\)
  • At each step, start with \(t = t_0\),
  • if \(f(x - tg(x)) > f(x) - \alpha t \|g(x)\|_2^2\), set \(t = \gamma t\)
  • otherwise, use \(t\) as step size

f0 <- function(x,y){
  4 * x^2 + y^2 + 2 * x * y - x - y 
}

f <- function(avec){
  x <- avec[1]
  y <- avec[2]
  4 * x^2 + y^2 + 2 * x * y - x - y 
}

g <- function(avec){
  x <- avec[1]
  y <- avec[2]  
  c(8*x + 2*y - 1, 2*y +2*x - 1)
}

x <- y <- seq(-20, 20, len=1000)
z <- outer(x, y, f0)
# Contour plots
contour(x,y,z, xlim = c(-10,10), ylim = c(-20,20), nlevels = 10, levels = seq(1,100,length.out = 10), main="Contour Plot")


x <- x0 <- 8; x_old <- 0; 
y <- y0 <- -10; y_old <- 0; 
curXY <- c(x,y)
preXY <- c(x_old, y_old)
t0 <- 1; k <- 0; error <- 1e-6
alpha = 1/3
beta = 1/2
trace <- list(curXY)
points(x0,y0, col=1, pch=1)

l2n <- function(avec){
  sqrt(sum(avec^2))
}
diffChange <- function(avec, bvec){
  deltaVec <- avec - bvec
  sumVec <- avec + bvec
  l2n(deltaVec)/l2n(sumVec)
}

while(diffChange(curXY, preXY) > error){
  k <- k + 1
  preXY <- curXY
  
  ## backtracking
  t <- t0
  while(f(preXY - t*g(preXY)) > f(preXY) - alpha * t * l2n(g(preXY))^2){
    t <- t * beta
  }
  
  curXY <- preXY - t*g(preXY)
  trace <- c(trace, list(curXY)) ## collecting results
  points(curXY[1], curXY[2], col=k, pch=19)
  segments(curXY[1], curXY[2], preXY[1], preXY[2])
  
}

print(k)

## [1] 25

Exercise, solve for logistic regression

library(ElemStatLearn)
glm_binomial_logit <- glm(svi ~ lcavol, data = prostate,  family = binomial())
summary(glm_binomial_logit)

## 
## Call:
## glm(formula = svi ~ lcavol, family = binomial(), data = prostate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.79924  -0.48354  -0.21025  -0.04274   2.32135  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0296     1.0429  -4.823 1.42e-06 ***
## lcavol        1.9798     0.4543   4.358 1.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.35  on 96  degrees of freedom
## Residual deviance:  64.14  on 95  degrees of freedom
## AIC: 68.14
## 
## Number of Fisher Scoring iterations: 6

• Optimize \(\beta_1\), the coefficient for lcavol and the intercept \(\beta_0\) simultaneously using gradient decent method.

Newton’s method for multivariate case

• For unconstrained, smooth univariate convex optimization \[\min f(\beta)\] where \(f\) is convex, twice differentiable, \(\beta \in \mathbb{R}^p\), \(f \in \mathbb{R}\). For gradient descent, start with initial value \(\beta^{(0)}\) and repeat the following \((k = 1,2,3,\ldots)\) until converge \[\beta^{(k)} = \beta^{(k - 1)} - t_k \nabla f(\beta^{(k - 1)})\]

• For Newton’s method, start with initial value \(\beta^{(0)}\) and repeat the following \((k = 1,2,3,\ldots)\) until converge \[\beta^{(k)} = \beta^{(k - 1)} - (\Delta f(\beta^{(k - 1)}))^{-1} \nabla f(\beta^{(k - 1)})\] Where \(\Delta f(\beta^{(k - 1)}) \in \mathbb{R}^{p\times p}\) is the Hessian matrix of \(f\) at \(\beta^{(k - 1)}\).

Implement multivariate Newton’s method

f0 <- function(x,y){
  4 * x^2 + y^2 + 2 * x * y - x - y 
}

f <- function(avec){
  x <- avec[1]
  y <- avec[2]
  4 * x^2 + y^2 + 2 * x * y - x - y 
}

g <- function(avec){
  x <- avec[1]
  y <- avec[2]  
  c(8*x + 2*y - 1, 2*y +2*x - 1)
}

h <- function(avec){
  x <- avec[1]
  y <- avec[2]  
  res <- matrix(c(8,2,2,2),2,2)  ## Hessian function
  return(res)  
} 

x <- y <- seq(-20, 20, len=1000)
z <- outer(x, y, f0)
# Contour plots
contour(x,y,z, xlim = c(-10,10), ylim = c(-20,20), nlevels = 10, levels = seq(1,100,length.out = 10), main="Contour Plot")


x <- x0 <- 8; x_old <- 0; 
y <- y0 <- -10; y_old <- 0; 
curXY <- c(x,y)
preXY <- c(x_old, y_old)
t0 <- 1; k <- 0; error <- 1e-6
trace <- list(curXY)
points(x0,y0, col=1, pch=1)

l2n <- function(avec){
  sqrt(sum(avec^2))
}
diffChange <- function(avec, bvec){
  deltaVec <- avec - bvec
  sumVec <- avec + bvec
  l2n(deltaVec)/l2n(sumVec)
}

while(diffChange(curXY, preXY) > error){
  k <- k + 1
  preXY <- curXY
  
  curXY <- preXY - solve(h(curXY)) %*% g(preXY)
  trace <- c(trace, list(curXY)) ## collecting results
  points(curXY[1], curXY[2], col=k, pch=19)
  segments(curXY[1], curXY[2], preXY[1], preXY[2])
}

k ## why k is only 2?

## [1] 2

Exercise

\[f(x,y) = \exp(xy) + y^2 + x^4\]

What are the x and y such that \(f(x,y)\) is minimized?

• Use Newton’s method to solve this problem

Backtracking line search for Newton’s method

• Usually pure Newton’s method works very well (do not necessary need to use backtracking)
• We have seen pure Newton’s method, which need not to decent at all iterations.
• In practice, we can use Newton’s method with backtracking search, which repeats \[\beta^{(k)} = \beta^{(k - 1)} - t (\Delta f(\beta^{(k - 1)}))^{-1} \nabla f(\beta^{(k - 1)})\]

• Note that the pure Newton’s method uses \(t = 1\)

• Algorithm:
  • Fix parameter \(0 < \gamma < 1\) and \(0 < \alpha \le 1/2\)
  • At each Newton’s iteration \(x\), start with \(t=1\), while \[f(\beta + t \Delta \beta) > f(\beta) + \alpha t \times \nabla f(\beta)^\top \Delta \beta\] Update \(t = \gamma t\)
  • In other word (set \(\Delta \beta = - (\Delta f(\beta))^{-1} \nabla f(\beta)\)), while \[f(\beta - t \Delta f(\beta)^{-1} \nabla f(\beta)) > f(\beta) - \alpha t (\nabla f(\beta))^\top (\Delta f(\beta))^{-1} \nabla f(\beta)\] Update \(t = \gamma t\)

Exercise, solve for logistic regression

library(ElemStatLearn)
glm_binomial_logit <- glm(svi ~ lcavol, data = prostate,  family = binomial())
summary(glm_binomial_logit)

## 
## Call:
## glm(formula = svi ~ lcavol, family = binomial(), data = prostate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.79924  -0.48354  -0.21025  -0.04274   2.32135  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0296     1.0429  -4.823 1.42e-06 ***
## lcavol        1.9798     0.4543   4.358 1.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.35  on 96  degrees of freedom
## Residual deviance:  64.14  on 95  degrees of freedom
## AIC: 68.14
## 
## Number of Fisher Scoring iterations: 6

• Optimize \(\beta_1\), the coefficient for lcavol and the intercept \(\beta_0\) simultaneously using Newton’s method.

Coordinate decent

• Suppose \(\beta \in \mathbb{R}^P\) and \(f(\beta)\) is a mapping from \(\mathbb{R}^P\) to \(R\).
• We can use coordinate descent to find a minimizer.
• Start with an initial guess \(\beta^{(0)} = (\beta^{(0)}_1, \beta^{(0)}_2, \ldots, \beta^{(0)}_p)^\top\) and repeat the following:
  • \(\beta_1^{(k)} = \arg \min_{\beta_1} f(\beta_1, \beta_2^{(k-1)}, \beta_3^{(k-1)}, \ldots, \beta_p^{(k-1)})\)
  • \(\beta_2^{(k)} = \arg \min_{\beta_2} f(\beta_1^{(k)}, \beta_2, \beta_3^{(k-1)}, \ldots, \beta_p^{(k-1)})\)
  • \(\beta_3^{(k)} = \arg \min_{\beta_3} f(\beta_1^{(k)}, \beta_2^{(k)}, \beta_3, \ldots, \beta_p^{(k-1)})\)
  • \(\ldots\)
  • \(\beta_p^{(k)} = \arg \min_{\beta_p} f(\beta_1^{(k)}, \beta_2^{(k)}, \beta_3^{(k)}, \ldots, \beta_p)\)
• Continue with \(k=1,2,3,\ldots\) until converge.
  • Note that after \(\beta_j^{(k)}\) has been updated, we use the new value for future update.

Example on linear regression

• Consider linear regression problem:

\[ \min_{\beta \in \mathbb{R}^p} f(\beta) = \min_{\beta \in \mathbb{R}^p} \frac{1}{2} \|y - X\beta\|_2^2\] where \(y \in \mathbb{R}^n\) and \(X \in \mathbb{R}^{n \times p}\) with columns \(X_1\) (intercept), \(X_2\), \(\ldots\), \(X_p\).

• Minimizing over \(\beta_j\) while fixing all \(\beta_i\), \(i \ne j\): \[0 = \frac{\partial f(\beta)}{\partial \beta_j} = X_j^\top (X\beta - y) = X_j^\top (X_j\beta_j + X_{-j}\beta_{-j} - y)\]

• We can solve for the coordinate descent updating rule:

\[\beta_j = \frac{X_j^\top (y - X_{-j}\beta_{-j})}{X_j^\top X_j}\]

implement coordinate decent using prostate cancer data

library(ElemStatLearn)

y <- prostate$lcavol
x0 <- prostate[,-match(c("lcavol", "train"), colnames(prostate))]
x <- cbind(1, as.matrix(x0))

beta <- rep(0,ncol(x))
beta_old <- rnorm(ncol(x))

error <- 1e-6
k <- 0
while(diffChange(beta, beta_old) > error){
  k <- k + 1
  beta_old <- beta
  for(j in 1:length(beta)){
    xj <- x[,j]
    xj_else <- x[,-j]
    beta_j_else <- as.matrix(beta[-j])
    beta[j] <- t(xj) %*% (y - xj_else %*% beta_j_else) / sum(xj^2)
  }
}

print(beta)

## [1] -2.419610178 -0.025954949  0.022114840 -0.093134994 -0.152983715
## [6]  0.367195579  0.197095411 -0.007114941  0.565822653

## compare with lm
lm(lcavol ~ . - train, data=prostate)

## 
## Call:
## lm(formula = lcavol ~ . - train, data = prostate)
## 
## Coefficients:
## (Intercept)      lweight          age         lbph          svi  
##   -2.420613    -0.025925     0.022113    -0.093138    -0.152966  
##         lcp      gleason        pgg45         lpsa  
##    0.367180     0.197249    -0.007117     0.565826

Reference

http://www.stat.cmu.edu/~ryantibs/convexopt/