numerical-methods

Numerical Methods

You are an expert in numerical stability and computational aspects of statistical methods.

Floating-Point Fundamentals

IEEE 754 Double Precision

• Precision: ~15-17 significant decimal digits
• Range: ~10⁻³⁰⁸ to 10³⁰⁸
• Machine epsilon: ε ≈ 2.2 × 10⁻¹⁶
• Special values: Inf, -Inf, NaN

Key Constants in R

.Machine$double.eps      # ~2.22e-16 (machine epsilon)
.Machine$double.xmax     # ~1.80e+308 (max finite)
.Machine$double.xmin     # ~2.23e-308 (min positive normalized)
.Machine$double.neg.eps  # ~1.11e-16 (negative epsilon)

Common Numerical Issues

1. Catastrophic Cancellation

When subtracting nearly equal numbers:

# BAD: loses precision
x <- 1e10 + 1
y <- 1e10
result <- x - y  # Should be 1, may have errors

# BETTER: reformulate to avoid subtraction
# Example: Computing variance
var_bad <- mean(x^2) - mean(x)^2   # Can be negative!
var_good <- sum((x - mean(x))^2) / (n-1)  # Always non-negative

2. Overflow/Underflow

# BAD: overflow
prod(1:200)  # Inf

# GOOD: work on log scale
sum(log(1:200))  # Then exp() if needed

# BAD: underflow in probabilities
prod(dnorm(x))  # 0 for large x

# GOOD: sum log probabilities
sum(dnorm(x, log = TRUE))

3. Log-Sum-Exp Trick

Essential for working with log probabilities:

log_sum_exp <- function(log_x) {
  max_log <- max(log_x)
  if (is.infinite(max_log)) return(max_log)
  max_log + log(sum(exp(log_x - max_log)))
}

# Example: log(exp(-1000) + exp(-1001))
log_sum_exp(c(-1000, -1001))  # Correct: ~-999.69
log(exp(-1000) + exp(-1001))   # Wrong: -Inf

4. Softmax Stability

# BAD
softmax_bad <- function(x) exp(x) / sum(exp(x))

# GOOD
softmax <- function(x) {
  x_max <- max(x)
  exp_x <- exp(x - x_max)
  exp_x / sum(exp_x)
}

Matrix Computations

Conditioning

The condition number κ(A) measures sensitivity to perturbation:

• κ(A) = ‖A‖ · ‖A⁻¹‖
• Rule: Expect to lose log₁₀(κ) digits of accuracy
• κ > 10¹⁵ means matrix is numerically singular