【Optimal Control (CMU 16-745)】Lecture 5 Optimization Part 3

Review:

• Line search
• Equality and inequality constraints
• KKT conditions

Lecture 5 Optimization Pt.3

Overview

• Algorithms for constrained minimization
• Augumented lagrangian method
• Quadratic programming
• More on regularization and line search

1. Inequality-constrained minimization

$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})\le 0\end{array}$$

KKT conditions:

$$\begin{array}{rl}\nabla f\left(\mathbf{x}\right)+{\left(\frac{\partial \mathbf{c}}{\partial \mathbf{x}}\right)}^{\top }\lambda =0& (\mathrm{stationarity})\\ \mathbf{c}\left(\mathbf{x}\right)\le 0& (\mathrm{primal}\mathrm{feasibility})\\ \lambda \ge 0& (\mathrm{dual}\mathrm{feasibility})\\ \lambda \odot \mathbf{c}\left(\mathbf{x}\right)={\lambda }^{\top }\mathbf{c}\left(\mathbf{x}\right)=0& (\mathrm{complementary}\mathrm{slackness})\end{array}$$

• ${\left(\frac{\partial \mathbf{c}}{\partial \mathbf{x}}\right)}^{\top }\lambda$ acts like a penalty term. Make sure its sign is consistent (makes cost worse for $\mathbf{c}\left(\mathbf{x}\right)$ infeasible)
• The notation $\odot$ is the Hadamard product (element-wise product)

2. Optimization algorithms

(1) Active-set methods

• Guess active/inactive constraints
• Solve equality constrained problem
• Used when you have a good heuristic for the active set

(2) Barrier/Interior-point methods

Replace inequalities with barrier function in objective
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})\le 0\end{array}$$

$$\Rightarrow \underset{\mathbf{x}}{\min }f(\mathbf{x})-\sum _{i=1}^m\frac{1}{\rho }\log (-c_i(\mathbf{x}))$$

• Gold standard for convex problems. (e.g. MPC)

Example:
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})=x_1^2+0.3x_2^2\\ & \text{s.t.}{c}_1(\mathbf{x})=x_1^2+x_2^2-2\le 0\\ & c_2(\mathbf{x})=x_1+x_2-1\le 0\end{array}$$

Plot the objective and constraints on 2D plane:

Plot the objective and constraints in 3D:

(3) Penalty

Replace constraints with penalty term that penalizes violations
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})\le 0\end{array}$$

$$\Rightarrow \underset{\mathbf{x}}{\min }f(\mathbf{x})+\frac{\rho }{2}{\left[\max (0,c(\mathbf{x}))\right]}^2$$

• easy to implement
• has issues with numerical ill-conditioning
• cannot achieve high accuracy
• rarely used in practice

(4) Augmented Lagrangian

Add Lagrange multiplier estimate to penalty method:
$$\underset{\mathbf{x}}{\min }{\mathcal{L}}_{\rho }(\mathbf{x},\tilde{\lambda })=f(\mathbf{x})+{\tilde{\lambda }}^{\top }\mathbf{c}(\mathbf{x})+\frac{\rho }{2}{\left[\max (0,c(\mathbf{x}))\right]}^2$$

${\mathcal{L}}_{\rho }(\mathbf{x},\tilde{\lambda })$ is called the augmented Lagrangian. $\tilde{\lambda }$ is the estimate of the Lagrange multiplier.

i. Update $\tilde{\lambda }$ by offloading the penalty term into $\tilde{\lambda }$ at each iteration:
$$\tilde{\lambda }\leftarrow \tilde{\lambda }+\rho \mathbf{c}(\mathbf{x})$$(for active constraints)

Insight:
$$\begin{array}{rl}& \frac{\partial f}{\partial \mathbf{x}}+{\tilde{\lambda }}^{\top }\frac{\partial \mathbf{c}}{\partial \mathbf{x}}++\rho {\mathbf{c}}^{\top }(\mathbf{x})\frac{\partial \mathbf{c}}{\partial \mathbf{x}}\\ =& \frac{\partial f}{\partial \mathbf{x}}+{\left[\tilde{\lambda }+\rho \mathbf{c}(\mathbf{x})\right]}^{\top }\frac{\partial \mathbf{c}}{\partial \mathbf{x}}\\ =& 0\end{array}$$ The term $\tilde{\lambda }+\rho c(\mathbf{x})$ looks like a Lagrange multiplier. It pushes at the same direction as the Lagrange multiplier, and pushes up against the constraint.

ii. Repeat until convergence:
(1) ${\min }_{\mathbf{x}}{\mathcal{L}}_{\rho }(\mathbf{x},\tilde{\lambda })$
(2) $\tilde{\lambda }\leftarrow \max (0,\tilde{\lambda }+\rho \mathbf{c}(\mathbf{x}))$ (clamping to guarantee $\tilde{\lambda }\ge 0$)
(3) $\rho \leftarrow \alpha \rho$ ($\alpha$ typically 10)

• Fixes ill-conditioning of penalty method
• Converges with finite $\rho$
• Works well on non-convex problems

Example: Quadratic program

$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Qx}+{\mathbf{q}}^{\top }\mathbf{x}(\mathbf{Q}\succ 0)\\ & \text{s.t.}\mathbf{Ax}\le \mathbf{b},\mathbf{Cx}=\mathbf{d}\end{array}$$

• Super useful in control
• Can be solved very fast (~kHz)

Parameters: $\mathbf{Q}=\left[\begin{array}{cc}0.5& 0\\ 0& 1\end{array}\right]$, $\mathbf{q}=0$, $\mathbf{A}=\left[\begin{array}{cc}1& 1\end{array}\right]$, $\mathbf{b}=-1$, $\mathbf{C}=0$, $\mathbf{d}=0$.

Plot the objective and constraints on 2D plane:

After 3 iterations:

Without $\rho$ update (slow convergence):

Example:

• Try with penalty method, full Augmented Lagrangian, and just $\lambda$ update

3. Regularization & Duality

(1) Equality constraints

Given
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})=0\end{array}$$

We might like to turn this into:
$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+p_{\infty }(\mathbf{c}(\mathbf{x})),p_{\infty }(\mathbf{x})=\{\begin{array}{ll}0& \mathbf{x}=0\\ +\infty & \text{otherwise}\end{array}$$

Practically it’s terrible, but we can get the same effect solving:
$$\underset{\mathbf{x}}{\min }\underset{\lambda }{\max }f(\mathbf{x})+{\lambda }^{\top }\mathbf{c}(\mathbf{x})$$

If the constraints cannot be satisfied ($\mathbf{c}(\mathbf{x})\ne 0$), then the cost is infinite.
The $\max$ acts like an infinite penalty.

(2) Inequality constraints

The above trick is also similar for inequality constraints:
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})\le 0\end{array}$$

$$\Rightarrow \underset{\mathbf{x}}{\min }f(\mathbf{x})+p_{\infty }(\mathbf{c}(\mathbf{x})),p_{\infty }(\mathbf{x})=\{\begin{array}{ll}0& \mathbf{x}\le 0\\ +\infty & \text{otherwise}\end{array}$$

Apply the “min-max” trick:
$$\Rightarrow \underset{\mathbf{x}}{\min }\underset{\lambda \ge 0}{\max }f(\mathbf{x})+{\lambda }^{\top }\mathbf{c}(\mathbf{x})$$

We denote
$$\mathcal{L}(\mathbf{x},\lambda )=f(\mathbf{x})+{\lambda }^{\top }\mathbf{c}(\mathbf{x})$$

The solution is a saddle point in $\left(\mathbf{x},\lambda \right)$ space.

Consider the eigenvalues of the KKT matrix: The eigenvalues related to $\mathbf{x}$ are all positive, and the eigenvalues related to $\lambda$ should be all negative.

(3) Interpretation

• KKT conditions define a saddle point in $\left(\mathbf{x},\lambda \right)$ space.
• KKT system should have $\dim (\mathbf{x})$ positive eigenvalues and $\dim (\lambda )$ negative eigenvalues at an optimum. It’s called a quasi-definite system.
• Regularize the KKT matrix as follows:

$$\left[\begin{array}{cc}H+\beta I& c^{\top }\\ c& -\beta I\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \end{array}\right]=\left[\begin{array}{c}-{\nabla }_{\mathbf{x}}\mathcal{L}\\ -\mathbf{c}\left(\mathbf{x}\right)\end{array}\right],\beta >0$$

(4) Duality

Question: Can we swap the $\min$ and $\max$?
Answer: If the problem is convex, yes! (Strong duality)

• For convex problems, the optimal value of the primal problem is equal to the optimal value of the dual problem.
• For non-convex problems, the optimal value of the primal problem is less than or equal to the optimal value of the dual problem. The difference is called the duality gap.

(5) Example

Without regularization, the solution will not converge to the optimal solution:

Check the eigenvalues of the KKT matrix:

H = ∇2f(xguess[:,end]) + ForwardDiff.jacobian(xn -> ∂c(xn)'*λguess[end], xguess[:,end])
C = ∂c(xguess[:,end])
K = [H C'; C 0]
eigvals(K)

We get:

3-element Vector{Float64}:
 -2.7504046059655027
 -0.5914642540616377
  1.6230587398907093

where two eigenvalues are negative and one is positive.

To fix this, we add regularization to the KKT matrix:

function regularized_newton_step(x,λ)
    β = 1.0
    H = ∇2f(x) + ForwardDiff.jacobian(xn -> ∂c(xn)'*λ, x)
    C = ∂c(x)
    K = [H C'; C 0]
    e = eigvals(K)
    while !(sum(e .> 0) == length(x) && sum(e .< 0) == length(λ)) # check if the KKT matrix is quasi-definite
        K = K + Diagonal([β*ones(length(x)); -β*ones(length(λ))]) # add regularization
        e = eigvals(K)
    end
    Δz = K$$-∇f(x)-C'*λ; -c(x)]
    Δx = Δz[1:2]
    Δλ = Δz[3]
    return Δx, Δλ
end

The solution converges to the optimal solution in 6 iterations:

• still overshoot (the 4th iteration) $\Rightarrow$ use line search

4. Merit function

(1) Question: How do we do a line search on a root-finding problem?

$\Rightarrow$ We find ${\mathbf{x}}^{\ast }$ such that $\mathbf{c}({\mathbf{x}}^{\ast })=0$.

We define a scalar merit function $p\left(\mathbf{x}\right)$ that measures the distance to solution.

Here are some standard choices:

• $L_2$ norm: $p\left(\mathbf{x}\right)=\frac{1}{2}{\mathbf{c}}^{\top }(\mathbf{x})\mathbf{c}(\mathbf{x})=\frac{1}{2}\Vert \mathbf{c}(\mathbf{x}){\Vert }_2^2$
• $L_1$ norm: $p\left(\mathbf{x}\right)=\Vert \mathbf{c}(\mathbf{x}){\Vert }_1$ (any norm works)

Now just do Armijo line search on $p\left(\mathbf{x}\right)$:
$\alpha =1$ (step length)
while $p\left(\mathbf{x}+\alpha \Delta \mathbf{x}\right)>p\left(\mathbf{x}\right)+b\alpha \nabla p{\left(\mathbf{x}\right)}^{\top }\Delta \mathbf{x}$
\quad $\alpha \leftarrow \theta \alpha$ ($\theta \in (0,1)$)
end
$\mathbf{x}\leftarrow \mathbf{x}+\alpha \Delta \mathbf{x}$

(2) Question: How about constrained minimization?

$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t.}\mathbf{c}(\mathbf{x})\le 0,\mathbf{d}(\mathbf{x})=0\end{array}$$

$$\Rightarrow \underset{\mathbf{x}}{\min }\mathcal{L}(\mathbf{x},\lambda ,\mu )=f(\mathbf{x})+{\lambda }^{\top }\mathbf{c}(\mathbf{x})+{\mu }^{\top }\mathbf{d}(\mathbf{x})$$

Pick a merit function:
$$p(\mathbf{x},\lambda ,\mu )=\frac{1}{2}\Vert \mathbf{r}(\mathbf{x},\lambda ,\mu ){\Vert }_2^2,\mathbf{r}(\mathbf{x},\lambda ,\mu )=\left[\begin{array}{c}{\nabla }_{\mathbf{x}}\mathcal{L}(\mathbf{x},\lambda ,\mu )\\ \min (0,\mathbf{c}(\mathbf{x}))\\ \mathbf{d}(\mathbf{x})\end{array}\right]$$ ($\mathbf{r}(\mathbf{x},\lambda ,\mu )$ is called the KKT residual)

or
$$p(\mathbf{x},\lambda ,\mu )=f(\mathbf{x})+\rho {\left|\left|\left[\begin{array}{c}\min (0,\mathbf{c}(\mathbf{x}))\\ \mathbf{d}(\mathbf{x})\end{array}\right]\right|\right|}_1$$

or the Augmented Lagrangian:
$$p(\mathbf{x},\lambda ,\mu )={\mathcal{L}}_{\rho }(\mathbf{x},\lambda ,\mu )$$

(3) Example

See the jupyter notebook.