<jupyter_start><jupyter_text>Constrained convex optimization*Selected Topics in Mathematical Optimization: 2017-2018***<NAME>** ([email](<EMAIL>))![](Figures/logo.png)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from teachingtools import plot_contour, add_path
# suppress warnings
import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Equality constrained convex optimization Problem outlineWe will start with convex optimization problems with linear equality constraints:$$\min_\mathbf{x} f(\mathbf{x}) \\\text{subject to } A\mathbf{x}=\mathbf{b}$$where $f : \mathbb{R}^n \rightarrow \mathbb{R}$ is convex and twice continuously differentiable and $A\in \mathbb{R}^{p\times n}$ with a rank $p < n$.The Lagrangian of this problem is$$L(\mathbf{x}, \boldsymbol{\nu}) = f(\mathbf{x}) + \boldsymbol{\nu}^\top(A\mathbf{x}-\mathbf{b})\,,$$with $\boldsymbol{\nu}\in\mathbb{R}^p$ the vector of Lagrange multipliers.A point $\mathbf{x}^\star\in$ **dom** $f$ is optimal for the above optimization problem only if there is a $\boldsymbol{\nu}^\star\in\mathbb{R}^p$ such that:$$A\mathbf{x}^\star = \mathbf{b}, \qquad \nabla f(\mathbf{x}^\star) + A^\top\boldsymbol{\nu}^\star = 0\,.$$We will reuse the same toy examples from the previous chapter, but add an equality constraint to both.- Simple quadratic problem:$$ \min_{\mathbf{x}} \frac{1}{2} (x_1^2 + 4 x_2^2)\\ \text{subject to }  x_1 - 2x_2 = 3$$- A non-quadratic function:$$ \min_{\mathbf{x}}\log(e^{x_1 +3x_2-0.1}+e^{x_1 -3x_2-0.1}+e^{-x_1 -0.1})\\ \text{subject to }  x_1 + 3x_2 = 0  $$<jupyter_code>from teachingtools import quadratic, grad_quadratic, hessian_quadratic
from teachingtools import nonquadratic, grad_nonquadratic, hessian_nonquadratic
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))
plot_contour(quadratic, (-11, 11), (-5, 5), ax1, [1.0, -2.0], 3)
ax1.set_title('Quadratic')
ax1.legend(loc=0)
plot_contour(nonquadratic, (-2, 2), (-1, 1), ax2, [1, 3], 0)
ax2.set_title('Non-quadratic')
ax2.legend(loc=0)
ax1.grid()<jupyter_output><empty_output><jupyter_text>Equality constrained convex quadratic optimizationConsider the following equality constrained convex optimization problem:$$\min_\mathbf{x}\frac{1}{2}\mathbf{x}^\top P \mathbf{x} + \mathbf{q}^\top \mathbf{x} + r  \\\text{subject to }  A\mathbf{x}=\mathbf{b}$$where $P$ is symmetric.The optimality conditions are$$A\mathbf{x}^\star = \mathbf{b}, \quad P\mathbf{x}^\star+\mathbf{q} +A^\top\boldsymbol{\nu}^\star=\mathbf{0}\,,$$which we can write as$$\begin{bmatrix}P & A^\top \\A & 0 \\     \end{bmatrix}     \begin{bmatrix}\mathbf{x}^\star\\\boldsymbol{\nu}^\star     \end{bmatrix}     =     \begin{bmatrix}-\mathbf{q} \\\mathbf{b}     \end{bmatrix}\,.$$Note that this is a block matrix.> If $P$ is positive-definite, the linearly constrained quadratic minimization problem has an unique solution.Solving this linear system gives both the constrained minimizer $\mathbf{x}^\star$ as well as the Lagrange multipliers. **Assignment 1**1. Complete the code to solve linearly constrained quadratic systems.2. Use this code to solve the quadratic toy problem defined above.<jupyter_code>from constrained import solve_constrained_quadratic_problem<jupyter_output><empty_output><jupyter_text>def solve_constrained_quadratic_problem(P, q, A, b):    """    Solve a linear constrained quadratic convex problem.    Inputs:        - P, q: quadratic and linear parameters of                the linear function to be minimized        - A, b: system of the linear constraints    Outputs:        - xstar: the exact minimizer        - vstar: the optimal Lagrange multipliers    """    p, n = ...   size of the problem     complete this code     HINT: use np.linalg.solve and np.bmat    solution = ...    xstar = solution[:n]    vstar = solution[n:]    return np.array(xstar), np.array(vstar)  complete the parametersP = ...    2 x 2 arrayA = ...    1 x 2 arrayb = ...   1 x 1 arrayq = ...    2 x 1 array xstar, vstar = solve_constrained_quadratic_problem(...   complete print('Minimizer:')print(xstar)fig, ax = plt.subplots()plot_contour(quadratic, (-11, 11), (-5, 5), ax, [1, -2], 3)ax.scatter(xstar[0,0], xstar[1,0], 50, 'r', label='minimum')ax.grid()ax.legend(loc=0)<jupyter_code>### Newton's method with equality constraints

To derive $\Delta \mathbf{x}_{nt}$ for the following equality constrained problem

$$
\min_\mathbf{x}  f(\mathbf{x}) \\
\text{subject to }  A\mathbf{x}=\mathbf{b}
$$

we apply a second-order Taylor approximation at the point $\mathbf{x}$, to obtain

$$
\min_\mathbf{x} \hat{f}(\mathbf{x} +\mathbf{v}) = f(\mathbf{x}) +\nabla f(\mathbf{x})^\top \mathbf{v}+ \frac{1}{2}\mathbf{v}^\top \nabla^2 f(\mathbf{x}) \mathbf{v} \\
\text{subject to } A(\mathbf{x}+\mathbf{v})=\mathbf{b}\,.
$$

Based on the solution of quadratic convex problems with linear constraints, the Newton $\Delta \mathbf{x}_{nt}$ step is characterized by

$$
\begin{bmatrix}
 \nabla^2 f(\mathbf{x})&  A^\top \\
A & 0 \\
     \end{bmatrix}
     \begin{bmatrix}
\Delta \mathbf{x}_{nt}\\
\mathbf{w}
     \end{bmatrix}
     =
     -\begin{bmatrix}
\nabla f(\mathbf{x}) \\
A\mathbf{x}-\mathbf{b}
     \end{bmatrix}
$$

- If the starting point $\mathbf{x}^{(0)}$ is chosen such that $A\mathbf{x}^{(0)}=\mathbf{b}$, the residual term vanishes and steps will remain in the feasible region. This is the **feasible start Newton method**.
- If we choose an arbitrary $\mathbf{x}^{(0)}\in$ **dom** $f$, not satisfying the constraints, this is the **infeasible start Newton method**. It will usually converge rapidly to the feasible region (check the final solution!).

Note that when we start at a feasible point, the residual vector $-(A\mathbf{x}-\mathbf{b})$ vanishes and the path will always remain in a feasible region. Otherwise we will converge to it.

In this chapter, we will use a fixed step size. For Newton's method this usually leads to only a few extra iterations compared to an adaptive step size.

>**input** starting point $\mathbf{x}\in$ **dom** $f$ (with $A\mathbf{x}=\mathbf{b}$ if using the feasible method), tolerance $\epsilon>0$.
>
>**repeat**
>
>>    1. Compute the Newton step $\Delta \mathbf{x}_{nt}$ and decrement $\lambda(\mathbf{x})$.
>>    2. *Stopping criterion*. **break** if $\lambda^2/2\leq \epsilon$.
>>    3. *Choose step size $t$*: either by line search or fixed $t$.
>>    4. *Update*. $\mathbf{x}:=\mathbf{x}+t \Delta \mathbf{x}_{nt}$.
>
>**output** $\mathbf{x}$

Again, the convergence can be monitored using the Newton decrement:

$$
\lambda^2(\mathbf{x}) = - \Delta \mathbf{x}_{nt}^\top \nabla f(\mathbf{x})\,.
$$

The algorithm terminates when

$$
\frac{\lambda(\mathbf{x})^2}{2} < \epsilon\,.
$$

The Newton decrement also indicates whether we are in or close to the feasible region.
**Assignment 2**

1. Complete the code for the linearly constrained Newton method.
2. Use this code to find the minimum of the non-quadratic toy problem, defined above (compare a feasible and infeasible start).<jupyter_output><empty_output><jupyter_text>from constrained import linear_constrained_newton<jupyter_code>def linear_constrained_newton(f, x0, grad_f,
              hess_f, A, b, stepsize=0.25, epsilon=1e-3,
              trace=False):
    '''
    Newton's method for minimizing functions with linear constraints.

    Inputs:
        - f: function to be minimized
        - x0: starting point (does not have to be feasible)
        - grad_f: gradient of the function to be minimized
        - hess_f: hessian matrix of the function to be minimized
        - A, b: linear constraints
        - stepsize: step size for each Newton step (fixed)
        - epsilon: parameter to determine if the algorithm is converged
        - trace: (bool) store the path that is followed?

    Outputs:
        - xstar: the found minimum
        - x_steps: path in the domain that is followed (if trace=True)
        - f_steps: image of x_steps (if trace=True)
    '''
    assert stepsize < 1 and stepsize > 0
    x = x0  # initial value
    p, n = A.shape
    if trace: x_steps = [x.copy()]
    if trace: f_steps = [f(x0)]
    while True:
        ddfx = hess_f(x)
        dfx = grad_f(x)
        # calculate residual
        Dx, _ = solve_constrained_quadratic_problem(... # complete!
        newton_decrement = ...
        if ...  # stopping criterion
            break  # converged
        # perform step
        if trace: x_steps.append(x.copy())
        if trace: f_steps.append(f(x))
    if trace: return x, x_steps, f_steps    
    else: return x
# to complete

A = ... #  1 x 2 array
b = ... #  1 x 1 array
x0 = np.array([[0.6], [-0.2]])  # feasible
#x0 = np.array([[0.6], [0.2]])  # unfeasible
xstar, x_steps, f_steps = linear_constrained_newton(..., trace=True)  # complete
print('Minimizer:')
print(xstar)

fig, ax = plt.subplots()
plot_contour(nonquadratic, (-2, 2), (-1, 1), ax, [1, 3], 0)
add_path(ax=ax, x_steps=x_steps)
ax.scatter(xstar[0,0], xstar[1,0], 50, 'r', label='minimum')
ax.grid()
ax.legend(loc=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15452636233005415, 'fraction_numerical': 0.017795954460042004, 'mean_word_length': 3.807651434643996, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '308311', 'license': 'MIT', 'n_tokens_mistral': 3089, 'n_tokens_neox': 2905, 'n_words': 920, 'path': 'Chapters/03.Constrained/Chapter3.ipynb', 'repo_name': 'MichielStock/SelectedTopicsOptimization', 'size': '9069'}