mtl5-python

Python bindings for MTL5 — NumPy/SciPy/JAX/PyTorch interop with hardware accelerator dispatch.

Built with nanobind for minimal overhead and zero-copy array interop.

Install

pip install mtl5

Prebuilt wheels are published for Linux, macOS, and Windows on CPython 3.10–3.12, so no compiler is required. Optional ecosystem integrations are available as extras:

pip install "mtl5[scipy]"   # SciPy sparse interop
pip install "mtl5[all]"     # scipy, torch, jax, pandas, scikit-learn

From source

Building the extension yourself — for contributing, or to enable a build option — needs Python 3.10+ and a C++20 compiler (GCC 12+, Clang 15+, MSVC 2022):

pip install .

See Development for an editable install.

Upgrading to 5.10.3

5.10.3 changes behaviour. mtl5.vector() and mtl5.matrix() no longer silently convert their input, so code that worked on 5.10.2 can raise TypeError. You are affected if you pass a non-contiguous array (a[::2], a column, a transpose) or a dtype outside float32, float64, int8, int16, int32, int64, uint8, complex64, complex128.

mtl5.vector(np.ascontiguousarray(a))  # keeps float64
mtl5.vector(a.astype(np.float64))  # or convert deliberately

That is a fix rather than a restriction: the old path silently downcast a float64 slice to float32 — losing precision — while reporting is_view=True and aliasing nothing. Full explanation in the v5.10.3 release notes.

Quick start

import numpy as np
import mtl5

# Vectors and norms
v = mtl5.vector(np.array([3.0, 4.0]))
print(mtl5.norm(np.array([3.0, 4.0])))  # 5.0

# Dot product
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
print(mtl5.dot(a, b))  # 32.0

# Solve Ax = b
A = np.array([[2.0, 1.0], [1.0, 3.0]])
b = np.array([5.0, 7.0])
x = mtl5.solve(A, b)
print(x)  # [1.6, 1.8]

Passing containers to NumPy and scikit-learn

The dense containers implement __array__, so anything that takes array-like input through np.asarray accepts them, scikit-learn included:

from sklearn.linear_model import LinearRegression

X, y = np.random.rand(100, 3), np.random.rand(100)
model = LinearRegression().fit(mtl5.matrix(X), mtl5.vector(y))
np.shares_memory(np.asarray(mtl5.matrix(X)), X)  # True: a view, not a copy

• Universal containers (posit16, fp16, ...) convert to float64 with a copy, so np.asarray(P, copy=False) raises. (That is NumPy 2; NumPy 1.x's np.asarray has no copy=.) scikit-learn then computes on the rounded values in float64.
• Sparse matrices have no __array__, so they aren't densified behind your back. Pass mtl5.sparse.to_scipy(S).
• A float32 container is upcast to float64 by an estimator that takes either float type, such as Ridge, because .dtype is the name "f32", not a NumPy dtype. np.asarray(M) first is zero-copy and keeps float32.

applications/sklearn_pipeline.py compares the test error of f64, f32, posit16 and fp16 feature storage.

Mixed precision

A mixed-precision operation has three independent precisions. MTL5 supplies the precision-generic kernels, Universal supplies the number systems, and this package composes them:

	chosen by
element (storage)	the container you pass in — mtl5.convert(x, "posit16")
accumulator (compute)	accumulator=
result (delivery)	result=

import numpy as np, mtl5

x = np.random.default_rng(0).standard_normal(4000)
v = mtl5.convert(x, "posit16")  # store narrow

mtl5.mixed.dot(v, v)  # accumulate in posit16 too
mtl5.mixed.dot(v, v, accumulator="f64")  # ...or in double
mtl5.mixed.dot(v, v, accumulator="quire")  # ...or exactly, in Universal's quire

Accumulating 4000 posit16 products, against the exact value of the same posit16 data:

accumulator	relative error
none (posit16)	1.8 × 10⁻¹
"f32"	3.8 × 10⁻⁷
"f64"	1.1 × 10⁻¹⁶
"fma"	1.1 × 10⁻¹⁶
"quire"	0 — bit-exact

mtl5.mixed.accumulators(dtype) lists what a given element type supports. The quire is available for f32, f64 and the posit, cfloat, lns and fixpnt families. Exactness varies by family. It is exact for f32, f64, posit, cfloat and fixpnt. On f32/f64 the result is the exact dot product rounded once, to nearest with ties to even, so it is correctly rounded. It is also slow: every product goes into a 584-bit (f32) or 4226-bit (f64) fixed-point register, and in mp-blas's level-3 benchmark a quire<double> accumulate runs about 3900× slower than a float64 one. On f32/f64, accumulator= is available on dot, norm and frobenius_norm. The lns quire's quire_mul forms each product through double and rounds it, so it is not exact and can do worse than the "f64" accumulator (universal#1203; see Universal's lns/fdp.hpp). Universal supplies the quire's accumulator functions itself; there is no bridge header in this package.

accumulator= is available on dot, norm (ord=2), frobenius_norm, matvec and matmul.

Integer operands: accumulator="i32"

The same idea with integers, which is where the hardware is. int8 and int16 operands accumulated in int32 map onto the widening multiply-accumulate — vpmaddwd / vpdpbusd on x86, SMLAL / SDOT on NEON:

a = np.random.randint(0, 256, 4096, dtype=np.uint8)  # activations
w = np.random.randint(-128, 128, 4096, dtype=np.int8)  # weights
mtl5.mixed.dot(a, w, accumulator="i32")

u8 × i8 is VNNI's native pairing on x86; ARM implements the symmetric ones first. Either order works — a dot product is symmetric, so MTL5 swaps the operands onto whatever the machine has.

The sum wraps, and sooner than vector length suggests. Products are always exact, but headroom goes as operand magnitude: about 2^(31-2b) terms at b bits. Measured at full range, one int16 × int16 product uses 2³⁰ of the int32 range so two overflow it, while int8 × int8 holds 131071 terms. That gap is why quantized inference is 8-bit. The wrap is two's complement, hence bit-identical across lane counts, backends and thread counts.

accumulator="i32" is required for these dtypes — the default is element precision, and an 8-bit accumulator overflows almost immediately. Build the arrays with NumPy: convert() does not target integers, because re-quantizing reals into 8 bits is a quantization scheme rather than a cast.

Matrix products come in two kernels, chosen by name:

A = np.random.randint(0, 256, (64, 512), dtype=np.uint8)  # activations
W = np.random.randint(-128, 128, (512, 32), dtype=np.int8)  # weights
mtl5.mixed.matmul(A, W, accumulator="i32")  # widen on load -> DenseMatrix_i32
mtl5.mixed.matmul_quad(A, W)  # quad multiply-accumulate (VNNI / SDOT)

An integer matmul returns a matrix of the accumulator's type, not the element type: i8 @ i8 is a DenseMatrix_i32. Both kernels give the same result bit for bit, equal to the naive triple loop at any size, because integer addition mod 2³² is associative. Only their speed differs. Upstream measures the quad kernel at 1.27×–3.64× over widen-on-load, against about 1.2× for the same instruction on a dot, because a GEMM is compute-bound.

pairing	matmul(..., accumulator="i32")	matmul_quad
i8 × i8, u8 × i8	DenseMatrix_i32	DenseMatrix_i32
u8 × u8	DenseMatrix_i32	DenseMatrix_u32
i16 × i16	DenseMatrix_i32	—
i8 × u8	DenseMatrix_i32	refused

matmul_quad refuses i8 × u8 rather than reordering it. The hardware op takes u8 × i8, and unlike a dot product a GEMM is not symmetric, so no swap of the operands fixes it for you. The quad kernel is never picked automatically, because nothing in a result shows which of the two kernels ran.

A released wheel is built at the x86-64 baseline, so it gets int8's bandwidth win — one byte per element against float64's eight — but not the VNNI instruction. mtl5.build_info() says which you have:

info = mtl5.build_info()
info["quad_dot"]  # 'none', 'partial' or 'all'
info["quad_dot_native"]  # {'u8xi8': ..., 'i8xi8': ..., 'u8xu8': ...}

"partial" is the usual answer on hardware that has the instruction at all. x86 implements u8 × i8 first (AVX3_DL) and the symmetric pairings only at AVX10.2. ARM is the other way round: i8 × i8 and u8 × u8 come with DotProd, and u8 × i8 only with I8MM. Build with -C cmake.define.MTL5_NATIVE_ARCH=ON to reach the instruction, on CPUs Highway will target. Alder Lake has AVX-VNNI but no AVX-512, so it still reports "none".

Iterative refinement

Factor cheaply in a low precision, then recover accuracy with a residual formed in float64:

x, info = mtl5.mixed.lu_iterative_refine(A, b, working="posit16", rel_tol=1e-14)
# info -> {'iters': 7, 'rel_residual': ..., 'converged': True}
mtl5.mixed.backward_error(A, x, b)

On a 100×100 system, forward error of the refined solution:

working precision	iterations	forward error
fp16	5	1.9 × 10⁻⁸
posit16	7	1.9 × 10⁻¹⁵
f32	2	3.8 × 10⁻¹⁶
f64	0	3.1 × 10⁻¹⁶

The result is always the best iterate found, so an over-long max_iter never degrades the answer.

mtl5.mixed.iterative_refine(A, M, b) is the sparse counterpart, refining through any factorization exposing solve() — the sparse direct factorizations below, or the ILU(0)/IC(0) preconditioners.

Arrays from universal_dtypes

universal_dtypes makes the same number systems real NumPy dtypes. Its arrays go into MTL5 and come back out bit for bit, without passing through float64:

import numpy as np, universal_dtypes as ud, mtl5  # pip install mtl5[universal]

a = np.array([1.0, 2.0, 3.0], dtype=ud.posit16)
v = mtl5.vector(a)  # DenseVector_posit16, the same 16-bit encodings
mtl5.norm(v)
v.to_universal()  # ndarray[posit16] again; to_numpy() still gives float64

mtl5.matrix() and the _copy forms work the same way. Supported: posit8, posit16, posit32, posit64, fixpnt8, fixpnt16, lns16, lns32, takum32, fp16 and the dd/td/qd cascades. These are the formats whose encoding matches this build. The test suite checks that every 8- and 16-bit pattern decodes the same in both libraries. Two dtypes are refused with a TypeError rather than misread:

• fp8e5m2, because MTL5's fp8 is e2m5;
• the formats MTL5 doesn't bind.

The data is copied, not viewed. The two packages are built on different Universal releases, and those releases share the encoding but not necessarily the C++ object layout. The named factories (mtl5.vector_posit16(x)) are unchanged: they round float64 values into a format, which is a different job.

Matrix Market I/O and sparsity pictures

A = mtl5.io.mm_read("circuit.mtx")  # -> CSR
D = mtl5.io.mm_read_dense("small.mtx")  # -> dense
mtl5.io.mm_write_sparse("out.mtx", A, "comment")

mtl5.io.spy(A, "pattern.png")  # binary non-zero pattern
mtl5.io.spy_magnitude(A, "mag.png", log_scale=True)
mtl5.io.spy_density(A, "dens.png", max_pixels=512)

Two things differ from what a SciPy user expects.

The function picks the container, not the file. scipy.io.mmread returns an ndarray for an array file and COO for a coordinate one. Here mm_read always gives CSR and mm_read_dense always gives dense, and both accept either file format — so reading a dense .mtx into CSR is a deliberate call rather than an error.

The PNGs are uncompressed. MTL5 writes them with a from-scratch encoder emitting DEFLATE stored blocks, which is what lets spy work with no image library and no plotting stack in the process. File size is therefore about width × height × channels — measured at the default max_pixels=1024, that is ~1.0 MB for spy (grayscale) and ~3.1 MB for spy_magnitude / spy_density (RGB). Pipe the output through any PNG optimizer if that matters.

.gz inputs are read transparently only when MTL5 is built with zlib (-C cmake.define.MTL5_WITH_ZLIB=ON); mtl5.build_info()["zlib"] reports it, and without it a .gz path raises rather than silently mis-parsing.

Test matrices

A named catalog of matrices with known pathologies — the inputs a mixed-precision experiment actually wants, and tedious to hand-roll correctly:

g = mtl5.generators

g.hilbert(8)  # cond ~1.5e10 — the canonical ill-conditioned matrix
g.clement(7)  # eigenvalues exactly -6,-4,-2,0,2,4,6
g.randspd(4, [1, 10, 100, 1000])  # SPD with exactly that spectrum
g.randsym(4, [-5, -1, 2, 7])  # controlled *indefinite* matrix
g.randsvd(20, 20, kappa=1e6)  # condition number exactly 1e6
g.laplacian_2d(64, 64)  # sparse, the usual solver benchmark

Also frank, pascal, wilkinson, rosser, magic, lehmer, lotkin, minij, ones, forsythe, kahan, moler, companion, vandermonde, randorth, laplacian_1d, poisson2d, and the published catalog via testsuite_names() / testsuite_kappa(name).

Every dense generator takes dtype=, which is what makes them useful here:

H = g.hilbert(8, dtype="posit16")  # correctly rounded posit16 Hilbert
mtl5.cholesky(H)  # ...then watch it fail

Generation happens in float64 and dtype= rounds. That is the right semantics for a test matrix — the definitions are over the reals, so you want the correctly rounded representation of the exact entry, not the result of evaluating the formula in low-precision arithmetic.

Range vectors follow NumPy, and take dtype= too:

mtl5.linspace(0, 1, 5)
mtl5.arange(0, 10, 3)
mtl5.logspace(0, 3, 4)
mtl5.geomspace(1, 1000, 4)  # 1, 10, 100, 1000 — a true geometric
# progression, not logspace's exponents

Dense factorizations

qr = mtl5.qr(A)  # Householder QR; tall or square
x = qr.solve(b)  # least squares
qr.Q, qr.R

lq = mtl5.lq(A)  # the row-space counterpart
lq.L, lq.Q

ld = mtl5.ldlt(A)  # symmetric, indefinite allowed
ld.solve(b)
ld.diagonal()  # D — its signs are the inertia

All accept an MTL5 matrix or a float32/float64 NumPy array, alongside the existing mtl5.lu and mtl5.cholesky.

Every dense factorization — lu, qr, lq, cholesky, ldlt and bunch_kaufman — is available for float32, float64 and all fifteen Universal dtypes, which is what makes a like-for-like comparison across number systems possible. The integer element types are not supported.

Availability is not usability, though, and the two diverge sharply at narrow widths. The orthogonalizing pair (qr, lq) fails first: below 16 bits every Householder reflector rounds to zero, so Q comes back as the exact identity and the factorization looks perfect while reconstructing nothing. The solving three (cholesky, ldlt, bunch_kaufman) have no orthogonality to lose and stay usable to 8 bits. tests/test_universal_factorizations.py records where each format stops.

Cholesky vs LDLᵀ across number systems

ldlt and cholesky are both available in every bound precision, which is what makes the interesting comparison possible. Cholesky takes square roots, so it refuses a matrix that has drifted out of positive-definiteness — the failure mode of a Kalman covariance update in low precision. LDLᵀ has no square roots, survives, and records what happened in D:

P = np.eye(6)
P[3, 3] = -1e-3  # covariance went indefinite

mtl5.cholesky(mtl5.convert(P, "posit16"))  # RuntimeError: not SPD
d = mtl5.ldlt(mtl5.convert(P, "posit16")).diagonal()
(d < 0).any()  # True — D names the bad direction

mtl5.bunch_kaufman(A) is the pivoting variant, for a symmetric matrix that plain ldlt rejects on a zero pivot.

Eigenvalues, BLAS 2/3, and matrix properties

The eigen entry points mirror numpy.linalg, and return NumPy arrays:

mtl5.eigvalsh(A)  # symmetric eigenvalues, real, ascending
w, Q = mtl5.eigh(A)  # ...with eigenvectors:  A = Q diag(w) Qᵀ
mtl5.eigvals(A)  # general spectrum, complex
w, V = mtl5.eig(A)  # ...with right eigenvectors
mtl5.spectral_radius(A)
mtl5.inertia(A)  # {'positive': …, 'negative': …, 'zero': …}

BLAS levels 2 and 3 write into a caller-supplied output, as BLAS does — that in-place accumulation is the point:

mtl5.ger(alpha, x, y, A)  # A += alpha x yᵀ
mtl5.symv(alpha, A, x, beta, y)  # y = alpha A x + beta y
mtl5.trsv(A, x, upper=True)  # x = A⁻¹x
mtl5.trmm(alpha, A, B, upper=True)  # B = alpha A B
mtl5.trsm(alpha, A, B, upper=True)  # solve A X = alpha B
mtl5.symm(alpha, A, B, beta, C)  # C = alpha A B + beta C
mtl5.syrk(alpha, A, beta, C)  # C = alpha A Aᵀ + beta C
mtl5.syr2k(alpha, A, B, beta, C)

Property predicates come in two cost classes. The docstrings are authoritative; the split is:

• O(n²) or cheaper — the structural checks (is_square, is_empty, is_symmetric, is_hermitian, is_triangular and the upper/lower variants, is_diagonal, is_banded, is_diagonally_dominant) and every vector predicate.
• O(n³) — anything that forms a product, factorizes, or runs an eigensolve: is_orthogonal, is_unitary, is_normal, is_spd, is_positive_definite, is_singular, is_nonsingular, is_invertible, spectral_radius, inertia, is_indefinite. Don't put these inside a loop.

SVD and the queries built on it:

U, s, V = mtl5.svd(A)  # s is the vector of singular values (NumPy's convention)
mtl5.svdvals(A)  # singular values only — cheaper
mtl5.condition_number(A)  # σ_max / σ_min
mtl5.rcond(A)  # σ_min / σ_max, safer near singular
mtl5.numerical_rank(A)
mtl5.nullity(A)

svd takes a tol. V and the reconstruction are accurate to machine precision regardless, but U's orthogonality is bounded by the iteration's tolerance rather than by eps. Over 160 matrices (n = 3..20, four tolerances) ‖UᵀU − I‖ had a median of about 1.4×tol and a worst case of 4.5×tol; the regression test asserts 20×tol to leave headroom. Treat tol as the knob: tighten it if you need an orthonormal U specifically, rather than relying on a particular multiple.

Complex numbers

mtl5.vector() and mtl5.matrix() accept complex64 and complex128 arrays and give the same zero-copy views as the real types, with dtypes c64/c128:

A = mtl5.matrix(np.array([[2 + 0j, 1 - 1j], [1 + 1j, 3 + 0j]]))
b = mtl5.vector(np.array([3 + 1j, 1 + 4j]))
x = mtl5.solve(A, b)  # complex LU with partial pivoting
mtl5.norm(x, 2)  # a real float, not a complex
A.real.to_numpy()  # real part, as a real matrix

Three things differ from the real case, and getting them wrong is quiet rather than loud, so they are worth stating.

dot is Hermitian. It computes sum(conj(a[i]) * b[i]), conjugating the first argument — that is NumPy's vdot, not NumPy's dot. The unconjugated bilinear product is mtl5.dot_real, which is what np.dot does for 1-D complex. Both exist because both are wanted; the names say which is which.

.T does not conjugate. MTL5's transpose is the plain one. .H (or mtl5.adjoint) is the conjugate transpose. For real elements the two coincide, which is exactly why the distinction has to be explicit here.

Hermitian and symmetric are different properties. mtl5.is_hermitian(A) tests A == Aᴴ and mtl5.is_symmetric(A) tests A == Aᵀ; for complex those are not the same matrix, and which one you have decides which solver is right.

What is available: containers and factories, all four norms, dot/dot_real, matmul/matvec, solve/lu/inv, transpose/adjoint/conj, ldlt_solve, cholesky, and qr/lq. What is not: bunch_kaufman, the eigen and SVD family, and the Krylov solvers — MTL5 has no complex implementation of any of them, and complex input raises a TypeError naming the alternative rather than silently taking a real part.

Complex least squares works through qr, which uses MTL5's complex Householder:

f = mtl5.qr(A)  # A complex, num_rows >= num_cols
x = f.solve(mtl5.vector(b))  # min ||Ax - b||_2
f.Q, f.R  # Q is unitary: Q^H Q = I

lq is the row-space counterpart, for the underdetermined case. Both were checked for the answer rather than the compile: the least-squares residual is orthogonal to range(A) to 3.5e-15, which is what distinguishes a solve applying Qᴴ from one applying Qᵀ — the latter would reconstruct A = QR perfectly and still solve the wrong problem.

Symmetric and Hermitian need different factorizations, and mtl5.ldlt_solve picks from the matrix: LDLᵀ for a complex symmetric one, LDLᴴ for a Hermitian one. They are not interchangeable — sending either input to the other's kernel gives a wrong answer — so a matrix that is neither is refused rather than guessed at.

S = mtl5.matrix(np.array([[2 + 1j, 1 - 1j], [1 - 1j, 3 + 2j]]))  # A == A^T
mtl5.ldlt_solve(S, b)  # LDL^T

H = mtl5.matrix(np.array([[2 + 0j, 1 - 1j], [1 + 1j, 3 + 0j]]))  # A == A^H
mtl5.ldlt_solve(H, b)  # LDL^H

mtl5.cholesky follows the same principle. For complex it computes A = L·Lᴴ, since the plain L·Lᵀ form and its diagonal ordering test mean nothing for complex elements — MTL5 static_asserts against that and this routes to its Hermitian variant. A non-real diagonal is reported as not Hermitian rather than as a definiteness failure, since those are different mistakes.

Complex is not in mtl5.dtypes(), which lists what mtl5.convert() accepts — the Universal number systems are real-only, so there is no complex target.

N-dimensional arrays

mtl5.array exposes MTL5's mtl/array layer — ranks 1 through 4, float32 and float64:

import mtl5

x = mtl5.array.asarray(a)  # zero-copy view of any strided NumPy array
x[1, 2]  # element access
x[:, 1:3]  # a view, never a copy
x.T.sum_axis(0)  # transpose is a view too
x.to_numpy()  # back to NumPy, strides included

Rank is a C++ template parameter, so the ranks are fixed at build time. Rank 5+ raises rather than silently flattening.

The point of the layer here is that it makes the dense containers sliceable without a copy:

M = mtl5.matrix(a)
col = mtl5.array.as_ndarray(M)[:, 2]  # a column of a DenseMatrix, no copy

asarray accepts any strided layout, so a NumPy transpose or slice comes through without materialising. Three things raise rather than quietly doing something else, all because the result aliases your buffer:

• negative strides — MTL5's strides are unsigned
• any dtype but float32/float64 — converting would return a view of a temporary, which is neither zero-copy nor the dtype you asked for
• read-only arrays — the view can write through. The error names the copying alternative; everything that only reads its input (norm, solve, convert, the *_copy factories, ...) takes a read-only array, such as a pandas 3 to_numpy(), as it is

Two details worth knowing:

.strides is in elements, not bytes. That is how MTL5 holds them. NumPy's .strides is in bytes, so the two differ by itemsize — [3, 1] here is (24, 8) there for float64.

Broadcasting is same-rank only. An extent of 1 stretches, so (2,3) + (2,1) works, but NumPy's rank promotion — (2,3) + (3,) — does not, because MTL5's broadcast_shape takes two shapes of equal rank. Reshape the operand first.

reshape and ravel return a view when the layout allows it and a copy otherwise, and never error; flatten always copies. All three give NumPy's element order even for a transposed or sliced source. MTL5's own reshape throws rather than copying in that case, which is right for a C++ caller but is not NumPy's contract — see docs/gap-analysis-2026-08.md §3.9.

Accumulator policy on the sparse factorizations

The four sparse direct factorizations that use a dense numeric workspace — splu, klu, supernodal_lu, supernodal_ldlt — take an accumulator= argument that types that workspace. A float32 factor can accumulate its updates in float64:

f = ms.splu(A32, accumulator="f64")  # narrow factor, wide arithmetic
f.accumulator  # 'f64'

The factor itself stays in the element type. Only the arithmetic that produced it widens: the accumulator removes the intermediate roundings inside each column's update chain, and each L/U entry is still rounded once on the way out. This is the mixed-precision knob a fixed-precision library structurally cannot offer.

On a badly scaled random sparse matrix (n = 400, cond ≈ 2.8e10) factored in float32, forward error of the direct solve — measured on x86-64 Linux:

ordering	nnz(L+U)	f32	f64	gain
colamd	81 712	8.19 × 10⁻³	6.39 × 10⁻³	1.28×
amd	92 026	5.35 × 10⁻²	3.17 × 10⁻²	1.69×
natural	98 729	3.52 × 10⁻²	1.39 × 10⁻²	2.54×

The direction is robust — a wider accumulator gives a smaller forward error in every ordering, matrix and platform tried, and the test suite asserts that. The magnitude is not: it ranges from roughly 1.3× to 3.3×, and which ordering benefits most changes between platforms, so read the column above as one measurement rather than a trend. (An earlier draft of this section claimed the gain grows with fill; that held on Linux and reversed on macOS.)

A 1.3–3× improvement is worth having but is not the headline. The headline is what it does to iterative refinement, which is what you would actually pair a narrow factor with — same matrix, natural ordering, refined against a float64 residual:

accumulator	iterations	forward error
f32	6	2.0 × 10⁻⁹
f64	4	1.4 × 10⁻¹⁰

A third fewer iterations and 14× closer, from a factor occupying the same memory.

Valid accumulators are "f32", "f64", "fma32", "fma64", or None for the element type. fma64 is measurably identical to f64 here. Two are refused rather than accepted quietly:

• narrower than the element type — f32 on a float64 matrix loses precision rather than gaining it
• quire — Universal defines quire<float> and quire<double>, but the factorizations keep one accumulator per row in their workspace, and a 4226-bit quire there needs its own design

One asymmetry worth knowing: refactor() replays the stored pivot sequence through an upstream entry point that takes no accumulator, so on splu, klu and supernodal_lu it refuses when a non-default accumulator is in effect rather than silently reverting to element precision — construct a new factorization instead. supernodal_ldlt.refactor() re-runs the full numeric factorization, so it carries the policy through.

Sparse storage formats

CSR is the default and what every solver here takes. Two other formats are available when the shape of the work suits them:

import mtl5.sparse as ms

C = ms.coo_matrix(row, col, data, shape=(n, n))  # triplets
C.insert(3, 7, 1.5)  # append; no pattern needed up front
A = C.tocsr()  # duplicates sum here

E = ms.ell_matrix(A)  # fixed-width, vectorisable
E.padding_ratio  # ...if this is small

COO is the format you build in when you don't know the pattern yet. It accumulates: a repeated (row, col) sums rather than overwrites, exactly as scipy.sparse.coo_matrix does. The duplicates stay separate — nnz counts them individually and a round trip through to_scipy preserves them — until tocsr() folds them into one entry.

ELL stores nrows × max_width slots and pads the short rows, which is what makes it regular enough to vectorise over. The cost is that one long row forces every other row to carry empty slots, so it is only a good idea when the row widths are near-uniform. padding_ratio is that decision as a number — 0.0 when every slot is occupied, and rising toward 1.0 as the matrix gets more ragged:

# a 2-D Laplacian: every row has 5 entries bar the boundary
ms.ell_matrix(mtl5.generators.laplacian_2d(32, 32)).padding_ratio  # 0.025

# an arrow matrix: one dense row and column, the rest diagonal
ms.ell_matrix(ms.from_scipy(arrow)).padding_ratio  # 0.985 — use CSR

ELL has no incremental build path: MTL5's ell_matrix has no element setter, so CSR is the only way in.

There is no CSC. MTL5's compressed2D used to accept a col_major orientation tag and ignore it — the inserter, element access and mult all treated the storage as row-major regardless, so genuine CSC arrays came back transposed with no error. That is now a compile-time rejection upstream (mtl5#355), but there is still no column-major container, so from_scipy converts a csc_matrix to CSR at the boundary — which is honest about the cost.

The thing CSC is usually wanted for here is a transpose product, and that is available directly, using MTL5's own transpose view rather than a scipy round trip:

A.rmatvec(x)  # A.T @ x, no second copy of the matrix

mtl5.sparse.as_linear_operator(A) uses it, so both directions of a SciPy LinearOperator now stay inside MTL5.

Tensor algebra

mtl5.tensor is index-notation tensor algebra — distinct from mtl5.array, which is NumPy-shaped N-D data. This one has fixed dimensions, Einstein summation and a metric:

import numpy as np, mtl5.tensor as mtt

A = mtt.asarray(np.arange(9.0).reshape(3, 3))
x = mtt.asarray(np.array([1.0, 2.0, 3.0]))

mtt.contract(A, "ij", x, "j")  # A @ x
mtt.contract(A, "ji", x, "j")  # A.T @ x
mtt.outer(x, x)  # rank 2

Rank and dimension are both compile-time in MTL5, so each pair is a separate type. Ranks 1, 2 and 4 over dimensions 2, 3 and 4 cover the module: nothing upstream uses rank 3, and rank 4 arises only as outer(rank2, rank2). Anything else raises rather than silently reshaping — use mtl5.array for general N-D data.

Index strings are read at runtime, which MTL5's own contract cannot do — it takes index names as compile-time template parameters. The space is small enough to enumerate, though: the repeated index sits in one of two positions on each side, so four instantiations cover every rank-2 contraction up to relabelling. Python therefore reaches MTL5's real contraction rather than a reimplementation. The letters are arbitrary; only which positions they share matters, so contract(A, "pq", B, "qr") and contract(A, "ij", B, "jk") agree.

A trace ("ii"), two shared indices, or none at all are each refused with a message naming which it was.

The metric operations cover raising and lowering. raise is a Python keyword, so the rank-1 pair is raise_index / lower_index, alongside upstream's raise_first / lower_second for rank 2:

g = mtt.minkowski_metric()  # signature (-, +, +, +)
v = mtt.asarray(np.array([1.0, 2.0, 3.0, 4.0]))
mtt.lower_index(v, g)  # [-1, 2, 3, 4]
mtt.raise_index(mtt.lower_index(v, g), g)  # back to v

is_symmetric and is_antisymmetric take a tolerance, and SymmetricTensor_d{dim}_{dtype} stores D(D+1)/2 components instead of D².

Matrix views

Eight ways to read part of a matrix:

mtl5.view.lower(A)  # and upper, strict_lower, strict_upper
mtl5.view.transposed(A)
mtl5.view.banded(A, 1, 1)  # tridiagonal
mtl5.view.map(A, [2, 0], [1, 2])  # out[i, j] = A[rows[i], cols[j]]
mtl5.view.hermitian(A)

These materialise — each returns an ordinary DenseMatrix that every other binding accepts. Upstream they are lazy accessors holding a const Matrix&, which buys nothing here (Python has no expression template to feed) and is a lifetime hazard: a caller passing a NumPy array or a converted scipy matrix hands over a temporary the view would outlive.

Two carry a trap worth knowing.

banded takes bandwidths, not signed diagonal offsets. banded(A, 1, 1) is tridiagonal and banded(A, 0, 0) the diagonal. Writing -1 for "one subdiagonal" reads plausibly but asks upstream for the wrong region entirely — negative values are refused rather than silently returning a band you did not mean.

hermitian is not the adjoint. It reads the upper triangle and mirrors it conjugated into the lower, discarding whatever was stored there — the usual "this matrix is Hermitian, only one triangle is meaningful" convention. For Aᴴ use mtl5.adjoint(A). The two agree only when the input is already Hermitian, which is exactly the case you would reach for first when checking:

Z = np.array([[1 + 0j, 2 - 1j], [9 + 9j, 4 + 0j]])  # lower triangle inconsistent on purpose
mtl5.view.hermitian(mtl5.matrix(Z))  # [[1, 2-1j], [2+1j, 4]] — lower discarded
mtl5.adjoint(mtl5.matrix(Z))  # [[1, 9-9j], [2+1j, 4]] — a real transpose

A non-real diagonal is refused: the view leaves the diagonal alone, so it would produce a matrix that is not in fact Hermitian.

Multigrid and smoothers

import mtl5, mtl5.sparse as ms

M = mtl5.mg.multigrid_1d(A, n_levels=5)  # builds the whole hierarchy
x = M.vcycle(mtl5.vector(x0), mtl5.vector(b), cycles=10)
M.level_sizes  # [127, 63, 31, 15, 7]

On 1-D Poisson the residual drops by a factor of about 0.03 per V-cycle, and that factor holds as the problem grows — mesh independence is what separates multigrid from a lone smoother, and the test suite asserts it rather than just checking the residual went down.

multigrid_1d builds the hierarchy for you: standard 1-D coarsening, Galerkin coarse operators, the chosen smoother at every level. n_levels is an upper bound — coarsening stops before a level would fall below 4 rows, so check level_sizes. wcycle is also available.

The seven smoothers are reachable directly, for use outside a hierarchy:

x = mtl5.mg.smooth(A, x, b, kind="symmetric_gauss_seidel", sweeps=5)
mtl5.mg.smoothers()
# ['jacobi', 'gauss_seidel', 'backward_gauss_seidel',
#  'symmetric_gauss_seidel', 'sor', 'backward_sor', 'symmetric_sor']

omega applies to the SOR variants; at omega=1.0 SOR is exactly Gauss-Seidel, which the tests use to confirm the parameter reaches the kernel.

Jacobi is a poor multigrid smoother, and that is a property of the method rather than of this binding. Multigrid needs a smoother that damps high-frequency error — the part the coarse grid cannot represent — and undamped Jacobi leaves the highest modes almost untouched. The usual remedy is damped Jacobi at omega ≈ 2/3, which MTL5's jacobi has no parameter for. Use a Gauss-Seidel or SOR variant in a hierarchy.

Grid transfer is exposed too, along with a sparse Galerkin product:

R = mtl5.mg.make_restriction_1d(127)  # 63 x 127, full weighting
P = mtl5.mg.make_prolongation_1d(63)  # 127 x 63, linear interpolation
Ac = mtl5.mg.galerkin(R, A, P)  # the coarse operator, still sparse

galerkin exists because spelling it R @ A @ P would go through a sparse-times-sparse product that returns a dense matrix — a fine-sized intermediate, which is exactly what you cannot afford in a hierarchy.

Krylov solvers and preconditioners

Ten solvers and eight preconditioners, in any combination except one, which is refused (see below):

import mtl5.sparse as ms

x, info = ms.gmres(A, b, M=ms.ilu0(A))
x, info = ms.cg(A, b, M=ms.ic0(A), rtol=1e-12)
x, info = ms.iterative_solve(A, b, solver="idr_s", M=ms.ssor(A, omega=1.4), s=8)

ms.solvers() and ms.preconditioners() list them:

symmetric	cg, minres
general	bicgstab, bicgstab_ell, cgs, gmres, idr_s, tfqmr, bicg, qmr
preconditioners	identity, diagonal, ic0, ildl, ilu0, ssor, ilut, block_diagonal

M= defaults to identity, so omitting it runs unpreconditioned. Solver-specific knobs are restart= (gmres), ell= (bicgstab_ell) and s= (idr_s).

The preconditioner is type-erased rather than being a template parameter of each solver, so the pairing is chosen at runtime. Binding the cross product directly would have meant 10 × 8 × 2 dtypes = 160 instantiations of a full iterative solver; this is 20.

bicg and qmr require a preconditioner with an exact adjoint. They are the only two that apply M⁻ᴴ. Every preconditioner here computes a real adjoint, except ssor: MTL5 implements its adjoint as its forward solve, which is exact only when A is symmetric. Each preconditioner reports this as .has_exact_adjoint, and .adjoint_solve(b) applies the adjoint directly. A pairing that would break down is refused rather than run:

ms.bicg(nonsymmetric_A, b, M=ms.ssor(nonsymmetric_A))  # ValueError, by design
ms.bicg(nonsymmetric_A, b, M=ms.ilu0(nonsymmetric_A))  # fine

One behaviour worth knowing: TFQMR can stagnate unpreconditioned — it plateaus rather than converging slowly, and more iterations do not help. Giving it any real preconditioner fixes it. That is a property of the method, not of this binding.

Sparse direct solvers

Seven factorizations, one interface — construct, .solve(b), .refactor(A2):

	for	notes
ms.splu	general square	Gilbert–Peierls, threshold pivoting
ms.klu	circuit matrices	block triangular form + per-block LU
ms.supernodal_lu	general square	dense block updates; .nsuper
ms.cholesky	symmetric positive definite	cheapest when it applies
ms.ldlt	symmetric, possibly indefinite	.diagonal() gives the inertia
ms.supernodal_ldlt	symmetric	dense block updates
ms.qr	least squares, rectangular	min ‖Ax − b‖₂

import mtl5.sparse as ms

lu = ms.splu(A, ordering="amd")  # analyze (ordering + symbolic) then factor
x = lu.solve(b)

lu.refactor(A2)  # numeric only — same pattern, new values
x2 = lu.solve(b2)

k = ms.klu(A)  # block triangular form + per-block LU
k.nblocks  # how reducible the matrix turned out to be

Two things scipy.sparse.linalg.splu cannot do.

Refactorization. A sequence of matrices sharing one sparsity pattern — the circuit-transient case — pays for the ordering and symbolic analysis once. On a 2-D Laplacian, n=3600, nnz(A)=17,760:

	nnz(factor)	factor	refactor	speedup
splu	205,636	12.5 ms	4.5 ms	2.8×
klu	119,530	6.7 ms	2.2 ms	3.0×
supernodal_lu	205,636	19.6 ms	4.0 ms	5.0×
cholesky	59,765	22.3 ms	20.8 ms	1.1×
ldlt	56,165	8.4 ms	7.0 ms	1.2×
supernodal_ldlt	56,165	4.7 ms	2.2 ms	2.1×

Two things to read off that table. Exploiting symmetry cuts the fill to about a quarter, and supernodal_ldlt is both the sparsest and the fastest option for a symmetric matrix. And the refactor win is small for cholesky/ldlt — they do not pivot, so their analysis is a cheap symbolic pass and there is little to skip; the LU-family factorizations, which must otherwise redo ordering and the pivot search, gain the most.

The saving is the analysis, so it scales with how much of the runtime that is. On a random pattern with catastrophic fill, numeric work dominates and refactor is no faster (we measured a slight loss). Structured sparsity is where it pays.

A factor narrower than the residual. The factor's precision is chosen by the dtype of A, independent of the precision you refine in:

lu32 = ms.splu(A.astype(np.float32), ordering="amd")
x, info = mtl5.mixed.iterative_refine(ms.from_scipy(A), lu32, b, rel_tol=1e-14)

On a 2-D Laplacian, n=1600:

	forward error
float64 factor, direct solve	2.0 × 10⁻¹⁵
float32 factor, direct solve	8.9 × 10⁻⁷
float32 factor + float64 refinement (3 iters)	7.0 × 10⁻¹⁶

Half the factorization memory and traffic, and the refined answer is better than the float64 direct solve. Every square factorization above works this way; qr is excluded because least squares is not the square system refinement corrects.

Orderings

ms.orderings() lists amd, colamd, rcm and natural. Each is also available standalone as a permutation, for inspection or for use on a scipy matrix directly:

p = ms.amd(A)  # or ms.colamd(A), ms.rcm(A), ms.ordering(A, name)
A[p][:, p]

The choice matters. Nonzeros in L+U for the 2-D Laplacian above (n=1600, 7840 nonzeros in A):

ordering	amd	colamd	rcm	natural
nnz(L+U)	41,542	62,944	90,040	128,078

colamd is the default because it suits unsymmetric matrices; amd is the better choice when the pattern is symmetric, as here.

Performance

Threading

MTL5's kernels are parallel but serial by default. The thread pool is sized once, on the first MTL5 call, and cannot be resized afterwards — so configure it before you do any work:

import mtl5

mtl5.set_num_threads(8)  # must precede the first MTL5 operation
print(mtl5.get_num_threads())

Equivalently, set MTL5_NUM_THREADS in the environment before importing:

MTL5_NUM_THREADS=8 python my_script.py

set_num_threads() raises RuntimeError rather than silently doing nothing if the pool has already been sized. The count is clamped to hardware concurrency.

All non-trivial kernels release the GIL, so MTL5 work overlaps with other Python threads instead of blocking the interpreter.

Double-precision matmul, 1000×1000, on a 20-core box:

threads	1	2	4	8
GFLOP/s	14.4	27.2	48.0	72.7

Build options

The wheel builds MTL5's blocked GEMM vectorised through Google Highway by default. Everything else is opt-in, since it needs a library on the build machine or makes the binary non-portable:

pip install . -C cmake.define.MTL5_WITH_BLAS=ON      # external BLAS
pip install . -C cmake.define.MTL5_WITH_LAPACK=ON    # geev/syev/gesdd dispatch
pip install . -C cmake.define.MTL5_NATIVE_ARCH=ON    # -march=native, local builds only

MTL5_NATIVE_ARCH=ON lets Highway target AVX2/AVX-512 rather than the x86-64 baseline — a large additional speedup, but the resulting binary only runs on machines like the one that built it.

Check what a given install actually has:

>>> mtl5.build_info()
{'blas': False, 'lapack': False, 'native_fast_gemm': True,
 'highway_simd': True, 'kpu': False, 'zlib': False, 'build_isa': 'SSE2'}
>>> mtl5.get_backend()
'native'

set_backend() validates against this build rather than silently accepting a backend that was never compiled in; backend selection itself is compile-time.

build_isa is the SIMD ISA this binary was compiled for, read from the compiler's own predefined macros. It is the field that makes the rest interpretable: highway_simd: True says Highway vectorised the micro-kernel, not how wide it went. Highway uses static dispatch — one ISA, fixed by the compiler flags — so a released wheel is x86-64 baseline whatever the machine running it supports. mtl5.system_info() reports that machine:

>>> mtl5.system_info()["cpu_simd"]      # what the CPU can do
'SSE2 AVX AVX2 FMA'
>>> mtl5.build_info()["build_isa"]      # what this binary may use
'SSE2'

A gap between those two is normal for a distributed wheel and is the reason MTL5_NATIVE_ARCH=ON is worth it for a local build. A gap you did not expect — after setting -march=native — means the flag never reached the compiler, which costs several x and shows up nowhere else. system_info() also carries the CPU brand, core count, OS, compiler and Release/Debug, which is what a benchmark report needs to be reproducible.

Benchmarking the number systems

The emulated formats are much slower than hardware float, and how much is the number that decides whether a mixed-precision experiment fits in a night. benchmarks/bench_blas.py measures it — dot, gemv, gemm, lu and qr across the bound number systems — and reports slowdown against float64 rather than raw GFLOP/s, because that is the figure experiment planning is built from.

python benchmarks/bench_blas.py            # default sweep
python benchmarks/bench_blas.py --quick    # smallest sizes, seconds
python benchmarks/bench_blas.py --json out.json

For scale, one recorded run (x86_64, single thread, no -march=native), GEMM at n=64: takum32 434x float64, cfloat32 570x, posit32 4024x. Those are an example from one host, not portable constants — every report embeds build_info(), thread count and architecture for that reason.

cfloat32 is the one to compare against rather than float: it is IEEE binary32 through the same emulation layer, so posit32 / cfloat32 separates what a format costs from what emulation costs.

benchmarks/README.md documents the methodology and how to read the output — including two traps: the small-n rows understate every ratio, and a narrow format can time perfectly normally while computing nothing useful.

Development

pip install -e ".[dev]"
pytest -v