Unit 5.8 · about 3 h
87 / 324Residuals, tolerances and iterative-solver error
A solver stops because of a rule, not because mathematics declares the answer correct. Residual scaling and conditioning tell you what that rule actually certifies.
- derive
- measure
You can now
- Derive the relationship between a linear-system residual, solution error and conditioning.
- Design norm-specific absolute and relative residual tests with meaningful scales.
- Distinguish residual convergence, stagnation, non-finite failure and iteration-cap termination.
Covers
- Residuals, tolerances and iterative-solver error
Study, in this order
3 steps
- Linear residuals and conditioningFundamentals of Numerical Computation
Section 2.8.1, Theorem 2.8.1; Section 2.8.2, Residual and backward error, Definition 2.8.2 and equation 2.8.11
Depth · Keep the normwise residual/error relationship. Do not reproduce the later least-squares solver comparison.
- Fixed-point errors versus iteration stoppingFundamentals of Numerical Computation
Sections 4.2.1 to 4.2.3, equations 4.2.2 to 4.2.7, Theorem 4.2.1 and Corollary 4.2.1; then Section 4.3.2's
newtonstopping implementationDepth · Observation 4.2.1 says the iteration diverges for all initial values when \(|g'(p)|>1\); read that as a local statement, since starting at \(p\) stays at \(p\). Theorem 4.2.1 as printed allows any invariant interval, but existence needs completeness: \(x/2\) maps \((0,1)\) into itself with no fixed point there. Read it with \(S\) a closed interval, then finish the existence step yourself in the self-test. Read the scalar stopping loop on the Newton page, not its graphical demos.
- Objective certificates rather than tiny stepsConvex Optimization
Section 5.5.1 (pp. 241 to 242) and Section 9.1.2, equations 9.9 to 9.10 (p. 460) and the strong convexity constants paragraph (p. 463)
Depth · Read only primal-dual gap stopping and the strong-convexity gradient-to-suboptimality bound, including why an unknown m makes that bound a conceptual rather than practical stopping rule.
Focus on
- For invertible \(A\), \(r=b-A\hat{x}=A(x^*-\hat{x})\). Hence \(\|\hat{x}-x^*\|_2\le\|A^{-1}\|_2\|r\|_2\); when \(b\) and \(x^*\) are nonzero, relative error \(\le\kappa_2(A)\|r\|_2/\|b\|_2\).
- A scale-aware residual test is \(\|r\|_2\le\mathit{atol}+\mathit{rtol}\|b\|_2\). State the norm, units, \(\mathit{atol}\) and \(\mathit{rtol}\). The normalized residual \(\|r\|_2/(\|A\|_2\|\hat{x}\|_2+\|b\|_2)\) is a useful normwise backward-error measure, not a forward-error certificate.
- For a \(\mu\)-strongly-convex differentiable unconstrained objective, \(\|\operatorname{grad} f(x)\|_2^2/(2\mu)\) bounds the objective gap. Constraints and nonsmooth penalties need KKT or proximal residuals instead; noisy gradients are not exact stationarity residuals.
- Check finite values first, then residual convergence, then stagnation and the iteration limit. Small successive steps can reflect poor scaling or rounding. Report the reason for exit rather than labeling every exit converged.
Reproduce closed-book
float64 set A=diag(1,1e-8), xstar=np.ones(2,dtype=np.float64), b=A@xstar and xhat=np.array([1,0],dtype=np.float64). Compute r=b-A@xhat: \(\|r\|_2=10^{-8}\), relative residual=1e-8, relative solution error=0.7071067812, \(\kappa_2(A)=10^8\), and normalized residual=5e-9 using the denominator \(\|A\|_2\|\hat{x}\|_2+\|b\|_2\). Allow relative scalar tolerance 1e-9 for these displayed nonzero quantities. The relative-error bound is 1, so the large observed error does not contradict it. An absolute-only threshold 1e-7 accepts this poor answer. With \(\mathit{atol}=10^{-12}\) and \(\mathit{rtol}=10^{-10}\), the Euclidean residual test rejects it. Multiply both \(A\) and \(b\) by \(10^8\): the residual becomes 1 and the relative residual remains about 1e-8; the same combined test rejects it again. Take one correction xhat+=np.linalg.solve(A,r) using the original A and r. The measured final Euclidean error and residual are both 0; require \(\le 10^{-12}\) rather than exact equality. Finally work the harmonic partial sums \(H_n\) by hand: successive increments tend to zero although \(H_n\) diverges. Use this to justify a stagnation status separate from residual convergence. No random draws are needed.Traps
- Treating a tiny residual as a tiny solution error without checking the inverse sensitivity or condition number.
- Using a relative-only tolerance when the target scale is zero, or comparing residuals in different units without rescaling
atol. - Declaring convergence at
maxiteror after a small update even though the residual test failed.
Practice
FNC Python: foil the successive-step stopping indicator
Exercise 4.3.8
The harmonic sums in Reproduce are one answer; find a different foiling sequence and prove both properties. Explain why a shrinking update alone is insufficient, then specify finite-value, residual, stagnation and iteration-cap statuses for your own loop.
Self-test
-
Which assumptions are used in the relative residual bound?
Invertibility gives \(x^*-\hat{x}=A^{-1}r\). Consistent operator and vector norms give \(\|\text{error}\|\le\|A^{-1}\|\|r\|\) and \(\|b\|\le\|A\|\|x^*\|\). Combine them to get the relative bound; \(b\) and \(x^*\) must be nonzero to divide by their norms. No random or convergence assumption is needed. -
Why does multiplying A and b by the same nonzero scalar leave relative residual unchanged?
Both \(\|r\|\) and \(\|b\|\) acquire the scalar's absolute magnitude. The solution is unchanged, but an absolute residual tolerance must be rescaled to preserve the same test exactly. -
Why do vanishing harmonic increments not establish convergence?
\(H_{n+1}-H_n=1/(n+1)\) tends to zero. But each block from \(2^{k-1}+1\) through \(2^k\) contributes at least \(1/2\), because it has \(2^{k-1}\) terms each at least \(1/2^k\). Thus \(H_{2^k}\ge 1+k/2\), including \(k=0\), and the sums diverge. -
What does a primal-dual gap certify, and what must be checked first?
For primal-feasible \(x\) and dual-feasible multipliers, weak duality bounds \(f(x)-p^*\) by \(f(x)-g\). Check both feasibilities before interpreting the gap. Strong duality is needed to make arbitrarily tight such certificates possible, not for this lower-bound argument. -
What makes the contraction-mapping existence guarantee valid?
A nonempty complete domain, such as a closed interval of \(\mathbb{R}\), an invariant map \(g(S)\) contained in \(S\), and a global Lipschitz constant \(q<1\). Invariance keeps all iterates in \(S\); the contraction bounds successive differences by a geometric sequence, hence the iterates are Cauchy. Completeness puts the limit in \(S\), and continuity gives \(g(p)=p\). The contraction also proves uniqueness. An open interval alone does not provide completeness.
Completion criterion
Stop when
You can derive the conditioning bound, explain the small-residual counterexample, and defend a norm-and-scale-specific stopping rule that never labels stagnation or an iteration cap as residual convergence.
Mark passed, then open 1.1 Only mark this when you can do it without notes.Leads to
- 2.1 Least squares: residuals, normal equations and stable solvers Use residual and solution-error checks when comparing linear solvers.
- 8.5 Overflow, underflow and finding the first invalid value Report convergence and numerical failure as separate statuses.
Resources for this unit
- Convex Optimization Read Chapters 1-5 and 9 with the course; use Chapter 10 only for equality-constrained Newton steps.
- 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.
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.
- 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. - PyTorch numerical type information Read
torch.finfo'seps,maxandsmallest_normalfields and its subnormal note; inspectfloat64,float32,float16andbfloat16on CPU. - 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.