Condition Numbers and Lost Digits

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve, hilbert, qr

def digits_correct(x_computed, x_true):
    """Number of correct decimal digits in x_computed."""
    rel_err = np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true)
    if rel_err == 0:
        return 16.0  # full double precision
    return -np.log10(rel_err)


def make_matrix(n, kappa):
    """Construct an n x n matrix with condition number exactly kappa.
    
    Uses A = U @ diag(sigma) @ V.T with random orthogonal U, V
    and singular values linearly spaced from 1 to 1/kappa.
    """
    U, _ = qr(np.random.randn(n, n))
    V, _ = qr(np.random.randn(n, n))
    sigma = np.linspace(1, 1.0 / kappa, n)
    return U @ np.diag(sigma) @ V.T

np.random.seed(42)
n = 50

kappas_ctrl = [10**k for k in range(0, 16, 2)]
digits_actual_ctrl = []
digits_pred_ctrl = []
rel_errors_ctrl = []

print(f"{'kappa':>12} | {'Predicted':>10} | {'Actual':>10} | {'Rel Error':>12}")
print("-" * 55)

for kappa in kappas_ctrl:
    A = make_matrix(n, kappa)
    x_true = np.random.randn(n)
    b = A @ x_true

    x_computed = solve(A, b)
    rel_err = np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true)
    d_actual = digits_correct(x_computed, x_true)
    d_pred = max(0, 16.0 - np.log10(kappa))

    digits_actual_ctrl.append(d_actual)
    digits_pred_ctrl.append(d_pred)
    rel_errors_ctrl.append(rel_err)

    print(f"{kappa:12.1e} | {d_pred:10.1f} | {d_actual:10.1f} | {rel_err:12.4e}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

log_kappas = [np.log10(k) for k in kappas_ctrl]

# Left: predicted vs actual digits
axes[0].plot(log_kappas, digits_pred_ctrl, 's--', label='Predicted (upper bound)', color='C1')
axes[0].plot(log_kappas, digits_actual_ctrl, 'o-', label='Actual', color='C0')
axes[0].set_xlabel(r'$\log_{10} \kappa(A)$')
axes[0].set_ylabel('Correct digits')
axes[0].set_title('Digits of accuracy vs condition number')
axes[0].set_ylim(-1, 17)
axes[0].axhline(0, color='r', ls='--', alpha=0.5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Right: relative error vs kappa
axes[1].loglog(kappas_ctrl, rel_errors_ctrl, 'o-', label='True relative error', color='C0')
eps = np.finfo(float).eps
kk = np.array(kappas_ctrl, dtype=float)
axes[1].loglog(kk, kk * eps, 's--', label=r'$\kappa \cdot \varepsilon_{\mathrm{mach}}$', color='C1')
axes[1].set_xlabel(r'$\kappa(A)$')
axes[1].set_ylabel('Relative error')
axes[1].set_title('Relative error vs condition number')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ns = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
kappas_hilbert = []
digits_actual_hilbert = []
digits_pred_hilbert = []
rel_errors_hilbert = []

print(f"{'n':>4} | {'kappa(H)':>12} | {'Predicted':>10} | {'Actual':>10} | {'Rel Error':>12}")
print("-" * 65)

for n in ns:
    H = hilbert(n)
    x_true = np.ones(n)
    b = H @ x_true

    x_computed = solve(H, b)
    kappa = np.linalg.cond(H)
    rel_err = np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true)
    d_actual = digits_correct(x_computed, x_true)
    d_pred = max(0, 16.0 - np.log10(kappa))

    kappas_hilbert.append(kappa)
    digits_actual_hilbert.append(d_actual)
    digits_pred_hilbert.append(d_pred)
    rel_errors_hilbert.append(rel_err)

    print(f"{n:4d} | {kappa:12.2e} | {d_pred:10.1f} | {d_actual:10.1f} | {rel_err:12.4e}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Left: condition number growth
axes[0].semilogy(ns, kappas_hilbert, 'o-')
axes[0].set_xlabel('Matrix size $n$')
axes[0].set_ylabel(r'$\kappa_2(H_n)$')
axes[0].set_title('Condition number of Hilbert matrix')
axes[0].grid(True, alpha=0.3)
axes[0].axhline(1 / np.finfo(float).eps, color='r', ls='--', alpha=0.5,
                label=r'$1/\varepsilon_{\mathrm{mach}}$')
axes[0].legend()

# Right: predicted vs actual digits
axes[1].plot(ns, digits_pred_hilbert, 's--', label='Predicted (upper bound)', color='C1')
axes[1].plot(ns, digits_actual_hilbert, 'o-', label='Actual', color='C0')
axes[1].set_xlabel('Matrix size $n$')
axes[1].set_ylabel('Correct digits')
axes[1].set_title('Digits of accuracy: Hilbert matrix')
axes[1].set_ylim(-1, 17)
axes[1].axhline(0, color='r', ls='--', alpha=0.5, label='No correct digits')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ns_vander = list(range(5, 25))
kappas_vander = []
digits_actual_vander = []
digits_pred_vander = []
rel_errors_vander = []

print(f"{'n':>4} | {'kappa(V)':>12} | {'Predicted':>10} | {'Actual':>10} | {'Rel Error':>12}")
print("-" * 65)

for n in ns_vander:
    # Equally spaced nodes on [0, 1]
    x_nodes = np.linspace(0, 1, n)
    V = np.vander(x_nodes, increasing=True)

    # True coefficients
    c_true = np.ones(n)
    b = V @ c_true

    c_computed = solve(V, b)
    kappa = np.linalg.cond(V)
    rel_err = np.linalg.norm(c_computed - c_true) / np.linalg.norm(c_true)
    d_actual = digits_correct(c_computed, c_true)
    d_pred = max(0, 16.0 - np.log10(kappa))

    kappas_vander.append(kappa)
    digits_actual_vander.append(d_actual)
    digits_pred_vander.append(d_pred)
    rel_errors_vander.append(rel_err)

    if n <= 12 or n % 4 == 0:
        print(f"{n:4d} | {kappa:12.2e} | {d_pred:10.1f} | {d_actual:10.1f} | {rel_err:12.4e}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Left: condition number growth
axes[0].semilogy(ns_vander, kappas_vander, 'o-')
axes[0].set_xlabel('Number of nodes $n$')
axes[0].set_ylabel(r'$\kappa_2(V)$')
axes[0].set_title('Condition number of Vandermonde matrix')
axes[0].grid(True, alpha=0.3)
axes[0].axhline(1 / np.finfo(float).eps, color='r', ls='--', alpha=0.5,
                label=r'$1/\varepsilon_{\mathrm{mach}}$')
axes[0].legend()

# Right: predicted vs actual digits
axes[1].plot(ns_vander, digits_pred_vander, 's--', label='Predicted (upper bound)', color='C1')
axes[1].plot(ns_vander, digits_actual_vander, 'o-', label='Actual', color='C0')
axes[1].set_xlabel('Number of nodes $n$')
axes[1].set_ylabel('Correct digits')
axes[1].set_title('Digits of accuracy: Vandermonde matrix')
axes[1].set_ylim(-1, 17)
axes[1].axhline(0, color='r', ls='--', alpha=0.5, label='No correct digits')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 5.5))

# Plot all three families
ax.loglog(kappas_ctrl, rel_errors_ctrl, 'o', ms=8, label='Controlled (SVD)', color='C0')
ax.loglog(kappas_hilbert, rel_errors_hilbert, 's', ms=8, label='Hilbert', color='C1')
ax.loglog(kappas_vander, rel_errors_vander, '^', ms=7, label='Vandermonde (uniform)', color='C2')

# Predicted bound
eps = np.finfo(float).eps
kk = np.logspace(0, 18, 200)
ax.loglog(kk, kk * eps, 'k--', lw=1.5, alpha=0.7, label=r'$\kappa \cdot \varepsilon_{\mathrm{mach}}$')

# Machine epsilon floor
ax.axhline(eps, color='gray', ls=':', alpha=0.5, label=r'$\varepsilon_{\mathrm{mach}}$')

# Numerically singular threshold
ax.axvline(1 / eps, color='r', ls=':', alpha=0.5, label=r'$\kappa = 1/\varepsilon_{\mathrm{mach}}$')

ax.set_xlabel(r'$\kappa(A)$', fontsize=12)
ax.set_ylabel('Relative error', fontsize=12)
ax.set_title('Backward stability in action: true error vs condition number', fontsize=12)
ax.set_xlim(1e-1, 1e19)
ax.set_ylim(1e-17, 1e2)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()