Feat/012 krylov lanczos

CHANGELOG.md

@@ -5,6 +5,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/).
 
+## [Unreleased]
+
+### Added
+
+- **Lanczos iteration.** `krylov::lanczos` tridiagonalizes a real symmetric matrix over a compile-time-sized Krylov
+  subspace, returning the projection as a new `TridiagonalMatrix<T, K>` and filling an orthonormal basis. Uses full
+  reorthogonalization, so the basis stays orthonormal to working precision.
+- `storage::Basis<T, K>`, a const-generic view over caller-provided memory holding the `K` basis vectors a Krylov
+  method builds up — the storage-layer piece Arnoldi and GMRES(m) will reuse.
+- `ConvergenceError::Breakdown`, reported when the Krylov subspace turns out to be invariant before reaching the
+  requested dimension (e.g. a repeated eigenvalue, or a starting vector inside a small invariant subspace).
+
 ## [Released]
 
 ## [0.4.0] - 2026-07-11

docs/SUMMARY.md

@@ -30,6 +30,7 @@
     - [SVD](book/06-decompositions/svd.md)
   - [Krylov Methods](book/08-krylov/README.md)
     - [Power Iteration](book/08-krylov/power-iteration.md)
+    - [Lanczos Iteration](book/08-krylov/lanczos.md)
   - [Eigenvalues & Eigenvectors](book/09-eigenvalues.md)
   - [Numerical Stability](book/12-numerical-stability/README.md)
     - [NaN/Inf Policy](book/12-numerical-stability/nan-inf-policy.md)

docs/book/08-krylov/README.md

@@ -1,8 +1,10 @@
 # Krylov Methods
 
 `rustebra` provides iterative eigenvalue methods in `rustebra::krylov`: power iteration for
-the dominant (largest-magnitude) eigenvalue, and inverse power iteration for the eigenvalue
-nearest an arbitrary shift. Unlike the direct decompositions in
-[Decompositions](../06-decompositions/README.md), these refine an estimate over many
-iterations and can fail to converge within a given budget, in addition to the usual
-dimension and non-finite-value failure modes.
+the dominant (largest-magnitude) eigenvalue, inverse power iteration for the eigenvalue
+nearest an arbitrary shift, and Lanczos iteration, which builds an orthonormal basis of a
+Krylov subspace and the symmetric tridiagonal matrix a symmetric operator projects onto
+within it. Unlike the direct decompositions in
+[Decompositions](../06-decompositions/README.md), these refine an estimate (or a basis) over
+many iterations and can fail to converge — or, for Lanczos, to extend the basis further —
+within a given budget, in addition to the usual dimension and non-finite-value failure modes.

docs/book/08-krylov/lanczos.md

@@ -0,0 +1,55 @@
+# Lanczos Iteration
+
+`lanczos` tridiagonalizes a symmetric `n x n` matrix `a` over a `K`-dimensional Krylov
+subspace: starting from a normalized `v0`, it builds an orthonormal basis `Q` of
+`span{v0, a*v0, ..., a^{K-1}*v0}` and returns the projection `T = Qᵗ * a * Q`, a symmetric
+`K x K` tridiagonal matrix. With `K == n` and no breakdown, `T` is orthogonally similar to
+`a` and therefore has exactly its spectrum; with `K < n`, `T`'s eigenvalues (the Ritz values)
+approximate the extreme eigenvalues of `a`, generally converging well before `K` reaches `n`.
+Unlike `power_iteration` and `inverse_power_iteration`, which each refine a single eigenpair,
+Lanczos builds a whole subspace at once — the basis for later extracting several eigenvalues,
+or for feeding other Krylov solvers (CG, MINRES) that need the same tridiagonal projection.
+
+Each step computes `w = a * q_j`, records the diagonal entry `α_j = q_jᵗ * w`, subtracts off
+the components along `q_j` and the previous basis vector `q_{j-1}` (the three-term
+recurrence), and normalizes what remains into `q_{j+1}`, recording its length as the
+off-diagonal entry `β_j`. In floating point, rounding error erodes the basis's orthogonality
+as the Ritz values converge, so every step also re-orthogonalizes `w` against *every* basis
+vector built so far (full reorthogonalization) rather than relying on the three-term
+recurrence alone.
+
+```rust
+{{#include ../../../examples/krylov/lanczos.rs}}
+```
+
+## Breakdown
+
+When the candidate for the next basis vector has (numerically) zero norm relative to `‖a *
+q_j‖`, the Krylov subspace is invariant: there's no new direction to extend the basis with,
+and the call fails with `ConvergenceError::Breakdown` rather than dividing by a vanishing
+norm. This is not a numerical failure to work around — it's a structural property of the
+pairing of `a` and `v0`. A repeated eigenvalue can contribute at most one basis vector to the
+Krylov subspace no matter how large `K` is (the identity matrix breaks down immediately for
+any `v0`, since `a * v0` never points anywhere new), and a `v0` that happens to lie in a
+proper invariant subspace of `a` breaks down as soon as that subspace is exhausted, even with
+a fully distinct spectrum. The remedy is different from a plain convergence failure: retry
+with a smaller `K`, or a different `v0`.
+
+## Gotchas
+
+- `a` is *assumed* symmetric, never verified. For a non-symmetric input the projection isn't
+  actually tridiagonal, and `T` silently misrepresents it — Lanczos has no equivalent of the
+  dimension or non-finite checks for this assumption, since checking symmetry itself would
+  cost as much as the decomposition it's protecting.
+- `tol` has no auto-computed default, the same as `power_iteration` and
+  `inverse_power_iteration` — see
+  [Krylov Tolerance and Convergence Criteria](../../specs/krylov-tolerance-and-convergence.md).
+  A `tol` of `0` detects only exact breakdown.
+- The basis size `K` is a `const` generic on the caller's `Basis` buffer, not a runtime
+  parameter — see
+  [Krylov Basis-Size Const-Generic Convention](../../specs/krylov-basis-size-const-generics.md).
+  `K > n` is a `DimensionMismatch`: an `n`-dimensional space has no `K` orthonormal directions
+  to find.
+- Both `ConvergenceError::ZeroVector` and `ConvergenceError::NonFinite` on `v0` are checked
+  up front, even when `K == 0` means no basis vector is ever written — a `K == 0` call is not
+  a shortcut around input validation, only around the iteration itself.

docs/specs/krylov-basis-size-const-generics.md

@@ -49,6 +49,6 @@ Krylov and other generic numerical code.
 
 ## Status
 
-Not yet implemented. Lanczos, Arnoldi, and GMRES(m) do not exist in the crate yet; this
-convention is decided ahead of their implementation so the first one written establishes the
-pattern correctly rather than needing a later rename.
+Partially implemented. Lanczos now exists and exposes its basis size as `const K: usize`,
+establishing the pattern this convention describes. Arnoldi and GMRES(m) do not exist in the
+crate yet.

docs/specs/krylov-tolerance-and-convergence.md

@@ -9,8 +9,8 @@ stabilization and eigenvector-residual stabilization — both fall within `tol`.
 
 ## Scope
 
-Applies to `power_iteration` and `inverse_power_iteration`, and to any future Krylov solver
-(CG, Lanczos, Arnoldi, GMRES(m)) that shares the same iterative-refinement shape. Does not
+Applies to `power_iteration`, `inverse_power_iteration`, and `lanczos`, and to any future
+Krylov solver (CG, Arnoldi, GMRES(m)) that shares the same iterative-refinement shape. Does not
 apply to the `algorithm::matrix` tolerance-taking functions (rank, SVD, condition-number
 estimation, Cholesky decomposition), which are covered by [[approximate-zero-tolerance]] and
 get an auto-computed default under [[auto-tolerance-defaults]] instead.
@@ -63,4 +63,7 @@ how much precision loss to expect.
 
 Implemented. `power_iteration` requires `tol`; `inverse_power_iteration` requires both `tol`
 and `singular_tol`. Both check eigenvalue and residual stabilization before declaring
-convergence, and neither exposes an auto-computed default.
+convergence, and neither exposes an auto-computed default. `lanczos` also requires `tol`,
+with no default, but as a basis-breakdown threshold rather than an eigenvalue/residual
+convergence check: it has no eigenvalue estimate to stabilize against, only a candidate basis
+vector's norm to compare against the local matrix-vector scale.

examples/krylov/lanczos.rs

@@ -0,0 +1,25 @@
+use rustebra::krylov::lanczos;
+use rustebra::storage::{Basis, StaticStorage};
+
+pub(crate) fn run() {
+    println!("\n== Lanczos iteration ==");
+    // [[4, 1, 2], [1, 3, 1], [2, 1, 5]], symmetric but not already tridiagonal.
+    let a = StaticStorage::new([4.0, 1.0, 2.0, 1.0, 3.0, 1.0, 2.0, 1.0, 5.0]);
+    let v0 = StaticStorage::new([1.0, 1.0, 1.0]);
+    let mut buffer = [0.0; 9];
+    let mut basis = Basis::<f64, 3>::new(&mut buffer, 3).unwrap();
+    let mut scratch = [0.0; 3];
+
+    let t = lanczos(&a, 3, &v0, 1e-12, &mut basis, &mut scratch).unwrap();
+    println!("diagonal = {:?}", t.diagonal());
+    println!("off_diagonal = {:?}", t.off_diagonal());
+
+    // Requesting fewer basis vectors than the matrix dimension (K < n) still produces the
+    // leading block of the same tridiagonal form, at a fraction of the memory: only `K`
+    // vectors of the basis are ever stored.
+    let mut partial_buffer = [0.0; 6];
+    let mut partial_basis = Basis::<f64, 2>::new(&mut partial_buffer, 3).unwrap();
+    let mut partial_scratch = [0.0; 3];
+    let partial_t = lanczos(&a, 3, &v0, 1e-12, &mut partial_basis, &mut partial_scratch).unwrap();
+    println!("partial diagonal (K = 2) = {:?}", partial_t.diagonal());
+}

examples/krylov/main.rs

@@ -0,0 +1,11 @@
+//! Tour `rustebra::krylov`'s subspace methods: Lanczos iteration, which builds an orthonormal
+//! basis of a Krylov subspace and the symmetric tridiagonal matrix a symmetric operator
+//! projects onto within it.
+//!
+//! Run with: `cargo run --example krylov`
+
+mod lanczos;
+
+fn main() {
+    lanczos::run();
+}