Bayesian Statistics and Computing

---

# Bayesian Statistics and Computing

## Lecture 1: Numerical Basics

#### *Yanfei Kang | BSC | Beihang University*

---
# Objectives

By the end of this lecture, you should be able to:

1. Explain double precision, machine precision, overflow, and underflow.
2. Recognize rounding error and catastrophic cancellation.
3. Distinguish conditioning from numerical stability.
4. Explain the numerical and statistical effects of centering and standardization.
5. Use big `$O$` notation for time and storage complexity.
6. Compare equivalent R implementations with `microbenchmark`.

---
# Machine parameters in R

```r
.Machine[c("double.eps", "double.xmin", "double.xmax")]
#> $double.eps
#> [1] 2.220446e-16
#> 
#> $double.xmin
#> [1] 2.225074e-308
#> 
#> $double.xmax
#> [1] 1.797693e+308
```

- `double.eps` is the difference between 1 and the next larger floating point number.
- The unit roundoff is approximately `double.eps / 2` under rounding to nearest.
- `double.xmin` is the smallest positive normal double. Subnormal numbers can be smaller but have reduced relative precision.
- `double.xmax` is the largest finite double.

---
# Decimal values in binary arithmetic

```r
0.1 + 0.2
#> [1] 0.3
0.1 + 0.2 == 0.3
#> [1] FALSE
isTRUE(all.equal(0.1 + 0.2, 0.3))
#> [1] TRUE
```

`0.1`, `0.2`, and `0.3` all require rounding in binary floating point arithmetic.

`all.equal()` returns `TRUE` for approximate equality and otherwise describes the difference. Use `isTRUE(all.equal(...))` inside a logical condition.

The tolerance must reflect the problem scale. Exact comparison with `==` can still be appropriate for integer counts, category codes, and exact flags.

---
# Absolute and relative error

If the exact target is `$z$` and the computed value is `$\hat z$`,

`$$e_{\mathrm{abs}}=|\hat z-z|,\qquad e_{\mathrm{rel}}=\frac{|\hat z-z|}{|z|}.$$`

Relative error is defined only when `$z\neq 0$`.

- Absolute error gives the error on the original scale.
- Relative error relates the error to the size of the target.
- When the exact value is unknown, use a high precision reference, an analytic result, a residual, or an error bound.

---
# Adding a small number to a large number

```r
x <- 1e16
(x + 1) - x
#> [1] 0
(x + 10000) - x
#> [1] 10000
```

In real arithmetic, the first expression equals 1.

Near `$10^{16}$`, adjacent double precision numbers are more than 1 apart. Therefore, `x + 1` and `x` can have the same floating point representation.

---
# Catastrophic cancellation

Subtracting two close values removes their shared leading digits. Small rounding errors in the inputs can then become large relative to the result.

For example,

`$$f(x)=\sqrt{x+1}-\sqrt{x}.$$`

For large `$x$`, the two square roots are close. Multiplying by the conjugate gives

`$$f(x)=\frac{1}{\sqrt{x+1}+\sqrt{x}}.$$`

The two formulas are mathematically equivalent but can behave differently in finite precision.

---
# Cancellation in R

```r
x <- 1e12
direct <- sqrt(x + 1) - sqrt(x)
stable <- 1 / (sqrt(x + 1) + sqrt(x))

c(
  direct = direct,
  stable = stable,
  relative_difference = abs(direct - stable) / abs(stable)
)
#>              direct              stable relative_difference 
#>        5.000038e-07        5.000000e-07        7.614494e-06
```

The second expression avoids direct subtraction of two close numbers.

This is a relative difference between two implementations, not a relative error against an unknown exact value.

---
# Overflow and underflow

```r
exp(1000)
#> [1] Inf
exp(-1000)
#> [1] 0
```

- Overflow occurs when the magnitude exceeds the largest finite floating point number.
- Underflow occurs when a nonzero result is too small to represent at the required precision.

`exp(1000)` returns `Inf`, while `exp(-1000)` underflows to zero.

---
# Likelihoods on the log scale

For independent observations,

`$$L(\theta)=\prod_{i=1}^n f(y_i\mid\theta).$$`

Products of many densities or probabilities can underflow. The log-likelihood is

`$$\ell(\theta)=\log L(\theta)=\sum_{i=1}^n\log f(y_i\mid\theta).$$`

It replaces a product by a sum and preserves the same maximizer.

Bayesian computations also combine priors, likelihoods, posterior kernels, MCMC acceptance probabilities, and importance weights on the log scale.

---
# A normal density in the tail

```r
dnorm(40)
#> [1] 0
dnorm(40, log = TRUE)
#> [1] -800.9189
```

The density underflows to zero, while the log-density remains representable.

The log scale prevents premature overflow or underflow. It does not solve problems caused by a poor model, poor data, or ill-conditioning.

---
# The log-sum-exp identity

Consider

`$$\log\left\{\sum_{i=1}^n\exp(a_i)\right\}.$$`

Let `$m=\max_{1\leq i\leq n}a_i$`. Then

`$$\log\left\{\sum_{i=1}^n\exp(a_i)\right\}=m+\log\left\{\sum_{i=1}^n\exp(a_i-m)\right\}.$$`

Every `$a_i-m$` is nonpositive, and at least one shifted exponential equals 1.

---
# Log-sum-exp in R

```r
log_sum_exp <- function(a) {
  m <- max(a)
  if (is.infinite(m)) {
    return(m)
  }
  m + log(sum(exp(a - m)))
}

a <- c(-1000, -1001, -1002)
c(
  direct = log(sum(exp(a))),
  stable = log_sum_exp(a)
)
#>  direct  stable 
#>    -Inf -999.59
```

The function assumes a nonempty input without missing values. The special case handles a maximum of `Inf` or an input containing only `-Inf`.

---
# Conditioning and numerical stability

| Concept | Meaning | Property of |
|:---|:---|:---|
| **Conditioning** | Sensitivity of the exact answer to small input perturbations | The mathematical problem |
| **Numerical stability** | Propagation or amplification of rounding error | The computational method |

An ill-conditioned problem can lose accuracy even with a stable algorithm. An unstable algorithm can lose accuracy on a well-conditioned problem.

For a linear system, a large `kappa(A)` indicates that small input changes may produce much larger changes in the solution.

---
# A nearly singular system

```r
A <- matrix(
  c(1, 1, 1, 1 + 1e-10),
  nrow = 2,
  byrow = TRUE
)
b1 <- c(2, 2 + 1e-10)
b2 <- b1 + c(0, 1e-12)

solution1 <- solve(A, b1)
solution2 <- solve(A, b2)
relative_change <- function(new, old) {
  sqrt(sum((new - old)^2)) / sqrt(sum(old^2))
}
```

Only the second component of `$\mathbf b$` changes, by `$10^{-12}$`.

---
# Sensitivity of the solution

```r
c(
  condition_number = kappa(A, exact = TRUE),
  relative_input_change = relative_change(b2, b1),
  relative_solution_change = relative_change(solution2, solution1)
)
#>         condition_number    relative_input_change relative_solution_change 
#>               4.0000e+10               3.5358e-13               1.0001e-02
cbind(solution1, solution2)
#>      solution1 solution2
#> [1,]         1      0.99
#> [2,]         1      1.01
```

The relative input change is about `$3.5\times10^{-13}$`, while the relative solution change is about one percent.

The main difficulty comes from the ill-conditioned linear system, not from whether `solve()` returns a result.

---
# Variance as a stability example

For observations `$x_1,\ldots,x_n$`, the second central moment with denominator `$n$` is

`$$v_n=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2=\frac{1}{n}\sum_{i=1}^n x_i^2-\bar{x}^2.$$`

The raw-moment formula is algebraically correct. When the mean is large and the spread is small, it subtracts two large, nearly equal numbers.

The formulas must use the same denominator before their numerical results are compared.

---
# Variance calculations in R

```r
set.seed(1)
x <- 1e8 + rnorm(100000)
n <- length(x)

variance_naive <- mean(x^2) - mean(x)^2
variance_centered <- mean((x - mean(x))^2)
variance_from_var <- var(x) * (n - 1) / n

c(
  naive_moment_formula = variance_naive,
  centered_formula = variance_centered,
  adjusted_var = variance_from_var
)
#> naive_moment_formula     centered_formula         adjusted_var 
#>                2.000                1.007                1.007
```

`var(x)` uses denominator `$n-1$`, so the code adjusts it to denominator `$n$` before comparing the results.

---
# Centering and standardization

Centering subtracts a location such as the sample mean. Standardization also divides by a scale such as the sample standard deviation.

For an unpenalized linear model with an intercept, nondegenerate centering and scaling preserve the column space when derived terms are handled consistently. Coefficients can be converted back to their original units, and fitted values can remain unchanged.

The scale and condition number faced by the numerical method can change substantially.

---
# Scaling a design matrix

```r
set.seed(2)
x_small <- rnorm(1000)
x_large <- 1e6 * rnorm(1000)

X_raw <- cbind(intercept = 1, x_small, x_large)
X_scaled <- cbind(
  intercept = 1,
  scale(cbind(x_small, x_large))
)

c(
  raw_design = kappa(X_raw),
  scaled_design = kappa(X_scaled)
)
#>    raw_design scaled_design 
#>    9.8420e+05    1.0448e+00
```

---
# Computational complexity

- Scanning a vector of length `$n$` usually costs `$O(n)$` operations.
- Multiplying two dense `$n\times n$` matrices by a classical algorithm costs `$O(n^3)$` operations.
- Storing a dense `$n\times p$` double matrix requires `$O(np)$` space.
- A dense QR factorization of an `$n\times p$` matrix with `$n\geq p$` usually costs `$O(np^2)$` operations.

Big `$O$` notation retains the dominant growth rate as the problem size increases.

It does not give an exact runtime or capture constants, caches, parallel execution, and numerical libraries.

---
# Object size and peak memory

```r
format(object.size(matrix(0, nrow = 1000, ncol = 1000)), units = "auto")
#> [1] "7.6 Mb"
format(object.size(matrix(0, nrow = 10000, ncol = 100)), units = "auto")
#> [1] "7.6 Mb"
```

Both matrices contain one million doubles. Their raw numerical values require about 8 MB.

A real analysis may also retain data copies, matrix factorizations, simulation draws, model objects, and plots. Peak memory can therefore be much larger than the input file.

---
# Three sum-of-squares methods

```r
library(microbenchmark)

set.seed(3)
x_bench <- rnorm(100000)

sum_squares_loop <- function(x) {
  out <- 0
  for (i in seq_along(x)) {
    out <- out + x[i]^2
  }
  out
}

sum_squares_vectorized <- function(x) {
  sum(x^2)
}

sum_squares_crossprod <- function(x) {
  drop(crossprod(x))
}
```

---
# Equality check

```r
stopifnot(
  isTRUE(all.equal(
    sum_squares_loop(x_bench),
    sum_squares_vectorized(x_bench)
  )),
  isTRUE(all.equal(
    sum_squares_vectorized(x_bench),
    sum_squares_crossprod(x_bench)
  ))
)
```

The benchmark should start only after the implementations have been shown to compute the same quantity.

---
# Sum-of-squares benchmark

```r
sumsq_benchmark <- suppressWarnings(
  microbenchmark(
    `Loop` = sum_squares_loop(x_bench),
    `Vectorized expression` = sum_squares_vectorized(x_bench),
    `Specialized function: crossprod` = sum_squares_crossprod(x_bench),
    times = 100L,
    unit = "ms",
    control = list(order = "random", warmup = 10L)
  )
)
```

---
# Benchmark results

```r
sumsq_summary <- summary(sumsq_benchmark, unit = "ms")
sumsq_timing <- data.frame(
  method = as.character(sumsq_summary$expr),
  median_ms = sumsq_summary$median,
  relative_time = sumsq_summary$median / min(sumsq_summary$median)
)
knitr::kable(
  sumsq_timing,
  row.names = FALSE,
  digits = 3,
  col.names = c("Method", "Median time (ms)", "Relative time")
)
```

|Method                          | Median time (ms)| Relative time|
|:-------------------------------|----------------:|-------------:|
|Loop                            |            4.357|        13.342|
|Vectorized expression           |            0.530|         1.621|
|Specialized function: crossprod |            0.327|         1.000|

---
# Interpreting vectorization

- `sum(x^2)` sends the loop to compiled code but creates the temporary vector `x^2`.
- `crossprod(x)` uses the inner-product structure directly.
- Actual timings depend on the R version, hardware, and numerical libraries.

Loops remain necessary for recursive algorithms, iterative optimization, and MCMC. Identify the actual bottleneck before choosing vectorization, a specialized function, or a different algorithm.

---
# Summary

- Double precision represents only a finite subset of the real numbers.
- Rounding, overflow, underflow, and cancellation can make a correct formula unreliable.
- Equivalent transformations, log-scale calculations, and specialized functions can improve evaluation.
- Conditioning describes input sensitivity. Stability describes error propagation by an algorithm.
- Centering and standardization can improve numerical scale but may change penalties, priors, and interpretation.
- Complexity describes growth with problem size. Benchmarking measures a specific implementation and environment.

---
# Lab

Choose one of the following textbook exercises:

1. Compare the direct and conjugate forms of `sqrt(x + 1) - sqrt(x)` for `$x=10^k$`, `$k=0,\ldots,20$`.
2. Improve `log_sum_exp()` so that its behavior is clear for empty inputs, `NA`, `Inf`, and inputs containing only `-Inf`.
3. Use `microbenchmark` to compare the loop, vectorized expression, and `crossprod()` for several vector lengths. Verify equality before timing.

Report the code, numerical output, and a short interpretation.

---
# References

Chapter 3 of our textbook.