The Binary Representation Problem
When you solve a system of linear equations by hand, you work with clean fractions like 1/3 or 2/5. However, computers do not store numbers as fractions by default. Instead, they use the IEEE 754 standard for floating-point arithmetic, which represents numbers in binary (base-2).
Just as the fraction 1/3 cannot be represented exactly in our decimal system (it becomes the repeating decimal 0.3333...), fractions like 1/10 or 1/5 cannot be represented exactly in binary. For example, 0.1 in binary is a repeating fraction: 0.00011001100110011.... Because computers have a finite number of bits (usually 64 bits for double-precision), they must truncate this infinite sequence, introducing a tiny rounding error of about 1 part in 1016.
How Row Operations Magnify Round-off Errors
In algorithms like numpy.linalg.solve or the least-squares fitting of numpy.polyfit, the computer performs Gaussian elimination or QR factorization. These processes involve thousands of arithmetic operations: multiplying rows, subtracting them, and dividing by pivot values.
During these steps, the tiny truncation errors from binary representation are added, subtracted, and multiplied. If an algorithm divides a row by a very small pivot, the round-off error is magnified exponentially. By the time the algorithm reaches the final solution, an entry that should be exactly 0 mathematically might be computed as -1.3877787807814457e-16. An entry that should be exactly 1 might become 1.0000000000000002. This decimal "junk" is a direct result of these accumulated errors.
Python Example: Visualizing the Floating-Point Gap
Consider this simple Python script that solves a basic system of equations:
import numpy as np
# System:
# 3x + 2y = 1
# 1x - 2y = 0
A = np.array([[3.0, 2.0], [1.0, -2.0]])
b = np.array([1.0, 0.0])
x = np.linalg.solve(A, b)
print(x) # Output: [0.25 0.125]
For this system, the values are clean (1/4 and 1/8) and can be represented exactly in binary. But if we change the coefficients slightly to introduce repeating binary fractions (like 1/3), NumPy will return floats with trailing junk. For exact algebraic results, you must use symbolic math libraries like SymPy or a dedicated fraction-based matrix calculator.
How to Deal with Floating-Point Junk
If you are coding a program that depends on matrix reduction, you can handle floating-point junk using these strategies:
- Thresholding (Tolerances): Instead of checking if an entry is exactly zero (
val == 0), check if its absolute value is smaller than a threshold (e.g.,abs(val) < 1e-12). NumPy providesnp.isclose()andnp.allclose()for this purpose. - Rational Math: Use a rational number class (like Python's
fractions.Fraction) that stores numbers as integer numerators and denominators. This eliminates binary rounding errors completely at the cost of execution speed.