Augmented Lagrangian method

Definitions

The nonlinear program we wish to solve is of the form:

$$\begin{array}{}\text{(P)}& \begin{array}{rlrlrlrl}& \underset{x}{\text{minimize}}& & f(x)& & & & f:{\mathrm{I}\mathrm{R}}^n\to \mathrm{I}\mathrm{R}\\ & \text{subject to}& & \underset{―}{x}\le x\le \bar{x}\\ & & & \underset{―}{z}\le g(x)\le \bar{z}& & & & g:{\mathrm{I}\mathrm{R}}^n\to {\mathrm{I}\mathrm{R}}^m\end{array}\end{array}$$

Define the convex sets $C$ and $D$ as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}C& =\{x\in {\mathrm{I}\mathrm{R}}^n\mid \underset{―}{x}\le x\le \bar{x}\}\\ D& =\{z\in {\mathrm{I}\mathrm{R}}^m\mid \underset{―}{z}\le z\le \bar{z}\}.\end{array}\end{array}$$

These rectangular boxes can be decomposed as Cartesian products of 1-dimensional closed intervals:

$$\begin{array}{}\text{(2)}& \begin{array}{rlrl}C& =C_1\times C_2\times \cdots \times C_n& & \text{where }{C}_i=[{\underset{―}{x}}_i,{\bar{x}}_i]\\ D& =D_1\times D_2\times \cdots \times D_m& & \text{where }{D}_i=[{\underset{―}{z}}_i,{\bar{z}}_i]\end{array}\end{array}$$

Using these definitions, problem $\text{(P)}$ can equivalently be expressed as

$$\begin{array}{}\text{(3)}& \begin{array}{rlrl}& \underset{x\in C}{\text{minimize}}& & f(x)\\ & \text{subject to}& & g(x)\in D.\end{array}\end{array}$$

After introduction of a slack variable $z$, problem $\text{(3)}$ can be stated as

$$\begin{array}{}\text{(P-ALM)}& \begin{array}{rlrl}& \underset{x\in C,z\in D}{\text{minimize}}& & f(x)\\ & \text{subject to}& & g(x)-z=0.\end{array}\end{array}$$

The Lagrangian function of problem $\text{(P-ALM)}$ is given by:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathcal{L}:{\mathrm{I}\mathrm{R}}^n\times {\mathrm{I}\mathrm{R}}^m\times {\mathrm{I}\mathrm{R}}^m\to \mathrm{I}\mathrm{R}:(x,z,y)& \mapsto \mathcal{L}(x,z,y)\\ & \triangleq f(x)+⟨g(x)-z,y⟩.\end{array}\end{array}$$

The vector $y\in {\mathrm{I}\mathrm{R}}^m$ is called the vector of Lagrange multipliers.

The augmented Lagrangian function with penalty factor $\mathrm{\Sigma }$ of the problem $\text{(P-ALM)}$ is defined as the sum of the Lagrangian function and a quadratic term that penalizes the constraint violation:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathcal{L}}_{\mathrm{\Sigma }}:{\mathrm{I}\mathrm{R}}^n\times {\mathrm{I}\mathrm{R}}^m\times {\mathrm{I}\mathrm{R}}^m\to \mathrm{I}\mathrm{R}:(x,z,y)& \mapsto {\mathcal{L}}_{\mathrm{\Sigma }}(x,z,y)\\ & \triangleq \mathcal{L}(x,z,y)+\frac{1}{2}{\Vert g(x)-z\Vert }_{\mathrm{\Sigma }}^2,\end{array}\end{array}$$

where $\mathrm{\Sigma }$ is a symmetric positive definite $m\times m$ matrix that defines a norm on ${\mathrm{I}\mathrm{R}}^m$, $\Vert z{\Vert }_{\mathrm{\Sigma }}^2\triangleq z^{\mathrm{\top }}\mathrm{\Sigma }z$.

The augmented Lagrangian method algorithm

The augmented Lagrangian method for solving Problem $\text{(P-ALM)}$ consists of the successive minimization of ${\mathcal{L}}_{\mathrm{\Sigma }}$ with respect to the decision variables $x$ and the slack variables $z$ (1), after which the Lagrange multipliers $y$ are updated (2), and the penalty factors ${\mathrm{\Sigma }}_{ii}$ corresponding to constraints with high violation are increased (3).

The augmented Lagrangian function is used as an exact penalty function for problem $\text{(P-ALM)}$, it is equivalent to the shifted quadratic penalty method with shift ${\mathrm{\Sigma }}^{-1}y$.

1. Minimization of the augmented Lagrangian

Using some algebraic manipulations, the augmented Lagrangian defined in $\text{(5)}$ can be expressed as

$$\begin{array}{}\text{(6)}& {\mathcal{L}}_{\mathrm{\Sigma }}(x,z,y)=f(x)+\frac{1}{2}\Vert g(x)-z+{\mathrm{\Sigma }}^{-1}y{\Vert }_{\mathrm{\Sigma }}^2-\frac{1}{2}\Vert y{\Vert }_{{\mathrm{\Sigma }}^{-1}.}^2\end{array}$$

At each iteration $\nu$ of the ALM algorithm, the following minimization problem is solved:

$$\begin{array}{}\text{(7)}& (x^{\nu },z^{\nu })=\underset{x\in C,z\in D}{\text{argmin}}{\mathcal{L}}_{{\mathrm{\Sigma }}^{\nu -1}}(x,z;y^{\nu -1})\end{array}$$

2. Update of the Lagrange multipliers

The update of the Lagrange multipliers corrects the shift ${\mathrm{\Sigma }}^{-1}y$ in $\text{(6)}$: if the constraint violation $g(x^{\nu })-z^{\nu }$ is positive, the shift is increased, in an attempt to drive the next iterate towards a smaller constraint violation $g(x^{\nu })-z^{\nu }$. The following update rule formalizes that idea:

$$\begin{array}{}\text{(8)}& y^{\nu }\leftarrow y^{\nu -1}+{\mathrm{\Sigma }}^{\nu -1}(g(x^{\nu })-z^{\nu })\end{array}$$

When the constraint violation becomes zero, the Lagrange multipliers are no longer updated.

As the penalty factors $\mathrm{\Sigma }$ tend towards infinity, the shift ${\mathrm{\Sigma }}^{-1}y$ has to vanish, because in that case, the quadratic penalty method without shifts solves the problem exactly. For ${\mathrm{\Sigma }}^{-1}y$ to vanish, the Lagrange multipliers must be bounded, which is achieved by the following projection:

Let $M>0$ be some large but finite bound.

$$\begin{array}{}\text{(9)}& {\underset{―}{y}}_i\triangleq \{\begin{array}{ll}0& {\underset{―}{z}}_i=-\mathrm{\infty }\\ -M& \text{otherwise},\end{array}{\bar{y}}_i\triangleq \{\begin{array}{ll}0& {\bar{z}}_i=+\mathrm{\infty }\\ +M& \text{otherwise}\end{array}\text{(10)}& Y\triangleq [{\underset{―}{y}}_1,{\bar{y}}_1]\times \cdots \times [{\underset{―}{y}}_m,{\bar{y}}_m]\end{array}$$

The result of $\text{(8)}$ is therefore clamped as follows:

$$\begin{array}{}\text{(11)}& y^{\nu }\leftarrow {\mathrm{\Pi }}_Y(y^{\nu -1}+{\mathrm{\Sigma }}^{\nu -1}(g(x^{\nu })-z^{\nu }))\end{array}$$

3. Update of the penalty factors

When the penalty factor for the $i$-th constraint, ${\mathrm{\Sigma }}_{ii}$ is increased, minimizing the violation of this constraint becomes more important in $\text{(7)}$. Therefore, if the constraint violation cannot be reduced by updating the shifts alone, the penalty factors are increased.

Selecting when and by how much each penalty factor should be increased is more of a heuristic. The strategy used here is to compare the violation at the current iterate with the violation at the previous iterate, it is the same strategy as used in QPALM. Denote the vector of constraint violations as $e^{\nu }\triangleq g(x^{\nu })-z^{\nu }$. Let $\theta \in (0,1)$.
If $|e_i^{\nu }|\le \theta |e_i^{\nu -1}|$, meaning that the constraint violation has decreased by at least a factor $\theta$ compared to the previous iteration, then the penalty factor is not updated.
If the constraint violation did not decrease sufficiently, then the penalty factor ${\mathrm{\Sigma }}_{ii}$ is increased by a factor

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }{\displaystyle \frac{|e_i^{\nu }|}{\Vert e^{\nu }{\Vert }_{\mathrm{\infty }}}},\end{array}$$

where $\mathrm{\Delta }>1$ is a tuning parameter. The violation of each individual constraint is scaled by the maximum violation of all constraints, such that the penalty factors of constraints with a large violation are increased more aggressively. If the factor in $\text{(12)}$ is less than one, the penalty factor is not updated (otherwise it would result in a reduction of the penalty).

PANOC

PANOC is an algorithm that solves optimization problems of the form:

$$\begin{array}{}\text{(P-PANOC)}& \begin{array}{rlrl}& \underset{x}{\text{minimize}}& & \psi (x)+h(x),\end{array}\end{array}$$

where $\psi :{\mathrm{I}\mathrm{R}}^n\to \mathrm{I}\mathrm{R}$ has Lipschitz gradient, and $h:{\mathrm{I}\mathrm{R}}^n\to \overline{\mathrm{I}\mathrm{R}}$ allows efficient computation of the proximal operator.

Recall the inner minimization problem $\text{(7)}$ in the first step of the ALM algorithm. It can be simplified to:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\underset{x\in C,z\in D}{min}{\mathcal{L}}_{\mathrm{\Sigma }}(x,z,y)& =-\frac{1}{2}\Vert y{\Vert }_{{\mathrm{\Sigma }}^{-1}}^2+\underset{x\in C}{min}\{f(x)+\underset{z\in D}{min}\left\{\frac{1}{2}{\Vert z-(g(x)+{\mathrm{\Sigma }}^{-1}y)\Vert }_{\mathrm{\Sigma }}^2\right\}\}\\ & =-\frac{1}{2}\Vert y{\Vert }_{{\mathrm{\Sigma }}^{-1}}^2+\underset{x\in C}{min}\{\underset{\triangleq {\psi }_{\mathrm{\Sigma }}(x;y)}{\underbrace{f(x)+\frac{1}{2}{\text{dist}}_{\mathrm{\Sigma }}^2(g(x)+{\mathrm{\Sigma }}^{-1}y,D)}}\}\end{array}\end{array}$$

Within the PANOC algorithm, the parameters $y$ and $\mathrm{\Sigma }$ remain constant, and will be omitted from the function names to ease notation:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\psi (x)& =f(x)+\frac{1}{2}{\text{dist}}_{\mathrm{\Sigma }}^2(g(x)+{\mathrm{\Sigma }}^{-1}y,D)\end{array}\end{array}$$

The inner problem in $\text{(13)}$ has the same minimizers as the following problem that will be solved using the PANOC algorithm:

$$\begin{array}{}\text{(15)}& \begin{array}{rlrl}& \underset{x\in C}{\text{minimize}}& & \psi (x),\end{array}\end{array}$$

This problem is an instance of problem $\text{(P-PANOC)}$ where the nonsmooth term $h$ is the indicator of the set $C$, $h(x)={\delta }_C(x)$.

Evaluation

The following is a list of symbols and formulas that are used in the implementation of the PANOC algorithm.

$$\begin{array}{rlr}y& \in {\mathrm{I}\mathrm{R}}^m& \text{Current Lagrange multipliers}\\ \mathrm{\Sigma }& \in \text{diag}({\mathrm{I}\mathrm{R}}_{>0}^m)& \text{Current penalty factor}\\ x^k& \in {\mathrm{I}\mathrm{R}}^n& \text{Current PANOC iterate}\\ {\gamma }_k& \in {\mathrm{I}\mathrm{R}}_{>0}& \text{Current proximal gradient step size}\\ {\zeta }^k& \triangleq g(x^k)+{\mathrm{\Sigma }}^{-1}y& \text{Shifted constraint value}\\ {\hat{z}}^k& \triangleq {\mathrm{\Pi }}_D(g(x^k)+{\mathrm{\Sigma }}^{-1}y)& \text{Closest feasible value for slack variable }z\\ & ={\mathrm{\Pi }}_D({\zeta }^k)\\ d^k& \triangleq {\zeta }^k-{\mathrm{\Pi }}_D({\zeta }^k)& \text{How far the shifted constraint value }\zeta \\ & ={\zeta }^k-{\hat{z}}^k& \text{is from the feasible set}\\ e^k& \triangleq g(x^k)-{\hat{z}}^k& \text{Constraint violation}\\ {\hat{y}}^k& \triangleq \mathrm{\Sigma }{d}^k& \text{Candidate Lagrange multipliers,}\\ & =\mathrm{\Sigma }(g(x^k)+{\mathrm{\Sigma }}^{-1}y-{\mathrm{\Pi }}_D(g(x^k)+{\mathrm{\Sigma }}^{-1}y))& \text{see }\text{(8)}\\ & =y+\mathrm{\Sigma }(g(x^k)-{\hat{z}}^k)\\ & =y+\mathrm{\Sigma }{e}^k\\ \psi (x^k)& ={\mathcal{L}}_{\mathrm{\Sigma }}(x^k,{\hat{z}}^k,y)+\frac{1}{2}\Vert y{\Vert }_{{\mathrm{\Sigma }}^{-1}}^2& \text{PANOC objective function}\\ & =f(x^k)+\frac{1}{2}{\mathbf{dist}}_{\mathrm{\Sigma }}^2(g(x^k)+{\mathrm{\Sigma }}^{-1}y,D)\\ & =f(x^k)+\frac{1}{2}{\Vert (g(x^k)+{\mathrm{\Sigma }}^{-1}y)-{\mathrm{\Pi }}_D(g(x^k)+{\mathrm{\Sigma }}^{-1}y)\Vert }_{\mathrm{\Sigma }}^2\\ & =f(x^k)+\frac{1}{2}{\Vert {\zeta }^k-{\hat{z}}^k\Vert }_{\mathrm{\Sigma }}^2\\ & =f(x^k)+\frac{1}{2}⟨d^k,{\hat{y}}^k⟩\\ \mathrm{\nabla }\psi (x^k)& =\mathrm{\nabla }f(x^k)+\mathrm{\nabla }g(x^k)\mathrm{\Sigma }(g(x^k)+{\mathrm{\Sigma }}^{-1}y-{\mathrm{\Pi }}_D(g(x^k)+{\mathrm{\Sigma }}^{-1}y))& \text{Gradient of the objective}\\ & =\mathrm{\nabla }f(x^k)+\mathrm{\nabla }g(x^k)\mathrm{\Sigma }({\zeta }^k-{\hat{z}}^k)\\ & =\mathrm{\nabla }f(x^k)+\mathrm{\nabla }g(x^k){\hat{y}}^k\\ {\hat{x}}^k& \triangleq T_{{\gamma }^k}\left(x^k\right)& \text{Next proximal gradient iterate}\\ & ={\mathrm{\Pi }}_C(x^k-{\gamma }^k\mathrm{\nabla }\psi (x^k))\\ p^k& \triangleq {\hat{x}}^k-x^k& \text{Proximal gradient step}\\ r^k& \triangleq \frac{1}{{\gamma }^k}{p}^k& \text{Fixed-point residual (FPR)}\\ {\phi }_{{\gamma }^k}(x^k)& =\psi (x^k)+h({\hat{x}}^k)+\frac{1}{2{\gamma }^k}\Vert {\hat{x}}^k-x^k{\Vert }^2+\mathrm{\nabla }\psi (x^k{)}^{\mathrm{\top }}({\hat{x}}^k-x^k)& \text{Forward-backward envelope (FBE)}\\ & =\psi (x^k)+\frac{1}{2{\gamma }^k}\Vert p^k{\Vert }^2+\mathrm{\nabla }\psi (x^k{)}^{\mathrm{\top }}{p}^k\\ q^k& \triangleq H_k{r}^k& \text{Quasi-Newton step}\\ x^{k+1}& =x^k+(1-\tau )p^k+\tau q^k& \text{Next PANOC iterate}\end{array}$$

Note that many of the intermediate values depend on the value of $x^k$, it is sometimes easiest to define them as functions:

$$\begin{array}{rl}\zeta (x)& \triangleq g(x)+{\mathrm{\Sigma }}^{-1}y\\ \hat{z}(x)& \triangleq {\mathrm{\Pi }}_D(g(x)+{\mathrm{\Sigma }}^{-1}y)\\ & ={\mathrm{\Pi }}_D(\zeta (x))\\ d(x)& \triangleq \zeta (x)-{\mathrm{\Pi }}_D(\zeta (x))\\ & =\zeta (x)-\hat{z}(x)\\ \hat{y}(x)& \triangleq \mathrm{\Sigma }d(x)\\ & =\mathrm{\Sigma }(g(x)+{\mathrm{\Sigma }}^{-1}y-{\mathrm{\Pi }}_D(g(x)+{\mathrm{\Sigma }}^{-1}y))\\ e(x)& \triangleq g(x)-\hat{z}(x)\end{array}$$

The result of the PANOC algorithm is the triple $({\hat{x}}^k,\hat{y}({\hat{x}}^k),\hat{z}({\hat{x}}^k))$.

The following graph visualizes the dependencies between the different values used in a PANOC iteration.