Skip to content

QP_weeds

Jump to bottom
Kiguli edited this page Oct 6, 2026 · 3 revisions

QP.py: solve_qp

Source: src/QP.py · User guide: Quadratic Programs · Theory: Scenario Theory

solve_qp is the quadratic-objective version of solve_lp. The constraints, the slack handling, the support-list computation and the return values are the same. What changes is the objective, which gains $\tfrac12 x^\top Q x$, and a validation step for Q. This page covers the whole interface and goes into detail only where the QP differs from the LP. The shared mechanics are on LP.py.

Problem solved

$$\begin{aligned} \min_{x \in \mathbb{R}^{n},\ \zeta \in \mathbb{R}^{N}} \quad & \tfrac{1}{2}\, x^\top Q\, x \;+\; c^\top x \;+\; \tau \lVert x - \bar{x} \rVert_p \;+\; \rho \sum_{i=1}^{N} \zeta_i \\\ \text{s.t.} \quad & A(\delta_i)\, x + b(\delta_i) \;\le\; \zeta_i \mathbf{1}, \qquad i = 1, \dots, N, \\\ & G x + h \;\le\; 0, \\\ & \zeta_i \;\ge\; 0, \qquad i = 1, \dots, N, \end{aligned}$$

where $Q = Q^\top \succeq 0$ and $\bar{x}$ is x_ref. The robust, regularization and relaxation formulations come from setting tau, rho and include_slack exactly as for the LP (see the mapping table).

Signature

from src.QP import solve_qp

x_out, zeta_out, cost_out, N, complexity, constraints, degeneracy = solve_qp(
    deltas, A_d, b_d, G, h, c, Q,
    tau=0.0, x_ref=np.array([0.0]), rho=0.0, norm_type=2,
    solver=None, include_slack=None,
)

Q is the seventh positional argument, after c. Everything before and after it matches solve_lp.

Parameters

Name Type / shape Default Meaning
deltas numpy.ndarray, (N, q) or (N,) required The scenarios, one per row. Each row goes unchanged to A_d and b_d. An empty deltas raises ValueError.
A_d callable, δ → ndarray (m_s, n) required Returns $A(\delta)$. n is taken from A_d(deltas[0]).shape[1].
b_d callable, δ → ndarray (m_s, 1) or (m_s,) required Returns $b(\delta)$; the result is reshaped to a column.
G numpy.ndarray, (p, n) required Hard-constraint matrix. Pass np.array([]) if there are none.
h numpy.ndarray, (p, 1) or (p,) required Hard-constraint offset, reshaped to a column. If G or h has size 0, no hard constraint is added.
c numpy.ndarray, (n, 1) or (n,) required Linear cost vector.
Q numpy.ndarray, (n, n) required Quadratic cost matrix. It must be symmetric positive semidefinite, up to a small tolerance (see Validation of Q). np.zeros((n, n)) is allowed and reduces the problem to an LP.
tau float 0.0 Regularization weight $\tau$.
x_ref scalar or array with 1 or n entries np.array([0.0]) Regularization centre $\bar{x}$. It is flattened, then either copied to every coordinate (1 entry) or reshaped to (n, 1) (n entries). Any other size raises ValueError.
rho float 0.0 Slack penalty $\rho$. With include_slack=None, slack exists exactly when rho != 0.
norm_type any number p >= 1, "inf" / np.inf, "fro" 2 Order of the vector norm of x - x_ref, as in the LP. Used only when tau > 0.
solver str or None None CVXPY solver for the main solve and for all re-solves. See Solver compatibility.
include_slack bool or None None True forces slack, False forbids it, and None means "slack if rho != 0". Same semantics as solve_lp.

Returns

Name Type / shape Meaning
x_out numpy.ndarray, (n, 1) Optimal $x$.
zeta_out numpy.ndarray, (N,) Optimal slacks, or np.zeros(N) without slack.
cost_out float prob.value, the full objective including $\tfrac12 x^\top Q x$, the $\tau$ term and the $\rho$ term. The web interface shows $c^\top x + \tfrac12 x^\top Q x$ instead.
N int len(deltas).
complexity int $k$, the number of scenarios in the support list.
constraints list of N cvxpy.constraints.Inequality Scenario constraints in the order of deltas. Hard constraints are excluded. Their variable values are overwritten by later re-solves, as described in the LP caveat.
degeneracy bool True when the lower risk bound is not certified.

Exceptions

Exception When
AssertionError Q is not symmetric within the tolerance: "Q must be symmetric". This check runs first.
AssertionError Q has an eigenvalue below the negative tolerance: "Q must be positive semi-definite (eigenvalues: ...)".
ValueError deltas is empty or None.
ValueError x_ref has neither 1 nor n entries.
Exception (from CVXPY) Q is not n × n: "Invalid dimensions for arguments to quad_form."
ValueError The status is not optimal or optimal_inaccurate: "QP did not solve to optimality. Status: ..."
cvxpy.error.SolverError The solver is missing or cannot handle the problem's cones.
ValueError get_support cannot recover a support list.

How it works

The pipeline is the same as solve_lp, with two additions.

Validation of Q

Before anything else, the function runs Q = check_psd(Q) (in src/Miscellaneous.py):

Q = np.asarray(Q, dtype=float)
scale = max(1.0, float(np.abs(Q).max()))
assert np.allclose(Q, Q.T, rtol=0, atol=tol * scale), "Q must be symmetric"       # tol = 1e-9
Q = (Q + Q.T) / 2
eigvals = np.linalg.eigvalsh(Q)
assert eigvals.min() >= -tol * scale, f"Q must be positive semi-definite (eigenvalues: {eigvals})"
  • Both tests use a relative tolerance of $10^{-9}$ times the largest absolute entry of Q (or $10^{-9}$ when all entries are below 1).
  • Rank-deficient matrices pass. Floating-point rounding often turns a zero eigenvalue into a tiny negative number, e.g. about -6e-18 for M @ M.T with M of shape 3 × 2. In a test with 20 such random matrices, all 20 were accepted.
  • Symmetry up to rounding passes, and the symmetrised matrix $\tfrac12 (Q + Q^\top)$ is what the solver uses.
  • The eigenvalues come from np.linalg.eigvalsh, the symmetric solver, applied after symmetrisation.
  • Q is not projected onto the PSD cone. A matrix whose most negative eigenvalue lies within the tolerance is passed on unchanged, and CVXPY is told to treat it as PSD (below). For matrices that are PSD up to rounding this is harmless.

Objective

objective = cp.Minimize((1/2) * cp.quad_form(x, cp.psd_wrap(Q))
                        + c.T @ x + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))

cp.psd_wrap(Q) tells CVXPY to trust that Q is PSD, so CVXPY does not repeat its own eigenvalue check (check_psd already did it). As with the LP, the norm term is left out when tau == 0.

Everything else

Variables, scenario constraints (with b(δ) reshaped to a column), hard constraints, the reshaping of x_ref, the status check and the call to get_support(constraints, non_risk_constraints, prob, objective, solver=solver) are line for line the same as in LP.py. This covers the dual screen at 1e-8, the value-based validation at relative tolerance 1e-6, greedy pruning, the degeneracy flag and the fallback path. The re-solves inside get_support reuse the same objective, so they are QPs too. All loops run in sequence, with no parallelism.

Solver compatibility

Case What CVXPY sends to the solver Solvers that can accept it (CVXPY 1.7.2 solver metadata)
tau == 0, or norm_type 1 or "inf" quadratic objective plus linear constraints conic solvers (CLARABEL, SCS, ECOS, CVXOPT, MOSEK) and QP solvers (OSQP, PIQP, PROXQP, DAQP, …)
tau > 0, norm_type 2 or "fro" quadratic objective plus a second-order cone conic solvers only
tau > 0, any other p quadratic objective plus a power cone CLARABEL, SCS, MOSEK

Tested on the example below (CVXPY 1.7.2):

  • CLARABEL gives $k = 1$, non-degenerate, without regularization and with tau=0.5 for norm_type 1, 2 and 3.
  • SCS gives $k = 1$ without regularization. With tau=0.5 its solutions are not accurate enough for the support-list screen: the fallback runs and reports degeneracy=True with $k = 77$ (norm_type=1) or $k = 274$ (norm_type=2).
  • OSQP gives $k = 1$ without regularization and with norm_type=1, and rejects norm_type=2 with tau > 0 ("Problem could not be reduced to a QP").
  • SCIPY is an LP solver and cannot solve QPs (SolverError); get_solvers offers it for LPs only.

Example

This example finds the minimum-norm point that satisfies $\delta_i^\top x \ge 1$ for 300 sampled $\delta_i$. Run it from the repository root.

import numpy as np
from src.QP import solve_qp
from src.Risk import quantify_risk

# Minimum-norm point that satisfies delta_i^T x >= 1 for every sampled delta_i:
#   min 1/2 x^T Q x   s.t.   1 - delta_i^T x <= 0
rng = np.random.default_rng(1)
deltas = rng.uniform(0.5, 1.5, size=(300, 2))      # N = 300 scenarios, q = 2

A_d = lambda d: -np.array([[d[0], d[1]]])          # (1, 2)
b_d = lambda d: np.array([[1.0]])                  # (1, 1)
G = h = np.array([])
c = np.zeros((2, 1))
Q = np.eye(2)                                      # symmetric PSD

x, zeta, cost, N, k, cons, degenerate = solve_qp(deltas, A_d, b_d, G, h, c, Q, solver="CLARABEL")
print("x =", x.ravel(), " cost =", round(cost, 6), " N =", N, " k =", k, " degenerate =", degenerate)
epsL, epsU = quantify_risk(k, N, 1e-6)
print(f"risk in [{epsL:.4f}, {epsU:.4f}] with confidence 1 - 1e-6")

# Same problem with regularisation towards x_ref = (1, 1) using the 1-norm
x, zeta, cost, N, k, cons, degenerate = solve_qp(deltas, A_d, b_d, G, h, c, Q,
                                                 tau=0.5, x_ref=np.array([1.0, 1.0]), norm_type=1,
                                                 solver="CLARABEL")
print("regularised x =", x.ravel(), " cost =", round(cost, 6), " k =", k, " degenerate =", degenerate)

Output (CVXPY 1.7.2, CLARABEL):

x = [0.85388612 0.95576728]  cost = 0.821306  N = 300  k = 1  degenerate = False
risk in [0.0000, 0.0662] with confidence 1 - 1e-6
regularised x = [0.88352631 0.92928663]  cost = 0.91569  k = 1  degenerate = False

Here a single half-space, the one closest to the origin, defines the optimum, so $k = 1$. In the second run cost includes $0.5 \lVert x - (1,1) \rVert_1$.

See also

Using The Tool

Home

Starting Up

Convex Optimization

Graphical Results

Getting into the Weeds

Scenario Theory

Implementation Details

Clone this wiki locally