Error Analysis

Error Analysis

Definitions

Formulas

Types of Errors

Truncation Error

Roundoff Error

Combined Error Model

The total error is approximately: total_error ≈ truncation_error + roundoff_error

For a method with truncation error O(h^p) and step size h:

• Decreasing h reduces truncation error but increases roundoff error (more operations)
• Optimal h balances both error sources
• Too small h leads to subtractive cancellation and roundoff dominates

Theorems

Python Implementation

import numpy as np
from typing import Callable, Tuple


class ErrorAnalyzer:
    """Comprehensive error analysis toolkit."""

    @staticmethod
    def condition_number_matrix(A: np.ndarray) -> float:
        """Compute 2-norm condition number of a matrix."""
        return np.linalg.cond(A, ord=2)

    @staticmethod
    def relative_error(computed: float, exact: float) -> float:
        """Compute relative error."""
        if exact == 0:
            return abs(computed)
        return abs(computed - exact) / abs(exact)

    @staticmethod
    def finite_difference_error(f: Callable, df: Callable, x: float,
                                h_range: Tuple[float, float] = (1e-16, 1)):
        """Analyze truncation vs roundoff error in finite differences."""
        h_values = np.logspace(np.log10(h_range[0]), np.log10(h_range[1]), 200)
        errors = []

        exact = df(x)
        for h in h_values:
            # Central difference
            approx = (f(x + h) - f(x - h)) / (2 * h)
            error = abs(approx - exact)
            errors.append(error)

        errors = np.array(errors)
        optimal_h = h_values[np.argmin(errors)]

        return h_values, errors, optimal_h

    @staticmethod
    def demonstrate_cancellation():
        """Demonstrate catastrophic cancellation."""
        print("\n--- Catastrophic Cancellation Demo ---")
        x = 1.0
        h_values = np.logspace(-1, -16, 16)

        print(f"Computing derivative of f(x) = sin(x) at x = 1")
        exact = np.cos(1.0)

        for h in h_values:
            # Forward difference (suffers from cancellation)
            fwd = (np.sin(1 + h) - np.sin(1)) / h
            # Central difference (better but still has issues)
            cen = (np.sin(1 + h) - np.sin(1 - h)) / (2 * h)

            fwd_err = abs(fwd - exact)
            cen_err = abs(cen - exact)

            print(f"  h={h:.0e}: Forward err={fwd_err:.2e}, "
                  f"Central err={cen_err:.2e}")

    @staticmethod
    def welford_online(data: np.ndarray) -> Tuple[float, float]:
        """Welford's online algorithm for mean and variance."""
        n = 0
        mean = 0.0
        M2 = 0.0

        for x in data:
            n += 1
            delta = x - mean
            mean += delta / n
            delta2 = x - mean
            M2 += delta * delta2

        variance = M2 / n if n > 1 else 0.0
        return mean, variance

    @staticmethod
    def naive_mean_variance(data: np.ndarray) -> Tuple[float, float]:
        """Naive two-pass mean and variance (unstable)."""
        n = len(data)
        mean = np.mean(data)
        variance = np.mean((data - mean) ** 2)
        return mean, variance

    @staticmethod
    def demonstrate_stability():
        """Compare stable vs unstable algorithms."""
        print("\n--- Algorithm Stability Demo ---")

        # Generate data with large mean and small variance
        np.random.seed(42)
        data = np.random.uniform(1e7, 1e7 + 1, 1_000_000)
        true_mean = np.mean(data)
        true_var = 1/12  # Variance of uniform(0,1)

        # Naive approach
        naive_mean, naive_var = ErrorAnalyzer.naive_mean_variance(data)

        # Welford
        welford_mean, welford_var = ErrorAnalyzer.welford_online(data)

        print(f"True mean:       {true_mean:.10f}")
        print(f"Naive mean:      {naive_mean:.10f}")
        print(f"Welford mean:    {welford_mean:.10f}")
        print(f"Mean error (naive):   {abs(naive_mean - true_mean):.2e}")
        print(f"Mean error (Welford): {abs(welford_mean - true_mean):.2e}")

        # Variance comparison
        naive_var_pass1 = np.mean(data**2) - np.mean(data)**2  # One-pass (unstable)
        print(f"\nNaive one-pass var: {naive_var_pass1:.6f}")
        print(f"Naive two-pass var: {naive_var:.6f}")
        print(f"Welford var:        {welford_var:.6f}")
        print(f"True var:           {true_var:.6f}")

    @staticmethod
    def demonstrate_condition_number():
        """Show effect of condition number on linear system solving."""
        print("\n--- Condition Number Effect on Linear Systems ---")

        # Well-conditioned
        A1 = np.array([[1, 0], [0, 2]])
        b1 = np.array([1, 1])
        x1 = np.linalg.solve(A1, b1)

        # Ill-conditioned
        A2 = np.array([[1, 1], [1, 1.0001]])
        b2 = np.array([2, 2.0001])
        x2 = np.linalg.solve(A2, b2)

        print(f"Well-conditioned (κ={np.linalg.cond(A1):.1f}):")
        print(f"  Solution: {x1}")

        print(f"Ill-conditioned (κ={np.linalg.cond(A2):.0f}):")
        print(f"  Solution: {x2}")

        # Perturbation experiment
        print(f"\nPerturbation experiment:")
        delta_b = np.array([1e-10, -1e-10])

        x1_pert = np.linalg.solve(A1, b1 + delta_b)
        x2_pert = np.linalg.solve(A2, b2 + delta_b)

        print(f"  Well-conditioned: Δx = {np.linalg.norm(x1_pert - x1):.2e}")
        print(f"  Ill-conditioned:  Δx = {np.linalg.norm(x2_pert - x2):.2e}")

    @staticmethod
    def demonstrate_backward_error():
        """Show backward error interpretation."""
        print("\n--- Backward Error Demonstration ---")

        # Solve Ax = b with known solution
        A = np.array([[1, 2], [3, 4]], dtype=float)
        x_true = np.array([1, 1])
        b = A @ x_true

        # Solve with perturbed matrix
        delta_A = np.array([[1e-10, 0], [0, 1e-10]])
        A_pert = A + delta_A
        x_computed = np.linalg.solve(A_pert, b)

        # Backward error: what perturbation to b would give this answer?
        delta_b = b - A @ x_computed
        backward_error = np.linalg.norm(delta_b) / np.linalg.norm(b)

        print(f"True solution:     {x_true}")
        print(f"Computed solution: {x_computed}")
        print(f"Forward error:     {np.linalg.norm(x_computed - x_true):.2e}")
        print(f"Backward error:    {backward_error:.2e}")
        print(f"Condition number:  {np.linalg.cond(A):.1f}")


def full_error_analysis():
    """Run comprehensive error analysis demonstrations."""
    print("=" * 60)
    print("NUMERICAL ERROR ANALYSIS")
    print("=" * 60)

    analyzer = ErrorAnalyzer()

    # 1. Condition numbers
    print("\n--- Condition Number Examples ---")
    matrices = {
        "Identity": np.eye(3),
        "Hilbert (5×5)": np.linalg.inv(np.vander(np.arange(5), increasing=True)),
        "Random 3×3": np.random.randn(3, 3),
        "Near-singular": np.array([[1, 1], [1, 1.0001]]),
    }

    for name, A in matrices.items():
        kappa = analyzer.condition_number_matrix(A)
        print(f"  {name:<20}: κ = {kappa:.4e}")

    # 2. Finite difference error
    print("\n--- Finite Difference Error Analysis ---")
    f = np.sin
    df = np.cos
    x = 1.0

    h_vals, errors, optimal_h = analyzer.finite_difference_error(f, df, x)
    print(f"  Optimal step size h: {optimal_h:.2e}")
    print(f"  Minimum error:       {np.min(errors):.2e}")
    print(f"  Theoretical optimal: {np.sqrt(np.finfo(float).eps):.2e}")

    # 3. Cancellation
    analyzer.demonstrate_cancellation()

    # 4. Stability comparison
    analyzer.demonstrate_stability()

    # 5. Condition number effect
    analyzer.demonstrate_condition_number()

    # 6. Backward error
    analyzer.demonstrate_backward_error()


full_error_analysis()

Applications in AI/ML

Error Sources in ML

Condition Numbers in Optimization

The condition number of the Hessian matrix determines convergence speed of gradient descent:

• κ(H) = 1: All directions have equal curvature → fastest convergence
• κ(H) ≫ 1: Elongated loss landscape → slow convergence, oscillation
• Preconditioning reduces effective κ, accelerating convergence

Numerical Precision in Deep Learning

Gradient Clipping as Stability Control

Gradient clipping (max norm or value clipping) is the ML equivalent of step-size control in ODE solvers. It prevents the instability that arises when the effective condition number of the optimization problem becomes too large during training.

Common Mistakes

Interview Questions

Q1: What is the difference between forward and backward error? A: Forward error measures how far the computed answer is from the true answer: |x̃ − x|. Backward error measures the smallest perturbation to the input that would produce the computed answer: find δ such that f(x + δ) = x̃. Backward error is often more meaningful because it tells you whether the computed answer is correct for a slightly different problem. A backward-stable algorithm produces the exact answer for a nearby problem.

Q2: Explain condition number and its practical significance. A: The condition number κ measures problem sensitivity. For a linear system Ax = b, κ(A) = ‖A‖·‖A⁻¹‖. If κ = 10^k, you lose about k digits of accuracy in the solution. For example, κ = 10^6 in double precision (15 digits) leaves about 9 correct digits. Large κ means the problem is ill-conditioned — no algorithm can produce accurate results without extra precision.

Q3: Why is Gaussian elimination with partial pivoting backward stable but the normal equations are not? A: Gaussian elimination with partial pivoting has a small backward error bound: ‖δA‖/‖A‖ ≤ n·ρ·ε. The normal equations (A^T A)x = A^T b square the condition number: κ(A^T A) = κ(A)². This means solving via normal equations loses twice as many digits, making it numerically unstable for ill-conditioned problems. Use QR factorization instead.

Q4: How do you choose the optimal step size for finite differences? A: The error has two components: truncation error O(h) and roundoff error O(ε/h). The total error is minimized when h_opt ≈ √(ε·|x|) for forward differences, or h_opt ≈ ε^(1/3)·|x| for central differences. For float64 (ε ≈ 10⁻¹⁶), h_opt ≈ 10⁻⁸ for forward and h_opt ≈ 10⁻⁵·|x| for central differences.

Q5: What is catastrophic cancellation and how do you detect it? A: Catastrophic cancellation occurs when subtracting nearly equal numbers, causing significant digits to cancel and amplifying relative error. Detect it by computing the condition number of the subtraction: κ = (|a| + |b|)/|a − b|. If κ ≫ 1, cancellation is occurring. Solutions: reformulate to avoid the subtraction, use higher precision, or apply algebraic identities (e.g., a² − b² = (a−b)(a+b)).

Q6: Explain the connection between condition number and floating point precision. A: If a problem has condition number κ, and the input has relative error ε, the output has relative error approximately κ·ε. In float64 with ε ≈ 10⁻¹⁶, a problem with κ = 10^10 loses 10 digits, leaving about 5 correct digits. If κ > 10^16, no digits are correct. This is why checking κ before solving is essential — it determines whether double precision is sufficient.

Practice Problems

import numpy as np


def naive_f(x):
    """Naive implementation — suffers from cancellation."""
    return (1 - np.cos(x)) / x**2


def stable_f(x):
    """Stable implementation using half-angle identity:
    1 - cos(x) = 2 sin²(x/2)
    """
    return 2 * np.sin(x / 2)**2 / x**2


# Test near x = 0 where cancellation is severe
print("f(x) = (1 - cos(x)) / x² near x = 0")
print(f"{'x':<14} {'Naive':<18} {'Stable':<18} {'True (Taylor)':<18}")
print("-" * 68)

# True value: f(0) = 1/2 (from Taylor expansion: cos(x) ≈ 1 - x²/2 + x⁴/24)
for x in [1, 0.1, 0.01, 1e-5, 1e-8, 1e-12, 1e-15]:
    naive = naive_f(x)
    stable = stable_f(x)
    # Taylor series: (1 - cos(x))/x² ≈ 1/2 - x²/24 + x⁴/720 - ...
    true_val = 0.5 - x**2 / 24 + x**4 / 720

    print(f"{x:<14.0e} {naive:<18.12f} {stable:<18.12f} {true_val:<18.12f}")

# Condition number of subtraction 1 - cos(x)
print("\nCondition number of 1 - cos(x):")
for x in [1, 0.1, 0.01, 1e-5, 1e-8]:
    a, b = 1, np.cos(x)
    kappa = (abs(a) + abs(b)) / abs(a - b)
    print(f"  x = {x:.0e}: κ = {kappa:.2e}")

import numpy as np


def naive_polynomial(coeffs, x):
    """Naive evaluation: Σ c_i * x^i."""
    result = 0.0
    for i, c in enumerate(coeffs):
        term = c * x**i  # Computing x**i from scratch each time
        result += term
    return result


def horner(coeffs, x):
    """Horner's method: (...((c_n * x + c_{n-1}) * x + c_{n-2}) * x + ...)"""
    result = 0.0
    for c in reversed(coeffs):
        result = result * x + c
    return result


# Test polynomial: (x-1)^10 expanded
coeffs = np.poly([1]*10)  # Coefficients of (x-1)^10

x = 1.00001  # Near the root
print("Polynomial: (x-1)^10 evaluated at x = 1.00001")
print(f"Naive:    {naive_polynomial(coeffs, x):.10f}")
print(f"Horner:   {horner(coeffs, x):.10f}")
print(f"Direct:   {(x - 1)**10:.10f}")

# Error analysis for large n
print(f"\nError growth analysis for x = 1 + h:")
for h in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]:
    x = 1 + h
    naive_val = naive_polynomial(coeffs, x)
    horner_val = horner(coeffs, x)
    true_val = h**10
    print(f"  h={h:.0e}: Naive err={abs(naive_val - true_val):.2e}, "
          f"Horner err={abs(horner_val - true_val):.2e}")

import numpy as np


def unstable_recurrence(n_terms):
    """Naive recurrence — unstable."""
    a = np.zeros(n_terms)
    a[0] = 0
    a[1] = 1

    for i in range(2, n_terms):
        a[i] = 2.25 * a[i-1] - 0.5 * a[i-2]

    return a


def stable_recurrence(n_terms):
    """Stable recurrence using correct characteristic roots.
    Roots of r² - 2.25r + 0.5 = 0 are r₁ = 2, r₂ = 0.25.
    General solution: aₙ = c₁·2ⁿ + c₂·(0.25)ⁿ
    With a₀=0, a₁=1: c₁ = 4/7, c₂ = -4/7
    """
    a = np.zeros(n_terms)
    c1 = 4/7
    c2 = -4/7

    for i in range(n_terms):
        a[i] = c1 * 2**i + c2 * (0.25)**i

    return a


print("Recurrence Relation Stability Analysis")
print(f"{'n':<6} {'Unstable':<18} {'Stable':<18} {'True':<18}")
print("-" * 60)

unstable = unstable_recurrence(20)
stable = stable_recurrence(20)

for i in range(min(15, len(unstable))):
    true_val = (4/7) * 2**i - (4/7) * (0.25)**i
    print(f"{i:<6} {unstable[i]:<18.6f} {stable[i]:<18.6f} {true_val:<18.6f}")

print(f"\nAt n=19: Unstable = {unstable[19]:.6f}, Stable = {stable[19]:.6f}")
print(f"True n=19: {(4/7)*2**19 - (4/7)*(0.25)**19:.6f}")

Quick Reference

Cross-References

• 086-Floating-Point — Roundoff error is defined by IEEE 754 properties; machine epsilon determines achievable accuracy
• 087-Root-Finding — Condition numbers determine sensitivity of root locations; Newton's method convergence depends on f''(x*)/f'(x*)
• 088-Integration — Quadrature error analysis uses truncation error bounds; adaptive quadrature controls error via subdivision
• 089-Interpolation — Interpolation error depends on derivatives of the function and node distribution
• Linear Algebra — Condition numbers of matrices determine accuracy of linear system solutions, eigenvalue computations, and least squares
• Optimization — Condition number of the Hessian determines convergence speed of gradient descent; preconditioning reduces effective condition number
• ODE Solvers: Stability regions determine maximum step size; stiff systems require implicit solvers