Quadratic program¶
A quadratic program is an optimization problem with a quadratic objective and affine equality and inequality constraints.
A common standard form is the following:
Here \(P \in \mathcal{S}^{n}_+\), \(q \in \mathcal{R}^n\), \(G \in \mathcal{R}^{m \times n}\), \(h \in \mathcal{R}^m\), \(A \in \mathcal{R}^{p \times n}\), and \(b \in \mathcal{R}^p\) are problem data and \(x \in \mathcal{R}^{n}\) is the optimization variable. The inequality constraint \(Gx \leq h\) is elementwise.
A simple example of a quadratic program arises in finance. Suppose we have \(n\) different stocks, an estimate \(r \in \mathcal{R}^n\) of the expected return on each stock, and an estimate \(\Sigma \in \mathcal{S}^{n}_+\) of the covariance of the returns. Then we solve the optimization problem:
to find a nonnegative portfolio allocation \(x \in \mathcal{R}^n_+\) that optimally balances expected return and variance of return.
When we solve a quadratic program, in addition to a solution \(x^\star\), we obtain a dual solution \(\lambda^\star\) corresponding to the inequality constraints. A positive entry \(\lambda^\star_i\) indicates that the constraint \(g_i^Tx \leq h_i\) holds with equality for \(x^\star\) and suggests that changing \(h_i\) would change the optimal value.
Quadratic Program Example¶
In the following code, we solve a quadratic program with CVXPY.
[1]:
# Import packages.
import cvxpy as cp
import numpy as np
rng = np.random.default_rng(1)
# generate a random feasible quadratic program.
m = 15
n = 10
p = 5
P = rng.standard_normal((n, n))
P = P.T @ P
q = rng.standard_normal(n)
G = rng.standard_normal((m, n))
# generate a feasible point first, then construct constraints around it
x_feasible = rng.standard_normal(n)
h = G @ x_feasible + rng.exponential(1, m)
A = rng.standard_normal((p, n))
b = A @ x_feasible
# define and solve the CVXPY problem.
x = cp.Variable(n)
prob = cp.Problem(
cp.Minimize((1 / 2) * cp.quad_form(x, P) + q.T @ x), [G @ x <= h, A @ x == b]
)
out = prob.solve(solver=cp.SCS, verbose=True)
# Print result.
print("The solver status is", prob.status)
print("\nThe optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("A dual solution corresponding to the inequality constraints is")
print(prob.constraints[0].dual_value)
(CVXPY) Aug 25 12:26:34 PM: Your problem has 10 variables, 20 constraints, and 0 parameters.
(CVXPY) Aug 25 12:26:34 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Aug 25 12:26:34 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Aug 25 12:26:34 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Aug 25 12:26:34 PM: Your problem is compiled with the CPP canonicalization backend.
(CVXPY) Aug 25 12:26:34 PM: Compiling problem (target solver=SCS).
(CVXPY) Aug 25 12:26:34 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS
(CVXPY) Aug 25 12:26:34 PM: Applying reduction Dcp2Cone
(CVXPY) Aug 25 12:26:34 PM: Applying reduction CvxAttr2Constr
(CVXPY) Aug 25 12:26:34 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Aug 25 12:26:34 PM: Applying reduction SCS
(CVXPY) Aug 25 12:26:34 PM: Finished problem compilation (took 2.980e-02 seconds).
(CVXPY) Aug 25 12:26:34 PM: Invoking solver SCS to obtain a solution.
(CVXPY) Aug 25 12:26:34 PM: Problem status: optimal
(CVXPY) Aug 25 12:26:34 PM: Optimal value: 1.678e+01
(CVXPY) Aug 25 12:26:34 PM: Compilation took 2.980e-02 seconds
(CVXPY) Aug 25 12:26:34 PM: Solver (including time spent in interface) took 7.689e-03 seconds
===============================================================================
CVXPY
v1.7.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
------------------------------------------------------------------
SCS v3.2.8 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 10, constraints m: 20
cones: z: primal zero / dual free vars: 5
l: linear vars: 15
settings: eps_abs: 1.0e-05, eps_rel: 1.0e-05, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 200, nnz(P): 55
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 6.47e+00 5.26e+00 3.10e+01 1.92e+01 1.00e-01 4.56e-03
100| 5.06e-05 4.41e-06 2.92e-04 1.68e+01 1.00e-01 4.74e-03
------------------------------------------------------------------
status: solved
timings: total: 4.76e-03s = setup: 3.29e-03s + solve: 1.48e-03s
lin-sys: 7.25e-05s, cones: 8.02e-06s, accel: 4.69e-06s
------------------------------------------------------------------
objective = 16.784126
------------------------------------------------------------------
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
The solver status is optimal
The optimal value is 16.78397968501814
A solution x is
[-0.72413284 0.46981846 -1.3871876 -1.09984742 -0.04464528 -1.02325643
0.35526475 -0.15215945 1.02594271 -0.99970021]
A dual solution corresponding to the inequality constraints is
[ 0.00000000e+00 3.35442583e-17 0.00000000e+00 0.00000000e+00
-5.11883259e-18 6.05443281e+00 1.64824687e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 0.00000000e+00]
[2]:
assert prob.status == cp.OPTIMAL, "Solver result not optimal"
assert prob.value is not None, "Problem value is None"