Algorithms#

An operation's own tolerance covers every implementation that rounds only where the operation must. Real implementations sometimes fix a particular order of additions, accumulate in another format, or, like fused attention kernels, round in extra places.

A named algorithm is one concrete way of computing an operation, with its own bounds measured against the operation's rounded. Algorithms in the Guide explains when to use one.

matmul algorithms#

The constructors live in er.linalg.products.

matmul_sequential(format, *, fma: bool = True, rounding: Rounding = Rounding.NEAREST_EVEN) -> MatmulAlgorithm#: Products added one at a time in index order, in format.

matmul_pairwise(format, *, fma: bool = True, rounding: Rounding = Rounding.NEAREST_EVEN) -> MatmulAlgorithm#: Products summed as a balanced binary tree, in format.

matmul_blocked(format, block_size: int, *, fma: bool = True, rounding: Rounding = Rounding.NEAREST_EVEN) -> MatmulAlgorithm#: Products summed sequentially in blocks of block_size, then the block sums sequentially, in format.

matmul_exact_accumulation(format, *, rounding: Rounding = Rounding.NEAREST_EVEN) -> MatmulAlgorithm#: Products summed exactly (a long accumulator), then rounded once to format.

matmul_slots16(format, *, rounding: Rounding = Rounding.NEAREST_EVEN) -> MatmulSlots16#: Products rounded to format and summed in 16 magnitude slots (MatmulSlots16).

Each returns an algorithm object with these methods:

class MatmulAlgorithm#

One fully specified matmul: each output element's k products, exact (fused multiply-add) or rounded to the accumulator format, summed by one accumulation method in that format, rounding to nearest (rounding: ties to even or away), then rounded to the output format. Real floating-point inputs only.

run(A, B, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Array#: This algorithm's output, emulated exactly (on the rounded values of operands that are earlier results).

error_model(k: int) -> ErrorModel#: The bound on the accumulated sum of k products, before the output rounding.

covers(output, A, B, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, standard: str = 'normal', flush_to_zero: bool = False) -> Array#: Where this algorithm can produce output, as Result.covers describes; the other arguments are tolerance's.

finite_range(A, B, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, standard: str = 'normal', flush_to_zero: bool = False) -> tuple#: The lowest and highest finite output of this algorithm per element, as Result.finite_range describes; the arguments are tolerance's.

may_be(kind: str, A, B, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, standard: str = 'normal', flush_to_zero: bool = False) -> Array#: Where this algorithm may output kind ("+inf", "-inf" or "nan"), as Result.may_be describes; the other arguments are tolerance's.

tolerance(A, B, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, standard: str = 'normal', flush_to_zero: bool = False) -> Array#

Worst-case |output - matmul(A, B).rounded| per element over this algorithm's finite outputs, as Result.tolerance describes.

A and B may be earlier results; their tolerances are then taken with standard and flush_to_zero (which change nothing for arrays).

class MatmulSlots16#: The 16-slot matmul (proofs/matmul.md, The 16-slot algorithm): per output element, each product rounded to the accumulator format (to nearest, rounding: ties to even or away), the products split by magnitude into 16 slots of one sixteenth of the format's precision each (in binades below the largest product; anything smaller goes to the last slot), each slot adding its values in increasing magnitude, alternating positive and negative values, and the slot sums then added from the smallest slot up. Real floating-point inputs only.

Code

import exactref as er

# Products rounded to float32 (no fused multiply-add), added in blocks of 4
alg = er.linalg.products.matmul_blocked("float32", 4, fma=False)
# Its name, and what it does in words
print("name:", alg.name)
print("description:")
print(alg.description)

# A 1x8 row of 0.1 and an 8x1 column of 1.0, in float32
A = er.asarray([[0.1] * 8], "float32")
B = er.asarray([[1.0]] * 8, "float32")

# The algorithm's output, emulated exactly
print("run:", alg.run(A, B).tolist())
# matmul's exact answer, rounded once to float32
print("rounded:", er.matmul(A, B).rounded.tolist())

# The bound on the error of a sum of 8 products, before the output rounding
print("error model:", alg.error_model(8))

Output

name: matmul_blocked
description:
Blocked summation of products rounded to float32 in float32 in blocks of 4
run: [[0.800000011920929]]
rounded: [[0.800000011920929]]
error model: 4.17e-07 * sum over k of |A[i,k] * B[k,j]| + 1.26e-44

Attention algorithms#

class AttentionAlgorithm(*, online_softmax: bool = False, block_size: int | None = None, exp2: bool = False, weights=None, accumulate=None, bias=None, reciprocal: bool = False)#

A fused-kernel way of computing attention, described by its extra rounding points.

The reference is unchanged (scaled_dot_product_attention); an algorithm only has a tolerance(): how far an implementation that rounds at the required boundaries and at these points can be from the correct answer. It has no run(), because kernels differ in details these methods leave open. Methods, any combination:

online_softmax: keys processed in blocks with a running maximum and sum, rescaling earlier partial sums by exp(old max - new max). With block_size None, any blocking and order; otherwise at most ceil(S / block_size) blocks of at most block_size keys (as kernels split the keys: every block full but one), in any order.
exp2: weights computed as exp2 of logits in base-2 units (the scale times log2(e), and mask values times log2(e), held in the compute format within ulps).
weights: the weights rounded to this format before the value product (and in the backward, P and dS before their products).
accumulate: every accumulation in this format instead of the compute format.
bias: float-mask values rounded to this format before they are added (a kernel computing ALiBi biases in lower precision).
reciprocal: normalization by multiplying with the row sum's reciprocal, rounded once, instead of dividing: for the output, and for the backward's probabilities (PyTorch's CPU softmax and its CPU flash kernel do this). That kernel's logsumexp for bfloat16 and float16 inputs uses an exp less accurate than the default 2 ulps: ulps=8 covers it.

backward_tolerance given the forward's rounded out and logsumexp describes a fused backward kernel: it recomputes P = exp(logit - logsumexp) and takes D = sum grad * out. Without them it describes a backward that recomputes the softmax, as PyTorch's math path does. online_softmax applies to the forward only.

Every input of tolerance and backward_tolerance, out and logsumexp included, may also be an earlier result, as Result describes. The bounds are proven in proofs/attention_algorithms.md.

tolerance(query, key, value, attn_mask=None, dropout_p: float = 0.0, is_causal: bool = False, scale=None, enable_gqa: bool = False, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, compute=None, ulps=None, flush_to_zero: bool = False, standard: str = 'normal') -> Outputs#

Worst-case |output - reference| per element of output and logsumexp, over this algorithm's finite outputs (Result.tolerance).

The arguments are scaled_dot_product_attention's; compute, ulps, flush_to_zero and standard are as in Result.tolerance. With online_softmax both standards give the any-order bound.

finite_range(*args, **kwargs) -> Outputs#: Per output, the lowest and highest finite output of this algorithm, as Result.finite_range describes; the arguments are tolerance's.

may_be(kind: str, *args, **kwargs) -> Outputs#: Per output, where this algorithm may output kind ("+inf", "-inf" or "nan"), as Result.may_be describes; the other arguments are tolerance's.

covers(outputs: dict, *args, **kwargs) -> Outputs#: For each output named in outputs ({"output": ..., "logsumexp": ...}, either or both), where this algorithm can produce that value, as Result.covers describes; the other arguments are tolerance's.

backward_tolerance(query, key, value, grad, out=None, logsumexp=None, attn_mask=None, dropout_p: float = 0.0, is_causal: bool = False, scale=None, enable_gqa: bool = False, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, compute=None, ulps=None, flush_to_zero: bool = False, standard: str = 'normal') -> Outputs#: Worst-case |gradient - reference| per element of grad_query, grad_key and grad_value, over this algorithm's finite outputs. Given the forward's rounded out and logsumexp (any values of the right shapes; a non-finite logsumexp or output gives an infinite bound where it is used), for a fused backward kernel; without them, for a backward that recomputes the softmax. compute, ulps, flush_to_zero and standard are as in Result.tolerance.

backward_finite_range(*args, **kwargs) -> Outputs#: Per gradient, the lowest and highest finite output of this algorithm, as Result.finite_range describes; the arguments are backward_tolerance's.

backward_may_be(kind: str, *args, **kwargs) -> Outputs#: Per gradient, where this algorithm may output kind ("+inf", "-inf" or "nan"), as Result.may_be describes; the other arguments are backward_tolerance's.

backward_covers(gradients: dict, *args, **kwargs) -> Outputs#: For each gradient named in gradients (grad_query, grad_key, grad_value), where this algorithm can produce that value, as Result.covers describes; the other arguments are backward_tolerance's.

Base classes#

class Algorithm#: Base for named algorithms (see the module docstring).

class ErrorModel(relative: Fraction, absolute: Fraction, scale: str)#: |error| <= relative * scale + absolute, with scale naming the input quantity.

Accumulation methods#

Accumulation methods: how an implementation may add up a list of terms.

Each method is a fully specified summation order. emulate reproduces it with a rounding function for the accumulator format. For tolerances, a method reports depth (the most roundings any one term passes through) and roundings (the total), which give the standard worst-case bound

|computed - exact| <= gamma(depth) * sum|t_i| + roundings * eta * (1 + gamma(depth))

with gamma(n) = n*u / (1 - n*u) (Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., chapters 3-4). The exact method sums exactly and rounds once, which no implementation can beat.

class AccumulationMethod(key: 'str', label: 'str', description: 'str', order_dependent: 'bool', parameters: 'tuple[str, ...]' = (), references: 'tuple[str, ...]' = ())#

emulate(terms: list, round_: Callable, add: Callable = operator.add, **params)#: Sum terms in this method's order, calling round_ wherever the accumulator rounds.

Code

from exactref.accumulation import get_method

# The pairwise summation method
pairwise = get_method("pairwise")
# The most roundings any one of 8 terms passes through
print("depth:", pairwise.depth(8))
# The number of roundings in a sum of 8 terms
print("roundings:", pairwise.roundings(8))

Output

depth: 4
roundings: 15

Pairwise summation rounds each term, then adds the two halves recursively, rounding after every addition. For 8 terms, a term passes through its own rounding and 3 levels of additions (a depth of 4), and there are 8 + 7 = 15 roundings in all.