What do multiples of machine epsilon mean?

Question

I'm trying to investigate the source of residuals of a matrix equation problem (Ax=b). To verify my answer, I subtract Ax-b, expecting 0. Instead of "pure" zeros, I obtain values in the same order of magnitude as machine epsilon, which is fine. The problem, is that these residuals seem to be multiples of each other, so I'm unsure how to interpret them.

I found some details here: Machine epsilon computation issue, this didn't clarify why multiples of the epsilon arise instead of just one or the other.

I checked on my system using np.finfo(float).eps which produced 2.220446049250313e-16. One of the residuals I get in the solution x is identical to this value, however, the other seems to be half the epsilon.

Here's the code I used:

# Arbitrary Matrix A and Vector b
A = np.array([[2,-1,0],[1,-2,1],[0,-1,2]])
b = np.array([[1],[0],[1]])

# Solve for Vector x
x = np.linalg.solve(A,b)

# Calculate difference, expected to be column of zeros
diff = A.dot(x) - b
print(diff)

Here's the output:

Output: 
[[ 0.00000000e+00]
 [-1.11022302e-16]  #-------> Is this machine epsilon...
 [-2.22044605e-16]] #-------> ...or this?

What is the explanation/interpretation for this? I understand that values smaller than epsilon can still be represented, but in that case why aren't both residuals -1.11022302e-16?

Thanks in advance!

Answer 1

The so-called machine epsilon is just the unit of least precision (ULP) at 1. That is, it is the position value of the least significant bit in the representation of 1. When there are 53 bits in the significand, 1 is represented by the binary numeral 1.000…000₂, where there are 52 zeros after the binary point. So the position value of the lowest digit is 2⁻⁵², and 2⁻⁵² is the ULP of 1.

In general, let ULP(x) stand for the unit of least precision for x. Generally, a floating-point format represents a number as (−1)^s • f • b^e, where b is a fixed base (two for the binary formats, ten for decimal, 16 for hexadecimal), s is a sign bit (0 for +, 1 for −), e is an exponent, and f is a significand with p digits, where p is a fixed amount for the format. For IEEE-754 binary32, p is 53, for 53 bits. The ULP is the position value of the lowest precision in the significand as scaled by the exponent, so, if some number x is represented in the floating-point format with sign bit s, significand f, and exponent e, its ULP is b^1−p • b^e. (I have assumed the format of the significand is one base-b digit before a radix point and p−1 digits after the radix point, and that is why its lowest digit has a position value of b^1−p. Such significands are in the interval [1, b). Sometimes significands are scaled differently, and the exponent is adjusted to compensate. For example, it may be useful in proofs for the significand to be an integer.)

In binary formats, ULP(2) = 2•ULP(1), ULP(½) = ½•ULP(1), ULP(¼) = ¼•ULP(1), and so on.

Suppose you have computed two values that are in the interval [1, 2) and that would, if computed with real-number arithmetic, be equal, but they have been computed with floating-point arithmetic and happen to differ slightly. Because of the format of the representation, they can only differ by multiples of ULP(1). When you subtract such numbers, you will often get 0, ULP(1), 2•ULP(1), or some other multiple of ULP(1), depending on the circumstances. When two numbers that would be identical if computed with real-number arithmetic are computed with floating-point arithmetic, they may experience different rounding errors in various parts of the computation.

If you compute two values that are in the interval [½, 1), they can differ only by multiples of ULP(½).

This is why you are seeing various multiples or binary fractions of ULP(1). It is simply an artifact of the the quantization of the floating-point format.