notebook.tex

\documentclass[8pt,letter]{article}
 
%% \usepackage[fleqn]{amsmath}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amsfonts,amsthm,bm}
\usepackage{breqn}
\usepackage{amsmath}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{tikz}
\usepackage[ruled,vlined,linesnumbered,lined,boxed,commentsnumbered]{algorithm2e}
\usepackage{siunitx}
\usepackage{graphicx}
\usepackage{subcaption}
%% \usepackage{datetime}
\usepackage{multirow}
\usepackage{multicol}
\usepackage{mathrsfs}
\usepackage{fancyhdr}
\usepackage{fancyvrb}
\usepackage{parskip} %turns off paragraph indent
\pagestyle{fancy}

\usepackage{xcolor}
\usepackage{mdframed}

\usetikzlibrary{arrows}

\DeclareMathOperator*{\argmin}{argmin}
\newcommand*{\argminl}{\argmin\limits}

\newcommand{\mathleft}{\@fleqntrue\@mathmargin0pt}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}} 
\newcommand{\N}{\mathbb{N}}
\newcommand{\ppartial}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\p}{\partial}
\newcommand{\te}[1]{\text{#1 }}
\newcommand{\norm}[1]{\|#1\|}

\setcounter{MaxMatrixCols}{20}

% remove excess vertical space for align equations
\setlength{\abovedisplayskip}{0pt}
\setlength{\belowdisplayskip}{0pt}
\setlength{\abovedisplayshortskip}{0pt}
\setlength{\belowdisplayshortskip}{0pt}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[theorem]
\newtheorem{lemma}[theorem]{Lemma}

% \newtheorem{mdtheorem}{Theorem}
% \newenvironment{theorem}
% {\begin{mdframed}[
%   backgroundcolor=green!10,
%   topline=false,
%   rightline=false,
%   bottomline=false,
%   leftline=false
%   ]\begin{mdtheorem}}
%     {\end{mdtheorem}\end{mdframed}}

\begin {document}

\lhead{Notes - Numerical Optimization, Bill (Yuan) Liu}

\begin{multicols*}{2}

  \section{General}

  \subsection{Wolfe conditions}

  Finding a step length $\alpha_k$ given a search direction $p_k$.
  
  \begin{itemize}
    
  \item Sufficient Decrease/ Armijo Condition
    
    Decrease proportional to $\alpha_k$ step length and directional derivative $\nabla f_k^T p_k$:
    \begin{align}
      f_{k+1} & \leq f_k + c_1 \alpha_k \nabla f_k^T p_k\\
      \nabla f_k^T p_k & \leq 0\ //by\ construction\\
      c_1 & \in (0,1)
    \end{align}
    
  \item Curvature Condition
    
    Reject extremely short steps:
    \begin{align}
      \nabla f_{k+1}^T p_k & \geq c_2 \nabla f_k^T p_k\\
      \frac{\partial \phi(\alpha_k)}{\partial \alpha_k} & \geq \frac{\partial \phi(0)}{\partial \alpha_k}\\
      &c_1,\alpha_k \in (0,1)\\
      &0<c_1<c_2<1
    \end{align}
    \begin{align*}
      \text{where}:& f_{k} = f(x_k)\\
      & f_{k+1} = f(x_k+\alpha_k p_k)\\
      & \phi(\alpha_k) = f(x_k + \alpha_k p_k)\\
      & \frac{\partial \phi(\alpha_k)}{\partial \alpha_k} = \nabla f_{k+1}^T p_k\\
      & \frac{\partial \phi(0)}{\partial \alpha_k} = \nabla f_{k}^T p_k //initial\ slope
    \end{align*}

  \end{itemize}
  
  \subsection{Strong Wolfe Condition}
  
  Modify curvature condition to reject large positive derivative:
  \begin{align}
    f_{k+1} & \leq f_k + c_1 \alpha_k \nabla f_k^T p_k\\
    | \nabla f_{k+1}^T p_k| & \leq c_2 |\nabla f_k^T p_k|\\
    0 & < c_1 < c_2 < 1
  \end{align}
  
  \subsection{Form of Search Direction}
  
  \begin{align}
    p_k = - B_k^{-1} \nabla f_k
  \end{align}
  
  $B_k$ symmetric, non-singular, postive definite $\implies$ $p_k$ is a descent direction:
  \begin{align}
    \nabla f_k^T(-B_k^{-1} \nabla f_k) < 0
  \end{align}

  \subsection{Goldstein Condition}
  
  \begin{align*}
    &f_k+(1-c) \alpha_k \nabla f_k^T p_k \leq f_{k+1} \leq f_k + c \alpha \nabla f_k^T p_k\\
    &c \in (0,\frac{1}{2})
  \end{align*}

  may miss minimizer of $f$, commonly used in Newtom algo

  \subsection{Back Tracking Line Search}

  Using sufficient decrease condition and back tracking for step length search:

  \begin{algorithm}[H]
    \SetKwInOut{Input}{$f, x, d, c_1, \alpha, \beta$}\SetKwInOut{Output}{$\alpha$}
    \Input{function, x, direction, gradient threshold, initial step length, contraction}
    \Output{found step length}
    \While{$f(x + \alpha d) > f(x) + c_1 \alpha \nabla f(x)^T d$}{
      $\alpha \leftarrow \alpha * \beta$\\
    }
    return $\alpha$
    \caption{Line Search\label{Algo_LineSearch}}
  \end{algorithm}
  

  \vfill\null
  \pagebreak
  
  \section{Conjugate Gradient}

  \subsection{linear method}
  Assuming unconstrained problem with strict convex quadratic objective function:
  \begin{align*}
    \frac{1}{2} x^TAx-b^Tx, A \succ 0, A^T=A
  \end{align*}


  $\nabla(\frac{1}{2} x^TAx-b^Tx) = Ax-b$, thus $\min_x x^TAx - b^Tx$ transformed to solving $Ax-b=0$. 

  let $x_{k+1} = x_k + \alpha_k p_k$, solve for $\alpha$:
  \begin{align*}
    %%     &\frac{\partial}{\partial \alpha} (\frac{1}{2}(x_k + \alpha_k p_k)^T A (x_k + \alpha_k p_k) - b^T(x_k + \alpha_k p_k)) = 0\\
    &Ax_{k+1}-b=0\\
    &A(x_k+\alpha_k p_k) - b=0\\
    &\alpha_k A p_k = b - Ax_k\\
    &r_k = Ax-b\\
    &\alpha_k A p_k = -r_k\\
    &\alpha_k p_k^T A p_k = - p_k^T r_k\\
    &\alpha_k = - \frac{p_k^T r_k}{p_k^T A p_k}
  \end{align*}
  
  \subsection{Conjugate Direction}
  
  Enforce by construction for search directions linearly independent wrt. A.:
  \begin{align*}
    (\forall i\not=j) p_i^T A p_j = 0
  \end{align*}
  Properties:
  \begin{itemize}
  \item Residual elimnated one direction at a time, resulting in max of n iterations.
  \item Optimal if Hessian is diagonal, if not can try preconditioning.
  \item Current residual is orthogonal to all previous search directions.
  \item Any set of conjugate directions can be used:
    \begin{itemize}
    \item eigenvectors
    \item Gram-Schmidt
    \item conjugate gradient algo.
    \end{itemize}
  \end{itemize}
  
  \subsubsection{Termination Steps}
  \begin{theorem}[Conjugate Direction Termination]\label{thm:CD_Term}
    Conjugate direction algorithm converges to solution $x \in \R^n$ of linear system in n steps. Theorem [T5.1] in Num. Opt. book.
  \end{theorem}
  \begin{proof}
    \begin{align*}
      &\text{given:}\\
      &x_{k+1} = x_k + \alpha_k p_k\\
      &\alpha_k = - \frac{r_k^T p_k}{p_k^T A p_k}\\
      &(\forall i \neq j) p_i^T A p_j = 0\ \text{[A conjugate]}\\
      &\text{then:}\\
      &x^*-x_0 = \sum_{i=0}^{n-1} \sigma_i p_i\\
      &p_k^T A(x^*-x_0) = p_k^T A(\sum_{i=0}^{n-1} \sigma_i p_i)\\
      &\sigma_i = \frac{p_k^T A(x^*-x_0)}{p_k^T A p_k}\ \text{(from conjugacy)}\\
      &x_k = x_0 + \sum_{i=0}^{k-1} \alpha_i p_i\\
      &p_k^T A x_k = p_k^T A (x_0+\alpha_0 p_0 + ..)\\
      &p_k^T A x_k = p_k^T A x_0\ \text{(from conjugacy)}\\
      &p_k^T A (x_k-x_0) = 0\\
      &p_k^T A(x^*-x_0) = p_k^T A(x^*-x_k)\\
      &p_k^T A(x^*-x_0) = p_k^T(b-Ax_k) = -p_k^T r_k\\
      &\sigma_i = \frac{-p_k^T r_k}{p_k^T A p_k} = \alpha_k\\
      &\alpha_k\ \text{is 1D minimizer of kth coordinate by construction}
    \end{align*}
    Conjugate direction algo. terminates in n steps.
  \end{proof}

  If A is not diagonal, can transform coordinates:
  \begin{align*}
    S\hat{x} &= x\\
    S &=[p_0\ p_1\ ..\ p_{n-1}],\ (\forall i \neq j)\ p_i A p_j =0\\
    \phi(x) &= \frac{1}{2} x^T A x - b^T x\\
    \phi(\hat{x}) &= \frac{1}{2} \hat{x}^T S^T A S \hat{x} - (S^T b)^T \hat{x}\\
    S^T A S\ diagonal & \implies \text{can optimize one coordinate at a time}\\
  \end{align*}
  
  \vfill\null
  \columnbreak
  
  \subsubsection{Expanding subspace minimizer}
  Using conjugate directions to generate sequence $\{x\}$, then:\\
  $r_k^T p_i = 0, \forall i < k$, $x_k$ is minimizer of $\frac{1}{2} x^TAx - b^Tx$ over $\{ x | x= x_0 + span\{p_0,...p_{k-1}\}$
  \begin{proof}
    \begin{align*}
      &\tilde{x} = x_0 + \sum_i \sigma_i p_i\\
      &\tilde{x} \text{ minimizes over }\{x_0 + span\{ p_0, ... p_{k-1}\}\} \iff r(\tilde{x})^Tp_i = 0\\
      &h(\sigma) = \phi(\tilde{x})\\
      &\phi(x) =\frac{1}{2}x^TAx-b^Tx\\
      &\text{h is also strictly convex quadratic,}\\
      &\text{with unique }\sigma^*\text{ satisfying:}\\
      &\frac{\partial h(\sigma^*)}{\partial \sigma_i} = 0, i=[0,k-1]\\
      &\frac{\partial h(\sigma^*)}{\partial \sigma_i} = \nabla \phi(\tilde{x})^T p_i = 0, i = [0,k-1]\\
      &\nabla \phi(x) = Ax-b=r\\
      &r(\tilde{x})^T p_i = 0, i=[0,k-1]\\
    \end{align*}  
  \end{proof}
  $p_i^T r_k=0, i=[0,k-1]$ via induction:
  \begin{proof}
    \begin{align*}
      &\text{base case of k=1, i=0}: x_1 = x_0 + \alpha_0 p_0\\
      &\text{ minimizes } \phi \text{ along } p_0\\
      &\implies r_1^T p_0 = (r_0 + \alpha_0  A p_0)^T p_0= 0 \\
      &\text{case: }r_{k-1}^T p_i = 0, i=[0,k-2]:\\
      &r_k = r_{k-1} + \alpha_{k-1} A p_{k-1}\\
      &(\forall i\in[0,k-2]) p_i^T r_k = p_i^T r_{k-1} + \alpha_{k-1} p_i^T A p_{k-1}\\
      &\text{(A-conjugacy by construction):}\\
      &(\forall i=[0,k-2]) p_i^T A p_{k-1}=0\\
      &(\forall i=[0,k-2]) p_i^T r_{k-1} = 0\ \text{(by induction hypothesis)}\\
      &p_i^T r_k = 0, i=[0,k-1]
    \end{align*}
  \end{proof}

  \vfill\null
  \columnbreak
  
  \subsection{Conjugate Gradient Method}
  Idea:
  \begin{itemize}
  \item uses only previous search direction to compute current search direction
  \item $p_k$ set to linear combination of $-r_k$ and $p_{k-1}$
  \item impose $p_k^T A p_{k-1}=0$
  \end{itemize}
  \begin{align*}
    p_k &= -r_k + \beta_k p_{k-1}\\
    p_{k-1}^T A p_k &= -p_{k-1}^T A r_k + \beta p_{k-1}^T A p_{k-1}\\
    0 &= -p_{k-1}^T A r_k + \beta p_{k-1}^T A p_{k-1}\\
    \beta &= \frac{p_{k-1}^T A r_k}{p_{k-1}^T A p_{k-1}}\\
    p_0& = -(Ax_0-b)=-r_0
  \end{align*}
  \begin{algorithm}[H]
    $x_0 = ..$\\
    $r_0 \leftarrow Ax_0-b$\\
    $p_0 = -r_0$\\
    \For{$k=[0,..n-1]$}{
      \uIf{$r_k==0$}{
        return $x_k$
      }
      \Else{
        $\alpha_k \leftarrow \frac{-r_k^Tp_k}{p_k^T A p_k} = \frac{r_k^T r_k}{p_k^T A p_k}$\\
        $x_{k+1} \leftarrow x_k + \alpha_k p_k$\\
        $r_{k+1} \leftarrow Ax_{k+1} - b = r_k + \alpha_k A p_k$\\
        $\beta_{k+1} \leftarrow \frac{p_k^T A r_{k+1}}{p_k^T A p_k}=\frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$\\
        $p_{k+1} \leftarrow -r_{k+1} + \beta_{k+1} p_k$
      }
    }
    \caption{Basic Conjugate Gradient Algorithm\label{Algo_CGBasic}}
  \end{algorithm}
  $p_k$ and $r_k$ is within krylov subspace:\\
  $K(r_0;k) =span\{r_0, Ar_0, ..A^k r_0\}$\\
  if $r_k \not= 0$:\\
  $r_k^T r_i = 0, i=[0,k-1]$\\
  $span\{r_0,..,r_k\} = span\{r_0,Ar_0,..,A^k r_0\}$\\
  $span\{p_0,..,p+k\}=span\{r_0,Ar_0,..,A^k r_0\}$\\
  $p_k^T A p_i=0, i=[0,k-1]$\\
  then, $\{x_k\} \rightarrow x^*$ in at most n steps.\\

  \vfill\null
  \columnbreak
  
  Simplification:
  \begin{align*}
    &p_{k+1} \leftarrow -r_{k+1} + \beta_{k+1} p_k\\
    &\alpha_k \leftarrow \frac{-r_k^Tp_k}{p_k^T A p_k}\\
    &\alpha_k \leftarrow \frac{-r_k^T(-r_{k} + \beta_{k} p_{k-1})}{p_k^T A p_k}\\
    &(\forall i=[0,k-1]) r_k^T p_i=0 \implies \beta_k r_k^T p_{k-1}=0\\
    &\alpha_k \leftarrow \frac{r_k^T r_k}{p_k^T A p_k} \ \text{(simplified)}
  \end{align*}
  \begin{align*}
    &r_{k+1}=r_k+\alpha_k A p_k\\
    &A p_k = \frac{r_{k+1}-r_k}{\alpha_k}\\
    &\beta = \frac{p_k^T A r_{k+1}}{p_k^T A p_k}\\
    &p_k^T A p_k = p_k^T \frac{r_{k+1}-r_k}{\alpha_k}= \frac{-p_k^T r_k}{\alpha} (\text{conjugacy})\\
    &p_k^T A p_k = -\frac{(- r_k + \beta_k p_{k-1})^T r_k}{\alpha}=\frac{r_k^T r_k}{\alpha} (\text{conjugacy})\\
    &p_k^T A r_{k+1} = r_{k+1}^T A p_k\\
    &p_k^T A r_{k+1} = r_{k+1}^T \frac{r_{k+1}-r_k}{\alpha_k}\\
    &r_k \in span\{p_k,p_{k-1}\}\ and\ r_{k+1}^T p_i = 0, i=[0,k] \implies\\
    &p_k^T A r_{k+1} = \frac{r_{k+1}^T r_{k+1}}{\alpha_k}\\
    &\beta_{k+1} \leftarrow \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k} \ \text{(simplified)}
  \end{align*}

  \subsection{Nonlinear Method}
  Minimize general convex function or nonlinear function. Variants: FR, PR.

  \subsubsection{FR (Fletcher Reaves)}
  Modify linear CG by:
  \begin{itemize}
  \item replace residual by gradient of nonlinear objective, $r_k \rightarrow \nabla f_k$
  \item replace $\alpha_k$ computation by a linear search to find approx. minimum along search direction of $f$
  \end{itemize}

  Equivalent to linear CG if objective is strongly convex quadratic.\\

  Linear search for $\alpha_k$ with strong Wolfe condition to ensure $p_k$'s are descent directions wrt. objective function.\\
  
  \subsubsection{PR}
  Replace $\beta_{k+1}$ computation in FR with:
  \begin{align*}
    \beta_{k+1}^{PR} & \leftarrow \frac{\nabla f_{k+1}^T (\nabla f_{k+1} - \nabla f_k)}{\nabla f_k^T \nabla f_k}\\
    \beta_{k+1}^+ & \leftarrow \max(\beta_{k+1}^{PR}, 0)
  \end{align*}


  CG Gradient Algo. Property [T.5.3]:
  \begin{align*}
    &(\forall i=[0,k-1])\ r_k^T r_i=0\\
    &span\{r_0,..,r_k\} = span\{r_0,Ar_0,..,A^kr_0\}\\
    &span\{p_0,..,p_k\} = span\{r_0,Ar_0,..,A^kr_0\}\\
    &(\forall i=[0,k-1])\ p_k^T A p_i = 0\\
    &\implies \{x_k\}\ \text{converges to } x^*\ \text{in at most n steps}    
  \end{align*}
  Proof by induction to show generated search direcitons are A-conjugate. Then apply Theorem \ref{thm:CD_Term} to conclude algo terminates within n steps.
  
  \vfill\null
  \pagebreak
  
  \section{Quasi Newton}

  \subsection{Concept}
  properties: $O(n^2)$, self correcting, slightly more iterations than Newton Method, linear convergence order and superlinear rate of convergence

  Derivation, using 2nd order model, with $B_k$ SPSD:
  \begin{align*}  
    m_k(p)=f_k + \nabla f_k^T p + \frac{1}{2} p^T B_k p
  \end{align*}
  taking gradient wrt. $p$ and solve, assuming minimum exists:
  \begin{align*}  
    0 & = \nabla f_k + B_k p\\
    p & = - B_k^{-1} \nabla f_k
  \end{align*}
  update equation:
  \begin{align*}
    x_{k+1} & = x_k + \alpha_k p_k\\
    x_{k+1} & = x_k - \alpha_k B_k^{-1} \nabla f_k
  \end{align*}
  pick $\alpha$ to satisfy Wolfe conditions

  updated model at next iterate $x_{k+1}$:
  \begin{align*}
    m_{k+1}(p)=f_{k+1} + \nabla f_{k+1}^T p + \frac{1}{2} p^T B_{k+1} p
  \end{align*}
  impose reasonable conditions:
  \begin{enumerate}
  \item gradient of $m_{k+1}$ matches gradient of objective function $f$ at iterate $x_k$
  \item gradient of $m_{k+1}$ matches gradient of objective function $f$ at iterate $x_{k+1}$
  \end{enumerate}
  \begin{align*}
    \nabla m_{k+1}(0) & = \nabla f_{k+1} \implies \text{ condition 2 satisfied}
  \end{align*}

  for condition 1 (gradient of $m_{k+1}$ match gradient of objective $f$ at iterate $x_k$):
  \begin{align*}
    \nabla m_{k+1}(-\alpha_k p_k) &= \nabla f_{k+1} - \alpha_k B_{k+1} p_k = \nabla f_k\\
    \nabla f_{k+1} - \nabla f_k &= \alpha_k B_{k+1} p_k\\
    y_k &= \nabla f_{k+1} - \nabla f_k\\
    s_k &= x_{k+1} - x_k = \alpha_k p_k\\
    y_k &= B_{k+1} s_k \text{ [secant equation]}
  \end{align*}

  $B_{k+1}$ SPD $\implies s_k^T B_{k+1} s_k = s_k^T y_k > 0$\\
  $f$ strongly convex $\implies s_k^T y_k > 0$ satisfied for any 2 points $x_k$ and $x_{k+1}$.\\
  $f$ nonconvex $\implies$ need to enforce $s_k^T y_k$ by Wolfe conditions, and use line search for step length $\alpha$.
  
  Using line search ensures curvature condition:
  \begin{proof}
    \begin{align*}
      &s_k=x_{k+1} - x_k = \alpha_k p_k\\
      &y_k = \nabla f_{k+1} - \nabla f_k\\
      &\text{Wolfe curvature condition:}\\
      &\nabla f_{k+1}^T p_k \geq c_2 \nabla f_k^T p_k, c_2 \in (0,1)\\
      &y_k^T s_k = (\nabla f_{k+1} - \nabla f_k)^T (\alpha_k p_k)\\
      &y_k^T s_k = \alpha_k (\nabla f_{k+1} - \nabla f_k)^T p^k\\
      &y_k^T s_k = \alpha_k (c_2 \nabla f_k - \nabla f_k)^T p^k\\
      &y_k^T s_k = \alpha_k (c_2-1) \nabla f_k^T p^k\\
      &\alpha > 0 \wedge c_2-1 < 0 \wedge \nabla f_k^T p_k < 0 (p_k \text{ is a descent dir})\\
      &\implies y_k^T s_k > 0
    \end{align*}
    Curvature condition holds when Wolfe line search is performed.
  \end{proof}

  Choosing approximate Hessian, $B$:
  
  Force a unique solution among infinite many due to extra degrees of freedom in the matrix, few (n) constraining conditions imposed by secand equation, few (n) constraining conditions of PD. One approach is solving a optimization problem to make row rank update to previous iterate:
  \begin{align*}
    \min_B & ||B-B_k||\\
    s.t.\ &B=B^T \\
    &Bs_k=y_k \text{ [secant equation]}\\
    &B \succ 0
  \end{align*}
  alternatively, constrain $B$'s inverse, $H$:
  \begin{align*}
    \min_H & ||H-H_k||\\
    s.t.\ &H=H^T\\
    &Hy_k=s_k
    &H_{k+1}y_k & = s_k \text{ [secant equation]}\\
  \end{align*}
  
  Different norms can be used. One choice: weighted Frobenius norm:
  \begin{align*}
    ||A||_W & := ||W^{1/2}AW^{1/2}||_F\\
    ||X||_F & := (\sum_i \sum_j (X_{ij})^2)^{1/2}
  \end{align*}

  \subsection{DFP update algorithm}
  let $\bar{G_k}$ be the average Hessian:
  \begin{align*}
    \bar{G_k} = \int_0^1 \nabla^2 f(x_k + \tau \alpha_k p_k) d \tau
  \end{align*}
  using Taylor's theorem:\\
  $y_k = \bar{G_k} s_k = \bar{G_k} \alpha_k p_k$\\
  let $W=\bar{G_k}^{-1}$\\
  solve optimization problem:
  \begin{align*}
    &\min_B \| B - B_k \|_W\\
    s.t.\ &W=\bar{G_k}^{-1}\\
    &Bs_k=y_k\\
    &B=B^T
  \end{align*}
  
  Solution:
  \begin{align*}
    B_{k+1} &= (I-\rho_k y_k s_k^T)B_k(I - \rho_k s_k y_k^T) + \rho_k y_k y_k^T\\
    \rho_k &= \frac{1}{y_k s_k}
  \end{align*}

  Computation simplification: optimize for inverse Hessian instead  which is used in update of search direction: $p_k = - B_k^{-1} \nabla f_k$\\
  Use Sherman-Morrison-Woodbury formula to obtain inverse.\\
  let $H_k = B_k^{-1}$, then DFP update becomes:\\
  \begin{align*}
    H_{k+1} = H_k - \frac{H_k y_k y_k^T H_k}{y_k^T H_k y_k} + \frac{s_k s_k^T}{y_k^T s_k}
  \end{align*}
  This is a rank 2 modification to previous iterate for efficiency.
  
  \subsection{BFGS update algorithm}

  Idea: impose condition on inverse of Hessian instead of Hessian approximation.\\

  let $H$ be approximate inverse of Hessian, then impose:
  \begin{align*}
    H_{k+1} \succ 0\\
    H_{k+1} = H_{k+1}^T\\
    H_{k+1} y_k = s_k
  \end{align*}

  solve for $H$ in the optimization problem:
  \begin{align*}
    &\min_H \|H-H_k\|_W\\
    s.t.\ &H=H^T\\
    &H y_k = s_k\\
  \end{align*}
  where $y_k=W s_k$\\
  
  let $W$ be the average Hessian:\\
  $\bar{G_k}=\int_0^1 \nabla^2 f(x_k + \tau \alpha_k p_k) d \tau$\\

  solve to obtain:
  \begin{align*}
    H_{k+1} & = (I-\rho_k s_k y_k^T) H_k (I - \rho_k y_k s_k^T) + \rho_k s_k s_k^T\\
    \rho_k & = \frac{1}{y_k^T s_k}
  \end{align*}
  
  %%   Weight matrix $W$ satisfies $Ws_k=y_k$.\\
  %%   $W$ is the average Hessian $\bar{G}$:
  %%   \begin{align*}
  %%     \bar{G} = \int_0^1 \nabla^2 f(x_k+\tau \alpha_k p_k) d\tau
  %%   \end{align*}

  %%   solution given by:
  %%   \begin{align*}
  %%     H_{k+1} & = (I-\rho_k s_k y_k^T)H_k(I-\rho_ky_ks_k^T) \\
  %%             &+ \rho_k s_k^T s_k\\
  %%     \rho_k & = \frac{1}{y_k^Ts_k}
  %%   \end{align*}

  Deduce from rank 2 update:\\
  $B_{k+1} = B_k + U_k + V_k$\\
  $U_k, V_k$ are rank 1 symmetric matrices:\\
  $U_k = \alpha uu^T, V_k = \beta vv^T$\\
  Let $u=y_k, v=B_ks_k$\\
  $B_{k+1} = B_k + \alpha y_k y_k^T + \beta B_k s_k (B_k s_k)^T$\\
  Impose secant condition:\\
  $B_{k+1} s_k = y_k = B_k s_k + \alpha y_k y_k^T s_k + \beta B_k s_k (B_k s_k)^T s_k$\\
  Solve for $\alpha, \beta$:\\
  $\alpha = \frac{1}{y_k^T s_k}$\\
  $\beta = \frac{-1}{(B_k s_k)^T s_k}$\\
  $B_{k+1} = B_k + \frac{y_k y_k^T}{y_k^T s_k} - \frac{B_k s_k (B_k s_k)^T}{s_k^T B_k^T s_k}$\\
  Use Sherman-Morrison formula to get $B^{-1}$ and optimize with inverse Hessian in order to avoid computing inverses in the actual algorithm:\\
  $B_{k+1}^{-1} = H_{k+1} = (I- \frac{s_k y_k^T}{y_k^T s_k})H_k(I- \frac{y_ks_k^T}{y_k^T s_k}) + \frac{s_k s_k^T}{y_k^T s_k}$\\
  
  initial $H_0$ can be chosen approximately:
  \begin{itemize}
  \item finite diferences
  \item $\alpha I$
  \end{itemize}
  
  \begin{algorithm}[H]
    \SetKwInOut{Input}{$H_0,x_0,\epsilon>0$}\SetKwInOut{Output}{$x$}
    \Input{inverse Hessian approx., initial point, convergence tolerance}
    \Output{solution}
    $k \leftarrow$ 0\\
    //until convergence\\
    \While{$||\nabla f_k|| > \epsilon$}{
      $p_k \leftarrow -B_k^{-1} \nabla f(x_k) = -H_k\nabla f(x_k)$\\
      $\alpha_k \leftarrow $ LineSearch($p_k, f, \nabla f$)\\
      $x_{k+1} \leftarrow x_k + \alpha_k p_k$\\
      $s_k \leftarrow x_{k+1} - x_k$\\
      $y_k \leftarrow \nabla f_{k+1} - \nabla f_k$\\
      $\rho_k \leftarrow \frac{1}{y_k^Ts_k}$\\
      //update inverse Hessian approximation for next iterate\\
      $H_{k+1} \leftarrow (I-\rho_k s_k y_k^T)H_k(I-\rho_k y_k s_k^T) + \rho_k s_k s_k^T$\\
      $k \leftarrow k+1$\\
    }
    return $x$
    \caption{BFGS Algorithm\label{Algo_BFGS}}
  \end{algorithm}

  Cost: $O(n^2)$ + cost of eval $f(.)$ + cost of eval $\nabla f(.)$.\\
  Order of convergence: superlinear, worse than Newton but more computationally efficient than Newton.\\   
  
  Using Sherman-Morrison-Woodbury formula to obtain Hessian update equation,
  \begin{align*}
    B_{k+1} & = B_k - \frac{B_ks_k s_k^T B_k^T}{s_k^T B_k^T s_k} + \frac{y_k y_k^T}{y_k^T s_k}
  \end{align*}
  , but is it more efficient to use the inverse version.\\
  
  Proper line search is required so that BFGS algo captures curvature information.\\
  
  Inaccurate line search can be used to reduce computation cost.

  \subsection{Sherman-Morrison}
  \begin{align*}
    & A^{-1}\ exists \implies\\
    & (A+uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u} \iff 1 + v^T A^{-1} u \neq 0
  \end{align*}

  \vfill\null
  \pagebreak
  
  \section{Trust Region Methods}
  idea:
  \begin{itemize}
  \item models local behaviour of the objective function (eg: 2nd order Taylor series)
  \item set local region to explore, then simultaneously find direction and step size to take
  \item region size adaptively set using results from previous iterations, using ratio of function value decrease vs. model value decrease
  \item step may fail due to inadequately set region, which need to be adjusted
  \item superlinear convergence when approximate model Hessian is equal to true Hessian
  \end{itemize}

  using 2nd order Taylor series model with symmetric matrix approximating Hessian
  \begin{align*}
    \min_{p \in \R^n} m_k(p) & = f_k + g_k^Tp + \frac{1}{2}p^TB_kp,\ st.\ ||p||\leq \Delta_k\\
    \Delta_k & :=\text{trust region radiius}\\
    g_k & =\nabla f(x_k)\\
    B_k & \succeq 0\\
    B_k^T & = B_k
  \end{align*}
  
  full step is ($p_k=-B_k^{-1}g_k)$ taken when $B\succ 0 $ and $||B_k^{-1} g_k|| \leq \Delta_k$

  evaluate goodness of model with $\frac{\text{actual function value change}}{\text{predicted model change}}$:
  \begin{align*}
    &\rho_k = \frac{f(x_k)-f(x_k+p_k)}{m_k(0)-m_k(p_k) \geq 0}\\
    &action \leftarrow
    \begin{cases}
      \text{expand trust region} &, \rho_k \approx 1\ (agreement)\\
      \text{shrink trust region} &, \rho_k  < 0 + \text{thresh}\\
      \text{keep trust region} &, o/w
    \end{cases}  
  \end{align*}

  \begin{algorithm}[H]
    $k \leftarrow$ 0\\
    \While{$||\nabla f_k|| > \epsilon$}{
      $p_k \leftarrow
      \argminl_{p} f_k + g_k^Tp + \frac{1}{2}p^TB_kp,\ st.\ ||p||\leq \Delta_k$\\
      $\rho_k = \frac{f(x_k)-f(x_k+p_k)}{m_k(0)-m_k(p_k)}$//ratio test\\
      //change trust region for next iterate\\
      \uIf{$\rho_k<\gamma (:\frac{1}{4})$}{
        $\Delta_{k+1} \leftarrow \alpha (:\frac{1}{4}) \Delta_k //shrink$
      }
      \uElseIf{$\rho_k> \beta (:\frac{3}{4})\ and\ ||p_k|| == \Delta_k$}{
        $\Delta_{k+1} \leftarrow min(2 \Delta_k, \hat{\Delta}) //expand$
      }
      \Else{
        $\Delta_{k+1} \leftarrow \Delta_k //keep$
      }

      //accept/reject for current iterate\\
      \uIf{$\rho_k > \eta (: \in [0, \frac{1}{4}))$}{
          $x_{k+1} \leftarrow x_k + p_k$
        }
        \Else{
          $x_{k+1} \leftarrow x_k$
        }
      }
      \caption{Trust Region Algorithm\label{Algo_TrustRegion}}
  \end{algorithm}

  Minimizer of the 2nd order Taylor series satisfy the following as a global solution for trust region iff:
  \begin{align*}
    (B+\lambda I)p^*=-g\\
    \text{complementary slackness:}\\
    \lambda(\Delta - ||p^*||)=0\\
    (B+\lambda I) \succeq 0\\
    \lambda \geq 0
  \end{align*}

  \begin{align*}
    \lambda \ge 0 \implies p \text{ parallel to negative gradient of model}
  \end{align*}
  
  Solving 2nd order Taylor series using approx methods, which reduce at least as much as Cauchy Point:
  \begin{itemize}
  \item dogleg (positive definite)
  \item 2-D subspace minimization (require smallest eigenvalue estimation if not positive definite)
  \item conjugate gradient based, effective when $B$ large, sparse
  \end{itemize}
  
  Cauchy Point (sufficient reduction in model value):\\
  Use 1st order approx. of model and steepest descent direciton to get next iterate, bounded within trust region.\\

  \vfill\null
  \pagebreak
  
  \subsection{Cauchy Point}
  Idea: find an approximate step that would still achieve convergence in a sequence of steps.\\
  Use linear version of the model:
  \begin{align*}
    & p_k^s = \argmin_p f_k + g_k^T p\\
    & s.t.\ \|p\| \leq \Delta_k
  \end{align*}

  $p_k^s = - \frac{\|\Delta_k\|}{\|g_k\|} g_k$\\

  Solve for step length in the original model with trust region:
  \begin{align*}
    & \tau_k = \argmin_{\tau \geq 0 } m_k(\tau p_k^s)\\
    & s.t.\ \|\tau p_k^s\| \leq \Delta_k
  \end{align*}

  Cauchy point: $p_k^c = \tau_k p_k^s$\\

  Solutions for $\tau_k$:
  \begin{align*}
    \begin{cases}
      g_k^T B_k g_k \leq 0 &, g_k \neq 0 \implies m_k(\tau_k p_k^s)\\
      & monotonically\ decreasing,\\
      & \tau_k = 1,\\
      g_k^T B_k g_k > 0 &, m_k(\tau_k p_k^s)\ convex\ \\
      & quadratic,\\
      & \tau_k=min(1, \frac{\|g_k\|^3}{\Delta_k g_k^T B_k g_k}),\\
      & either: boundary\ val\ of\ \Delta_k,\\
      & or\ unconstrained\ minimizer\ of\ \\
      & quadratic
    \end{cases}
  \end{align*}

  Cauchy point $p_k^c$ inexpensive to calculate.\\

  [TODO] derive $\tau_k$ for case of $g_k^T B_k g_k > 0$.
  
  \subsection{Dogleg method}

  Idea: improve on Cauchy point in considering curvature information. Requires $B \succ 0$.\\
  
  if $B \succ 0$:\\
  $p^B=-B^{-1}g$\\
  $p^*=p^B$ if $\Delta \geq \|p^B\|$\\
  $p^U=\frac{-g^Tg}{g^TBg}g$ (intermediate point along direction of steepest descent)\\
  Resulting algo, interpolate between $p^U$ and $p^B$:\\
  $\tilde{p}(\tau) =
  \begin{cases}
    \tau p^U &, \tau_\in [0,1]\\
    p^U + (\tau-1)(p^B-p^U) &, \tau \in [1,2]
  \end{cases}
  $
  \begin{align*}
    B \succ 0 \implies & \|\tilde{p}(\tau)\| \text{increases wrt. }\tau\ \land\\
    & m(\tilde{p}(\tau)) \text{decreases wrt. } \tau
  \end{align*}
  
  if $\|p^B\| \leq \Delta$: $p$ chosen at $p^B$\\
  else $p$ chosen at intersection of $\tilde{p}(\tau)$ and trust region boudnary by solving:\\
  $\|p^U+(\tau-1)(p^B-p^U)\|^2 = \|\Delta^2\|$

  $p_k^S=\argminl_p f_k+g_k^Tp, \|p\| \leq \Delta_k$\\
  $\tau_k = \argminl_{\tau\geq0} m_k(\tau p_k^S), \|\tau p_k^S\| \leq \Delta_k$\\
  $p_k^S=\frac{-\Delta_k g_k|}{\|g_k\|}$\\
  $p_k^C=\tau_k p_k^S$\\
  $p_k^C=-\tau_k \frac{\Delta_k g_k}{\|g_k\|}$\\
  $\tau_k=\begin{cases}
  1 &, g_k^T B_k g_k \leq 0\\
  min(\frac{\|g_k|^3}{\Delta_k g_k^T B_k g_k}, 1) &, o/w
  \end{cases}$

  [TODO: add section on 2D subspace minimization method]\\
  
  \vfill\null
  % \pagebreak

  \subsection{Iterative Solution}
  Idea: solve subproblem $\min_{\|p\| \leq \Delta} m(p)$ by applying Newton's method to find $\lambda$ that matches trust region radius. This is slightly more accurate per step compared to Dogleg. Use $(B+\lambda I)p^*=-g$ to solve $\min_{\|p\| \leq \Delta} m(p)$ for $\lambda$.

  If $\lambda=0$ and $(B+\lambda I)p^*=-g, \|p^*\| \leq \Delta$ and $(B+\lambda I) \succeq 0$: return $\lambda$\\
  Else: find $\lambda$ s.t. $(B+\lambda I) \succeq 0$ and $\|p(\lambda)\| = \Delta, p(\lambda)=-(B+\lambda I)^{-1}g$. Solve and return $\lambda$.
  
  Solve $\|p(\lambda)\|-\Delta=0, \lambda > \lambda_1$ via Newton's method (root finding). Approx. this to nearly a linear problem for easy solving:\\

  \begin{algorithm}[H]
    \For{$l=0,1,..$}{
      solve $B+\lambda ^l I = R^T R$\\
      $R^TR p_l = -g$\\
      $R^T q_l = p_l$\\
      $\lambda^{l+1} \leftarrow \lambda^l + (\frac{\|p_l\|}{\|q_l\|})^2(\frac{\|p_l\|-\Delta)}{\Delta})$
      check $\lambda \geq \lambda_1$
    }
    \caption{Subproblem Algo\label{Algo_TrustRegionIterativeSubproblem}}
  \end{algorithm}

  \vfill\null
  \columnbreak
  
  \subsection{Trust Region Newton CG Method}

  \label{section_trust_region_newton_cg}
  Idea: use trust region algo for the outer iteration, use iterative CG based algorithm to solve for inner optimization problem.

  Outer iteraion: Algo \ref{Algo_TrustRegion}.\\

  Inner optimization problem:
  \begin{align*}
    \min_{p \in \R^n} m(p) & = f + g^T p + \frac{1}{2} p^T B p\\
    & s.t.\ \|p\| \leq \Delta
  \end{align*}
  
  \begin{algorithm}[H]
    $\epsilon_k = \eta_k \| \nabla f_k\|$\\
    $z_0 = 0$\\
    $r_0 = \nabla f_k$\\
    $d_0 = -r_0 = \nabla f_k$\\
    \uIf{$\|r_0\| < \epsilon_k$}{
      return $p_k=z_0=0$
    }
    \For{$j=0,1,..$}{
      //dir of nonpositive curvature test\\
      \uIf{$d_j^T B_k d_j \leq 0$}{ 
        return $\argmin_{p_k} m_k(p_k=z_j + \tau d_j)$ s.t. $\|p_k\| = \Delta_k$
      }
      $\alpha_j = \frac{r_j^T r_j}{d_j^T B_k d_j}$\\
      $z_{j+1}=z_j + \alpha_j d_j$\\
      //trust region check\\
      \uIf{\|z_{j+1}\| \geq \Delta_k}{
        return $p_k : p_k = z_j + \tau d_j$ s.t. $\|p_k\| = \Delta_k, \tau \geq 0$
      }
      $r_{j+1} = r_j + \alpha_j B_k d_j$\\
      \uIf{$\|r_{j+1}\| < \epsilon_k=\eta_k \|\nabla f_k\|, \eta_k=min(\frac{1}{2},\|\nabla f_k\|^{\frac{1}{2}})$}{
        return $p_k=z_{j+1}$
      }
      $B_{j+1} = \frac{r_{j+1}^T r_{j+1}}{r_j^T r_j}$\\
      $d_{j+1} = -r_{j+1} + B_{j+1} d_j$\\
    }\\
    \caption{Trust Region Newton-CG Subproblem (CG Steinhaug) \label{Algo_CG_Steihaug}}
  \end{algorithm}
  
  $\epsilon_j$ can be chosen similarly as in Line Search Newton-CG method, where $\{\eta_k\}$ is the forcing sequence.
  
  \vfill\null
  \pagebreak
  
  \section{Proximal Algorithm}
  Idea:
  \begin{itemize}
  \item reliance on easy to evaluate proximal operators
  \item separability allows parallel evaluation
  \item generalization of projection based algorithms
  \end{itemize}

  \begin{align*}
    prox_{\lambda f}(v) = \argminl_{x} f(x) + \frac{1}{2 \lambda} \|x-v\|_2^2
  \end{align*}

  Resolvent of subdifferential operator:
  \begin{align*}
    &z=prox_{\lambda f}(x) \implies z \in (I+\lambda \partial f)^{-1}(x)\\
    &(I+\lambda \partial f)^{-1} :=\text{resolvent of operator }\partial f
  \end{align*}
  
  \subsection{Proximal Gradient Method}
  Solve $\min_{x} g(x)+f(x)$, where $f,g$ are closed, convex functions and $f$ differentiable.\\
  
  Idea: combine Proximal point algorithm with gradient projection method. Gradient step on $f$, and proximal point algo step on $g$.
  \begin{align*}
    &x^* = prox_{\lambda^k g}(x^k-\lambda^k \nabla f(x^k)) \\
    &=\argminl_{x} g(x) + \frac{1}{2 \lambda^k} \|x - (x^k - \lambda^k \nabla f(x^k))\|_2^2\\
    &\text{tradeoff between g and and gradient step}\\
    &g=I_C(x) \implies \text{projected gradient step}\\
    &g=0 \implies \text{gradient descent}\\
    &f=0 \implies \text{proximal minimization}
  \end{align*}

  Relation to Pixed Point:\\
  $x^*$ is a fixed point solution of $\min_{x} g(x)+f(x)$ iff $0 \in \nabla f(x^*) + \partial g(x^*)$ iff $x^* = (I+\lambda \partial g)^{-1}(I-\lambda \nabla f)(x^*)$\\
  
  Forward Euler, Backward Euler stepping is same as the proximal gradient iteration, $prox_{\lambda^k g}(x^k-\lambda^k \nabla f(x^k))$\\

  \subsection{Accelerated Proximal Gradient Method}
  Introduce extrapolation:
  \begin{align*}
    &y^{k+1} = x^k + w^k(x^k-x^{k-1})\\
    &x^{k+1} = prox_{\lambda^k g}(y^{k+1}-\lambda^k \nabla f(y^{k+1}))\\
    &w^k \in [0,1)
  \end{align*}

  Example: $w^k = \frac{k}{k+3}, w^0 = 0, \lambda^k \in (0,1/L], L:= \text{Lipschitz constant of }\nabla f,\ or\ \lambda^k$ found via line search.
    Line search for $\lambda^k$ (Beck and Teboulle):

    \begin{algorithm}[H]
      $\hat{f}(x,y):=f(y) + \nabla f(y)^T(x-y) + \frac{1}{2\lambda}\|x-y\|^2_2$\\
      \While{True}{
        $z = prox_{\lambda g}(y^k-\lambda \nabla f(y^k))$\\
        \uIf{$f(x) \leq \hat{f}(z,y^k)$}{
          break
        }
        $\lambda = \beta \lambda$
      }
      return $\lambda^k := \lambda, x^{k+1} := z$\\
      \caption{Proximal Gradient Algorithm\label{Algo_ProximalGradient}}
    \end{algorithm}

    \subsection{Types of Proximal Operators}
    \begin{itemize}
    \item quadratic functions
      \begin{align*}
        f&=\frac{1}{2} \|.\|_x^2 \implies prox_{\lambda f}(v) = (\frac{1}{1+\lambda})v\\
        f&=\frac{1}{2}x^T A x + b^T x + c, A \in S_+^n \implies\\
        &prox_{\lambda f}(v)=(I+\lambda)^{-1}(v-\lambda b)
      \end{align*}
    \item unconstrained problem: use gradient methods such as Newton, Quasi-Newton
    \item constrained: use projected subgradient for non-smooth, projected gradient or interior method for smooth
    \item separable function: if scalar, may be solved analytically, eg:
      L1 norm separable to:
      \begin{align*}
        f(x) = |x| \implies prox_{\lambda f}(v) =
        \begin{cases}
          v-\lambda ,& v \geq \lambda\\
          0,& |v| \leq \lambda\\
          v+\lambda, & v \leq -\lambda
        \end{cases}
      \end{align*}
      \begin{align*}
        f(x) = -log(x) \implies prox_{\lambda f}(v) = \frac{v+\sqrt{v^2+4\lambda}}{2}
      \end{align*}
    \item general scalar function
      \begin{itemize}
      \item localization: using a subgradient oracle and bisection algorithm
      \item twice continuously differentiable: guarded Newton method
      \end{itemize}
    \item polyhedra constraint, quadratic objective: solve as QP problem
      \begin{itemize}
      \item duality to reduce number of variables to solve if possible
      \item gram matrix caching
      \item affine constraint($Ax=b$): use pseudo-inverse, $A^{+}$:
        \begin{align*}
          \Pi_C(v) &= v - A^{+}(Av-b)\\
          &A \in R^{m \times n}, m < n \implies A^{+} = A^T(AA^T)^{-1}\\
          &A \in R^{m \times n}, m > n \implies A^{+} = (A^TA)^{-1}A^T
        \end{align*}
      \item hyperplane constaint($a^Tx=b$):
        \begin{align*}
          \Pi_C(v) = v + (\frac{b-a^Tb}{\|a\|_2^2})a
        \end{align*}
      \item halfspace
        \begin{align*}
          \Pi_C(v) = v - \frac{max(a^T v -b,0)}{\|a\|_2^2}a
        \end{align*}
      \item box($l\leq x \leq u$)
        \begin{align*}
          \Pi_C(v)_k = min(max(v_k, l_k), u_k)
        \end{align*}
      \item probability simplex($1^Tx= 1, x \geq 0$)\\
        bisection algo on $\nu$:
        \begin{align*}
          \Pi_C(v) = (v- \nu 1)_+\\
          \text{intial }[l_k,u_k] = [\max_i v_i-1, \max_i vi]
        \end{align*}
        analytically solve when bounded in between 2 adjacent v_i's      
      \end{itemize}
    \item cones ($\kappa$: proper cone)\\
      problem of the form:
      \begin{align*}
        min_x \|x-v\|_2^2\\
        s.t.: x\in \kappa
      \end{align*}
      \begin{align*}
        x \in \kappa\\
        v=x-\lambda\\
        \lambda\in\kappa^*\\
        \lambda^Tx=0
      \end{align*}
      \begin{itemize}
      \item cone $C=\R_+^n$
        \begin{align*}
          \Pi_C(v)=v_+
        \end{align*}
      \item 2nd order cone $C=\{(x,t)\in \R^{n+1}: \|x\|_2 \leq t\}$
        \begin{align*}
          \Pi_C(v,s) = \begin{cases}
            0,& \|v\|_2 \leq -s\\
            (v,s),& \|v\|_2 \leq s\\
            \frac{1}{2}(1+\frac{s}{\|v\|_2})(v,\|v\|_2),& \|v\|_2 \geq |s|
          \end{cases}
        \end{align*}
      \item PSD cone $S_+^n$
        \begin{align*}
          \Pi_C(V)&=\sum_i (\lambda_i)_+ u_i u_i^T\\
          V&=\sum_i \lambda_i u_i u_i^T\ (eigendecomp)
        \end{align*}
      \item exponential cone\\
        Todo
      \end{itemize}
    \item pointwise supremum
      \begin{itemize}
      \item max function
      \item support function
      \end{itemize}
    \item norms
      \begin{itemize}
      \item L2
      \item L1
      \item L-$\inf$
      \item elastic net
      \item sum of norms
      \item matrix norm
      \end{itemize}
    \item sublevel set
    \item epigraph
    \item matrix functions
      Todo
    \end{itemize}

    \vfill\null
    \pagebreak
    
    \section{Smoothness}
    $\beta$ smoothness:
    \begin{align*}
      & \frac{\beta}{2} \|x\|_2^2 - f(x) \text{ is convex (fit quadratic on top)}\\
      & f(y) \leq f(x) + \nabla f(x)^T (y-x) + \frac{\beta}{2} \|y-x\|_2^2\\
      & \implies \text{f is smooth}
    \end{align*}

    $\alpha$ strong convexity:
    \begin{align*}
      & f(x) - \frac{\alpha}{2} \|x\|_2^2 \text{ is convex (fit quadratic below)}\\
      & f(y) \geq f(x) + \nabla f(x)^T (y-x) + \frac{\alpha}{2} \|y-x\|_2^2\\
      & \text{f may not be smooth}
    \end{align*}

    Smoothness near optimal point: there always exists a step (length $<\frac{2}{\beta}$) that is a descent direction.

    Property of f wrt. optimality value:\\
    f is $\beta$-smooth $\implies$
    \begin{align*}
      \frac{1}{2 \beta} \|\nabla f(x) \|_2^2 \leq f(x) -f(x^*) \leq \frac{\beta}{2} \| x-x^* \|_2^2
    \end{align*}
    f is $\alpha$-strongly convex $\implies$
    \begin{align*}
      \frac{\alpha}{2} \| x-x^* \|_2^2 \leq f(x) -f(x^*) \leq \frac{1}{2 \alpha} \|\nabla f(x) \|_2^2
    \end{align*}

    Co-coercivity for $\beta$-smooth $f$:
    \begin{align*}
      \left(\nabla f(x) - \nabla f(y)\right)^T(x-y) \geq \frac{1}{\beta} \|\nabla f(x) - \nabla f(y)\|_2^2
    \end{align*}
    Coercivity for $\alpha$-strongly convex $f$:
    \begin{align*}
      \left(\nabla f(x) - \nabla f(y)\right)^T(x-y) \geq \alpha \|x - y\|_2^2
    \end{align*}
    
    \subsection {Oracle based lower bounds}
    Lipschitz convex function: error=$\mathcal{O}(\frac{1}{\sqrt{T}})$\\
    Smooth convex function: $\epsilon=\mathcal{O}(\frac{1}{T^2})$\\
    Smooth + strongly convex function: $\epsilon=\mathcal{O}(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1})$

    \subsection {Gradient descent rate of convergence}
    $\beta$-smooth $f$: $\mathcal{O}(\frac{1}{T})$\\
    $\beta$-smooth + $\alpha$-strongly convex $f$: $\mathcal{O}((\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1})^T)$

    \subsection {Gradient descent with momentum}
    $f$ $\beta$-smooth, $\alpha$-s.c.:
    \begin{align*}
      & x_1=y_1 = x_{init}\\
      & y_{t+1} = x_t - \frac{1}{\beta} \nabla f(x_t)\\
      & x_{t+1} = \left(1+\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right) y_{t+1} - \left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right) y_t
    \end{align*}

    $f$ $\beta$-smooth:
    \begin{align*}
      & x_{t+1}=x_t-\eta \nabla f(x_t) \text{ [normal GD]}\\
      & \text{momentum, nesterov accel.:}\\
      & d_{t+1} = \gamma_{t+1}(x_{t+1} - x_t)\\
      & x_{t+1} = (x_t + d_t) - \eta \nabla f(x_t + d_t)\\
      & \epsilon=\mathcal{O}\left(\frac{1}{T^2}\right)
    \end{align*}

    \vfill\null
    \pagebreak
    
    \section {Common Algorithms}

    \subsection {Projected Gradient}
    \begin{align*}
      & y_{t+1} = x_t - \eta g_t, g_t \in \partial f(x_t)\\
      & x_{t+1} = proj_{X}(y_{t+1})\\
      & x^*=proj_{X}(x^*) \text{ [fixed point at optimality]}
    \end{align*}
    
    \subsection {Proximal Gradient}
    \begin{align*}
      prox_{\eta h}(y) = \argmin_x h(x) + \frac{1}{2 \eta} \|x-y\|_2^2\\
      x_{t+1} = prox_{\eta h}(x_t - \eta \nabla f(x_t))
    \end{align*}
    
    Special case of projected gradient:\\
    let $h(x)=I_X(x) \implies prox_{\eta h} = proj_X$\\
    where $I_X(x) =
    \begin{cases}
      0 &, x \in X\\
      +\infty &, o/w
    \end{cases}$

    \subsubsection{Common proximal operators}
    $h(x) = \|x\|_1$:
    \begin{align*}
      (prox_{\eta h}(x))_i & = \begin{cases}
        x_i - \eta &, x_i \geq \eta\\
        0 &, |x_i| \leq \eta\\
        x_i + \eta &, x_i \leq -\eta
      \end{cases}\\
      & = \max(|x_i|-\eta,0)\ sign(x_i)
    \end{align*}
    
    $h(x)=\frac{1}{2} x^T Q x + q^T x + q_0, Q \succeq 0$:
    \begin{align*}
      prox_{\eta h} = (I + \eta Q)^{-1} (x- \eta q)
    \end{align*}

    $h(x)=\sum_i h_i(x_i)$:
    \begin{align*}
      (prox_{\eta h}(x))_i = prox_{\eta h_i}(x_i) \text{ [parallelism]}
    \end{align*}

    Composite function:\\
    $\min_x f(x) = g(x) + h(x)$, $g$ smooth, $h$ not smooth but has prox. operator. Convergence due to non-smooth function: $\mathcal{O}(\frac{1}{\sqrt{T}})$

    \subsection{L1 regularization with subgradient descent}
    
    \begin{align*}
      & \min_x \|Ax-y\|_2^2 + \labda \|x\|_1\\
      & g=\|Ax-y\|_2^2\\
      & h=\lambda \|x\|_1\\
      & \partial_x h(x) = \begin{cases}
        -1 &, x_i < 0\\
        [-1, 1] &, x_i = 0\\
        1 &, x_i > 0
      \end{cases}\\
      & x_{t+1} = x_t - \eta (\nabla g(x_t) + \lambda \partial h(x))\\
      & x_{t+1} = x_t - \eta (2A^T(Ax-y) + \lambda z)\\
      & z=\begin{cases}
      -1 &, x_i < 0\\
      [-1,1] &, x_i = 0\\
      1 &, x_i > 0
      \end{cases}
    \end{align*}
    Convergence: $\epsilon=\mathcal{O}\left(\frac{1}{\sqrt{T}}\right)$
    
    \subsection{ISTA (with proximal gardient)}
    \begin{align*}
      & \min_x g(x) + h(x), g\ \beta\text{-smooth}\\
      & g(x) = \|Ax-y\|_2^2\\
      & h(x) \lambda \|x\|_1\\
      & x_{t+1} = prox_{\eta h}(x_t - \eta \nabla g(x_t))\\
      & g\ \beta \text{-smooth, select } \eta = \frac{1}{\beta}\\
      & x_{t+1} = prox_{\frac{1}{\beta} \lambda \|\dot\|_1}(x_t - \frac{1}{\beta} \nabla g(x_t))\\
      & x_{t+1} = \argmin_x \lambda \|x\|_1 + \frac{\beta}{2} \|x-(x_t-\frac{1}{\beta} \nabla g(x_t))\|_2^2\\
      & \nabla g(x_t) = 2A^T(Ax-y)
    \end{align*}
    Convergence: $\mathcal{O}\left(\frac{1}{\sqrt{T}}\right)

    \subsection{FISTA}
    \begin{align*}
      & \min_x g(x) + h(x), g\ \beta\text{-smooth}\\
      & \mu_0=0\\
      & \mu_t = \frac{1+\sqrt{1+4\mu_t^2-1}}{2}\\
      & \gamma_t = \frac{1-\mu_t}{\mu_{t+1}}\\
      & x_1=z_1=\text{arbitrary}\\
      & x_{t+1} = (1-\gamma_t) x_{t+1} + \gamma_t z_t\\
      & z_{t+1} = prox_{\frac{\lambda}{\beta}\| \dot \|_1}(x_t - \frac{1}{\beta} \nabla g(x_t))
    \end{align*}
    Convergence: $\epsilon=\mathcal{O}\left(\frac{1}{T^2}\right)

    \vfill\null
    \pagebreak
    
    \section{Subgradient Method}
    TODO

    \vfill\null
    \pagebreak
    
    \section{Large Scale Optimization}
    \begin{enumerate}
    \item memory reduction techniques with quasi-Newton methods.
    \item Inexact Newton methods with superlinear local convergence properties, indefinite Hessian, Hessian-free impl.
    \item function that is partial separable to use smaller subspaces.
    \end{enumerate}

    
    \subsection{Inexact Newton Methods}
    
    Iterative method with local convergence properties.\\
    Residual: $r_k = \nabla^2 f_k p_k + \nabla f_k$\\
    Termination condition for iterative solver: $\|r_k\| \leq \eta_k \|\nabla f_k\|, \eta_k \in (0,1)$, where ${\eta_k}$ is the forcing sequence.\\
    Convergence properties: see Thm 7.1, 7.2 (Nocedal)\\
    
    2 such algorithms, using CG as the iterative solver:
    \begin{enumerate}
    \item Line Search Newton-CG
    \item Trust Region Newton-CG (Section \ref{section_trust_region_newton_cg})
    \end{enumerate}

    \subsubsection{Line Search Newton-CG}
    Augment original CG algorithm with line search, curvature check and forcing function check.\\
    
    \begin{algorithm}[H]
      $x_0 = ..$\\
      \For{$k=0,1,..$}{
        \eta_k = \min(\frac{1}{2}, \| \nabla f_k\|^{\frac{1}{2}})\\
        $z_0 = 0, r_0=\nabla f_k, d_0=-r_0=-\nabla f_k$\\
        \For{$j=0,1,..$}{
          //choose direction\\
          \uIf{$d_j^T B_k d_j \leq 0$}{
            \uIf{$j=0$}{
              break $p_k = - \nabla f_k$
            }
            \Else{
              break $p_k = z_j$
            }
          }
          $\alpha_j = \frac{r_j^T r_j}{d_j^T B_k d_j}$\\
          $z_{j+1}  = z_j + \alpha_j d_j$\\
          $r_{j+1} = r_j + \alpha_j B_k d_j$\\\
          //forcing function constraint\\
          \uIf{$\|r_{j+1}\| < \epsilon_k = \eta_k \|\nabla f_k\|$}{
            break $p_k = z_{j+1}$
          }
          $B_{j+1} = \frac{r_{j+1}^T r_{j+1}}{r_j^T r_j}$\\
          $d_{j+1} = -r_{j+1} + B_{j+1} d_j$\\
        }
        $\alpha_k = $ Wolfe/Goldstein/Armijo step search\\
        $x_{k+1} = x_k + \alpha_k p_k$\\
      }\\
      \caption{Line Search Newton-CG Algo \label{Algo_line_search_newton_cg}}
    \end{algorithm}

    \subsection{Partially Separable Functions}
    Decompose into sum of element functions with a large number of variables being linearly independent.\\

    Ideally: $\nabla f(x) = \sum_i \nabla f_i(x), \nabla^2 f(x) = \sum_i \nabla^2 f_i(x)$.\\

    Use Hessian approximations for individual element functions.\\ 

    Use chain rule to transform into msaller spaces and back to original space via compactifying matrices.\\

    eg: $\nabla^2 f(_1(x) \approx U_i^T B_{[1]} U_i$ mapping $B_{[1]}$ back into correct position in the full Hessian approximation.

    \begin{align*}
      x_{[i]} &= U_i x\\
      f(x) &= \sum_i \phi_i(U_i x)\\
      \nabla f_i(x) &= U_i^T \nabla \phi_i(U_i x)\\
      \nabla^2 f_i(x) &= U_i^T \nabla^2 \phi_i(U_i x) U_i\\
    \end{align*}

    Use a Quasi-Newton to approximate $\nabla^2 \phi_i$ which require smaller resources by updating a new approximation $B_{[i]}$ iteratively:
    \begin{align*}
      s_{[i]} & = x_{[i]}^+ - x_{[i]}\\
      y_{[i]} & = \nabla \phi_i(x_{[i]}^+) - \nabla \phi_i(x_{[i]})\\
      B_{[i]} & \approx \nabla^2 \phi_i(U_i x) = \nabla^2 \phi_i(x_{[i]})
    \end{align*}

    Map back to place in original matrix:\\

    $\nabla^2 f_i(x) \approx U_i^T B_{[i]} U_i$\\

    Complete Hessian approximation then becomes:\\

    $B = \sum_i U_i^T B_{[i]} U_i$\\

    which can be used in trust region algo to get $p_k$:\\
    
    $B_k p_k = -\nabla f_k$.

    \vfill\null
    \pagebreak

    \section{Quadratic Programming}

    TODO

    \vfill\null
    \pagebreak

    \section{Sequential QP}

    TODO

    \vfill\null
    \pagebreak
    
    \section{Nonlinear Penalty Based Methods}
    
    \subsection{Quadratic Penalty}
    
    Inexact, iterative method. Idea: add a weighted quadratic penalty term per constraint to the objective. Form a sequence of probems with increasing weights to penalty terms.\\
    
    $\min_x f(x)\ s.t.\ c_i(x)=0, i\in \mathcal{E}$\\
    $Q(x;\mu) = f(x) + \frac{\mu}{2} \sum_i c_i(x)^2$\\
    $\mu_k \rightarrow +\infty\ as\ k \rightarrow +\infty$\\

    Use unconstrained optimization techniques and use previous iterate as starting point for current iteration.\\
    
    May diverges when penalty weight if not big enough.\\

    Ill conditioning of Hessian when $\mu_k$ becomes large at the minimizer. CG and Quasi-Newton not suitable. Newton's method may still have problems: numerical inaccuracies when solving linear equations to calculate step, and Taylor series 2nd order approximation onyl accurate in a small neighbourhood of x so it may not generate useful step unless $\mu_{k+1}$ is chosen slightly larger than $\mu_k$.


    let $A(x)^T = [\nabla c_i(x)]_{i \in \mathcal{E}}$\\
    $Q(x; \mu_k) = f(x) + \sum_{i \in \mathcal{E}} \mu_k c_i(x_k)$\\
    $\nabla_x Q(x; \mu_k) = \nabla f(x) + \sum_{i \in \mathcal{E}} \mu_k c_i(x_k) \nabla c_i(x_k)$\\
    $\nabla_{xx}^2 Q(x; \mu_k) = \nabla^2 f(x) + \sum_{i \in \mathcal{E}} \mu_k c_i(x) \nabla^2 c_i(x) + \mu_k A(x)^TA(x)$\\
    $\nabla_{xx}^2 Q(x; \mu_k) \approx \nabla_{xx}^2 L(x; \lambda^*) + \mu_k A(x)^T A(x)$ near minimizer where $\mu_k A(x)^T A(x)$ may be ill conditioned as $\mu_k \rightarrow \infty$\\
    
    For the case of equality and inequality constrained problem:
    \begin{align*}
      & \min_x f(x)\\
      & s.t.\ c_i(x)=0, i\in \mathcal{E}\\
      & s.t.\ c_i(x)\geq 0, i\in \mathcal{I}\\
      & Q(x;\mu) = f(x) + \frac{\mu}{2} \sum_i c_i(x)^2 + \frac{\mu}{2} \sum_i max(-c_i(x), 0)^2
    \end{align*}

    \vfill\null
    \columnbreak
    
    \begin{algorithm}[H]
      $x_0^s = ..$\\
      $\mu_0 > 0$\\
      $\{\tau_k\} \rightarrow 0$\\
      \For{$k=0,1,..$}{
        solve $x_k = \argmin_x Q(x; \mu_k)$ with initial guess $x_k^s$ with termination condition $\|\nabla_x Q(x;\mu_k)\| \leq \tau_k$\\
        \uIf{final convergence criteria passes}{
          return approx solution
        }
        update $\mu_{k+1}: \mu_{k+1} > \mu_k$ (eg: chosen adaptively)\\
        update $x_{k+1}^s \leftarrow x_k$
      }\\
      \caption{Quadratic Penalty Algo \label{Algo_quadratic_penalty}}
    \end{algorithm}

    \subsection{Nonsmooth Penalty Methods}

    \subsubsection{L1 Penalty}
    Exact method, not rely on updating stratefor penalty parameter. Penalty function:\\

    $\phi(x;\mu) = f(x) + \mu \sum_{i \in \mathcal{E}}\|c_i(x)\| + \mu \sum_{i \in \mathcal{I}} \max(0, -c_i(x))$\\

    If $x^*$ is a strict local minimizer then $x^*$ is a local minimizer of $\phi(x;\mu)$ for $\forall \mu > \mu^*$, where $\mu^*=||\lambda^*||_{\infty} = \max_{i \in \mathcal{I} \cup \mathcal{E}} | \lambda_i^*|$, $\lambda_i$s being lagrange multiplers. This generates enough penalization when going into infeasible resion to keep $x^*$ as local minimizer.

    Nonsmooth function, but has directional derivatives $D(\phi(x;\mu);p)$ along $p$.\\

    $\hat{x}$ is a stationary point for $\phi(x;\mu)$ if $D(\phi(\hat{x}, \mu); p) \geq 0, \forall p$\\

    $\hat{x}$ is a stationary point of measure of infeasibility $h(x) = \sum_{i \in \mathcal{E}}|c_i(x)| + \sum_{i \in \mathcal{I}} max(0, -c_i(x))$ if $D(\phi(\hat{x}, \mu); p) \geq 0, \forall p$\\

    \begin{algorithm}[H]
      $x_0^s = ..$\\
      $\mu_0 > 0$\\
      $\tau = ..$\\
      \For{$k=0,1,..$}{
        solve approximately $x_k = \argmin_x \phi(x; \mu_k)$ with initial guess $x_k^s$
        \uIf{$h(x_k) \leq \tau$}{
          return $x_k$
        }
        update $\mu_{k+1}: \mu_{k+1} > \mu_k$\\
        update new starting point: $x_{k+1}^s$
      }\\
      \caption{Classic L1 Penalty Algo \label{Algo_l1_classic_penalty}}
    \end{algorithm}

    Practical L1 penalty method: linearize constraints, replace nonlinear objective by a quadratic. Convert to a smooth QP problem by new slack variables and solve with standard QP solver.
    \begin{align*}
      & \min_{p,r,s,t} f(x) + \frac{1}{2} p^T W p + \nabla f(x)^T p\\
      & + \mu \sum_{i \in \mathcal{E}} (r_i + s_i) + \mu \sum_{i \in \mathcal{I}} t_i\\
      s.t. & \nabla c_i(x)^T p + c_i(x) = r_i - s_i, i \in \mathcal{E}\\
      & \nabla c_i(x)^T p + c_i(x) \geq -t_i, i \in \mathcal{I}\\
      & r,s,t \geq 0
    \end{align*}    

    Possible issue: $\mu_k$ is critical to success of iteration. Strategy: choose $\mu$ slightly larger than $\|\lambda^*\|_{\infty}$ but that depends on lagrange multipler estimate.
    
    \vfill\null
    \columnbreak
    
    \subsection{Augmented Lagrangian}
    Add estimated Lagrange multipliers term to quadratic penalty formulation. Differs from standard Lagrangian with the extra quadratic terms.\\
    $Q(x, \lambda; \mu) = f(x) - \sum_i \lambda_i c_i(x) + \frac{\mu}{2} \sum_i c_i(x)^2, \forall i \in \mathcal{E}$\\
    
    Idea: fix $\mu_k, \lambda_k$ at current iteration, minimized wrt. $x$, update estimate of Lagrangian variables.

    $x_k = \argmin_x L(x, \lambda_k; \mu_k)$\\
    Optimality condition for uncstained minimization:\\
    $0 \approx \nabla_x L(x_k, \lambda_k; \mu_k) = \nabla f(x_k) - \sum_{i \in \mathcal{E}}(\lambda_{k,i} - \mu_k c_i(x_k)) \nabla c_i(x_k)$\\
    $\lambda_i^* \approx \lambda_{k,i} - \mu_k c_i(x_k), \forall i \in \mathcal{E}$\\
    Update of Lagrangian variables follow naturally.\\
    $\lambda_{k+1,i} \leftarrow \lambda_{k,i} - \mu_k c_i(x_k), \forall i \in \mathcal{E}$\\
    $c_i(x_k) \approx \frac{\lambda_{k,i} - \lambda_i^*}{\mu_k}$\\
    If $\lambda_{k,i}$ close to $\lambda_i^*$, then $c_i(x_k) \ll \frac{1}{\mu_k}$, this is better than quadratic penalty method's $c_i(x_k)\approx \frac{-\lambda_i^*}{\mu_k}, \forall i \in \mathcal{E}$\\

    Even with modest penalty weight, $\mu_k$, $c_i(x_k)$ is more like 0 which is require for optimality condition, this allows convergence without having $\mu_k \rightarrow +\infty$. Less ill conditioning and initial starting point is less critical.\\

    Tolerance $\tau_k$ can be chosen on infeasibility, $\sum_i |c(x_k)|$, and can be increased if infeasibility reduction in insufficient.
    
    \begin{algorithm}[H]
      $x_0^s = ..$\\
      $\lambda_0 = ..$\\
      $\tau_0 > 0$\\
      $\mu_0 > 0$\\
      \For{$k=0,1,..$}{
        solve $x_k = \argmin_x L(x, \lambda_k; \mu_k)$ starting at $x_k^s$ with termination condition $\|\nabla_x Q(x, \lambda_k;\mu_k)\| \leq \tau_k$\\
        \uIf{$c_i(x)=0, \forall i \in \mathcal{E}$}{
          return $x_k$
        }
        $\lambda_{k+1,i} \leftarrow \lambda_{k,i} - \mu_k c_i(\lambda_{k,i})\\
        update $\mu_{k+1}: \mu_{k+1} > \mu_k$ (eg: chosen adaptively)\\
        update $x_{k+1}^s \leftarrow x_k$\\
        update $\tau_{k+1}^s \leftarrow ...$
      }\\
      \caption{Augmented Lagrangian Algo \label{Algo_augmented_lagrangian}}
    \end{algorithm}

    Augmented Lagrangian Approaches:
    \begin{itemize}
    \item Bound Constained
    \item Linearly Constrained Lagrangian
    \item Unconstrained
    \end{itemize}

    \subsubsection{Bound Constrained}
    Introduce slack variables for inequalities, transform into equality constrained problem with variable bound constraints.\\

    $\min_x f(x),\ s.t.\ c_i(x) = 0\ \forall i,\ l \leq x \leq u$\\
    
    Solve the subprobem using augmented Lagrangian with nonlinear gradient probjection method by fixing $\lambda, \mu$.\\
    $\min_x L(x, \lambda_k; \mu_k) = f(x) - \sum_i \lambda_i c_i(x) + \frac{\mu}{2}\sum_i c_i(x)^2$\\
    
    1st order necessary conditions for $x$:\\
    $x - Proj_{[l,u]}(x-\nabla_x L(x,\lambda; \mu), l, u) = 0$\\

    TODO: [insert algo 17.4 nocedal]

    \subsubsection{Linearly Constrained Lagrangian}
    Idea: lienarize constraints and generate a step by minimizing Lagragian.
    \begin{align*}
      & \min_x F_k(x) = f(x) - \sum_i \lambda_i^k \bar{c}_i^k(x) + \frac{\mu}{2} \sum_i [\bar{c}_i^k(x)]^2\\
      s.t.\ & c(x_k) + A_k(x-x_k)=0\\
      & \bar{c}_i^k(x) = c_i(x) - (c_i(x_k) + \nabla c_i(x_k)^T (x-x_k))\\
      & l \leq x \leq u
    \end{align*}

    where $A_k$ is the constraint Jacobian.\\

    \subsubsection{Unconstrained}
    Idea: use proximal point approach for inequality constrained problems.
    \begin{align*}
      \min_x F(x) & = \max_{\lambda \geq 0} \{f(x) - \sum_{i \in \mathcal{I}} \lambda_i c_i(x) \} \\
      &= \begin{cases}
        f(x) &, if\ x\ feasible\\
        \infty &, o/w
      \end{cases}
    \end{align*}
    
    \begin{align*}
      \begin{cases}
        x\ \text{feasible} &, (c_i(x)<0),\ select\ \lambda_i\ large,\ \lambda_j = 0, \forall j \neq i\\
        & \implies F(x) \rightarrow \infty\\
        x\ \text{infeasible} &, (c_i(x) \geq 0), \forall i \in \mathcal{I},\ select\ (\forall i \in \mathcal{I})\lambda_i= 0\\
        & \implies F(x) = f(x)
      \end{cases}
    \end{align*}

    Then $\min_x F(x) = \min_{x\ feasible} f(x)$\\

    Use an approximation $\hat{F}(x; \lambda^k, \mu_k)$ to $F$ that is smooth for practical optimization:

    $\hat{F}(x; \lambda_^, \mu_k) = \max_{\lambda \geq o} \{ f(x) - \sum_{i \in \mathcal{I}} \lambda_i c_i(x) - \frac{1}{2 \mu_k} \sum_{i \in \mathcal{I}}(\lambda_i - \lambda_i^k)^2 \}$\\

    This is a bound constrained quadratic program in $lambda$, separable in individual components $\lambda_i$. The additional proximal term helps staying close to $\lambda_i^k$.\\

    We perform maximization explicitly:\\
    
    $\lambda_i = \begin{cases}
      0 &, -c_i(x) + \frac{\lambda_i^k}{\mu_k} \leq 0\\
      \lambda_i^k - \mu_k c_i(x) &, o/w
    \end{cases}$

    Substitute $\lambda_i$ back in $\hat{F}$:\\
    
    $\hat{F}(x; \lambda^k, \mu_k) = f(x) + \sum_{i \in \mathcal{I}} \psi(c_i(x), \lambda_i^k; \mu_k)$\\
    
    $\psi(t,\sigma; \mu) = \begin{cases}
      - \sigma t + \frac{\mu}{2} t^2 &, t - \frac{\sigma}{\mu} \leq 0\\
      - \frac{1}{2} \mu} \sigma^2} &, o/w
    \end{cases}$

    New iterate $x_k$ is found by minimizing $\hat{F}$:\\

    $x_K = \argmin_x \hat{F}(x; \lambda_k, \mu_k)$\\

    New Lagrangian multiplier update rule:\\
    
    $\lambda_i^{k+1} = \begin{cases}
      0 &, -c_i(x) + \frac{\lambda_i^k}{\mu_k} \leq 0\\
      \lambda_i^k - \mu_k c_i(x) &, o/w
    \end{cases}$
    
    \vfill\null
    \pagebreak

    \section{Linear System, Iterative}

    \subsection{Idea}
    Solve $Ax=b$\\
    Idea: split $A$ into $M+N$ and use back/forward substitutions to solve:\\
    $Mx^{k+1}=Nx^k+b$
    
    \subsection{Jacobi}
    \subsubsection{Splitting}
    \begin{align*}
      A &= M-N\\
      M &= D = Diagonal\\
      N &= M - A = D - A = -(L+U)\\
      L &= \text{Lower left}\\
      U &= \text{Upper right}\\
      (M-N)x &= b\\
      Mx &= Nx + b\\
      Dx &= -(L+U)x + b\\
      x &= D^{-1}(b -(L+U)x)
    \end{align*}

    \subsubsection{Algo}
    \begin{align*}    
      x^{k+1} &= D^{-1}(b -(L+U)x^k)
    \end{align*}
    
    \subsection{Gauss-Seidel}
    Use most recently updated values
    \subsubsection{Splitting}
    \begin{align*}
      A &= M-N\\
      M &= D + L\\
      N &= -U\\
      Mx &= Nx +b\\
      (D+L)x &= -Ux + b\\
      x &=  (D+L)^{-1}(b -Ux)
    \end{align*}

    \vfill\null
    \columnbreak
        
    \subsubsection{Algo}
    \begin{align*}
      x^{k+1} &=  (D+L)^{-1}(b -Ux^k)\\
      (D+L) x^{k+1} &= (b -Ux^k)\\
      (I+D^{-1}L) x^{k+1} &= D^{-1} (b -Ux^k)\\
      x^{k+1} &= D^{-1} (b -Ux^k) - D^{-1}L x^{k+1}\\
      x^{k+1} &= D^{-1} (b - L x^{k+1} - Ux^k)\\
      (\forall i)\ x_i^{k+1} & \leftarrow \frac{1}{A_{ii}} \left(b_i - \sum_{j<i} A_{ij} x_j^{k+1} - \sum_{j>i} A_{ij} x_j^k\right)
    \end{align*}
    
    \subsection{SOR(Successive Over-Relaxation)}
    \begin{align*}
      & x^{k+1} \leftarrow x^k + w (x_{GS}^{k+1} - x^k)\\
      & x^{k+1} \leftarrow (1-w) x^k + w x_{GS}^{k+1}\\
      &\begin {cases}
        w > 1 &, \text{over-relaxation}\\
        w \in (0,1) &, \text{under-relaxation}\\
        w = 1 &, \text{Gauss-Seidel}
      \end {cases}
    \end{align*}

    \subsubsection{Splitting}
    \begin{align*}
      & A = M - N\\
      & M = \frac{1}{w} D + L\\
      & N = (\frac{1}{w} - 1)D - U\\
      & Mx = Nx + b\\
      & x = M^{-1}(Nx + b)\\
      & (\frac{1}{w} D + L)x^{k+1} = ((\frac{1}{w}-1) D - U) w^k +b\\
      & x^{k+1} = (D+wL)^{-1} (((1-w)D - wU)x^k + wb)
    \end{align*}

    \subsubsection{Algo}
    \begin{align*}
      & (\frac{1}{w} D + L)x^{k+1} = ((\frac{1}{w}-1) D - U) w^k +b\\
      & (\frac{1}{w}I + D^{-1}L)x^{k+1} = ((\frac{1}{w} - 1)I - D^{-1}U)x^k + D^{-1}b\\
      & \frac{1}{w}x^{k+1} = -D^{-1}L x^{k+1} + ((\frac{1}{w} - 1)I - D^{-1}U)x^k + D^{-1}b\\
      & x^{k+1} = w(D^{-1}(b-Lx^{k+1}-Ux^k)+(\frac{1}{w}-1)x^k)\\
      & x^{k+1} = wD^{-1}(b-Lx^{k+1}-Ux^k)+(1-w)x^k\\
      & x^{k+1} = w x^{k+1}_{Gauss-Seidel} + (1-w)x^k
    \end{align*}
    
    \vfill\null
    \pagebreak
    
    \section{Factorization}
    
    \subsection{LU}
    
    By design:
    \begin{align*}
      U&=MA\\
      M^{-1} U &= A\\
      LU &= A\\
      L &= M^{-1} = (M_{n-1} .. M_1)^{-1}\\
      L &= M_1^{-1}..M_{n-1}^{-1} = L_1 .. L_{n-1}
    \end{align*}

    With partial pivoting in rows:
    \begin{align*}
      PAx &= Pb\\
      LU &= PA\\
      LUx &= Pb\\
      y &= Ux\ //solve\ for\ y\\
      Ly &= Pb\\
      Ux &= y\ //Solve\ for\ x
    \end{align*}

    \begin{algorithm}[H]
      $P \leftarrow \begin{bmatrix} 0 & .. & n-1 \end{bmatrix}$//permute\ rows\\
      $S \leftarrow 0$//row\ swaps\\
      \For{$col\ in\ 0\ ..\ n$}{
        let $rowmax \leftarrow $col\\
        \For{$row\ in\ col+1\ ..\ n$}{
          \uIf{$A_{row, col} > A_{rowmax,col}$}{
            $rowmax \leftarrow row$\\
          }
        }
        
        swap($A_{rowmax, :}$, $A_{col, :}$)\\
        swap($P_{rowmax}$, $P_{col}$)\\
        $S \leftarrow S +  rowmax\ !=\ col\ ?\ 1\ :\ 0$\\

        \uIf{$A_{col, col} < \epsilon_{tol}$}{
          raise\ degenerate\ matrix\ error
        }
        \For{$row\ in\ col+1\ ..\ n$}{
          $r \leftarrow A_{row, col} / A_{col,col}$\\
          $A_{row, col} \leftarrow r$ //save value for $L$\\
          //reduction for $U$\\
          \For{$c\ in\ col+1\ ..\ n$}{
            $A_{row,c} \leftarrow A_{row,c} - r A_{col, c}$\\
          }
        }
      }
      return $A, P, S$
      \caption{LU with Partial Pivot \label{Algo_LU_partial_pivot}}
    \end{algorithm}

    $L$ stored in lower left triangle of $A$ with implicit diagonal ones.\\
    
    $U$ stored in upper right triangle of $A$.\\
    
    \vfill\null
    \columnbreak

    Determinant:\\
    
    $LU=PA$\\
    $det(A) = det(P^{-1} LU)$\\
    $det(A) = det(P^{-1}) det(L) det(U)$\\
    $det(P^{-1}) = (-1)^S$\\
    where $S$ is the number of row exchanges from the permutation matrix $P$\\
    $det(A) = (-1)^S \prod_i L_{i,i} \prod_i U_{i,i}$\\
    With implicit diagonal ones of $L$ in Algo \ref{Algo_LU_partial_pivot} and storage of $U$ in returned $A$:\\
    $det(A) = (-1)^S \prod_i A_{i,i}$\\
    
    \vfill\null
    \pagebreak
    
    \subsection{QR}
    $A = QR$\\
    $A = \begin{bmatrix} Q_1 &Q_2 \end{bmatrix} \begin{bmatrix} R\\0 \end{bmatrix}$\\
    $Q$: m x m, orthogonal\\
    $R$: n x n, upper triangular\\
    m \geq n\\

    Reduced QR factorization of A:\\
    $A = Q_1 R$\\
    $Q_1$: m x n\\
    $R$: n x n, upper triangular\\

    $Ax = b$\\
    $Q_1 R x = b$\\
    $Rx = Q_1^{-1} b = Q_1^T b$\\
    solve via subsubstitution for x\\

    \subsection{Least Squares Problem}
    Orthogonal matrix is norm preserving for $L2$:\\
    $\|Qx\|_2^2 = x^T Q^T Q x = x^T x = \|x\|_2^2$
    \begin{align*}
      & Ax=b\\
      & r = b-Ax\\
      & \|r\|_2^2 = \|b-Ax\|_2^2\\
      & \|r\|_2^2 = \|b-Q \begin{bmatrix} R \\ 0\end{bmatrix}x\|_2^2\\
      & \|r\|_2^2 = \|Q^Tb - Q^T Q \begin{bmatrix} R \\ 0\end{bmatrix}x\|_2^2\\
      & \|r\|_2^2 = \|Q^Tb - I\begin{bmatrix} R \\ 0\end{bmatrix}x\|_2^2\\
      & Q^T b = \begin{bmatrix} c_1 \\ c_2 \end{bmatrix}\\
      & \|r\|_2^2 = \|c_1 - Rx\|_2^2 + \|c_2\|_2^2\\
    \end{align*}

    Solve $c_1 - Rx=0 \implies \|r\|_2^2 = \|c_2\|^2$

    \vfill\null
    \columnbreak
        
    \subsubsection{Householder Reflection Method}
    Idea: relfect vector so all coordinates but one disappear.\\
    let $| \alpha | = \|x\|_2$\\
    $u \leftarrow x - \alpha e_1$\\
    $v = \frac{u}{\|u\|}$\\
    $Q = I - 2 vv^T$
    \begin{align*}
      & (I - 2 \frac{vv^t}{v^Tv}) (I - 2 \frac{vv^t}{v^Tv})^T\\
      & = (I - 2 \frac{vv^t}{v^Tv}) (I - 2 \frac{vv^t}{v^Tv})\\    
      & = I - 4 \frac{vv^t}{v^Tv} + 4 \frac{vv^t vv^T}{v^Tv v^Tv}\\
      & = I - 4 \frac{vv^t}{v^Tv} + 4 \frac{vv^t}{v^Tv}\\
      & = I\\
      & \implies (I - 2 \frac{vv^t}{v^Tv})\ \text{orthongal, symmetric}
    \end{align*}
    
    $Q^T = Q^{-1} = Q$\\
    
    $Qx = \begin{bmatrix} \alpha \\ 0 \\ ..\\ 0 \end{bmatrix}$\\

    $A$: m x n, $m \geq n$\\
    
    $ Q_1 A =
    \begin{bmatrix}
      \alpha_1 & .. \\
      0 & A'
    \end{bmatrix}$\\

    $ Q_x =
    \begin{bmatrix}
      I_{k-1} & 0 \\
      0 & Q_k'
    \end{bmatrix}$\\    

    Transform to upper triangular matrix:\\
    $R = Q_t .. Q_1 A$\\
    $(Q_t .. Q_1)^T R = (Q_t .. Q_1)^T (Q_t .. Q_1) A$\\
    let $Q = (Q_t .. Q_1)^T = (Q_1 .. Q_t)$\\
    $Q R = Q Q^T A$\\
    $QR = (Q_1 .. Q_t) R = A$\\

    Solve $Ax=b$\\
    $QRx = b$\\
    $Rx = Q^{-1}b$\\
    $Rx = Q^T b$\\
    $Rx = (Q_1 .. Q_t)^T b$\\    
    $Rx = (Q_t^T .. Q_1^T) b$\\
    $Rx = (Q_t .. Q_1) b$\\
    Solve with backsubstitution\\

    \vfill\null
    \columnbreak
    
    \section{Eigenvalues & Eigenvectors}

    \subsection{Power Iteration}

    $x_k = A x_{k-1} = A^k x_0$\\
    $x_k = A^k \sum_j \alpha_j v_j$\\
    $x_k = \sum_j \alpha_j A^k v_j$\\
    $x_k = \sum_j \alpha_j \lambda_j^k v_j$\\
    $x_k = \alpha_1 \lambda_1^k v_1 + \sum_{j\neq1} \alpha_j \lambda_j^k v_j$\\
    $x_k = \lambda_1^k (\alpha_1 v_1 + \sum_{j\neq1} \alpha_j \frac{\lambda_j^k}{\lambda_1^k} v_j$\\
    $\lim_{k \rightarrow \infty} |(\frac{\lambda_j}{\lambda_1})^k| \rightarrow 0$\\
    $\lim_{k \rightarrow \infty} x_k = \lambda_1^k \alpha_1 v_1$

    \begin{algorithm}[H]
      let $x_0 \leftarrow ...$\\
      \For{$k=0..$}{
        $y_{k+1} = A x_k$\\
        $x_{k+1} = \frac{y_{k+1}}{\|y_{k+1}\|_{\infty}}$ //normalize
      }\\
      \caption{Power Iteration Algo \label{Algo_power_iteration}}
    \end{algorithm}
    Issue: failure to converge to linear comb. of eigenvectors corresponding to same eigenvalue.\\

    \subsection{Inverse Iteration}
    Inverse iteration: use factorization and $eig(A^{-1}) = \frac{1}{eig(A)}$\\
    
    For computing smallest eigenvalue:
    
    \begin{algorithm}[H]
      let $x_0 \leftarrow ...$\\
      \For{$k=0..$}{
        Solve for $y_{k+1}$ with factorization:\\
        $A y_{k+1} = x_k$\\
        $x_{k+1} = \frac{y_{k+1}}{\|y_{k+1}\|_{\infty}}$ //normalize
      }\\
      \caption{Inverse Iteration \label{Algo_inverse_iteration}}
    \end{algorithm}
    converge to eigenvector corresponding to smallest eigenvalue of $A$ = largest eigenvalue of $A^{-1}$\\

    \vfill\null
    \columnbreak
        
    \subsection{Shifting}
    Rayleigh Quotient for better eigenvalue estimation:\\
    $\lambda = \frac{x^T A x}{x^T x}$\\
    $\sigma_{smallest}(A-\sigma I) = \lambda - \sigma$, $\lambda$: eigenvalue of $A$ closest to $\sigma$\\
    use eigenvalue estimate from Rayleigh Quotient in inverse iteration to get better eigenvector\\
    
    \begin{algorithm}[H]
      let $x_0 \leftarrow$ guess of eigenvector\\
      \For{$k=0..$}{
        $\sigma_{k+1} = \frac{x_k^T A x_k}{x_k^T x_k}$//shift\\
        solve $y_{k+1}$: $(A-\sigma_{k+1}I)y_{k+1} = x_k$//next vector iterate\\
        $x_{k+1} = \frac{y_{k+1}}{\|y_{k+1}\|_{\infty}}$//normalize\\
      }\\
      \caption{Inverse Iteration with Shifting \label{Algo_inverse_iteration_shifting}}
    \end{algorithm}        
        
    \subsection{Simultaneous Iteration}
    
    \begin{algorithm}[H]
      let $X_0$: $n \times p$ matrix of rank $p$ (columns of linearly independent vectors)\\
      \For{$k=1..$}{
        $X_k = A X_{k-1}$\\
      }\\
      \caption{Simultaneous Iteration Algo \label{Algo_simultaneous_iteration}}
    \end{algorithm}    

    let $S_0 = span(X_0)$\\
    $S \coloneqq$ invariant subspace spanned by $p$ eigenvectors associated with largest eigenvalues of A.\\
    $S_k = A^k S_0$\\
    columns of $X_k$ form a basis for $S_k$\\
    $S_k \rightarrow S$\\

    Issue: $X_k$ beomes ill-conditioned bassi for subspace $S_k$.\\
    One solution: apply normalization such as QR factorization of columns of $X_k$:
    
    \begin{algorithm}[H]
      let $X_0$: $n \times p$ matrix of rank $p$ (columns of linearly independent vectors)\\
      \For{$k=1..$}{
        $\hat{Q}_k R_k = X_{k-1}$ //normalize\\
        $X_k = A \hat{Q}_k$ //gen next iterate\\
      }\\
      \caption{Orthogonal Iteration Algo \label{Algo_orthogonal_iteration}}
    \end{algorithm}

    $\hat{Q}_k:  n \times p$ orthogonal\\
    $R_k:  p \times p$ upper triangle\\

    \vfill\null
    \pagebreak
    
    \subsection{QR Iteration}
    Reference: section 4.5.6 from Heath\\
    
    $A\hat{Q} = \hat{Q}B$\\
    $B$ triangular\\

    Define $A_k = \hat{Q}_k^H A \hat{Q}_k$\\
    $\hat{Q}_k$ orthogonal $\implies$ $\sigma(A) = \sigma(\hat{Q}_k^H A \hat{Q}_k)$\\
    $\hat{Q}_k A_k = A \hat{Q}_k$\\
    $A \hat{Q}_k = \hat{Q}_k A_k$\\
    Shur Form of $A$:\\
    $A_k$ triangular / block triangular when it converges\\

    Avoiding explicit computation of $X_k$ and its factorization $\hat{Q}_k R_k = X_k$:\\
    let $X_0 = I$\\
    factorize: $\hat{Q}_0 R_0 = X_0$\\
    $\hat{Q}_0 = R_0 = I$\\
    $X_1 = A \hat{Q}_0 = A I = A$\\
    factorize: $\hat{Q}_1 R_1 = X_1 = A$\\
    using $A_k = \hat{Q}_k^H A \hat{Q}_k$:\\
    $A_1 = \hat{Q}_1^H A \hat{Q}_1$\\
    $A_1 = \hat{Q}_1^H \hat{Q}_1 R_1 \hat{Q}_1$, where $\hat{Q}_1 R_1 = A$\\
    $A_1 = R_1 \hat{Q}_1$ (notice reverse of $\hat{Q}_1 R_1 = A$)\\
    $X_2 = A \hat{Q}_1$\\
   
    $X_2 = \hat{Q}_1 \hat{Q}_1^H A \hat{Q}_1$\\
    $X_2 = \hat{Q}_1 A_1$\\

    instead of factorizing: $\hat{Q}_2 R_2 = X_2 = A \hat{Q}_1$\\
    factorize this: $Q_2 R_2 = A_1$:\\
    $X_2 = \hat{Q}_1 Q_2 R_2$\\
    let $\hat{Q}_2 = \hat{Q}_1 Q_2$\\
    $X_2 = \hat{Q}_2 R_2$\\

    $A_2 = \hat{Q}_2^H A \hat{Q}_2$\\
    $A_2 = \hat{Q}_2^H \hat{Q}_1 R_1 \hat{Q}_2$, where $A = \hat{Q}_1 R_1$\\
    since $\hat{Q}_2 = \hat{Q}_1 Q_2$\\
    $A_2 = Q_2^H \hat{Q}_1^H \hat{Q}_1 R_1 \hat{Q}_1 Q_2$\\
    $A_2 = Q_2^H R_1 \hat{Q}_1 Q_2$\\    
    since $R_1 \hat{Q}_1 = A_1$\\
    $A_2 = Q_2^H A_1 Q_2$\\
    since $Q_2 R_2 = A_1$\\
    $A_2 = Q_2^H Q_2 R_2 Q_2$\\
    $A_2 = R_2 Q_2$ (notice reverse of $Q_2 R_2 = A_1$)\\

    $A_3 = \hat{Q}_3^H A \hat{Q}_3$\\
    since $\hat{Q}_1 R_1 = A$\\
    $A_3 = \hat{Q}_3^H \hat{Q}_1 R_1 \hat{Q}_3$\\
    since $\hat{Q}_3 = \hat{Q}_1 Q_2 Q_3$\\
    $A_3 = Q_3^H Q_2^H \hat{Q}_1^H \hat{Q}_1 R_1 \hat{Q}_1 Q_2 Q_3$\\
    $A_3 = Q_3^H Q_2^H  R_1 \hat{Q}_1 Q_2 Q_3$\\
    since $A_1 = R_1 \hat{Q}_1$\\
    $A_3 = Q_3^H Q_2^H  A_1 Q_2 Q_3$\\
    since $Q_2 R_2 = A_1$\\
    $A_3 = Q_3^H Q_2^H Q_2 R_2 Q_2 Q_3$\\
    $A_3 = Q_3^H R_2 Q_2 Q_3$\\
    since $R_2 Q_2 = A_2$\\
    $A_3 = Q_3^H A_2 Q_3$\\
    factorize $Q_3 R_3 = A_2$\\
    $A_3 = Q_3^H Q_3 R_3 Q_3$\\
    $A_3 = R_3 Q_3$ (notice reverse of $Q_3 R_3 = A_2$)\\
    ...\\

    Thus $X_k$, and factoriaztion $\hat{Q}_k R_k = X_k$ is not calculated.\\
    Convergence is checked on $A_k$ for (block) triangular form, in which the eigenvalues can be read off.\\
    
    \begin{algorithm}[H]
      let $A_0 = A$\\
      \For{$k=1..$}{
        $Q_k R_k = A_{k-1}$ //normalize\\
        $A_k = R_k Q_k$ //gen next iterate\\
      }\\
      \caption{QR Iteration Algo \label{Algo_QR_iteration}}
    \end{algorithm}

    TODO 2
    
\end{multicols*}
\end {document}