Numerical Methods

---


name: numerical-methods
description: Numerical algorithms and computational techniques for statistics


---

# 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
```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:
```r
# 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
```r
# 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:
```r
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
```r
# 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

```r
# Check condition number
kappa(X, exact = TRUE)

# For regression: check X'X conditioning
kappa(crossprod(X))
```

### Solving Linear Systems

**Prefer**: Decomposition methods over explicit inversion

```r
# BAD: explicit inverse
beta <- solve(t(X) %*% X) %*% t(X) %*% y

# GOOD: QR decomposition
beta <- qr.coef(qr(X), y)

# BETTER for positive definite: Cholesky
R <- chol(crossprod(X))
beta <- backsolve(R, forwardsolve(t(R), crossprod(X, y)))

# For ill-conditioned: SVD/pseudoinverse
beta <- MASS::ginv(X) %*% y
```

### Symmetric Positive Definite Matrices
Always use specialized methods:
```r
# Cholesky for SPD
L <- chol(Sigma)

# Eigendecomposition
eig <- eigen(Sigma, symmetric = TRUE)

# Check positive definiteness
all(eigen(Sigma, symmetric = TRUE, only.values = TRUE)$values > 0)
```

## Optimization Stability

### Gradient Computation
```r
# Numerical gradient (for verification)
numerical_grad <- function(f, x, h = sqrt(.Machine$double.eps)) {
  sapply(seq_along(x), function(i) {
    x_plus <- x_minus <- x
    x_plus[i] <- x[i] + h
    x_minus[i] <- x[i] - h
    (f(x_plus) - f(x_minus)) / (2 * h)
  })
}

# Central difference is O(h²) accurate
# Forward difference is O(h) accurate
```

### Hessian Stability
```r
# Check Hessian is positive definite at optimum
check_hessian <- function(H, tol = 1e-8) {
  eigs <- eigen(H, symmetric = TRUE, only.values = TRUE)$values
  min_eig <- min(eigs)

  list(
    positive_definite = min_eig > tol,
    min_eigenvalue = min_eig,
    condition_number = max(eigs) / min_eig
  )
}
```

### Line Search
For gradient descent stability:
```r
backtracking_line_search <- function(f, x, d, grad, alpha = 1, rho = 0.5, c = 1e-4) {
  # Armijo condition
  while (f(x + alpha * d) > f(x) + c * alpha * sum(grad * d)) {
    alpha <- rho * alpha
    if (alpha < 1e-10) break
  }
  alpha
}
```

## Integration and Quadrature

### Numerical Integration Guidelines
```r
# Adaptive quadrature (default choice)
integrate(f, lower, upper)

# For infinite limits
integrate(f, -Inf, Inf)

# For highly oscillatory or peaked functions
# Increase subdivisions
integrate(f, lower, upper, subdivisions = 1000)

# For known singularities, split the domain
```

### Monte Carlo Integration
```r
mc_integrate <- function(f, n, lower, upper) {
  x <- runif(n, lower, upper)
  fx <- sapply(x, f)

  estimate <- (upper - lower) * mean(fx)
  se <- (upper - lower) * sd(fx) / sqrt(n)

  list(value = estimate, se = se)
}
```

## Root Finding

### Newton-Raphson Stability
```r
newton_raphson <- function(f, df, x0, tol = 1e-8, max_iter = 100) {
  x <- x0
  for (i in 1:max_iter) {
    fx <- f(x)
    dfx <- df(x)

    # Check for near-zero derivative
    if (abs(dfx) < .Machine$double.eps * 100) {
      warning("Near-zero derivative")
      break
    }

    x_new <- x - fx / dfx

    if (abs(x_new - x) < tol) break
    x <- x_new
  }
  x
}
```

### Brent's Method
For robust root finding without derivatives:
```r
uniroot(f, interval = c(lower, upper), tol = .Machine$double.eps^0.5)
```

## Statistical Computing Patterns

### Safe Likelihood Computation
```r
# Always work with log-likelihood
log_lik <- function(theta, data) {
  # Compute log-likelihood, not likelihood
  sum(dnorm(data, mean = theta[1], sd = theta[2], log = TRUE))
}
```

### Robust Standard Errors
```r
# Sandwich estimator with numerical stability
sandwich_se <- function(score, hessian) {
  # Check Hessian conditioning
  H_inv <- tryCatch(
    solve(hessian),
    error = function(e) MASS::ginv(hessian)
  )

  meat <- crossprod(score)
  V <- H_inv %*% meat %*% H_inv

  sqrt(diag(V))
}
```

### Bootstrap with Error Handling
```r
safe_bootstrap <- function(data, statistic, R = 1000) {
  results <- numeric(R)
  failures <- 0

  for (i in 1:R) {
    boot_data <- data[sample(nrow(data), replace = TRUE), ]
    result <- tryCatch(
      statistic(boot_data),
      error = function(e) NA
    )
    results[i] <- result
    if (is.na(result)) failures <- failures + 1
  }

  if (failures > 0.1 * R) {
    warning(sprintf("%.1f%% bootstrap failures", 100 * failures / R))
  }

  list(
    estimate = mean(results, na.rm = TRUE),
    se = sd(results, na.rm = TRUE),
    failures = failures
  )
}
```

## Debugging Numerical Issues

### Diagnostic Checklist
1. **Check for NaN/Inf**: `any(is.nan(x))`, `any(is.infinite(x))`
2. **Check conditioning**: `kappa(matrix)`
3. **Check eigenvalues**: For PD matrices
4. **Check gradients**: Numerically vs analytically
5. **Check scale**: Variables on similar scales?

### Debugging Functions
```r
# Trace NaN/Inf sources
debug_numeric <- function(x, name = "x") {
  cat(sprintf("%s: range [%.3g, %.3g], ", name, min(x), max(x)))
  cat(sprintf("NaN: %d, Inf: %d, -Inf: %d\n",
              sum(is.nan(x)), sum(x == Inf), sum(x == -Inf)))
}

# Check relative error
rel_error <- function(computed, true) {
  abs(computed - true) / max(abs(true), 1)
}
```

## Best Practices Summary

1. **Always work on log scale** for products of probabilities
2. **Use QR or Cholesky** instead of matrix inversion
3. **Check conditioning** before solving linear systems
4. **Center and scale** predictors in regression
5. **Handle edge cases** (empty data, singular matrices)
6. **Use existing implementations** (LAPACK, BLAS) when possible
7. **Test with extreme values** (very small, very large, near-zero)
8. **Compare analytical and numerical gradients**
9. **Monitor convergence** in iterative algorithms
10. **Document numerical assumptions** and limitations

## Key References

- Higham
- Golub & Van Loan