Worksheet Chapter 12, Part 3

Worksheet Chapter 12, Part 3#

Convex Quadratic Optimization#

Notebook Info#

Topics#

Strictly Convex Quadratic Optimization
Convex Quadratic Optimization

Learning goals#

Use cvxpy to code both the primal and dual formulations of a strictly convex quadratic optimization problem.
Obtain the primal and dual solutions from a single run.
Use cvxpy to solve a (non-strict) convex quadratic program.

import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp

%matplotlib inline
np.set_printoptions(precision=4, suppress=True)

1. Strictly Convex Quadratic Programming#

Problem Setup#

Consider the following strictly convex quadratic program:

\[\begin{split}\begin{aligned} \text{(P)} \quad & f^* = \min_{\mathbf{x}} \quad & \mathbf{x}^T Q \mathbf{x} + 2\mathbf{c}^T \mathbf{x} \\ & \text{such that} & A\mathbf{x} \leq \mathbf{b} \end{aligned}\end{split}\]

Data:

\[\begin{split}Q = \begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \mathbf{c} = \begin{pmatrix} 1 \\ 0 \\ -1 \end{pmatrix}, \quad A = \begin{pmatrix} 1 & 1 & 0 \\ 0 & 1 & 1 \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} 2 \\ 1 \end{pmatrix}\end{split}\]

Equivalently in scalar form:

\[\begin{split}\begin{aligned} \min \quad & 2x_1^2 + 2x_2^2 + x_3^2 + 2x_1 - 2x_3 \\ \text{such that} \quad & x_1 + x_2 \leq 2 \\ & x_2 + x_3 \leq 1 \end{aligned}\end{split}\]

# Set up the data for the problem
Q_prob = np.array([[2, 0, 0], 
                   [0, 2, 0], 
                   [0, 0, 1]], dtype=float)
c_prob = np.array([1, 0, -1], dtype=float)
A_prob = np.array([[1, 1, 0], 
                   [0, 1, 1]], dtype=float)
b_prob = np.array([2, 1], dtype=float)

❓❓❓ Question ❓❓❓: Check that this problem is strictly convex by ensuring that \(Q\) is positive definite.

# your code here

Solve the Primal Problem using CVXPY#

❓❓❓ Question ❓❓❓: Use CVXPY to solve the primal problem as written above. Store the optimal value in f_star_primal and the optimal solution in x_star.

Hints:

The objective is: x @ Q @ x + 2 * c @ x
The constraint is: A @ x <= b. This needs to get passed to cp.Problem as the only entry in a list, so save it as something like constraints = [A@x<=b].
Use cp.Problem(cp.Minimize(...), [...]) to formulate the problem.
Call .solve() on the problem to solve it
Access the optimal value with problem.value and the solution with variable.value

# Your code here

# x = 
# objective = 
# constraints = [...]
# primal_problem = 

# primal_problem.solve()

# f_star_primal = 
# x_star = 

Solve the Dual Problem#

From the lecture, recall that for a strictly convex QP with positive definite \(Q\), the dual problem is:

\[\begin{split}\begin{aligned} \text{(D)} \quad & q^* = \max_{\boldsymbol{\lambda}} \, q(\boldsymbol{\lambda}) \\ & \text{such that} \quad \boldsymbol{\lambda} \geq 0 \end{aligned}\end{split}\]

where the dual objective is

\[q(\boldsymbol{\lambda}) = -\boldsymbol{\lambda}^T A Q^{-1} A^T \boldsymbol{\lambda} - 2(AQ^{-1}\mathbf{c} + \mathbf{b})^T \boldsymbol{\lambda} - \mathbf{c}^T Q^{-1} \mathbf{c}.\]

❓❓❓ Question ❓❓❓: We’re going to break up the dual objective functions into pieces.

Compute \(Q^{-1}\) and store it in Q_inv.
Compute the three components of the dual objective:
- \(H = A Q^{-1} A^T\) (store as H_dual)
- \(g = AQ^{-1}\mathbf{c} + \mathbf{b}\) (store as g_dual)
- \(\mathbf{c}^T Q^{-1} \mathbf{c}\) (store as constant_term)

# Your code here

❓❓❓ Question ❓❓❓: We will now set up the objective function, which with our replacements above is

\[q(\boldsymbol{\lambda}) = -\boldsymbol{\lambda}^T H \boldsymbol{\lambda} - 2 g^T \boldsymbol{\lambda} - e,\]

where \(H =\) H_dual, \(g =\) g_dual, and \(e = \) constant_term.

Setup and solve the problem using cvxpy. To do this:

Set up the variables \(\lambda\) as lam. Remeber that Python reserves the word lambda so do not save the variables as this!
Formulate the dual objective. Be sure to use cp.quad_form(lam, H_dual) to represent \(\boldsymbol{\lambda}^T H \boldsymbol{\lambda}\) for the first term so that cvxpy recognizes it as a convex function.
Add the constraints \(\boldsymbol{\lambda} \geq 0\).
Create and solve the dual problem using cp.Problem(cp.Maximize(...), [...]). Save the problem as dual_problem so we can use it later.
Store the dual optimal value in q_star_dual and the optimal dual variables in lam_star

# your code here

Comparison#

Strong duality tells us that for a strictly convex problem, the primal and dual optimal values must be equal (\(f^* = q^*\)).

❓❓❓ Question ❓❓❓: Check that your solutions from above, f_star_primal and q_star_dual are the same.

# Your code here

It turns out that cvxpy was actually solving both the primal and dual for you in the background when you were solving the primal problem. From the primal problem, you saved your constraints as a list constraints = [...].

constraints

[Inequality(Expression(AFFINE, UNKNOWN, (2,)))]

These constraints save the dual value.

# Extract dual variables from primal solve
lam_primal = np.array([constraints[i].dual_value for i in range(len(constraints))])
print(f"λ from primal solve: {lam_primal}")

λ from primal solve: [[0. 0.]]

❓❓❓ Question ❓❓❓: Check that the \(\lambda\) found in the primal problem and the ones you calculated directly using the dual problem are the same.

# Your code or notes here

2. Convex Quadratic Programming#

Now consider a convex (not necessarily strictly convex) quadratic program in the same notation as the slides:

\[\begin{split}\begin{aligned} \text{(P2)}\quad & \min_{\mathbf{x}} \quad & \mathbf{x}^T Q\mathbf{x} + 2\mathbf{c}^T\mathbf{x} \\ & \text{such that} & A\mathbf{x} \leq \mathbf{b} \end{aligned}\end{split}\]

For this example, use

\[\begin{split}Q = \begin{pmatrix}2&0&0\\0&0&0\\0&0&1\end{pmatrix},\quad \mathbf{c} = \begin{pmatrix}-1\\1\\0\end{pmatrix},\end{split}\]

\[\begin{split}A = \begin{pmatrix}1&1&0\\-1&0&1\\0&1&1\\0&-1&0\end{pmatrix},\quad \mathbf{b}=\begin{pmatrix}2\\1\\2\\1\end{pmatrix}.\end{split}\]

Notice: \(Q\succeq 0\) (positive semidefinite), so this is convex, but not strictly convex.

# Data for convex QP (Q is PSD, not PD)
Q2_prob = np.array([[2, 0, 0],
                    [0, 0, 0],
                    [0, 0, 1]], dtype=float)
c2_prob = np.array([-1, 1, 0], dtype=float)
A2_prob = np.array([[1, 1, 0],
                    [-1, 0, 1],
                    [0, 1, 1],
                    [0, -1, 0]], dtype=float)
b2_prob = np.array([2, 1, 2, 1], dtype=float)

Solving original problem#

❓❓❓ Question ❓❓❓: Setup and solve this original convex QP using cvxpy.

\[\min_{\mathbf{x}}\ \mathbf{x}^TQ\mathbf{x}+2\mathbf{c}^T\mathbf{x}\quad\text{s.t. }A\mathbf{x}\le\mathbf{b}.\]

Save the value/optimizer as f_star_convex, x_star_convex. What is the dual solution?

#Your code here

Solving the modified original problem#

Ok, so you already know from above that the dual problem has been solved, but let’s try to set up the dual problem directly to see the results. To do this, first we need to find a matrix \(D\) such that \(D^\top D = Q\) where we’ve saved \(Q =\) Q2_prob. We can do this using spectral factorization.

eigvals2 = np.linalg.eigvals(Q2_prob)
print(f"Eigenvalues of Q2: {eigvals2}")
print(f"Q2 is PSD (convex): {np.all(eigvals2 >= -1e-10)}")
print(f"Q2 is PD (strictly convex): {np.all(eigvals2 > 1e-10)}")

# Build D so that Q = D^T D using spectral factorization
w, V = np.linalg.eigh(Q2_prob)
w = np.clip(w, 0.0, None)
D2 = np.diag(np.sqrt(w)) @ V.T

Eigenvalues of Q2: [2. 0. 1.]
Q2 is PSD (convex): True
Q2 is PD (strictly convex): False

Check ||D^T D - Q||_F: 4.44e-16

❓❓❓ Question ❓❓❓: Check that \(Q=D^T D\).

# Your code here

❓❓❓ Question ❓❓❓: Now, set up and solve the reformulated primal problem

\[\begin{split} \begin{aligned} & \text{min} & & \|\mathbf{z}\|^2 + 2 \mathbf{c}^{\top} \mathbf{x} \\ & \text{such that} & & A \mathbf{x} \leq \mathbf{b}, \\ & & & \mathbf{z} = D \mathbf{x}. \end{aligned} \end{split}\]

Save the results as f_star_reform, x_star_reform, z_star_reform. Did you get the same answer as the original problem version?

# Your code here

Solving the dual problem#

Ok, so finally, we can solve the dual problem for the reformulated version. ❓❓❓ Question ❓❓❓: Set up and solve the convex-QP dual using cvxpy.

\[\max_{\lambda,\mu}\ -\|\mu\|^2-2\mathbf{b}^T\lambda\]

\[\text{s.t. }\mathbf{c}+A^T\lambda-D^T\mu=0,\ \lambda\ge0,\ \mu\in\mathbb{R}^n.\]

Save the dual value as q_star_convex_dual.

# Your code here

Compare#

Now that you’ve solved all the different versions, run the following code. Do you get the same answer every time?

print("="*60)
print("CONVEX QP VIA REFORMULATION + DUAL")
print("="*60)

print(f"\nOriginal primal value f*:   {f_star_convex:.8f}")
print(f"Reformulated value f*:      {f_star_reform:.8f}")
print(f"Dual value q*:              {q_star_convex_dual:.8f}")

print(f"\n|f_original - f_reform|:    {abs(f_star_convex - f_star_reform):.2e}")
print(f"|f_reform - q_dual|:        {abs(f_star_reform - q_star_convex_dual):.2e}")

print(f"\nOriginal x*:                {x_star_convex}")
print(f"Reformulated x*:            {x_star_reform}")
print(f"Reformulated z* (D x*):     {z_star_reform}")
print(f"D x* from reformulated x*:  {D2 @ x_star_reform}")

============================================================
CONVEX QP VIA REFORMULATION + DUAL
============================================================

Original primal value f*:   -2.50000000
Reformulated value f*:      -2.50000000
Dual value q*:              -2.50000000

|f_original - f_reform|:    0.00e+00
|f_reform - q_dual|:        0.00e+00

Original x*:                [ 0.5 -1.   0. ]
Reformulated x*:            [ 0.5 -1.   0. ]
Reformulated z* (D x*):     [0.     0.     0.7071]
D x* from reformulated x*:  [0.     0.     0.7071]