Chapter 039 — Numerical Stability

# --- Visualization: Stable vs unstable recurrences ---
# Classic example: computing I_n = integral_0^1 x^n * e^x dx
# Exact recurrence (integration by parts): I_n = e - n * I_{n-1}
# This is UNSTABLE in the forward direction (n growing).
# Stable direction: compute backwards from a high starting estimate.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

e = np.e
N_TERMS = 30

# True values via scipy-free numerical integration (Gaussian quadrature estimate)
# We use the stable backward recurrence as the 'true' reference
def stable_backward(n_max, start_estimate=0.5):
    """Compute I_0..I_{n_max} using backward stable recurrence."""
    # Start from I_{n_max} = start_estimate, recurse backward: I_{n-1} = (e - I_n) / n
    I = np.zeros(n_max + 1)
    I[n_max] = start_estimate  # rough guess — error will wash out going backward
    for k in range(n_max, 0, -1):
        I[k-1] = (e - I[k]) / k
    return I

def unstable_forward(n_max):
    """Compute I_0..I_{n_max} using forward recurrence: I_n = e - n * I_{n-1}."""
    I = np.zeros(n_max + 1)
    I[0] = e - 1  # exact: I_0 = integral_0^1 e^x dx = e - 1
    for n in range(1, n_max + 1):
        I[n] = e - n * I[n-1]
    return I

I_stable = stable_backward(N_TERMS, start_estimate=1 / (N_TERMS + 1))  # rough guess
I_unstable = unstable_forward(N_TERMS)

ns = np.arange(N_TERMS + 1)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(ns, I_stable, 'o-', color='steelblue', markersize=4, label='Stable (backward)')
axes[0].plot(ns, I_unstable, 's-', color='crimson', markersize=4, label='Unstable (forward)', alpha=0.8)
axes[0].set_ylim(-2, 3)
axes[0].set_xlabel('n', fontsize=11)
axes[0].set_ylabel('Iₙ', fontsize=11)
axes[0].set_title('I_n = ∫₀¹ xⁿeˣ dx: stable vs unstable recurrence', fontsize=10)
axes[0].legend()
axes[0].axhline(0, color='gray', linewidth=0.5)

# Error growth in the forward direction
# Relative error: compare to stable solution as reference
valid = I_stable != 0
rel_err = np.where(valid, np.abs(I_unstable - I_stable) / np.abs(I_stable), np.nan)
axes[1].semilogy(ns[1:], rel_err[1:], 'o-', color='crimson', markersize=4)
axes[1].set_xlabel('n', fontsize=11)
axes[1].set_ylabel('Relative error |forward - stable| / |stable|', fontsize=11)
axes[1].set_title('Error grows exponentially in the unstable forward direction', fontsize=10)

plt.tight_layout()
plt.show()

# --- Visualization: Horner's method vs naive polynomial evaluation ---
# Horner's method evaluates p(x) = a_0 + a_1*x + ... + a_n*x^n
# using n multiplications and n additions — and with better numerical stability.
# Naive method computes each power x^k separately, accumulating rounding errors.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

# Wilkinson's polynomial: (x-1)(x-2)...(x-20), roots at 1,2,...,20
# Expanded coefficients (computed exactly)
roots = list(range(1, 21))
coeffs = np.poly(roots)  # numpy gives coefficients from highest to lowest degree

def horner(coeffs_high_to_low, x):
    """Evaluate polynomial using Horner's method."""
    result = coeffs_high_to_low[0]
    for c in coeffs_high_to_low[1:]:
        result = result * x + c
    return result

def naive_eval(coeffs_high_to_low, x):
    """Evaluate polynomial naively: compute each power separately."""
    n = len(coeffs_high_to_low) - 1
    result = 0.0
    for i, c in enumerate(coeffs_high_to_low):
        result += c * x**(n - i)
    return result

# Evaluate near a root (x=10.5 should give a clear positive value;
# at x=10, should give 0 but floating-point makes it noisy)
x_vals = np.linspace(0.5, 20.5, 2000)
y_horner = np.array([horner(coeffs, x) for x in x_vals])
y_naive  = np.array([naive_eval(coeffs, x) for x in x_vals])

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Full view
axes[0].semilogy(x_vals, np.abs(y_horner) + 1e-10, color='steelblue', lw=1, label="Horner's method")
axes[0].semilogy(x_vals, np.abs(y_naive) + 1e-10, color='crimson', lw=1, alpha=0.7, label='Naive evaluation')
axes[0].set_xlabel('x', fontsize=11)
axes[0].set_ylabel('|p(x)|', fontsize=11)
axes[0].set_title("Wilkinson's polynomial: |p(x)| near roots", fontsize=10)
axes[0].legend()

# Zoom near a root
x_zoom = np.linspace(14.9, 15.1, 500)
axes[1].plot(x_zoom, [horner(coeffs, x) for x in x_zoom], color='steelblue', lw=2, label="Horner's")
axes[1].plot(x_zoom, [naive_eval(coeffs, x) for x in x_zoom], color='crimson', lw=2, alpha=0.8, label='Naive')
axes[1].axhline(0, color='gray', linewidth=0.8)
axes[1].set_xlabel('x', fontsize=11)
axes[1].set_ylabel('p(x)', fontsize=11)
axes[1].set_title('Zoom near root x=15: noise from naive method', fontsize=10)
axes[1].legend()

plt.tight_layout()
plt.show()

# --- Implementation: Numerical Stability Toolkit ---

import numpy as np

def horner_eval(coeffs_low_to_high, x):
    """
    Evaluate polynomial using Horner's method.

    Args:
        coeffs_low_to_high: coefficients [a_0, a_1, ..., a_n]
                            where polynomial = a_0 + a_1*x + ... + a_n*x^n
        x: evaluation point (scalar or array)

    Returns:
        p(x) computed with O(n) operations and O(n*eps_m) backward error
    """
    result = coeffs_low_to_high[-1]  # start from highest-degree coefficient
    for c in reversed(coeffs_low_to_high[:-1]):
        result = result * x + c
    return result


def condition_number_scalar(f, df, x):
    """
    Compute the relative condition number of scalar function f at point x.

    Args:
        f: callable, the function
        df: callable, the derivative f'
        x: evaluation point

    Returns:
        kappa = |x * f'(x) / f(x)|
    """
    fx = f(x)
    if abs(fx) < 1e-300:
        return np.inf  # near a root: condition number diverges
    return abs(x * df(x) / fx)


def measure_algorithm_error(algorithm, exact, x_values):
    """
    Measure relative error of an algorithm over a range of inputs.

    Args:
        algorithm: callable(x) -> approximate result
        exact: callable(x) -> exact result
        x_values: array of test points

    Returns:
        array of relative errors
    """
    errors = []
    for x in x_values:
        approx = algorithm(x)
        true = exact(x)
        if abs(true) > 1e-300:
            errors.append(abs(approx - true) / abs(true))
        else:
            errors.append(abs(approx - true))
    return np.array(errors)


# --- Validation ---
# Horner vs naive for p(x) = 1 + x + x^2 + x^3 + x^4
coeffs = [1, 1, 1, 1, 1]  # a_0, ..., a_4 (low to high)
x = 2.0
horner_result = horner_eval(coeffs, x)
# Geometric sum: (x^5 - 1)/(x - 1) = 31 for x=2
print(f"p(2) = 1+2+4+8+16 = {horner_result}  (expected: 31)")

# Condition number of f(x) = x^2 - 1 at x = 1 (near root)
f  = lambda x: x**2 - 1
df = lambda x: 2*x
for xval in [0.5, 1.0, 1.001, 2.0, 10.0]:
    k = condition_number_scalar(f, df, xval)
    print(f"  kappa at x={xval:.3f}: {k:.4f}  {'← ILL-CONDITIONED near root' if k > 100 else ''}")

# --- Experiment 1: Error growth in the unstable integral recurrence ---
# Hypothesis: starting from I_0 = e - 1 with a tiny perturbation,
# the error grows as n! through the forward recurrence.
# Try changing: PERTURBATION to see how even machine epsilon causes divergence.

import numpy as np

e = np.e
N_TERMS = 25    # <-- modify this
PERTURBATION = 1e-15  # <-- modify this (try 0, 1e-15, 1e-10)

def forward_recurrence(n_max, I0_error=0.0):
    """I_n = e - n * I_{n-1}, starting from I_0 = (e-1) + error."""
    I = np.zeros(n_max + 1)
    I[0] = e - 1 + I0_error
    for n in range(1, n_max + 1):
        I[n] = e - n * I[n-1]
    return I

I_exact = forward_recurrence(N_TERMS, I0_error=0.0)
I_perturbed = forward_recurrence(N_TERMS, I0_error=PERTURBATION)

errors = np.abs(I_perturbed - I_exact)

# Theoretical prediction: error_n = PERTURBATION * n!
import math
predicted = np.array([abs(PERTURBATION) * math.factorial(n) for n in range(N_TERMS + 1)])

import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

fig, ax = plt.subplots(figsize=(9, 5))
ns = np.arange(N_TERMS + 1)
ax.semilogy(ns, errors + 1e-20, 'o-', color='crimson', markersize=5, label='Actual error')
ax.semilogy(ns, predicted + 1e-20, '--', color='steelblue', lw=2, label='Predicted: ε₀ × n!')
ax.set_xlabel('n', fontsize=11)
ax.set_ylabel('|error|', fontsize=11)
ax.set_title(f'Error growth: perturbation = {PERTURBATION:.0e}, amplified by n!', fontsize=11)
ax.legend()
plt.tight_layout()
plt.show()

# --- Experiment 2: Stable backward recurrence — error suppression ---
# Hypothesis: starting from a WRONG estimate at high n and recursing backward,
# the error shrinks by 1/n at each step — total error is O(eps_m).
# Try changing: START_ERROR to see how robustly backward recurrence suppresses it.

import numpy as np

e = np.e
N_MAX = 30         # <-- modify this
START_ERROR = 1.0  # <-- modify this (try 0.1, 1.0, 10.0 — the backward method absorbs all of these)

def backward_recurrence(n_max, start_estimate_error=0.0):
    """I_{n-1} = (e - I_n) / n, starting from rough estimate at n_max."""
    I_true_at_n = 1 / (n_max + 1)  # rough but reasonable estimate
    I = np.zeros(n_max + 1)
    I[n_max] = I_true_at_n + start_estimate_error  # intentionally wrong
    for k in range(n_max, 0, -1):
        I[k-1] = (e - I[k]) / k
    return I

I_exact   = backward_recurrence(N_MAX, start_estimate_error=0.0)
I_wrong   = backward_recurrence(N_MAX, start_estimate_error=START_ERROR)

errors = np.abs(I_wrong - I_exact)

print(f"Starting error at n={N_MAX}: {START_ERROR}")
print(f"Resulting error at n=0: {errors[0]:.4e}")
print(f"Error suppression factor: {START_ERROR / (errors[0] + 1e-300):.2e}")
print("\nErrors across n:")
for n in [N_MAX, N_MAX-5, N_MAX//2, 5, 0]:
    print(f"  n={n:2d}: error = {errors[n]:.4e}")

# --- Experiment 3: Condition number vs algorithmic error ---
# Hypothesis: for well-conditioned problems, both naive and stable algorithms agree.
# For ill-conditioned problems, unstable algorithms amplify error far beyond κ * eps_m.
# Try changing: POWER to control how ill-conditioned the problem is.

import numpy as np

# Evaluate f(x) = x^POWER - 1 near its root x=1
# Condition number near x=1: kappa = |x * f'(x) / f(x)| = |x * POWER * x^(POWER-1) / (x^POWER - 1)|
# As x -> 1, kappa -> POWER / |x-1|  (grows as we approach the root)

POWER = 20   # <-- modify this

def f(x): return x**POWER - 1
def df(x): return POWER * x**(POWER-1)

x_test_points = 1.0 + np.logspace(-15, -1, 100)

kappas = np.array([abs(x * df(x) / f(x)) if abs(f(x)) > 1e-300 else np.inf for x in x_test_points])
floor_errors = kappas * np.finfo(float).eps  # minimum achievable relative error

import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

fig, ax = plt.subplots(figsize=(9, 5))
ax.loglog(x_test_points - 1, kappas, color='darkorange', lw=2, label='Condition number κ')
ax.loglog(x_test_points - 1, floor_errors, '--', color='steelblue', lw=2, label='Error floor = κ × ε_m')
ax.axhline(np.finfo(float).eps, color='green', linestyle=':', lw=1.5, label='ε_m (perfect algorithm)')
ax.set_xlabel('x − 1  (distance from root)', fontsize=11)
ax.set_ylabel('Relative quantity', fontsize=11)
ax.set_title(f'Condition number of x^{POWER}−1 near root x=1', fontsize=11)
ax.legend(fontsize=10)
plt.tight_layout()
plt.show()

# --- Mini Project: Comprehensive Summation Benchmark ---
# Problem:
#   Summation of n floating-point numbers is a fundamental operation.
#   Different algorithms trade accuracy against speed.
#   Build a benchmark that compares 5 methods across multiple distributions
#   and produces a 'stability vs speed' report.
#
# Methods to compare:
#   1. Python built-in sum() — left-to-right
#   2. numpy sum (default) — may use SIMD/chunked but not compensated
#   3. Kahan compensated sum
#   4. Sorted sum (ascending order)
#   5. math.fsum — exact (uses extended precision internally)

import numpy as np
import math
import time
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

def kahan_sum(values):
    s, c = 0.0, 0.0
    for x in values:
        y = x - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

# Test scenarios
N = 100_000  # <-- modify this
rng = np.random.default_rng(42)

scenarios = {
    'Uniform [0,1]':         rng.uniform(0, 1, N),
    'Mixed large+small':     np.concatenate([rng.uniform(1e6, 2e6, N//2), rng.uniform(0, 1, N//2)]),
    'Alternating sign':      np.where(np.arange(N) % 2 == 0, rng.uniform(1, 2, N), -rng.uniform(1, 2, N)),
    'Exponential':           rng.exponential(1.0, N),
}

methods = [
    ('Built-in sum',   lambda arr: sum(arr)),
    ('numpy sum',      lambda arr: float(np.sum(arr))),
    ('Kahan sum',      kahan_sum),
    ('Sorted sum',     lambda arr: float(np.sum(np.sort(arr)))),
    ('math.fsum',      math.fsum),
]

print(f"{'Scenario':<25} | {'Method':<15} | {'Rel. Error':>12} | {'Time (ms)':>10}")
print('-' * 72)

for scenario_name, arr in scenarios.items():
    exact = math.fsum(arr)  # best available reference
    for method_name, fn in methods[:-1]:  # skip fsum in comparison (it's the reference)
        REPS = 3
        t0 = time.perf_counter()
        result = None
        for _ in range(REPS):
            result = fn(arr)
        t1 = time.perf_counter()
        rel_err = abs(result - exact) / (abs(exact) + 1e-300)
        time_ms = (t1 - t0) / REPS * 1000
        print(f"{scenario_name:<25} | {method_name:<15} | {rel_err:>12.4e} | {time_ms:>10.3f}")
    print()