tolerances

Numerical tolerances derived from backward error analysis.

This module provides principled tolerance functions for numerical testing, based on standard backward error bounds from numerical linear algebra.

Error Bounds (Golub & Van Loan, Matrix Computations, 4th ed.)¶

QR Factorization (§5.2): The computed Q̂, R̂ satisfy A + δA = Q̂R̂ where ‖δA‖ ≤ O(ε)‖A‖

Triangular Solve (§3.1): (R + δR)x̂ = b where ‖δR‖ ≤ O(nε)‖R‖

Forward error: ‖x - x̂‖/‖x‖ ≤ O(κ(A) × ε) Forward error: ‖x - x̂‖/‖x‖ ≤ O(κ(A) × ε)

At OLS solution: X’e = 0 theoretically At OLS solution: X’e = 0 theoretically Computed: ‖X’ê‖ ≤ O(ε × κ(X) × ‖X‖ × ‖y‖ × n)

Design Principles¶

1. Filter pathological cases with assume() at generation time

2. Tolerances derived from theory, not empirical tuning

3. Include safety factors for implementation variations

4. Document error bound source for each function

References¶

.. [golub2013] Golub, Van Loan. Matrix Computations (2013), 4th ed. .. [higham2002] Higham. Accuracy and Stability of Numerical Algorithms (2002).

Attributes:

Functions:

Attributes¶

DEFAULT_SAFETY¶

DEFAULT_SAFETY: float = 10.0

EPS¶

EPS: float = np.finfo(np.float64).eps

MAX_COND¶

MAX_COND: float = 10000000000.0

MAX_COND_INVERSE¶

MAX_COND_INVERSE: float = 100000.0

Functions¶

algorithm_comparison_atol¶

algorithm_comparison_atol(scale: float, cond: float, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for comparing different algorithms.

When comparing results from two different numerical methods (e.g., QR vs SVD), each method has its own rounding pattern. The difference between them is bounded by the sum of their individual errors.

Error bound: O(2 × ε × κ × scale)

Parameters¶

scale : float Characteristic scale of the result (e.g., max(|coef|)). cond : float Condition number of the problem. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for algorithm comparison.

Examples¶

result_qr = qr_solve(X, y) result_svd = svd_solve(X, y) cond = np.linalg.cond(X) scale = max(np.abs(result_qr.coef).max(), 1.0) atol = algorithm_comparison_atol(scale, cond) np.testing.assert_allclose(result_qr.coef, result_svd.coef, atol=atol)

algorithm_comparison_rtol¶

algorithm_comparison_rtol(cond: float, safety: float = DEFAULT_SAFETY) -> float

Relative tolerance for comparing different algorithms.

Parameters¶

cond : float Condition number of the problem. safety : float, default=10.0 Safety factor.

Returns¶

float Relative tolerance for algorithm comparison.

decomposition_atol¶

decomposition_atol(A: ArrayLike, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for decomposition properties.

For verifying properties like Q’Q = I, R is upper triangular, etc.

Based on backward error: ‖δA‖ ≤ O(ε × ‖A‖ × n)

Parameters¶

A : array-like Matrix being decomposed. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for decomposition checks.

fitted_atol¶

fitted_atol(X: ArrayLike, y: ArrayLike, cond: float, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for fitted value comparisons.

For verifying ŷ = Xβ̂ and y = ŷ + e.

Parameters¶

X : array-like Design matrix. y : array-like Response vector. cond : float Condition number of X. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for fitted value checks.

glm_score_atol¶

glm_score_atol(X: ArrayLike, y: ArrayLike, cond: float, weights: ArrayLike | None = None, *, iterative: bool = False, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for GLM score equation checks.

For verifying X’W(y - μ) ≈ 0 at IRLS convergence.

GLMs use iterative fitting (IRLS), which accumulates more numerical error than closed-form OLS. Non-Gaussian families need looser tolerances because the score function is nonlinear in the coefficients.

Error bound: O(ε × κ(X) × ‖X‖ × ‖y‖ × ‖W‖ × n × safety × irls_factor)

Parameters¶

X : array-like Design matrix (n × p). y : array-like Response vector (n,). cond : float Condition number of X. weights : array-like, optional Weight vector or None for unweighted. iterative : bool, default=False If True, uses looser tolerance suitable for iterative algorithms like IRLS (e.g., Poisson, Gamma). Set True for non-Gaussian families. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for score equation checks.

Examples¶

result = fit_glm(X, y, family=‘binomial’) score = X.T @ (y - result.mu) atol = glm_score_atol(X, y, np.linalg.cond(X)) np.testing.assert_allclose(score, 0, atol=atol)

has_full_rank¶

has_full_rank(A: ArrayLike) -> bool

Check if matrix has full column rank.

Parameters¶

A : array-like Matrix to check (n × p where n ≥ p).

Returns¶

bool True if rank(A) == min(n, p).

hat_matrix_atol¶

hat_matrix_atol(X: ArrayLike, cond: float, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for hat matrix property checks.

For verifying H² = H (idempotence) and H’ = H (symmetry).

The hat matrix H = X(X’X)^{-1}X’ involves inverting X’X, which has condition number κ(X)². Additionally, checking H² = H involves matrix multiplication that accumulates error.

Error bound: O(ε × κ(X)² × n²)

Parameters¶

X : array-like Design matrix (n × p). cond : float Condition number of X. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for hat matrix checks.

Notes¶

For well-conditioned tests, use assume(is_well_conditioned(X, MAX_COND_INVERSE)) to filter cases where κ(X) > 1e5 (meaning κ(X’X) > 1e10).

Examples¶

H = X @ np.linalg.inv(X.T @ X) @ X.T cond = np.linalg.cond(X) np.testing.assert_allclose(H @ H, H, atol=hat_matrix_atol(X, cond))

inference_atol¶

inference_atol(coef: ArrayLike, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for inference result comparisons.

For verifying CI midpoints, t-statistics, and other inference quantities that depend on coefficient estimates.

Error bound: O(ε × max(|β|) × safety)

Parameters¶

coef : array-like Coefficient estimates (the scale reference). safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for inference comparisons.

Examples¶

result = infer(coef, vcov, df) ci_midpoint = (result.ci_lower + result.ci_upper) / 2 np.testing.assert_allclose(ci_midpoint, coef, atol=inference_atol(coef))

is_well_conditioned¶

is_well_conditioned(A: ArrayLike, threshold: float = MAX_COND) -> bool

Check if matrix is well-conditioned for stable computation.

Use with hypothesis.assume() to filter pathological test cases.

Parameters¶

A : array-like Matrix to check. threshold : float, default=MAX_COND Maximum acceptable condition number.

Returns¶

bool True if condition number is below threshold.

Examples¶

from hypothesis import assume assume(is_well_conditioned(X))

orthogonality_atol¶

orthogonality_atol(n: int, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for orthogonality checks.

For verifying Q’Q = I or P’P = I where Q/P is an orthonormal matrix.

Based on backward error for orthogonal matrices: O(ε × n).

Parameters¶

n : int Dimension of the orthogonal matrix (number of rows). safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for orthogonality checks.

Examples¶

Q, R = np.linalg.qr(X) n = Q.shape[0] np.testing.assert_allclose(Q.T @ Q, np.eye(Q.shape[1]), atol=orthogonality_atol(n))

residual_atol¶

residual_atol(X: ArrayLike, y: ArrayLike, cond: float, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for residual orthogonality checks.

For verifying X’e ≈ 0 (normal equations).

Based on bound: ‖X’ê‖ ≤ O(ε × κ(X) × ‖X‖ × ‖y‖ × n)

Parameters¶

X : array-like Design matrix (n × p). y : array-like Response vector (n,). cond : float Condition number of X. safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance for residual orthogonality.

Examples¶

result = qr_solve(X, y) score = X.T @ result.residuals np.testing.assert_allclose(score, 0, atol=residual_atol(X, y, cond))

solve_atol¶

solve_atol(scale: float, cond: float, n: int = 1, safety: float = DEFAULT_SAFETY) -> float

Absolute tolerance for linear solve operations.

For near-zero solutions where relative tolerance is meaningless.

Parameters¶

scale : float Characteristic scale of the problem (e.g., max(|X|, |y|)). cond : float Condition number of the system matrix. n : int, default=1 Problem dimension (for accumulated rounding). safety : float, default=10.0 Safety factor.

Returns¶

float Absolute tolerance with minimum floor.

solve_rtol¶

solve_rtol(cond: float, safety: float = DEFAULT_SAFETY) -> float

Relative tolerance for linear solve operations.

Based on forward error bound for QR/SVD solve: ‖x - x̂‖/‖x‖ ≤ O(κ(A) × ε)

Parameters¶

cond : float Condition number of the system matrix. safety : float, default=10.0 Safety factor for implementation variations.

Returns¶

float Relative tolerance for comparing solutions.

Examples¶

cond = 1e6 rtol = solve_rtol(cond) np.testing.assert_allclose(computed, expected, rtol=rtol)