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);
}

Keyboard shortcuts

rustebra