第 3 章数值计算基础

在上一章中，我们已经看到：统计公式写在纸上是一回事，把它稳定地交给计算机执行是另一回事。本章进一步讨论计算统计中最基本的数值问题。它们看起来有些“底层”，但会反复出现在后面的线性代数、优化、Monte Carlo、Bootstrap 和贝叶斯计算中。

本章不试图完整讲授数值分析，而是帮助读者建立几种必要的计算直觉：

计算机中的数不是数学中的实数；
等价的数学公式可能具有完全不同的数值表现；
数据规模会改变算法选择；
写出能运行的代码不够，还要知道它是否稳定、是否高效。

3.1 计算机中的数字

数学中的实数有无限多个，但计算机只能表示有限多个数字。R 默认使用双精度浮点数（double precision floating-point numbers）来存储数值型对象。双精度浮点数大致可以理解为用有限位数记录一个数的有效数字和数量级。因此，它能表示很大或很小的数，但不能表示所有实数，也不能保证每一次四则运算都完全精确。

R 中 .Machine 保存了当前平台上浮点数的一些重要信息。

.Machine[c("double.eps", "double.xmin", "double.xmax")]
#> $double.eps
#> [1] 2.22e-16
#> 
#> $double.xmin
#> [1] 2.225e-308
#> 
#> $double.xmax
#> [1] 1.798e+308

其中，.Machine$double.eps 可以理解为机器能分辨的相对精度量级；.Machine$double.xmin 是标准双精度浮点数中最小的正规格化数；.Machine$double.xmax 是最大有限数。它们提醒我们：计算机能表示的数字是有边界、有精度限制的。

一个很经典的例子是：

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

这里并不是 R 算错了，而是 0.1、0.2 和 0.3 在二进制浮点系统中不能都被精确表示。比较浮点数时，通常不应直接使用 ==，而应使用允许小误差的比较方式，例如 all.equal()。

3.2 舍入误差与有效数字

浮点数只有有限精度。当一个很大的数和一个很小的数相加时，小数可能完全“看不见”。

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

从数学上看，第一行应该得到 1，但计算机返回的结果可能是 0。这是因为在 $10^{16}$ 这个数量级上，浮点数已经无法分辨相差 1 的两个数。

舍入误差还会在相减时被放大。当两个非常接近的数相减时，前面的有效数字会相互抵消，剩下的结果可能只包含很少的有效信息。这种现象称为灾难性抵消（catastrophic cancellation）。

例如，考虑函数

\[ f(x)=\sqrt{x+1}-\sqrt{x}. \]

当 $x$ 很大时，两个平方根非常接近，直接相减容易损失精度。利用代数恒等式，我们可以把它改写为

\[ f(x)=\frac{1}{\sqrt{x+1}+\sqrt{x}}. \]

两种公式在数学上等价，但数值表现不同。

x <- 1e12
direct <- sqrt(x + 1) - sqrt(x)
stable <- 1 / (sqrt(x + 1) + sqrt(x))
c(direct = direct, stable = stable, relative_error = abs(direct - stable) / stable)
#>         direct         stable relative_error 
#>      5.000e-07      5.000e-07      7.614e-06

这个例子很重要：统计计算中常常不是“公式对不对”的问题，而是“这个公式是否适合直接计算”的问题。

3.3 上溢、下溢与对数尺度

如果计算结果超出了浮点数可表示范围，就会出现上溢（overflow）或下溢（underflow）。例如：

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

exp(1000) 超过了双精度浮点数的表示范围，返回 Inf；exp(-1000) 太接近 0，返回 0。在统计计算中，这类问题尤其常见，因为似然函数经常包含指数函数和大量概率密度的乘积。

解决这类问题的基本策略是尽量在对数尺度上计算。以正态密度为例：

dnorm(40)
#> [1] 0
dnorm(40, log = TRUE)
#> [1] -800.9

直接密度值发生下溢，但对数密度仍然可以稳定表示。后面的极大似然估计、重要性抽样和 MCMC 中都会频繁使用对数概率。

另一个常用技巧是 log-sum-exp。假设我们要计算

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

如果 $a_i$ 都是很小的负数，先计算 exp(a_i) 再求和可能会全部下溢。

a <- c(-1000, -1001, -1002)
log(sum(exp(a)))
#> [1] -Inf

m <- max(a)
m + log(sum(exp(a - m)))
#> [1] -999.6

第二种写法利用了恒等式

\[ \log\left(\sum_i \exp(a_i)\right) = m + \log\left(\sum_i \exp(a_i-m)\right), \]

其中 $m=\max_i a_i$。这样做不会改变数学结果，却显著改善数值稳定性。

3.4 数值稳定性

如果一个算法在小的输入扰动或舍入误差下仍能给出可靠结果，我们说它具有较好的数值稳定性。稳定性不是抽象概念，它直接影响统计结果。

考虑样本方差。数学上，

\[ \operatorname{Var}(X)=E(X^2)-[E(X)]^2. \]

因此，有人可能会用 mean(x^2) - mean(x)^2 计算方差。但当数据的均值很大、方差很小时，这种写法容易发生有效数字抵消。

set.seed(1)
x <- 1e8 + rnorm(100000)

naive_var <- mean(x^2) - mean(x)^2
stable_var <- var(x)

c(naive = naive_var, stable = stable_var)
#>  naive stable 
#>  2.000  1.007

var() 的实现比上面的朴素公式更稳定。这个例子也说明，R 中很多基础函数已经封装了较好的数值算法。学习计算统计并不是要求我们所有算法都从头写，而是要知道什么时候可以信任现有函数、什么时候需要理解其背后的计算原理。

另一个常见问题是矩阵接近奇异。在线性回归中，如果自变量之间高度相关，矩阵 $X^\top X$ 会变得病态，直接求逆可能不可靠。

A <- matrix(c(1, 1, 1, 1 + 1e-12), 2, 2)
kappa(A)
#> [1] 3.999e+12
solve(A, c(1, 1))
#> [1] 1 0

kappa() 给出矩阵条件数的估计。条件数越大，问题越病态，计算结果对输入误差越敏感。第 4 章会进一步讨论矩阵分解与稳定计算。

3.5 中心化与标准化

很多统计模型在计算前都会对变量进行中心化或标准化。中心化是减去均值，标准化是在中心化后再除以标准差。它们并不会改变变量中包含的信息，却能改善数值尺度，使优化和矩阵计算更稳定。

set.seed(2)
income <- exp(rnorm(1000, mean = 10, sd = 0.8))
income_scaled <- scale(log(income))

summary(log(income))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    7.82    9.49   10.04   10.05   10.62   12.41
summary(as.vector(income_scaled))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -2.7428 -0.6831 -0.0117  0.0000  0.6986  2.9034

在经济统计中，收入、资产、销售额等变量常常具有很大的尺度差异，并且分布右偏。先取对数，再标准化，往往能让模型估计更稳定，也让优化算法更容易收敛。当然，是否进行这些变换还要结合研究问题和解释需求。

3.6 计算复杂度与存储

同一个统计问题，在小数据上可能几乎瞬间完成，在大数据上却可能变得不可行。算法所需的计算时间和存储空间会随着样本量、变量个数和模型复杂度增长。

我们通常用复杂度的大致量级来描述算法规模。例如：

遍历一个长度为 $n$ 的向量通常是 $O(n)$；
两个 $n\times n$ 矩阵相乘通常约为 $O(n^3)$；
存储一个 $n\times p$ 的数据矩阵需要 $O(np)$ 的空间。

这些符号不追求精确秒数，而是帮助我们判断当数据规模扩大 10 倍或 100 倍时，程序是否还能运行。

存储也同样重要。下面看几个矩阵对象大致需要多少内存。

object.size(matrix(0, nrow = 1000, ncol = 1000))
#> 8000216 bytes
object.size(matrix(0, nrow = 10000, ncol = 100))
#> 8000216 bytes

一个双精度数通常占 8 字节。一个 $10000 \times 100$ 的矩阵已经有一百万个数字。真实数据分析中，还会有多个中间对象、模型对象和图形对象，因此内存占用往往比原始数据大得多。

3.7 R 中的向量化

R 是面向向量的语言。很多时候，向量化写法不仅更简洁，也更快。下面比较用循环和向量化方式计算平方和。

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

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

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

c(loop = sum_loop(x), vectorized = sum_vectorized(x))
#>       loop vectorized 
#>     100818     100818
system.time(sum_loop(x))
#>    user  system elapsed 
#>   0.001   0.000   0.002
system.time(sum_vectorized(x))
#>    user  system elapsed 
#>   0.001   0.000   0.000

这并不意味着循环永远不能用。循环在模拟、迭代算法和 MCMC 中非常常见。关键在于：当一个操作本质上可以交给底层向量化函数完成时，通常应优先使用已有的向量化实现；当算法本身是递推或迭代过程时，清晰的循环反而更自然。

3.8 R 中的性能评估

在优化代码前，应该先测量。R 中最基本的计时工具是 system.time()。如果只想快速比较几种实现，它已经足够。

set.seed(4)
X <- matrix(rnorm(5000 * 20), 5000, 20)
y <- rnorm(5000)

system.time(crossprod(X, y))
#>    user  system elapsed 
#>       0       0       0
system.time(t(X) %*% y)
#>    user  system elapsed 
#>   0.001   0.000   0.000

crossprod(X, y) 和 t(X) %*% y 在数学上等价，但前者专门为交叉乘积设计，通常更直接。类似地，tcrossprod()、rowSums()、colMeans() 等函数都值得优先考虑。

不过，性能评估也要谨慎。一次运行的时间会受系统状态影响，短小代码的计时尤其不稳定。正式比较时，可以多运行几次，也可以使用专门的基准测试包。本书中只要求读者掌握一个基本习惯：不要凭感觉判断代码快慢，先测量，再优化。

3.9 本章小结

本章讨论了统计计算中的几个基础问题。浮点数有限精度意味着计算机不能完全等同于数学实数；舍入误差、上溢、下溢和灾难性抵消会影响看似简单的计算；对数尺度、稳定公式、中心化、标准化和矩阵分解可以显著提高计算可靠性；数据规模决定了算法的时间和存储成本；R 中的向量化和专用矩阵函数常常能带来更好的性能。

后续章节会不断回到这些原则。线性代数章节会讨论矩阵分解和病态问题；优化章节会讨论迭代算法和收敛；Monte Carlo 与 MCMC 章节会讨论模拟误差；统计模型实现章节会把这些工具放回实际建模流程中。

3.10 练习

运行 0.1 + 0.2 == 0.3 和 all.equal(0.1 + 0.2, 0.3)，解释二者结果为什么不同。
对不同的 $x$ 值比较 sqrt(x + 1) - sqrt(x) 和 1 / (sqrt(x + 1) + sqrt(x))，观察何时差异变得明显。
令 a <- c(-100, -500, -1000)，分别用直接方法和 log-sum-exp 方法计算 log(sum(exp(a)))。
构造两个高度相关的变量，计算其设计矩阵的 kappa()，并尝试拟合一个线性回归模型。
用 system.time() 比较 sum(x)、手写循环求和、mean(x) * length(x) 三种写法在大向量上的运行时间。

第 3 章 数值计算基础