add and multiply#

Elementwise addition and multiplication combine arrays element by element: x[i] + alpha * y[i] and x[i] * y[i]. Broadcasting first stretches arrays of different shapes to a common shape, as NumPy and PyTorch do, so a row can be added to every row of a matrix. They are the most basic operations of all, and they still round: the exact sum or product of two floats usually needs more bits than the format has. A single product, or a sum with alpha 1, rounds once, so there is usually exactly one correct answer; a sum with another alpha, and a complex product, round in more places.

In ExactRef, er.add(x, y, alpha=1) computes x + alpha * y and er.multiply(x, y) computes x * y, element by element with NumPy broadcasting. They are also exactref.arithmetic.basic.add and exactref.arithmetic.basic.multiply.

Each output element is the exact value rounded once. Both are matmul's exact sums underneath (x * y a one-term product, x + alpha * y the two-term sum [x, y] . [1, alpha]), so they share matmul's tolerances and normal answer. With at most two exact terms, the tolerance is the exact worst case, and often 0.

Arguments#

Argument

Meaning

Default

x, y

the operands, broadcast together; a Python scalar is taken exactly

required

alpha (add only)

the factor on y, taken exactly

1

dtype

the output format, of the same kind as the inputs

the operands' promoted format

rounding

the rounding of the reference and of the final step into the output format

er.Rounding.NEAREST_EVEN

overflow

what a result past the output format's range becomes: "ieee", "saturate" or "nan" for floats; "wrap", "error" or "saturate" for integers

each format's own rule; integers wrap

Examples#

Code
import exactref as er

# A shorthand for these examples: the values as a float32 array
def f32(values):
    return er.asarray(values, "float32")

# A float32 vector
x = f32([1.0, 2.0, 3.0])

# x + [0.1, 0.2, 0.3], each element the exact sum rounded once to float32
print("add:", er.add(x, f32([0.1, 0.2, 0.3])).rounded.tolist())
# Each product rounds once, so the tolerance here is 0
print("multiply tolerance:",
      er.multiply(f32([0.1, 0.2]), f32([0.3, 0.7])).tolerance().tolist())
Output
add: [1.100000023841858, 2.200000047683716, 3.299999952316284]
multiply tolerance: [0.0, 0.0]

Broadcasting#

Code
# A column of 2 times a row of 3: both stretch to the shape (2, 3)
print("product:",
      er.multiply(f32([[1.0], [2.0]]), f32([10.0, 20.0, 30.0])).rounded.tolist())
Output
product: [[10.0, 20.0, 30.0], [20.0, 40.0, 60.0]]

alpha and scalar operands#

alpha and a scalar operand (an int, float, complex, a Fraction with a power-of-two denominator, a FloatValue or a ComplexValue) are taken exactly and do not change the output format:

Code
# 1 + 3 * y with y the float32 nearest 0.1, rounded once
print("alpha=3:", er.add(f32([1.0]), f32([0.1]), alpha=3).rounded.tolist())
# x times the Python float 0.5; the output stays float32
print("x * 0.5:", er.multiply(x, 0.5).rounded.tolist())
Output
alpha=3: [1.2999999523162842]
x * 0.5: [0.5, 1.0, 1.5]

Mixed formats#

Operands of different formats are promoted the way PyTorch promotes them. An operand of shape () takes part in the promotion only when every operand has shape (), and one in a wider format is first rounded once to the promoted format, as PyTorch casts it:

Code
# The output format of x * y
def promoted(x, y):
    return er.multiply(x, y).rounded.dtype.name

# Two arrays of shape (1,): float32 is the narrowest format that holds both
print("bfloat16 and float16:",
      promoted(er.asarray([1.5], "bfloat16"), er.asarray([2.0], "float16")))
# A float64 of shape () beside an array of shape (1,): only the array counts
print("float64 of shape () and float16:",
      promoted(er.asarray(1.5, "float64"), er.asarray([2.0], "float16")))
# Both of shape (): both formats take part
print("both of shape ():",
      promoted(er.asarray(1.5, "float64"), er.asarray(2.0, "float16")))
Output
bfloat16 and float16: float32
float64 of shape () and float16: float16
both of shape (): float64

Signs of zero#

An exact zero takes the sign IEEE arithmetic gives, as in PyTorch: -0.0 + -0.0 is -0, and a cancellation is +0, or -0 when rounding toward -inf.

Code
# A float32 -0
neg = f32([-0.0])
# value((0,)) is element 0's exact value; negative is True for a minus sign
print("-0 + -0 negative:", er.add(neg, neg).rounded.value((0,)).negative)

# Float32 1 and -1
one, minus_one = f32([1.0]), f32([-1.0])
# 1 + -1 cancels to an exact zero
print("1 + -1 negative:", er.add(one, minus_one).rounded.value((0,)).negative)
# The same sum, rounded toward -inf
down = er.Rounding.TOWARD_NEGATIVE
print("1 + -1 toward -inf negative:",
      er.add(one, minus_one, rounding=down).rounded.value((0,)).negative)
Output
-0 + -0 negative: True
1 + -1 negative: False
1 + -1 toward -inf negative: True

Special values#

Complex special values follow IEEE rules in each part, as in PyTorch. A coefficient that is exactly 0, such as alpha=0, turns an infinity into NaN, as PyTorch's does. Where mathematics gives an answer, ExactRef gives it even when PyTorch differs: a zero part of a complex alpha adds no term.

Code
# A shorthand: the values as a complex128 array
def c128(values):
    return er.asarray(values, "complex128")

# The imaginary part is inf * 0 + 0 * 1, and inf * 0 is NaN
print("(inf+0j) * (1+0j):",
      er.multiply(c128([complex("inf")]), c128([1 + 0j])).rounded.tolist())
# Each part is its own sum: inf + 1 and 1 + 1
print("(inf+1j) + (1+1j):",
      er.add(c128([complex("inf+1j")]), c128([1 + 1j])).rounded.tolist())
# alpha=0 gives the term 0 * inf, which is NaN
print("1 + 0 * inf:",
      er.add(f32([1.0]), f32([float("inf")]), alpha=0).rounded.tolist())
# alpha's real part is 0, so the term 0 * inf is left out
print("(1+1j) + 1j * (inf+1j):",
      er.add(c128([1 + 1j]), c128([complex("inf+1j")]), alpha=1j).rounded.tolist())
Output
(inf+0j) * (1+0j): [(inf+nanj)]
(inf+1j) + (1+1j): [(inf+2j)]
1 + 0 * inf: [nan]
(1+1j) + 1j * (inf+1j): [infj]

PyTorch gives nan+infj for the last one, because it multiplies alpha * y out in complex arithmetic.

Checking PyTorch#

Code
import torch

# Two float32 tensors
a = torch.tensor([0.1, 0.2, 0.3])
b = torch.tensor([0.7, 0.8, 0.9])

# Could a correct implementation produce PyTorch's a + 0.5 * b? covers() answers
# per element, and all() asks whether every element is covered.
print("add covered:",
      er.add(a, b, alpha=0.5).covers(torch.add(a, b, alpha=0.5)).all())
# The same question for PyTorch's a * b
print("multiply covered:", er.multiply(a, b).covers(a * b).all())
Output
add covered: True
multiply covered: True

PyTorch converts alpha before using it (on the CPU, to the tensor's dtype), and in bfloat16 and float16 it may store alpha * y first. Pass the converted alpha, and use any_of when either form may run; PyTorch equivalents has the details and an example.

Backward#

er.add.backward(x, y, grad, alpha=) returns grad_x = grad and grad_y = conj(alpha) * grad; er.multiply.backward(x, y, grad) returns grad_x = grad * conj(y) and grad_y = grad * conj(x). Each is summed over its operand's broadcast copies and rounded once.

Code
# A column of 2 and a row of 3, as in the broadcasting example
col, row = f32([[1.0], [2.0]]), f32([10.0, 20.0, 30.0])
# The reference backward of col * row, for an output gradient of ones
g = er.multiply.backward(col, row, f32([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]))
# grad times row, summed over the 3 copies of each element of col: 10 + 20 + 30
print("grad_x:", g.rounded.grad_x.tolist())
# grad times col, summed over the 2 copies of each element of row: 1 + 2
print("grad_y:", g.rounded.grad_y.tolist())
Output
grad_x: [[60.0], [60.0]]
grad_y: [3.0, 3.0, 3.0]

Each gradient is summed over the copies broadcasting made of its operand. A scalar operand gets no gradient, so the result then holds only the other operand's gradient:

Code
# The scalar operand 3.0 gets no gradient, so rounded is grad_x alone
print("grad_x:",
      er.multiply.backward(f32([1.0, 2.0]), 3.0, f32([1.0, 1.0])).rounded.tolist())
Output
grad_x: [3.0, 3.0]

Signs of zero in gradients#

A gradient element with one copy, of an operand with the output's number of dimensions, is a single product, so an exact zero takes that product's IEEE sign, as in PyTorch. Every other gradient element is a sum, and its exact zero is +0 (-0 when rounding toward -inf), as PyTorch's sums start from +0; PyTorch also sums a gradient over the dimensions broadcasting adds in front. In add's backward, a complex alpha other than 1 follows the coefficient rule, where PyTorch multiplies the gradient by conj(alpha) in complex arithmetic, so a zero part can differ in sign. gemm's grad_C with a complex beta other than 1 is the same case.

Full specification#

er.add#

add(x, y, *, alpha=1, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Elementwise x + alpha * y with broadcasting (NumPy rules).

Each output element is the exact value rounded once to dtype (default: the operands' promoted format, with operands of shape () as in multiply). alpha is taken exactly as given, like a scalar operand. tolerance() budgets alpha * y rounded in the compute format, the sum rounded in it, and the final rounding; with alpha 1 or -1 nothing is rounded but the sum. An implementation that stores alpha * y in a narrower format first is the chain add(x, multiply(y, alpha, dtype=...)); any_of covers both when either may run. PyTorch's CPU add also rounds alpha to the dtype first, so to check it pass that rounded alpha: then float32 matches the plain form, bfloat16 the stored form, and float16 either one, element by element (any_of). Integer results wrap on overflow unless overflow says otherwise. An exact zero takes the sign IEEE arithmetic gives x + alpha * y, as in PyTorch: -0.0 + -0.0 is -0, and a cancellation is +0 (-0 when rounding toward -inf). In a complex sum, a zero part of alpha (the imaginary part of 1, say) adds no term, by the rule below, where PyTorch multiplies alpha * y out: (-0+1j) + (-0-1j) has the real part -0, PyTorch's +0. A complex non-finite part reaches only the parts it is added to, as in PyTorch: (inf + 1j) + (1 + 1j) is inf + 2j. A complex alpha follows the same rule, which is the mathematical answer: add(1+1j, inf+1j, alpha=1j) is infj, where PyTorch, multiplying alpha * y in complex arithmetic, gives nan+infj. x and y may also be earlier results.

er.multiply#

multiply(x, y, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Elementwise product with broadcasting (NumPy rules).

Each output element is the exact product rounded once to dtype (default: the operands' promoted format: the one that holds both, else the narrowest standard format that does, as PyTorch promotes). As in PyTorch, an operand of shape () takes part in the promotion only when both are of shape (); one in a wider format is first rounded once to the promoted format. A scalar operand (int, float, complex, a Fraction with a power-of-two denominator, a FloatValue or ComplexValue) is taken exactly and does not change the format. tolerance() budgets the product's rounding in the compute format and the final rounding; a complex product's parts are two-term sums, budgeted as in matmul. Integer results wrap on overflow unless overflow says otherwise. A complex product with a non-finite part follows IEEE rules in each part, as PyTorch's does: (inf + 0j) * (1 + 0j) is inf + nanj. An exact zero takes the sign IEEE arithmetic gives the product (-0.0 * 1 is -0), as in PyTorch. x and y may also be earlier results, as Result describes.

er.add.backward#

backward(x, y, grad, *, alpha=1, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Gradients of add(x, y, alpha=alpha) for the output gradient grad.

grad_x = grad and grad_y = conj(alpha) * grad, each summed over the copies that broadcasting made of its operand, in that operand's shape and dtype; each is one product, so every element is one exact sum rounded once. An element with one copy, of an operand with the output's number of dimensions, is a single product, so an exact zero takes that product's IEEE sign. Any other element is a sum, as in PyTorch, which also sums the dimensions broadcasting adds in front: its exact zero is +0 (-0 when rounding toward -inf). A zero part of a complex conj(alpha) adds no term, as in add, where PyTorch multiplies grad by conj(alpha) in complex arithmetic unless alpha is 1: with alpha=2, a gradient 1-0j gives grad_y 2-0j, PyTorch's 2+0j. A scalar operand gets no gradient. grad must have the output's shape.

er.multiply.backward#

backward(x, y, grad, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Gradients of multiply(x, y) for the output gradient grad.

grad_x = grad * conj(y) and grad_y = grad * conj(x), each summed over the copies that broadcasting made of its operand, in that operand's shape and dtype. Each is one product (the copies are folded into its sum), so every element is one exact sum rounded once. A scalar operand gets no gradient. Where an operand is broadcast, PyTorch rounds each product to the dtype before summing the copies, an extra rounding the default tolerance does not budget. An element with one copy, of an operand with the output's number of dimensions, is a single product, so an exact zero takes that product's IEEE sign, as in PyTorch. Any other element is a sum, as in PyTorch, which also sums the dimensions broadcasting adds in front: its exact zero is +0 (-0 when rounding toward -inf). grad must have the output's shape.

See also#