Optimization¶

Unconstrained optimization¶

• Problem

• Maximization: $\min_{\mathbf{x}\in\mathbb{R}^d} -f(\mathbf{x})$

Golden section search¶

• $d=1$.

• Applicable to functions defined on an interval of the real line.

• Assumption: $f$ is continuous and strictly unimodal on the interval

\

A function $f: [a,b] \to \mathbb{R}$ is strictly unimodal if 
$$
f(\lambda x + (1-\lambda) y) < \max\{ f(x), f(y) \},
\quad
x, y \in [a, b],~x\neq y, ~\lambda \in (0, 1).
$$

• Idea: keep track of three points $a < x_1 < b$ such that $f(x_1) < \min\{f(a), f(b)\}$

• Choose a new point $x_0$ such that $x_0$ belongs to the longer of intervals $(a, x_1)$ and $(x_1, b)$.

• WLOG, suppose $a < x_0 < x_1$.

  • If $f(x_0) > f(x_1)$, then the minimum lies in the interval $(x_0, b)$. Hence set $(a, x_1, b) := (x_0, x_1, b)$.
  • Otherwise, the minimum lies in the interval $(a, x_1)$. Hence set $(a, x_1, b) := (a, x_0, x_1)$.

• Above rules do not say how to choose $x_0$. A possibility is to keep the lengths of the two possible candidate intervals the same: $b - x_0 = x_1 - a$, and the relative position of the intermediate points ($x_1$ to $(x_0, b)$; $x_0$ to $(a, x_1)$) the same.

• Let $\alpha = \frac{x_0 - a}{b - a}$ and $y = \frac{x_1 - x_0}{b -a}$. Then we require $\frac{b - x_1}{b - a} = \alpha$ and

By solving this equation, we have $\alpha = \frac{3 - \sqrt{5}}{2}$ or $1 - \alpha = \frac{\sqrt{5} - 1}{2}$, the golden section.

Convergence of golden section search is slow but sure to find the global minimum (if the assumptions are met).

Newton's Method¶

• First-order optimality condition

• Apply Newton-Raphson to $g(\mathbf{x})=\nabla f(\mathbf{x})$.

• Idea: iterative quadratic approximation.

• Second-order Taylor expansion of the objective function around the current iterate $\mathbf{x}^{(t)}$

and then minimize the quadratic approximation.

• To maximize the quadratic appriximation function, we equate its gradient to zero

which suggests the next iterate $$ \begin{eqnarray*} \mathbf{x}^{(t+1)} &=& \mathbf{x}^{(t)} - [\nabla^2 f(\mathbf{x}^{(t)})]^{-1} \nabla f(\mathbf{x}^{(t)}), \end{eqnarray*} $$ a complete analogue of the univariate Newton-Raphson for solving $g(x)=0$.

• Considered the gold standard for its fast (quadratic) convergence: $\mathbf{x}^{(t)} \to \mathbf{x}^{\star}$ and

In words, the estimate gets accurate by two decimal digits per iteration.

We call this naive Newton's method.

• Stability issue: naive Newton's iterate is not guaranteed to be a descent algorithm, i.e., that ensures $f(\mathbf{x}^{(t+1)}) \le f(\mathbf{x}^{(t)})$. It's equally happy to head uphill or downhill. Following example shows that the Newton iterate converges to a local maximum, converges to a local minimum, or diverges depending on starting points.

Practical Newton¶

• A remedy for the instability issue:

  1. approximate $\nabla^2 f(\mathbf{x}^{(t)})$ by a positive definite $\mathbf{H}^{(t)}$ (if it's not), and
  2. line search (backtracking) to ensure the descent property.

• Why insist on a positive definite approximation of Hessian? First define the Newton direction:

By first-order Taylor expansion, $$ \begin{align*} f(\mathbf{x}^{(t)} + s \Delta \mathbf{x}^{(t)}) - f(\mathbf{x}^{(t)}) &= \nabla f(\mathbf{x}^{(t)})^T s \Delta \mathbf{x}^{(t)} + o(s) \\ &= - s \cdot \nabla f(\mathbf{x}^{(t)})^T [\mathbf{H}^{(t)}]^{-1} \nabla f(\mathbf{x}^{(t)}) + o(s). \end{align*} $$ For $s$ sufficiently small, $f(\mathbf{x}^{(t)} + s \Delta \mathbf{x}^{(t)}) - f(\mathbf{x}^{(t)})$ is strictly negative if $\mathbf{H}^{(t)}$ is positive definite.
The quantity $\{\nabla f(\mathbf{x}^{(t)})^T [\mathbf{H}^{(t)}]^{-1} \nabla f(\mathbf{x}^{(t)})\}^{1/2}$ is termed the Newton decrement.

• In summary, a practical Newton-type algorithm iterates according to

where $\mathbf{H}^{(t)}$ is a positive definite approximation to $\nabla^2 f(\mathbf{x}^{(t)})$ and $s$ is a step size.

• For strictly convex $f$, $\nabla^2 f(\mathbf{x}^{(t)})$ is always positive definite. In this case, the above algorithm is called damped Newton. However, line search is still needed (at least for a finite number of times) to guarantee convergence.

Line search¶

• Line search: compute the Newton direction and search $s$ such that $f(\mathbf{x}^{(t)} + s \Delta \mathbf{x}^{(t)})$ is minimized.

• Note the Newton direction only needs to be calculated once. Cost of line search mainly lies in objective function evaluation.

• Full line search: $s = \arg\min_{\alpha} f(\mathbf{x}^{(t)} + \alpha \Delta \mathbf{x}^{(t)})$. May use golden section search.

• Approximate line search: step-halving ($s=1,1/2,\ldots$), Amijo rule, ...

• Backtracking line search (Armijo rule)

# Backtracking line search
    # given: descent direction ∆x, x ∈ domf, α ∈ (0,0.5), β ∈ (0,1).
    t <- 1.0
    while (f(x + t * delx) > f(x) + alpha * t * sum(gradf(x) * delx) {
        t <- beta * t
    }

The lower dashed line shows the linear extrapolation of $f$, and the upper dashed line has a slope a factor of α smaller. The backtracking condition is that $f$ lies below the upper dashed line, i.e., $0 \le t \le t_0$.

• How to approximate $\nabla^2 f(\mathbf{x})$? More of an art than science. Often requires problem specific analysis.

• Taking $\mathbf{H}^{(t)} = \mathbf{I}$ leads to the gradient descent method (see below).

Fisher scoring¶

• Consider MLE in which $f(\mathbf{x}) = -\ell(\boldsymbol{\theta})$, where $\ell(\boldsymbol{\theta})$ is the log-likelihood of parameter $\boldsymbol{\theta}$.

• Fisher scoring method: replace $- \nabla^2 \ell(\boldsymbol{\theta})$ by the expected Fisher information matrix

which is true under exchangeability of tne expectation and the differentiation (true for most common distributions).

Therefore we set $\mathbf{H}^{(t)}=\mathbf{I}(\boldsymbol{\theta}^{(t)})$ and obtain the Fisher scoring algorithm: 

• Combined with line search, a descent algorithm can be devised.

Example: logistic regression¶

• Binary data: response $y_i \in \{0,1\}$, predictor $\mathbf{x}_i \in \mathbb{R}^{p}$.

• Model: $y_i \sim $Bernoulli$(p_i)$, where

• MLE: density

• Log likelihood of the data $(y_i,\mathbf{x}_i)$, $i=1,\ldots,n$, and its derivatives are

(why the last line?)

• Therefore for this problem Newton's method == Fisher scoring:

where $$ \mathbf{z}^{(t)} = \mathbf{X} \beta^{(t)} + s[\mathbf{W}^{(t)}]^{-1} (\mathbf{y} - \mathbf{p}^{(t)}) $$ are the working responses. A Newton iteration is equivalent to solving a weighed least squares problem $$ \min_{\beta} \sum_{i=1}^n w_i (z_i - \mathbf{x}_i^T \beta)^2 $$ for which we know how to solve well. Thus the name IRLS (iteratively re-weighted least squares).

• Implication: if a weighted least squares solver is at hand, then logistic regression models can be fitted.
  • IRLS == Fisher scoring == Newton's method

Example: Poisson regression¶

• Count data: response $y_i \in \{0, 1, 2, \dotsc \}$, predictor $\mathbf{x}_i \in \mathbb{R}^{p}$.

• Model: $y_i \sim $Poisson$(\lambda_i)$, where

• MLE: density

• Log likelihood of the data $(y_i,\mathbf{x}_i)$, $i=1,\ldots,n$, and its derivatives are

(why the last line?)

• Therefore for this problem Newton's method == Fisher scoring:

where $$ \mathbf{z}^{(t)} = \mathbf{X} \beta^{(t)} + s[\mathbf{W}^{(t)}]^{-1} (\mathbf{y} - \boldsymbol{\lambda}^{(t)}) $$ are the working responses. A Newton iteration is again equivalent to solving a weighed least squares problem $$ \min_{\beta} \sum_{i=1}^n w_i (z_i - \mathbf{x}_i^T \beta)^2 $$ which leads to the IRLS.

• Implication: if a weighted least squares solver is at hand, then Poisson regression models can be fitted.
  • IRLS == Fisher scoring == Newton's method

Homework: repeat this using the Armijo rule.

Generalized Linear Models (GLM)¶

That IRLS == Fisher scoring == Newton's method for both logistic and Poisson regression is not a coincidence. Let's consider a more general class of generalized linear models (GLM).

Exponential families¶

• Random variable $Y$ belongs to an exponential family if the density

* $\eta$: natural parameter.  
* $\phi>0$: dispersion parameter.  
* Mean: $\mu= b'(\eta)$. When $b'(\cdot)$ is invertible, function $g(\cdot)=[b']^{-1}(\cdot)$ is called the canonical link function.
* Variance $\mathbf{Var}{Y}=b''(\eta)a(\phi)$.

• For example, if $Y \sim \text{Ber}(\mu)$, then

Hence $$ \eta = \log \frac{\mu}{1-\mu}, \quad \mu = \frac{e^{\eta}}{1+e^{\eta}}, \quad b(\eta) = -\log (1-\mu) = \log(1+e^{\eta}) $$ Hence $$ b'(\eta) = \frac{e^{\eta}}{1+e^{\eta}} = g^{-1}(\eta). $$ as above.

Generalized linear models¶

GLM models the conditional distribution of $Y$ given predictors $\mathbf{x}$ through the conditional mean $\mu = \mathbf{E}(Y|\mathbf{x})$ via a strictly increasing link function $$ g(\mu) = \mathbf{x}^T \beta, \quad \mu = g^{-1}(\mathbf{x}^T\beta) = b'(\eta) $$

From these relations we have (assuming no overdispersion, i.e., $a(\phi)\equiv 1$) $$ \mathbf{x}^T d\beta = g'(\mu)d\mu, \quad d\mu = b''(\eta)d\eta, \quad b''(\eta) = \mathbf{Var}[Y] = \sigma^2, \quad d\eta = \frac{1}{b''(\eta)}d\mu = \frac{1}{b''(\eta)g'(\mu)}\mathbf{x}^T d\beta. $$ Therefore $$ d\mu = \frac{1}{g'(\mu)}\mathbf{x}^T d\beta . $$

Then, after some workout with matrix calculus, we have for $n$ samples:

• Score, Hessian, information

• Fisher scoring method:

IRLS with weights $w_i = [1/g'(\mu_i)]^2/\sigma_i^2$ and some working responses $z_i$.

• For canonical link, $\mathbf{x}^T\beta = g(\mu) =[b']^{-1}(\mu) = \eta$. The second term of Hessian vanishes because $d\eta=\mathbf{x}^Td\beta$ and $b''(\eta)=1/g'(\mu)$. The Hessian coincides with Fisher information matrix. IRLS == Fisher scoring == Newton's method. Hence MLE is a convex optimization problem.

• Non-canonical link, Fisher scoring != Newton's method, and MLE is in general a non-convex optimization problem.

  Example: Probit regression (binary response with probit link).

where $\Phi(\cdot)$ is the cdf of a standard normal. Exceptionally, this problem is a convex optimization problem.

• R implements the Fisher scoring method, aka IRLS, for GLMs in function glm().

Nonlinear regression - Gauss-Newton method¶

• Now we finally get to the problem Gauss faced in 1801!

Relocate the dwarf planet Ceres https://en.wikipedia.org/wiki/Ceres_(dwarf_planet) by fitting 24 observations to a 6-parameter (nonlinear) orbit. - In 1801, Jan 1 -- Feb 11 (41 days), astronomer Piazzi discovered Ceres, which was lost behind the Sun after observing its orbit 24 times. - Aug -- Sep, futile search by top astronomers; Laplace claimed it unsolvable. - Oct -- Nov, Gauss did calculations by method of least squares, sent his results to astronomer von Zach. - Dec 31, von Zach relocated Ceres according to Gauss’ calculation.

• Nonlinear least squares (curve fitting):

For example, $y_i =$ dry weight of onion and $x_i=$ growth time, and we want to fit a 3-parameter growth curve $$ \mu(x, \beta_1,\beta_2,\beta_3) = \frac{\beta_3}{1 + \exp(-\beta_1 - \beta_2 x)}. $$

• If $\mu_i$ is a linear function of $\beta$, i.e., $\mathbf{x}_i^T\beta$, then NLLS reduces to the usual least squares.

• "Score" and "information matrices"

where $\mathbf{J}(\beta)^T = [\nabla \mu_1(\mathbf{x}_1,\beta), \ldots, \nabla \mu_n(\mathbf{x}_n,\beta)] \in \mathbb{R}^{p \times n}$.

• Gauss-Newton (= "Fisher scoring method") uses $\mathbf{I}(\beta)$, which is always positive semidefinite.

• Justification
  1. Residuals $y_i - \mu_i(\mathbf{x}_i, \beta)$ are small, or $\mu_i$ is nearly linear;
  2. Model the data as $y_i \sim N(\mu_i(\mathbf{x}_i, \beta), \sigma^2)$, where $\sigma^2$ is assumed to be known. Then,

Thus
$$
    \mathbf{E}[-\nabla^2 \ell(\beta)] = \sum_{i=1}^n \nabla \mu_i(\beta) \nabla \mu_i(\beta)^T.
$$

• In fact Gauss invented the Gaussian (normal) distribution to justify this method!

Gradient descent method¶

• Iteration:

for some step size $\gamma_t$. This is a special case of the Netwon-type algorithms with $\mathbf{H}_t = \mathbf{I}$ and $\Delta\mathbf{x} = \nabla f(\mathbf{x}^{(t)})$.

• Idea: iterative linear approximation (first-order Taylor series expansion).

and then minimize the linear approximation within a compact set: let $\Delta\mathbf{x} = \mathbf{x}^{(t+1)}-\mathbf{x}^{(t)}$ to choose. Solve $$ \min_{\|\Delta\mathbf{x}\|_2 \le 1} \nabla f(\mathbf{x}^{(t)})^T\Delta\mathbf{x} $$ to obtain $\Delta\mathbf{x} = -\nabla f(\mathbf{x}^{(t)})/\|\nabla f(\mathbf{x}^{(t)})\|_2 \propto -\nabla f(\mathbf{x}^{(t)})$.

• Step sizes are chosen so that the descent property is maintained (e.g., line search).

  • Step size must be chosen at any time, since unlike Netwon's method, first-order approximation is always unbounded below.

• Pros

  • Each iteration is inexpensive.
  • No need to derive, compute, store and invert Hessians; attractive in large scale problems.

• Cons

  • Slow convergence (zigzagging).
  
  • Do not work for non-smooth problems.

Convergence¶

• If

  1. $f$ is convex and differentiable over $\mathbb{R}^d$;
  2. $f$ has $L$-Lipschitz gradients, i.e., $\|\nabla f(\mathbf{x}) - \nabla f(\mathbf{y})\|_2 \le L\|\mathbf{x} - \mathbf{y} \|_2$;
  3. $p^{\star} = \inf_{\mathbf{x}\in\mathbb{R}^d} f(x) > -\infty$ and attained at $\mathbf{x}^{\star}$,

  then

for constant step size ($\gamma_t = \gamma \in (0, 1/L]$), a *sublinear* convergence. Similar upper bound with line search.

• In general, the best we can obtain from gradient descent is linear convergence (cf. quadratic convergence of Newton).

Examples¶

• Least squares: $f(\beta) = \frac{1}{2}\|\mathbf{X}\beta - \mathbf{y}\|_2^2 = \frac{1}{2}\beta^T\mathbf{X}^T\mathbf{X}\beta - \beta^T\mathbf{X}^T\mathbf{y} + \frac{1}{2}\mathbf{y}^T\mathbf{y}$

  • Gradient $\nabla f(\beta) = \mathbf{X}^T\mathbf{X}\beta - \mathbf{X}^T\mathbf{y}$ is $\|\mathbf{X}^T\mathbf{X}\|_2$-Lipschitz. Hence $\gamma \in (0, 1/\sigma_{\max}(\mathbf{X})^2]$ guarantees descent property and convergence.

• Logistic regression: $f(\beta) = -\sum_{i=1}^n \left[ y_i \mathbf{x}_i^T \beta - \log (1 + e^{\mathbf{x}_i^T \beta}) \right]$

  • Gradient $\nabla f(\beta) = - \mathbf{X}^T (\mathbf{y} - \mathbf{p})$ is $\frac{1}{4}\|\mathbf{X}^T\mathbf{X}\|_2$-Lipschitz.

  • May be unbounded below, e.g., all the data are of the same class.

  • Adding a ridge penalty $\frac{\rho}{2}\|\beta\|_2^2$ makes the problem identifiable.

• Poisson regression: $f(\beta) = -\sum_{i=1}^n \left[ y_i \mathbf{x}_i^T \beta - \exp(\mathbf{x}_i^T \beta) - \log(y_i!) \right]$

  • Gradient

is *not* Lipschitz (HW). What is the consequence?

Recall

Coordinate descent¶

• Idea: coordinate-wise minimization is often easier than the full-dimensional minimization

• Similar to the Gauss-Seidel method for solving linear equations.

• Block descent: a vector version of coordinate descent

• Q: why objective value converges?

  • A: descent property

• Caution: even if the objective function is convex, coordinate descent may not converge to the global minimum and converge to a suboptimal point.