LP_weeds

LP.py: solve_lp

Source: src/LP.py · User guide: Linear Programs · Theory: Scenario Theory

solve_lp builds a scenario linear program in CVXPY, solves it, and then works out the complexity $k$ (the size of a support list) and a degeneracy flag. It does not compute the risk. To get the risk, pass the returned complexity and N to quantify_risk.

This page is the developer reference. For a walkthrough of the four formulations, see Linear Programs.

---

Problem solved

$$\begin{aligned} \min_{x \in \mathbb{R}^{n},\ \zeta \in \mathbb{R}^{N}} \quad & 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 $\bar{x}$ is x_ref. All inequalities hold element by element. Each scenario contributes $m_s$ rows, and all of them share the same scalar slack $\zeta_i$. When slack is disabled, $\zeta \equiv 0$. The term $\tau \lVert x - \bar{x} \rVert_p$ is added only when $\tau > 0$.

The four formulations used in the web interface map onto the arguments like this:

Formulation	tau	rho	include_slack
Robust	0	0	False (or None)
Regularization	> 0	0	False (or None)
Relaxation	0	> 0	True (or None)
Relaxation + regularization	> 0	> 0	True (or None)

---

Signature

from src.LP import solve_lp

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

Parameters

Name	Type / shape	Default	Meaning
deltas	numpy.ndarray, shape (N, q) or (N,)	required	The scenarios $\delta_1,\dots,\delta_N$, one per row. N = len(deltas). Each row deltas[i] goes unchanged to A_d and b_d. If deltas is 1-D, each call gets a scalar. An empty deltas (or None) raises ValueError.
A_d	callable, δ → ndarray (m_s, n)	required	Returns $A(\delta)$. It must return a 2-D array. The number of decision variables is taken from n = 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 (m_s, 1), so a 1-D array gives the same constraints.
G	numpy.ndarray, (p, n)	required	Hard-constraint matrix. Pass np.array([]) if there are no hard constraints.
h	numpy.ndarray, (p, 1) or (p,)	required	Hard-constraint offset, reshaped to a column. If either G or h has size == 0, no hard constraint is added.
c	numpy.ndarray, (n, 1) or (n,)	required	Linear cost vector.
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 first. One entry is copied to every coordinate. n entries are reshaped to (n, 1). Any other size raises ValueError.
rho	float	0.0	Slack penalty $\rho$. When include_slack=None, it also decides whether slack variables exist.
norm_type	any number p >= 1, "inf" / np.inf, "fro"	2	Order of the vector norm of x - x_ref (the column is flattened first). "fro" is the 2-norm. Used only when tau > 0.
solver	str or None	None	CVXPY solver name. It is used for the main solve and for every re-solve inside get_support. None lets CVXPY choose. The web interface sends the solver picked in its dropdown, and falls back to "CLARABEL" if none is sent. See Solver compatibility.
include_slack	bool or None	None	True always adds $\zeta \ge 0$, even when rho == 0. The problem is then normally unbounded, and ValueError is raised. False never adds slack: $\zeta \equiv 0$ and the $\rho$ term is the constant 0. None adds slack exactly when rho != 0.

Returns

solve_lp returns a 7-tuple:

Name	Type / shape	Meaning
x_out	numpy.ndarray, (n, 1)	Optimal $x$, read straight after the main solve.
zeta_out	numpy.ndarray, (N,)	Optimal slacks when slack is enabled, otherwise np.zeros(N). A positive entry marks a violated scenario.
cost_out	float	prob.value, the full objective, including the $\tau$ and $\rho$ terms. The web interface shows $c^\top x$ instead.
N	int	Number of scenarios, len(deltas).
complexity	int	$k$, the number of scenarios in the support list found by get_support. It counts scenarios, not rows.
constraints	list of N cvxpy.constraints.Inequality	The scenario constraints, one per scenario, in the same order as deltas. Hard constraints are not included. See the caveat below.
degeneracy	bool	True when get_support could not certify that the support list is the unique one. The lower risk bound is then not certified (see Risk.py).

Note. The support list itself is computed but not returned; only its cardinality is. To find out which scenarios are of support, build the CVXPY problem yourself and call get_support directly (example in Miscellaneous.py).

Caveat. When solve_lp returns, the CVXPY variables behind constraints hold the solution of the last re-solve done inside get_support, not x_out. Dual values of constraints used in those re-solves are overwritten too. So constraints[i].violation() and constraints[i].dual_value do not describe the returned optimum. To check a constraint at the optimum, compute A_d(deltas[i]) @ x_out + b_d(deltas[i]). In the example below, x_out = [-0.198, 2.200], but afterwards the variable held [2.002, 0.0]. After the relaxed run, 10 of the 200 constraints had dual_value equal to None.

Exceptions

Exception	When
ValueError	deltas is empty or None: The input `deltas` must have at least one row.
ValueError	x_ref has neither 1 nor n entries.
ValueError	The solve status is not "optimal" or "optimal_inaccurate" (for example "infeasible", "unbounded" or "user_limit"). The message gives the status, the objective value and x.
cvxpy.error.SolverError	The solver is not installed, or it cannot handle the cones in the problem (see Solver compatibility).
ValueError	Raised by get_support when no valid support list can be recovered: "Error calculating support list! Perhaps try another solver?"

---

How it works

1. Dimensions. N = check_scenarios(deltas) (that is, len(deltas), after checking that deltas is not empty) and n = A_d(deltas[0]).shape[1].

2. Variables. x = cp.Variable((n, 1)). If slack is enabled (see include_slack), zeta = cp.Variable(N, nonneg=True). Otherwise zeta = np.zeros(N).

3. Scenario constraints. For each i, in a plain Python loop:

   b_i = np.asarray(b_d(deltas[i]), dtype=float).reshape(-1, 1)   # b(δ) as a column
   constraints.append(A_d(deltas[i]) @ x + b_i <= zeta[i])

   Each scenario gives one CVXPY Inequality with $m_s$ rows, and the scalar zeta[i] is broadcast across those rows. This is why the complexity counts scenarios, not rows. A_d and b_d are each called once per scenario. A_d is also called on deltas[0] to find the dimensions.

4. Hard constraints. non_risk_constraints = [G @ x + h <= 0] (with G made 2-D and h reshaped to a column), or [] if G or h is empty. They live in a separate list, so every re-solve enforces them but they are never support-list candidates and never count towards $k$.

5. Regularization centre. x_ref is flattened and reshaped to (n, 1), either by copying a single value or by reshaping n values.

6. Objective.

   objective = cp.Minimize(c.T @ x + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))

   regularization_term returns tau * cp.norm(x - x_ref, p), with x - x_ref flattened to a vector and "fro" mapped to 2. When tau == 0 it returns 0, so the norm is not in the model and a robust or relaxed LP is a pure LP.

7. Solve. prob = cp.Problem(objective, constraints + non_risk_constraints) and then prob.solve(solver=solver), with the solver's default settings (no options are passed). Any status other than optimal or optimal_inaccurate raises ValueError.

8. Support list and degeneracy. The function calls get_support(constraints, non_risk_constraints, prob, objective, solver=solver) with its defaults threshold=1e-8 and sol_tol=1e-6. In short:

   • Screen. Keep each scenario whose dual has $\max\lvert\lambda\rvert > 10^{-8}$. Such a scenario is active or violated. Under relaxation, the duals of a violated scenario sum to $\rho$ (stationarity in $\zeta_i$). This was checked numerically on the example problem with rho=0.1: every violated scenario had dual sum 0.1.
   • Record spurious screens. Note which screened constraints are strictly slack at the optimum. These are solver noise.
   • Validate. Re-solve with only the screened set plus the hard constraints. If the optimal value is reproduced (relative tolerance $10^{-6}$), greedily prune the set. degeneracy is set to True only if a removed constraint was active or violated.
   • Fallback. Otherwise, prune the full list of N scenarios and set degeneracy = True as a precaution.

   The full algorithm is on Miscellaneous.py → get_support.

9. Return. zeta_out and cost_out are read, and the 7-tuple is returned. The support list is thrown away.

Cost of the support computation

Each re-solve calls Problem.solve on a reduced problem with the same objective.

Path	Number of re-solves	Size of each re-solve
Normal (the screen reproduces the optimum)	$1 + s$, where $s$ is the number of screened scenarios	at most $s$ scenario constraints plus the hard constraints
Fallback (the screen fails)	$N + 2$	up to $N - 1$ scenario constraints

The fallback path usually dominates the runtime for large N. It also prints a message (see Misc_weeds).

Parallelism

There is none. The scenario loop, the dual screen, the pruning loop and all re-solves run one after another in a single Python thread. Only the solver's own internal threading, if it has any, is parallel.

Numerical tolerances

Where	Value	Purpose
Accepted solve status	"optimal", "optimal_inaccurate"	Anything else raises ValueError.
Dual screen (get_support, threshold)	1e-8 on $\max \lvert \text{dual} \rvert$	Picks the candidate support scenarios.
Value reproduction (test_support)	np.isclose(ref, new, rtol=1e-6, atol=1e-9)	Decides whether removing a scenario changes the optimum.
Strictly slack test (_is_slack)	margin $> 10^{-6} + 10^{-6}\cdot\text{scale}$	Separates spurious screens from genuinely active constraints.
Solver tolerances	solver defaults	No solver options are passed.

Solver compatibility

The cones that reach the solver depend on tau and norm_type:

Case	Cones produced	Solvers that can accept it (CVXPY 1.7.2 solver metadata)
tau == 0 (robust or relaxed LP)	linear only	any LP solver, e.g. SCIPY, GLPK, HIGHS, as well as the conic solvers
tau > 0, norm_type 1 or "inf"	linear only	as above
tau > 0, norm_type 2 or "fro"	linear + second-order cone	conic solvers such as CLARABEL, SCS, ECOS, CVXOPT, MOSEK
tau > 0, any other p	linear + power cone	CLARABEL, SCS, MOSEK

Tested with CVXPY 1.7.2 on the example below:

• SCIPY solves the robust LP ($k = 2$) and the relaxed LP ($k = 10$), and with tau > 0 it works for norm_type=1 but fails for norm_type=2 (SolverError).
• OSQP solves the robust LP ($k = 2$), but on the relaxed LP its duals are too inaccurate for the screen: the fallback path runs and reports degeneracy=True with $k = 164$. This is why the web interface offers OSQP for QPs only.
• SCS solves the robust LP but returns $k = 31$ with degeneracy=True, where the true $k$ is 2.

An interior-point solver such as MOSEK (recommended) or CLARABEL is the safe choice when you need the complexity.

---

Pitfalls

• Strings are not accepted directly. A_d and b_d must be Python callables. To build them from strings such as "delta[0]*delta[2]", use generate_matrix_function.
• cost_out includes the penalty terms. Compute float(c.ravel() @ x_out.ravel()) if you only want $c^\top x$.

How the web interface calls it

The web interface's solve code, solve_form in src/solve_form.py, calls

solve_lp(scenarios, A_d, b_d, G, h_vec, c, tau, theta_bar, rho, p, run_solver, include_slack=include_slack)

once for every $(\tau, \rho)$ pair in the comma-separated sweep fields:

• In symbolic mode, A_d and b_d come from generate_matrix_function. In numeric mode they come from generate_numeric_A_b.
• The selected option decides which parameters apply: rho is set to 0 for robust and regularization, and tau to 0 for robust and relaxation, whatever the hidden fields contain.
• include_slack is True for the options relaxation and regularization-relaxation, and False for robust and regularization.
• p is parsed from the form: "inf" becomes np.inf, "fro" stays a string, and anything else goes through float(...) and must be at least 1.
• The cost shown is $c^\top x$, not cost_out.
• If a MOSEK licence is uploaded, the same solve_form runs in a subprocess with solver="MOSEK" (see mosek_solve.py).

---

Example

This example finds the smallest interval $[x_1 - x_2,\ x_1 + x_2]$ that contains 200 samples of a scalar $\delta$. Run it from the repository root (python example.py), so that src can be imported.

import numpy as np
from src.LP import solve_lp
from src.Risk import quantify_risk

# Smallest interval [x1 - x2, x1 + x2] containing N samples of a scalar delta:
#   min x2   s.t.   delta_i - x1 - x2 <= 0,   x1 - x2 - delta_i <= 0
rng = np.random.default_rng(0)
deltas = rng.standard_normal((200, 1))              # N = 200 scenarios, q = 1

A_d = lambda d: np.array([[-1.0, -1.0], [1.0, -1.0]])   # (m_s, n) = (2, 2)
b_d = lambda d: np.array([[d[0]], [-d[0]]])             # (m_s, 1) = (2, 1)
c = np.array([[0.0], [1.0]])                            # (n, 1)
G = h = np.array([])                                    # no hard constraints

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

# Relaxation: slack zeta_i >= 0 per scenario, penalised by rho
x, zeta, cost, N, k, cons, degenerate = solve_lp(deltas, A_d, b_d, G, h, c, rho=0.1, solver="CLARABEL")
print("relaxed x =", x.ravel(), " k =", k, " violated =", int((zeta > 1e-7).sum()), " degenerate =", degenerate)
epsL, epsU = quantify_risk(k, N, 1e-6)
print(f"risk in [{epsL:.4f}, {epsU:.4f}] with confidence 1 - 1e-6")

# Degenerate instance: duplicate the largest sample, so two support lists exist
deltas_dup = np.vstack([deltas, deltas[[deltas.argmax()]]])
x, zeta, cost, N, k, cons, degenerate = solve_lp(deltas_dup, A_d, b_d, G, h, c, solver="CLARABEL")
print("duplicated max: N =", N, " k =", k, " degenerate =", degenerate)

Output (CVXPY 1.7.2, CLARABEL):

x = [-0.19792014  2.20031272]  cost = 2.200313  N = 200  k = 2  degenerate = False
support scenarios: argmin = 159  argmax = 79
risk in [0.0000, 0.1101] with confidence 1 - 1e-6
relaxed x = [-0.10880633  1.86414658]  k = 10  violated = 10  degenerate = False
risk in [0.0000, 0.1852] with confidence 1 - 1e-6
PRUNE removed 1 active/violated support scenario(s) (3 screened -> 2 irreducible); flagging degeneracy (lower risk bound not certified).
duplicated max: N = 201  k = 2  degenerate = True

What the output shows:

• Robust run. The two extreme samples are the support list, so $k = 2$. The interval's centre and half-width are exactly $(\max\delta + \min\delta)/2$ and $(\max\delta - \min\delta)/2$.
• Relaxed run. The 10 violated samples form the support list, so $k = 10$, which equals the number of positive slacks.
• Duplicated maximum. The two copies make the support list non-unique. Pruning removes one copy, so degeneracy=True, and only the upper risk bound should be used.

See also

• QP.py and SDP.py: the same pipeline with a quadratic objective or LMI constraints.
• Miscellaneous.py: get_support, test_support, the solver table and the parsing helpers.
• Risk.py: turning (complexity, N, beta) into risk bounds.