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 matmul stacks rank >= 3 on its last two axes (numpy's own np.matmul rule) and dot/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.

la is a built-in nickname for linalg, so every linalg: call can also be written with the shorter la: prefix:

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 four also answer to their CL operator spellings -- linalg:+, linalg:-, linalg:* and linalg:/ -- which are n-ary left folds of add / sub / mul / div, so (linalg:+ a b c) broadcasts step by step and each step is still the accelerated kernel. The degenerate arities follow CL: no argument gives the identity (0 / 1), a single argument to + and * is itself, and a single argument to - and / is the negation / reciprocal. (The vec: package's operator aliases are strictly binary instead -- its kernels are fixed-arity by design.)

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, the binary linalg:power, 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. linalg:softmax, linalg:log-softmax and linalg:erf sit here for the same reason relu does -- they are not in numpy proper (they are scipy's / torch's), but they are the array-level primitive an activation layer needs. The two softmaxes are the max-subtracted, numerically stable forms; erf is what the exact GELU is built from, and it sums the all-positive-term series rather than the alternating Maclaurin series, so it stays accurate to a double's last ulps where the naive form loses every digit by |x| ~ 3.

Reductions along an axis

The reductions linalg:sum, linalg:mean, linalg:amax and linalg:amin take numpy's keyword arguments: an integer :axis (negative counts from the end) reduces along that axis instead of over the whole array -- the axis is dropped from the result, or kept as extent 1 when :keepdims 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 keyword and return per-slice indices (a packed double array for matrices, since linalg arrays have no integer width). linalg:var and linalg:std take the same :axis / :keepdims plus a :ddof divisor correction (0 by default -- numpy's np.var and torch's unbiased=False; 1 gives the sample variance), which with mean along the same axis is the LayerNorm normalizer. linalg:reshape accepts one -1 extent and infers it from the element count.

Rank-N shapes: joins, slices and stacked products

Everything above is rank-generic, and a handful of operations exist to reshape that rank. linalg:expand-dims inserts an extent-1 axis and linalg:squeeze removes one (numpy's expand_dims / squeeze, torch's unsqueeze / squeeze); linalg:concatenate joins a list of arrays along an axis that already exists and linalg:stack along a new one; linalg:slice is basic numpy slicing, and linalg:triu / linalg:tril keep one triangle of a matrix.

rontolisp has no x[:, :n] syntax, so slice spells it as a list with one spec per axis: nil leaves that axis whole, (start end) or (start end step) selects along it, a negative index counts from the end, nil in the start or end position means "from the beginning" / "to the end", and a missing trailing spec leaves the remaining axes whole. Every axis is kept, exactly as numpy's x[:, 0:3] keeps both -- dropping an axis is what linalg:row does.

linalg:matmul is rank-generic too: at rank >= 3 on either side it is the stacked product (torch's bmm / matmul), where the last two axes are the matrix and every leading axis broadcasts. That is the shape a batched attention score has, so (batch heads n d) times (batch heads d n) gives (batch heads n n) in one call. (linalg:dot stays rank <= 2 on purpose: numpy's np.dot contracts against a different axis at higher rank, so passing it a stack signals an error pointing here rather than returning a wrong answer.)

Indexing, selection and masks

linalg:take-rows selects axis-0 slices by an index vector (numpy's x[mask], any rank) and keeps axis 0, while linalg:row takes one slice by an integer and drops it (numpy's x[i], so one image of a batch arrives at a forward pass as a plain vector). 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, or -- better -- pass it to linalg:where, which selects between two operands on a non-zero mask (numpy's np.where). Selecting rather than multiplying is what lets a -infinity mask reach linalg:softmax as a weight of exactly zero: multiplying an infinity by zero would give a NaN. 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 (default 1) along :axis (default the last axis): each step shortens that axis by one, so a matrix differences within each row by default and down each column with :axis 0. 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 :element-type keyword (the default is 'double-float; pass :element-type '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 precision rule, which covers the reductions and the matrix product; element-wise results stay bit-identical. That last sentence is --simd's alone: --gpu also takes the element-wise transcendentals, and a device computes them with its own library, so under that flag they are close to the portable definition rather than equal to it.

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.