Unit 5.6 · about 4 h
85 / 324Floating point, rounding, overflow and cancellation
Mathematical real arithmetic is not a machine's arithmetic. A dtype's range and spacing predict different failures, and rewriting a formula can matter more than increasing precision.
- derive
- measure
You can now
- Derive spacing and unit roundoff from a binary significand and distinguish precision from exponent range.
- Compare
float64,float32,float16andbfloat16, including normal and subnormal values. - Construct overflow, underflow and cancellation failures and check an algebraically stable alternative.
Covers
- Floating-point representation, rounding and machine precision
- Overflow, underflow and catastrophic cancellation
Study, in this order
3 steps
- Binary representation and roundingFundamentals of Numerical Computation
Section 1.1: Definitions 1.1.1 to 1.1.2, equations 1.1.1 to 1.1.6, and Sections 1.1.2 to 1.1.3
Depth · Use the Python examples. The main-text minimum-range discussion concerns normal values; read footnote 2 and the
dtypereference below before treating \(\mathit{tiny}\) as the smallest positive representable value. Skip the Julia-specific Exercise 1.1.5. - Conditioning and unstable cancellationFundamentals of Numerical Computation
Sections 1.2.1 to 1.2.2, condition numbers and error estimates, on the preceding Problems and conditioning page; Section 1.4.1 and Examples 1.4.1 to 1.4.2 on Stability
Depth · Focus on subtractive cancellation and the stable quadratic-root rewrite; skip the later polynomial-root perturbation demo.
- Inspect low-precision formatsPyTorch numerical type information
torch.finfo:bits,eps,max,tiny,smallest_normaland the note about subnormal numbers
Focus on
- With \(p\) total significand bits including the hidden leading one, spacing in \([2^e,2^{e+1})\) is \(2^{e-(p-1)}\). \(\mathit{eps}\) is the gap above 1; unit roundoff for round-to-nearest is \(\mathit{eps}/2\). The relative bound needs normal-range results and no overflow.
float64,float32,float16andbfloat16have \(p=53,24,11,8\) and exponent fields of 11,8,5,8 bits respectively. bfloat16 has wide range but coarser spacing than float16.- Overflow produces infinities; gradual underflow first reaches subnormals with reduced relative precision. \(\mathit{tiny}\) is the smallest positive normal, not the smallest positive representable value.
- Subtracting nearby rounded operands can destroy relative accuracy even when the subtraction itself is exact. Separate conditioning of the requested difference from avoidable instability in how its operands were formed.
Reproduce closed-book
np.finfo for np.float64, np.float32 and np.float16, and torch.finfo(torch.bfloat16) on CPU for bfloat16. Print eps, tiny and max to three significant digits: float64 (2.22e-16,2.23e-308,1.80e308), float32 (1.19e-7,1.18e-38,3.40e38), float16 (9.77e-4,6.10e-5,6.55e4), bfloat16 (7.81e-3,1.18e-38,3.39e38). Derive \(\mathit{eps}=2^{-(p-1)}\) and check each \(\mathit{eps}\) within relative scalar tolerance 1e-6. Inspect np.nextafter(dtype(0),dtype(1)), and torch.nextafter on CPU bfloat16 scalars, to see values below \(\mathit{tiny}\). In NumPy float64 set x=np.float64(1e8): \(\sqrt{x^2+1}-x\) returns 0, while \(1/(\sqrt{x^2+1}+x)\) returns about 5.00e-9. Rationalize the expression by hand and bound the exact result's difference from \(1/(2x)\) by \(1/(8x^3)\); allow absolute 1e-20 for the stable scalar. Finally evaluate np.exp(np.float64(1000)) and np.exp(np.float64(-1000)) inside np.errstate(over='ignore',under='ignore'): outputs are inf and 0. These are range failures, not exact evaluations of \(\exp\). No random draws are needed.Traps
- Calling \(\mathit{eps}\) the smallest representable number, or calling
finfo.tinythe smallest subnormal. - Assuming bfloat16 is more precise because its exponent range is larger.
- Applying a relative error tolerance at a true zero, or claiming that a stable rewrite recovers information already rounded out of its inputs.
Practice
FNC Python: enumerate a toy format and derive single precision
Exercises 1.1.1 (a) to (c) and 1.1.4 (a) to (b)
Attempt by hand before using finfo or nextafter. Translate the single-precision checks to NumPy float32 rather than the exercise's Julia type name.
Self-test
-
Why is round-to-nearest error bounded by eps/2 in the normal range?
In each binade the absolute error is at most half its spacing. Dividing by \(|x|\ge 2^e\) bounds relative error by \(2^{-(p-1)}/2\). The lower bound on \(|x|\) uses normal binade representation; subnormals instead have fixed absolute spacing. Overflow is excluded. -
What is lost in \(\sqrt{x^2+1}-x\) at \(x=10^8\)?
\(x^2\) is about \(10^{16}\), wherefloat64spacing is 2, so adding 1 rounds away under ties-to-even. The square root returns \(x\) and subtraction gives zero. Rationalization avoids subtracting nearly equal operands; it does not make \(x^2+1\) exact. -
Is the computed zero an underflow of the true answer?
No. The true answer is about 5e-9, well within the normalfloat64range. It is cancellation after an earlier rounding loss. -
Can \((1+\mathit{eps}/2)-1\) and \(1+(\mathit{eps}/2-1)\) differ in
float64?Yes. Under round-to-nearest ties-to-even, the first evaluates to zero and the second to about 1.11e-16. The real-arithmetic associative law does not transfer to rounded operations.
Completion criterion
Stop when
You can derive each format's spacing, distinguish \(\mathit{eps}\) from normal and subnormal limits, and explain why the stable cancellation rewrite passes its absolute-error check while the naive formula fails.
Mark passed, then open 5.7 Only mark this when you can do it without notes.Leads to
- 8.5 Overflow, underflow and finding the first invalid value Use constructed failure mechanisms to interpret non-finite diagnostics.
- 8.4 Precision formats, autocast and gradient scaling Apply range-versus-precision reasoning when selecting dtypes.
Resources for this unit
- Fundamentals of Numerical Computation Study §§1.1-1.4, 2.7-2.8, 3.2-3.4 and 4.2-4.3 for floating point, conditioning, stability, least squares and iteration error, in the Python edition.
- PyTorch numerical type information Read
torch.finfo'seps,maxandsmallest_normalfields and its subnormal note; inspectfloat64,float32,float16andbfloat16on CPU.
Also in this section
- Stanford EE364a - Convex Optimization Follow lectures 1-9 and 15-16; use lecture 17 only for KKT systems. Leave full interior-point treatment for later.
- Convex Optimization Read Chapters 1-5 and 9 with the course; use Chapter 10 only for equality-constrained Newton steps.
- D2L - Optimization Algorithms Use the convexity, GD, SGD, minibatch, momentum, AdaGrad, RMSProp, Adam and scheduler sections for ML derivations and code.
- CVXPY basic and DCP tutorials Use as the solver oracle after deriving a constrained problem and its KKT conditions by hand.
- Proximal Algorithms Read Sections 1.1, 4.2 and the elementwise soft-thresholding formula in 6.5.2, then the lasso formulation on p. 196. Skip 7.1.1's penalty-scaling and Gram-rewrite errors on p. 197; derive the update from 4.2 and 6.5.2.
- PyTorch optimizer and gradient-clipping contracts Use AdamW's algorithm and
clip_grad_norm_'s parameters to locate weight decay and clipping precisely. Inspectstate_dictfor optimizer restart; no network training is assigned. - SciPy minimize: BFGS Read the BFGS method description and
gtol,maxiter,c1andc2options. Supply an analytic gradient and compare residuals and evaluation counts, not just iterations. - The Accuracy of Floating Point Summation Read Section 2's recursive-sum error bound and Section 3's pairwise and compensated summation algorithms. Leave the statistical error models and historical machine details aside.
- Accurately computing running variance Read the Welford mean and centered-sum recurrences. Prove the centered-sum invariant and translate the recurrence to NumPy; skip the C++ class.