Approximate Value Functions

Chapter 1: Parametric Representations

A parametric approximation U_θ(s) represents the value function as a fixed functional form controlled by a parameter vector θ. Instead of storing one number per state, we store θ and compute U_θ(s) on demand for any query state.

The most important categorization is local vs global:

In both cases, the fitting procedure finds θ by minimizing a loss between U_θ(s_i) and the backed-up target values y_i. The loss is typically mean squared error:

Fitting the Approximator: Regression on Bellman Targets

In approximate value iteration, we compute Bellman targets y_i at each sample state and then fit the approximator to them. The full loop is:

1. Initialize: U_θ(s) = 0 for all s (or some prior)
2. Sample states: Draw m states s₁, ..., s_m from S (uniform, grid, or policy rollouts)
3. Backup: y_i = max_a[R(s_i,a) + γ ∑_s' T(s'|s_i,a) U_θ(s')] for each i
4. Fit: Find θ that minimizes ∑_i(U_θ(s_i) − y_i)²
5. Repeat steps 3–4 until convergence or budget exhausted

The backup in step 3 requires evaluating U_θ(s') at the next state s', which uses the current approximation. This creates the same "moving target" instability as neural network regression (Chapter 7). With linear and k-NN approximators, the instability is less severe but still exists — there is no general proof of convergence for approximate value iteration.

One exception: If the approximation family is rich enough to contain the true value function, and the sample distribution covers S well, approximate value iteration converges to a small error ball. The formal guarantee: ||U_θ* − U*||_∞ ≤ ε/(1−γ) where ε is the approximation error of the best-fit U* in the function class. This bound degrades with discount γ → 1 and improves as ε → 0.

Chapter 2: Nearest Neighbor

Nearest neighbor (k-NN) is the simplest local approximator. To compute U_θ(s) at a query state s:

1. Find the k closest sample states to s (by Euclidean distance or any metric)
2. Return the average of their stored values

where N_k(s) is the set of k-nearest sample indices and θ_i = y_i is the backed-up value at sample s_i. With k=1, the approximation is piecewise constant — every query returns the value of its single nearest neighbor. As k increases, the estimate becomes smoother but less responsive to local variation.

A critical design choice: distance metric. Raw Euclidean distance treats all dimensions equally, which fails when state variables have very different scales (e.g., angle in radians vs. velocity in m/s). Always normalize state variables to comparable ranges, or use the Mahalanobis distance which accounts for scale and correlation.

python
import numpy as np

class KNNValueFunction:
    """k-Nearest Neighbor value function approximation."""
    def __init__(self, k=3):
        self.k = k
        self.states = []   # list of sample states (arrays)
        self.values = []   # list of target values (floats)

    def update(self, states, values):
        """Replace stored samples with new (state, value) pairs."""
        self.states = [np.array(s) for s in states]
        self.values = list(values)

    def __call__(self, s):
        """Query value at state s."""
        s = np.array(s)
        dists = [np.linalg.norm(s - si) for si in self.states]
        idx = np.argsort(dists)[:self.k]
        return np.mean([self.values[i] for i in idx])

# Usage in approximate value iteration
approx_U = KNNValueFunction(k=3)
sample_states = [np.array([s]) for s in np.linspace(0, 10, 20)]
for iteration in range(50):
    targets = [bellman_backup(s, approx_U) for s in sample_states]
    approx_U.update(sample_states, targets)

Efficient k-NN with k-d Trees

Naively finding the k nearest neighbors requires scanning all m samples: O(m) per query. For large m, this is too slow. The k-d tree (k-dimensional tree) is a binary tree that partitions the state space recursively, enabling nearest-neighbor queries in O(log m) average time.

In Python: scipy.spatial.cKDTree builds a k-d tree in milliseconds for m=10,000 points, and queries are 100× faster than brute force. For the approximate value iteration setting where we make thousands of queries per iteration, this speedup is critical.

python
from scipy.spatial import cKDTree
import numpy as np

class FastKNNValueFunction:
    def __init__(self, k=3):
        self.k = k; self.tree = None; self.values = None

    def update(self, states, values):
        self.tree = cKDTree(np.array(states))  # O(m log m) build
        self.values = np.array(values)

    def __call__(self, s):
        dists, idx = self.tree.query(np.array(s), k=self.k)  # O(log m)
        return np.mean(self.values[idx])

Chapter 3: Kernel Smoothing

Kernel smoothing is a weighted average over all (or nearby) sample states, where the weight is a decreasing function of distance. Common kernels include the Gaussian and the uniform (box) kernel:

Here K is the kernel function and h is the bandwidth — the single most important hyperparameter. Small h: tight kernel, high local sensitivity, noisy estimate. Large h: wide kernel, smooth estimate, may miss local structure.

K-NN is actually a special case of kernel smoothing with the "top-k" kernel: w=1/k for the k nearest points and 0 for the rest. The Nadaraya-Watson estimator (kernel smoothing) can be derived from a nonparametric regression model where U(s) is locally constant in a neighborhood of size h. Maximizing the weighted log-likelihood under a Gaussian noise model gives exactly the kernel smoothing formula above. This probabilistic interpretation means the bandwidth h has a natural interpretation: it is the standard deviation of the assumed "spread" in the value function, not just an arbitrary smoothing parameter.

For value function approximation specifically, a key practical issue is extrapolation. Both k-NN and kernel smoothing extrapolate poorly: they return the nearest sample's value (k-NN) or a weighted average of nearby samples (kernel) for query states outside the convex hull of the training data. This means that states that are never visited during training have poorly estimated values. For problems with bounded state spaces, this is manageable; for problems where the optimal policy may visit novel states (exploration), this creates a chicken-and-egg problem: you need good values to explore, but you need exploration to get good values. This is why on-policy sampling (Chapter 9 of the textbook) helps: it focuses value estimation where the current policy actually goes, ensuring the approximation is good where it matters.

Kernel smoothing generalizes naturally to multivariate states using a product kernel: w(s, s_i) = ∏_j K((s_j − s_ij)/h_j) where h_j is the bandwidth for dimension j. Using separate bandwidths per dimension is important when state variables have different scales or different degrees of variation. For example, in a pendulum MDP, the angle variable may need a narrow bandwidth (value changes quickly near the unstable equilibrium) while the angular velocity may need a wider bandwidth.

Choosing the Bandwidth: Cross-Validation and the Bias-Variance Tradeoff

The optimal bandwidth minimizes the mean squared error between the approximation and the true value function. This has two competing components:

The bias-variance tradeoff is the fundamental problem of statistical estimation. For the value function approximation setting, we typically use leave-one-out cross-validation to select h: fit the approximation on all but one sample state and measure the prediction error at the held-out point. Average over all held-out points and choose the h that minimizes this cross-validation error.

An alternative: Silverman's rule of thumb sets h proportional to the sample standard deviation and inversely to m^1/(n+4) where n is the state dimension. This rule is optimal for Gaussian-distributed data and works well in practice even for non-Gaussian cases. For a 1D problem with standard deviation σ and m samples: h = 1.06 σ m^−1/5.

python
import numpy as np

def silverman_bandwidth(states):
    """Silverman's rule of thumb for kernel bandwidth (1D)."""
    m = len(states)
    sigma = np.std(states)
    return 1.06 * sigma * m ** (-0.2)

def gaussian_kernel_smooth(query_states, sample_states, values, h):
    """Vectorized Gaussian kernel smoothing."""
    diffs = query_states[:, None] - sample_states[None, :]  # (nq, m)
    weights = np.exp(-0.5 * (diffs / h) ** 2)                 # (nq, m)
    weights /= weights.sum(axis=1, keepdims=True) + 1e-10     # normalize
    return weights @ values                                       # (nq,)

# Usage: smooth a noisy value estimate over 30 sample states
sample_s = np.linspace(0, 10, 30)
values = np.sin(sample_s) * 5 + np.random.randn(30) * 0.5  # noisy
h = silverman_bandwidth(sample_s)
query_s = np.linspace(0, 10, 200)
smooth_v = gaussian_kernel_smooth(query_s, sample_s, values, h)

Chapter 4: Linear Interpolation

Linear interpolation works on a grid: a regular rectangular discretization of the state space. Given a query state s between grid points, the value is computed as a weighted combination of the surrounding grid values, with weights proportional to how close s is to each grid point.

In 1D with grid points at x_l and x_r (the two neighbors of query x), the linear interpolation is:

In 2D, bilinear interpolation uses the four surrounding grid corners. In n dimensions, multilinear interpolation uses 2ⁿ surrounding corners. This exponential scaling is why linear interpolation is practical in low dimensions (1–4) but problematic in higher ones.

Choosing a Grid Resolution

For grid-based interpolation methods, the resolution per dimension is a critical design parameter. Choosing it requires balancing approximation accuracy, memory, and computation:

A practical heuristic: start with a coarse 5×5 (or 5ⁿ) grid, solve approximately, inspect the value function, and double the resolution in dimensions where the function varies most rapidly. ACAS X used a non-uniform grid with finer resolution near zero relative altitude (where collision risk is highest) and coarser resolution at large separations (where the optimal action is always "do nothing").

Linear interpolation is also the foundation of the Q-table with linear interpolation method used in classic RL papers. Sutton & Barto's mountain car problem (1998) used a tile-coded linear approximation that is equivalent to piecewise linear interpolation on a non-uniform grid. The tile coding creates a grid implicitly: each "tile" is a rectangular region in state space with its own value estimate. When multiple tilings overlap, the interpolated value is a weighted sum of tiles, approximating the same multilinear formula as regular grid interpolation. The advantage: tile coding adapts naturally to the state-space geometry (tiles can be offset and rotated) while maintaining the O(1) query time of grid methods.

Despite this limitation, piecewise linear interpolation on a grid is the most widely deployed approximate value function method in safety-critical aviation. The original TCAS (Traffic Collision Avoidance System, deployed in virtually all commercial aircraft since 1993) used a lookup table. ACAS X (2020s) uses a 5D grid with piecewise linear interpolation, achieving far better performance than TCAS while maintaining hard safety guarantees. The key insight: in the aviation application, 5 state dimensions × modest resolution is manageable; the smoothness assumption holds because physics is continuous.

The multilinear interpolation formula in n dimensions:

where λ_j is the fractional offset in dimension j, and corner_b is the grid corner indexed by binary vector b. This formula has 2ⁿ terms and requires 2ⁿ table lookups per query. For n=5 that is 32 lookups; for n=10 that is 1024. Simplex interpolation replaces this with just n+1 = 6 or 11 lookups respectively.

python
import numpy as np

def linear_interp_1d(grid_x, grid_U, query_x):
    """Piecewise linear interpolation in 1D."""
    # Find the two surrounding grid points
    idx = np.searchsorted(grid_x, query_x) - 1
    idx = np.clip(idx, 0, len(grid_x) - 2)
    xl, xr = grid_x[idx], grid_x[idx + 1]
    Ul, Ur = grid_U[idx], grid_U[idx + 1]
    lam = (xr - query_x) / (xr - xl)  # weight for left point
    return lam * Ul + (1 - lam) * Ur

# Example: value function on [0, 1] with 5 grid points
grid_x = np.linspace(0, 1, 5)
grid_U = np.array([0.0, 0.8, 1.2, 0.9, 0.1])
print(linear_interp_1d(grid_x, grid_U, 0.35))  # 0.96

Chapter 5: Simplex Interpolation

Simplex interpolation avoids the exponential 2ⁿ corner problem of multilinear interpolation. Instead of requiring all 2ⁿ corners of a hypercube, it decomposes each grid cell into simplices (triangles in 2D, tetrahedra in 3D, etc.) using only n+1 vertices each.

A simplex in n dimensions has n+1 vertices. Any point inside a simplex can be written as a convex combination of the vertices with non-negative barycentric coordinates λ₀, ..., λ_n that sum to 1:

The interpolated value is then U_θ(s) = ∑_i λ_i U(s_i). This uses only n+1 lookups per query instead of 2ⁿ. The tradeoff: the interpolation is linear within each simplex but may have discontinuous derivatives at simplex boundaries.

Worked Numerical Example: 2D Simplex Interpolation

Suppose our 2D grid has spacing 1 in each dimension. Grid corners at (0,0), (1,0), (0,1), (1,1) with values U(0,0)=2, U(1,0)=5, U(0,1)=3, U(1,1)=8. We query at s=(0.7, 0.4).

Step 1: Fractional coordinates. Already in [0,1]: d_x=0.7, d_y=0.4.

Step 2: Determine triangle. Since d_x = 0.7 ≥ d_y = 0.4, we're in the lower triangle with vertices (0,0), (1,0), (1,1).

Step 3: Barycentric coordinates. λ₀ = 1 − d_x = 0.3, λ₁ = d_x − d_y = 0.3, λ₂ = d_y = 0.4. Verify: 0.3 + 0.3 + 0.4 = 1.0 ✓.

Step 4: Interpolate. U(0.7, 0.4) = 0.3×U(0,0) + 0.3×U(1,0) + 0.4×U(1,1) = 0.3×2 + 0.3×5 + 0.4×8 = 0.6 + 1.5 + 3.2 = 5.3.

Compare: bilinear interpolation at (0.7, 0.4) uses all 4 corners. λ_x=0.7, λ_y=0.4: U = (1−0.7)(1−0.4)×2 + 0.7(1−0.4)×5 + (1−0.7)0.4×3 + 0.7×0.4×8 = 0.36 + 2.1 + 0.36 + 2.24 = 5.06. The two methods agree closely but use different numbers of lookups (3 vs 4). In 10D: simplex uses 11 lookups vs bilinear's 1024.

Here is the full simplex interpolation algorithm for a 2D grid with state (x, y) ∈ [0,1]²:

python
import numpy as np

def simplex_interp_2d(grid_U, query_x, query_y, x_min, x_max, y_min, y_max):
    """
    Simplex interpolation in 2D.
    grid_U: 2D array of values at grid points
    query: (x, y) continuous state
    Returns interpolated value.
    """
    nx, ny = grid_U.shape
    # Map to fractional grid coordinates [0, nx-1] x [0, ny-1]
    fx = (query_x - x_min) / (x_max - x_min) * (nx - 1)
    fy = (query_y - y_min) / (y_max - y_min) * (ny - 1)
    # Integer and fractional parts
    ix, iy = int(np.clip(fx, 0, nx-2)), int(np.clip(fy, 0, ny-2))
    dx, dy = fx - ix, fy - iy  # fractional parts in [0, 1)
    # Sort fractional parts: determines which triangle we're in
    if dx >= dy:  # lower triangle: (0,0)-(1,0)-(1,1)
        lam = [1 - dx, dx - dy, dy]  # barycentric coordinates
        pts = [(ix, iy), (ix+1, iy), (ix+1, iy+1)]
    else:         # upper triangle: (0,0)-(0,1)-(1,1)
        lam = [1 - dy, dy - dx, dx]
        pts = [(ix, iy), (ix, iy+1), (ix+1, iy+1)]
    return sum(lam[k] * grid_U[pts[k]] for k in range(3))

# Only 3 lookups (n+1 = 3 for 2D) vs 4 for bilinear (2^2 = 4)

Chapter 6: Linear Regression

Linear regression approximation represents the value function as a linear combination of basis functions β₁(s), ..., β_p(s):

The parameters θ are found by minimizing the mean squared error over sample states. This has a closed-form solution: the ordinary least squares estimator θ^* = (B^TB)⁻¹B^Ty, where B is the matrix with rows β(s_i)^T and y is the vector of target values y_i.

Basis function choices for value function approximation:

• Polynomials: β_k(s) = s₁^a₁s₂^a₂... (monomials up to some degree)
• Radial basis functions: β_k(s) = exp(−||s − c_k||²/(2σ²)) (Gaussians centered at c_k)
• Fourier features: sin/cos of random projections (useful for periodic value functions)
• Tile coding: binary features from overlapping tilings (a classic RL trick)

The OLS estimator θ^* = (B^TB)⁻¹B^Ty requires matrix inversion. For p basis functions and m sample states, the matrix B^TB is p×p. Inverting it costs O(p³). For p=100 basis functions this is fast; for p=10,000 (e.g., with many tile coding tilings) it becomes expensive. In practice, use numpy.linalg.lstsq or iterative solvers (conjugate gradient) for large p. Alternatively, use stochastic gradient descent to update θ online rather than solving the full OLS every iteration — this is the basis of TD-learning with linear function approximation, analyzed by Sutton (1988) and Tsitsiklis & Van Roy (1997).

Tile coding is a classic trick from RL (Sutton & Barto) that creates binary features by overlaying multiple offset grids ("tilings") over the state space. Each tiling is like a low-resolution grid; the feature for a state is 1 in the bin it falls into on each tiling, 0 everywhere else. Using k overlapping tilings at different offsets gives k features per state. The advantage: linear regression on tile-coded features can represent nonlinear functions, while remaining fully linear in parameters.

Chapter 7: Neural Network Regression

A neural network approximator U_θ(s) passes the state through a stack of learned transformations. Each layer computes a linear combination of the previous layer's activations, followed by a nonlinear activation function:

The parameters θ = {W₁, b₁, W₂, b₂, ..., W_out, b_out} are learned by stochastic gradient descent on the MSE loss. Common activation functions: ReLU σ(x) = max(0,x), tanh, or sigmoid.

Neural networks are universal approximators: with enough width or depth, a ReLU network can approximate any continuous function to arbitrary precision. In practice, they are the default for high-dimensional state spaces — they power DQN, A3C, PPO, and virtually all modern deep RL.

For the approximate value iteration setting (offline, batch), neural networks are trained by minimizing the mean squared Bellman error over sample transitions (s, a, r, s'). The key data pipeline is: (1) sample a batch of transitions from a replay buffer, (2) compute Bellman targets y_i = r_i + γ · max_a' U_θ−(s'_i, a') using the frozen target network, (3) gradient step on ∑(U_θ(s_i, a_i) − y_i)². The target network weights θ

Architecture Design for Value Function Networks

The architecture of the value network matters enormously. Key design decisions:

For the MDP in the showcase (1D continuous, 2 actions), a minimal network is 2 hidden layers × 32 units. For Atari (image input, 18 actions), DQN used 3 convolutional layers followed by 2 fully connected layers with 512 units. For robotics with 50D state spaces, typical architectures are 3 fully connected layers with 256 units each.

python
import torch
import torch.nn as nn

class ValueNetwork(nn.Module):
    def __init__(self, state_dim, hidden=64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(state_dim, hidden), nn.ReLU(),
            nn.Linear(hidden, hidden),  nn.ReLU(),
            nn.Linear(hidden, 1)            # scalar value output
        )
    def forward(self, s):
        return self.net(s).squeeze(-1)  # (batch,) scalar per state

# Approximate value iteration step
def approx_value_iter_step(net, target_net, optimizer, states, rewards, next_states, gamma=0.99):
    with torch.no_grad():
        targets = rewards + gamma * target_net(next_states)  # Bellman backup
    preds = net(states)
    loss = nn.functional.mse_loss(preds, targets)
    optimizer.zero_grad(); loss.backward(); optimizer.step()
    return loss.item()

Why Neural Networks Are Unstable in Value Iteration

The fundamental problem: the Bellman backup computes y_i = R + γ · max_a U_θ(s'). But U_θ uses the current network weights θ. So when we take a gradient step to make U_θ(s) closer to y_i, we also change U_θ(s'), which changes y_i on the next iteration. We are chasing a moving target. This can cause oscillations or divergence.

DQN (Mnih et al., 2015, DeepMind) introduced two key fixes that made neural value functions stable for the first time:

Together these two fixes turned a theoretically unstable algorithm into a practical powerhouse. DQN achieved superhuman performance on 49 Atari games using this architecture: a CNN for image features, two hidden layers with 512 ReLU units, output of Q-values for each action, trained with replay buffer of 1M transitions and target network update every 10,000 steps.

Data Flow in Neural Approximate Value Iteration

Double DQN and Other Improvements

Since the original DQN paper (2015), several improvements have become standard. For the approximate value iteration context:

All of these improvements operate on the same approximate value iteration loop from the textbook. The core insight is unchanged: sample, backup, fit. The improvements address specific failure modes of the basic loop when scaled to complex environments. For the textbook's simple MDP examples, the basic DQN loop without improvements is sufficient and often preferred for pedagogical clarity.

Concretely, the tensor shapes for a 1D continuous state problem with 2 actions:

The forward pass: s → hidden layers → Q-values (batch, 2) → gather Q(s,a) at chosen action → compare with target y. Gradient flows back through gather → Q-head → hidden layers only (not through the target network, which is detached from the computation graph).

Chapter 8: Showcase — Approximate Value Iteration

Here is the full approximate value iteration loop, visualized in real time on a 1D continuous MDP. The state is a position s ∈ [0, 10]. The reward is high at the goal (s=8) and zero elsewhere. The discount is γ=0.9. There is no analytical solution — we must approximate.

You choose the approximator: k-NN, kernel smoothing, or linear regression with RBF features. The loop samples states randomly, computes Bellman backups, refits the approximator, and shows how the value function evolves. Watch the value function grow toward the goal.

From Approximate Value Function to Policy Extraction

Once we have a good approximation U_θ, extracting a policy requires evaluating the Bellman maximization at each decision point:

For discrete action spaces (e.g., left/right/stop), this is a finite argmax — cheap. For continuous actions, this becomes a non-convex optimization problem. Common approaches:

In the showcase's 1D MDP with 2 discrete actions (left/right), policy extraction is O(2). For ACAS X with 5 actions (alerts), O(5). The simplicity of discrete, small action spaces is one reason they are preferred in safety-critical applications.

Chapter 9: Connections & What's Next

Convergence Theory for Approximate Value Iteration

Exact value iteration converges because the Bellman operator T is a γ-contraction: ||TU − TU'||_∞ ≤ γ ||U − U'||_∞. The fixed point of T is U*, and repeated application of T shrinks any error by γ per iteration.

With approximate value iteration, the error does not shrink to zero. Each iteration first applies T (contraction) then applies the approximation operator Π (projection onto the function class). The composition ΠT is no longer a contraction in general. Three cases:

The error bound ε/(1−γ) has two terms: ε is the best-in-class approximation error (how well the function family can represent U*), and 1/(1−γ) is the discount amplification factor. As γ → 1, even small approximation errors get amplified arbitrarily. This is why deep RL with γ=0.99 is so finicky: the 1/(1−0.99) = 100 amplification magnifies any approximation error by 100x in the value bound.

On-Policy vs Off-Policy Sampling

The choice of sample states is closely related to the concept of model bias in approximate dynamic programming. If the sample states are not representative of the states visited by the optimal policy, the approximation will be accurate in the wrong region of state space, leading to a suboptimal policy. This is a form of distribution mismatch: we minimize loss on the training distribution (sample states) but the policy deploys in a different distribution (states the optimal policy visits). The solution is iterative: use the current approximate policy to generate new sample states, refit, and repeat. This is called approximate policy iteration (API) and is the algorithmic framework underlying SARSA, on-policy actor-critic, and TRPO in modern deep RL.

The sampling distribution for sample states is a crucial design choice that the textbook addresses explicitly. Two main strategies:

For the grid-based methods (linear interpolation, simplex), a uniform grid is natural — the grid points are the sample states. For k-NN and kernel methods, on-policy sampling often works best after the first few iterations. For neural networks, experience replay (off-policy) is standard: it reduces correlation between consecutive samples and allows reuse of data, which is critical when environment interaction is expensive.

The textbook also notes an important practical consideration: the sampling distribution of states in S matters enormously. If your sample states cluster in low-value regions but the optimal policy mostly visits high-value regions, the approximation will be poor where it matters most. Common solutions: (1) follow the current approximate policy to generate samples (on-policy sampling), (2) importance-weight samples from a fixed distribution, or (3) use a uniform grid if the state space is bounded and low-dimensional. The interplay between sample distribution and approximation error is the core technical challenge of approximate dynamic programming and motivates most of the innovations in deep RL (experience replay, on-policy methods, distributional RL) that have appeared since the DQN paper.

Further Reading

The textbook's Chapter 8 is a strong survey. For deeper coverage:

• Sutton & Barto, "Reinforcement Learning: An Introduction" (2nd ed., 2018), Chs. 9–12. The canonical treatment of function approximation in RL: linear methods, neural networks, eligibility traces.
• Mnih et al. (2015), "Human-level control through deep reinforcement learning." The DQN paper. Shows fitted value iteration with neural networks at scale for the first time.
• Ernst, Geurts & Wehenkel (2005), "Tree-Based Batch Mode Reinforcement Learning." Original fitted value iteration paper; shows any regression oracle can be used.
• Kochenderfer et al., "Algorithms for Decision Making" Ch. 8. The textbook chapter; most comprehensive single-chapter survey of approximate value function methods.

Computing Approximate Bellman Backups in Continuous State Spaces

The Bellman backup y_i = max_a[R(s_i,a) + γ ∑_s' T(s'|s_i,a) U_θ(s')] requires integrating U_θ(s') over the transition distribution T. For deterministic transitions (s' = f(s,a) exactly), this is trivial: y_i = max_a[R(s_i,a) + γ U_θ(f(s_i,a))]. For stochastic transitions (s' ~ T(·|s,a)), we must approximate the expectation:

• Monte Carlo: Sample K next states s'_k ~ T(·|s,a); approximate E[U(s')] ≈ (1/K) ∑_k U_θ(s'_k). Simple but requires K environment rollouts per backup.
• Gaussian quadrature: If T is Gaussian, use exact quadrature points and weights. O(Kⁿ) points for n-dimensional state but exact for polynomial approximators up to degree 2K−1.
• Reparameterization: Write s' = g(s,a,ε) for ε ~ N(0,I); estimate E[U(s')] ≈ U(g(s,a,E[ε])) (deterministic equivalent, exact for linear U_θ).

For the showcase's 1D MDP, the transitions are deterministic (s' = clamp(s + action, 0, 10)), so backups are exact single-step lookups. In the real-world ACAS X problem, the aircraft dynamics are stochastic (due to wind and sensor noise), and Monte Carlo backup with K=10–50 samples per state per action is the standard approach.

Real-World Deployment: ACAS X as a Case Study

The Airborne Collision Avoidance System X (ACAS X) is the premier real-world deployment of approximate value function methods. Developed at MIT Lincoln Laboratory and the FAA, it uses a 5D piecewise linear approximation over the states: relative altitude, relative bearing, intruder closure rate, own-ship altitude rate, and intruder altitude rate.

Key engineering choices and why they were made:

The key lesson: for safety-critical aviation systems, neural networks are not acceptable — regulators require deterministic behavior that can be formally verified. This does not mean neural networks are never used in aviation: they are used for perception (camera-based object detection), speech recognition, and anomaly detection in maintenance. But for the decision-making layer that determines pilot advisories, only certifiable methods (look-up tables, piecewise linear functions) are acceptable under current FAA standards. The push toward neural network verification (using formal methods like Marabou or α,β-CROWN) may change this in the future.

The ACAS X value table is publicly released and has been analyzed extensively by the formal verification community. Katz et al. (2017) verified properties of small ACAS X networks using the Reluplex satisfiability checker. Julian et al. (2019) showed how to compress the 5D grid into a neural network that is 1000x smaller, while Kouvaros et al. (2021) formally verified certain safety properties of those compressed networks. This research direction — compressing safe grid-based policies into neural networks and then formally verifying the neural version — represents the state of the art in bridging the certifiability gap between classical and neural approaches. The textbook's approximate value function methods are the starting point for this entire research thread.

How This Chapter Connects to the Rest of the Book

Chapter 8's approximate value functions feed directly into every subsequent chapter. Chapter 9 (Online Planning) replaces offline precomputation with real-time search using rollout policies — the rollout policy is exactly the policy extracted from an approximate value function via the greedy policy improvement step. Chapter 10 (Policy Search) parameterizes π_θ(a|s) directly and optimizes θ with gradient ascent; approximate value functions appear here as critics in actor-critic methods, providing low-variance gradient estimates. Chapters 12–17 (Model-Free RL) learn U_θ or Q_θ without knowing T or R — but the Bellman operator being approximated is identical; the only difference is whether samples come from a simulator or from real environment interactions.

The critical design axis that runs through all of these chapters is the approximation-exploration tradeoff: the more compact U_θ is (fewer parameters, simpler basis), the faster it can be fit and the less data it needs — but the worse it approximates the true U* in complex regions of state space. The practitioner's job is to choose an approximator that is rich enough to capture the value structure that matters for decision-making, while being simple enough to be trained reliably from the available data. Chapter 8 gives you the vocabulary and toolkit to make this choice deliberately rather than by default.

Generalization and the Credit Assignment Problem

One of the deepest challenges in approximate value functions is credit assignment: when the agent takes action a in state s and eventually receives a reward r, how much of r is attributable to a versus to later actions? The Bellman equation solves this by definition — Q(s,a) = R(s,a) + γ max_a' Q(s',a') assigns exactly the right credit to action a by recursion. But with function approximation, each gradient update that adjusts θ to reduce the Bellman error at one (s,a) pair also changes U_θ(s') for all other states s' where the basis functions overlap. This generalization is both the power and the hazard of function approximation: a good approximator generalizes correctly (nearby states have similar values), while a bad approximator corrupts good value estimates when it over-generalizes.

Tile coding (Section 8.3) manages this tradeoff explicitly: the tiles are designed so that generalization happens only within each tile, never across tiles. This gives exact local estimates at the cost of coarse approximation at boundaries. Neural networks generalize more smoothly but without explicit control. The field of representation learning (Bengio et al., 2013; Lesort et al., 2018) studies which learned state representations support accurate value generalization — this is the same as asking what features make a good basis for U_θ.