Vectors & Matrices (linalg)

The linalg package provides a numpy-style API for vectors and matrices: constructors, shape manipulation, elementwise arithmetic, products, reductions, discrete calculus (differences and numerical derivatives), and linear algebra (determinant, inverse, linear solving).

Like the JSON library, linalg is implemented once in Lisp source (linalg.lisp): the interpreter loads the definitions lazily on the first use of a linalg: function, and the compile path splices them into the program when it references the package. There is no per-backend code, so every function behaves identically on the interpreter, the JVM compiler, WASM Preview 1 and the WASM component.

Data representation

linalg constructors build packed float arrays: unboxed (array double-float) values, the same representation as an #d(...) literal. A vector is a rank-1 array, printed #d(1.0 2.0 ...), and a matrix is a rank-2 array, printed with the nested #d((...) ...) form -- the #d marks the unboxed packed representation, so the printed form reads back as a packed array. Individual elements are read and written with aref, and any array built elsewhere -- packed or a general boxed array -- can be handed to a linalg function. Arrays of higher rank work too: the elementwise operations, the reductions, reshape/flatten and array-equal walk the elements in flat row-major order and accept any rank, while dot/matmul/outer/det/inv/solve/trace/transpose stay defined for vectors and matrices (rank <= 2), like numpy's specialized routines. linalg:from-list / linalg:to-list convert between arrays and lists.

linalg computes in floating point, prioritizing speed: every constructor and array-building operation returns a packed double-float array by default (single-float is available too -- see Single-float precision), and linalg:det, linalg:inv and linalg:solve run in floating point (like numpy), so a general inverse carries the usual rounding and a nearly singular determinant can be a small epsilon rather than exactly 0. A reduction follows the element type numpy-style: a reduction over a packed or float array is a double, while a reduction over a plain integer array (a bare #(1 2 3) literal) stays an integer or exact ratio; linalg:norm is always a float because sqrt is. One cross-backend caveat: the WASM backends print a non-terminating float at fewer significant digits than the interpreter and JVM, so a rounded inverse or an irrational norm can look different across backends even though the underlying double is identical.

A worked example

The inv and solve matrices above are chosen so their float results are exact and print identically on every backend; a general inverse such as (linalg:inv #2A((1 2) (3 4))) computes the same values but carries floating-point rounding.

Elementwise arithmetic and broadcasting

linalg:add, linalg:sub, linalg:mul and linalg:div operate elementwise and broadcast by numpy's rules: a scalar operand on either side is broadcast over the other operand's shape, and two arrays of different shapes align their trailing axes -- each aligned pair of extents must be equal or contain a 1 (a missing leading axis counts as 1), and the axis of extent 1 is stretched over the other operand's extent. A pair that fits neither rule signals a shape-mismatch error. The result keeps the first array operand's element type, matching the mixed-width rule. Note that mul is the Hadamard (elementwise) product -- the matrix product is linalg:matmul (or the rank-dispatching linalg:dot). Arbitrary per-element transformations go through linalg:emap.

The frequent per-element operations also exist under their numpy ufunc names: linalg:exp, linalg:log, linalg:tanh, linalg:sin, linalg:cos, linalg:tan, linalg:asin, linalg:acos, linalg:atan, linalg:sinh, linalg:cosh, linalg:sqrt, linalg:abs, linalg:square, linalg:negative, linalg:sign and linalg:reciprocal, plus the comparison selects linalg:maximum, linalg:minimum, linalg:clip and linalg:relu (defined by the strict comparison (if (> x y) x y) and its mirrors, so the second operand or the bound wins any false comparison -- ties and NaN included, identically on every backend). Each is equivalent to the obvious emap (or mul / div / maximum / minimum call), but as named functions they are accelerated under --simd, which emap with an arbitrary callback never is.

Reductions along an axis

The reductions linalg:sum, linalg:mean, linalg:amax and linalg:amin take an optional integer axis (negative counts from the end, numpy style) and reduce along that axis instead of over the whole array: the axis is dropped from the result, or kept as extent 1 when the optional keepdims flag is non-nil -- the shape that broadcasts back over the input, which is how a batch softmax subtracts its row maxima. linalg:argmax and linalg:argmin take the same axis argument and return per-slice indices (a packed double array for matrices, since linalg arrays have no integer width). linalg:reshape accepts one -1 extent and infers it from the element count.

Indexing, selection and masks

linalg:take-rows selects axis-0 slices by an index vector (numpy's x[mask], any rank), linalg:gather picks one element per row (y[np.arange(n), t]), and linalg:one-hot builds a label matrix. The elementwise comparisons linalg:equal, linalg:greater, linalg:greater-equal, linalg:less and linalg:less-equal return 0.0/1.0 masks (with scalar operands and broadcasting): multiply by a mask where numpy would boolean-index. linalg:zeros-like allocates a zero array of the same shape and width.

Random numbers

The np.random analog is seeded and cross-backend deterministic: linalg:seed resets a Wichmann-Hill generator whose draws are exact integer and IEEE double arithmetic, so a seeded sequence of linalg:rand, linalg:randn, linalg:uniform, linalg:choice and linalg:permutation is bit-identical on the interpreter, the JVM and both WASM targets -- weight initialization and mini-batch sampling reproduce exactly everywhere. randn uses the Irwin-Hall sum of twelve uniforms rather than Box-Muller (whose log/cos would diverge on WASM), so its tails clip at six standard deviations; fine for initialization, but not a distribution-exact np.random.randn.

Discrete calculus

linalg:diff and linalg:gradient are numpy's discrete-calculus pair (np.diff / np.gradient). diff takes the n-th discrete difference along the last axis (default 1): each step shortens that axis by one, and a matrix differences within each row. gradient estimates the derivative of a vector of samples with second-order central differences (first-order one-sided at the two ends), so the result keeps the input's length; the optional second argument is either a uniform sample spacing (a number, default 1) or a coordinate vector of the same length for non-uniformly spaced samples. Both preserve the input's width like every other linalg transform. The arithmetic is floating point as usual, but sample values that differentiate exactly -- polynomials read at integer coordinates, like every example below -- print identically on every backend.

The gradient of #(0 1 4 9 16) -- the parabola y = x^2 sampled at x = 0..4 -- recovers the true derivative 2x exactly at the interior points (central differences are exact for quadratics; the two ends are first-order estimates), and the coordinate-vector form stays exact even for the unevenly spaced samples on the last line. examples/ml/numerical-calculus.lisp works these ideas through a projectile-motion walkthrough.

Single-float precision

linalg computes in double-float by default, but it is width-polymorphic: it accepts and preserves packed single-float (#f) arrays, which use half the memory and twice the SIMD lane count. Every constructor takes an optional trailing element-type (the default is 'double-float; pass 'single-float for a #f result), and every transform -- add/sub/mul/div/emap, transpose/reshape, dot/matmul/outer, inv/solve -- preserves its input's width. A single-float value therefore stays single-float all the way through: a functional weight update (linalg:sub W grad) keeps W's width rather than silently widening it back to double (which, on the JVM --simd path, would force a mixed-width error on the following vec:matvec). Reach for single-float when you want the speed and memory of f32 and can accept its lower precision, and keep the double-float default for precision-critical work such as det/inv/solve.

SIMD acceleration

linalg needs no flag to be correct anywhere, but the --simd flag accelerates it: thirty-two functions -- add, sub, mul, div, sum, norm, amax, amin, argmax, argmin, trace, transpose, reshape, dot, outer, the unary ufuncs exp, log, tanh, sin, cos, tan, asin, acos, atan, sinh, cosh, sqrt, abs, negative, sign, and the comparison selects maximum, minimum -- are routed to native vector kernels (jdk.incubator.vector on the interpreter and the JVM, WebAssembly v128 on wasm-GC), and mean, matmul, flatten, solve, square, reciprocal, clip and relu are accelerated with them because they are written in terms of them. Nothing changes in what a program accepts or rejects: an input a kernel cannot handle (a general boxed array, mixed widths, a plain number) simply runs the portable linalg.lisp definition instead, with the same result and the same error messages. The only observable difference is the single-float reduction precision rule; element-wise results and the full matrix product stay bit-identical.

There is no reason to switch packages for speed: under --simd, linalg and vec land on the same kernels. See Choosing between vec and linalg -- the short version is: write against linalg by default, and reach for vec only for its -into destination-passing loops, a --no-gc target (where linalg does not compile), or its fail-fast width strictness.

First-class functions

linalg functions are ordinary defuns, so #'linalg:norm and friends work as first-class values wherever a function is expected:

Because arrays compare by identity (eq) only, results are compared with linalg:array-equal, which checks shape and numeric equality (1 and 1.0 compare equal).

Packages

linalg is a package of its own and does not use cl: inside (in-package linalg) the standard functions are not visible under their bare names and would need cl: qualification (cl:print, cl:mapcar, ...). Most programs should therefore stay in the default cl-user package and call the qualified linalg: names, as every example on this page does.