Expand description
§Linear Algebra Engine
A custom, lightweight linear algebra module optimized for the specific structures arising in these control problems.
§1. Compressed Sparse Row (CSR) Matrix (CsrMatrix)
A memory-efficient storage format for sparse matrices.
- Problem: The finite difference discretization of the HJB equation produces a highly sparse matrix (tridiagonal or pentadiagonal structure) where most elements are zero.
- Solution: Stores only non-zero elements and their column indices, reducing memory usage from $O(N^2)$ to $O(N)$.
Source: src/linalg/csr.rs
§2. Successive Over-Relaxation (SOR) Solver (solve_sor)
An iterative method for solving linear systems of equations $Ax = b$.
- Problem: The “Evaluation Step” of Policy Iteration requires solving a large, sparse linear system at every time step. Direct solvers (like Gaussian elimination) are too slow ($O(N^3)$) and memory-intensive.
- Solution: SOR is an extension of the Gauss-Seidel method that converges faster by extrapolating the new value (using a relaxation factor $\omega$). It exploits the sparsity of the CSR matrix to solve the system in $O(k \cdot N)$ time, where $k$ is the number of iterations.
Math derivation: docs/math/sor.md
Source: src/linalg/sor.rs
§3. Symmetric Tridiagonal Eigenvalue Solver (eigen_symmetric_tridiagonal)
A specialized solver for computing eigenvalues and eigenvectors of symmetric tridiagonal matrices.
- Problem: The “Exact” analytical solution requires computing the matrix exponential $\exp(-M(T-t))$. For the Avellaneda model, $M$ is symmetric and tridiagonal.
- Solution: Implements the Implicit QL algorithm (a variant of QR). It diagonalizes the matrix $M = Q \Lambda Q^T$, allowing the matrix exponential to be computed efficiently as $\exp(M) = Q \exp(\Lambda) Q^T$.
Math derivation: docs/math/eigen_symmetric_tridiagonal.md
Source: src/linalg/eigen.rs
§4. Benchmarks
Benchmarks run on a symmetric tridiagonal matrix of size 5000×5000.
§4.1 Eigenvalue Decomposition
Computing all eigenvalues of a symmetric tridiagonal matrix.
| Implementation | Method | Time (ms) | Notes |
|---|---|---|---|
| Rust LAPACK | dstev (Intel MKL) | 261.32 | Fastest |
| Python Scipy | eigh_tridiagonal | 297.45 | Uses LAPACK/MKL |
| Rust Faer | selfadjoint_eigenvalues | 13,630 | Dense solver on sparse data |
| Rust Custom | Implicit QL | 101,548 | Pure Rust implementation |
§4.2 Linear System Solve ($Ax=b$)
Solving a tridiagonal linear system.
| Implementation | Method | Time (ms) | Notes |
|---|---|---|---|
| Rust LAPACK | dgtsv (Intel MKL) | 0.088 | Fastest |
| Python LAPACK | dgtsv | 0.118 | Direct wrapper |
| Python Banded | solve_banded | 0.142 | General banded solver |
| Rust SOR | Custom Iterative | 1.08 | ~8x slower than Direct |
| Python Sparse | spsolve | 2.29 | General sparse solver |
§4.3 Running Benchmarks
-
Generate Data:
cargo run --release --bin generate_matrix -- 5000 bench_matrix.txt -
Run Comparison:
cargo run --release --bin benchmark_runner -- bench_matrix.txt
Python equivalents are in benches/python/.