Validating Numerical Precision for Low-Precision Ops

This document describes methods and tools to detect numerical precision issues when implementing mathematical operations for low-precision floating-point types like bfloat16.

Overview

When implementing ops like softplus, log, exp, etc. in bfloat16, precision issues can arise from:

1. Rounding errors - Limited mantissa bits (7 for bf16 vs 23 for fp32)
2. Catastrophic cancellation - Subtracting nearly equal numbers
3. Absorption - Adding small values to large values (e.g., 1 + small → 1)
4. Error amplification - Operations like exp() and log() can amplify errors

Method 1: Sensitivity Analysis with Perturbation Testing

Check if small input changes cause disproportionate output changes.

import torch
import numpy as np

def detect_precision_loss(fn, x, dtype, epsilon=1e-4):
    """
    Check if small input changes cause disproportionate output changes.
    
    Args:
        fn: Function to test
        x: Input value or tensor
        dtype: torch dtype to test (e.g., torch.bfloat16)
        epsilon: Threshold for acceptable relative change
        
    Returns:
        bool: True if precision loss detected
    """
    x_tensor = torch.tensor(x, dtype=dtype)
    y = fn(x_tensor)
    
    # Perturb by 1 ULP (Unit in Last Place)
    x_perturbed = x_tensor + torch.finfo(dtype).eps * x_tensor.abs()
    y_perturbed = fn(x_perturbed)
    
    # If output changes more than expected, precision issue
    relative_change = (y_perturbed - y).abs() / y.abs()
    return relative_change > epsilon

When to use: Quick check for input sensitivity at specific values.

Method 2: Compare Against Higher Precision Reference

The gold standard - compare bf16 results against fp32 or fp64.

def find_precision_boundary(fn, x_range, beta=1.0, threshold=0.01):
    """
    Find where bf16 diverges from fp32 reference.
    
    Args:
        fn: Function to test (takes tensor and beta)
        x_range: Tuple of (min, max) for input sweep
        beta: Parameter for functions like softplus
        threshold: Relative error threshold
        
    Returns:
        List of (x, error) tuples where precision loss occurs
    """
    problem_points = []
    
    for x in np.linspace(x_range[0], x_range[1], 1000):
        bf16_result = fn(torch.tensor(x, dtype=torch.bfloat16), beta).item()
        fp32_result = fn(torch.tensor(x, dtype=torch.float32), beta).item()
        
        if fp32_result != 0:
            rel_error = abs(bf16_result - fp32_result) / abs(fp32_result)
            if rel_error > threshold:
                problem_points.append((x, rel_error))
                
    return problem_points

# Example usage for softplus
def softplus_ref(x, beta):
    return torch.nn.functional.softplus(x, beta=beta)

problems = find_precision_boundary(softplus_ref, (-10, 5), beta=2.0)
for x, err in problems[:10]:
    print(f"x={x:.2f}: error={err:.2%}")

When to use: Systematic sweep to find problem regions.

Method 3: ULP (Unit in Last Place) Analysis

Count how many ULPs apart two values are - more meaningful than relative error for floating point.

def ulp_difference(expected, actual, dtype):
    """
    Count how many ULPs apart two values are.
    
    Args:
        expected: Reference value (tensor or float)
        actual: Computed value (tensor or float)
        dtype: The dtype being tested
        
    Returns:
        Number of ULPs difference
    """
    if not isinstance(expected, torch.Tensor):
        expected = torch.tensor(expected, dtype=dtype)
    if not isinstance(actual, torch.Tensor):
        actual = torch.tensor(actual, dtype=dtype)
        
    eps = torch.finfo(dtype).eps
    ulp = (actual - expected).abs() / (expected.abs() * eps + 1e-30)
    return ulp

def analyze_ulp_distribution(expected, actual, dtype):
    """
    Analyze distribution of ULP differences.
    
    Returns dict with statistics.
    """
    ulps = ulp_difference(expected, actual, dtype)
    
    return {
        'max_ulp': ulps.max().item(),
        'mean_ulp': ulps.mean().item(),
        'median_ulp': ulps.median().item(),
        'pct_within_1_ulp': (ulps <= 1).float().mean().item(),
        'pct_within_2_ulp': (ulps <= 2).float().mean().item(),
        'pct_within_4_ulp': (ulps <= 4).float().mean().item(),
    }

When to use: Understanding the magnitude of errors in floating-point terms.

Method 4: Condition Number Analysis

For a function $f(x)$, the condition number is: $$\kappa(x) = \left| \frac{x \cdot f’(x)}{f(x)} \right|$$

High condition number = sensitive to input perturbations.

def condition_number_softplus(x, beta):
    """
    Estimate condition number for softplus.
    
    softplus(x) = log(1 + exp(beta*x)) / beta
    softplus'(x) = sigmoid(beta*x)
    
    kappa = |x * sigmoid(beta*x) / softplus(x)|
    """
    import numpy as np
    
    beta_x = beta * x
    exp_bx = np.exp(np.clip(beta_x, -500, 500))  # Prevent overflow
    sigmoid = exp_bx / (1 + exp_bx)
    softplus = np.log1p(exp_bx) / beta
    
    if softplus == 0:
        return float('inf')
    return abs(x * sigmoid / softplus)

# Sweep to find high-condition regions
for x in np.linspace(-10, 10, 50):
    kappa = condition_number_softplus(x, beta=2.0)
    if kappa > 10:
        print(f"x={x:.1f}: condition number = {kappa:.1f} (unstable!)")

When to use: Theoretical analysis to identify inherently unstable regions.

Method 5: Error Propagation Analysis

Track error bounds through each operation step.

def error_propagation_softplus(x, beta, dtype=torch.bfloat16):
    """
    Track error accumulation through softplus computation.
    
    Shows where precision is lost in the computation chain.
    """
    eps = torch.finfo(dtype).eps.item()
    
    # Step 1: beta * x
    beta_x = beta * x
    err_mul = abs(beta_x) * eps
    
    # Step 2: exp(beta_x)
    exp_bx = np.exp(beta_x)
    err_exp = exp_bx * (err_mul + eps)  # exp amplifies relative error
    
    # Step 3: 1 + exp_bx (CRITICAL for small exp_bx!)
    sum_val = 1 + exp_bx
    if exp_bx < eps:
        # This is where precision is lost!
        err_sum = eps  # Rounds to 1.0, loses all information about exp_bx
        absorption_warning = True
    else:
        err_sum = err_exp
        absorption_warning = False
    
    # Step 4: log(sum)
    log_val = np.log(sum_val)
    if log_val != 0:
        err_log = err_sum / sum_val + eps * abs(log_val)
        relative_error = err_log / abs(log_val)
    else:
        err_log = float('inf')
        relative_error = float('inf')
    
    # Step 5: divide by beta
    final = log_val / beta
    
    return {
        'x': x,
        'beta_x': beta_x,
        'exp_bx': exp_bx,
        'sum_val': sum_val,
        'log_val': log_val,
        'final': final,
        'estimated_relative_error': relative_error,
        'absorption_warning': absorption_warning,
    }

# Example: analyze the problem region
print("Error propagation analysis for softplus with beta=2.0:")
print("-" * 60)
for x in [-1, -2, -3, -4, -5, -6, -8, -10]:
    result = error_propagation_softplus(x, beta=2.0)
    warning = " ⚠️ ABSORPTION!" if result['absorption_warning'] else ""
    print(f"x={x:3d}: exp(βx)={result['exp_bx']:.2e}, "
          f"err={result['estimated_relative_error']:.2%}{warning}")

When to use: Understanding WHY precision loss occurs at each step.

Method 6: Automated Threshold Finder

Sweep to find where Taylor approximation beats standard formula.

def find_optimal_neg_threshold(beta=1.0, target_error=0.05):
    """
    Find optimal negative threshold for switching to Taylor approximation.
    
    For softplus: when should we use exp(beta*x)/beta instead of log(1+exp(beta*x))/beta?
    
    Args:
        beta: Beta parameter
        target_error: Maximum acceptable relative error
        
    Returns:
        Optimal threshold value
    """
    for neg_thresh in np.arange(-1, -15, -0.5):
        max_error_standard = 0
        max_error_taylor = 0
        
        # Test in the region around this threshold
        for x in np.linspace(neg_thresh - 2, neg_thresh + 2, 200):
            beta_x = beta * x
            
            # Reference (fp64)
            ref = float(np.log1p(np.exp(beta_x))) / beta
            if ref == 0:
                continue
            
            # Standard formula simulated in bf16
            exp_bx = float(torch.tensor(np.exp(beta_x), dtype=torch.bfloat16))
            one_plus = float(torch.tensor(1.0 + exp_bx, dtype=torch.bfloat16))
            log_val = float(torch.tensor(np.log(one_plus), dtype=torch.bfloat16))
            std_bf16 = log_val / beta
            
            # Taylor approximation: log(1+z) ≈ z for small z
            taylor = exp_bx / beta
            
            err_std = abs(std_bf16 - ref) / abs(ref)
            err_taylor = abs(taylor - ref) / abs(ref)
            
            max_error_standard = max(max_error_standard, err_std)
            max_error_taylor = max(max_error_taylor, err_taylor)
        
        # Find where Taylor becomes better than standard
        if max_error_taylor < max_error_standard:
            print(f"neg_thresh={neg_thresh}: "
                  f"standard_err={max_error_standard:.2%}, "
                  f"taylor_err={max_error_taylor:.2%}")
            
            if max_error_taylor < target_error:
                print(f"\n✓ Optimal threshold: {neg_thresh} "
                      f"(Taylor error < {target_error:.0%})")
                return neg_thresh
    
    return None

# Find threshold for different beta values
for beta in [1.0, 2.0, 4.0]:
    print(f"\n{'='*60}")
    print(f"Finding threshold for beta={beta}")
    print('='*60)
    find_optimal_neg_threshold(beta=beta, target_error=0.05)

When to use: Determining optimal cutoffs for approximation switches.

Method 7: Visual Debugging

Plot error as a function of input to immediately see problem regions.

import matplotlib.pyplot as plt

def visualize_precision_loss(beta=2.0, save_path=None):
    """
    Visualize where precision loss occurs for softplus in bf16.
    """
    x = np.linspace(-10, 5, 1000)
    
    # Reference (fp64)
    ref = np.log1p(np.exp(beta * x)) / beta
    
    # bf16 standard formula (simulated)
    bf16_std = []
    bf16_taylor = []
    for xi in x:
        beta_xi = beta * xi
        exp_bx = float(torch.tensor(np.exp(beta_xi), dtype=torch.bfloat16))
        one_plus = float(torch.tensor(1.0 + exp_bx, dtype=torch.bfloat16))
        log_val = float(torch.tensor(np.log(max(one_plus, 1e-30)), dtype=torch.bfloat16))
        bf16_std.append(log_val / beta)
        bf16_taylor.append(exp_bx / beta)
    
    bf16_std = np.array(bf16_std)
    bf16_taylor = np.array(bf16_taylor)
    
    # Calculate errors
    err_std = np.abs(bf16_std - ref) / (np.abs(ref) + 1e-30)
    err_taylor = np.abs(bf16_taylor - ref) / (np.abs(ref) + 1e-30)
    
    # Plot
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
    
    # Top: Function values
    ax1.plot(x, ref, 'b-', label='Reference (fp64)', linewidth=2)
    ax1.plot(x, bf16_std, 'r--', label='bf16 standard', alpha=0.7)
    ax1.plot(x, bf16_taylor, 'g:', label='bf16 Taylor', alpha=0.7)
    ax1.set_ylabel('softplus(x)')
    ax1.set_title(f'Softplus Computation (beta={beta})')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Bottom: Relative error
    ax2.semilogy(x, err_std, 'r-', label='Standard formula error', linewidth=1.5)
    ax2.semilogy(x, err_taylor, 'g-', label='Taylor approx error', linewidth=1.5)
    ax2.axvline(x=-3/beta, color='purple', linestyle='--', label=f'neg_thresh=-3 (x={-3/beta:.1f})')
    ax2.axvline(x=-5/beta, color='orange', linestyle='--', label=f'neg_thresh=-5 (x={-5/beta:.1f})')
    ax2.axhline(y=0.06, color='gray', linestyle=':', label='6% threshold')
    ax2.axhline(y=0.08, color='gray', linestyle='-.', label='8% threshold')
    ax2.set_xlabel('x')
    ax2.set_ylabel('Relative Error')
    ax2.set_ylim(1e-6, 1)
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    
    if save_path:
        plt.savefig(save_path, dpi=150, bbox_inches='tight')
        print(f"Saved plot to {save_path}")
    
    plt.show()

# Generate plots for different beta values
# visualize_precision_loss(beta=1.0, save_path='softplus_precision_beta1.png')
# visualize_precision_loss(beta=2.0, save_path='softplus_precision_beta2.png')

When to use: Quick visual identification of problem regions.

Recommended Workflow

1. Start with visual analysis - Plot errors to get a quick overview
2. Use error propagation - Understand WHERE in the computation precision is lost
3. Compare against fp64 reference - Quantify the actual error
4. Find optimal thresholds - Use automated sweep for approximation cutoffs
5. Validate with ULP analysis - Ensure errors are within acceptable ULP bounds
6. Test edge cases - Verify at boundaries and extreme values

Quick Checklist for New Ops

When implementing a new op for bfloat16:

• Test against fp32/fp64 reference across full input range
• Identify operations that can cause absorption (adding small to large)
• Identify operations that amplify errors (exp, log, division by small numbers)
• Check condition number at boundaries
• Consider Taylor series approximations for edge cases
• Determine if intermediate promotion to fp32 is needed
• Document numerical stability considerations in code comments

Example: Applying to Softplus

The softplus precision issue was detected and resolved using these methods:

1. Visual analysis showed errors spiking for negative inputs
2. Error propagation revealed 1 + exp(beta*x) rounds to 1.0 for small exp values
3. Threshold finder determined -3 is optimal cutoff for Taylor switch
4. Reference comparison validated the fix achieves < 8% relative error

See the companion Softplus kernel implementation notes for the full analysis.