v0.3.1

numax v0.3.1 moves the library to MAX 26.7 and Mojo 1.2.0.

• max-core is pinned to ==26.7 and Mojo to >=1.2.0,<2.0.0.
• CI uses prefix-dev/setup-pixi v0.11.2 with pixi v0.83.0.
• numax follows the MAX 26.7 kernel APIs: RawFileDescriptor, TileTensor subscripts in place of slice, std.sort taking a function, Coord-based matmul epilogues and matrix_band_part input functions, top_k's KEngine parameter, and the builtin_kernels homes of the reduce monoids and rowwise_softmax.
• Mojo 1.2's where prover no longer discharges a few equalities between a function's own constraints and the callee's, so those constraints are restated in the order the callee expects or the operands are retyped. The conditions themselves are unchanged.

There are no public API changes.

Full Changelog: v0.3.0...v0.3.1

numax v0.3.0

0.3.0 makes every Tensor-tier routine run on the accelerator, adds four SciPy areas numax never had, and renames a lot of the surface so one word means one thing. It is a hard break: the renames below land without deprecation aliases. 173 commits since v0.2.0, 1,824 tests in 157 suites, and the whole gate (tests, examples, accuracy, type-check, compile-fail checks, format, CI) ran on Linux with an NVIDIA A10G as well as on the Apple M3 Pro the release was developed on.

Breaking changes

If you are upgrading from 0.2.0, these are the edits you will make.

Old	New
View, ViewOver	TensorView, TensorViewOver (now exported)
Tensor.view()	Tensor.tile()
static_view, dynamic	as_static, as_dynamic
numax.core.array	numax.core.tensor
numax.core.tensor (the map/reduce drivers)	numax.core.functional
TensorLU, TensorQR, TensorSVD, TensorEigh, TensorSchur, TensorMinimizeResult, TensorIVPResult and the other Tensor-prefixed result structs	LU, QR, SVD, Eigh, Schur, MinimizeResult, IVPResult, ...
scalar-tier IVPResult, register-tier MinimizeResult	ScalarIVPResult, ArrayMinimizeResult
Chebyshev, CubicSpline in the register tier	ArrayChebyshev, ArrayCubicSpline
norm selectors inf, fro, neg_inf	NORM_INF, NORM_FRO, NORM_NEG_INF
signal modes full, same, valid	MODE_FULL, MODE_SAME, MODE_VALID
count-based arange[num](start, step)	arange_n; arange now takes NumPy's stop or start, stop, step

Behavior changes that are not renames:

• Tensor's constructors take the DeviceContext last and optional, like the factories. A wrong-length initializer raises "Tensor: N values for M elements" instead of reading past the end.
• Tensor.buffer and layout, and TensorView's tile and device, are private. Reach the data through to_host and copy_from_host.
• solve_triangular's upper now names the triangle that is stored.
• The six register-tier subpackages (linalg, fft, signal, integrate, optimize, interpolate) no longer have a public .array module. See "One name per operation".
• sort, ravel and diagflat on a run-time-shaped tensor no longer ask MAX for a buffer of -1 elements; a run-time-shaped input gives a run-time-shaped result.

Everything runs on the accelerator

Every public Tensor-tier routine now has a device path. Where MAX ships a kernel numax calls it. Where MAX does not, numax writes one in MAX's idiom and records the gap in docs/parity.md. Host-only is no longer something a routine may ship with.

What changed on the device:

• The spectral chain (sytrd, eigvalsh, eigh, hessenberg, eigvals, schur, gebrd, svdvals, svd, matrix_rank, and everything built on them: the matrix functions, pinv, orth, null_space, polar, cond, lstsq(method="svd")) works at gpu=True. In 0.2 that was a compile error, because sytrd[gpu=True] disagreed with the host at n = 4, block = 1. That is fixed and each routine has a device test.
• The operators (a + b and the rest) follow the tensor to wherever it lives. So do the shape routines (reshape, concatenate, stack, split, tile, repeat, roll, flip, broadcast_to, diag, rot90, and so on) and the creation functions (zeros, ones, full, eye, linspace, arange).
• Scans: cumsum, cumprod, cumulative_trapezoid. MAX's nn.cumsum is CPU-only, so this is numax's own blocked scan.
• Sorting and selection: sort, argsort, median, quantile, percentile, iqr, searchsorted, unique, rankdata. The sort is numax's own bitonic sort, because MAX's GPU argsort returns wrong results past 256 elements on Metal. Axis-wise median and mode sort every slice at once.
• Statistics: histograms (histogram, histogram2d, histogramdd, bincount, digitize) by atomics into per-block bins, the descriptive statistics, pearsonr, linregress, spearmanr, kendalltau, and every hypothesis test except ks_1samp, whose Float64 cdf parameter Metal cannot compile.
• Data-dependent output length: nonzero, argwhere, extract, compress, find_peaks, peak_prominences, by mask, scan and scatter. The one scalar read back sizes the result.
• Recurrences: lfilter, filtfilt, sosfilt, decimate, resample, hilbert, istft, savgol_filter, and spline construction.
• Banded solvers: solve_banded and solveh_banded for tridiagonal systems, by parallel cyclic reduction.
• Iterative drivers: least_squares and curve_fit keep residuals and the Jacobian on the device, and solve_ivp no longer stages state through the host.
• The Array tier's tier-2 routines (the optimize drivers, pinv, firwin) no longer need float64 or raises in their bodies, so they run inside a kernel.

The exceptions are listed per routine in docs/parity.md: a single banded Cholesky or Toeplitz system, the narrow-band iteration inside the eigen and SVD routines (its O(n^2) accumulation is on the device), and ks_1samp.

How it is checked: tests/ still never names gpu=True, because a GPU-less CI runner cannot compile the device path. A separate tests_gpu/ tree compares every gpu=True routine against its host answer at float32. It passed on Metal and, in this release, on an A10G (all 16 areas). CI compiles it for sm_80 on every push.

Four new SciPy areas

numax.spatial has distance.cdist, pdist and squareform; Euclidean goes through the Gram-matrix identity on linalg.matmul, the other metrics are one elementwise pass each. KDTree with query and query_ball_point, built as a linear BVH over Morton-ordered points by a device sort and queried one point per thread. Delaunay, ConvexHull, Voronoi and Rotation are not here.

numax.ndimage has convolve, correlate, gaussian_filter, uniform_filter, median_filter, map_coordinates, shift, zoom, label, binary_erosion, binary_dilation and distance_transform_edt. It is imported as numax.ndimage, not from the root, because its convolve and correlate are numax.signal's names with a different meaning.

numax.cluster has vq.vq, kmeans and kmeans2 as Lloyd's iteration over cdist and MAX's ArgMin, and hierarchy.linkage and fcluster. linkage runs the primitive algorithm with the distance matrix and each merge's row update on the device; it is not the nearest-neighbor-chain version.

numax.optimize gains linprog (SciPy's former homogeneous self-dual interior point, with the normal matrix formed and factored on the device), milp (branch and bound over it), and minimize(method="slsqp") with LinearConstraint and NonlinearConstraint. trust-constr and the global optimizers are not here.

Examples are clustering.mojo (spatial and cluster), image_processing.mojo and constrained_optimization.mojo under examples/intermediate/.

NumPy and SciPy gaps closed

• Core: unique_counts and unique_inverse (NumPy 2's split of unique), isin, intersect1d, union1d, setdiff1d, nan_to_num, logaddexp, logaddexp2, sinc, heaviside, gcd, lcm, the bitwise operations and shifts, nanargmax, nanargmin, ediff1d, unwrap, and keepdims on the axis reductions.
• fft: fftn, ifftn, rfftn, irfftn, irfft2, hfft, ihfft, dctn, idctn, dstn, idstn, and norm= on the transforms. The engine needed no axis argument: each pass transforms the last axis and writes through a rotated view. rfftn and irfftn stop at rank 3.
• signal: Second-order sections (output=OUTPUT_SOS on every design, plus tf2sos, zpk2sos, tf2zpk, sos2tf, sosfreqz, sosfilt_zi, sosfiltfilt), the order estimators buttord, cheb1ord, cheb2ord, ellipord, iirnotch, iirpeak, bilinear, group_delay, lfiltic and lfilter(zi), upfirdn, resample_poly, five more windows (tukey, gaussian, flattop, nuttall, chebwin), chirp, sawtooth, square, peak_widths, argrelextrema, convolve2d and correlate2d.
• stats: .rvs() on every distribution by inverse CDF over a shared Generator; Generator.choice, permutation, shuffle, gamma, beta, poisson, binomial, lognormal, multivariate_normal and spawn; mean, var, std, interval and entropy on every distribution; new distributions uniform_dist, lognorm, weibull_min, cauchy, laplace, rayleigh, logistic, pareto, bernoulli, geom, nbinom, hypergeom and multivariate_normal; tests kruskal, levene, bartlett, fisher_exact, binomtest, chi2_contingency, skewtest, kurtosistest, normaltest, jarque_bera, bootstrap and permutation_test; and median_abs_deviation, ecdf and gaussian_kde.
• special: expit, log_expit, ndtr, log_ndtr, ndtri, gammaincinv, betaincinv, polygamma, ellipkinc, ellipeinc, elliprf, elliprd and wofz. elliptic_k and elliptic_e are now computed through Carlson's forms. betainc's lower tail reaches below 1e-30.
• linalg: ldl, generalized eigh(a, b) and eigvalsh(a, b), solve_sylvester, solve_continuous_lyapunov, solve_discrete_lyapunov, solve_continuous_are, and lu returning P, L, U.
• integrate: solve_bvp, cumulative_simpson, romb, and t_eval, dense output and events on solve_ivp.
• interpolate: General-degree BSpline and make_interp_spline, interpn in up to eight dimensions, RBFInterpolator, NearestNDInterpolator and griddata (nearest-neighbor only, since there is no Delaunay).

Masked arrays, datetime64 and structured dtypes, the numpy.polynomial classes, and sparse and Krylov solvers are recorded as absent, each with its reason, in docs/parity.md.

Three accuracy fixes came with this: log1p, log2, exp2 and cosh at float64 now compute to double precision, max and min of all-infinite input return the infinity (MAX's `Red...

v0.2.0

The big things, and why you'd want this release

Dense linear algebra on Tensor is now most of scipy.linalg, not a third of it. 0.1 gave you cholesky, lu_factor, qr_factor and solve, blocked and running on MAX's GEMM. 0.2 adds the eigendecomposition, the SVD, the Schur form, and the matrix functions built on top of it, all with the same shape: device-resident, blocked, cubic term in linalg.matmul. If you were reaching for NumPy or SciPy to diagonalize a matrix, take an SVD, or compute expm/logm/sqrtm, that call now has a Mojo counterpart that runs on your GPU.

You can now borrow instead of copying. Every subpackage moved onto a TensorLike bound this release, so a routine written against it takes an owned Tensor or a borrowed View of a sub-block interchangeably. Factoring a 4x4 corner of a 1000x1000 matrix, or filtering a slice of a recording, no longer costs a copy first.

Thousands of small linear-algebra problems, one GPU thread each. map_blocks is new: hand it a batch of small problems in a structure-of-arrays layout and it runs one whole factorization, solve, or eigenproblem per SIMD lane or per GPU thread. examples/advanced/batched_solve.mojo runs 4096 independent 4x4 Cholesky solves this way, plus the derivative of every one of them at Dual, from one call.

scipy.signal's filter design is complete. butter, cheby1, cheby2, ellip, and the iirfilter front door over all four, plus freqz, are new. If you design filters, this closes the one part of numax.signal that wasn't here yet.

Two real bugs are fixed. maximum/minimum silently dropped NaN instead of propagating it (a NumPy-visible difference, now fixed, with fmax/fmin added as the NaN-skipping alternative). And the new spectral routines' gpu=True path, which answered wrong rather than raising, is now a compile-time error instead, so a silent wrong answer is no longer reachable.

Plus 38 new NumPy/SciPy names picked up in a full audit, and a documented reason for every name still missing (docs/parity.md).

Improvements and the detailed API surface

Spectral decompositions (new), over Tensor:

Function	What it gives you
eigh, eigvalsh	Symmetric eigendecomposition: ascending eigenvalues, orthonormal eigenvectors
svd, svdvals	Singular value decomposition, any m >= n
eigvals	General eigenvalues (re, im), Francis double-shift QR
schur	Real Schur form A = Z T Z^T, the base every matrix function past expm builds on
sytrd, gebrd, hessenberg	The tridiagonal, bidiagonal, and Hessenberg reductions underneath the four above, each its own TensorTridiagonal/TensorSVD-style result with a .q()

Every one takes a block size (block == n recovers the unblocked algorithm, pinned by tests) and a gpu: Bool, currently False-only: gpu=True is now a compile error on these six plus pinv/cond/matrix_rank, guarding the Metal correctness bug mentioned above rather than papering over it.

Matrix functions (new), over Tensor, on the Schur form: expm (scaling and squaring, Padé), logm (inverse scaling and squaring), sqrtm (Björck-Hammarling, handles repeated/defective eigenvalues), funm(a, f) (any FloatLike scalar function via block Parlett), cosm, sinm, fractional_matrix_power.

TensorLike and View (new): the trait bound behind owned/borrowed genericity, in numax/core/tensorlike.mojo. Every one of core, linalg, stats, signal, fft, integrate, interpolate, optimize, special and io was moved onto it this release. dim[T, i] reads a compile-time extent off either type; is_row_major[T] refuses a strided block at compile time where a routine needs to flatten.

map_blocks (new): numax/core/tensor.mojo. One small problem (a factorization, a solve, an eigenproblem) per lane instead of one scalar, structure-of-arrays so a lane's load is contiguous.

IIR filter design (new), over Tensor: butter, cheby1, cheby2, ellip (lowpass/highpass/bandpass/bandstop), iirfilter (dispatch by name), freqz. SciPy's exact route: analog prototype, frequency warping, lp2lp/lp2hp/lp2bp/lp2bs, bilinear transform.

The 38-name surface sweep, grouped by area, all built from routines already in the library:

• Elementwise: sign, square, reciprocal, rint, degrees, radians, fmax, fmin
• Array manipulation: flatten, dstack, atleast_3d, rot90, array_split
• Sorting/selection: compress, partition, argpartition
• Polynomials: polyval, polyder, polyint, roots, polyfit
• Linear algebra: orth, null_space, polar, tanm, rq
• Signal: csd, coherence, istft, decimate, oaconvolve, peak_prominences, zpk2tf
• Special/integration: betaln, fixed_quad, dblquad
• Statistics: ks_2samp, wilcoxon

Correctness fixes: maximum/minimum now propagate NaN on float dtypes (previously followed SIMD's maxNum/minNum, which silently skip it); gpu=True on any spectral routine or pinv/cond/matrix_rank is a compile error, not a wrong answer.

Faster, measured on the same machine before and after: quantile/percentile by three-way quickselect instead of a sort (1,096 ms to 151 ms at 2^24 elements); cov/corrcoef by Welford means plus one GEMM (125 ms to 4.1 ms); the FFT engine fused into one bit-reversal-plus-six-stages kernel (3x behind pocketfft on a 2^20 transform, down from 14x); lfilter/filtfilt reading through the tensor's own host mapping instead of copying into a List (filtfilt now ahead of SciPy).

Platform: bumped to max-core ==26.6 and mojo >=1.0.0,<2.0.0 (from 26.5/1.0).

Performance: NVIDIA A10G and Apple M3 Pro, and where numax wins

Every number below is measured on one processor at a time and never mixed with the other (docs/performance.md has the full harness for each).

Where numax wins outright:

• The GEMM every factorization runs on beats the vendor path on both platforms. On the A10G (Ampere), MAX's matmul reaches 20,459 GFLOP/s against cuBLAS's FP32 path at 15,342 (PyTorch) and 14,995 (CuPy), a 33-36% lead. On the CPU side (EPYC 7R32), it now reaches 93% of OpenBLAS, up from 79% before this release's MAX bump.
• The elementwise kernel matches hand-written CUDA and beats every ML framework's dispatch. On the A10G, a FloatLike Gaussian kernel with no device name and no dtype in its signature hits 61,013 M elem/s, within 1% of a hand-written cupy.ElementwiseKernel (60,352) and 10% ahead of torch.compile (53,670), at 83% of the card's 600 GB/s bandwidth spec. Same story on the M3 Pro's Metal GPU: 14,465 against torch.compile's 13,831, and about 3x ahead of eager PyTorch and MLX (4,807 and 4,866).
• solve and qr_factor overtake PyTorch's MPS backend as matrices get bigger. On the M3 Pro, solve passes PyTorch MPS at n = 2048 (59.7 against 50.4 GFLOP/s) and qr_factor is 1.48x ahead at n = 512.
• BLAS-1 beats the vendor CPU libraries outright. dot and nrm2 on Tensor beat OpenBLAS's own sdot/snrm2 by 2.5-5x on the EPYC and Accelerate's by 1.7-3.2x on the M3 Pro, because routing them through MAX's rowwise scaffolder threads them where the vendor routines don't.
• A run of NumPy/SciPy-shaped CPU wins on the M3 Pro (against Accelerate): cov/corrcoef 3.3x ahead, welch 8.9x, norm.cdf 6.3x, interp 9.4x, filtfilt now ahead of SciPy outright.

Where numax is still behind, stated rather than hidden:

• The new spectral decompositions ship CPU-only this release. eigh reaches 0.45 of LAPACK at n = 1024, svd 0.056; gpu=True is a compile error rather than a live path (see the correctness fix above). Closing the device path is the headline item for 0.3.
• The factorizations' device solvers still trail cuSOLVER's fused kernels, even though the GEMM underneath wins: on the A10G, cholesky is 62.5 GFLOP/s against cuSOLVER's 595.0, lu_factor 30.6 against 282.9. The gap is a single-block, effectively serial panel kernel and its launch count, not the multiply. cuSOLVER's advantage is a parallel panel; that's the next lever.
• quantile is 2.85x behind NumPy, because it still downloads the whole tensor before selecting on it.

What 0.2.0 represents

0.1 proved the two-axis idea (one kernel, several conformers; one tensor, any device) on a narrow slice of dense linear algebra. 0.2 is where the Tensor tier stops being three factorizations and a solve and becomes close to the whole of scipy.linalg, with genuine wins on the primitives underneath it (the GEMM, the elementwise dispatch, BLAS-1) and an honestly stated gap on the parts still catching up (the device spectral path, the factorization panels).

Requires max-core ==26.6 and mojo >=1.0.0,<2.0.0.

[dependencies]
numax = { git = "https://github.com/ehsanmok/numax.git", tag = "v0.2.0" }

Full diff: v0.1.0...v0.2.0

v0.1.0

First release of νMAX, a numerical computing library built on MAX.

One kernel, several meanings

Every kernel in the library is written once, against the FloatLike trait.
The conformer you instantiate it with decides what the call returns:

Type	What a call gives you
Plain[dtype, width]	Ordinary SIMD at hardware speed
Dual[Inner]	f(x) and f'(x) from one call, forward-mode
Gradient[Inner, n_vars]	Every ∂f/∂xᵢ from one call
Compensated[dtype, width]	Roughly double dtype's precision, same algorithm
Decimal	Exact base-10 fixed point
Complex[Inner]	Complex results
Interval[Inner]	Rigorous enclosures

They nest. Gradient[Dual[Plain]] is a Hessian, Complex[Dual[Plain]]
differentiates holomorphically, Dual[Compensated] takes a derivative in
double-double. Nobody implemented autodiff for the incomplete beta function;
the chain rule lives in Dual's arithmetic rather than in any kernel, so
differentiating one costs nothing extra.

from numax import Dual, FloatLike, Plain, f32

def g[T: FloatLike](x: T) -> T:              # written once, no `dtype`
    return (-(x * x)).exp()

def main():
    print(g(Plain[f32](0.5)).v)              # 0.7788008
    var d = g(Dual[Plain[f32]].seed(0.5))
    print(d.value.v, d.deriv.v)              # 0.7788008 -0.7788008

One tensor, every device

Tensor owns a MAX DeviceBuffer, so the DeviceContext handed to a
factory decides host or device memory. Nothing else in the program changes,
and .view() yields the TileTensor every MAX kernel takes. The tensor
drivers (map, reduce, reduce_axis, broadcast_op_rows,
map_threaded) take one gpu: Bool compile-time parameter instead of
shipping two implementations, and map[gpu=True] launches through
DeviceContext.enqueue_function on Metal and CUDA with no per-backend
branch anywhere in the library.

The two-tier rule is what pays for that. Tier-1 kernels run a fixed
iteration count with no per-lane branching, doing selection arithmetically
instead of with an if, which is what makes them GPU-launchable unmodified.
Tier 2 is Plain-only and may loop to a tolerance. Every module and every
function declares which tier it is in, so a reader at a call site never has
to audit a body.

What is in it

Subpackage	Contents
numax.core	FloatLike and its conformers, Tensor, elementwise math, comparisons, sorting, the tensor engine
numax.special	Γ and B, erf, Bessel J/Y, Lambert W, elliptic K/E, orthogonal polynomials, activations
numax.linalg	cholesky, lu, qr, eigh, svd, solve, inverse, pinv, det, cond, norms, matmul
numax.optimize	newton/halley/bisection at fixed iterations; newton_tol, brentq, bfgs to a tolerance
numax.integrate	Gauss-Legendre/Simpson/trapezoid, rk4/dopri5; quad, quad_vec, solve_ivp
numax.interpolate	Horner, cubic splines, Chebyshev fits
numax.fft	fft/ifft, rfft/irfft, fft2, fftfreq, fftshift
numax.signal	convolve, correlate, Hann/Hamming/Blackman windows
numax.stats	Descriptive statistics, nine scipy.stats-shaped distribution namespaces, sampling
numax.io	.npy read/write, byte-identical to numpy.save, plus numax's own NMX1 format

Subpackages carry NumPy and SciPy names, so an import has an obvious
counterpart. from numax.prelude import * is the common surface in one
line. It leaves out every name that would shadow a Mojo builtin, because a
module-level definition replaces the builtin rather than overloading it.

The .npy support is small but it removes a conversion step: a program
ported from NumPy can read the files it already has and hand results back
the same way.

Performance

Every number below was measured, and CPU and GPU results are never combined
into one comparison. On an A10G's GPU, numax reaches 60,955 M elem/s against
torch.compile's 53,607, about 83% of the card's bandwidth spec. On the same
machine's CPU, map over a Plain kernel runs at 0.998x a hand-written
raw-SIMD loop and 4.5x NumPy. Every row in docs/performance.md compares
implementations on one processor.

Accuracy

Every approximation documents an error bound, and pixi run accuracy checks
each one against mpmath references computed at 50 digits and checked into
the repository.

Limits worth knowing before you start

Shapes are compile-time parameters. There is no broadcasting, no strided
views, no slicing and no fancy indexing. Those wait on a runtime-shape
owner, because a view over a slice is not row-major and every driver in
numax.core.tensor requires that it is.

Above the tensor drivers the surface is not uniformly rank-generic yet.
numax.stats reduces over every element with no axis=; sorting flattens;
transpose is 2-D; concatenate, split and stack are rank-1. The
SciPy-shaped algorithms are fixed-size Array kernels, so linalg works on
matrices and fft2 on a square transform, while quad, solve_ivp, the
splines and the distributions are 1-D.

Anything whose length depends on a value rather than a parameter has to say
so in its signature. unique returns a full-length tensor plus a count
rather than a right-sized array, and its docstring tells you to read only
the first count entries.

.npz is out of scope. It is a zip container, and while stored-only entries
would be easy, deflate is not.

Running the GPU examples needs a real device. CI compiles them against a
pinned target so their call sites stay honest, but execution is a local
check, verified by hand on both Metal and CUDA.

This is a young library and the APIs may change.

Install

[workspace]
channels = ["https://conda.modular.com/max", "conda-forge"]
preview = ["pixi-build"]

[dependencies]
numax = { git = "https://github.com/ehsanmok/numax.git", tag = "v0.1.0" }

Requires pixi. mojo ==1.0.0 and max ==26.5 come in
transitively. numax is distributed as source rather than a .mojopkg,
because a package embeds the compiler version it was built with and refuses
to load against any other.

Verified before release on linux-64: the git dependency resolves from this
tag, builds, and runs the example above with no import path flag. The
manifest declares linux-64, linux-aarch64 and osx-arm64; CI runs on linux-64
and osx-arm64.

Verification

44 test suites and 664 tests, 14 CPU examples, the accuracy harness against
the mpmath references, the CPU benchmark, mojo doc over the whole public
generic surface, and a formatting check. All green on the tagged commit, on
both ubuntu-latest and macos-15.

Full documentation: features ·
architecture · parity ·
performance