Numerical Precision Challenges in ATen Operations

This document covers precision loss concerns in PyTorch ATen operations and best practices for validating numerically sensitive operations.

Table of Contents

1. Overview of Precision Loss
2. Common Numerically Sensitive Operations
3. Understanding ULP (Units in Last Place)
4. Validation Methods
5. Best Practices for Testing
6. PyTorch-Specific Considerations

Overview of Precision Loss

Floating-point arithmetic has inherent limitations due to finite precision. Precision loss becomes particularly problematic in certain mathematical operations:

Catastrophic Cancellation

When subtracting two nearly-equal numbers, most significant digits cancel out, leaving only the less accurate trailing digits:

# Example: computing log(1 + x) for small x
x = 1e-15
naive = math.log(1 + x)      # 1 + x rounds to 1.0, result is 0.0
correct = math.log1p(x)      # ≈ 1e-15 (preserves precision)

Absorption

When adding numbers of vastly different magnitudes, the smaller value gets “absorbed”:

big = 1e16
small = 1.0
result = big + small - big   # Should be 1.0, but equals 0.0

Overflow/Underflow

Intermediate results exceed representable range:

# Computing sqrt(x² + y²) naively
x, y = 1e200, 1e200
naive = math.sqrt(x**2 + y**2)   # Overflow! x² = inf
correct = math.hypot(x, y)        # ≈ 1.414e200 (uses scaling internally)

Common Numerically Sensitive Operations

1. Companion Functions (Avoid Precision Loss)

2. Log-Space Operations (Avoid Overflow)

3. Normalization Operations

4. Special Mathematical Functions

Many special functions have regions where naive implementations fail:

• Bessel functions near zeros
• Gamma/lgamma for large arguments
• Error function (erf/erfc) at extremes
• Trigonometric functions for large arguments (argument reduction)

Understanding ULP (Units in Last Place)

What is ULP?

ULP measures the spacing between adjacent floating-point numbers. The ULP error is how many representable floats exist between the computed and true values.

For float32 near 1.0:  1 ULP ≈ 1.19e-7
For float32 near 1e10: 1 ULP ≈ 1024
For float64 near 1.0:  1 ULP ≈ 2.22e-16

Why ULP Can Be Misleading

Problem 1: ULP varies with magnitude

# Near 1.0, 1 ULP ≈ 1e-7 (float32)
# Near 0.0, 1 ULP ≈ 1e-45 (denormalized range)
# An "acceptable" 1 ULP error means very different things!

Problem 2: Numerically unstable operations can have huge ULP errors

# log1p(1e-10) with naive log(1+x):
# true value:  1.0e-10
# computed:    0.0 (1+1e-10 rounds to 1.0 in float32)
# ULP error:   ~670 million ULPs!

Problem 3: ULP doesn’t account for mathematical conditioning Some functions are inherently ill-conditioned (small input changes → large output changes). High ULP error may be unavoidable for ill-conditioned inputs.

When ULP is Useful

• Comparing implementations of the same algorithm
• Verifying against high-precision reference (mpfr)
• Testing well-conditioned operations
• Regression testing (has ULP error increased?)

Validation Methods

1. Absolute + Relative Tolerance (Recommended Default)

The standard approach used by PyTorch and NumPy:

def allclose(computed, expected, rtol=1e-5, atol=1e-8):
    return |computed - expected| <= atol + rtol * |expected|

Advantages:

• Handles values near zero (atol dominates)
• Handles large values (rtol dominates)
• Intuitive to tune

PyTorch defaults by dtype:

2. ULP-Based Testing

import numpy as np

def ulp_diff(a, b):
    """Count ULPs between two floats."""
    # Convert to integer representation
    a_int = np.float32(a).view(np.int32)
    b_int = np.float32(b).view(np.int32)
    
    # Handle sign differences
    if (a_int < 0) != (b_int < 0):
        return abs(a_int) + abs(b_int)
    return abs(a_int - b_int)

# Usage: assert ulp_diff(computed, expected) <= max_ulp

When to use:

• Testing against arbitrary-precision reference
• Verifying mathematical library functions
• Typical acceptable ranges: 0.5-4 ULP for basic ops

3. High-Precision Reference Testing

Use arbitrary-precision libraries to compute “ground truth”:

from mpmath import mp
mp.dps = 50  # 50 decimal places

def test_log1p_precision():
    test_values = [1e-15, 1e-10, 1e-5, 0.1, 0.5]
    for x in test_values:
        expected = float(mp.log1p(x))
        computed = torch.tensor(x).log1p().item()
        
        rel_error = abs(computed - expected) / abs(expected)
        assert rel_error < 1e-6, f"log1p({x}): rel_error={rel_error}"

4. Domain-Specific Testing

Test in regions where precision loss is expected:

def test_log1p_small_values():
    """log1p should be accurate for small x where log(1+x) fails."""
    # Critical region: |x| < epsilon where 1+x rounds to 1
    small_x = torch.tensor([1e-7, 1e-8, 1e-10, 1e-15], dtype=torch.float32)
    
    # In this region, log1p(x) ≈ x (Taylor series)
    result = torch.log1p(small_x)
    
    # Check relative error against Taylor approximation
    # log1p(x) = x - x²/2 + x³/3 - ...
    taylor_approx = small_x - small_x**2/2
    torch.testing.assert_close(result, taylor_approx, rtol=1e-5, atol=0)

5. Gradient Checking for Autograd

Numerical gradients can also suffer from precision issues:

def test_gradient_log1p():
    x = torch.tensor([1e-10], requires_grad=True, dtype=torch.float64)
    
    # Analytical gradient: d/dx log1p(x) = 1/(1+x)
    y = torch.log1p(x)
    y.backward()
    analytical = x.grad.item()
    
    # For small x, gradient ≈ 1
    expected = 1.0 / (1.0 + x.item())
    
    assert abs(analytical - expected) / expected < 1e-10

6. Special Value Testing

Always test edge cases:

def test_special_values():
    special = torch.tensor([
        0.0,           # Exact: log1p(0) = 0
        -0.0,          # Signed zero
        float('inf'),  # log1p(inf) = inf
        float('-inf'), # log1p(-inf) = nan (domain error for x < -1)
        float('nan'),  # log1p(nan) = nan
        -1.0,          # log1p(-1) = -inf
        -1.0 - 1e-10,  # Slightly below domain: should be nan
    ])
    
    result = torch.log1p(special)
    
    assert result[0] == 0.0
    assert result[1] == 0.0  # or -0.0
    assert result[2] == float('inf')
    assert torch.isnan(result[3])
    assert torch.isnan(result[4])
    assert result[5] == float('-inf')
    assert torch.isnan(result[6])

7. Monotonicity and Consistency Testing

For monotonic functions, verify ordering is preserved:

def test_monotonicity():
    x = torch.linspace(-0.99, 100, 10000)
    y = torch.log1p(x)
    
    # log1p is strictly increasing
    diffs = y[1:] - y[:-1]
    assert (diffs > 0).all(), "log1p should be strictly increasing"

Best Practices for Testing

1. Choose Appropriate Tolerances by Region

def get_tolerances_for_log1p(x):
    """Different tolerances for different input regions."""
    x_abs = abs(x)
    
    if x_abs < 1e-7:
        # Near zero: result ≈ x, use relative tolerance to x
        return {'rtol': 1e-6, 'atol': 0}
    elif x_abs < 1:
        # Normal range: standard tolerances
        return {'rtol': 1e-6, 'atol': 1e-12}
    else:
        # Large x: log1p(x) ≈ log(x), relative tolerance fine
        return {'rtol': 1e-6, 'atol': 0}

2. Test Multiple Dtypes

@pytest.mark.parametrize("dtype", [
    torch.float16, torch.bfloat16, torch.float32, torch.float64
])
def test_log1p_dtype(dtype):
    x = torch.tensor([1e-5, 0.1, 1.0, 100.0], dtype=dtype)
    result = torch.log1p(x)
    
    # Use dtype-appropriate tolerances
    tols = {
        torch.float16: {'rtol': 1e-3, 'atol': 1e-4},
        torch.bfloat16: {'rtol': 2e-2, 'atol': 1e-3},
        torch.float32: {'rtol': 1e-5, 'atol': 1e-6},
        torch.float64: {'rtol': 1e-10, 'atol': 1e-12},
    }[dtype]
    
    expected = torch.tensor([float(mp.log1p(v)) for v in x.tolist()], dtype=dtype)
    torch.testing.assert_close(result, expected, **tols)

3. Comparative Testing Against NumPy/SciPy

def test_against_numpy():
    x = np.random.uniform(-0.99, 100, 10000).astype(np.float32)
    
    numpy_result = np.log1p(x)
    torch_result = torch.log1p(torch.from_numpy(x)).numpy()
    
    np.testing.assert_allclose(torch_result, numpy_result, rtol=1e-6, atol=1e-6)

4. Stress Testing for Numerical Stability

def test_logsumexp_stability():
    """logsumexp should not overflow even with extreme values."""
    # Values that would overflow with naive exp()
    x = torch.tensor([1000.0, 1000.0, 1000.0])
    result = torch.logsumexp(x, dim=0)
    
    # log(3 * exp(1000)) = 1000 + log(3)
    expected = 1000.0 + math.log(3)
    assert abs(result.item() - expected) < 1e-4
    
    # Very negative values should not underflow to -inf
    x_neg = torch.tensor([-1000.0, -1000.0, -1000.0])
    result_neg = torch.logsumexp(x_neg, dim=0)
    expected_neg = -1000.0 + math.log(3)
    assert abs(result_neg.item() - expected_neg) < 1e-4

5. Test Across the Full Domain

def test_full_domain():
    """Test across entire valid domain with logarithmic spacing."""
    # Positive values: from smallest denormal to largest finite
    pos_exp = torch.linspace(-45, 38, 1000)  # float32 range
    pos_values = 10 ** pos_exp
    
    # Negative values: from -1+ε to 0
    neg_values = -torch.logspace(-7, -1, 100)
    
    all_values = torch.cat([neg_values, pos_values])
    
    # Should not produce NaN (except for x < -1)
    result = torch.log1p(all_values)
    assert not torch.isnan(result).any()

PyTorch-Specific Considerations

1. CUDA Fast-Math Intrinsics

From aten/src/ATen/NumericUtils.h:

template <typename T>
C10_HOST_DEVICE inline T log1p(T x) {
#if defined(__CUDA_ARCH__) || defined(__HIP_ARCH__)
  // NOTE: There is no __log1pf so unfortunately we lose precision.
  return __logf(1.0f + x);  // Fast but imprecise!
#else
  return ::log1p(x);  // Precise
#endif
}

Implication: CUDA float32 log1p loses precision for small values. Tests may need looser tolerances for CUDA.

2. Device-Specific Testing

@pytest.mark.parametrize("device", ['cpu', 'cuda'])
def test_log1p_device(device):
    if device == 'cuda' and not torch.cuda.is_available():
        pytest.skip("CUDA not available")
    
    x = torch.tensor([1e-10], device=device, dtype=torch.float32)
    result = torch.log1p(x)
    
    # CUDA needs looser tolerance due to __logf approximation
    rtol = 1e-3 if device == 'cuda' else 1e-6
    expected = torch.tensor([1e-10], device=device)
    torch.testing.assert_close(result, expected, rtol=rtol, atol=0)

3. Half-Precision (float16/bfloat16)

These formats have very limited precision:

• float16: 10-bit mantissa, ~3 decimal digits
• bfloat16: 7-bit mantissa, ~2 decimal digits

def test_half_precision():
    # float16 can't represent 1 + 1e-5 distinctly from 1
    x = torch.tensor([1e-5], dtype=torch.float16)
    
    # This will have significant error
    result = torch.log1p(x)
    
    # Use very loose tolerances
    torch.testing.assert_close(
        result, 
        torch.tensor([1e-5], dtype=torch.float16),
        rtol=0.1, atol=1e-4  # 10% relative tolerance!
    )

4. Using torch.testing.assert_close

PyTorch’s recommended testing function:

import torch.testing

# Basic usage
torch.testing.assert_close(actual, expected, rtol=1e-5, atol=1e-8)

# With custom message
torch.testing.assert_close(
    actual, expected, 
    rtol=1e-5, atol=1e-8,
    msg=f"Failed for input {input_val}"
)

# Check NaN equality (by default, NaN == NaN is True)
torch.testing.assert_close(
    torch.tensor([float('nan')]),
    torch.tensor([float('nan')]),
    equal_nan=True  # Default
)

Summary: Validation Strategy for Numerically Sensitive Ops

1. Use combined atol + rtol tolerances as the primary method
2. Adjust tolerances by:
   • Input region (near zero, normal, large)
   • Data type (float16 needs 10-100x looser)
   • Device (CUDA may have reduced precision)
3. Test comprehensively:
   • Special values (0, ±inf, nan, domain boundaries)
   • Critical regions (where precision loss would occur with naive impl)
   • Full domain with logarithmic spacing
4. Use high-precision reference (mpfr) for ground truth when needed
5. ULP testing is useful for comparing implementations, but be aware it varies with magnitude
6. Verify mathematical properties (monotonicity, symmetry, identities)
7. Test gradients separately since they can have different precision characteristics

References

• What Every Computer Scientist Should Know About Floating-Point Arithmetic
• IEEE 754 Floating Point Standard
• PyTorch Testing: torch.testing.assert_close
• NumPy Testing: numpy.testing.assert_allclose
• Arbitrary precision: mpmath library