Skip to content

Risk_weeds

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

Risk.py: quantify_risk

Source: src/Risk.py · Theory: Scenario Theory · Where $k$ comes from: get_support

quantify_risk(k, N, beta) turns the complexity $k$ of a solved scenario program, the number of scenarios $N$ and a confidence parameter $\beta$ into a lower bound $\underline{\varepsilon}$ and an upper bound $\overline{\varepsilon}$ on the risk $V(x^\ast_N)$, the probability that a new scenario is violated. If the problem is non-degenerate, then with probability at least $1-\beta$:

$$\underline{\varepsilon}(k) \;\le\; V(x^\ast_N) \;\le\; \overline{\varepsilon}(k).$$

These are the two-sided bounds of Garatti and Campi, Risk and complexity in scenario optimization, Math. Programming 191 (2022). They are computed by bisection on a reformulation that uses the regularized incomplete beta function, following Campi and Garatti, Compression, generalization and learning, JMLR 24 (2023), Appendix B. The function has no notion of degeneracy; see Degeneracy.

Signature

from src.Risk import quantify_risk

epsL, epsU = quantify_risk(k, N, beta)

Parameters

Name Type Meaning
k int, $0 \le k \le N$ Complexity, the cardinality of the support list. This is the complexity returned by solve_lp, solve_qp or solve_sdp. Integral floats such as 2.0 give the same result.
N int, $N \ge 1$ Number of scenarios.
beta float, $0 < \beta < 1$ Confidence parameter. The bounds hold with probability at least $1 - \beta$. Typical values are 1e-6 to 1e-3.

The inputs are validated before any computation: ValueError is raised unless N is an integer $\ge 1$, k is an integer between 0 and N, and $0 < \beta < 1$ (see Edge cases).

Returns

Name Type Meaning
epsL float Lower bound $\underline{\varepsilon}(k)$. It is 0.0 whenever the theoretical lower bound is $\le 0$, and always when k == 0.
epsU float Upper bound $\overline{\varepsilon}(k)$. It is 1.0 when k == N, and also when the upper root lies within $10^{-10}$ of 1 (for example k=99, N=100, beta=1e-6).

Both values are plain Python floats.

What is computed

For $k < N$, the theorem uses the polynomial equation in $t$:

$$\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 roots $\underline{t}(k) \le \overline{t}(k)$ in $[0, \infty)$. For $k = N$ the equation is $1 - \frac{\beta}{6N}\sum_{i=N+1}^{4N}\binom{i}{N} t^{i-N} = 0$, which has one root $\overline{t}(N)$, and one sets $\underline{t}(N) = 0$. Then

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

The equation the code actually solves

Instead of summing the polynomial's binomial terms (up to $4N$ of them), the code works in $\varepsilon = 1 - t$ and solves $L(\varepsilon) = R(\varepsilon)$ with

$$\begin{aligned} L(\varepsilon) &= \frac{\beta}{3}\, I_\varepsilon(k+1,\ N-k) \;+\; \frac{\beta}{6}\, I_\varepsilon(k+1,\ 4N+1-k), \\\ R(\varepsilon) &= \Big(1 + \frac{\beta}{6N}\Big)\, N\varepsilon\, \Big( I_\varepsilon(k,\ N-k+1) - I_\varepsilon(k+1,\ N-k) \Big), \end{aligned}$$

where $I_x(a,b)$ is the regularized incomplete beta function (jax.scipy.special.betainc(a, b, x)).

Why the two equations agree. Write $p = \varepsilon$ and $\mathrm{Bin}(M, p)$ for a binomial count. Three standard identities are needed:

$$\begin{aligned} p^{k+1} \sum_{i=k}^{M-1} \binom{i}{k} (1-p)^{i-k} &= I_p(k+1,\ M-k) = \mathbb{P}\big[\mathrm{Bin}(M,p) \ge k+1\big], \\\ I_p(k,\ N-k+1) - I_p(k+1,\ N-k) &= \binom{N}{k} p^k (1-p)^{N-k} = \mathbb{P}\big[\mathrm{Bin}(N,p) = k\big], \\\ I_p(k+1,\ N+1-k) &= I_p(k+1,\ N-k) + p\,\mathbb{P}\big[\mathrm{Bin}(N,p) = k\big]. \end{aligned}$$

The first is the negative-binomial identity. Multiply the polynomial by $N\varepsilon^{k+1}$, then rewrite $\sum_{i=N+1}^{4N}$ as $\sum_{i=k}^{4N} - \sum_{i=k}^{N}$. The result is exactly $L = R$. In words:

$$\Big(1+\tfrac{\beta}{6N}\Big) N\varepsilon\, \mathbb{P}[\mathrm{Bin}(N,\varepsilon)=k] \;=\; \tfrac{\beta}{3}\,\mathbb{P}[\mathrm{Bin}(N,\varepsilon)\ge k+1] + \tfrac{\beta}{6}\,\mathbb{P}[\mathrm{Bin}(4N+1,\varepsilon)\ge k+1].$$

So $\underline{\varepsilon}$ and $\overline{\varepsilon}$ are the two roots of one equation in $\varepsilon$. In the upper-bound loop the code writes the first coefficient as (beta/2 - beta/6). This equals $\beta/3$, so both loops solve the same equation and differ only in the bracket they search.

The same equation also covers the edge cases. When $k = N$, the term $I_\varepsilon(N+1, 0)$ must be 0. When $k = 0$, the term $I_\varepsilon(0, N+1)$ must be 1. JAX's betainc returns exactly these limits (betainc(11, 0, 0.3) = 0.0 and betainc(0, 11, 0.3) = 1.0), whereas scipy.special.betainc returns nan for both. Replacing the JAX function with SciPy's would break the k = 0 and k = N cases.

Bisection

threshold = 1e-10

# lower bound: search [0, k/N]
t1, t2 = 0.0, k/N
while t2 - t1 > threshold:
    t = (t1 + t2) / 2
    if L(t) > R(t): t1 = t        # still left of the lower root
    else:           t2 = t
epsL = t1                         # left end, so rounded down

# upper bound: search [k/N, 1], skipped when k == N (then epsU = 1)
t1, t2 = k/N, 1
while t2 - t1 > threshold:
    t = (t1 + t2) / 2
    if L(t) > R(t): t2 = t        # right of the upper root
    else:           t1 = t
epsU = t2                         # right end, so rounded up
  • Brackets. The search assumes that the empirical violation frequency $k/N$ lies between the two roots, that is $R(k/N) \ge L(k/N)$. This held in every case checked below.
  • Rounding. Both results are rounded outwards, by at most $10^{-10}$: epsL is the left end of its final interval and epsU the right end.
  • No lower root. If $R \ge L$ on all of $(0, k/N]$, which happens when the theoretical lower bound is $\le 0$, the lower loop shrinks to 0 and returns epsL = 0.0. This reproduces the $\max\lbrace 0, \cdot\rbrace$.
  • Iteration count. The loops run $\lceil \log_2((k/N)/10^{-10}) \rceil$ and $\lceil \log_2((1-k/N)/10^{-10}) \rceil$ times, so at most 34 each. For example, $k=2$, $N=200$ takes 27 + 34 iterations. Each iteration makes four betainc calls.

Floating-point precision

quantify_risk runs the bisection inside jax.experimental.enable_x64(), so JAX computes in float64. Without it, JAX computes in float32. According to the source comment, $L$ and $R$ differ by only about $10^{-8}$ near the upper root, which float32 cannot resolve. The upper bound could then come out up to about 0.01 too low at $\beta = 10^{-6}$, whereas float64 matches a SciPy reference to about $10^{-10}$. The setting is local to the call: importing src.Risk does not change JAX's global configuration, and other JAX code in the same process keeps the float32 default.

Verification against the polynomial

The table compares quantify_risk with roots of the polynomial computed independently. Those roots were found by bisection in 80-digit decimal arithmetic with exact integer binomials. The last column is the one-sided bound $\varepsilon(k) = 1 - t(k)$, where $t(k)$ solves $\frac{\beta}{N}\sum_{i=k}^{N-1}\binom{i}{k}t^{i-k} = \binom{N}{k}t^{N-k}$. This bound needs no non-degeneracy assumption, and quantify_risk does not compute it.

$k$ $N$ $\beta$ epsL exact $\underline{\varepsilon}$ epsU exact $\overline{\varepsilon}$ one-sided $\varepsilon(k)$
0 100 1e-6 0 0 0.1586324059 0.1586324059 0.15244298
5 100 1e-6 0 0 0.2653136026 0.2653136027 0.25876431
20 100 1e-6 0.0432916579 0.0432916580 0.4718639313 0.4718639313 0.46530794
60 100 1e-6 0.3194619689 0.3194619689 0.8398057520 0.8398057520 0.83581827
99 100 1e-6 0.7925930397 0.7925930398 1 0.9999999999 ≈ 1
100 100 1e-6 0.8186299158 0.8186299159 1 1 1
0 1000 1e-6 0 0 0.0172040952 0.0172040952 0.01648037
10 1000 1e-6 0 0 0.0393671338 0.0393671338 0.03843961
50 1000 1e-3 0.0258867719 0.0258867720 0.0851435092 0.0851435092 0.08354206
200 1000 1e-6 0.1333601984 0.1333601985 0.2792148898 0.2792148898 0.27742067

In every row the code agrees with the exact roots to within $10^{-10}$.

Degeneracy

quantify_risk does not take the degeneracy flag and always returns both bounds. The two-sided certificate needs non-degeneracy, meaning the support list is unique. So when solve_lp, solve_qp or solve_sdp returns degeneracy=True:

  • Do not use epsL. It is not certified.
  • epsU can still be used as an upper bound. In every case in the table above, the two-sided epsU was larger, and therefore more conservative, than the one-sided bound $\varepsilon(k)$ that holds without non-degeneracy.
  • k itself may be larger than the true complexity. Pruning finds an irreducible support list, but not necessarily a smallest one. A larger $k$ only makes epsU more conservative.

The web interface's results table has a row "Degeneracy Detected?" (Yes/No). When it says Yes, the "Risk Bounds" row shows the lower bound as N/A (e.g. [N/A, 0.1414]) and the ρ/τ graphs leave it out. The downloaded risk field still holds both numbers.

Edge cases

Input Behaviour
k == 0 The lower interval is [0, 0], so the loop is skipped and epsL = 0.0. The upper search uses $I_\varepsilon(0, N+1) = 1$.
k == N epsU = 1.0 with no search. The lower search uses $I_\varepsilon(N+1, 0) = 0$.
Upper root within $10^{-10}$ of 1 epsU stays at its starting value, 1.0.
Integral floats, e.g. k=2.0 Accepted and converted to int.
k < 0, k > N or non-integral k ValueError, e.g. k must be an integer between 0 and N = 3 (got 5).
N < 1 or non-integral N ValueError, e.g. N must be a positive integer (got 0).
beta outside $(0, 1)$ ValueError, e.g. beta must lie strictly between 0 and 1 (got 0).

Timing

Step Measured cost
import jax together with src.Risk about 0.3 s
First call in a process about 200 ms (JAX warm-up)
Later calls about 20–35 ms each

The run time does not depend on $N$: quantify_risk(5000, 1_000_000, 1e-6) took 32 ms. This is because betainc is a special-function evaluation rather than a sum over $4N$ terms, and the iteration count depends only on $k/N$. The JAX calls run eagerly and are not JIT-compiled. The web interface reports this time separately from the solve time, as risk_time.

How the web interface uses beta

For each $(\tau, \rho)$ pair in a sweep, the web interface's solve code (solve_form) calls quantify_risk(complexity, N, beta) with the $\beta$ the user entered. It reports an overall confidence of $1 - q\beta$, where $q$ is the number of $(\tau,\rho)$ pairs. This is a union bound over the $q$ attempts. The MOSEK licence-upload path runs the same code, so it reports the same values. A missing $\beta$, or one outside $(0, 1)$, is rejected before any solve (in the browser, and by the server with HTTP 400).

Example

Run from the repository root:

import time
from src.Risk import quantify_risk

for k, N, beta in [(2, 200, 1e-6), (0, 200, 1e-6), (10, 200, 1e-6), (200, 200, 1e-6), (25, 1000, 1e-3)]:
    t0 = time.perf_counter()
    epsL, epsU = quantify_risk(k, N, beta)
    dt = time.perf_counter() - t0
    print(f"k={k:>3} N={N:>4} beta={beta:.0e}:  epsL={epsL:.6f}  epsU={epsU!r:<22} ({dt*1e3:.0f} ms)")

Output:

k=  2 N= 200 beta=1e-06:  epsL=0.000000  epsU=0.11005086670164022    (223 ms)
k=  0 N= 200 beta=1e-06:  epsL=0.000000  epsU=0.0829441822716035     (21 ms)
k= 10 N= 200 beta=1e-06:  epsL=0.000000  epsU=0.18516078560496682    (33 ms)
k=200 N= 200 beta=1e-06:  epsL=0.901218  epsU=1.0                    (18 ms)
k= 25 N=1000 beta=1e-03:  epsL=0.009714  epsU=0.05192692162963794    (33 ms)

What the output shows:

  • The first call includes the JAX warm-up, which is why it takes longer.
  • For k == N, epsU is 1.0, while epsL is still informative.
  • The first and third rows match the robust and relaxed runs of the LP example, which have $k = 2$ and $k = 10$ with $N = 200$.

See also

Using The Tool

Home

Starting Up

Convex Optimization

Graphical Results

Getting into the Weeds

Scenario Theory

Implementation Details

Clone this wiki locally