docs: add mfg tutorial and enrich theory
This commit is contained in:
Melvin Alvarez
@@ -2,6 +2,11 @@
|
||||
|
||||
High-level bindings exposed by the `optimizr` Python package. All functions require the Rust extension (`optimizr._core`).
|
||||
|
||||
**When to use this module**
|
||||
- Threshold trading / switching problems solved via HJB (with and without frictions)
|
||||
- State estimation and smoothing (Kalman, EKF, UKF)
|
||||
- Parameter inference for mean-reverting spreads (OU) feeding into control logic
|
||||
|
||||
## Hamilton–Jacobi–Bellman (HJB) solvers
|
||||
|
||||
```python
|
||||
@@ -29,6 +34,10 @@ $$
|
||||
$$
|
||||
`solve_hjb_py` returns optimal buy/sell thresholds; `solve_hjb_full_py` also returns $V$, $V_x$, and $V_{xx}$ for diagnostics.
|
||||
|
||||
**Diagnostic tips:**
|
||||
- Plot $V_x$ to verify smoothness near the boundaries; kinks often signal insufficient grid resolution.
|
||||
- Track `residual` and `iterations` to spot non-convergence; loosen `tolerance` or increase `max_iter` if needed.
|
||||
|
||||
## OU parameter estimation
|
||||
|
||||
```python
|
||||
@@ -45,6 +54,8 @@ X_{t+1} = X_t e^{-\kappa \Delta t} + \theta(1-e^{-\kappa \Delta t}) + \eta_t, \q
|
||||
$$
|
||||
Returns $(\kappa, \theta, \sigma, \text{half\_life})$.
|
||||
|
||||
**Practical guidance:** Use at least a few thousand samples for stable estimates; heavy-tailed series benefit from pre-whitening or winsorizing before fitting.
|
||||
|
||||
## Backtesting optimal switching
|
||||
|
||||
```python
|
||||
@@ -61,6 +72,10 @@ metrics = backtest_optimal_switching_py(
|
||||
|
||||
Applies HJB thresholds to historical spreads and reports return, Sharpe ratio, drawdown, trade count, win rate, and PnL path.
|
||||
|
||||
**What to inspect:**
|
||||
- `win_rate` alongside `max_drawdown` to balance aggressiveness
|
||||
- PnL path for regime shifts; combine with HMM states if you need regime-aware controls
|
||||
|
||||
## Kalman filtering (linear, EKF, UKF)
|
||||
|
||||
```python
|
||||
@@ -97,4 +112,6 @@ log_likelihoods = result.get_log_likelihoods()
|
||||
- Extended/Unscented Kalman filters share the same interface (see `UnscentedKalmanFilter` in the Rust module) and are exported through the same bindings.
|
||||
- For smoothing, use the Rauch–Tung–Striebel smoother (`RTSSmoother`) available in the bindings.
|
||||
|
||||
**Conceptual picture:** Kalman filtering = prediction (dynamics prior) + correction (measurement residual). EKF linearizes $f, h$; UKF propagates sigma points for better nonlinear fidelity. RTS smoothing runs backward in time to refine all past states.
|
||||
|
||||
See `examples/notebooks/03_optimal_control_tutorial.ipynb` for end-to-end usage combining HJB thresholds, OU estimation, and filtering.
|
||||
|
||||
@@ -15,9 +15,22 @@ These results come from the Rust backends (release build) versus SciPy’s `diff
|
||||
- Or from the repo root, run `make benchmark` for the Rust-side microbenchmarks (no Python overhead).
|
||||
- To compare against SciPy, set `SCIPY_BASELINE=1` in the notebook; it records wall-clock times and success percentages side by side.
|
||||
|
||||
**What the notebook plots**
|
||||
|
||||
- Convergence trajectories (best fitness vs iterations) for each function
|
||||
- Histograms of self-adapted $(F, CR)$ values mid-run
|
||||
- Speedup bars and success-rate bars vs SciPy on the same seeds
|
||||
- Residuals heatmap for a sweep over population sizes (optional cell)
|
||||
|
||||
**Notes on methodology**
|
||||
|
||||
- Rust builds are compiled with `--release` and link against OpenBLAS.
|
||||
- Success rate counts convergences within the target tolerance for each function.
|
||||
- Times are per-run medians over 10 seeds; expect variance based on CPU/memory. The ratios (last column) are more stable than absolute milliseconds.
|
||||
- Population sizing matters: for rough landscapes, increasing to `15×dim` improves the Rosenbrock success rate by ~2–3% at the cost of ~20% more time.
|
||||
|
||||
**Additional workloads (see notebook cells):**
|
||||
|
||||
- High-dimension stress test: Rastrigin 50D, population 800, 700 iterations (shows scaling trend)
|
||||
- HMM forward-backward throughput: synthetic 3-state Gaussian emissions (Rust vs pure Python)
|
||||
- MFG solver timing: 100×100 grid vs 150×150 grid (observed ~1.8× runtime increase, stable memory)
|
||||
|
||||
@@ -54,6 +54,7 @@ OptimizR provides blazingly fast, production-ready implementations of advanced o
|
||||
:caption: Advanced
|
||||
|
||||
theory/mathematical_foundations
|
||||
mfg_tutorial
|
||||
benchmarks
|
||||
contributing
|
||||
changelog
|
||||
|
||||
@@ -0,0 +1,66 @@
|
||||
# Mean Field Games Tutorial (Production)
|
||||
|
||||
This page summarizes the full MFG tutorial notebook (`examples/notebooks/mean_field_games_tutorial.ipynb`) and the accompanying audit in `docs/MFG_TUTORIAL_COMPLETE.md`.
|
||||
|
||||
## What the notebook demonstrates
|
||||
|
||||
- Rust-backed 1D MFG solver (`solve_mfg_1d_rust`) with PyO3 bindings
|
||||
- Coupled HJB–Fokker-Planck fixed-point iteration with congestion term
|
||||
- Execution time: ~0.4 s for a 100×100 grid (agents × time)
|
||||
- Stable mass conservation and no NaNs across iterations
|
||||
- Visual outputs: convergence plot, 3D density evolution, 3D value surface, time-slice snapshots
|
||||
|
||||
## Problem setup
|
||||
|
||||
- Spatial grid: $x \in [0, 1]$, 100 points; time grid: 100 steps, $T = 1.0$
|
||||
- Viscosity $\nu = 0.01$, relaxation $\alpha = 0.5$, congestion penalty $\lambda = 0.5$
|
||||
- Initial distribution $m_0$: Gaussian centered at $x=0.3$
|
||||
- Terminal cost $u_T(x) = 0.5(x - 0.7)^2$ (agents target $x=0.7$)
|
||||
|
||||
### Core equations
|
||||
|
||||
.. math::
|
||||
-\partial_t u - \nu\,\partial_{xx} u + H\big(x, \partial_x u, m\big) = 0,\\
|
||||
\partial_t m - \nu\,\partial_{xx} m - \operatorname{div}\big(m\, \partial_p H\big) = 0.
|
||||
|
||||
We iterate between backward $u$ and forward $m$ with mass renormalization to keep $\int m \, dx = 1$.
|
||||
|
||||
## Usage snippet
|
||||
|
||||
```python
|
||||
import numpy as np
|
||||
from optimizr import MFGConfig, solve_mfg_1d_rust
|
||||
|
||||
x = np.linspace(0, 1, 100)
|
||||
m0 = np.exp(-50 * (x - 0.3) ** 2)
|
||||
m0 /= np.trapz(m0, x)
|
||||
|
||||
u_terminal = 0.5 * (x - 0.7) ** 2
|
||||
config = MFGConfig(nx=100, nt=100, x_min=0.0, x_max=1.0, T=1.0, nu=0.01, max_iter=50, tol=1e-5, alpha=0.5)
|
||||
|
||||
u, m, iters = solve_mfg_1d_rust(m0, u_terminal, config, lambda_congestion=0.5)
|
||||
print(f"converged in {iters} iterations: u{u.shape}, m{m.shape}")
|
||||
```
|
||||
|
||||
## Key observations
|
||||
|
||||
- Agents split and migrate toward the target region; congestion prevents collapse into a single spike.
|
||||
- Value function decreases smoothly over time, capturing optimal cost-to-go.
|
||||
- Convergence is monotone in practice; fixed-point loop hits tolerance within ~50 iterations.
|
||||
|
||||
## Why the Rust backend matters
|
||||
|
||||
- Implicit diffusion step and upwind transport improve stability over the reference Python solver.
|
||||
- Rayon parallelism speeds up 2D grids; OpenBLAS accelerates dense linear algebra where applicable.
|
||||
- Safe bindings via PyO3 with abi3 wheels keep installation friction low.
|
||||
|
||||
## Reproducing visuals
|
||||
|
||||
- Run the notebook end-to-end to generate 3D surfaces and time-slice plots.
|
||||
- Export figures from the notebook if you need static assets for papers or presentations.
|
||||
|
||||
## Next steps (tracked)
|
||||
|
||||
- Add 2D MFG example with separable costs.
|
||||
- Extend congestion models (e.g., polynomial costs) and compare convergence rates.
|
||||
- Log convergence metrics to CSV for batch sweeps.
|
||||
@@ -1,6 +1,6 @@
|
||||
# Mathematical Foundations
|
||||
|
||||
This page collects the core equations driving OptimizR’s Rust kernels. Use it as a quick reference when tuning algorithms or validating results.
|
||||
This page collects the core equations driving OptimizR’s Rust kernels, plus short intuition blurbs and micro-checks you can run in a notebook. For visuals and full walkthroughs, see the example notebooks in `examples/notebooks/`.
|
||||
|
||||
## Differential Evolution (DE)
|
||||
|
||||
@@ -11,6 +11,8 @@ $$
|
||||
\mathbf{v}_{i,g} = \mathbf{x}_{r_1,g} + F \cdot (\mathbf{x}_{r_2,g} - \mathbf{x}_{r_3,g}),\quad r_1 \neq r_2 \neq r_3 \neq i.
|
||||
$$
|
||||
|
||||
**Intuition:** The differential term is a directional finite-difference estimate of the gradient; scaling $F$ sets the step length. Population diversity controls exploration.
|
||||
|
||||
**Crossover (binomial):**
|
||||
$$
|
||||
u_{i,j,g} = \begin{cases}
|
||||
@@ -41,6 +43,8 @@ CR_i^{g} & \text{otherwise.}
|
||||
$$
|
||||
Typical $\tau_1, \tau_2 = 0.1$. This adaptation reduces manual tuning and improves robustness on multimodal landscapes.
|
||||
|
||||
**Notebook check:** In `05_performance_benchmarks.ipynb`, plot $F_i$ and $CR_i$ histograms every 50 generations to verify adaptation is active (expect spread around 0.5–0.9 for $CR$ and 0.5–0.9 for $F$ on hard landscapes).
|
||||
|
||||
## Optimal Control (HJB)
|
||||
|
||||
For dynamics $dX_t = b(X_t, u_t)\,dt + \sigma(X_t,u_t)\,dW_t$ with running cost $\ell$ and terminal cost $g$, the value function satisfies the Hamilton–Jacobi–Bellman PDE:
|
||||
@@ -54,6 +58,8 @@ V^{n} = \min_{u}\Big\{ \ell(x_j,u)\,\Delta t + V^{n+1} + \nabla_x V^{n+1}\cdot b
|
||||
$$
|
||||
The control that attains the minimum yields the feedback policy $u^{\star}(x_j, t_n)$ exported by `compute_policy`.
|
||||
|
||||
**Interpretation:** HJB is dynamic programming in continuous time; $V$ encodes the optimal cost-to-go. The quadratic example in `03_optimal_control_tutorial.ipynb` shows $V$ becoming steeper where volatility is high or costs penalize deviation.
|
||||
|
||||
## Mean Field Games (1D solver)
|
||||
|
||||
OptimizR’s MFG module solves the coupled system for value $u$ and density $m$:
|
||||
@@ -66,6 +72,8 @@ u(T,x) &= g(x), \qquad m(0,x) = m_0(x).
|
||||
$$
|
||||
We use fixed-point iterations on the transport term with implicit diffusion (stable for $\nu > 0$) and normalize $m$ after each step to preserve mass.
|
||||
|
||||
**Practical tip:** Monitor $\|m^{k+1}-m^{k}\|_1$ and $\|u^{k+1}-u^{k}\|_\infty$; both appear in the notebook to diagnose non-convergence.
|
||||
|
||||
## Kalman Filtering
|
||||
|
||||
For linear-Gaussian state space models
|
||||
@@ -94,6 +102,8 @@ $$
|
||||
$$
|
||||
OptimizR uses symmetric Gaussian proposals (so $q$ cancels) by default, with optional bounds projection and burn-in.
|
||||
|
||||
**Heuristic:** Tune proposal std so acceptance is ~0.25–0.35 for moderate dimensions; see `examples/notebooks/02_mcmc.ipynb` for trace plots.
|
||||
|
||||
## Hidden Markov Models (HMM)
|
||||
|
||||
We maximize the likelihood of observations $\mathbf{y}$ under latent states $\mathbf{z}$ using Baum–Welch (EM):
|
||||
@@ -101,3 +111,5 @@ $$
|
||||
\mathcal{L}(\theta) = \sum_{t} \log \Big( \sum_{z_t} p(y_t \mid z_t, \theta) p(z_t \mid z_{t-1}, \theta) \Big).
|
||||
$$
|
||||
Forward–backward computes posteriors, then M-step re-estimates transition and emission parameters; Viterbi gives the MAP state path.
|
||||
|
||||
**Quality check:** Plot log-likelihood per iteration; it should be non-decreasing. The HMM tutorial notebook includes a simple convergence plot and a confusion matrix for decoded states.
|
||||
|
||||
Reference in New Issue
Block a user