Scenario_Theory

Scenario Theory

This page covers the parts of scenario theory that Scen-Opt implements. It explains what problem is solved and what the reported quantities mean: the complexity $k$, the risk bounds $[\underline{\varepsilon}, \overline{\varepsilon}]$ and the confidence $1-\beta$. It also says when the lower bound can be trusted. For the code, see Risk.py and Miscellaneous.py. For the three program classes, see Linear Programs, Quadratic Programs and Semidefinite Programs.

The running example is the smallest enclosing interval (benchmarks/LP_half_width_2d). It takes $N = 100$ points $\delta_i \in [0,1]$ and finds the narrowest interval that contains them. The decision is $x = [o, w]^\top$, with centre $o$ and half-width $w$. It is also the example on First Example: Python Calls.

1. Scenarios

The uncertainty is a vector $\delta \in \Delta \subseteq \mathbb{R}^q$ with a probability distribution $\mathbb{P}$. Scen-Opt never needs $\mathbb{P}$. It only needs $N$ samples $\delta_1, \ldots, \delta_N$, called scenarios, drawn independently from that same $\mathbb{P}$ (i.i.d.).

In Scen-Opt, each scenario is one row of the scenarios file:

• In symbolic mode the row is $\delta_i$ itself. Its components are delta[0], ..., delta[q-1].
• In numeric mode the row holds the constraint data of that scenario, already evaluated. See First Example: Web Interface.

The i.i.d. assumption is the only statistical assumption, and it matters. The rows must be independent draws from the conditions the decision will face. If the data were filtered, hand-picked or correlated in time, the guarantees below do not apply.

2. The scenario program

All formulations in Scen-Opt are special cases of

$$\begin{aligned} \min_{x \in \mathbb{R}^d,\ \zeta_1, \ldots, \zeta_N \ge 0} \quad & c(x) + \tau \lVert x - \bar{x} \rVert_p + \rho \sum_{i=1}^{N} \zeta_i \\\ \text{subject to} \quad & f(x, \delta_i) \le \zeta_i, \qquad i = 1, \ldots, N, \\\ & \text{hard (scenario-independent) constraints on } x . \end{aligned}$$

The table shows what $c(x)$ and the constraints are for each program class:

Program	Cost $c(x)$	Scenario constraint	Hard constraint
LP	$c^\top x$	$A(\delta_i)x + b(\delta_i) \le \zeta_i \mathbf{1}$	$Gx + h \le 0$
QP	$c^\top x + \tfrac12 x^\top Q x$	$A(\delta_i)x + b(\delta_i) \le \zeta_i \mathbf{1}$	$Gx + h \le 0$
SDP	$c^\top x + \tfrac12 x^\top Q x$	$F_0(\delta_i) + \sum_j x_j F_j(\delta_i) \preceq \zeta_i I$	$E_0 + \sum_j x_j E_j \preceq 0$

So $f(x,\delta)$ is the largest entry of $A(\delta)x + b(\delta)$ for LP and QP. For SDP it is the largest eigenvalue of $F_0(\delta) + \sum_j x_j F_j(\delta)$.

There is one slack per scenario. It is shared by all $m_s$ rows of $A(\delta_i)$, or by the whole LMI. A "scenario constraint" is therefore the whole block for one $\delta_i$, and $k$ counts scenarios, not rows.

Hard constraints do not depend on $\delta$:

• They are always enforced.
• They are not scenarios.
• They never count towards $k$.
• They play no part in the risk.

The four formulations

Web interface option	$\tau$	$\rho$	Slacks $\zeta_i$
Robust	$0$	–	none ($\zeta_i = 0$)
Regularization	$> 0$	–	none
Relaxation	$0$	$> 0$	yes
Regularization + Relaxation	$> 0$	$> 0$	yes

• Regularization pulls the solution towards a reference point $\bar{x}$ (x_ref) in the $p$-norm.
• Relaxation lets the program violate a scenario constraint by $\zeta_i$, at a cost of $\rho \zeta_i$.
• In the web interface, the selected option decides which of $\tau$ and $\rho$ are used; the fields of the other options are ignored.
• In the Python API, the formulation follows from the arguments tau, rho and include_slack. With the default include_slack=None, slacks are created exactly when rho != 0. The regularization term is added only when tau > 0.
• In the relaxation formulations $\rho$ must be positive. With $\rho = 0$ the slacks cost nothing, so the scenario constraints have no effect. For the interval example, the solver then reports the problem as unbounded.

The theory assumes that the program has a solution and that the solution is unique. The usual ways to secure this are:

• With relaxation, the scenario constraints alone cannot make the program infeasible.
• Regularization ($\tau > 0$) is the usual way to break ties between equally good solutions.

Scen-Opt decides support-list membership by comparing optimal values, not optimizers, so that test still works when the optimizer is not unique.

Below, $x^\ast_N$ is the $x$-part of the optimal solution for the chosen formulation.

3. Appropriateness and risk

A decision $x$ is appropriate for $\delta$ if $f(x,\delta) \le 0$, that is, if $x$ satisfies the original, unrelaxed scenario constraint for $\delta$. With relaxation, $x^\ast_N$ can be inappropriate for some of the scenarios used to compute it: those with $\zeta_i > 0$.

The risk of $x$ is the probability that a new, unseen scenario makes $x$ inappropriate:

$$V(x) = \mathbb{P}\big\lbrace \delta \in \Delta : f(x, \delta) > 0 \big\rbrace .$$

For LP and QP this means at least one row of $A(\delta)x + b(\delta) \le 0$ fails. For SDP it means the matrix is not negative semidefinite.

The quantity of interest is $V(x^\ast_N)$. Three facts shape how it is handled:

• $V(x^\ast_N)$ cannot be computed, because $\mathbb{P}$ is unknown.
• $V(x^\ast_N)$ is itself random, because $x^\ast_N$ depends on the sample.
• Scenario theory bounds $V(x^\ast_N)$ using only $N$, the complexity and the confidence parameter $\beta$. No validation data and no model of $\mathbb{P}$ are needed.

4. Support lists and complexity

A support list is a sublist of $(\delta_1, \ldots, \delta_N)$ with three properties:

1. It contains every scenario that $x^\ast_N$ violates ($f(x^\ast_N, \delta_i) > 0$), together with some of the scenarios at which $x^\ast_N$ is active ($f(x^\ast_N, \delta_i) = 0$).
2. Solving the program with only these scenarios (plus the hard constraints) gives back $x^\ast_N$.
3. It is irreducible: removing any element changes the solution.

More than one support list can exist. The smallest size among them is the complexity $s^\ast_N$. Scen-Opt reports its computed value as $k$:

• In Python, it is the complexity return value of solve_lp, solve_qp and solve_sdp.
• In the web interface, it is the "Complexity (support list size)" row.

In the interval example, the two extreme points fix the interval, so $k = 2$ whatever the value of $N$.

How relaxation affects $k$

A violated scenario belongs to every support list, so $k$ is at least the number of violated scenarios. The penalty $\rho$ therefore trades cost against risk:

• A small $\rho$ makes violations cheap. More scenarios are violated, $k$ grows, and so do the risk bounds.
• A large $\rho$ makes the program behave like the robust one.

The table shows this for the interval example. The settings are: Relaxation formulation, $N = 100$, $\beta = 10^{-6}$, solver CLARABEL. The code is on First Example: Python Calls.

$\rho$	Half-width $c^\top x^\ast$	Violated scenarios	$k$	$[\underline{\varepsilon}, \overline{\varepsilon}]$
0.02	0.2736	50	50	[0.2367, 0.7638]
0.05	0.4205	20	20	[0.0433, 0.4719]
0.1	0.4644	10	10	[0, 0.3432]
0.25	0.4911	4	4	[0, 0.2476]
1.0	0.4983	0	2	[0, 0.2085]

A narrower interval (lower cost) comes with a higher certified risk. The ρ and τ graphs in the web interface plot exactly this trade-off.

Regularization changes which solution is selected, and so it can change $k$ too. In the same example, with $\tau = 0.1$, $\bar{x} = 0$ and $p = 2$ (Regularization + Relaxation), $\rho = 0.05, 0.1, 0.25, 1$ give $k = 22, 11, 5, 2$ instead of $20, 10, 4, 2$.

5. The risk certificate

The confidence parameter $\beta \in (0,1)$ is the probability that the certificate itself is wrong. This probability is taken over the random draw of the $N$ scenarios. Typical values are $10^{-6}$ or smaller; all benchmarks use $\beta = 10^{-6}$. The statements below hold for every distribution $\mathbb{P}$.

5.1 Upper bound (valid in every case)

For $k = 0, \ldots, N-1$, let $t(k)$ be the unique solution in $(0,1)$ of

$$\frac{\beta}{N} \sum_{i=k}^{N-1} \binom{i}{k} t^{\,i-k} - \binom{N}{k} t^{\,N-k} = 0,$$

and set $\varepsilon(k) = 1 - t(k)$, with $\varepsilon(N) = 1$. Then

$$\mathbb{P}^N \big\lbrace V(x^\ast_N) > \varepsilon(s^\ast_N) \big\rbrace \le \beta .$$

This one-sided result needs no assumption beyond those in Section 2 (Garatti & Campi, 2025).

5.2 Two-sided bound (non-degenerate case)

For $k = 0, \ldots, N-1$, consider the polynomial equation

$$\binom{N}{k} t^{\,N-k} - \frac{\beta}{2N} \sum_{i=k}^{N-1} \binom{i}{k} t^{\,i-k} - \frac{\beta}{6N} \sum_{i=N+1}^{4N} \binom{i}{k} t^{\,i-k} = 0 .$$

It has exactly two solutions $\underline{t}(k) \le \overline{t}(k)$ in $[0, \infty)$.

For $k = N$, the equation is instead

$$1 - \frac{\beta}{6N} \sum_{i=N+1}^{4N} \binom{i}{N} t^{\,i-N} = 0 .$$

It has one solution $\overline{t}(N)$, and $\underline{t}(N) = 0$.

The bounds are

$$\underline{\varepsilon}(k) = \max\lbrace 0,\ 1 - \overline{t}(k) \rbrace, \qquad \overline{\varepsilon}(k) = 1 - \underline{t}(k).$$

If the problem is non-degenerate (Section 6), then

$$\mathbb{P}^N \big\lbrace \underline{\varepsilon}(s^\ast_N) \le V(x^\ast_N) \le \overline{\varepsilon}(s^\ast_N) \big\rbrace \ge 1 - \beta$$

(Garatti & Campi, 2022).

5.3 What quantify_risk computes

src.Risk.quantify_risk(k, N, beta) returns (eps_lower, eps_upper), which are $\underline{\varepsilon}(k)$ and $\overline{\varepsilon}(k)$ from Section 5.2. It follows the numerical procedure of Campi & Garatti (2023, Appendix B):

1. It rewrites the polynomial of Section 5.2 in the variable $\varepsilon = 1 - t$, using the regularized incomplete beta function $I_\varepsilon(a,b)$ (jax.scipy.special.betainc(a, b, eps)).
2. The two bounds are then the two zeros in $[0,1]$ of

$$g(\varepsilon) = \frac{\beta}{3} I_\varepsilon(k+1, N-k) + \frac{\beta}{6} I_\varepsilon(k+1, 4N+1-k) - \Big(1 + \frac{\beta}{6N}\Big) N \varepsilon \big[ I_\varepsilon(k, N-k+1) - I_\varepsilon(k+1, N-k) \big].$$

The zeros are found by bisection:

1. Lower bound. Bisect on $[0, k/N]$ until the interval is narrower than $10^{-10}$. Return the left end of the final interval, which is the conservative side. For $k = 0$ the interval is empty and $\underline{\varepsilon} = 0$.
2. Upper bound. If $k = N$, return $\overline{\varepsilon} = 1$. Otherwise, bisect on $[k/N, 1]$ in the same way and return the right end, which is again the conservative side.

On test cases, these values match the roots of the polynomial in Section 5.2 found by direct root-finding, to about $10^{-10}$.

quantify_risk computes in JAX 64-bit floats (it switches jax_enable_x64 on only while it runs). In float32, rounding noise could understate $\overline{\varepsilon}$ by up to about 0.01 at $\beta = 10^{-6}$.

The reported upper bound is also valid without non-degeneracy. Scen-Opt reports $\overline{\varepsilon}(k)$ from Section 5.2, not $\varepsilon(k)$ from Section 5.1. On a grid of cases ($N$ from 5 to 300, $k$ from 0 to $N$, $\beta$ from $10^{-9}$ to $10^{-1}$), $\overline{\varepsilon}(k)$ was never smaller than $\varepsilon(k)$. For example, with $N = 100$, $k = 2$ and $\beta = 10^{-6}$, $\varepsilon(k) = 0.2021$ and $\overline{\varepsilon}(k) = 0.2085$. So $\overline{\varepsilon}$ is a valid, slightly conservative upper bound in the degenerate case too.

$\overline{\varepsilon}$ also increases with $k$. In a degenerate case the computed $k$ can exceed $s^\ast_N$, and the bound then errs on the safe side.

Some example values of quantify_risk:

$N$	$k$	$\beta$	$\underline{\varepsilon}$	$\overline{\varepsilon}$
100	2	$10^{-6}$	0	0.2085
1000	2	$10^{-6}$	0	0.0230
10000	2	$10^{-6}$	0	0.0023
1000	20	$10^{-6}$	0.0039	0.0561
1000	100	$10^{-6}$	0.0544	0.1629
100	2	$10^{-3}$	0	0.1404
100	2	$10^{-9}$	0	0.2680

The bounds bracket $k/N$ and tighten as $N$ grows. $\beta$ has a much weaker effect than $N$ or $k$: making $\beta$ a million times smaller ($10^{-3} \to 10^{-9}$) raises $\overline{\varepsilon}$ from 0.140 to 0.268 when $N = 100$.

5.4 Reading the certificate

"Risk $[0, 0.2085]$ with confidence $1 - 10^{-6}$" means the following. Unless the 100 scenarios were an exceptionally unlucky draw (probability at most $10^{-6}$), at most 20.85% of future $\delta$ fall outside the computed interval. It does not mean the solution fails with probability $\beta$. $\beta$ is the probability that the bound itself is wrong.

6. Degeneracy

The two-sided bound needs the following property.

Non-degeneracy. For every $N$, and with probability one over the scenarios, there is exactly one support list.

Equivalently, the scenarios that are active or violated at the optimum already form a support list: none of them is redundant. Degeneracy means more scenarios are active or violated than are needed to fix the solution. Typical causes:

• Repeated scenarios. Each copy can be dropped on its own, so several support lists exist. This can happen with discrete distributions or rounded data.
• Many constraints through one point. In the Growth Bound benchmark, essentially all of the 3127 scenario constraints are active at the optimum, but only 6 form a support list. In Radiation Therapy, 20 are active and 5 form a support list.

What Scen-Opt does

get_support (in src/Miscellaneous.py, details on Miscellaneous.py) computes $k$ and tests for degeneracy at the same time:

1. Screen. Keep each scenario constraint whose dual value exceeds $10^{-8}$. These are the active or violated ones, called the support scenarios.
2. Check. Test whether these candidates alone, plus the hard constraints, reproduce the optimal value, to a relative tolerance of $10^{-6}$.
3. Prune. If they do, remove candidates greedily while the optimal value stays the same. If the pruning removes a scenario that was genuinely active or violated at the optimum, raise the degeneracy flag.
   • Interior-point solvers can leave a small dual on a constraint that is strictly slack. Removing such a constraint is treated as solver noise and does not raise the flag.
   • When this happens, the console prints a "Discarded ... spuriously screened constraint(s)" message.
4. Fallback. If the screened set does not reproduce the optimum, prune the full scenario list instead and raise the flag as a precaution. The resulting list is valid, but it cannot be certified to be the smallest one.

How degeneracy is reported

• Python API: the last return value of solve_lp, solve_qp and solve_sdp is degeneracy. quantify_risk does not know about it and always returns both bounds. If degeneracy is True, discard eps_lower.
• Web interface: the results table shows Degeneracy Detected? Yes, and the "Risk Bounds" cell shows the lower bound as N/A, e.g. [N/A, 0.1414]. The ρ and τ graphs leave the lower bound of that run out. The downloaded results still contain the numerical lower bound in risk; do not use it when degeneracy is true.
• Benchmarks: run.py prints the bounds as [unreliable, <upper>]. The benchmark table in the README marks the lower bound as N/A (Growth Bound and Radiation Therapy).

In every case the upper bound $\overline{\varepsilon}$ remains valid (Section 5.3).

Limits of the check

The check uses only the scenarios you supplied. It cannot prove that non-degeneracy holds with probability one. If the flag is not raised, it is still the user's responsibility to judge whether the lower bound applies.

Solver accuracy also matters. On a 200-point version of the interval example (the example on LP.py), the first-order solver SCS gives $k = 31$ with the flag raised, while CLARABEL gives $k = 2$, non-degenerate. Use a high-accuracy solver (MOSEK, which is recommended, or CLARABEL) whenever you rely on $k$.

7. Sweeping several ρ or τ values: the union bound

The web interface accepts comma-separated lists for $\rho$ and $\tau$. It solves one program for every pair, $n_\tau \cdot n_\rho$ programs in total, and certifies each one at the entered $\beta$.

• Each certificate fails with probability at most $\beta$, so the probability that any of them fails is at most $n_\tau n_\rho \beta$ (union bound).
• All certificates therefore hold simultaneously with confidence at least

$$1 - \beta \, n_\tau \, n_\rho .$$

This is what allows you to look at all the results and then pick a $(\rho, \tau)$ pair. The certificate of the chosen pair holds at the overall level, even though the choice was made after seeing the data.

The interface reports this overall value:

• The confidence label changes to $1 - \beta \cdot n_\rho$, $1 - \beta \cdot n_\tau$ or $1 - \beta \cdot n_\tau \cdot n_\rho$, depending on the formulation.
• The "Confidence" row of the results table shows the same value with all its digits (e.g. 0.999992), clamped at 0.

For example, 4 values of $\rho$ and 2 of $\tau$ with $\beta = 10^{-6}$ give $1 - 8 \times 10^{-6} = 0.999992$.

For an overall confidence of $1 - \beta_{\text{tot}}$, enter $\beta = \beta_{\text{tot}} / (n_\tau n_\rho)$. Because $\beta$ has only a weak effect on the bounds (Section 5.3), the bounds widen only slightly.

The Python API does no such bookkeeping. If you loop over several values yourself and compare the results, apply the same union bound.

References

The theory on this page is summarized from the works below.

• M. C. Campi and S. Garatti. Introduction to the Scenario Approach. SIAM, 2018.
• G. C. Calafiore and M. C. Campi. The scenario approach to robust control design. IEEE Transactions on Automatic Control, 51(5):742–753, 2006.
• M. C. Campi and S. Garatti. The exact feasibility of randomized solutions of uncertain convex programs. SIAM Journal on Optimization, 19(3):1211–1230, 2008.
• M. C. Campi and S. Garatti. Wait-and-judge scenario optimization. Mathematical Programming, 167(1):155–189, 2018.
• M. C. Campi and S. Garatti. Scenario optimization with relaxation: a new tool for design and application to machine learning problems. 2020 59th IEEE Conference on Decision and Control (CDC), pp. 2463–2468, 2020. doi:10.1109/CDC42340.2020.9303914
• M. C. Campi, A. Carè and S. Garatti. The scenario approach: A tool at the service of data-driven decision making. Annual Reviews in Control, 52:1–17, 2021.
• S. Garatti and M. C. Campi. Risk and complexity in scenario optimization. Mathematical Programming, 191(1):243–279, 2022.
• M. C. Campi and S. Garatti. Compression, generalization and learning. Journal of Machine Learning Research, 24(339):1–74, 2023.
• S. Garatti and M. C. Campi. Non-convex scenario optimization. Mathematical Programming, 209(1):557–608, 2025.