Nonlinear Reachability

Chapter 1: Interval Arithmetic

Before we can propagate sets through nonlinear functions, we need a calculus for intervals. Interval arithmetic defines how to add, subtract, multiply, and divide intervals so the result is guaranteed to contain every possible outcome.

Start with the simplest case. Suppose x ∈ [1, 3] and y ∈ [2, 5]. What values can x + y take? The smallest sum is 1 + 2 = 3. The largest is 3 + 5 = 8. So x + y ∈ [3, 8]. That is interval addition:

Each operation follows the same principle: consider all possible combinations and take the tightest enclosing interval.

The Four Basic Operations

Subtraction is where it gets interesting. If x ∈ [1, 3] and y ∈ [2, 5], then x − y ranges from 1 − 5 = −4 to 3 − 2 = 1. Notice the endpoints swap: we subtract the upper bound of y from the lower bound of x to get the new lower bound.

Multiplication is messier because signs matter. If one interval is entirely positive and the other entirely negative, the extremes swap. The safest rule: compute all four products of the endpoints and take the min and max.

Elementary Functions

For monotone functions, interval evaluation is straightforward. If f is monotonically increasing on [a, b], then f([a, b]) = [f(a), f(b)]. If decreasing, f([a, b]) = [f(b), f(a)]. For non-monotone functions like sin or x², we must find the extrema within the interval.

The squaring rule for intervals containing zero is important: the minimum of x² is 0 (achieved at x = 0), not min(a², b²).

python
class Interval:
    def __init__(self, lo, hi):
        self.lo, self.hi = lo, hi

    def __add__(self, other):
        return Interval(self.lo + other.lo, self.hi + other.hi)

    def __sub__(self, other):
        return Interval(self.lo - other.hi, self.hi - other.lo)

    def __mul__(self, other):
        prods = [self.lo*other.lo, self.lo*other.hi,
                 self.hi*other.lo, self.hi*other.hi]
        return Interval(min(prods), max(prods))

# Example: [1,3] + [2,5] = [3,8]
x = Interval(1, 3)
y = Interval(2, 5)
z = x + y  # Interval(3, 8)

Chapter 2: Inclusion Functions

Interval arithmetic gives us operations on intervals. But what we really need is to evaluate an entire function on an interval input. An inclusion function is a function that maps intervals to intervals and is guaranteed to contain the true image.

Here [x] denotes an interval, and [f̅]([x]) is the inclusion function's output. The guarantee is: no matter which x you pick from [x], f(x) will land inside [f̅]([x]). This is the fundamental enclosure property.

The simplest way to build an inclusion function is natural interval evaluation: take the expression for f, replace every real operation with its interval counterpart, and evaluate. If f(x) = x² + 2x + 1, then [f̅]([x]) = [x]² + 2[x] + [1, 1].

Properties of Good Inclusion Functions

A useful inclusion function should satisfy three properties:

Thin convergence says that when the input interval shrinks to a point, the output interval also shrinks to a point — the true function value. This is critical: it means that by subdividing the input into small enough pieces, we can make the overapproximation as tight as we want.

Inclusion isotonicity says that a wider input produces a wider output. This is a monotonicity property of the inclusion function itself (not of f). Natural interval evaluation always satisfies this because each interval operation is isotone.

Convergence Order

How fast does the overapproximation shrink as the input interval narrows? The excess width is w([f̅]([x])) − w(f([x])), where w denotes interval width. If this excess is O(w([x])^k), we say the inclusion function has convergence order k.

Natural interval evaluation has order 1 in general: halving the input width roughly halves the excess. The mean value form (Chapter 4) achieves order 2: halving the input width reduces the excess by a factor of 4. Higher-order Taylor models can achieve order k for any desired k.

Chapter 3: The Dependency Effect

Here is the central weakness of natural interval evaluation: it treats every occurrence of a variable as independent. When the same variable appears more than once in an expression, the method does not track the correlation between the occurrences. This leads to overapproximation called the dependency problem (or wrapping effect for set propagation).

A Concrete Example: x² − x on [−1, 2]

Let f(x) = x² − x. What is the true range of f on [−1, 2]?

Taking the derivative: f'(x) = 2x − 1 = 0 gives x = 0.5. We check the critical point and endpoints:

So the true image is [−0.25, 2].

Now let's try natural interval evaluation. We compute [−1, 2]² − [−1, 2]:

Natural interval evaluation gives [−4, 5]. The true answer was [−0.25, 2]. That is a width of 9 instead of 2.25 — a 4× overestimate.

When is the Dependency Problem Worst?

The dependency problem is worst when:

1. The same variable appears many times. Each additional occurrence introduces a new "phantom" degree of freedom that inflates the result.

2. The interval is wide. On a tiny interval, all occurrences of x are near the same value, so treating them as independent barely matters. On a wide interval, the discrepancy is large.

3. The function is far from monotone. Monotone functions have no dependency problem at all (the range is just [f(a), f(b)] or [f(b), f(a)]). It is the turning points inside the interval that create the trouble.

Chapter 4: The Mean Value Theorem Trick

Natural interval evaluation is crude: it replaces every real operation with an interval operation and hopes for the best. The mean value form is smarter. It exploits calculus — specifically, the mean value theorem — to build a tighter inclusion function.

The idea: linearize f around the center of the interval, then bound the linearization error using interval evaluation of the derivative. For a function f: R → R and interval [x] with center c = (a + b) / 2:

We do not know ξ, but we know it lies in [x]. So f'(ξ) ∈ f'([x]) (evaluated using interval arithmetic on the derivative). The mean value inclusion function is:

Here f(c) is a real number (the function evaluated at the center), f'([x]) is the interval-evaluated derivative, and ([x] − c) is the interval [a − c, b − c] = [−w/2, w/2] where w is the interval width.

Worked Example

Let f(x) = x² − x on [−1, 2]. The center is c = 0.5. f(c) = 0.25 − 0.5 = −0.25. The derivative is f'(x) = 2x − 1, so f'([−1, 2]) = 2[−1, 2] − 1 = [−2, 4] − 1 = [−3, 3]. And [x] − c = [−1.5, 1.5].

Width 9.0 — barely better than the natural evaluation's width 9.0 for this example. Why? Because f'([x]) = [−3, 3] is a wide interval centered at 0, and multiplying by the half-width gives a huge spread.

Convergence Order

The mean value form has quadratic convergence: if the interval width is w, the excess width is O(w²). This means halving the input interval reduces the excess by a factor of 4, compared to only 2 for natural evaluation. On small intervals, this makes a huge difference.

For multivariable functions f: Rⁿ → R^m, the mean value form uses the Jacobian matrix:

where J_f([x]) is the interval-evaluated Jacobian. Each entry of the Jacobian is an interval, and the matrix-vector product uses interval arithmetic.

Chapter 5: Taylor Expansions

The mean value form uses a first-order Taylor expansion and bounds the remainder with interval arithmetic on the derivative. The natural next step: use higher-order Taylor expansions for tighter bounds.

Second-Order Form

Expand f around the center c of interval [x]:

The first two terms are exact (evaluated at the point c). Only the last term needs bounding. Since ξ ∈ [x], we bound f''(ξ) ∈ f''([x]) using interval arithmetic:

Notice that f'(c) is a real number (no interval). This means the linear term contributes no overestimation — it shifts the interval exactly. All overestimation is confined to the quadratic remainder term.

Worked Example

f(x) = x² − x on [−1, 2], c = 0.5.

f(c) = −0.25, f'(c) = 0 (the derivative at the center is 2(0.5) − 1 = 0), f''(x) = 2 (constant).

This is the exact answer! Because f is a quadratic polynomial, the second-order Taylor expansion captures it perfectly — there is no remainder.

The General Pattern

The first k−1 terms are exact point evaluations. Only the k-th term is bounded with interval arithmetic. This structure is key: the more terms you compute exactly, the less overestimation you incur.

The Hessian in Multiple Dimensions

For f: Rⁿ → R, the second-order form involves the Hessian matrix H_f:

The Hessian H_f([x]) is an n × n matrix of intervals. Computing it requires n(n+1)/2 second partial derivatives, each evaluated on the input interval. For n = 100 (a small neural network layer), that is 5050 interval Hessian entries — expensive, but still cheaper than naive Monte Carlo sampling for formal guarantees.

Chapter 6: Taylor Models

Taylor expansions bound the function output as an interval. But an interval is a box — it cannot capture correlations between output dimensions. Taylor models go further: the output is a polynomial (which can represent correlations) plus a small interval remainder.

Here p(x) is a polynomial of degree k (the Taylor polynomial), and [R] is an interval bounding the approximation error. The polynomial captures the smooth part of f, and the remainder captures everything else.

Why Polynomials Beat Boxes

Consider a 2D function f: R² → R² where f₁ and f₂ are correlated (e.g., f₁(x) = x, f₂(x) = x²). The true image is a parabolic curve. An interval box bounds this curve with a rectangle, wasting huge area in the corners. A Taylor model represents the output as (x, x² + [R]), which traces the parabola with only a thin error band — much tighter.

Arithmetic on Taylor Models

Just as interval arithmetic defines +, −, · for intervals, Taylor model arithmetic defines operations on polynomial-plus-remainder pairs:

Addition is straightforward: add polynomials, add remainders. Multiplication is trickier: the product of two degree-k polynomials is degree 2k. We truncate to degree k, sweeping the higher-order terms into the remainder interval. This is where overapproximation enters.

Composing Taylor Models

A major advantage of Taylor models: they compose naturally. If you have a Taylor model for g(x) and want to compute f(g(x)), you substitute the Taylor model for g into the Taylor expansion of f. The polynomial part remains polynomial (after truncation), and the remainder absorbs all the error terms.

This is exactly what happens when propagating sets through a deep network layer by layer. Each layer's Taylor model feeds into the next. The remainder grows with each composition — but much slower than with plain interval arithmetic, because the polynomial captures the dominant behavior.

Chapter 7: Symbolic vs. Concrete

We have seen two strategies: concrete methods (interval arithmetic, mean value form) that evaluate the function on intervals and produce interval outputs, and symbolic methods (Taylor models) that maintain a polynomial expression in the input variables. Each has distinct tradeoffs.

Concrete (Interval) Methods

Concrete methods evaluate greedily: at each operation, they produce an interval and forget the structure that produced it. This is fast and simple, but information is lost at every step — especially the dependency between variables.

Symbolic (Taylor Model) Methods

Symbolic methods carry forward the algebraic relationship between inputs and outputs. A Taylor model knows that "output dimension 1 is approximately 2x + 3y" and "output dimension 2 is approximately x − y". This structural knowledge means that when you later compute dimension 1 + dimension 2, the x and y terms combine correctly, respecting dependencies.

The Composition Tradeoff

Symbolic methods shine when composing many functions (deep networks, long time horizons). But the polynomial complexity grows with each composition. After L layers of a neural network with n neurons each, a degree-k Taylor model has O((kL + n choose n)) terms — exponential in L for fixed k. In practice, you must periodically concretize: evaluate the polynomial on its input intervals, collapse it to a new interval, and start a fresh Taylor model. This resets the polynomial complexity at the cost of introducing one round of wrapping.

Chapter 8: Optimization & MILP

So far we have propagated sets forward using arithmetic. There is an entirely different approach: formulate the reachability question as an optimization problem. Instead of asking "what interval does f([x]) land in?", ask "what is the maximum (or minimum) of f over [x]?"

If we can solve these two optimization problems, the range of f on [x] is exactly [f_min, f_max] — no overapproximation at all. The catch: for general nonlinear f, global optimization is NP-hard.

Support Functions

A support function generalizes the idea to any direction d. Instead of finding the range along coordinate axes, find the extreme of the dot product d^Tf(x):

By evaluating ρ in many directions, you can reconstruct a tight polytope bounding the reachable set. This is the support function approach to reachability, and it reduces set propagation to a sequence of optimization problems.

ReLU Networks and MILP

For neural networks with ReLU activations, there is a beautiful trick. The ReLU function max(0, x) is piecewise linear. Its behavior at each neuron can be encoded as a binary variable: z = 1 if the neuron is active (x ≥ 0), z = 0 if inactive (x < 0).

Here L and U are lower and upper bounds on x (computed by interval arithmetic), and z ∈ {0, 1}. These four constraints, together with the binary variable z, exactly encode the ReLU. The entire network becomes a mixed-integer linear program (MILP).

The MILP approach gives exact bounds on the network output — no overapproximation. But it is exponential in the number of unstable neurons. For a network with 1000 unstable neurons, 2¹⁰⁰⁰ branches is intractable. In practice, modern MILP solvers use branch-and-bound with LP relaxations to prune the tree, making problems with hundreds of unstable neurons tractable.

Chapter 9: Partitioning & Neural Nets

All the methods so far propagate a single set through the nonlinear function. There is a universal trick to tighten any overapproximation: subdivide the input.

Split the input interval [x] into smaller pieces [x]₁, [x]₂, ..., [x]_k. Compute the inclusion function on each piece. The union of the results still contains the true image, but each piece's overapproximation is tighter (because the dependency problem is less severe on narrow intervals).

Neural Network Verification

A neural network with L layers and n_l neurons per layer is a composition f = f_L ∘ f_L−1 ∘ … ∘ f₁, where each f_l is an affine map followed by a nonlinear activation. The verification question: given input set X, does the output ever enter an unsafe region?

The Combinatorial Explosion

For ReLU networks, the number of possible activation patterns is 2^N where N is the total number of neurons. Each pattern corresponds to a linear region of the network. Verification requires reasoning about which patterns are reachable from the input set.

Partitioning the input helps by reducing the number of unstable neurons in each partition. If a partition is small enough that every neuron's sign is determined, the network is linear on that partition and verification is trivial (LP, not MILP). The challenge: finding the right partition granularity. Too coarse = too many unstable neurons. Too fine = too many partitions.

What comes next: This chapter handled nonlinear continuous systems. Chapter 10 tackles discrete systems, where the state space is finite and reachability becomes a graph search problem. The tools change from calculus to combinatorics, but the core question remains: where can the system go?