Home

Linear algebra for embedded systems, microcontrollers, and real-time applications.

A hybrid no_std/alloc library. Stack-first by default. Scales to sparse matrices and Krylov subspace solvers when a heap is available.

API Reference · GitHub

Why rustebra exists

Rust currently lacks a linear algebra library that is simultaneously serious about no_std support and complete enough to cover sparse matrices and iterative solvers. Existing options either assume a heap is always available, or only provide a partial set of operations for constrained environments.

rustebra closes that gap.

Design principles

No allocator required by default — The core works entirely on the stack using const generics to fix sizes at compile time.
Allocation is opt-in — Heap-backed structures and algorithms are available behind the alloc feature flag.
Generic over numeric precision — Works across floating-point types, from microcontrollers without double-precision units to desktop systems.
Explicit error handling — Recoverable failures are reported through Result, not panics.

rustebra vs. the competition

Feature	rustebra	ndarray	nalgebra
no_std support	Yes Full	Partial Optional	Partial Optional
Stack-only (no heap required)	Yes Default	No No	Yes For fixed-size
Sparse matrices	Yes v0.3.0+ (COO, CSR, CSC)	No Separate crate	Partial Limited
Krylov solvers	Yes v0.4.0+ (power iteration)	Partial Via ndarray-linalg	No Not in core
3D math/graphics	No Not focused	No Not provided	Yes Excellent
BLAS/LAPACK integration	No No	Yes Excellent	No Pure Rust
Maturity	Early v0.4.0	Yes Mature	Yes Mature
Embedded systems	Yes Best choice	No Poor fit	Partial For fixed-size only

When to use rustebra

Use rustebra if:

You need linear algebra without dynamic allocation (embedded, real-time, microcontroller)
You’re working with sparse matrices in an embedded context
You want no_std + optional alloc (best of both worlds)
You need predictable stack-only memory

Use ndarray if:

You need production BLAS/LAPACK routines (scientific computing, data science)
You’re comfortable with heap allocation and want optimal performance
You need large matrices with sophisticated solvers
Building NumPy-like workflows in Rust

Use nalgebra if:

You need 3D graphics, robotics, or game engine math (Points, Isometries, Rotations)
You want optional no_std support with fixed-size matrices
Building low-level geometric transformations

Getting started

[dependencies]
rustebra = "0.4.0"

# Optional: heap-backed structures and Krylov solvers
rustebra = { version = "0.4.0", features = ["alloc"] }

# no_std build (default)
cargo build
cargo test

# with alloc feature
cargo build --features alloc
cargo test --features alloc

Explore the docs

API Reference — generated from cargo doc
Algorithms — mathematical reference for every algorithm
Design & Specs — design decisions behind the library
Contributing Guide — how to get involved

Licensed under the Apache License 2.0.

Algorithms

Mathematical reference for every algorithm implemented in this project. Each document covers the theory, the method, and when to use it. Language-agnostic.

v0.2.0 — Matrix decompositions and structural properties

Cofactor expansion
Rank
LU decomposition
QR decomposition
Cholesky decomposition
Singular Value Decomposition (SVD)
Condition number

Cholesky Decomposition

What it computes

A factorization of a symmetric positive-definite matrix A into:

A = L · Lᵀ

where L is lower triangular with positive diagonal entries. It is the “square root” of a matrix in the sense that multiplying L by its own transpose recovers A.

Intuition

LU decomposition works for any invertible square matrix, but it doesn’t exploit symmetry. If you know A is symmetric and positive-definite, you can do roughly half the work: instead of finding two different triangular matrices L and U, you find one (L) and observe that U is simply Lᵀ.

The positive-definite condition (all eigenvalues positive) guarantees that the square root exists and that all diagonal entries of L are real and positive — which is what makes the algorithm work without pivoting.

What “symmetric positive-definite” means

Symmetric: Aᵢⱼ = Aⱼᵢ for all i, j (the matrix equals its own transpose).
Positive-definite: for every non-zero vector x, xᵀAx > 0. Intuitively, A doesn’t flip any vector to point in an opposite direction.

Common sources: covariance matrices in statistics, stiffness matrices in structural engineering, kernel matrices in machine learning.

Method

Compute L column by column. For column j:

Lⱼⱼ = √(Aⱼⱼ − Σ Lⱼₖ²)        (k from 1 to j−1)

Lᵢⱼ = (Aᵢⱼ − Σ LᵢₖLⱼₖ) / Lⱼⱼ  (k from 1 to j−1, for i > j)

If at any point the expression inside the square root is negative or zero, the matrix is not positive-definite and the decomposition does not exist.

Example:

A = | 4   2   2 |       L = | 2    0    0  |
    | 2   5   3 |           | 1    2    0  |
    | 2   3   6 |           | 1    1    √3 |

Verify: L · Lᵀ = A.

Computational cost

O(n³/3) — approximately half the cost of LU decomposition for the same matrix size, because the symmetry means only the lower triangle needs to be computed.

When to use it

Solving Ax = b when A is symmetric positive-definite — the most efficient general method in that case.
Sampling from multivariate normal distributions (used in statistics and Monte Carlo methods).
As a fast check for positive-definiteness: if Cholesky succeeds, the matrix is positive- definite; if it fails (negative under the square root), it is not.

Limitations

Only applicable to symmetric positive-definite matrices — will fail or produce incorrect results otherwise.
Requires a square root operation at each diagonal step.
Does not need pivoting (unlike LU), which simplifies the implementation but also means there is no fallback if the positive-definite condition is violated.

Cofactor Expansion

What it computes

The determinant of a square matrix.

Intuition

The determinant measures how much a matrix scales areas (in 2D) or volumes (in 3D) when applied as a linear transformation. A determinant of 2 means areas double; a determinant of 0 means the transformation collapses space into a lower dimension, losing information irreversibly.

Method

Pick any row or column. For each element in that row or column, multiply it by the determinant of the submatrix obtained by removing that element’s row and column (called the minor), then alternate signs following the checkerboard pattern:

+ - + - ...
- + - + ...
+ - + - ...

The sum of those signed products is the determinant.

Base cases:

For a 1×1 matrix:

det([a]) = a

For a 2×2 matrix (the formula most people remember):

det | a  b | = a·d − b·c
    | c  d |

General case — expanding along the first row:

det(A) = Σ (−1)^(1+j) · a₁ⱼ · det(M₁ⱼ)
          j

where M₁ⱼ is the minor obtained by deleting row 1 and column j.

3×3 example:

det | a  b  c |
    | d  e  f | = a·(e·i − f·h) − b·(d·i − f·g) + c·(d·h − e·g)
    | g  h  i |

Computational cost

O(n!) — grows factorially with matrix size. This is impractical for large matrices. For matrices larger than approximately 4×4, LU decomposition computes the same result far more efficiently.

When to use it

Small matrices (up to approximately 4×4) where the direct formula is exact and cheap.
Educational contexts where the explicit formula is needed.
As a building block for computing the inverse via the adjugate matrix.

Limitations

Impractical for large matrices due to O(n!) cost.
Accumulates floating-point rounding errors as matrix size grows.
Only defined for square matrices.

Condition Number

What it computes

A single non-negative number κ(A) that measures how sensitive the solution of a linear system Ax = b is to small changes in A or b.

κ(A) = σ₁ / σₙ

where σ₁ is the largest singular value of A and σₙ is the smallest.

Equivalently, for invertible matrices:

κ(A) = ‖A‖ · ‖A⁻¹‖

Intuition

Suppose you solve Ax = b and get a solution x. Now perturb b by a tiny amount δb. How much does the solution change?

The condition number bounds the answer:

‖δx‖ / ‖x‖  ≤  κ(A) · ‖δb‖ / ‖b‖

A condition number of 1 means the system is perfectly well-conditioned: a 1% error in b causes at most a 1% error in x. A condition number of 10⁶ means a tiny relative error in b can cause an error up to a million times larger in x — the system is nearly impossible to solve accurately with floating-point arithmetic.

Geometrically: a large condition number means the matrix A is “almost singular” — it maps vectors in some direction to nearly zero, making it nearly impossible to tell which input produced a given output.

Scale and interpretation

κ(A)	Interpretation
1	Perfect — no amplification of errors
10 – 100	Well-conditioned — safe for most computations
10³ – 10⁶	Moderately ill-conditioned — results may lose
	several digits of precision
> 10⁸	Severely ill-conditioned — floating-point results
	may be essentially meaningless
∞	Singular matrix — no solution or infinitely many

As a rule of thumb: if κ(A) ≈ 10^k, you lose approximately k digits of precision in the solution.

Method

The most reliable method is via SVD:

κ(A) = σ_max / σ_min

For symmetric positive-definite matrices, the eigenvalues equal the squared singular values, so:

κ(A) = λ_max / λ_min

A cheaper but less reliable estimate uses the LU decomposition — several condition number estimators (LAPACK-style) exist that avoid a full SVD while giving a good approximation.

When to use it

Before solving a linear system, to anticipate how accurate the solution can be.
When comparing different formulations of the same problem — a better-conditioned formulation gives more accurate results with the same arithmetic.
After computing a decomposition, as a sanity check on the result.
In iterative methods (Krylov), a large condition number is why preconditioning is needed: it transforms the problem into a better-conditioned one.

Relationship to other algorithms

Computed exactly via SVD (most reliable).
Estimated cheaply after LU decomposition (less reliable but faster).
Directly motivates preconditioning in Krylov methods.
A condition number of ∞ (or very large) is the formal definition of a matrix being singular (or nearly singular).

LU Decomposition

What it computes

A factorization of a square matrix A into two triangular matrices:

A = L · U

where L is lower triangular (zeros above the diagonal) and U is upper triangular (zeros below the diagonal).

Intuition

Solving a linear system Ax = b directly is expensive. But if A is already triangular, solving is trivial — you just substitute forward or backward. LU decomposition transforms the general problem into two easy triangular ones.

Once you have L and U, solving Ax = b becomes:

Solve Ly = b for y (forward substitution — trivial because L is lower triangular).
Solve Ux = y for x (back substitution — trivial because U is upper triangular).

The decomposition pays for itself when you need to solve multiple systems with the same A but different right-hand sides b — compute L and U once, reuse them many times.

Method

Gaussian elimination, recorded as a matrix factorization.

The same elimination steps used for row reduction can be captured as multiplications by elementary lower-triangular matrices. Collecting those elimination steps gives L; the result of the elimination is U.

At each step k, for each row i below row k:

multiplier mᵢₖ = aᵢₖ / aₖₖ
row i ← row i − mᵢₖ · row k

The multipliers mᵢₖ fill the lower triangle of L; the remaining matrix is U.

Example:

A = | 2  1  1 |       L = | 1    0   0 |    U = | 2  1   1  |
    | 4  3  3 |           | 2    1   0 |        | 0  1   1  |
    | 8  7  9 |           | 4    3   1 |        | 0  0   2  |

Verify: L · U = A.

Pivoting

If a diagonal entry (pivot) is zero or very small, the division aᵢₖ / aₖₖ becomes undefined or numerically unstable. Partial pivoting solves this by swapping rows before each step to put the largest available entry on the diagonal:

P · A = L · U

where P is a permutation matrix recording the row swaps. Partial pivoting is the standard in practice — LU without pivoting can fail even when the matrix is technically invertible.

Computing the determinant via LU

Once U is available, the determinant is the product of its diagonal entries, adjusted for the sign of the row permutations:

det(A) = (−1)^s · u₁₁ · u₂₂ · ... · uₙₙ

where s is the number of row swaps performed during pivoting. This is O(n²) — far cheaper than cofactor expansion.

Computational cost

O(n³) for the decomposition. O(n²) for each subsequent solve once L and U are known.

When to use it

Solving general square linear systems Ax = b, especially when solving multiple times with the same A.
Computing the determinant efficiently for matrices larger than approximately 4×4.
Computing the matrix inverse (though direct solve is usually preferable).
As a building block for other algorithms.

Limitations

Only directly applicable to square matrices.
Requires partial pivoting for numerical stability.
Not the best choice when A has known special structure (symmetric positive-definite → Cholesky; orthogonal → QR).

QR Decomposition

What it computes

A factorization of any m×n matrix A (with m ≥ n) into:

A = Q · R

where Q is orthogonal (its columns are unit vectors, all perpendicular to each other, and Qᵀ·Q = I) and R is upper triangular.

Intuition

An orthogonal matrix Q is a pure rotation or reflection — it doesn’t distort lengths or angles. R is triangular, so easy to work with. Together they let you decompose any linear transformation into “rotate/reflect first, then scale and shear”.

Because Qᵀ = Q⁻¹ (a property unique to orthogonal matrices), solving Ax = b via QR becomes:

QRx = b  →  Rx = Qᵀb

which is a triangular system — solved immediately by back substitution.

Methods

Gram-Schmidt orthogonalization

Take the columns of A one by one. For each new column, subtract its projections onto all previously processed columns to remove any component in their direction, then normalize to unit length. The normalized columns form Q; the projection coefficients fill R.

For columns a₁, a₂, …, aₙ:

u₁ = a₁
u₂ = a₂ − proj_{u₁}(a₂)
u₃ = a₃ − proj_{u₁}(a₃) − proj_{u₂}(a₃)
...

qᵢ = uᵢ / ‖uᵢ‖

where proj_u(a) = (a·u / u·u) · u

Simple to understand, but accumulates floating-point errors for nearly dependent columns. Modified Gram-Schmidt (projecting against already-orthogonalized columns rather than the originals) reduces this error significantly.

Householder reflections

Rather than building Q column by column, apply a sequence of reflection matrices Hₖ that zero out the entries below the diagonal one column at a time:

Hₙ · ... · H₂ · H₁ · A = R
Q = H₁ · H₂ · ... · Hₙ

Each Hₖ is a Householder reflector of the form:

H = I − 2·vvᵀ / (vᵀv)

chosen so that it maps a given vector onto a multiple of a coordinate axis. More numerically stable than Gram-Schmidt and preferred in practice.

Computational cost

O(mn²) for an m×n matrix — more expensive than LU for square systems, but applicable to non-square matrices and more numerically stable.

When to use it

Least-squares problems: solving overdetermined systems (more equations than unknowns), where no exact solution exists and the best approximate solution is sought.
Eigenvalue algorithms: QR iteration is the basis of the standard method for computing all eigenvalues of a matrix.
Arnoldi and Lanczos iterations: QR is used internally to orthogonalize the Krylov basis.
When numerical stability matters more than raw speed.

Limitations

More expensive than LU for square systems when stability is not a concern.
The full Q matrix is large; often only Qᵀb is needed, not Q itself.

Rank

What it computes

The rank of a matrix is the number of linearly independent rows (or equivalently, linearly independent columns). It tells you the effective dimensionality of the information the matrix carries.

Intuition

If a matrix has 4 rows but one of them is a linear combination of the others (e.g. it equals the sum of the first two rows), that row adds no new information. The rank counts only the rows that are genuinely independent — the rest are redundant.

A matrix of size m×n has rank at most min(m, n). When the rank equals min(m, n), the matrix is said to have full rank. When it is lower, the matrix is rank-deficient.

Method

Gaussian elimination to row echelon form.

Apply a sequence of elementary row operations to reduce the matrix:

Find the leftmost column with a non-zero entry (the pivot column).
Swap rows if needed to bring a non-zero entry to the top of that column.
Subtract multiples of that row from all rows below it to zero out every entry below the pivot.
Repeat for the submatrix below and to the right of the current pivot.

The result is a matrix in row echelon form: each row either starts with a leading non-zero entry (a pivot) further to the right than the row above it, or is entirely zero.

The rank is the number of non-zero rows in the result.

Example:

| 1  2  3 |       | 1  2  3 |
| 2  4  6 |  →    | 0  0  0 |   rank = 1
| 3  6  9 |       | 0  0  0 |

Every row was a multiple of the first; after elimination, only one non-zero row remains.

| 1  0  2 |       | 1  0  2 |
| 0  1  3 |  →    | 0  1  3 |   rank = 2
| 0  0  0 |       | 0  0  0 |

Two independent rows survive.

Computational cost

O(m·n·min(m,n)) — efficient for matrices of practical size.

When to use it

Checking whether a linear system has a unique solution, infinitely many, or none.
Detecting redundancy in a set of equations or vectors.
As a prerequisite for understanding SVD and the condition number.

Relationship to other algorithms

A square matrix has full rank if and only if its determinant is non-zero.
The rank is the number of non-zero singular values (see SVD).
A rank-deficient matrix has no inverse.

Singular Value Decomposition (SVD)

What it computes

A factorization of any m×n matrix A into three matrices:

A = U · Σ · Vᵀ

where:

U is an m×m orthogonal matrix (the left singular vectors).
Σ is an m×n diagonal matrix with non-negative entries σ₁ ≥ σ₂ ≥ … ≥ 0 on the diagonal (the singular values).
V is an n×n orthogonal matrix (the right singular vectors).

Intuition

Every linear transformation — no matter how complex — can be broken into three simple steps:

Rotate/reflect the input space (Vᵀ).
Scale each axis independently (Σ — the singular values are the scale factors).
Rotate/reflect the output space (U).

The singular values tell you how much A stretches or compresses space in each “direction”. Large singular values correspond to directions where A amplifies; small ones where it barely affects the input; zero singular values correspond to directions that get completely collapsed (those are the nullspace of A).

Visually: SVD reveals the “principal axes” of the transformation.

Singular values and their meaning

The number of non-zero singular values equals the rank of A.
The largest singular value σ₁ is the maximum amount A can stretch any unit vector.
The smallest non-zero singular value σᵣ is the minimum stretch in the column space.
If any singular value is zero, the matrix is rank-deficient (same as having determinant 0 for square matrices).

Method

Computing SVD directly is a multi-step process:

Form AᵀA (an n×n symmetric positive-semidefinite matrix).
Compute its eigendecomposition: AᵀA = V · D · Vᵀ, where D is diagonal with the eigenvalues on the diagonal.
Singular values: σᵢ = √λᵢ where λᵢ are the eigenvalues of AᵀA.
Left singular vectors: uᵢ = A·vᵢ / σᵢ for each non-zero singular value.

In practice, the eigendecomposition in step 2 is computed via the QR algorithm (iterative QR decompositions applied repeatedly until the matrix converges to diagonal form), not by solving the characteristic polynomial directly.

Truncated SVD

Often only the k largest singular values are needed (k ≪ min(m,n)). Computing only those (truncated SVD) is far cheaper than the full decomposition and sufficient for most applications.

Computational cost

O(min(m,n) · m · n) — the most expensive decomposition in this library. For large matrices, iterative methods (Lanczos, Arnoldi) are used to approximate only the needed singular values.

When to use it

Least-squares problems: the most robust method, handling rank-deficient cases where QR would struggle.
Pseudoinverse: A⁺ = V · Σ⁺ · Uᵀ, where Σ⁺ inverts the non-zero diagonal entries.
Low-rank approximation: keep only the k largest singular values and vectors to get the best rank-k approximation of A (used in data compression, PCA, and dimensionality reduction).
Condition number: σ₁ / σₙ (see condition number document).
Rank determination: count non-zero singular values (more numerically reliable than Gaussian elimination for ill-conditioned matrices).

Limitations

The most computationally expensive decomposition in this library.
Forming AᵀA explicitly can amplify floating-point errors — bidiagonalization-based algorithms (Golub-Reinsch) avoid this but are more complex to implement.

Book

This book documents rustebra, a hybrid no_std/alloc linear algebra crate for Rust. It starts with the foundations — scalars and vectors — and builds up through matrices, sparse storage, decompositions, and Krylov subspace solvers.

Part I covers the foundations: the Scalar trait that algorithms in this crate are generic over, and the StaticVector/DynamicVector types built on top of it. Later parts assume familiarity with these building blocks.

Getting Started

rustebra is a hybrid no_std/alloc linear algebra crate for Rust, scaling from embedded targets to dynamic Krylov subspace solvers. The crate is no_std by default (#![cfg_attr(not(test), no_std)] in lib.rs), so it has no dependency on the standard library and no allocator requirement unless you opt in.

Adding the dependency

Add rustebra to Cargo.toml:

[dependencies]
rustebra = "0.3"

`no_std` vs. `alloc`

By default, only stack-allocated, compile-time-sized types are available — for example StaticVector<T, N>, whose length N is a const generic. This works on embedded targets with no heap.

Enabling the alloc feature pulls in extern crate alloc and unlocks heap-allocated, runtime-sized types such as DynamicVector<T>:

[dependencies]
rustebra = { version = "0.3", features = ["alloc"] }

A first example

The core building block is Scalar, a trait implemented for f32 and f64 that defines the arithmetic surface (add, sub, mul, div, sqrt, sin, cos) the rest of the crate is generic over. Vector types like StaticVector are built on top of it:

#![allow(unused)]
fn main() {
use rustebra::vector::StaticVector;

let a = StaticVector::new([1.0, 2.0, 3.0]);
let b = StaticVector::new([4.0, 5.0, 6.0]);

println!("a + b = {:?}", a.add(&b));
println!("||a|| = {:.4}", a.norm());
}

StaticVector::new infers N from the array literal’s length, so a and b here are both StaticVector<f64, 3>. See Vectors for the full construction and operations surface.

Scalars & Numeric Types

Every algorithm in rustebra is generic over a Scalar type rather than hard-coded to f32 or f64. Scalar (defined in rustebra::scalar) is the minimal arithmetic surface the rest of the crate depends on: the additive and multiplicative identities (zero, one), the four basic operations (add, sub, mul, div), and the transcendental functions needed by the algorithm layer (sqrt, sin, cos). f32 and f64 are the two implementors provided today.

#![allow(unused)]
fn main() {
use rustebra::scalar::Scalar;

fn double<T: Scalar>(x: T) -> T {
    x.add(x)
}
}

Note the call style for the transcendental methods: Scalar::sqrt(4.0f64) rather than 4.0f64.sqrt(). f64/f32 already have an inherent sqrt from std, which would shadow the trait method if called with dot syntax, so code that needs to stay generic over T: Scalar (or that wants the trait’s specific implementation) calls it as an associated function instead.

How `sqrt`, `sin`, and `cos` are computed

Scalar is implemented for f32/f64 without relying on std’s floating-point intrinsics for these three methods, since the crate is no_std by default:

sqrt uses a fixed-iteration Newton-Raphson (Babylonian) iteration: x_{n+1} = (x_n + self / x_n) / 2. self == 0 short-circuits to 0 rather than iterating (the formula would otherwise divide by zero on the first step), and negative inputs return 0 rather than NaN or a panic, since Scalar is an infallible trait that must also support future non-float implementors with no NaN representation.
sin and cos use a fixed-iteration Taylor series expansion around zero, after first reducing the input to [-pi, pi] by subtracting the nearest multiple of 2*pi.

In all three cases the iteration count is fixed rather than convergence-checked, so the amount of work done is predictable and independent of the input — important in no_std contexts where there’s no easy way to bound worst-case iteration count otherwise.

`FloatTolerance`

Comparing floating-point results for approximate equality requires a tolerance, and rustebra doesn’t hard-code one into Scalar itself. Instead, FloatTolerance is a separate trait (also implemented for f32/f64) that adds a single method, epsilon(), reporting the type’s machine epsilon — the smallest positive value e such that T::one().add(e) != T::one(). It’s kept separate from Scalar so that a hypothetical future Scalar implementor with no meaningful notion of machine epsilon (e.g. a fixed-point or exact-rational type) isn’t forced to implement it just to satisfy Scalar.

#![allow(unused)]
fn main() {
use rustebra::scalar::FloatTolerance;

assert_eq!(f64::epsilon(), f64::EPSILON);
}

Gotchas

Call the transcendental methods as Scalar::sqrt(x) / Scalar::sin(x) / Scalar::cos(x) when x’s type isn’t pinned to a concrete float — dot-call syntax (x.sqrt()) resolves to std’s inherent method on concrete f32/f64 values, not the trait method, and won’t compile at all in a generic T: Scalar context.
Scalar::sqrt of a negative number returns 0, not NaN — don’t rely on NaN propagation to detect a negative input; check the sign yourself if that matters to your code.

Vectors

rustebra provides two vector types: StaticVector<T, N>, stack-allocated with a compile-time length, and DynamicVector<T>, heap-allocated behind the alloc feature with a runtime length. This section covers how to construct them, the arithmetic and dot-product operations they support, and how to compute their norm.

Construction

rustebra vectors come in two flavors, distinguished by where their elements live and how their length is tracked.

StaticVector<T, N> is stack-allocated: N is a const generic, so the length is part of the type itself and is known at compile time. Two StaticVectors with different N are different types — the compiler rejects mismatched-length operations before your code ever runs. This is the default choice in a no_std context, since it requires no allocator.

DynamicVector<T> is heap-allocated (it wraps a Vec<T> internally) and its length is only known at runtime. It’s available behind the alloc feature. Because two DynamicVectors can’t be proven to have the same length at compile time, operations that combine two of them return a Result instead of a bare value.

Both types are constructed directly from their backing collection — StaticVector::new takes a fixed-size array, DynamicVector::new takes a Vec:

#![allow(unused)]
fn main() {
use rustebra::vector::StaticVector;

let v = StaticVector::new([1.0, 2.0, 3.0]);
}

Here N is inferred from the array literal’s length (3), so it rarely needs to be written out explicitly.

Gotchas

StaticVector<T, 3> and StaticVector<T, 4> are unrelated types. If you find yourself needing to add vectors whose lengths aren’t known until runtime, reach for DynamicVector instead — trying to force different Ns together is a compile error, not something you can work around with a cast.
DynamicVector needs the alloc feature enabled. Without it, rustebra::vector only exports StaticVector, and code referencing DynamicVector won’t compile.

Operations

StaticVector and DynamicVector support the same core set of operations: element-wise addition and subtraction, scalar scaling, the dot product, and the Euclidean norm. The example below constructs two StaticVectors and runs each of them in turn.

#![allow(unused)]
fn main() {
use rustebra::vector::StaticVector;

pub(crate) fn run() {
    println!("== StaticVector ==");
    let a = StaticVector::new([1.0, 2.0, 3.0]);
    let b = StaticVector::new([4.0, 5.0, 6.0]);

    println!("a = {a:?}");
    println!("b = {b:?}");
    println!("a + b = {:?}", a.add(&b));
    println!("a - b = {:?}", a.sub(&b));
    println!("a scaled by 2 = {:?}", a.scale(2.0));
    println!("a . b = {}", a.dot(&b));
    println!("||a|| = {:.4}", a.norm());
}
}

Gotchas

add, sub, and dot on StaticVector return a plain value, not a Result — the const generic N guarantees both operands have the same length at compile time, so there’s nothing to fail at runtime. The DynamicVector equivalents return Result<_, LengthMismatch> instead, since two DynamicVectors aren’t guaranteed to match in length.
scale takes the factor by value (T, not &T), matching Scalar’s Copy bound — pass the number directly rather than a reference.

Norms & Distances

Both StaticVector and DynamicVector expose a single norm() method: the Euclidean (L2) norm, ‖v‖ = sqrt(v . v), computed as the square root of the vector’s dot product with itself.

#![allow(unused)]
fn main() {
use rustebra::vector::StaticVector;

let v = StaticVector::new([3.0, 4.0]);
assert_eq!(v.norm(), 5.0);
}

There’s no separate distance method. The Euclidean distance between two vectors is just the norm of their difference, and falls out of the existing sub and norm operations:

#![allow(unused)]
fn main() {
use rustebra::vector::StaticVector;

let a = StaticVector::new([1.0, 2.0, 3.0]);
let b = StaticVector::new([4.0, 6.0, 3.0]);

let distance = a.sub(&b).norm();
assert_eq!(distance, 5.0);
}

For DynamicVector, sub returns a Result (two DynamicVectors aren’t guaranteed to match in length at compile time), so the equivalent distance computation is a.sub(&b)?.norm().

Gotchas

norm() itself never fails or returns a Result — even on DynamicVector — since computing the norm of a single vector has no dimension-mismatch case to report. It’s sub/add/dot between two DynamicVectors that can fail.
sqrt under the hood is Scalar::sqrt, a fixed-iteration approximation rather than a hardware intrinsic (see Scalars & Numeric Types) — norms of extremely large or small magnitude vectors may lose a little precision relative to f64::sqrt.

Matrices

rustebra provides two matrix types: StaticMatrix<T, R, C>, stack-allocated with compile-time-known shape, and DynamicMatrix<T>, heap-allocated behind the alloc feature with a runtime-known shape. This section covers how to construct them and the operations — arithmetic, matrix-vector and matrix-matrix products, transpose, rank, and the LU/QR/SVD decompositions — they support.

Construction

Like vectors, rustebra matrices come in a stack-allocated and a heap-allocated flavor.

StaticMatrix<T, R, C> stores its shape as const generics, so R and C are part of the type and known at compile time. It’s constructed from an array of rows — [[T; C]; R] — and requires no allocator:

#![allow(unused)]
fn main() {
use rustebra::matrix::StaticMatrix;

let m = StaticMatrix::new([[1.0, 2.0], [3.0, 4.0]]);
}

DynamicMatrix<T> is heap-allocated (available behind the alloc feature) and stores its shape — rows and cols — as runtime fields rather than type parameters. It’s constructed from a flat, row-major Vec<T> plus the explicit row and column counts, and returns a Result since the data’s length isn’t statically guaranteed to match rows * cols:

#![allow(unused)]
fn main() {
use rustebra::matrix::DynamicMatrix;

let m = DynamicMatrix::new(2, 2, vec![1.0, 2.0, 3.0, 4.0]).unwrap();
assert_eq!(m.rows(), 2);
assert_eq!(m.cols(), 2);
}

Gotchas

StaticMatrix<T, 2, 3> and StaticMatrix<T, 3, 2> are unrelated types, just like differently-sized StaticVectors — the compiler rejects shape-mismatched operations at compile time rather than at runtime.
DynamicMatrix::new returns Err(DimensionMismatch) rather than panicking if data.len() != rows * cols — check the result rather than assuming construction always succeeds.

Operations

StaticMatrix and DynamicMatrix support the same core set of operations: element-wise addition and subtraction, scalar scaling, matrix-vector and matrix-matrix products, transpose, rank, and the LU, QR, Cholesky, and SVD decompositions, plus condition number. The example below runs each of them on a StaticMatrix.

#![allow(unused)]
fn main() {
use rustebra::matrix::StaticMatrix;
use rustebra::vector::StaticVector;

pub(crate) fn run() {
    println!("== StaticMatrix ==");
    let a = StaticMatrix::new([[1.0, 2.0], [3.0, 4.0]]);
    let b = StaticMatrix::new([[5.0, 6.0], [7.0, 8.0]]);
    let v = StaticVector::new([1.0, 1.0]);

    println!("a = {a:?}");
    println!("b = {b:?}");
    println!("a + b = {:?}", a.add(&b));
    println!("a - b = {:?}", a.sub(&b));
    println!("a scaled by 2 = {:?}", a.mul_scalar(2.0));
    println!("a * v = {:?}", a.mul_vector(&v));
    println!("a * b = {:?}", a.mul_matrix(&b));
    println!("a^T = {:?}", a.transpose());
    println!("det(a) = {:?}", a.determinant());
    println!("rank(a) = {:?}", a.rank());

    let (l, u, swap_count) = a.lu();
    println!("lu(a) = (l = {l:?}, u = {u:?}, swap_count = {swap_count})");

    let (q, r) = a.qr().expect("a has at least as many rows as columns");
    println!("qr(a) = (q = {q:?}, r = {r:?})");

    let spd = StaticMatrix::new([[4.0, 2.0], [2.0, 2.0]]);
    println!("cholesky(spd) = {:?}", spd.cholesky());

    let mut svd_scratch = [0.0; 5 * 2 * 2 + 2 + 2];
    let (svd_u, sigma, v) = a.svd(&mut svd_scratch).expect("scratch is correctly sized");
    println!("svd(a) = (u = {svd_u:?}, sigma = {sigma:?}, v = {v:?})");

    let mut condition_scratch = [0.0; 7 * 2 * 2 + 3 * 2];
    println!(
        "condition_number(a) = {:?}",
        a.condition_number(&mut condition_scratch)
    );
}
}

Gotchas

On StaticMatrix, operations between two matrices (add, sub, mul_matrix, mul_vector) return a plain value, not a Result — the const generics guarantee compatible shapes at compile time. The DynamicMatrix equivalents return Result<_, DimensionMismatch> instead, since two DynamicMatrixs aren’t guaranteed to match in shape.
qr returns Result<_, DimensionMismatch> on both types, because R >= C isn’t something the type system can enforce even for StaticMatrix — stable Rust has no way to bound one const generic against another.
svd and condition_number on StaticMatrix take a caller-provided scratch buffer sized by an explicit formula (5 * C * C + C + R for svd, 7 * N * N + 3 * N for condition_number) — get the size wrong and you get Err(DimensionMismatch) back rather than a panic. DynamicMatrix’s svd/condition_number allocate their own scratch space internally and don’t take this parameter.
determinant on StaticMatrix<T, N, N> returns Err(DeterminantError::MatrixTooLargeWithoutAlloc) if the alloc feature is disabled and N > 4 — for larger matrices without alloc, use rustebra::algorithm::matrix::determinant_lu directly with your own scratch buffer.
cholesky requires self to be symmetric positive-definite; it returns Err(CholeskyError::NotPositiveDefinite) rather than a garbage result or a panic when that doesn’t hold.

Sparse Matrices

rustebra also supports sparse matrices — matrices where most entries are zero and only the non-zeros are stored. This section covers the storage formats (COO, CSR, CSC), how to construct and convert between them, and the operations (scaling, matrix-vector and matrix-matrix products, addition, pruning) available on them. Sparse support requires the alloc feature.

CSR/CSC

Compressed sparse row (CSR) and compressed sparse column (CSC) are the two primary sparse storage formats rustebra operates on. Both store only the non-zero entries, using three parallel arrays instead of a dense grid.

CSR (CsrMatrix<T>) stores:

row_ptr — length rows + 1; the range row_ptr[i]..row_ptr[i + 1] gives the slice of col_indices/values belonging to row i.
col_indices — the column index of each stored entry.
values — the value of each stored entry.

CSC (CscMatrix<T>) is the transpose layout: a col_ptr of length cols + 1 plays the role row_ptr plays for CSR, and row_indices plays the role col_indices plays.

#![allow(unused)]
fn main() {
use rustebra::sparse::CsrMatrix;

// 3x3 identity: row i's single entry is at column i, value 1.0.
let eye = CsrMatrix::new(3, 3, vec![0, 1, 2, 3], vec![0, 1, 2], vec![1.0_f64, 1.0, 1.0])
    .unwrap();
assert_eq!(eye.row_ptr(), &[0, 1, 2, 3]);
}

Within a row (CSR) or column (CSC), the stored indices don’t have to be sorted — but every index must be in-bounds, and the pointer array must be non-decreasing, starting at 0 and ending at the total non-zero count. CsrMatrix::new/CscMatrix::new validate these invariants and return Err rather than constructing a malformed matrix.

`SortedCsrMatrix` / `SortedCscMatrix`

SortedCsrMatrix<T> and SortedCscMatrix<T> wrap a CsrMatrix/CscMatrix and add the additional guarantee that indices within each row/column are in ascending order. This enables O(log(nnz/rows)) binary-search lookup of a specific entry, and is a precondition some algorithms (sparse triangular solves, certain preconditioners) require. Both types implement Deref to their unsorted counterpart, so every read-only accessor is available without unwrapping. Construct one directly with SortedCsrMatrix::from_csr/ SortedCscMatrix::from_csc (which sorts, paying the cost up front), or get one as the output of any operation that produces sorted results as a side effect of its own algorithm — coo_to_csr, csr_to_csc, csc_to_csr, spmm_csr, add_csr, and add_csc.

Gotchas

A CSR matrix with unsorted column indices is still a valid CsrMatrix — sortedness is an opt-in stronger guarantee via SortedCsrMatrix, not a base invariant of CsrMatrix itself.
Storing an explicit zero is legal in CsrMatrix/CscMatrix::new (it doesn’t reject zero-valued entries), but validate_csr-style validation used elsewhere treats an explicit zero as a violation — don’t assume every zero has been pruned just because construction succeeded.

Construction & Conversion

Sparse matrices can be built directly in CSR or CSC format, or assembled from unordered (row, col, value) triplets via CooMatrix and converted afterward. The three examples below construct a matrix in each of the three formats:

#![allow(unused)]
fn main() {
use rustebra::sparse::CooMatrix;

pub(crate) fn run() {
    println!("\n== CooMatrix ==");

    let eye = CooMatrix::new(3, 3, vec![0, 1, 2], vec![0, 1, 2], vec![1.0_f64, 1.0, 1.0])
        .expect("valid COO triplets");

    println!("3×3 identity:");
    println!(
        "  rows={}, cols={}, nnz={}",
        eye.rows(),
        eye.cols(),
        eye.nnz()
    );
    println!("  row_indices = {:?}", eye.row_indices());
    println!("  col_indices = {:?}", eye.col_indices());
    println!("  values      = {:?}", eye.values());

    // Duplicate (row, col) entries are legal — they are summed on conversion.
    let dup = CooMatrix::new(2, 2, vec![0, 0, 1], vec![0, 0, 1], vec![3.0_f64, 4.0, 5.0])
        .expect("valid COO with duplicate at (0,0)");
    println!(
        "with duplicate at (0,0): nnz={} (summed on coo_to_csr)",
        dup.nnz()
    );
}
}

#![allow(unused)]
fn main() {
use rustebra::sparse::CsrMatrix;

pub(crate) fn run() {
    println!("\n== CsrMatrix ==");

    // row_ptr has rows+1 entries; row_ptr[i]..row_ptr[i+1] spans row i.
    let eye = CsrMatrix::new(
        3,
        3,
        vec![0, 1, 2, 3],
        vec![0, 1, 2],
        vec![1.0_f64, 1.0, 1.0],
    )
    .expect("valid CSR arrays");

    println!("3×3 identity:");
    println!(
        "  rows={}, cols={}, nnz={}",
        eye.rows(),
        eye.cols(),
        eye.nnz()
    );
    println!("  row_ptr     = {:?}", eye.row_ptr());
    println!("  col_indices = {:?}", eye.col_indices());
    println!("  values      = {:?}", eye.values());
    println!("  row_range(1) = {:?}", eye.row_range(1));
}
}

#![allow(unused)]
fn main() {
use rustebra::sparse::CscMatrix;

pub(crate) fn run() {
    println!("\n== CscMatrix ==");

    // col_ptr has cols+1 entries; col_ptr[j]..col_ptr[j+1] spans column j.
    let eye = CscMatrix::new(
        3,
        3,
        vec![0, 1, 2, 3],
        vec![0, 1, 2],
        vec![1.0_f64, 1.0, 1.0],
    )
    .expect("valid CSC arrays");

    println!("3×3 identity:");
    println!(
        "  rows={}, cols={}, nnz={}",
        eye.rows(),
        eye.cols(),
        eye.nnz()
    );
    println!("  col_ptr     = {:?}", eye.col_ptr());
    println!("  row_indices = {:?}", eye.row_indices());
    println!("  values      = {:?}", eye.values());
    println!("  col_range(2) = {:?}", eye.col_range(2));
}
}

CooMatrix tolerates duplicate (row, col) triplets — they’re logically summed when the matrix is later used, e.g. by conversion to CSR. Converting between formats is a distinct step from construction; the example below covers all four conversion functions (coo_to_csr, csr_to_coo, csr_to_csc, csc_to_csr):

#![allow(unused)]
fn main() {
use rustebra::sparse::{
    CooMatrix, CscMatrix, CsrMatrix, SortedCsrMatrix, coo_to_csr, csc_to_csr, csr_to_coo,
    csr_to_csc,
};

pub(crate) fn run() {
    println!("\n== convert ==");

    // coo_to_csr: duplicate (row, col) entries are summed; columns sorted per row.
    let coo = CooMatrix::new(2, 2, vec![0, 0, 1], vec![0, 0, 1], vec![3.0_f64, 4.0, 5.0])
        .expect("valid COO");
    let csr: SortedCsrMatrix<f64> = coo_to_csr(coo).expect("dimensions fit within limits");
    println!("coo_to_csr (duplicate (0,0) summed: 3+4=7):");
    println!(
        "  row_ptr={:?}  col_indices={:?}  values={:?}",
        csr.row_ptr(),
        csr.col_indices(),
        csr.values()
    );

    // csr_to_coo: expands row_ptr into explicit row indices.
    let csr_eye = CsrMatrix::new(
        3,
        3,
        vec![0, 1, 2, 3],
        vec![0, 1, 2],
        vec![1.0_f64, 1.0, 1.0],
    )
    .expect("valid CSR");
    let coo_out = csr_to_coo(csr_eye).expect("dimensions fit within limits");
    println!(
        "\ncsr_to_coo (3×3 identity): row_indices={:?}",
        coo_out.row_indices()
    );

    // csr_to_csc: transposes storage layout, sorted by (col, row).
    let csr_rect = CsrMatrix::new(
        2,
        3,
        vec![0, 2, 4],
        vec![0, 2, 1, 2],
        vec![1.0_f64, 2.0, 3.0, 4.0],
    )
    .expect("valid 2×3 CSR");
    let csc_out = csr_to_csc(csr_rect).expect("dimensions fit within limits");
    println!("\ncsr_to_csc (2×3):");
    println!(
        "  col_ptr={:?}  row_indices={:?}  values={:?}",
        csc_out.col_ptr(),
        csc_out.row_indices(),
        csc_out.values()
    );

    // csc_to_csr: transposes back, sorted by (row, col).
    let csc_eye = CscMatrix::new(
        3,
        3,
        vec![0, 1, 2, 3],
        vec![0, 1, 2],
        vec![1.0_f64, 1.0, 1.0],
    )
    .expect("valid CSC");
    let csr_back = csc_to_csr(csc_eye).expect("dimensions fit within limits");
    println!(
        "\ncsc_to_csr (3×3 identity): row_ptr={:?}",
        csr_back.row_ptr()
    );
}
}

Gotchas

coo_to_csr sums duplicate (row, col) entries and sorts column indices within each row as a side effect, which is why it returns a SortedCsrMatrix rather than a plain CsrMatrix — you get the stronger guarantee “for free” from the conversion algorithm.
CsrMatrix::new/CscMatrix::new/CooMatrix::new all return Result rather than panicking on malformed input (mismatched array lengths, out-of-bounds indices, an invalid pointer array) — check the result rather than assuming construction always succeeds.

Sparse Operations

Sparse matrices support scaling, matrix-vector and matrix-matrix products against a dense buffer, sparse-sparse addition, sparse-sparse multiplication, and pruning near-zero entries. Each operation is a free function taking the sparse matrix by value or reference, rather than a method — and each is provided in both a CSR and a CSC variant.

Scaling every stored value by a constant:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CscMatrix, CsrMatrix, scale_csc, scale_csr};

pub(crate) fn run() {
    println!("\n== scale ==");

    let csr =
        CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![4.0_f64, 6.0]).expect("valid CSR");
    let scaled = scale_csr(csr, 0.5);
    println!("scale_csr by 0.5: values = {:?}", scaled.values());

    let csc =
        CscMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![4.0_f64, 6.0]).expect("valid CSC");
    let scaled_csc = scale_csc(csc, 0.5);
    println!("scale_csc by 0.5: values = {:?}", scaled_csc.values());
}
}

Multiplying a sparse matrix by a dense vector:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CscMatrix, CsrMatrix, matvec_csc, matvec_csr};

pub(crate) fn run() {
    println!("\n== matvec ==");

    // [2  0]   [1]   [ 2]
    // [0  5] × [3] = [15]
    let m_csr =
        CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 5.0]).expect("valid CSR");
    let y = matvec_csr(&m_csr, &[1.0, 3.0]).expect("dimensions match");
    println!("matvec_csr: [2 0; 0 5] × [1, 3] = {y:?}");

    let m_csc =
        CscMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 5.0]).expect("valid CSC");
    let y = matvec_csc(&m_csc, &[1.0, 3.0]).expect("dimensions match");
    println!("matvec_csc: [2 0; 0 5] × [1, 3] = {y:?}");
}
}

Multiplying a sparse matrix by a dense matrix:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CscMatrix, CsrMatrix, matmat_csc, matmat_csr};

pub(crate) fn run() {
    println!("\n== matmat ==");

    // [2  0]   [1  0]   [2  0]
    // [0  3] × [0  4] = [0 12]   (x is row-major)
    let x = [1.0_f64, 0.0, 0.0, 4.0];

    let m_csr =
        CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 3.0]).expect("valid CSR");
    let y = matmat_csr(&m_csr, &x, 2).expect("dimensions match");
    println!("matmat_csr: [2 0; 0 3] × [1 0; 0 4] = {y:?}");

    let m_csc =
        CscMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 3.0]).expect("valid CSC");
    let y = matmat_csc(&m_csc, &x, 2).expect("dimensions match");
    println!("matmat_csc: [2 0; 0 3] × [1 0; 0 4] = {y:?}");
}
}

Adding two sparse matrices of the same shape:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CscMatrix, CsrMatrix, add_csc, add_csr};

pub(crate) fn run() {
    println!("\n== add ==");

    // [2  0]   [1  3]   [3  3]
    // [0  5] + [0  4] = [0  9]
    let a = CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 5.0]).expect("valid CSR");
    let b = CsrMatrix::new(2, 2, vec![0, 2, 3], vec![0, 1, 1], vec![1.0_f64, 3.0, 4.0])
        .expect("valid CSR");
    let c = add_csr(&a, &b).expect("shapes match");
    println!(
        "add_csr: col_indices={:?}  values={:?}",
        c.col_indices(),
        c.values()
    );

    let a = CscMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![2.0_f64, 5.0]).expect("valid CSC");
    let b = CscMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![1.0_f64, 4.0]).expect("valid CSC");
    let c = add_csc(&a, &b).expect("shapes match");
    println!(
        "add_csc: row_indices={:?}  values={:?}",
        c.row_indices(),
        c.values()
    );
}
}

Sparse-sparse matrix multiplication (spmm_csr), which multiplies two CsrMatrix operands and returns a sorted result:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CsrMatrix, spmm_csr};

pub(crate) fn run() {
    println!("\n== spmm ==");

    // [1  0]   [1  2]   [1  2]
    // [0  2] × [3  4] = [6  8]
    let a = CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![1.0_f64, 2.0]).expect("valid CSR");
    let b = CsrMatrix::new(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![1.0_f64, 2.0, 3.0, 4.0],
    )
    .expect("valid CSR");
    let c = spmm_csr(&a, &b).expect("inner dimensions match");
    println!("spmm_csr: [1 0; 0 2] × [1 2; 3 4]:");
    println!(
        "  row_ptr={:?}  col_indices={:?}  values={:?}",
        c.row_ptr(),
        c.col_indices(),
        c.values()
    );
}
}

Pruning entries within a tolerance of zero:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CsrMatrix, prune_csr};

pub(crate) fn run() {
    println!("\n== prune ==");

    // Entries within [-tolerance, tolerance] are dropped.
    let m = CsrMatrix::new(
        1,
        3,
        vec![0, 3],
        vec![0, 1, 2],
        vec![1e-15_f64, 2.0, -1e-15],
    )
    .expect("valid CSR");
    let pruned = prune_csr(m, 1e-10).expect("dimensions fit within limits");
    println!("prune_csr (tolerance=1e-10, two near-zeros removed):");
    println!("  nnz={}", pruned.nnz());
    println!(
        "  col_indices={:?}  values={:?}",
        pruned.col_indices(),
        pruned.values()
    );
}
}

Sparse matrices also implement SparseLinearOp, an abstraction for “apply this matrix to a dense vector, writing into a caller-supplied buffer” — the interface Krylov solvers are written against so they work with any sparse format:

#![allow(unused)]
fn main() {
use rustebra::sparse::{CsrMatrix, SparseLinearOp};

let eye = CsrMatrix::new(2, 2, vec![0, 1, 2], vec![0, 1], vec![1.0_f64, 1.0]).unwrap();
let mut y = [0.0; 2];
eye.apply(&[3.0, 5.0], &mut y).unwrap();
assert_eq!(y, [3.0, 5.0]);
}

Gotchas

Every operation here returns Result (e.g. DimensionMismatch for shape-incompatible operands) rather than panicking — sparse matrices carry their shape as runtime fields, so there’s no type-level guarantee two operands are compatible the way there is for StaticMatrix.
apply (from SparseLinearOp) writes into a caller-supplied output buffer instead of allocating a new one, so solver loops built on top of it can reuse the same buffer across iterations without allocating per call.
validate_csr/validate_csc (in rustebra::sparse) check the canonical format invariants — including that no stored value is an explicit zero — separately from construction. CsrMatrix::new/CscMatrix::new accept explicit zeros; only validate_csr/validate_csc flag them.

Decompositions

rustebra factors matrices via LU (Gaussian elimination with partial pivoting), QR (orthogonal-triangular factorization), Cholesky (for symmetric positive-definite matrices), and SVD (singular value decomposition). Each is available both as a low-level function in rustebra::algorithm::matrix operating on Storage, and as an ergonomic method on StaticMatrix/DynamicMatrix. This section covers LU, QR, and SVD in detail.

LU

LU decomposition factors a square matrix a as l * u, where l is unit lower triangular (1s on the diagonal) and u is upper triangular, up to a row permutation: l * u == p * a, where p is the permutation built by applying the reported row swaps, in order, to the identity.

rustebra computes it via Gaussian elimination with partial pivoting: at each step, the row with the largest-magnitude entry in the current column is swapped into the pivot position before elimination proceeds, which avoids dividing by a small or zero pivot. The low-level lu function delegates to lu_partial_pivot, which documents the pivoting strategy in more detail; both are available directly, alongside the ergonomic StaticMatrix/DynamicMatrix::lu() method used in Matrix Operations.

#![allow(unused)]
fn main() {
use rustebra::algorithm::matrix::{lu, lu_partial_pivot};
use rustebra::storage::StaticStorage;

pub(crate) fn run() {
    println!("\n== LU decomposition ==");
    // [[4, 3], [6, 3]].
    let a = StaticStorage::new([4.0, 3.0, 6.0, 3.0]);
    let mut l = [0.0; 4];
    let mut u = [0.0; 4];

    let swaps = lu(&a, 2, 2, &mut l, &mut u).unwrap();
    println!("l = {l:?}, u = {u:?} (row swaps: {swaps})");

    let mut l_explicit = [0.0; 4];
    let mut u_explicit = [0.0; 4];
    lu_partial_pivot(&a, 2, 2, &mut l_explicit, &mut u_explicit).unwrap();
    println!("l (explicit partial pivot) = {l_explicit:?}, u = {u_explicit:?}");
}
}

Gotchas

The permutation p isn’t returned as its own matrix — only the number of row swaps performed (swap_count) is reported. This is enough to recover the permutation’s parity (used by determinant), but not the row ordering itself; there’s no direct way to recover p as a matrix from the public API.
lu/lu_partial_pivot are only defined for square matrices — they return Err(DimensionMismatch) for a non-square input, rather than a rectangular LU variant.

QR

QR decomposition factors a rows x cols matrix a (with rows >= cols) as q * r, where q has orthonormal columns and r is upper triangular. rustebra provides two algorithms:

qr_householder — builds q as a full rows x rows orthogonal matrix via Householder reflections, zeroing out the sub-diagonal of each column in turn.
qr_gram_schmidt — modified Gram-Schmidt orthogonalization, producing a rows x cols q (only as many orthonormal columns as a has, not a full square orthogonal matrix). Projecting each column against the running, already-orthogonalized vector (rather than the original column) is what keeps rounding error from compounding as badly as classical Gram-Schmidt does.

qr is the general-purpose entry point and currently delegates to qr_householder. Both algorithms, plus the entry point, are demonstrated below, alongside the ergonomic StaticMatrix/DynamicMatrix::qr() method used in Matrix Operations.

#![allow(unused)]
fn main() {
use rustebra::algorithm::matrix::{qr, qr_gram_schmidt, qr_householder};
use rustebra::storage::StaticStorage;

pub(crate) fn run() {
    println!("\n== QR decomposition ==");
    // [[3, 5], [4, 0]].
    let a = StaticStorage::new([3.0, 5.0, 4.0, 0.0]);
    let mut scratch = [0.0; 2];

    let mut q = [0.0; 4];
    let mut r = [0.0; 4];
    qr(&a, 2, 2, &mut q, &mut r, &mut scratch).unwrap();
    println!("q = {q:?}, r = {r:?}");

    let mut q_householder = [0.0; 4];
    let mut r_householder = [0.0; 4];
    qr_householder(
        &a,
        2,
        2,
        &mut q_householder,
        &mut r_householder,
        &mut scratch,
    )
    .unwrap();
    println!("q (explicit householder) = {q_householder:?}, r = {r_householder:?}");

    let mut q_gram_schmidt = [0.0; 4];
    let mut r_gram_schmidt = [0.0; 4];
    qr_gram_schmidt(
        &a,
        2,
        2,
        &mut q_gram_schmidt,
        &mut r_gram_schmidt,
        &mut scratch,
    )
    .unwrap();
    println!("q (explicit gram-schmidt) = {q_gram_schmidt:?}, r = {r_gram_schmidt:?}");
}
}

Gotchas

qr_householder’s q is rows x rows; qr_gram_schmidt’s q is only rows x cols. Don’t assume both produce a same-shaped q if you swap between them.
If a column of a lies entirely in the span of the previous columns (or is zero), Gram- Schmidt has no direction left to normalize into that column of q — rather than erroring, it leaves that column as 0, since linear dependence is a property of the input, not a malformed call.
Both algorithms require rows >= cols; this can’t be checked by the type system on StaticMatrix (stable Rust can’t bound one const generic against another), so qr() returns Result on both matrix types even though most other StaticMatrix operations don’t need to.

SVD

Singular value decomposition factors any rows x cols matrix a as u * diag(sigma) * vᵗ, where u has orthonormal columns, sigma is a length-cols vector of non-negative singular values sorted in descending order, and v is a cols x cols orthogonal matrix. Unlike LU or Cholesky, a doesn’t need to be square — the SVD exists for any matrix.

rustebra computes it via svd_qr_iteration: it eigendecomposes aᵗ * a using a fixed-iteration QR algorithm (100 iterations, rather than convergence-checked, so behavior stays predictable in no_std contexts) to get v and sigma, then recovers u from a * v scaled by sigma. svd is the general-purpose entry point, delegating to svd_qr_iteration with an automatically-computed tolerance so callers don’t have to pick one. Both are shown below, alongside the ergonomic StaticMatrix/DynamicMatrix::svd() method used in Matrix Operations.

#![allow(unused)]
fn main() {
use rustebra::algorithm::matrix::{svd, svd_qr_iteration};
use rustebra::storage::StaticStorage;

pub(crate) fn run() {
    println!("\n== Singular Value Decomposition (SVD) ==");
    // [[2, 0], [0, 1]].
    let a = StaticStorage::new([2.0, 0.0, 0.0, 1.0]);
    let mut scratch = [0.0; 5 * 2 * 2 + 2 + 2];

    let mut u = [0.0; 4];
    let mut sigma = [0.0; 2];
    let mut v = [0.0; 4];
    svd(&a, 2, 2, &mut u, &mut sigma, &mut v, &mut scratch).unwrap();
    println!("sigma = {sigma:?}, u = {u:?}, v = {v:?}");

    let mut sigma_explicit = [0.0; 2];
    svd_qr_iteration(
        &a,
        2,
        2,
        &mut u,
        &mut sigma_explicit,
        &mut v,
        &mut scratch,
        1e-9,
    )
    .unwrap();
    println!("sigma (explicit QR iteration) = {sigma_explicit:?}");
}
}

Gotchas

svd/svd_qr_iteration need a caller-provided scratch buffer sized 5 * cols * cols + cols + rows — get the length wrong and you get Err(DimensionMismatch) back, not a panic. DynamicMatrix::svd() allocates this internally and doesn’t expose the parameter.
The QR iteration count inside svd_qr_iteration is fixed (100), not convergence-checked — for matrices with closely clustered singular values, more iterations might be needed for full precision than this budget provides.

Krylov Methods

rustebra provides iterative solvers and eigenvalue methods in rustebra::krylov:

Power iteration for the dominant (largest-magnitude) eigenvalue
Inverse power iteration for the eigenvalue nearest an arbitrary shift
Conjugate Gradient (CG) for solving symmetric positive-definite linear systems
Lanczos iteration, which builds an orthonormal basis of a Krylov subspace and the symmetric tridiagonal matrix a symmetric operator projects onto within it
Arnoldi iteration, the non-symmetric counterpart to Lanczos, which builds an orthonormal basis of a Krylov subspace and the upper Hessenberg matrix a general operator projects onto within it
GMRES(m), a restarted iterative solver for general (possibly non-symmetric) linear systems, built on top of Arnoldi iteration

Unlike the direct decompositions in Decompositions, 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.

Power Iteration

power_iteration computes the dominant (largest-magnitude) eigenvalue of an n x n matrix a, along with a corresponding unit-length eigenvector. Starting from an initial guess v0 (normalized first), each iteration multiplies the current estimate by a and renormalizes: v_{k+1} = a * v_k / ‖a * v_k‖. The eigenvalue estimate is the Rayleigh quotient λ_k = v_kᵗ * (a * v_k). Renormalizing every iteration is what keeps repeated multiplication from overflowing or underflowing — only the iterate’s direction evolves, its magnitude is reset to 1 each round.

#![allow(unused)]
fn main() {
use rustebra::krylov::power_iteration;
use rustebra::storage::StaticStorage;

// [[2, 0], [0, 1]]; dominant eigenvalue 2, eigenvector (1, 0).
let a = StaticStorage::new([2.0, 0.0, 0.0, 1.0]);
let v0 = StaticStorage::new([1.0, 1.0]);
let mut eigenvector = [0.0; 2];
let mut scratch = [0.0; 2];

let eigenvalue = power_iteration(&a, 2, &v0, 100, 1e-10, &mut eigenvector, &mut scratch)
    .unwrap();
assert!((eigenvalue - 2.0).abs() < 1e-9);
}

Iteration stops once both the eigenvalue estimate and the eigenpair residual (‖a * v_k - λ_k * v_k‖) fall within the requested tolerance, relative to the eigenvalue’s magnitude. Convergence speed depends on how well-separated the dominant eigenvalue is from the second-largest: the eigenvector error shrinks by roughly |λ2 / λ1| per iteration.

Inverse power iteration

Plain power iteration can only find the largest-magnitude eigenvalue. inverse_power_iteration targets the eigenvalue nearest an arbitrary shift instead, by applying the same direction-refinement loop to the operator (a - shift * I)⁻¹ — computed via one partial-pivoted LU factorization up front, then a triangular solve per iteration rather than a plain matrix-vector product. With shift == 0, this finds the smallest-magnitude eigenvalue. Convergence rate becomes |λ1 - shift| / |λ2 - shift| (nearest and second-nearest the shift), so a shift close to a good eigenvalue estimate converges very fast.

Gotchas

Power iteration doesn’t converge — and returns ConvergenceError::MaxIterationsExceeded — when the two largest-magnitude eigenvalues are exactly tied in magnitude but not value (e.g. λ2 == -λ1): the iterate oscillates forever rather than settling.
inverse_power_iteration returns ConvergenceError::SingularShift if shift coincides with (or lands within singular_tol of) an actual eigenvalue, rather than solving against amplified numerical noise. A shift merely near an eigenvalue but past that check is not a problem — it converges faster, not worse.
Both functions return ConvergenceError::ZeroVector if v0 has zero norm, and ConvergenceError::NonFinite if an iterate goes NaN/infinite — errors, not panics or silent garbage output, per the crate’s no-unwrap-in-library-code policy.

Conjugate Gradient (CG)

conjugate_gradient solves the symmetric positive-definite (SPD) linear system A x = b via iterative refinement: building a sequence of conjugate directions from residuals, then stepping along each direction to get closer to the solution. CG exploits the SPD structure for rapid convergence — in exact arithmetic, it converges to the true solution in at most n iterations for an n × n matrix, but in practice finite-precision arithmetic and early stopping via residual tolerance reach an approximate solution much faster.

Each iteration does one matrix-vector product (A * p_k) and two dot products, with no factorization and only 4 fixed-size working vectors (r, p, Ap, x). This makes CG ideal for large systems where factorization is too expensive.

#![allow(unused)]
fn main() {
use rustebra::krylov::conjugate_gradient;
use rustebra::storage::StaticStorage;

// Symmetric positive-definite 2x2 matrix: [[2, 1], [1, 2]].
// System: 2x + y = 1, x + 2y = 2. Solution: x = [0, 1].
let a = StaticStorage::new([2.0_f64, 1.0, 1.0, 2.0]);
let b = [1.0, 2.0];
let x0 = [0.0, 0.0];
let mut x = [0.0; 2];
let mut r = [0.0; 2];
let mut p = [0.0; 2];
let mut ap = [0.0; 2];

conjugate_gradient(&a, 2, &b, &x0, 100, 1e-10, &mut x, &mut r, &mut p, &mut ap)
    .unwrap();

// Solution satisfies A x = b
assert!((x[0] - 0.0).abs() < 1e-8);
assert!((x[1] - 1.0).abs() < 1e-8);
}

Convergence

The convergence rate depends on the condition number κ(A) = λ_max / λ_min of a: the iteration converges to within factor (√κ - 1) / (√κ + 1) per step in the energy norm ‖e‖_a = √(e^T A e). This means:

Well-conditioned matrices (κ ≈ 1): Fast convergence, often within a handful of iterations.
Ill-conditioned matrices (κ >> 1): Slow convergence, but still graceful degradation.

For example, a matrix with κ = 100 converges at rate roughly (10 - 1) / (10 + 1) ≈ 0.82 per iteration, so even pathological cases eventually reach a solution with patience.

SPD Detection

A genuinely SPD matrix has all positive eigenvalues. CG detects a non-SPD input operationally: if the dot product p_k · (A p_k) ever becomes zero, negative, or NaN mid-iteration, the algorithm returns ConvergenceError::NonFinite. This happens automatically for:

Negative-definite or indefinite matrices
Singular or numerically singular matrices
NaN/Inf from overflow or ill-conditioning

#![allow(unused)]
fn main() {
use rustebra::krylov::{conjugate_gradient, ConvergenceError};
use rustebra::storage::StaticStorage;

// Indefinite matrix: [[1, 0], [0, -1]] has eigenvalues 1 and -1.
let a = StaticStorage::new([1.0_f64, 0.0, 0.0, -1.0]);
let b = [1.0, 1.0];
let x0 = [0.0, 0.0];
let mut x = [0.0; 2];
let mut r = [0.0; 2];
let mut p = [0.0; 2];
let mut ap = [0.0; 2];

let result = conjugate_gradient(&a, 2, &b, &x0, 100, 1e-10, &mut x, &mut r, &mut p, &mut ap);
assert_eq!(result, Err(ConvergenceError::NonFinite));
}

Gotchas

Initial guess matters: A good initial guess x0 can reduce iterations significantly. If you have no prior estimate, x0 = [0, ..., 0] is standard.
Tolerance choice: Requesting tol = 1e-15 on f64 (machine epsilon ~2.2e-16) can lead to MaxIterationsExceeded. Use tol relative to the scale of b.
Matrix must be SPD: Non-symmetric or indefinite matrices cause immediate failure. If unsure, use a general-purpose solver like QR or SVD instead.
Nonfinite iterates: If A x_k contains NaN or Inf, CG returns ConvergenceError::NonFinite on the next residual check — not a silent failure or panic, per the crate’s policy.

Example: Least-Squares via Normal Equations

While CG solves A x = b directly, you can solve overdetermined least-squares min ‖A x - b‖ via the normal equations: solve A^T A x = A^T b. The system is symmetric and positive-definite (assuming full column rank), so CG applies:

#![allow(unused)]
fn main() {
use rustebra::krylov::conjugate_gradient;
use rustebra::storage::StaticStorage;

// Overdetermined system: 3 equations, 2 unknowns.
// A = [[1, 1], [1, 2], [1, 3]], b = [2, 3, 5]
// A^T A = [[3, 6], [6, 14]] (symmetric, positive-definite)
// A^T b = [10, 24]
let ata = StaticStorage::new([3.0_f64, 6.0, 6.0, 14.0]);
let atb = [10.0, 24.0];
let x0 = [0.0, 0.0];
let mut x = [0.0; 2];
let mut r = [0.0; 2];
let mut p = [0.0; 2];
let mut ap = [0.0; 2];

conjugate_gradient(&ata, 2, &atb, &x0, 100, 1e-10, &mut x, &mut r, &mut p, &mut ap)
    .unwrap();

// Verify solution approximately
assert!(x[0] > 0.5 && x[0] < 1.5);
assert!(x[1] > 1.0 && x[1] < 2.0);
}

Note: Normal equations can amplify rounding error for ill-conditioned matrices (κ(A^T A) ≈ κ(A)²). For numerical robustness, use QR decomposition or SVD instead.

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.

#![allow(unused)]
fn main() {
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());
}
}

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. 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. 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.

Arnoldi Iteration

arnoldi reduces a general, row-major n x n matrix a to upper Hessenberg form 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 H = Qᵗ * a * Q, which is upper Hessenberg. Unlike Lanczos Iteration, a need not be symmetric — there is no three-term recurrence to exploit, so every step orthogonalizes the candidate vector against the entire basis built so far, by modified Gram-Schmidt, rather than just the two previous vectors.

#![allow(unused)]
fn main() {
use rustebra::krylov::arnoldi;
use rustebra::storage::{Basis, StaticStorage};

pub(crate) fn run() {
    println!("\n== Arnoldi iteration ==");
    // Non-symmetric on purpose: Arnoldi handles general operators, unlike Lanczos.
    let a = StaticStorage::new([4.0, 1.0, 2.0, 3.0, 3.0, 1.0, 5.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 (h, reached) = arnoldi(&a, 3, &v0, 1e-12, &mut basis, &mut scratch).unwrap();
    println!("reached = {reached}");
    for r in 0..3 {
        let row: Vec<f64> = (0..3).map(|c| h.entry(r, c).unwrap()).collect();
        println!("h[{r}] = {row:?}");
    }

    // Requesting fewer basis vectors than the matrix dimension (K < n) still produces the
    // leading block of the same Hessenberg 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_h, partial_reached) =
        arnoldi(&a, 3, &v0, 1e-12, &mut partial_basis, &mut partial_scratch).unwrap();
    println!("partial reached (K = 2) = {partial_reached}");
    println!("partial h[0][0] = {}", partial_h.entry(0, 0).unwrap());

    // A rank-1 matrix started from its own range breaks down after one vector — a good
    // outcome (the invariant subspace was found), reported as `Ok`, not an error.
    let rank_one = StaticStorage::new([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]);
    let v_rank_one = StaticStorage::new([1.0, 1.0, 1.0]);
    let mut breakdown_buffer = [0.0; 9];
    let mut breakdown_basis = Basis::<f64, 3>::new(&mut breakdown_buffer, 3).unwrap();
    let mut breakdown_scratch = [0.0; 3];
    let (_, breakdown_reached) = arnoldi(
        &rank_one,
        3,
        &v_rank_one,
        1e-10,
        &mut breakdown_basis,
        &mut breakdown_scratch,
    )
    .unwrap();
    println!("breakdown reached (rank-1 input) = {breakdown_reached}");
}
}

Orthogonalization

Modified Gram-Schmidt (subtracting each projection immediately, rather than computing all projections against the original candidate and subtracting them at the end) is the standard trade for Arnoldi: markedly more stable than classical Gram-Schmidt at the same O(K * n) per-step cost, though still less stable than Householder Arnoldi, which trades that extra stability for losing the explicit basis vectors the Krylov projection needs. No reorthogonalization pass is added on top, unlike Lanczos’s full reorthogonalization — Arnoldi’s every-vector orthogonalization does not erode as quickly as Lanczos’s three-term recurrence does, so a second pass doesn’t earn its keep the same way.

Breakdown is success, not failure

When the candidate vector’s norm falls to (numerically) zero relative to ‖a * q_j‖ after orthogonalization, q_0, ..., q_j already span an invariant subspace of a: there’s no new direction to extend the basis with, but the vectors and the leading block of H already built are exact and useful. This is reported as Ok((h, reached)) with reached < K, not an error — unlike Lanczos, which reports the analogous condition as ConvergenceError::Breakdown. The difference is what the caller does next: a Lanczos caller that requested K vectors and got fewer has nothing it can use without changing its request, while GMRES, built on top of Arnoldi, can solve directly in the smaller subspace Ok reports — the exact subspace is often exactly where the true solution already lives. Callers that do need the full K vectors distinguish this case from a complete run by checking reached < K.

Gotchas

K is a const generic on the caller’s Basis buffer, not a runtime parameter — see Krylov Basis-Size Const-Generic Convention. K > n is a DimensionMismatch.
Both ConvergenceError::ZeroVector and ConvergenceError::NonFinite on v0 are checked up front, even when K == 0 means no basis vector is ever written.
tol has no auto-computed default — see Krylov Tolerance and Convergence Criteria. A tol of 0 detects only exact breakdown.

GMRES(m)

gmres solves the general (possibly non-symmetric) linear system A x = b via restarted GMRES, GMRES(M): unlike Conjugate Gradient, A need not be symmetric positive-definite. Each restart cycle runs Arnoldi Iteration from the current residual to build an M-dimensional Krylov basis, solves the resulting small least-squares problem via Givens rotations, and updates x — restarting from the improved iterate until either the residual meets tol or max_restarts cycles are exhausted.

A is supplied as a sparse linear operator rather than a dense matrix: applying it never allocates, so restart cycles reuse the same workspace (out_x, basis, scratch) throughout.

#![allow(unused)]
fn main() {
use rustebra::krylov::gmres;
use rustebra::sparse::CsrMatrix;
use rustebra::storage::Basis;

pub(crate) fn run() {
    println!("\n== GMRES(m) ==");

    // Non-symmetric, non-SPD system, unlike Conjugate Gradient's requirements:
    // [[4, 1], [2, 3]] x = [1, 2]. Solution: x = [0.1, 0.6].
    let a = CsrMatrix::new(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![4.0_f64, 1.0, 2.0, 3.0],
    )
    .expect("valid CSR");
    let b = [1.0, 2.0];
    let x0 = [0.0, 0.0];
    let mut out_x = [0.0; 2];
    let mut buffer = [0.0; 4];
    let mut basis = Basis::<f64, 2>::new(&mut buffer, 2).unwrap();
    let mut scratch = [0.0; 2];

    gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).expect("converges");
    println!("full-basis solve: x = {out_x:?}");

    // Restart size smaller than the problem dimension (M = 1 < n = 3): the solution still
    // emerges, just carried forward across several restart cycles instead of one.
    let a3 = CsrMatrix::new(
        3,
        3,
        vec![0, 2, 5, 7],
        vec![0, 1, 0, 1, 2, 1, 2],
        vec![5.0_f64, 1.0, 1.0, 4.0, 1.0, 1.0, 3.0],
    )
    .expect("valid CSR");
    let b3 = [6.0, 6.0, 4.0];
    let x0_3 = [0.0, 0.0, 0.0];
    let mut out_x3 = [0.0; 3];
    let mut buffer3 = [0.0; 3];
    let mut basis3 = Basis::<f64, 1>::new(&mut buffer3, 3).unwrap();
    let mut scratch3 = [0.0; 3];

    gmres(
        &a3,
        &b3,
        &x0_3,
        500,
        1e-10,
        &mut out_x3,
        &mut basis3,
        &mut scratch3,
    )
    .expect("converges across restarts");
    println!("restarted GMRES(1) solve: x = {out_x3:?}");
}
}

Algorithm

Each restart cycle:

Computes the residual r = b - A x and its norm β = ‖r‖, returning Ok immediately if β <= tol.
Runs Arnoldi iteration from q_0 = r / β, building an orthonormal basis Q of up to M vectors and the upper Hessenberg projection H, stopping early (before M steps) on breakdown — an invariant subspace found before the basis filled up, the same “success, not failure” case documented on Arnoldi Iteration.
Solves min_y ‖β e_1 - H y‖ via incremental Givens rotations, then updates x <- x + Q y.

Convergence

GMRES’s residual norm decreases monotonically within a cycle (each additional basis vector can only improve the least-squares fit) and never increases across a restart, because restarting recomputes the same residual the next cycle continues from. It is not guaranteed to decrease strictly every cycle, though: a starting vector aligned with an invariant subspace the operator doesn’t expand (breakdown on the very first Arnoldi step) leaves x unchanged, and the iteration stagnates. Restarting also discards the larger Krylov subspace full (non-restarted) GMRES would have kept building, so GMRES(M) can converge slower, or stagnate on problems full GMRES would resolve — the restart budget max_restarts bounds the cost of that risk rather than eliminating it.

Gotchas

The restart size M is a const generic on the caller’s Basis buffer, not a runtime parameter — see Krylov Basis-Size Const-Generic Convention. M > n (the operator’s dimension) is a DimensionMismatch.
M == 0 never reduces a nonzero residual: no basis vector can be built, so every restart cycle is a no-op check against tol, and a nonzero residual exhausts max_restarts without ever moving x.
tol has no auto-computed default — see Krylov Tolerance and Convergence Criteria.
Non-finite residuals or Arnoldi iterates (NaN or infinite, including values that overflow f64 mid-computation) return ConvergenceError::NonFinite rather than silently producing a wrong answer.
A zero pivot in the small least-squares system built from H — a coincidental exact singularity in the projected system that Arnoldi’s own breakdown test didn’t already catch — returns ConvergenceError::Breakdown.

Eigenvalues & Eigenvectors

An eigenvector of a square matrix a is a non-zero vector v whose direction a leaves unchanged: a * v = λ * v for some scalar λ, its corresponding eigenvalue. Geometrically, a only stretches or shrinks v (by a factor of λ) rather than rotating it into a different direction. Eigenvalues show up throughout numerical linear algebra: a matrix’s condition number (see Condition Number) is the ratio of its largest to smallest singular value — themselves the square roots of the eigenvalues of aᵗ * a — and stability analysis of iterative systems generally comes down to whether relevant eigenvalues stay inside or outside the unit circle.

rustebra doesn’t compute eigenvalues via a direct, closed-form decomposition (the way it computes LU or QR). Instead, Krylov Methods provides two iterative methods: power_iteration, which converges to the dominant (largest-magnitude) eigenvalue and a corresponding eigenvector, and inverse_power_iteration, which converges to whichever eigenvalue lies nearest a caller-chosen shift — useful for targeting the smallest eigenvalue, or any eigenvalue you already have a rough estimate for. See that section for the algorithms, their convergence behavior, and worked examples.

Gotchas

These are iterative approximations, not exact decompositions — they return an error (MaxIterationsExceeded) rather than a result if convergence criteria aren’t met within the given iteration budget, and convergence speed depends on how well-separated the target eigenvalue is from its neighbors. See Power Iteration for the specifics.

Numerical Stability

This section documents two aspects of how rustebra behaves numerically: its policy on NaN/Inf propagation, and the precision trade-offs behind its fixed-iteration approximations for sqrt, sin, and cos.

NaN/Inf Policy

Outside of rustebra::krylov, the crate does not check for NaN or infinite values — it lets IEEE-754 arithmetic propagate them the way plain f32/f64 operations naturally do. An operation that receives or produces a non-finite value doesn’t return an error; the non-finite value just flows through.

#![allow(unused)]
fn main() {
let x = 0.0_f64 / 0.0;
assert!(x.is_nan());
// rustebra's arithmetic (add/sub/mul/div, dot products, norms, ...) doesn't intercept
// this — a NaN input produces a NaN output, the same as raw f64 arithmetic would.
}

Krylov methods (power_iteration, inverse_power_iteration) are the deliberate exception: they loop, so a poisoned iterate doesn’t just produce one bad result — it burns the entire remaining iteration budget computing on garbage. On embedded targets that budget is bounded and can’t be spent elsewhere, so these functions detect a non-finite iterate and stop early with ConvergenceError::NonFinite, rather than paying for the full iteration count on a chain of NaN arithmetic.

Gotchas

Don’t rely on a NaN result from most rustebra operations to signal an error condition the way Result does elsewhere in the crate — only the Krylov methods actively check for and report non-finite values. Everywhere else, a NaN/Inf input or intermediate result is expected to propagate silently, same as raw floating-point arithmetic.
Scalar::sqrt of a negative number is a deliberate exception to this exception: it returns 0, not NaN — see Scalars & Numeric Types.

Precision Guarantees

rustebra implements sqrt, sin, and cos itself rather than depending on std’s floating-point intrinsics, since the crate is no_std by default. Each uses a fixed-iteration approximation rather than a convergence-checked loop:

sqrt uses Newton-Raphson (Babylonian) iteration, x_{n+1} = (x_n + value / x_n) / 2, run for a fixed 50 iterations.
sin/cos use a Taylor series expansion around zero (after range-reducing the input to [-pi, pi]), run for a fixed 20 iterations.

In both cases the iteration count is fixed instead of stopping once successive iterates stop changing, so the amount of work done is predictable and independent of the input — important in no_std contexts, where there’s no easy way to bound worst-case iteration count for a convergence-checked loop otherwise. The trade-off is that precision isn’t uniform across the input domain: iteration counts were chosen generously enough to converge for the normal range of inputs these functions expect (vector norms and angles typical of linear algebra work), but may lose precision at the extreme ends of a type’s exponent range, where more iterations would be needed to leave the initial guess’s slow-convergence region.

#![allow(unused)]
fn main() {
use rustebra::scalar::Scalar;

let result = Scalar::sqrt(2.0_f64);
assert!((result - core::f64::consts::SQRT_2).abs() < 1e-9);
}

Comparing two floating-point results for approximate equality — accounting for the rounding error these iterative approximations (and floating-point arithmetic generally) introduce — is what FloatTolerance::epsilon() is for; see Scalars & Numeric Types.

Gotchas

These functions don’t document a specific worst-case error bound (e.g. “accurate to N ULPs”) — the guarantee is the mechanism (a fixed, generously-sized iteration budget tuned for typical inputs), not a numeric precision figure. Don’t assume bit-for-bit parity with std’s f64::sqrt/sin/cos for values at the edges of the valid range.

no_std & Embedded

rustebra is no_std by default: src/lib.rs declares #![cfg_attr(not(test), no_std)], so the crate has no dependency on the standard library (outside test builds) and no allocator requirement unless you opt in.

Without any feature flags, only stack-allocated, compile-time-sized types are available: StaticVector<T, N> and StaticMatrix<T, R, C>, whose lengths/shapes are const generics. These work on embedded targets with no heap at all.

Enabling the alloc feature pulls in extern crate alloc and unlocks heap-allocated, runtime-sized types: DynamicVector<T>, DynamicMatrix<T>, and every sparse matrix type in rustebra::sparse (CooMatrix, CsrMatrix, CscMatrix, and their sorted variants), all of which are gated behind #[cfg(feature = "alloc")] and simply don’t exist in the crate’s public API without it.

# no_std, no allocator: only StaticVector/StaticMatrix are available.
[dependencies]
rustebra = "0.3"

# adds DynamicVector/DynamicMatrix and the sparse module.
[dependencies]
rustebra = { version = "0.3", features = ["alloc"] }

Gotchas

Code that references DynamicVector, DynamicMatrix, or anything in rustebra::sparse without the alloc feature enabled won’t compile — these items aren’t hidden behind a runtime check, they’re absent from the crate entirely.
StaticMatrix::determinant on matrices larger than 4x4 returns Err(DeterminantError::MatrixTooLargeWithoutAlloc) when the alloc feature is disabled — cofactor expansion beyond that size needs heap-allocated scratch space that isn’t available in a pure no_std build.

Performance

The main performance-relevant design choice in rustebra is the split between Static* and Dynamic* types (see no_std & Embedded). StaticVector<T, N> and StaticMatrix<T, R, C> are stack-allocated with their size fixed at compile time — no heap allocation, no indirection through a Vec, and their length/shape checks (where the type system can rule out a mismatch) are eliminated entirely rather than deferred to runtime. DynamicVector<T> and DynamicMatrix<T> trade that for runtime-flexible sizing, at the cost of a heap allocation per construction and a Result-returning runtime shape check on every operation between two of them.

For algorithms with a caller-provided scratch buffer (svd, condition_number, qr on StaticMatrix), the caller controls whether that buffer lives on the stack or the heap — sizing it as a stack array keeps the whole computation allocation-free.

Beyond this structural choice, rustebra doesn’t publish benchmark numbers or make specific throughput/complexity claims — the underlying algorithms (Gaussian elimination for LU, Householder/Gram-Schmidt for QR, QR iteration for SVD) use standard, textbook time complexity for their respective operations.

Cookbook

This section collects short, task-oriented recipes for common jobs — the kind of thing you’d otherwise piece together from the reference chapters above. There are no recipes here yet; check back as the cookbook grows.

Specs Index

This index lists every spec in docs/specs/, ordered from foundational decisions to the ones that build on them. Read top to bottom for the full picture; jump directly to a spec if you already know the area you’re touching.

Crate shape and layering

The decisions everything else sits on top of: how the crate is distributed, how it’s layered internally, and what the core scalar and storage abstractions look like.

Single Crate with Incremental Feature Flags
Trait-Based Generic Layered Architecture
Generic Scalar Type
Hybrid Memory Allocation Strategy
Minimal Shared Storage Trait Between Static and Dynamic Types

Cross-cutting conventions

Conventions that apply across every module, independent of which mathematical domain they show up in.

Per-Module Error Handling with Result
NaN/Inf Policy

Scalar functions and numerical tolerance

How elementary functions are implemented and how the library decides when a computed value is small enough to treat as zero.

Public Elementary Scalar Functions, Scoped to Internal Needs
Elementary Function Precision
Numerical Tolerance for Approximate-Zero Comparisons
Auto-Computed Tolerance Defaults

Dense matrix API

Two-Level Public API for Matrix Operations

Sparse matrix API

Sparse-specific decisions, ordered from the underlying storage representation up to the operations built on it.

Sparse Index Type (u32)
Sparse Matrix Public API Shape
SortedCsrMatrix — Type-Level Column-Sort Invariant
Sparse Entry Pruning and Explicit Zero Semantics

Krylov subspace methods

Krylov Basis-Size Const-Generic Convention
Krylov Tolerance and Convergence Criteria

Single Crate with Incremental Feature Flags

Summary

The library ships as one crate, not a multi-crate workspace, and feature flags are added one at a time, only when the code they gate actually exists.

Scope

Applies to the crate’s overall distribution shape — how the mathematical domains it covers (dense linear algebra, sparse matrices, Krylov methods, and other domains added over time) are packaged, versioned, and gated behind Cargo features.

Decision

All mathematical domains live in a single published crate rather than being split across a workspace of per-domain crates. New feature flags are introduced incrementally, in the same change that introduces the code they gate — never added ahead of time in anticipation of a module that doesn’t exist yet.

A workspace split is not ruled out permanently; it remains available later if a specific domain grows large enough on its own to justify independent versioning and release cycles. It is simply not the starting point, because most domains this crate will eventually cover don’t exist yet, and designing a multi-crate boundary around code that hasn’t been written is exactly the kind of upfront speculation the project avoids elsewhere.

Constraints

A single manifest, changelog, and version number must be enough to describe the whole crate’s public surface at any point in time — callers should never need to reason about which of several crate versions is compatible with which.
A feature flag may only be added in the same change that introduces the functionality it gates. No flag is reserved or stubbed in for a domain that is still unimplemented.
Splitting into a workspace later must remain structurally possible; nothing in the current module layout should make a future split disproportionately expensive.

Status

Implemented. At present the crate carries a single feature flag, gating the optional heap-backed storage tier; further flags are expected to arrive one per domain as those domains are implemented, each introduced alongside its own code rather than in advance.

Trait-Based Generic Layered Architecture

Summary

The crate is organized into four dependency-ordered layers — scalar, storage, algorithm, and public API — each depending only on the layers before it, following the architecture pattern used by established Rust linear algebra libraries.

Scope

Applies to the overall organization of the entire crate: where a new type, trait, or function belongs, and which other parts of the codebase it is allowed to depend on.

Decision

The crate has four layers:

Scalar layer. Defines what it means to be a number the library can operate on — arithmetic and the relevant elementary functions. Depends on nothing else in the crate.
Storage layer. A minimal trait describing how to access the elements and dimensions of a vector or matrix, independent of whether the underlying memory is stack- or heap-based. Implemented by both storage strategies the crate offers.
Algorithm layer. Mathematical algorithms — arithmetic operations, norms, decompositions, Krylov methods, and so on — written generically against the scalar and storage layers. Each algorithm is implemented once and works over any type satisfying the relevant traits, regardless of which storage strategy backs it.
Public API layer. The concrete, ergonomic types a caller actually uses. These wire a specific storage strategy together with the algorithm layer, so that callers interact with concrete types rather than generics or trait bounds directly.

Each layer depends only on the layers listed before it; there is no dependency running in the other direction.

Constraints

The storage layer’s trait must stay minimal, growing only as real algorithms need new capability from it — a broad, leaky storage abstraction would undermine the entire layering.
Every mathematical algorithm belongs in the algorithm layer and must be written against the scalar and storage traits, never against a concrete storage type directly, so a single implementation always covers every storage strategy.
The public API layer must remain the only layer an ordinary caller needs to understand; generics and trait bounds from the lower layers should not leak into it unnecessarily.

Status

Implemented. The four layers are in place: the scalar and storage traits at the base, a generic algorithm layer built on top of them, and the concrete public-facing vector and matrix types assembled from that foundation.

Generic Scalar Type

Summary

Every operation in the library is written generically over a scalar type from the start, rather than hardcoded against a single concrete numeric type such as f64.

Scope

Applies to every function, algorithm, and data structure in the crate that operates on numeric values — vectors, matrices, sparse formats, and the algorithms defined over them.

Decision

Operations are written against a scalar type parameter rather than a fixed concrete type. The set of capabilities that type parameter must provide — basic arithmetic, whichever elementary functions are actually needed, identity values, and so on — is defined incrementally, as each capability is actually required by the mathematics being implemented, rather than specified upfront in full.

Constraints

The library must remain usable with f32 on resource-constrained targets that lack double-precision floating-point hardware, and with f64 elsewhere, without maintaining two separate codebases for the two cases.
Generality must be established before substantial algorithm code exists, so that switching or adding a numeric type later does not force a disruptive rewrite of function signatures across the crate.
The scalar abstraction should not grow capability the current mathematics doesn’t need; each new requirement is added only when an algorithm actually needs it, in keeping with the same discipline applied to the crate’s other core traits.

Status

Implemented. All operations across the crate are generic over the scalar type; f32 and f64 are the concrete types currently in use, with room for further scalar types (such as fixed-point) if a real need for one arises.

Hybrid Memory Allocation Strategy

Summary

The library is built around two tiers of memory use: a strict, allocation-free core that is the default build, and an optional feature flag that unlocks heap-backed dynamic data structures and algorithms for callers who have an allocator available.

Scope

Applies to every data layout and elementary operation in the crate. The baseline tier covers fixed-size vectors, buffers, and matrices whose dimensions are known at compile time. The optional tier covers growable, Vec-backed matrices and vectors, and any algorithm whose workspace size cannot be known until runtime — sparse matrix formats and iterative solvers being the primary examples.

Decision

By default, the crate builds against no allocator at all. Baseline types use const generics (<const N: usize>) to fix the size of every vector, buffer, and matrix at compile time, and live entirely on the stack.

A single conditional-compilation feature unlocks a second tier: dynamic, Vec-backed data structures and the algorithms that need a workspace whose size depends on runtime input rather than a compile-time constant. Enabling this feature is opt-in and additive — it does not change the behavior or requirements of the baseline tier.

Static (compile-time-sized) and dynamic (heap-backed) types are kept as distinct, separate type families rather than unified behind one shared public API. How much internal behavior the two families share at the implementation level is treated as its own, separate design question, not settled by this decision.

Constraints

The baseline build must run unmodified on microcontrollers with only a few kilobytes of RAM and no global memory allocator — this is the floor the strict tier is designed against, not an aspirational target.
In the baseline tier, internal memory overflow and heap-fragmentation failures must be impossible by construction (no allocator exists to fragment or exhaust), not merely avoided by discipline or testing.
Baseline, stack-allocated operations must have deterministic, bounded latency, with no allocator-dependent jitter — this rules out any hidden heap use inside the default build.
Some algorithms — sparse matrix multiplication and iterative solvers among them — fundamentally need a workspace whose size is not known until runtime; a purely stack-only design cannot express them, which is why the optional tier exists rather than trying to force every algorithm into the const-generic model.

Status

Implemented (Revised). The baseline crate remains allocator-free with const-generic, stack-only storage; the optional feature gates every Vec-backed dynamic type and every workspace-heavy algorithm. The accepted trade-offs are additional signature complexity in the baseline tier from const generics, and the maintenance of static and dynamic types as two separate, non-unified families.

Minimal Shared Storage Trait Between Static and Dynamic Types

Summary

Stack-based (compile-time-sized) and heap-based (runtime-sized) storage share a small, narrowly scoped trait so that generic algorithms can be written once and run over either kind of storage, without unifying the two families’ public-facing types.

Scope

Applies to the internal seam between the algorithm layer and the two storage strategies the crate offers. Does not apply to the public types themselves (the stack-allocated and heap-allocated vector/matrix types remain separate, independently designed types for direct use by callers) — this decision only governs how generic algorithm code accesses storage internally.

Decision

A minimal storage trait exposes only what generic algorithms actually need — element access and dimensions — and nothing more. It is not designed upfront in full; new capability is added to the trait only at the point some concrete algorithm needs it, not speculatively.

This trait exists purely as an internal seam for the algorithm layer. It does not replace or merge the public stack-allocated and heap-allocated types, which remain distinct types with their own ergonomics, aimed at their own use cases. A caller using either type directly sees no difference in their public API as a result of this trait existing.

Expensive algorithms — Krylov methods among them — are written once against this trait and work unmodified over both storage strategies, rather than requiring a duplicate implementation per strategy.

Constraints

The trait’s surface must only grow when a specific, real algorithm needs a capability it doesn’t yet expose — no capability is added speculatively ahead of a concrete caller.
The trait must not leak the dynamic-size assumptions of heap-backed storage into the compile-time-size guarantees the stack-backed storage relies on, or vice versa.
The public-facing stack-allocated and heap-allocated types must remain unaffected by this trait’s existence — callers using them directly should observe no API change.

Status

Implemented. The trait is scoped to element access and dimensions today; expensive algorithms such as Krylov methods are already written once against it and operate over both storage strategies without duplicated implementations.

Per-Module Error Handling with Result

Summary

Recoverable failures are reported through Result, never through panics, and each module defines its own error type scoped to the failures that are actually possible within it, rather than funneling every domain into one shared, crate-wide error enum.

Scope

Applies library-wide, to every module capable of failing: operations on dynamically-sized data that discover incompatible dimensions at runtime, sparse format conversions given malformed input, iterative solvers that fail to converge, and any future module with its own class of recoverable failure.

Decision

Every module that can fail defines its own error type, containing only the variants that can actually occur within that module’s own operations. A shared module holds error-related code that is genuinely common across domains — such as conversion traits or formatting helpers — but it does not define a single catch-all error type that every other module must funnel into.

No panics are used for recoverable failure. An operation that can fail for a reason the caller might reasonably want to handle returns Result, full stop.

Constraints

The library targets bare-metal and other failure-sensitive environments where an uncontrolled abort may be unacceptable to the calling application — panicking on invalid, merely-unexpected input is not an option.
A module’s error type must not force otherwise-independent mathematical domains to depend on each other’s error vocabulary; a change to one domain’s error representation must not ripple into an unrelated domain.
Error types must stay precise: a module’s error type should only contain variants that can actually occur in that module, not a large shared enum with variants irrelevant to most call sites.

Status

Implemented. Each domain module (dense operations, sparse formats, iterative solvers, and others) defines and exports its own error type; no crate-wide catch-all error type exists.

NaN/Inf Policy

Summary

The library propagates non-finite values through single arithmetic operations and elementary functions by default, following IEEE 754 semantics. Iterative Krylov solvers are the one documented exception: they detect a non-finite iterate mid-loop and return an error instead of continuing to spend iterations on a poisoned state.

Scope

Applies library-wide to scalar arithmetic and elementary functions (sqrt, sin, cos, and similar) on every scalar type. Separately, applies to the Krylov iterative solvers — power iteration, inverse power iteration, and future additions such as CG, Lanczos, Arnoldi, and GMRES(m) — where the exception described below takes over.

Decision

For a single arithmetic operation or elementary function call, a NaN or Inf operand produces a NaN or Inf result per standard IEEE 754 rules, with no special detection, no branch, and no documentation obligation at each call site — the behavior is exactly what the underlying floating-point hardware already does.

Krylov iterative solvers do not follow this rule. Each iteration explicitly checks whether the current iterate has become non-finite; if it has, the solver stops immediately and reports the failure through its error type rather than continuing to loop. This is a deliberate carve-out, not an inconsistency: a single operation always costs exactly one evaluation no matter what value it produces, but an iterative solver that keeps looping on a poisoned (NaN/Inf) state can burn its entire iteration budget converging toward nothing, producing a wrong or meaningless answer only after paying the full cost of every remaining iteration.

Constraints

The common path (a single scalar operation) must not carry extra branching cost for non-finite detection — propagation is free precisely because it does nothing beyond what hardware floating-point already does.
On resource-constrained targets, the compute budget spent inside an iteration loop is not negligible; a solver that silently iterates to max_iter on poisoned state wastes a real, bounded resource rather than just wall-clock time.
Failure reporting must go through the existing Result-based error model — no panics, no silent truncation of the iteration, and no ambiguity between “converged” and “gave up on garbage input.”
The non-finite failure must be distinguishable from other convergence failures (for example, an iterate that is finite but has no direction to normalize, or a shift that makes an inner matrix singular) so callers can tell why a solve stopped early.

Status

Implemented. Single operations and elementary functions propagate non-finite values with no extra checks. Krylov solvers detect non-finite iterates mid-loop and return a dedicated error variant distinct from other convergence-failure variants.

Public Elementary Scalar Functions, Scoped to Internal Needs

Summary

Elementary functions on the scalar type — sqrt, sin, cos, and others — are implemented from scratch, are exposed publicly, and are scoped only to what the library’s own operations actually need, rather than aiming to be a complete replacement for an external math library.

Scope

Applies to every elementary function required by the scalar trait and its implementations, and to the decision of whether and how those functions are exposed to callers of the crate.

Decision

Elementary functions are implemented internally rather than pulled in from an external dependency, consistent with the library’s preference for a minimal dependency footprint. Because these functions already have to exist for the library’s own internal use, they are also exposed as public API — a caller depending on this crate for its math operations does not need a separate math dependency for any function this crate already implements for itself.

Coverage is deliberately incomplete: only the functions the library’s own vector, matrix, and future group operations actually require are implemented. Functions with no internal caller — special functions, for instance — are out of scope and are not added purely to serve external callers who might want a general-purpose function library.

Constraints

No function is added to this surface solely for external convenience; every elementary function implemented here must have at least one internal caller within the crate’s own operations.
Precision and edge-case behavior are whatever this crate’s own implementations provide; there is no promise of matching any particular external math library’s behavior bit for bit.

Status

Implemented. Elementary functions required internally are implemented and exposed publicly; no functions exist in this surface purely for external use.

Elementary Function Precision

Summary

Precision targets for elementary functions are tiered by scalar type, not by target device: f64 targets a relative error below 1e-14, f32 targets a relative error below 1e-6, both measured over the primary domain [-2π, 2π].

Scope

Applies to every elementary function the scalar type provides — trigonometric functions, exponential, logarithm, and similar — over the domain [-2π, 2π]. Behavior outside that domain is explicitly out of scope for any precision guarantee.

Decision

Two fixed precision tiers exist, keyed only to which concrete floating-point type is in use, with no per-device variation: every target using f64 gets the same precision target as every other target using f64, regardless of whether that target is a desktop machine or a microcontroller, and likewise for f32.

No separate, reduced-precision fast path is provided for smaller or more constrained devices. That option was considered and rejected: there is no concrete evidence of demand for it, and building it would compete directly with higher-priority work (Krylov solver development) for the same engineering time.

Constraints

Outside the primary domain [-2π, 2π], elementary function behavior is documented as degraded and untested — not bounded by any stated precision guarantee. Callers operating outside this domain get no promise about the result’s accuracy.
Precision targets are a function of scalar type alone; no code path may vary an elementary function’s precision based on which target device it is compiled for.
A reduced-precision fast path must not be added without a concrete, evidenced use case — the absence of demand today is a real constraint on scope, not an oversight to be filled in speculatively.

Status

Implemented. Elementary function implementations are verified against a high-precision reference over [-2π, 2π], meeting the stated per-type targets. No fast-path variant exists.

Numerical Tolerance for Approximate-Zero Comparisons

Summary

Functions whose entire purpose is judging whether a computed floating-point value is “small enough not to count” — numerical rank, singular-value negligibility, and positive-definite rounding noise — take an explicit tolerance rather than comparing against exact zero; functions that are instead asking “is there truly nothing left, given the best available candidate was already chosen” keep exact-zero comparisons and get their real fix elsewhere.

Scope

Applies to every function that must decide whether a computed value counts as zero: rank computation, singular value decomposition and condition-number estimation (deciding whether a singular value is negligible), and Cholesky decomposition (deciding whether a value that should be non-negative before a square root is negative only due to rounding noise, or genuinely violates positive-definiteness). Also applies, by contrast, to pivot selection and reflector/normalization zero-norm checks in LU and QR decomposition, which are addressed directly rather than given a tolerance parameter. How a default tolerance is computed for callers who don’t want to supply their own is a separate concern.

Decision

Two categories of exact-zero comparison exist, and they are treated differently.

The first category — “is there truly nothing left to act on” — is an algebraic fact about the matrix, correctly decided by comparing to exact zero, provided the algorithm is already choosing the best available candidate at each step rather than merely the first nonzero one. Pivot selection and zero-norm checks fall here: their real defect was picking the first nonzero candidate instead of the largest-magnitude one, and that defect is fixed by switching to true largest-magnitude selection, after which the existing exact-zero check is correct on its own terms. No tolerance parameter is added to these functions.

The second category — “is this small enough not to count” — has no exact answer once real, non-toy floating-point data is involved. Rank, singular-value negligibility, and positive-definiteness-under-rounding-noise all fall here. Functions answering this kind of question gain an explicit tolerance parameter, compared either relative to an already-computed scale (where one exists as a natural side effect of the algorithm) or as an absolute threshold the caller chooses relative to their own data (where no such natural scale exists).

Constraints

A decomposition function (LU, QR) must remain mathematically faithful to its actual input; it must never silently approximate its result differently depending on a threshold the caller did not ask for and may not know was applied.
The “is there truly nothing left” category must never gain a tolerance parameter as a substitute for fixing genuine pivot/candidate-selection defects — thresholding would mask the real bug instead of fixing it.
The two tolerance semantics in use — relative to an existing computed scale, versus an absolute threshold — must be documented clearly enough at each function that a caller cannot mistake one for the other.

Status

Implemented. Rank, singular value decomposition, condition-number estimation, and Cholesky decomposition each take an explicit tolerance parameter, relative or absolute as appropriate. LU and QR pivot/normalization checks use true largest-magnitude selection with exact-zero comparison and take no tolerance parameter. See also [[auto-tolerance-defaults]] for how a default tolerance is offered to callers who don’t want to supply their own.

Auto-Computed Tolerance Defaults

Summary

Operations that require an explicit tolerance (rank, singular value decomposition, condition-number estimation, Cholesky decomposition) each expose a second, general-user entry point that computes a sensible default tolerance internally, so a caller who doesn’t know what tolerance to pick can still call the operation.

Scope

Applies to every operation covered by [[approximate-zero-tolerance]]’s second category — the functions that already take an explicit, caller-supplied tolerance — and to how a default value for that tolerance is computed and exposed for callers who don’t want to supply one.

Decision

Each tolerance-taking operation gets two entry points rather than a single function with an optional tolerance: a general-user name with no tolerance parameter, and the existing explicit name that takes a caller-supplied tolerance. A single function with one signature would need one trait bound covering every call, which would force the auto-computing capability’s trait requirement onto callers who already supply their own tolerance and never needed it — splitting the two paths into separate functions avoids that coupling.

Computing a default requires machine epsilon for the concrete scalar type in use. Rather than adding an epsilon method to the crate’s core scalar abstraction — which would force every scalar implementation, including hypothetical future ones with no meaningful notion of machine epsilon, to answer a question that may not apply to it — a separate, narrower trait holds it, implemented only by the concrete floating-point types. The general-user entry points require both the core scalar trait and this narrower one; the explicit, caller-supplied-tolerance entry points require only the core scalar trait, exactly as before. A future scalar type with no meaningful epsilon simply never implements the narrower trait, and loses only the auto-computing convenience — every explicit function remains callable by supplying a tolerance directly.

The default itself follows the same relative-where-it’s-free, absolute-elsewhere split used for caller-supplied tolerance: where an operation already computes a natural scale as a side effect (such as the largest singular value, or the largest-magnitude diagonal entry), the default is derived relative to that scale; where no such scale exists without adding work that wouldn’t make the comparison more correct, the default is an absolute value derived from machine epsilon and the problem size alone.

Constraints

No auto-computing capability may be added to the core scalar trait itself; it must live behind a separate, narrower trait that only floating-point-like types need implement.
A caller who already supplies their own tolerance must never be forced to satisfy a trait bound that exists only to support the auto-computing default path.
The default tolerance for a given operation must be derived either relative to a scale the operation already computes, or as an absolute value tied to machine epsilon and problem size — never an arbitrary constant unrelated to either.
A scalar type that cannot meaningfully provide machine epsilon must still be able to call every explicit, tolerance-taking function by supplying its own value; it only loses access to the general-user, default-computing entry point.

Status

Implemented. Rank, singular value decomposition, condition-number estimation, and Cholesky decomposition each expose a general-user entry point with an auto-computed default, backed by a narrow epsilon-exposing trait implemented by the floating-point scalar types, alongside their existing explicit, caller-supplied-tolerance entry points.

Two-Level Public API for Matrix Operations

Summary

Any matrix operation with more than one available algorithm exposes two levels of public API: a high-level function that picks an algorithm automatically, and one explicit function per algorithm that always runs the named algorithm regardless of input.

Scope

Applies to every matrix operation where multiple algorithms exist to compute the same result — solving a linear system, computing a determinant, decomposing a matrix — and where those algorithms differ in cost, numerical stability, or precondition (for example, one algorithm requiring symmetric positive-definite input while another has no such requirement).

Decision

For each such operation, two entry points are published:

A high-level function (for example, a general determinant or solve) that selects an algorithm automatically based on observable properties of the input — size, symmetry, detectable positive-definiteness. Which specific algorithm runs underneath is not part of the function’s public contract; it may change between versions without being a breaking change, as long as the mathematical result stays correct.
One explicit algorithm function per algorithm (for example, a cofactor-expansion determinant and an LU-based determinant) that always runs the named algorithm, regardless of what the input looks like. These functions are stable in the stronger sense: the named algorithm never changes underneath them.

Both levels are public and documented; neither is treated as internal or discouraged.

Constraints

General users who don’t want to choose an algorithm must have a single, ergonomic entry point per operation that always returns a correct result.
Mathematical users who need to study or rely on a specific algorithm’s exact behavior must be able to call it by name and trust that name never silently starts running a different algorithm.
The heuristics a high-level function uses to select an algorithm must be documented well enough that a caller can predict what will run, even though that choice is not a formal, version-stable guarantee.

Status

Implemented. Operations with multiple available algorithms — determinant and solve among them — expose both a high-level automatic entry point and named explicit-algorithm functions.

Sparse Index Type (u32)

Summary

Every sparse matrix storage format — COO, CSR, and CSC — represents row indices, column indices, and row/column pointer arrays as u32, not usize. Any value that would not fit in a u32 is rejected through the format’s own error type at the point it would be stored, rather than being truncated or silently wrapped.

Scope

Applies to the internal index storage of CooMatrix, CsrMatrix, and CscMatrix (row_indices/col_indices, row_ptr/col_ptr), and to every operation that produces or consumes those arrays: format construction, COO/CSR/CSC conversion, spmm, matmat, add, entry pruning, and validation. It does not apply to matrix dimensions themselves (rows, cols), which remain usize, or to scalar values stored in the matrix.

Decision

Index and pointer arrays are Vec<u32> (owned) and &[u32] (borrowed) throughout the sparse module. Any point where a usize value — a dimension, a running count, a computed offset — must be written into one of these arrays goes through a checked conversion (u32::try_from) or checked arithmetic (checked_add, checked_mul), and a failure there is surfaced as the format’s DimensionOverflow error variant. Constructors and conversions never truncate a usize into a u32 with an as cast where the value could plausibly exceed u32::MAX.

Constraints

Index arrays, not values, dominate memory use for sparse formats at typical sparsity levels — a CSR matrix with more rows/columns than nonzeros spends more memory on row_ptr and col_indices than on the stored values themselves. Halving the width of every index entry (u32 versus a 64-bit usize on the platforms this library primarily targets) is a direct, unconditional reduction in that footprint, which matters most on RAM-constrained targets.
No operation may panic or silently wrap on a value that doesn’t fit — every narrowing conversion from usize to u32 must be checked and reported through the crate’s existing Result-based error model.
u32::MAX (4,294,967,295) rows or columns is far beyond what a sparse matrix on a memory-constrained target could hold in the first place, so the narrower range is not a realistic ceiling in practice, while the memory savings apply to every matrix regardless of size.

Status

Implemented. All format constructors, conversions, and sparse operations (spmm, matmat, add, prune, validate) that write into an index or pointer array go through a checked path and return DimensionOverflow on failure instead of wrapping.

Sparse Matrix Public API Shape

Summary

Sparse matrix operations (COO, CSR, CSC) are exposed as free functions rather than methods, COO is a build-only format with no arithmetic of its own, shape failures across every sparse operation share one error type, matvec/matmat return owned allocations, and a narrow trait lets iterative solvers call into any sparse format without knowing which one they’re holding.

Scope

Applies to the three sparse storage formats — coordinate (COO), compressed sparse row (CSR), and compressed sparse column (CSC) — and to every operation defined over them: construction, conversion between formats, matrix-vector and matrix-matrix products, addition, scaling, and the interface iterative solvers use to call into a sparse matrix without depending on its concrete format.

Decision

Free functions, not methods. All sparse operations are free functions, following the same convention used throughout the algorithm layer elsewhere in the crate. This keeps operation logic decoupled from the format type definitions, so operations can be added incrementally without reopening format modules — and stays consistent with the crate’s avoidance of trait-object-based dynamic dispatch, which a method surface on the types themselves would not by itself require, but which the free-function convention sidesteps entirely regardless.

COO is build-only. COO is the natural way to build a sparse matrix incrementally — it requires no ordering invariant, and duplicate entries at the same position are valid and expected while building. Arithmetic directly on COO is either expensive or semantically ambiguous (adding two COO matrices without first sorting both has no efficient implementation, and naively concatenating entry lists produces a result with unresolved duplicate positions that behaves differently from the equivalent CSR operation). COO supports construction and entry-by-entry manipulation; arithmetic requires converting to CSR or CSC first.

One shared shape-mismatch error. Every sparse operation that can fail due to incompatible dimensions — matvec’s vector-length-versus-column-count check, binary operations’ equal-dimension check, matmat’s inner-dimension-agreement check — returns the same single error type rather than a module-specific one. All of these failures represent the same underlying fact from a caller’s perspective, and a single type means one import and one ?-target regardless of which sparse operation is being called.

Owned output. Matrix-vector and matrix-matrix products return an owned, newly allocated vector or matrix. Every sparse format already requires an allocator to exist at all, so there is no no_std-without-alloc path through this API to preserve, and an owned return keeps call sites simple.

Solver interface without dynamic dispatch. Iterative solvers need only a matrix-vector product and the operand dimensions, and should not be forced to pick a specific format at the call site. A trait exposing rows, cols, and an apply-style matrix-vector product provides this indirection through generic, monomorphized dispatch — no trait objects, no heap allocation for the indirection itself. Every sparse format type implements it, and solver functions are written generically against it rather than against any one format.

CSR/CSC parity. Both formats expose the same operation set (scale, matvec, add, matmat), because a caller who receives a matrix in one format should not be forced to convert before using it just because an operation happened to be implemented for the other format first.

Constraints

A shape-mismatch failure anywhere in the sparse module must surface through the same error type, regardless of which operation detected it.
The solver-facing trait must add no heap allocation of its own beyond what the concrete sparse format already requires, and must not rely on trait objects.
Any per-iteration matrix-vector product called from inside a solver’s loop must be able to write into a caller-supplied output buffer rather than allocate fresh on every call — a per-iteration allocation is both a throughput cost and, on real-time embedded targets, a determinism problem.
COO must never be extended with arithmetic operations that produce ambiguous results (such as unresolved duplicate entries after an add); any operation with that risk requires converting to CSR or CSC first.

Status

Implemented. All sparse operations are free functions; COO supports construction and conversion only; a single shared error type covers every shape mismatch; the solver-facing trait writes into a caller-supplied buffer for its matrix-vector product; CSR and CSC expose matching operation sets.

SortedCsrMatrix — Type-Level Column-Sort Invariant

Summary

A newtype wrapper around the CSR format carries a compile-time guarantee that every row’s column indices are stored in ascending order, so algorithms that require this ordering can state the requirement in their signature instead of checking or re-sorting at runtime.

Scope

Applies to any CSR matrix produced by an operation that already naturally yields sorted column indices per row, and to any current or future algorithm — sparse triangular solve and sparse Cholesky factorization among them — that requires that ordering as a precondition.

Decision

A wrapper type around the plain CSR format restricts construction so that only genuinely sorted matrices can exist as values of the wrapped type. Two construction paths exist: a public constructor that sorts an arbitrary CSR matrix’s rows in place before wrapping it, safe for any input; and a second, module-internal constructor that asserts the invariant without independently verifying it, used only by the specific conversion and multiplication operations that already guarantee sorted output as a natural consequence of how they work (so no separate sort pass is needed in those two cases).

The wrapper transparently exposes every read-only method of the underlying plain CSR type through automatic dereferencing, so existing code that reads fields or calls accessors on a result that happens to be sorted needs no changes. A consuming conversion lets a caller discard the sorted guarantee explicitly when they need to pass the value to code that only accepts the plain format.

Operations that already produce sorted output by construction return the sorted wrapper type directly, communicating the guarantee to their callers at the type level rather than by convention or documentation alone. Ordered COO → CSR → CSC by format/flow rather than by release history (all six ship together, so no chronological ordering applies), these are: coo_to_csr, csr_to_csc, csc_to_csr, spmm, add_csr, and add_csc.

Constraints

The module-internal, invariant-skipping constructor must only be reachable from code that can actually guarantee sorted output as a structural consequence of its own algorithm — never exposed as a general-purpose “trust me” escape hatch.
The wrapper must add no runtime cost beyond what plain construction already costs; the guarantee is enforced once at construction, not re-checked on every subsequent read.
Existing code operating on the plain format’s read-only accessors must continue to compile unchanged when handed a value of the wrapped type.

Status

Implemented. The wrapper type is used as the return type of the six operations that produce sorted output as a natural consequence of their algorithm — coo_to_csr, csr_to_csc, csc_to_csr, spmm, add_csr, and add_csc — and a public, verifying constructor is available for wrapping an arbitrary CSR or CSC matrix from any other source.

Sparse Entry Pruning and Explicit Zero Semantics

Summary

A pruning function drops sparse-matrix entries whose magnitude falls at or below a caller-supplied threshold, using only ordering comparisons rather than requiring the scalar type to provide an absolute-value method, and the threshold is always supplied by the caller — no auto-computed default exists for this operation.

Scope

Applies to CSR sparse matrices accumulating numerically negligible entries over the course of iterative computation — fill-in from products, cancellation from additions, ordinary floating-point rounding — where those entries are structurally present but too small to be worth the memory and arithmetic cost of keeping.

Decision

An entry with value v is dropped when its magnitude is at or below a threshold τ > 0 (equivalently, it is kept when v is strictly greater than τ or strictly less than -τ). This matches the standard definition of threshold-based pruning used elsewhere in numerical software.

Because the scalar type deliberately exposes no absolute-value or ordering methods of its own — keeping its required surface minimal rather than assuming every possible scalar implementation has a meaningful notion of magnitude — the pruning function adds the ordering bound it needs directly to its own signature instead. The magnitude comparison is then expressed without an absolute-value method at all, by comparing against both the threshold and its negation using only ordering and the scalar type’s existing arithmetic.

The threshold is always supplied by the caller. Unlike operations where a machine-epsilon- derived default is a sensible baseline, the “right” pruning threshold depends entirely on the caller’s own domain — the physical scale being modeled, the acceptable accuracy loss downstream, whether the goal is memory compression or cancellation-noise removal — and no neutral default is correct across those cases. A caller who wants no pruning at all simply passes a zero threshold.

The function returns the plain CSR format, not the sorted-invariant wrapper type, even though pruning only removes entries and never reorders the ones that remain: if the input was already known to be sorted, the output is too, but expressing that in the return type would require a second function or an overload, which is not justified without a concrete caller needing it. A caller who needs the sorted wrapper type after pruning applies the public sorting constructor to the result.

Only the row-oriented format is covered. The column-oriented equivalent is not implemented, because no current caller needs it — it would be structurally identical with rows and columns swapped, and is deferred until a real use case exists rather than added speculatively.

Constraints

The pruning function must not require any capability beyond ordering comparisons and the scalar type’s existing arithmetic — no absolute-value method is added to the core scalar abstraction to support it.
No auto-computed default threshold is offered; every call must supply an explicit value, because no domain-independent default is correct.
Pruning must never reorder surviving entries — only remove entries, preserving whatever ordering the input already had.

Status

Implemented for the row-oriented format. The column-oriented equivalent is deferred until a concrete caller needs it; callers working in that format today convert to the row-oriented format, prune, and convert back.

Krylov Basis-Size Const-Generic Convention

Summary

Krylov subspace methods that build up a basis of a fixed maximum size expose that size as a const K: usize (Lanczos, Arnoldi) or const M: usize (GMRES(m) restart size) type parameter directly on the public function, following the notation standard in the numerical linear algebra literature.

Scope

Applies to any current or future Krylov subspace method whose basis (or restart) size is known at compile time — Lanczos and Arnoldi iterations, and the restart parameter of GMRES(m). Does not apply to power_iteration/inverse_power_iteration, which have no basis of growing size to bound. Does not introduce a new trait; it extends the existing minimal [[static-dynamic-storage-trait]] the same way every other generic algorithm in the crate already does.

Decision

const K: usize names the Lanczos/Arnoldi basis size, and const M: usize names the GMRES(m) restart size — the same symbols the literature uses for each, rather than a single shared name across both algorithm families or a crate-invented one.

This size is a direct const-generic type parameter on the public function signature, not hidden behind a wrapper type. Matrix types already expose their dimensions as const generics throughout the public API (StaticMatrix<T, const R: usize, const C: usize>), so a basis size expressed the same way is consistent with how callers already read every other compile-time size in this crate; wrapping it would introduce a second convention for the same kind of information with no offsetting benefit.

The const generic is an extension of the existing narrow Storage trait used to write generic algorithms once against both static and dynamic storage, not a new trait — consistent with the crate’s one-pattern-reused-across-algorithms approach already established for Krylov and other generic numerical code.

Constraints

K and M name the basis/restart size specifically for Lanczos/Arnoldi and GMRES(m) respectively; a different Krylov method with a differently-named parameter in its own literature uses that name instead, rather than being forced onto K/M for consistency’s own sake.
The basis/restart size is exposed as a direct type parameter on the public function, never hidden behind a wrapper type, matching how matrix dimensions are already exposed elsewhere in the public API.
No new trait is introduced to support this; any capability the algorithms need beyond what the existing storage trait exposes is added to that trait only when a concrete algorithm needs it, not speculatively.

Status

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.

Krylov Tolerance and Convergence Criteria

Summary

Krylov iterative methods take tol (convergence) and, for inverse iteration, singular_tol (shifted-matrix singularity) as required caller-supplied parameters, with no auto-computed default. Convergence is only declared once two independent measures — eigenvalue stabilization and eigenvector-residual stabilization — both fall within tol.

Scope

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.

Decision

tol, and singular_tol for inverse iteration, are required parameters with no general-user, default-computing entry point, unlike the algorithm::matrix category-2 functions. Those functions can derive a default from a scale the algorithm already computes (the largest singular value, the largest-magnitude diagonal entry) or from machine epsilon and problem size alone. A Krylov solver has no equivalent problem-independent quantity before it starts: its natural scale is the eigenvalue estimate being refined, which doesn’t exist until at least one iteration has run. Callers choose tol relative to their own problem scale and the accuracy their use case needs.

Convergence requires two measures to independently fall within tol (relative to the current eigenvalue estimate) before an iterate is accepted:

eigenvalue stabilization: |λ_k - λ_{k-1}| <= tol * |λ_k|
eigenvector-residual stabilization: ‖A·v_k - λ_k·v_k‖ <= tol * |λ_k|

Neither measure is sufficient alone. An eigenvalue estimate can plateau prematurely while the eigenvector direction is still rotating toward its limit, and a small residual computed from an estimate that hasn’t yet stabilized can understate how much the eigenvalue itself is still moving. Requiring both closes each measure’s blind spot with the other. Because the eigenvalue-stabilization check needs a prior estimate to compare against, at least one full iteration must run before convergence can be declared: max_iter < 2 can never converge.

inverse_power_iteration additionally takes singular_tol, governing how close the shifted matrix a - shift * I may be to singular before the solver reports SingularShift instead of continuing to iterate on amplified noise. See [[nan-inf-policy]] for how SingularShift relates to the NonFinite and ZeroVector failure modes of the same solvers.

condition_number is not consumed internally by the Krylov solvers as a preconditioning signal; it exists purely as a caller-facing diagnostic for deciding, before or after a solve, how much precision loss to expect.

Constraints

tol and singular_tol never get an auto-computed default; every Krylov entry point requires the caller to supply them directly.
Convergence must never be declared on a single measure — both eigenvalue stabilization and residual stabilization must hold before an iterate is accepted.
max_iter < 2 must never report convergence, since the eigenvalue-stabilization check requires a prior estimate that doesn’t exist before the first iteration completes.
condition_number remains a caller-facing diagnostic only; no Krylov solver may consume it internally as a preconditioning signal.

Status

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. 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.

Keyboard shortcuts

rustebra