第 3 章 数值计算基础

第 1 章已经说明，统计公式规定了“要计算什么”，但数学上正确的表达式并不自动给出可靠的计算方法；第 2 章则介绍了用 R 组织和复现统计计算的基本方式。本章进一步讨论计算机数制、数值误差、问题的条件性、算法的稳定性以及计算成本。这些概念构成后续数值线性代数、优化、Monte Carlo、Bootstrap 和贝叶斯计算的共同基础。

本章不以完整讲授数值分析为目标，而是建立分析统计计算方法所需的基本语言和判断标准。读者应逐步形成一个习惯：在确认公式和代码正确之后，还要检查有限精度、输入敏感性、算法实现以及时间和存储成本。

3.1 本章学习目标

完成本章后，读者应当能够：

1. 说明双精度浮点数与数学实数的差异，并正确解释机器精度、上溢和下溢；
2. 识别舍入误差与灾难性抵消，并利用等价变换或对数尺度改进求值方法；
3. 区分问题的条件性与算法的数值稳定性，解释二者对计算误差的不同影响；
4. 说明中心化和标准化如何改善变量尺度，同时识别这些变换对模型解释、惩罚项和先验分布的影响；
5. 使用大 \(O\) 记号粗略分析算法的时间复杂度和空间复杂度；
6. 使用 microbenchmark 对若干等价实现进行可重复的计时比较，并结合准确性和存储需求解释结果。

3.2 计算机中的数字

数学中的实数有无限多个，而计算机只能用有限位数表示有限多个数值。在当前主流平台上，R 中不带 L 后缀的普通数值常量以及大多数连续型统计计算通常采用双精度浮点数；整数型和复数型对象另有相应的存储类型。双精度浮点数的表示方式通常对应 IEEE 754 的 binary64 格式。一个非零正规格化浮点数由符号、有效数字和指数共同表示。因此，双精度浮点数可以覆盖很宽的数量级，却不能精确表示所有实数，也不能保证每次四则运算都等同于实数运算 (Goldberg 1991; Higham 2002)。

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 是 1 与下一个更大的可表示浮点数之差；在舍入到最近数的规则下，单位舍入误差（unit roundoff）约为它的一半。.Machine$double.xmin 是最小的正正规格化双精度数；次正规数（subnormal numbers）还可以表示更接近 0 的数，但其相对精度会降低。.Machine$double.xmax 是最大的有限双精度数。

设 \(a\) 和 \(b\) 已经是可表示的浮点数，并且精确运算结果处于正规数范围内。一次正确舍入的基本运算常用下式近似描述：

\[ \operatorname{fl}(a\circ b) = (a\circ b)(1+\delta), \qquad |\delta|\leq u, \]

其中 \(\circ\) 表示一种基本算术运算，\(\operatorname{fl}(\cdot)\) 表示浮点计算结果，\(u\) 是单位舍入误差。这个模型说明，单次运算的误差通常很小，但误差可能在长序列运算中累积，也可能被问题本身或不稳定算法显著放大。

十进制小数在二进制浮点系统中未必能够精确表示。例如：

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

这里并不是 R 进行了错误的实数运算，而是 0.1、0.2 和 0.3 的二进制浮点表示均含有舍入。当两个数在理论上应当相等、但由浮点运算得到时，通常应比较二者之差是否小于与问题尺度相适应的容许误差。all.equal() 在近似相等时返回 TRUE，否则返回描述差异的字符串；因此，在程序的逻辑判断中通常写成 isTRUE(all.equal(...))。容许误差的选择仍应结合具体问题，而不能机械地套用同一个阈值。对于整数计数、类别编码或专门设置的精确标记值，直接使用 == 仍可能是适当的。

3.3 舍入误差、相对误差与灾难性抵消

若目标量的精确值为 \(z\)，计算得到的近似值为 \(\hat z\)，则绝对误差和相对误差分别为

\[ e_{\mathrm{abs}}=|\hat z-z|, \qquad e_{\mathrm{rel}}=\frac{|\hat z-z|}{|z|}, \]

其中相对误差只在 \(z\neq 0\) 时定义。绝对误差反映误差的原始量级，相对误差则把误差与目标量的尺度联系起来。实际统计问题中的精确值往往未知，此时只能借助高精度参考值、解析结果、残差或误差上界评估准确性。因此，下面代码报告的是两种实现之间的相对差异，而不是相对于未知真值的严格相对误差。

浮点数的有效位数有限。当两个数的数量级相差很大时，把较小的数加到较大的数上，结果可能不再发生变化。

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

在实数算术中，第一个表达式等于 1；在双精度浮点算术中，\(10^{16}\) 附近相邻可表示数之间的距离已经大于 1，因此 x + 1 可能与 x 具有相同的浮点表示。

舍入误差在两个非常接近的数相减时尤其值得注意。相减会消去二者共同的前导有效数字，使输入中原本很小的舍入误差在结果中占据较大比例。这一现象称为灾难性抵消（catastrophic cancellation）(Wilkinson 1963; Higham 2002)。

例如，考虑

\[ 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_difference = abs(direct - stable) / abs(stable)
)
#>              direct              stable relative_difference 
#>           5.000e-07           5.000e-07           7.614e-06

第二种写法避免了两个相近数的直接相减。这个例子揭示了数值计算的一个基本原则：数学等价并不意味着数值等价；表达式的书写形式和实际求值方法需要分别考察。

3.4 上溢、下溢与对数尺度

当计算结果的绝对值超过最大有限浮点数时，会发生上溢（overflow）；当非零结果的绝对值过小而不能以所需精度表示时，会发生下溢（underflow）。

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

exp(1000) 超出双精度数的表示范围并返回 Inf，而 exp(-1000) 下溢为 0。在统计计算中，这类问题尤其常见。若观测相互独立，似然通常写为

\[ L(\theta)=\prod_{i=1}^n f(y_i\mid\theta), \]

大量小于 1 的密度值或概率相乘很容易下溢。相应的对数似然为

\[ \ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log f(y_i\mid\theta), \]

它把乘积转化为求和，并保留同一个极大值点。在贝叶斯计算中，先验密度、似然与后验核也通常在对数尺度上组合；MCMC 的接受概率和重要性权重同样需要避免直接计算极小概率的乘积 (Gentle 2009; Monahan 2011)。

第 1 章曾用正态密度初步展示下溢现象。这里从对数尺度计算的角度再次考察该例：

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

直接密度值已经下溢为 0，但对数密度仍能稳定表示。需要注意，对数尺度避免了过早的上溢或下溢，却不能消除模型不恰当、数据质量不佳或问题病态所造成的困难。

另一个常用技巧是 log-sum-exp。考虑

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

若某些 \(a_i\) 很大，exp(a_i) 可能上溢；若所有 \(a_i\) 都是绝对值很大的负数，则各项又可能下溢。令\(m=\max_{1\leq i\leq n}a_i\)，可利用恒等式

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

由于所有 \(a_i-m\leq 0\)，指数项不会上溢，并且至少有一项等于 1。

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.6

函数中的特殊判断处理了最大元素为 Inf 或所有元素均为 -Inf 的边界情形；这里暂定输入非空且不含缺失值。后续的混合模型、重要性抽样、序贯 Monte Carlo 和贝叶斯模型比较都会反复使用这一变换。

3.5 条件性与数值稳定性

数值误差的来源至少包含两个层面。条件性（conditioning）描述问题本身对输入扰动的敏感程度；数值稳定性（numerical stability）描述算法如何传播或放大舍入误差。条件性是数学问题的属性，稳定性是计算方法的属性，二者不能混为一谈 (Higham 2002)。

3.5.1 问题的条件性

设问题把输入 \(z\) 映射为输出 \(g(z)\)。如果很小的输入扰动可能造成很大的输出变化，则该问题是病态的（ill-conditioned）；反之则是良态的（well-conditioned）。病态问题即使使用稳定算法，也可能无法从有限精度或含有测量误差的输入中得到高精度答案。

在线性方程组 \(\mathbf A\mathbf x=\mathbf b\) 中，给定一种相容的矩阵范数，非奇异矩阵\(\mathbf A\) 的条件数可定义为

\[ \kappa(\mathbf A)=\lVert\mathbf A\rVert\,\lVert\mathbf A^{-1}\rVert. \]

当只扰动非零右端向量，并且原解非零时，在相容范数下有

\[ \frac{\lVert\Delta\mathbf x\rVert}{\lVert\mathbf x\rVert} \leq \kappa(\mathbf A) \frac{\lVert\Delta\mathbf b\rVert}{\lVert\mathbf b\rVert}. \]

因此，较大的条件数意味着输入中的相对误差可能被显著放大。下面构造一个接近奇异的矩阵，并只对\(\mathbf b\) 的第二个分量施加 \(10^{-12}\) 的扰动。代码中的 kappa(A, exact = TRUE) 计算基于 2 范数的条件数；若使用默认设置，kappa() 返回的是条件数估计。

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))
}

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.000e+10                3.536e-13                1.000e-02
cbind(solution1, solution2)
#>      solution1 solution2
#> [1,]         1      0.99
#> [2,]         1      1.01

右端向量的相对变化约为 \(3.5\times 10^{-13}\)，但解的相对变化可达到百分之一的量级。这里的主要困难来自线性方程组本身的病态性，而不是 solve() 是否成功返回结果。第 4 章将进一步讨论条件数、数值秩和矩阵分解 (Golub and Van Loan 2013)。

3.5.2 算法的数值稳定性

稳定算法应当控制运算过程中舍入误差的传播。一种常见判断方式是后向误差分析：如果计算结果可以解释为某个轻微扰动后的输入所对应的精确解，则称算法具有后向稳定性。后向稳定并不保证前向误差一定很小；当问题病态时，轻微的输入扰动本身就可能对应很大的输出变化。

样本方差说明了算法形式的重要性。对于观测值 \(x_1,\ldots,x_n\)，分母为 \(n\) 的二阶中心矩为

\[ 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. \]

右端的“原始二阶矩减均值平方”公式在代数上正确，但当均值很大而离散程度很小时，需要把两个很大的近似数相减，容易发生灾难性抵消。为了进行口径一致的比较，下面把 R 的样本方差 var(x) 从分母 \(n-1\) 调整为分母 \(n\)。

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

中心化公式避免了两个约为 \(10^{16}\) 的量直接相减；var() 返回分母为 \(n-1\) 的样本方差，经口径调整后应与中心化结果接近。这个例子还说明，比较两种计算方法之前，必须先确认它们计算的是同一个统计量。

3.6 中心化、标准化与尺度问题

中心化是从变量中减去某个位置量，标准化通常是在中心化后再除以一个尺度量。对于含截距、未加惩罚的线性模型，如果对数值自变量实施非退化的中心化和尺度变换，并一致处理交互项等派生项，那么设计矩阵的列空间保持不变，系数可以换算回原单位，拟合值也可以保持不变；但计算所面对的尺度和条件数可能显著改善。

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.842e+05     1.045e+00

这个例子中的两个变量具有相近的随机结构，但数值尺度相差约 \(10^6\) 倍。标准化后的设计矩阵具有更均衡的列尺度，因此条件数明显减小。在梯度法等迭代算法中，尺度差异还会造成目标函数在不同方向上的曲率差异，使统一步长难以选择。

数据变换仍然首先是建模决定。例如，对收入、资产或销售额取对数会改变模型的函数形式和参数解释，不能仅被视为数值技巧。对于岭回归、Lasso 等带惩罚方法，变量尺度会改变相同惩罚参数对各系数的实际约束；在贝叶斯模型中，尺度也会影响先验分布的含义和后验几何形状。因此，中心化、标准化和其他变换的规则应当在分析中明确记录，并在预测新数据时保持一致。

3.7 计算复杂度与存储

同一个统计问题在小规模数据上可能瞬间完成，在大规模数据上却可能不可行。算法所需的计算时间和存储空间会随观测数、变量数、参数维数以及模型结构增长。大 \(O\) 记号保留问题规模增长时的主导项，用于描述这种增长速度。例如：

• 遍历一个长度为 \(n\) 的向量通常需要 \(O(n)\) 次运算；
• 按经典算法相乘两个稠密的 \(n\times n\) 矩阵通常需要 \(O(n^3)\) 次运算；
• 存储一个稠密的 \(n\times p\) 双精度数据矩阵需要 \(O(np)\) 的空间；
• 当 \(n\geq p\) 时，对一个稠密的 \(n\times p\) 设计矩阵进行 QR 分解通常需要 \(O(np^2)\) 次运算。

大 \(O\) 记号不提供精确运行时间，也会忽略常数、缓存、并行化和底层线性代数库等因素。它的主要用途是判断规模扩大后某一方法是否仍然可行，而不是代替实际计时。

下面比较两个元素总数相同、形状不同的矩阵所占用的对象空间。

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"

双精度数通常占 8 字节，因此一百万个数的原始数值存储约为 8 MB。真实分析还可能同时保留数据副本、矩阵分解、模拟样本、模型对象和图形对象，峰值内存往往远高于原始数据文件的大小。判断算法能否运行时，应当估计中间对象，而不能只查看输入数据。

3.8 向量化与专用函数

R 是面向向量的语言。向量化表达式通常比逐元素解释执行的 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))
}

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)
  ))
)

sumsq_benchmark <- suppressWarnings(
  microbenchmark(
    `手写循环` = sum_squares_loop(x_bench),
    `向量化表达式` = sum_squares_vectorized(x_bench),
    `专用函数 crossprod` = sum_squares_crossprod(x_bench),
    times = 100L,
    unit = "ms",
    control = list(order = "random", warmup = 10L)
  )
)

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("方法", "耗时中位数（毫秒）", "相对耗时")
)

方法	耗时中位数（毫秒）	相对耗时
手写循环	5.995	12.971
向量化表达式	0.742	1.605
专用函数 crossprod	0.462	1.000

sum(x^2) 把循环交给底层实现，但会生成临时向量 x^2；crossprod(x) 直接利用内积结构计算平方和。具体速度取决于 R 版本、硬件和底层线性代数库，因此表中的数值只描述当前计算环境。这个例子也不意味着循环应当被一概避免：递推算法、迭代优化和 MCMC 天然包含循环。优化代码时，应优先识别真正的性能瓶颈，再判断是否存在等价的向量化表达式、专用函数或更合适的算法。

3.9 R 中的性能评估

性能判断应当建立在测量之上。system.time() 适合测量一次较长任务的用户时间、系统时间和总耗时；对于运行时间很短的若干表达式，microbenchmark() 通过重复执行提供更稳定的比较。下面比较两种计算\(\mathbf X^\top\mathbf y\) 的方法。

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

crossprod_result <- drop(crossprod(X, y))
explicit_result <- drop(t(X) %*% y)
stopifnot(isTRUE(all.equal(crossprod_result, explicit_result)))

matrix_benchmark <- suppressWarnings(
  microbenchmark(
    `crossprod(X, y)` = crossprod(X, y),
    `t(X) %*% y` = t(X) %*% y,
    times = 200L,
    unit = "ms",
    control = list(order = "random", warmup = 10L)
  )
)

matrix_summary <- summary(matrix_benchmark, unit = "ms")
matrix_timing <- data.frame(
  method = as.character(matrix_summary$expr),
  median_ms = matrix_summary$median,
  relative_time = matrix_summary$median / min(matrix_summary$median)
)
knitr::kable(
  matrix_timing,
  row.names = FALSE,
  digits = 3,
  col.names = c("方法", "耗时中位数（毫秒）", "相对耗时")
)

方法	耗时中位数（毫秒）	相对耗时
crossprod(X, y)	0.415	1.000
t(X) %*% y	0.664	1.599

crossprod(X, y) 和 t(X) %*% y 在数学上等价，前者则直接表达交叉乘积结构，并避免显式构造转置矩阵。这里的计时把构造 t(X) 的成本包含在第二种方法中；如果应用场景允许预先计算并多次复用转置矩阵，基准测试的边界就应相应改变。

可靠的性能比较至少应满足四项要求：先验证不同实现给出同一目标量；在相同输入和计算环境中进行比较；重复运行并报告中位数或其他稳健汇总；明确哪些数据准备和中间对象构造被计入时间。运行时间也不是唯一标准，数值误差、峰值内存、可扩展性和代码可维护性都可能改变最终选择。

3.10 进一步阅读

浮点表示和浮点运算的经典综述可参见 Goldberg (1991)；舍入误差的系统讨论可参见 Wilkinson (1963) 和 Higham (2002)。统计问题中的数值方法可参见 Gentle (2009) 与 Monahan (2011)。矩阵条件数、数值秩和稳定矩阵分解可参见 Golub and Van Loan (2013)。

3.11 本章小结

本章介绍了计算统计中的基本数值概念。双精度浮点数只能表示实数的有限子集，其他实数只能近似表示；舍入误差、上溢、下溢和灾难性抵消可能使数学上正确的表达式产生不可靠的计算结果；等价变换、对数尺度和专用数值函数可以改善求值过程。条件性刻画问题对输入扰动的敏感程度，数值稳定性刻画算法传播误差的方式，稳定算法不能消除病态问题固有的输入敏感性。中心化和标准化可以改善尺度，但其统计含义必须与模型、惩罚项和先验分布一并考虑。复杂度分析用于判断规模增长的影响，基准测试则用于评估给定环境中的实际性能；二者都不能取代对准确性、存储和可复现性的检查。

后续章节会反复使用这些原则。数值线性代数章节将讨论矩阵分解、条件数与数值秩；优化章节将讨论目标函数尺度、迭代误差与停止准则；Monte Carlo 和 MCMC 章节将进一步区分数值误差与模拟误差，并使用对数尺度计算概率和权重。

3.12 思考题

1. .Machine$double.eps 与单位舍入误差有什么关系？为什么它不能被简单理解为“计算结果允许出现的绝对误差”？
2. 在什么情形下可以直接用 == 比较数值？在什么情形下应采用带容许误差的比较？容许误差为什么应当考虑问题尺度？
3. sqrt(x + 1) - sqrt(x) 的两种计算形式在数学上完全等价，为什么有限精度下会得到不同结果？请用有效数字抵消解释。
4. 区分病态问题与不稳定算法。稳定算法用于病态问题时，结果为什么仍可能不准确？不稳定算法用于良态问题时又可能发生什么？
5. 为什么似然、后验核、重要性权重和 MCMC 接受概率通常在对数尺度上计算？取对数不能解决哪些问题？
6. 标准化为什么可能改善优化和矩阵计算？在带惩罚回归或贝叶斯模型中，为什么标准化又不仅是一个数值技巧？
7. 大 \(O\) 复杂度与 microbenchmark 得到的运行时间分别回答什么问题？为什么不能仅依据一次基准测试断言某种写法在所有平台和规模下都更快？

3.13 上机实验（Lab）

1. Lab 1：浮点间距与灾难性抵消。 逐步增大整数 \(k\)，考察 1 + 2^(-k) == 1 何时成立，并把结果与 .Machine$double.eps 联系起来。再对 \(k=0,1,\ldots,20\) 令 \(x=10^k\)，比较 sqrt(x + 1) - sqrt(x) 与共轭变换后的公式，画出两种结果相对差异随 \(k\) 的变化。
2. Lab 2：稳定的对数计算。 先明确空向量和含 NA 输入的处理约定，再完善本章的 log_sum_exp()，使其同时正确处理这些输入以及含 Inf 或全部元素均为 -Inf 的情形。分别对包含很大有限正数和绝对值很大的有限负数的向量比较直接算法与稳定算法，并编写单元测试说明函数的边界行为。
3. Lab 3：方差算法与平移。 生成一组固定噪声\(Z_i\overset{\mathrm{iid}}{\sim}N(0,1)\)，观测值记为 \(z_i\)，并以mean((z - mean(z))^2) 作为参考量；依次令 \(x_i=c+z_i\)，其中\(c=10^2,10^4,\ldots,10^{12}\)。比较原始矩公式、中心化公式和经过分母调整的 var() 相对于参考量的误差，并讨论误差为何随 \(c\) 改变。
4. Lab 4：条件数、扰动与标准化。 构造两个相关程度逐渐升高的变量，记录设计矩阵的条件数。对响应向量施加若干小扰动，比较最小二乘系数的相对变化；随后对自变量标准化，说明标准化能改善什么、不能解决什么。
5. Lab 5：可重复的性能比较。 使用 microbenchmark 比较手写循环、向量化表达式和专用函数在向量求和或矩阵运算中的耗时。至少考察五种数据规模，先验证结果的一致性，再报告耗时中位数和相对耗时，并绘制对数耗时关于对数规模的图形。讨论为什么单一规模下的速度比较不能代替复杂度分析，以及临时对象和硬件环境如何影响结论。

参考文献

Gentle, James E. 2009. Computational Statistics. New York: Springer.

Goldberg, David. 1991. “What Every Computer Scientist Should Know about Floating-Point Arithmetic.” ACM Computing Surveys 23 (1): 5–48. https://doi.org/10.1145/103162.103163.

Golub, Gene H., and Charles F. Van Loan. 2013. Matrix Computations. 4th ed. Baltimore: Johns Hopkins University Press.

Higham, Nicholas J. 2002. Accuracy and Stability of Numerical Algorithms. 2nd ed. Philadelphia: SIAM.

Monahan, John F. 2011. Numerical Methods of Statistics. 2nd ed. Cambridge: Cambridge University Press.

Wilkinson, J. H. 1963. Rounding Errors in Algebraic Processes. London: Her Majesty’s Stationery Office.