Optimiz-rs 🚀

High-performance optimization algorithms in Rust with Python bindings

Optimiz-rs provides blazingly fast, production-ready implementations of advanced optimization and statistical inference algorithms. Built with Rust for maximum performance and exposed to Python through PyO3, it delivers up to 86× speedup over pure-Python references on intrinsically loopy / sequential workloads.

800-particle mean-reverting McKean–Vlasov flow simulated by optimizr.mean_reverting_mckean_vlasov — two clouds at x = ±2 fuse under dX_t = θ(m̄_t − X_t) dt + σ dW_t. Source: examples/animate_mckean_vlasov.py.

✨ What's New in v2.0.0

v2 ships eight brand-new CPU-only generic numerical primitive groups with full Python bindings, on top of every v1.x algorithm (which remain available). The Python module name is unchanged: import optimizr as opt.

Rough volatility & integral equations

• solve_fractional_ode(h0, alpha, t_horizon, n_steps, rhs) — Caputo fractional ODE Adams scheme.
• solve_volterra(g, kernel, t_horizon, n_steps) — second-kind Volterra integral equation by trapezoidal product integration.
• geometric_grid_lift(kernel, t_samples, n_factors, gamma_min, gamma_max) — multi-exponential approximation of a kernel by NNLS on a geometric rate grid (Markovian lift à la Abi Jaber–El Euch).
• fourier_invert(char_fn, t_grid, x_grid) — characteristic-function → density inversion (Carr–Madan style).
• mittag_leffler_py(z, alpha, beta) — generalised Mittag-Leffler reference function.

Backward SDEs & PDEs

• linear_bsde_constant_coeffs(a, b, c, terminal, n_steps, t_horizon, theta=0.5) — backward SDE θ-scheme (closed-form analytic test against dY = -ρ Y dt).
• fokker_planck_constant(...) — 1-D forward Fokker–Planck solver with conservative central differences.
• hjb_quadratic_2d(...) — explicit upwind solver for 2-D HJB on a Cartesian grid.
• poisson_2d_zero_boundary(...) — 2-D Poisson −Δu = f SOR solver.

Stochastic & quadratic-impact control

• optimal_switching_dp(...) — discrete-time optimal switching by dynamic programming.
• pontryagin_lqr(...) — Pontryagin maximum principle for LQ control.
• two_sided_intensities(...) — bilateral intensity-controlled jump process.
• quadratic_impact_control_py(...) — convex quadratic-cost control on a controlled SDE.

Mean-field & agent-based dynamics

• mean_reverting_mckean_vlasov(initial, theta, sigma, n_steps, t_horizon, seed) — N-particle McKean–Vlasov simulator (returns paths_flat, n_particles, n_steps, time_grid).
• consensus_dynamics(...) — synchronous opinion-dynamics consensus on a graph.
• solve_mfg_1d_rust(MFGConfig) — 1-D mean-field game (HJB ↔ Fokker–Planck fixed-point).

Propagation of chaos

For an interacting N-particle system

dX^{i,N}_t \;=\; b\!\bigl(X^{i,N}_t,\; \mu^N_t\bigr)\, dt \;+\; \sigma\, dW^i_t,
\qquad
\mu^N_t \;=\; \frac{1}{N}\sum_{j=1}^{N}\delta_{X^{j,N}_t},

Sznitman's theorem (1991) states that whenever b is Lipschitz in both arguments, the empirical measure \mu^N_t converges in Wasserstein-2 to the law \mu_t of the McKean–Vlasov limit at rate \mathcal{O}(1/\sqrt{N}), and any finite $k$-tuple of particles becomes asymptotically independent — chaos propagates from t=0 to all later times:

\operatorname{Law}\!\bigl(X^{1,N}_t,\dots,X^{k,N}_t\bigr) \;\xrightarrow[N\to\infty]{w}\; \mu_t^{\otimes k}.

The animation above runs four parallel simulations with N\in\{20,100,500,4000\} sharing the same bimodal initial law, the same drift -\theta(x-\bar{x}) and the same noise \sigma\, dW. The bottom panel tracks the Wasserstein-2 distance W_2(\mu^N_t, \mu_t) to a high-resolution reference and visibly decays as 1/\sqrt{N}. Source: examples/animate_propagation_of_chaos.py. Companion notebook: examples/notebooks/14_mckean_vlasov.ipynb.

Topology, graphs & path signatures

• vietoris_rips_filtration, persistent_homology, bottleneck_distance — TDA primitives.
• combinatorial_laplacian_py, normalised_laplacian_py, random_walk_laplacian_py, spectral_cluster_py — graph spectral analysis.
• path_signature, path_log_signature, random_signature, signature_kernel, shuffle_product, concatenate_signatures — Chen–Strichartz iterated integrals and signature kernels.

Risk, robust inference & calibration

• historical_var_py(losses, alpha), parametric_var_py(...), cvar_value_py(...), minimize_cvar_py(...) — coherent risk measures.
• robust_drift(...), estimate_hurst(...), scale_dependent_hurst(...) — robust drift / Hurst estimation.
• mmd_gaussian(...), f_alpha_lambda_py(...) — generative-calibration hooks (MMD, fractional kernels).

Point processes & Kalman filtering

• simulate_hawkes, simulate_bivariate_hawkes, simulate_fbm, simulate_mixed_fbm — order-flow simulators.
• LinearKalmanFilter, UnscentedKalmanFilter, RTSSmoother — state-space inference.

Quality bar

• 20-test analytic non-regression suite for the v2 public API: tests/test_v2_api.py.
• ABI3 wheels, Python ≥ 3.8.
• Crate name: optimiz-rs (Rust); distribution name: optimiz-rs (PyPI); module name: optimizr (Python import).

pip install --upgrade optimiz-rs
python -c "import optimizr; print(optimizr.__version__)"   # 2.0.0

v2 benchmark (single-threaded, best-of-3, Apple M2)

Generated by examples/benchmark_v2.py:

Workload	Pure Python / NumPy	optimiz-rs (Rust)	Speedup
HMM Baum-Welch (2 states, 5 000 obs, 10 iters)	970.94 ms	14.34 ms	67.7×
Differential evolution (Rastrigin d=5, 50 iters × 20 pop)	417.45 ms	30.03 ms	13.9×
Path signature (T=300, d=3, depth=3)	11.07 ms	0.99 ms	11.2×
Hawkes process (T=100, μ=1, α=0.6, β=1.2)	2.75 ms	0.83 ms	3.3×
MCMC random-walk MH (5 000 samples, d=2)	35.41 ms	20.24 ms	1.7×

Workloads that are fully vectorisable in NumPy (e.g. drift updates for an N-particle SDE without callback) are not in this table: a tight NumPy loop on contiguous arrays is hard to beat from Rust through a PyO3 callback boundary. Use optimiz-rs for the algorithms above and numpy for the rest — both are first-class citizens.

✨ What's New in v1.1.0

This release adds a broad collection of CPU-only generic numerical primitives, all purely additive:

• optimal_control::matrix_riccati — backward RK4 solver for the matrix Riccati ODE.
• timeseries_utils::nonsync_covariance — Hayashi--Yoshida asynchronous covariance estimator.
• timeseries_utils::wavelet — DWT and MODWT (Haar, Daubechies 2--10).
• risk_measures — empirical / parametric VaR and CVaR, plus convex CVaR minimisation.
• graph::laplacian + graph::spectral_clustering — combinatorial / normalised / random-walk Laplacians and spectral clustering with Jacobi diagonalisation.
• topology — Vietoris--Rips persistent homology and bottleneck distance.
• volterra — fractional Caputo Adams solver, Markovian lift by NNLS on a geometric grid, second-kind Volterra solver, direct Fourier inversion of characteristic functions.
• signatures — truncated tensor signatures, log-signatures, random-reservoir projection (Cuchiero--Schmocker--Teichmann), Salvi--Cass--Lyons signature kernel, shuffle product / Chen concatenation.

All new modules are exposed via the Rust API only in this release; Python bindings will follow in a subsequent minor release. The previously stable Python API is unchanged.

✨ What's New in v1.0.0

🎉 Production Ready - First stable release with comprehensive documentation
📚 ReadTheDocs - Full documentation at https://optimiz-r.readthedocs.io
🏗️ Published to crates.io - Install with cargo add optimiz-rs
🐍 Published to PyPI - Install with pip install optimiz-rs
🔒 Stable API - Semantic versioning from v1.0.0 forward

Features

✨ Algorithms Included:

• Mean Field Games: 1D MFG solver, HJB-Fokker-Planck coupling, agent population dynamics
• Differential Evolution: 5 strategies (rand/1, best/1, current-to-best/1, rand/2, best/2), adaptive jDE, convergence tracking
• Optimal Control: HJB solvers, regime switching, jump diffusion, MRSJD framework
• Hidden Markov Models: Baum-Welch training, Viterbi decoding, Gaussian emissions
• MCMC Sampling: Metropolis-Hastings, adaptive proposals, Bayesian inference
• Sparse Optimization: Sparse PCA, Box-Tao decomposition, Elastic Net, ADMM
• Risk Metrics: Hurst exponent, half-life estimation, time series analysis
• Information Theory: Mutual information, Shannon entropy, feature selection
• Mathematical Toolkit: Gradient, Hessian, Jacobian, statistics, linear algebra

🚀 Performance:

• 50-100× faster than pure Python implementations
• 95% memory reduction vs NumPy/SciPy
• Parallel-ready with Rayon infrastructure
• Production-tested on multi-dimensional problems

🐍 Python-First API:

• Clean, intuitive NumPy-based interface
• Rich result objects with convergence diagnostics
• Type hints and comprehensive documentation
• Jupyter notebook integration

Installation

From PyPI (Python)

pip install optimiz-rs

The PyPI distribution name is optimiz-rs (with dash). The Python import name is optimizr (no dash): import optimizr as opt.

From crates.io (Rust)

cargo add optimiz-rs

The Rust crate is also optimiz-rs; its library name is optimizr (no dash) — matching the Python module.

From Source

# Clone the repository
git clone https://github.com/ThotDjehuty/optimiz-r.git
cd optimiz-r

# Install with maturin
pip install maturin
maturin develop --release

# Or install in editable mode
pip install -e .

Using Docker

# Start Jupyter notebook server with examples
docker-compose up dev
# Access at http://localhost:8888

# Run all tests
docker-compose run test

# Build distribution wheels
docker-compose run build

Quick Start

Differential Evolution (Enhanced in v0.2.0)

import numpy as np
from optimizr import differential_evolution

# Rosenbrock function (challenging non-convex problem)
def rosenbrock(x):
    return sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

# Optimize with adaptive jDE (self-tuning parameters)
result = differential_evolution(
    objective_fn=rosenbrock,
    bounds=[(-5, 5)] * 10,
    maxiter=1000,
    strategy='best1',  # 5 strategies: rand1, best1, currenttobest1, rand2, best2
    adaptive=True,     # Adaptive F and CR parameters (jDE algorithm)
    atol=1e-6
)

print(f"Optimum: {result.x}")
print(f"Value: {result.fun} (expected: 0.0)")
print(f"Converged: {result.converged}, Iterations: {result.nit}")
# Typical speedup: 74-88× faster than SciPy

Mathematical Toolkit (New in v0.2.0)

from optimizr import maths_toolkit as mt
import numpy as np

# Numerical differentiation
f = lambda x: x[0]**2 + 2*x[1]**2 + x[0]*x[1]
x = np.array([1.0, 2.0])

gradient = mt.gradient(f, x)        # ∇f(x)
hessian = mt.hessian(f, x)          # H(f)(x)
jacobian = mt.jacobian(f, x)        # J(f)(x)

# Statistics
data = np.random.randn(1000)
stats = {
    'mean': mt.mean(data),
    'var': mt.variance(data),
    'std': mt.std_dev(data),
    'skew': mt.skewness(data),
    'kurt': mt.kurtosis(data)
}

# Linear algebra
A = np.random.randn(5, 5)
norm_l1 = mt.norm_l1(A)
norm_l2 = mt.norm_l2(A)
A_norm = mt.normalize(A)

Mean Field Games (New in v0.3.0)

from optimizr import MFGConfig, solve_mfg_1d_rust
import numpy as np

# Configure MFG problem for population dynamics
config = MFGConfig(
    nx=100, nt=100,           # 100 spatial × 100 temporal grid points
    x_min=0.0, x_max=1.0,     # Spatial domain [0, 1]
    T=1.0,                     # Time horizon
    nu=0.01,                   # Viscosity (diffusion coefficient)
    max_iter=50,               # Fixed-point iteration limit
    tol=1e-5,                  # Convergence tolerance
    alpha=0.5                  # Relaxation parameter
)

# Initial distribution (agents start at x=0.3)
x = np.linspace(0, 1, 100)
m0 = np.exp(-50 * (x - 0.3)**2)
m0 /= np.sum(m0) * (x[1] - x[0])

# Terminal cost (agents want to reach x=0.7)
u_terminal = 0.5 * (x - 0.7)**2

# Solve coupled HJB-Fokker-Planck system
u, m, iterations = solve_mfg_1d_rust(
    m0, u_terminal, config,
    lambda_congestion=0.5
)

print(f"✓ Converged in {iterations} iterations")
print(f"Solution: u{u.shape}, m{m.shape}")
# Typical time: 0.4s for 10,000 space-time points

Hidden Markov Model

from optimizr import HMM
import numpy as np

# Fit HMM with regime switching
returns = np.random.randn(1000)
hmm = HMM(n_states=3)
hmm.fit(returns, n_iterations=100)

# Decode most likely state sequence
states = hmm.predict(returns)
print(f"Transition Matrix:\n{hmm.transition_matrix_}")
print(f"Detected states: {states}")

MCMC Sampling

from optimizr import mcmc_sample

# Define log-posterior
def log_likelihood(params, data):
    mu, sigma = params
    return -0.5 * np.sum(((data - mu) / sigma) ** 2) - len(data) * np.log(sigma)

# Sample from posterior
data = np.random.randn(100) + 2.0
samples = mcmc_sample(
    log_likelihood_fn=log_likelihood,
    data=data,
    initial_params=[0.0, 1.0],
    param_bounds=[(-10, 10), (0.1, 10)],
    n_samples=10000,
    burn_in=1000,
    proposal_std=0.1
)

print(f"Posterior mean: {np.mean(samples, axis=0)}")

Optimal Control (New in v0.2.0)

from optimizr import optimal_control
import numpy as np

# Hamilton-Jacobi-Bellman equation solver
# For stochastic control problem: dX_t = μ dt + σ dW_t

# Define problem parameters
grid = np.linspace(-5, 5, 100)
dt = 0.01
horizon = 1.0

# Solve HJB equation
value_function = optimal_control.solve_hjb(
    grid=grid,
    drift=lambda x: -0.1 * x,      # Mean reversion
    diffusion=lambda x: 0.2,        # Constant volatility
    cost=lambda x, u: x**2 + u**2,  # Quadratic cost
    dt=dt,
    horizon=horizon
)

# Compute optimal control policy
policy = optimal_control.compute_policy(value_function, grid)
print(f"Value at origin: {value_function[len(grid)//2]:.4f}")

Information Theory

from optimizr import mutual_information, shannon_entropy

# Calculate mutual information between two variables
x = np.random.randn(1000)
y = 2 * x + np.random.randn(1000) * 0.5

mi = mutual_information(x, y, n_bins=10)
print(f"Mutual Information: {mi:.4f}")

# Calculate entropy
entropy = shannon_entropy(x, n_bins=10)
print(f"Shannon Entropy: {entropy:.4f}")

Algorithm Details

Hidden Markov Models

Implementation of the Baum-Welch algorithm (Expectation-Maximization) for learning HMM parameters:

• Forward-Backward Algorithm: Efficient computation of state probabilities
• Viterbi Decoding: Find most likely state sequence
• Gaussian Emissions: Continuous observation models
• Normalization: Numerical stability for long sequences

Use Cases:

• Regime detection in time series
• Speech recognition
• Biological sequence analysis
• Financial market state identification

MCMC Sampling

Metropolis-Hastings algorithm for sampling from arbitrary probability distributions:

• Adaptive Proposals: Gaussian random walk
• Burn-in Period: Discard initial samples
• Bounded Parameters: Constraint handling
• Convergence Diagnostics: Track acceptance rates

Use Cases:

• Bayesian parameter estimation
• Posterior inference
• Integration of complex distributions
• Uncertainty quantification

Differential Evolution (Enhanced in v0.2.0)

Advanced global optimization for non-convex, multimodal, high-dimensional problems:

5 Mutation Strategies:

• rand/1/bin: Random base vector (exploration)
• best/1/bin: Best individual base (exploitation)
• current-to-best/1/bin: Balanced exploration/exploitation
• rand/2/bin: Two difference vectors (diversity)
• best/2/bin: Best with two differences (aggressive)

Adaptive jDE Algorithm:

• Self-tuning mutation factor (F) and crossover rate (CR)
• Parameter adaptation per individual
• τ₁, τ₂ control adaptation speed
• Eliminates manual parameter tuning

Convergence Features:

• Early stopping with tolerance detection
• Convergence history tracking
• Best fitness evolution monitoring
• Rich diagnostic information

Performance:

• 74-88× faster than SciPy (Python)
• Efficient for 10-1000 dimensional problems
• Memory-efficient population management
• Parallel-ready architecture

Use Cases:

• Hyperparameter optimization (ML/DL)
• Engineering design problems
• Inverse problems and calibration
• Non-smooth, noisy objectives
• Constrained optimization with penalties

Grid Search

Exhaustive search over parameter space:

• Complete Coverage: Evaluate all grid points
• Parallel Ready: Independent evaluations
• Flexible Bounds: Per-parameter ranges
• Best Score Tracking: Return optimal parameters

Use Cases:

• Small parameter spaces
• Benchmark comparisons
• Hyperparameter tuning
• Global optima verification

Information Theory Metrics

Quantify information content and dependencies:

• Mutual Information: I(X;Y) = H(X) + H(Y) - H(X,Y)
• Shannon Entropy: H(X) = -∑ p(x) log p(x)
• Binning Strategy: Histogram-based estimation
• Normalized Variants: Available through Python API

Use Cases:

• Feature selection
• Dependency detection
• Time series analysis
• Causality testing

Mathematical Toolkit (New in v0.2.0)

Centralized mathematical utilities for all algorithms:

Numerical Differentiation:

• gradient(): ∇f(x) with central differences
• hessian(): H(f)(x) second-order derivatives
• jacobian(): J(f)(x) for vector functions
• Configurable step size (h)

Statistics:

• mean(), variance(), std_dev()
• skewness(), kurtosis() for distribution shape
• correlation(), covariance() for dependencies
• Efficient single-pass algorithms

Linear Algebra:

• norm_l1(), norm_l2(), norm_frobenius()
• normalize() for vector/matrix normalization
• trace(), outer_product()
• ndarray-linalg integration

Integration:

• trapz(): Trapezoidal rule
• simpson(): Simpson's rule

Special Functions:

• sigmoid(), softmax()
• soft_threshold() for proximal methods

Use Cases:

• Algorithm development
• Sensitivity analysis
• Statistical inference
• Custom optimization methods

Optimal Control (New in v0.2.0)

Hamilton-Jacobi-Bellman equation solvers for stochastic control:

Features:

• HJB PDE solver with finite difference schemes
• Regime-switching models (Markov chains)
• Jump diffusion processes (Poisson jumps)
• MRSJD (Markov Regime Switching Jump Diffusion)

Components:

• Value function computation
• Optimal policy extraction
• Boundary conditions handling
• Grid-based discretization

Use Cases:

• Portfolio optimization under uncertainty
• Resource management with regime changes
• Risk-sensitive control
• Dynamic programming problems

Performance Benchmarks

Comparison against pure Python/NumPy/SciPy implementations (v0.2.0):

Algorithm	Problem Size	Optimiz-rs (Rust)	NumPy/SciPy	Speedup
DE - rand/1	50D Rosenbrock	285ms	21.2s	74×
DE - best/1	50D Rosenbrock	270ms	23.8s	88×
DE - adaptive jDE	50D Rosenbrock	310ms	24.5s	79×
HMM Fit	10k samples	45ms	3.2s	71×
MCMC Sample	100k iterations	120ms	8.5s	71×
Sparse PCA	1000×100 matrix	180ms	12.5s	69×
Mutual Information	50k points	12ms	380ms	32×
Gradient (numerical)	100D function	8ms	145ms	18×
Hessian (numerical)	50D function	95ms	4.2s	44×

Benchmarks run on Apple M1 Pro, 10 cores, 32GB RAM

Documentation

API Reference

Full API documentation is available in the docs/ directory:

• HMM API
• MCMC API
• Differential Evolution API
• Grid Search API
• Information Theory API

Examples & Tutorials

Complete Jupyter notebook tutorials in examples/notebooks/ (all validated in v0.3.0):

1. Hidden Markov Models - Regime detection, Baum-Welch, Viterbi ✅
2. MCMC Sampling - Metropolis-Hastings, Bayesian inference ✅
3. Differential Evolution - 5 strategies, adaptive jDE, convergence ✅
4. Optimal Control - HJB, regime switching, jump diffusion (theory) ℹ️
5. Real-World Applications - Complete workflows ✅
6. Performance Benchmarks - Rust vs Python comparisons ✅
7. Mean Field Games - Population dynamics, HJB-FP coupling ✅ NEW in v0.3.0

All notebooks tested and production-ready! See NOTEBOOK_AUDIT_REPORT.md for validation details.

Python script examples:

• HMM Regime Detection
• Parallel DE Benchmark
• Polarway-Optimizr Integration
• Timeseries Integration

Mathematical Background

Detailed mathematical descriptions and references:

• HMM Theory
• MCMC Theory
• Differential Evolution Theory - Updated for v0.2.0
• Information Theory

📚 Documentation & Getting Started

Comprehensive documentation is available on ReadTheDocs:

👉 https://optimiz-r.readthedocs.io/en/latest/

The documentation includes:

• 🚀 Quick Start Guide - Get up and running in minutes
• 📖 Installation - Detailed setup instructions for all platforms
• 🎓 Tutorials - Step-by-step guides for each algorithm
• 📚 API Reference - Complete function and class documentation
• 🔬 Theory & Math - Mathematical foundations and references
• 💡 Examples - Real-world use cases and code samples
• ⚡ Performance - Benchmarks and optimization tips

New to Optimiz-rs? Start with the Quick Start Guide or try the Mean Field Games Tutorial.

Development

Building from Source

# Setup development environment
git clone https://github.com/ThotDjehuty/optimiz-r.git
cd optimiz-r

# Install development dependencies
pip install -e ".[dev]"

# Build Rust extension
maturin develop

# Run tests
pytest tests/ -v

# Run Rust tests
cargo test

# Run benchmarks
cargo bench

Code Quality

# Format code
black python/
cargo fmt

# Lint
ruff check python/
cargo clippy

# Type checking
mypy python/

Contributing

Contributions are welcome! Please see CONTRIBUTING.md for guidelines.

Areas for Contribution

• Advanced DE variants (JADE, SHADE, L-SHADE)
• GPU acceleration via CUDA/ROCm (see Roadmap)
• Additional optimization algorithms (PSO, CMA-ES, NES)
• More probability distributions for HMM
• Additional language bindings (R, Julia, JavaScript)
• Documentation improvements and tutorials
• Benchmark comparisons and case studies

License

MIT License - see LICENSE file for details.

Citation

If you use Optimiz-rs in your research, please cite:

@software{optimizr2024,
  title = {Optimiz-rs: High-Performance Optimization Algorithms in Rust},
  author = {HFThot Research Lab},
  year = {2024},
  version = {1.0.0},
  url = {https://github.com/ThotDjehuty/optimiz-r}
}

Acknowledgments

Built with:

• Rust - Systems programming language
• PyO3 - Rust bindings for Python
• Maturin - Build and publish Rust crates as Python packages
• NumPy - Numerical computing in Python

Inspired by:

• scipy.optimize
• scikit-learn
• hmmlearn
• emcee

Contact

• Issues: GitHub Issues
• Discussions: GitHub Discussions
• Website: HFThot Research Lab
• Email: contact@hfthot-lab.eu

---

Optimiz-rs - Fast optimization for data science and machine learning 🚀