2 误差
2.1 误差的种类
统计计算的算法要得到正确的结果， 就需要尽可能减少误差。 统计问题中的误差有模型误差、实验误差和数值计算误差， 在统计计算研究中主要解决的是如何减少数值计算误差的问题。
统计计算的算法通常是用来求解某种统计模型。 任何用来解决实际问题的数学模型都或多或少地简化了实际问题， 忽略掉一些细节，从而模型误差不可避免。 如果模型不合适，其它误差控制得再完美， 问题也不能得到解决； 更糟的是， 良好的计算结果会给使用者以错误的信心。 比如，我们使用的回归模型要求观测是独立的， 而实际数据观测有不可忽略的序列相关性， 尽管我们用软件算出了很完美的结果， 这个结果也是错误的。 我们应当仔细选择模型， 尽可能减少模型误差。
建立统计模型所需的数据来自实验、观测、抽样调查等过程， 在这样的过程中会出现实验误差， 包括随机误差、系统误差、过失误差。
随机误差是试验过程中由一系列随机因素引起的不易控制的误差， 可以通过多次重复试验或改进模型设计来减小随机误差。 随机误差可能来自物理量本身的波动， 比如测量风速，就是在测量一个随时变化的物理量， 不可避免地受到随机误差影响。 随机误差可能来自不可控制的随机因素影响， 比如，在用雷达测量飞机的方位和速度时， 可能受到地磁、气温、地形的影响。 由于测量仪器精度的限制也会产生随机误差， 比如用最小刻度是1度的温度计测量温度， 只能把不足1度的值四舍五入或者估计小数点后一位数字。 随机误差也可能来自特定条件下才发生的程序错误。
系统误差是多次测量持续偏高或偏低的误差， 多次重复测量不能消除或减少系统误差。 系统误差可能来自仪器本身的误差， 比如用不锈钢直尺测量家具高度， 直尺本身在温度不同时长度有细微变化。 系统误差也可能来自仪器使用不当， 比如用天平测量质量时天平没有配准。 当发现有系统误差时， 必须找出引起误差的原因并消除。
在记录实验数据时由于人的过失可以导致误差发生, 这样的误差称为过失误差。 比如，在记录仪表（如水表、电表）的读数时看错数字， 在记录数值时写错小数点位置， 在上传数据时报告了过时的或错误的数据，等等。 统计数据分析必须甄别并改正这样的过失误差， 否则会对分析结果产生严重影响。 在使用计算机软件处理数据时程序设计的质量问题也会导致误差发生。 比如，当输入条件不满足模型要求而程序中未进行检查时， 可能给出错误的结果。
2.2 数值误差
数值误差是用电子计算机进行数据存储和计算时产生的误差， 设计算法时必须了解并尽可能避免这种误差。
设\(A\)为结果的精确值，\(a\)为\(A\)在计算机中的近似值，则 \(\Delta = a - A\)称为绝对误差, 简称误差。 \(\delta = \frac{a-A}{A}= \frac{\Delta}{A}\) 称为相对误差。 相对误差没有量纲，常常用百分数表示。 绝对误差和相对误差常常仅考虑其绝对值。 实际中如果能估计绝对误差和相对误差的取值区间或绝对值的上限， 则可以确定误差大小。 如果绝对误差\(\Delta\)的绝对值上限估计为\(\tilde \Delta\), 则相对误差的绝对值上限可以用\(\delta = \frac{\tilde\Delta}{a}\)估计， 因为\(A\)的真实值一般是未知的。 有\(p\)位有效数字的一个十进制近似数的相对误差控制在\(5 \times 10^{-p}\)以下。
没有数值计算经验的读者往往会有一个错误认识， 认为计算机得到的结果总是准确的、可信的。 即使不考虑模型误差、随机误差和系统误差的影响， 用计算机进行数值计算也不可避免地受到误差影响， 只要我们把误差控制在可接受的范围就可以了。
2.2.1 计算机二进制
为了解数值计算误差， 我们首先需要了解计算机中数值的存储方式。 计算机中的数都用二进制表示，包括定点表示和浮点表示两种。
2.2.1.1 定点表示
定点数主要用于保存整数，是精确表示的， 定点数的四则运算可以得到精确结果， 但整数除法需要特别注意余数问题， 另外定点表示的整数范围有限， 比如用32个二进制位可以表示\(-2^{31} \sim 2^{31}-1\)之间的整数（约10位有效数字）， 定点数的计算结果有可能溢出， 即超出能表示的整数范围。
例2.1 (Julia中的定点整数(*)) 显示Julia语言中整数的界限与二进制表示。
整数1的表示：
x=1
typeof(x)
## Int64
bitstring(x)
## "0000000000000000000000000000000000000000000000000000000000000001"上面的x是64位定点整数，首位的0是符号位，表示非负，
后面有62个0和1个1。
上面的x是8位定点整数，
首位的0是符号位，表示非负，后面有6个0和1个1。
上面的x是8位无符号定点整数，有7个0和1个1, 没有符号位。
上面的x是8位定点整数，首位1表示负号，后面有7个1。
Julia中各种整数类型的最大可表示值：
Int8约有2-3位十进制有效数字。
Int16约有4-5位十进制有效数字。
Int32约有9-10位十进制有效数字。
Int64约有18-19位十进制有效数字。
Int64的加法超出存储界限时发生溢出， 而且不提供错误信息：
typemax(Int64) + 1
## -9223372036854775808
bitstring(typemax(Int64))
## "0111111111111111111111111111111111111111111111111111111111111111"
bitstring(typemax(Int64)+1)
## "1000000000000000000000000000000000000000000000000000000000000000"变成了负数。
但是Int8, Int16, Int32可以自动提升到更大储存的类型:
注意Julia中整数加法、减法、乘法结果还是整数， 除法按实数除法计算， 整数之间的乘法也是整数。 两个整数之间进行乘方计算时结果也是整数， 这是Julia语言的特殊规定， 其它语言中一般不这样规定， R语言中的四则运算和乘方运算基本都是按浮点数计算。
Julia程序：
溢出了。写成实数的乘方：
※※※※※
2.2.1.2 浮点表示
数值计算时主要使用浮点表示, 比如用64个二进制位能表示绝对值在\(10^{-308} \sim 10^{308}\)之间的实数， 有效数字大约有\(16 \sim 17\)位。 二进制浮点数的表示包括\((S,E,F)\)三个部分， 即正负号、指数部分、小数部分。 小数部分的位数\(d\)决定了精度。一般存储\(2^{-1} \leq F < 1\), 这样\(F\)中首位一般是1，有些计算机不保存这一位。 例如，二进制数\(101.101_2=0.101101_2 \times 2^3\) 可以表示为\((+, +11_2, .101101_2)\)。
不同的计算机硬件可能采用不同的浮点数表示，
在现代常用的IEEE 64位浮点数的表示中，
一个浮点数\(x\)被表示成
\[
x = \pm 2^{q - 1023} \times (1.b_1 b_2 \dots b_{52})
\]
其中符号\(\pm\)用64位的首位表示，
随后11位存储指数部分，
看成一个11位无符号定点整数\(q\)，
指数部分\(E\)代表\(\times 2^{q-1023}\)，
但是\(q=0\)和\(q=2^{11}-1\)有特殊解释。
最后的52位是二进制小数部分\(F\)，
约定小数点前面总是1并且不保存这个1。
当\(q=0\)时，小数部分约定小数点前面是0，
保存所有的二进制小数位，
这时小数部分的有效位数比正常情况少。
\(q > 0\)时绝对值最小的浮点数约为2.2E-308。
当\(q=2^{11}-1\)(所有的11位指数位都为1)时，
代表某些不正常的实数，
这时如果小数位都是零，
用来代表Inf(无穷大);
如果小数位有1，
用来代表NaN(比如0.0/0.0的结果)。
\(q < 2^{11}-1\)时绝对值最大的浮点数约为1.80E308。
由于\(E\)和\(F\)两部分的位数限制， 二进制浮点数只有有限个， 这些数的集合记为\(\mathcal F\)， 在这个集合中可以制定一个次序， \((S,E,F)\)的下一个数是\((S, E, F+2^{-d})\)。 对于实数域\(\mathbb R\)中的实数\(x\)， 当\(|x|\)超出浮点数范围时，我们不能表示； 当\(|x|\)没有超出浮点数范围时， 一般也只能用\(\mathcal F\)中的数来近似表示。 比如，十进制数\(0.1\)如果表示成二进制数， 是一个无限循环小数\(0.000\dot{1}10\dot{0}_2\), 只能存储其中有限位。 在近似表示时对于不能保存的位数计算机中可能用了舍入法或者截断法处理, 这样造成的误差叫做舍入误差。
例2.2 (Julia中的浮点数的二进制表示(*)) 用Julia语言演示64位浮点数\(0.1\)的二进制存储。
这是Julia语言中Float64类型的浮点数，按IEEE 64位浮点数格式存储。 首位0是符号位， 后面11位是指数位，存储\(q\)的无符号整数值，恰好等于1019， 表示\(2^{1019 - 1023} = 2^{-4}\)， 最后52位是小数位， 都在二进制小数点后面， 而小数点前面有一个1， 所以近似表示\(1.\dot 100 \dot 1\)， 因此表示为 \[\begin{aligned} 0.1 \approx & + 2^{-4} \times \left[ 1 + (1 \times 2^{-1} + 0 \times 2^{-2} + 0 \times 2^{-3} + 1 \times 2^{-4}) \right. \\ & \left. + (1 \times 2^{-5} + 0 \times 2^{-6} + 0 \times 2^{-7} + 1 \times 2^{-8}) + \cdots \right] . \end{aligned}\]
※※※※※
设实数\(z\)在浮点数可表示范围内， 将其浮点表示记为\(\mbox{fl}(z)\)， 则绝对舍入误差满足 \[\begin{aligned} |\mbox{fl}(z) - z| \leq U |z|, \forall z. \end{aligned}\] 其中\(U\)称为机器单位(machine unit)， \(\mbox{fl}\)用截断近似时\(U=2^{-(d-1)}\), \(\mbox{fl}\)用四舍五入时\(U = 2^{-d}\)。 一个浮点数\(z\)用浮点表示\(\mbox{fl}(z)\)后误差范围如下 \[\begin{align*} \mbox{fl}(z) = z(1+u), \ |u| \leq U. \end{align*}\]
和\(U\)类似的一个数是机器\(\varepsilon_{\text{m}}\)， \(\varepsilon_{\text{m}}\)是使得\(1+\varepsilon_{\text{m}}\) 和\(1\)的表示不相同的最小正浮点数。 可以认为\(U \approx \varepsilon_{\text{m}}/2\)。 典型的双精度计算\(\varepsilon_{\text{m}}\)约为\(10^{-16}\)数量级， 而单精度的相应值与双精度相差约\(10^9\)倍, 即单精度计算与双精度计算的精度相差约9位有效数字。 使用双精度实数可以减少表示误差和计算误差， 但是不能消除这两种误差， 所以不能盲目相信使用了双精度后就总是能准确地进行数值计算。
在R软件中用变量.Machine$double.eps保存双精度计算的机器\(\varepsilon_{\text{m}}\)值。
在Julia语言中用eps(Float64)获得双精度浮点数的机器\(\varepsilon_{\text{m}}\)值。
2.2.2 计算误差
由于计算机存储的有限性， 浮点数的四则运算不符合结合律，计算结果与计算次序有关。 不同的算法可能导致不同的计算误差， 应该尽可能选用计算精度高的数学公式， 另外在设计算法时需要注意避免一些损失精度的做法。
例2.3 (结合律的反例) 浮点数的四则运算不符合结合律，计算结果与计算次序有关。
\(1.1 + 0.1 - 1.2\)的精确结果是0， 但是在R中计算，得
而
两个计算次序不同，结果都不等于零而且不同计算次序结果不同。 Julia计算结果与R的结果相同。
※※※※※
例2.4 (累加的简化例子) 计算 \[ \sum_{n=1}^{999} \frac{1}{n(n+1)} \]
可以直接累加计算，造成很大的累积误差。 只要把公式变换成 \[\begin{aligned} \sum_{n=1}^{999} \frac{1}{n(n+1)} = \sum_{n=1}^{999}\left( \frac{1}{n} - \frac{1}{n+1} \right) = 1 - \frac{1}{1000} = 0.999 \end{aligned} \] 就只要计算一个除法和一个减法。
nmax <- 999
exact <- 1 - 1/(nmax + 1)
direct <- 0.0
for(n in 1:nmax){
direct <- direct + 1/(n*(n+1))
}
direct
## [1] 0.999
direct - 0.999
## [1] 6.661338e-16
exact - 0.999
## [1] 0在R中计算， 后一算法与精确值0.999的差等于0， 前一算法与精确值0.999的差等于6.66133814775094e-16。
Julia中计算：
let
nn = 999
s = 0.0
for n=1:nn
s += 1.0/(n*(n+1))
end
@show s;
@show s - 0.999;
end;
## s = 0.9990000000000007
## s - 0.999 = 6.661338147750939e-16累加计算的结果有\(10^{-16}\)级别的误差。
※※※※※
在进行浮点数加减时，两个绝对值相差很大的数的加减严重损失精度。 设\(|a| \gg |b|\)(\(\gg\)表示“远大于”), 则\(a+b\)会被舍入为\(a\)。 比如，计算 \[\begin{align*} 0.1234 + 0.00001234 \end{align*}\] 如果仅允许保留4位有效数字，则结果为\(0.1234\)。 为避免这样的问题，应该只对大小相近的数相加减， 比如，可以把小的数组合后分别相加，再加起来。
两个相近数相减损失有效数字。如: \[\begin{align*} 0.8522 - 0.8511 = 0.0011 \end{align*}\] 有效数字个数从4位减低到2位。 统计中如下的方差计算公式： \[\begin{align*} \sum_{i=1}^n (x_i - \bar x)^2 = \sum_{i=1}^n x_i^2 - \frac{1}{n} (\sum_{i=1}^n x_i)^2 \end{align*}\] 就存在这样的问题。 在对分布函数\(F(x)\)计算右尾部概率\(1 - F(x)\)时， 如果\(x\)很大则\(F(x)\)很接近1，计算\(1 - F(x)\)会造成结果有效数字损失。 为此，统计软件中计算分布函数时常常可以直接计算右尾概率， 以及概率的对数。
例2.5 逻辑斯蒂分布函数为\(F(x) = \frac{1}{1 + e^{-x}}\), 当\(x\)很大时如果直接计算 \(1 - F(x) = 1 - \frac{1}{1 + e^{-x}}\)，就会出现有效位数损失； 只要改用 \[\begin{aligned} 1 - F(x) = \frac{1}{1 + e^x} \end{aligned}\] 计算，就可以避免两个相近数相减，提高精度。
例如， 用\(1 - \frac{1}{1 + e^{-x}}\)计算\(1 - F(8)\)结果为:
用\(\frac{1}{1 + e^x}\)计算\(1 - F(8)\)结果为
第一种计算公式损失了4位以上的有效数字:
cat(sprintf("%20.16E", y1), "\n", sprintf("%20.16E", y2), "\n", sep="")
## 3.3535013046637197E-04
## 3.3535013046647811E-04用Julia语言进行上述计算：
2-element Array{String,1}:
"0.00033535013046637197"
"0.0003353501304664781" 第一种计算公式损失了4位以上的有效数字。
※※※※※
例2.6 计算\(\sqrt{x^2 + 1} - |x|\)， 直接计算在\(|x|\)很大时会有较大误差， 改成\(1/(\sqrt{x^2 + 1} + |x|)\)则不损失精度。
例2.7 (数值微分中的浮点运算误差(*)) 如果某函数自变量相对变动在机器精度U附近时函数值也有变化， 则函数值的计算结果被舍入误差严重影响，不能得到有效计算结果。
对一元可微函数\(f(x)\), \[f'(x) = \lim_{h\to 0} \frac{f(x+h) - f(x)}{h}\] 当\(h\)很小时可以用\(\frac{f(x+h) - f(x)}{h}\)近似计算\(f'(x)\)。
取\(f(x) = e^x\), 用这种方法近似计算\(x=0\)处\(f'(x)\)的值。真值等于1。
令\(g(h) = |\frac{e^h - 1}{h} - 1|\)为近似计算的误差。 用浮点数计算，分别取\(h=10^{-1}, 10^{-2}, \dots, 10^{-20}\)。
hvec <- 10^seq(-1, -20, by=-1)
yp <- abs((exp(hvec) - 1) / hvec - 1.0)
plot(hvec, yp,
xlab="h", main="数值方法估计的导数的误差",
ylab="误差", log="xy")
abline(h=0, col="gray", lty=3)
cbind(hvec, yp)
## hvec yp
## [1,] 1e-01 5.170918e-02
## [2,] 1e-02 5.016708e-03
## [3,] 1e-03 5.001667e-04
## [4,] 1e-04 5.000167e-05
## [5,] 1e-05 5.000007e-06
## [6,] 1e-06 4.999622e-07
## [7,] 1e-07 4.943368e-08
## [8,] 1e-08 6.077471e-09
## [9,] 1e-09 8.274037e-08
## [10,] 1e-10 8.274037e-08
## [11,] 1e-11 8.274037e-08
## [12,] 1e-12 8.890058e-05
## [13,] 1e-13 7.992778e-04
## [14,] 1e-14 7.992778e-04
## [15,] 1e-15 1.102230e-01
## [16,] 1e-16 1.000000e+00
## [17,] 1e-17 1.000000e+00
## [18,] 1e-18 1.000000e+00
## [19,] 1e-19 1.000000e+00
## [20,] 1e-20 1.000000e+00发现\(g(0.1) = 0.052\)有较大误差，逐步减小\(h\)的值， 同时误差变小， \(g(10^{-8}) = 6.1 \times 10^{-9}\)， 这时误差已经比较小。 为了进一步减小误差， 取更小的\(h\)值，发现并不能减小误差， \(g(10^{-11}) = 8.3\times 10^{-8}\)与\(h=10^{-8}\)时误差相近。 继续减小\(h\)的值，发现\(g(10^{-15}) = 0.11\)， 误差反而变得很大， 这是因为\(\exp(10^{-15}) - 1\)是两个十分接近的数相减， 有效位数损失很多，除以\(10^{-15}\)又放大了误差。 进一步减小\(h\)，\(g(10^{-16})=1\), 这是因为\(\exp(10^{-16}) - 1\)的浮点结果等于0。
此例的Julia语言版本：
20×2 Matrix{Float64}:
0.1 0.0517092
0.01 0.00501671
0.001 0.000500167
0.0001 5.00017e-5
1.0e-5 5.00001e-6
1.0e-6 4.99962e-7
1.0e-7 4.94337e-8
1.0e-8 6.07747e-9
1.0e-9 8.27404e-8
1.0e-10 8.27404e-8
1.0e-11 8.27404e-8
1.0e-12 8.89006e-5
1.0e-13 0.000799278
1.0e-14 0.000799278
1.0e-15 0.110223
1.0e-16 1.0
1.0e-17 1.0
1.0e-18 1.0
1.0e-19 1.0
1.0e-20 1.0※※※※※
例2.8 误差的累积会放大误差。
把0.1累加一千万次，与真实值1e6比较。
Julia语言编程：
误差约\(10^{-4}\)。
注意双精度数\(\varepsilon_m\)在\(10^{-16}\)左右，
1e6的表示误差在\(10^{-10}\)左右，
这个累加误差比表示误差大了6个数量级。
※※※※※
浮点数做乘法运算时，结果有效位主要取决于精度较低的数。 做除法运算时，如果分母绝对值很小则会产生较大的绝对误差。
在比较两个浮点数是否相等时， 因为浮点数表示和计算的不精确性， 程序中不应该直接判断两个浮点数是否相等， 而应该判断两者差的绝对值是否足够小。
设计计算机数值计算的算法时一定要注意浮点算法的局限。 要了解能保存的最大和最小实数。例如，浮点数的范围有时不够用。 超过最大值会发生向上溢出, 绝对值小于最小正浮点数会发生向下溢出。 向上溢出一般作为错误；向下溢出经常作为结果零。 还要了解相对误差的范围， 用机器精度\(U\)或\(\varepsilon_m\)表示。
例2.9 (计算同生日概率) 设一个班有\(n\)个人， 则至少有两人的生日的月、日重合的概率为 \[\begin{aligned} p_n = 1 - \frac{P_{365}^n}{365^n} \end{aligned}\] 其中\(P_{365}^n = 365\times 364 \times\cdots \times (365-n+1)\)。
如果直接计算分式的分子和分母， 则当\(n\)达到121时分母的计算溢出。
只要改用如下计算次序就可以解决这个问题： \[\begin{aligned} P_n = 1 - \prod_{j=0}^{n-1} \frac{365-j}{365} \end{aligned}\]
days <- 365
nvec <- 1:days
pvec <- numeric(days)
pvec[1] <- 0
for(n in 2:days){
pvec[n] <- 1 - prod((365 - (0:(n-1)))/365)
}
plot(nvec, pvec, main="n个人中存在生日同月日的概率",
pch=16, cex=0.2,
xlab="人数", ylab="概率", ylim=c(0,1))
事实上，以上的对\(n\)的循环还可以用cumprod()函数简化为：
Julia语言版本：
days = 365
nvec = 1:days
pvec = zeros(days)
pcum = 1
for n=1:days
pcum *= (365 - (n-1))/365
pvec[n] = 1 - pcum
end
pvec注意Julia语言的循环不损失效率， 不需要费尽心思改成向量化的程序风格； 而R如果用循环则会慢一到二个数量级。
Julia也支持“广播”（类似于向量化运算）和cumsum，上述Julia程序也可改写为：
※※※※※
例2.10 (softmax函数的稳定性问题(*)) 在机器学习中常用如下的softmax函数： \[ \text{softmax}(x_1, \dots, x_n) = \left(\frac{e^{x_i}}{\sum_{j=1}^n e^{x_j}}, \ i=1,2,\dots,n \right) . \] 讨论其由于数值范围限制造成的误差。
当各自变量\(x_i\)取值范围适中时， 这个函数没有什么问题。 但是当某个\(x_i\)取值很大时， \(e^{x_i}\)会发生向上溢出错误； 当所有\(x_i\)都是绝对值很大的负数时， \(e^{x_j}\)会向下溢出变成0或者只有很少的有效位数， 使得分母等于零或接近于零， 造成除零错误或者产生很大相对误差。 为此，将公式修改为 \[ \text{softmax}(x_1, \dots, x_n) = \left(\frac{\exp\{x_i - \max_{k} x_k\}}{\sum_{j=1}^n \exp\{x_j - \max_{k} x_k\}}, \ i=1,2,\dots,n \right), \] 因为\(x_i - \max_{k} x_k \leq 0\)， 所以不会向上溢出； 因为分母中恰好有一个\(e^0 = 1\)， 所以不会出现分母绝对值很小导致大的相对误差的情形。
※※※※※
2.2.3 截断误差
除了数值计算中浮点数表示和计算带来的误差， 另一类更大的误差是截断误差。 比如，计算机中的超越函数如\(e^x\)通常用级数逼近， 而级数只能计算到有限项， 舍去部分造成误差。 在确定需要计算的项数时， 机器精度\(U\)或\(\varepsilon_\text{m}\)给出了一个精度的界限， 任何使得结果相对变化小于\(U\)或\(\varepsilon_\text{m}\)的尾项都不需要再加进来。 数值计算中方程求根、求最小值点经常使用迭代法求解， 而迭代法也不能无限执行， 一般预定一个精度，达到精度时计算停止， 会造成误差。 所以，即使是编程语言内置的函数或公认的程序包的结果也是有数值计算误差的, 不一定总能达到机器精度\(U\)， 见例2.7。
2.3 计算复杂度
例2.11 计算多项式 \[\begin{aligned} P_n(x) = a_0 + a_1 x_1 + \dots + a_n x^n \end{aligned}\]
若直接计算，\(a_k x^k\)需要计算\(k\)次乘法， 总共需要计算\(1 + 2 + \dots + n = \frac12 n(n+1)\)次乘法和\(n\)次加法, 我们称这样的算法需要\(O(n^2)\)次浮点运算。
如果用如下递推算法（秦九韶法）：
| 秦九韶法 |
|---|
| \(u \leftarrow a_n\) |
| for (\(k\) in \(1:n\)) { |
| \(\qquad\) \(u \leftarrow u \cdot x + a_{n-k}\) |
| } |
| 输出\(u\) |
就只需要\(n\)次乘法和\(n\)次加法, 即\(O(n)\)次浮点运算， 在提高了效率的同时也会提高计算精度。
例2.11展示了不同算法的运算次数可以有重大的差别。 一个算法的计算开销主要是计算时间长短和占用的存储空间大小。 计算时间取决于算法的运算次数， 计算机中各种运算的种类较多， 数值计算的算法主要计算浮点运算次数(Floating Point Operations, FLOPs)。 浮点运算次数将算法中每一次四则运算计为一次浮点运算。 浮点运算次数一般表示为\(O(n)\), \(O(n^2)\)等随问题规模\(n\)而变化的形式， 最好能精确到具体的比例系数，如\(O(\frac43 n^3)\)。
按照数学分析中的记号习惯， 称浮点运算次数\(T_n = O(n^k)\)，\(k>0\)， 是指存在\(C>0\)使得当\(n\)充分大时\(T_n \leq C n^k\)。 例如， 例2.11的直接计算方法的浮点运算次数为 \(\frac12 n(n+1) + n = \frac12 n^2 + \frac32 n\)， 令\(C = \frac32\), 则\(n\geq 1\)时\(\frac12 n(n+1) + n \leq C n^2\)， 所以可以称此算法的浮点运算次数为\(O(n^2)\)。 若\(\lim_{n\to\infty} \frac{T_n}{C n^k} = 1\)， 称\(T_ n \sim C n^k\)， 这时也称浮点运算次数为\(O(C n^k)\)， 比如例2.11的直接计算方法的浮点运算次数为\(O(\frac12 n^2)\)。 事实上，若\(T_n\)为\(n\)的\(k\)次多项式， 最高次项系数为\(C\)， 则\(T_n \sim C n^k\)， 而且对任意小的\(\delta > 0\)都有\(n\)充分大时\(T_n \leq (C+\delta) n^k\)， 所以将\(T_n \sim C n^k\)记作\(T_n = O(C n^k)\)是合理的。 称浮点次数为\(O(n^k)\)的算法具有\(k\)阶时间复杂度。
若\(\lim_{n\to\infty} \frac{T_n}{n^k} = 0\)，记作\(T_n = o(n^k)\)。 \(o(n^k)\)一定是\(O(n^k)\)。 \(T_n \sim Cn^k\)可以记作\(T_n = C n^k + o(n^k)\)。
2.4 随机误差的度量
随机误差来自于观测样本或随机模拟中的随机因素， 一般随着样本量增大而减小， 但随机误差是随机变量， 不能给出严格误差界限， 只能用概率方法描述随机误差大小。
为了方便讨论当样本量\(n\)趋于无穷时随机误差大小的变化规律， 引入如下的\(O_p\)和\(o_p\)的记号。
\(O_p(\cdot)\)是概率论中表示依概率有界的记号。 如果\(\{ \xi_n \}\)和\(\{ \eta_n \}\)是两个随机变量序列， 满足 \[\begin{aligned} \lim_{M \to \infty} \sup_{n \geq 1} P\left( \left| \frac{\xi_n}{\eta_n} \right | > M \right) = 0, \end{aligned}\] 则称\(\{ \frac{\xi_n}{\eta_n} \}\) 依概率有界， 记为\(\frac{\xi_n}{\eta_n} = O_p(1)\)或 \(\xi_n = O_p(\eta_n)\)。
如果对\(\forall \delta>0\)都有 \[\begin{align*} \lim_{n\to\infty} P\left( \left| \frac{\xi_n}{\eta_n} \right | > \delta \right) = 0 \end{align*}\] 称\(\xi_n / \eta_n\) 依概率趋于零， 记为\(\xi_n / \eta_n = o_p(1)\)或\(\xi_n = o_p(\eta_n)\)。
\(O_p\)和\(o_p\)是我们用来估计随机误差幅度的有效工具， 有如下性质:
- \(o_p(1)\) = \(O_p(1)\);
- 如果\(\{ \eta_n \}\)是趋于零的常数列， 则\(\xi_n = O_p(\eta_n)\)时必有\(\xi_n = o_p(1)\)。
- \(o_p(1) + o_p(1) = o_p(1)\);
- \(o_p(1) + O_p(1) = O_p(1)\);
- \(o_p(1) \cdot o_p(1) = o_p(1)\);
- \(o_p(1) \cdot O_p(1) = o_p(1)\);
- 如果\(\{ \xi_n \}\) 以概率1趋于零， 或\(\{ \xi_n \}\)是趋于零的非随机实数列， 则\(\xi_n = o_p(1)\);
- 单个随机变量\(\xi = O_p(1)\);
- 如果随机变量序列\(\xi_n\)依分布收敛到分布\(F\), 则\(\xi_n=O_p(1)\)， 且\(\xi_n + o_p(1)\)仍依分布收敛到分布\(F\)。
如果\(\xi_n\)满足中心极限定理 \[\begin{aligned} \frac{\xi_n - E \xi_n}{\sqrt{\text{Var}(\xi_n)}} \stackrel{d}{\rightarrow} \text{N}(0,1) \end{aligned}\] 设\(n \text{Var}(\xi_n)\)有上界\(D\)， 则 \[\begin{aligned} \xi_n - E \xi_n = \frac{\sqrt{n \text{Var}(\xi_n)}}{\sqrt{n}} \cdot O_p(1) = O_p\left( \frac{\sqrt{D}}{\sqrt{n}} \right) = o_p(1) . \end{aligned}\] 例如，设\(X_1, X_2, \dots\) 独立同分布，有共同的方差\(\sigma^2\), \(\bar X_n = \frac{1}{n} \sum_{i=1}^n X_i\), 则\(\bar X_n - EX_1 = O_p(\frac{\sigma}{\sqrt{n}})\)。
2.5 问题的适定性与算法稳定性
把算法粗略看成 \[\begin{aligned} \text{输出} = f(\text{输入}) \end{aligned}\] 问题的适定性可以看成是输入微小变化时输出变化的大小， 如果输出连续地依赖于输入并且输入的微小变化引发的输出变化也是微小的， 则问题是适定的。
为了问题的适定性，定义条件数如下： \[\begin{aligned} \frac{|f(\text{输入}+\delta) - \text{输出}|}{|\text{输出}|} = \text{条件数} \cdot \frac{|\delta|}{|\text{输入}|}. \end{aligned}\] 即条件数(condition index)是输出的相对变化与输入的相对变化的比值。 当\(f(x)\)为一元可微函数时， 计算\(f(x)\)的条件数可以用微分表示为 \[\begin{aligned} \kappa(f,x) = |x f'(x) / f(x)|. \end{aligned}\] 条件数较小，比如\(\kappa < 10\)的时候， 可以设计算法给出问题的精确解。 条件数很大的问题称为病态问题(ill-conditioned), 可能会有很大误差。 条件数等于无穷或不存在的问题称为不适定问题(ill-posed)。
即使问题是适定的， 同一个问题也可以有多种不同算法实现， 不同算法结果的精度可能有很大的差别， 只有设计良好的算法才能使得结果达到条件数代表的精度。 设计不当的算法可能因为浮点运算的不良影响使得结果精度远低于条件数所代表的精度， 这样的算法称为不稳定的算法。 对于适定问题， 应该设计稳定的算法使得浮点运算的不良影响受控， 使得结果能够基本达到条件数所代表的精度。
例2.12 (解一元二次方程的稳定性问题(*)) 考虑二次方程 \[\begin{aligned} z^2 - x_1 z + x_2 = 0, \ (x_1>0, x_2>0) . \end{aligned}\] 设其中\(x_2\)很小，这时两个根都是正根，求其中较小根\(z_2\)。
\(z_2\)定义为 \[\begin{align} z_2 = z_2(x_1, x_2) = \frac{x_1 - \sqrt{x_1^2 - 4 x_2}}{2} \tag{2.1} \end{align}\] 在\(x_2\)变化时条件数 \[\begin{aligned} C =& \left| x_2 \frac{ \partial z_2(x_1, x_2)}{\partial x_2} / z_2(x_1, x_2) \right| \\ =& x_2 \frac{1}{x_1^2 - 4 x_2} \frac{2}{x_1 - \sqrt{x_1^2 - 4 x_2}} = \frac{z_1}{z_1 - z_2} \end{aligned}\] 当\(z_1\)较大，\(z_2\)很小时\(C\)接近1， 问题适定性很好。 当判别式\(x_1^2 - 4 x_2\)很接近于零时两个根\(z_1\)和\(z_2\)差很小， 条件数很大， 这时\(x_2\)的一个微小变化可能引起根\(z_2\)的很大变化。
即使在\(C=1\)时不适当的算法也可能会造成结果不稳定。 设判别式\(x_1^2 - 4 x_2\)不接近于零， 两个根\(z_1\)和\(z_2\)差距很大， 但是当\(x_2\)很小时公式(2.1)的分子中\(x_1\)和\(\sqrt{x_1^2 - 4 x_2}\)很接近， 相减会造成精度损失， 改用公式 \[\begin{align} z_2 = \frac{2}{x_1 + \sqrt{x_1^2 - 4 x_2}} \tag{2.2} \end{align}\] 则避免了相近数相减， 提高了精度。
下面程序中x1>0, x2>0，方程的左边看成z, x1, x2的函数：
给定x1, x2, 按(2.1)计算z1(较大根)和z2(较小根)的函数为：
z1_1 <- function(x1, x2) (x1 + sqrt(x1^2 - 4*x2))/2
z2_1 <- function(x1, x2) (x1 - sqrt(x1^2 - 4*x2))/2给定x1, x2, 按(2.2)计算z2(较小根)的函数为：
给定x1, x2, 计算条件数的函数为：
取正常的x1=1，
取x2为很小的值，这里取\(0.5\varepsilon_m\)：
这时，第一种方法解出的z2为：
第二种方法的解出的z2为：
两个解的差为：
对x2增加一个扰动\(\delta\)：
这时两种方法解出的z2及其差为：
z2_1(x1, x2 + dx2)
## [1] 1.665335e-16
z2_2(x1, x2 + dx2)
## [1] 1.332268e-16
z2_1(x1, x2 + dx2) - z2_2(x1, x2 + dx2)
## [1] 3.330669e-17两种方法的结果差别很大。 增加扰动后的变化与条件数给出的相对误差界限比较， 条件数为：
条件数等于1，问题适定；
但是两种方法在输入的x2加了扰动之后的相对变化分别为：
(z2_1(x1, x2 + dx2) - z2_1(x1,x2))/z2_1(x1, x2)
## [1] 0.5
(z2_2(x1, x2 + dx2) - z2_2(x1,x2))/z2_2(x1, x2)
## [1] 0.2这些相对变化应该都近似等于\(C \delta / x_2 = 0.2\)，
但第一种方法的相对变化超出了应有的值，
其原因是(2.1)的分子中x2很小，
基本相当于x1 - x1，
是两个很相近的数相减，会严重损失精度。
此例计算的Julia语言版本：
f0(z, x1, x2) = z^2 - x1*z + x2 # 方程左边
# 较小的根，第一种公式，直接计算
z2v1(x1,x2) = (x1 - sqrt(x1^2 - 4*x2))/2.0 # 直接求较小根的公式
# 较小的根，第二种公式，避免减法损失精度
z2v2(x1,x2) = 2*x2/(x1 + sqrt(x1^2 - 4*x2)) # 改进的求较小根的公式
# 计算z2的条件数，在x2变化时z2的条件数
f1cond(x1,x2) = 2*x2/(x1^2 - 4*x2)/(x1 - sqrt(x1^2 - 4*x2))
f1cond(1.0, 0.5*eps(1.0)) # 取x1=1, x2=0.5ε, 求条件数，约为1
## 1.0000000000000004取x1=1, x2=0.5ε, 直接算较小根:
取x1=1, x2=0.6ε, 直接求较小根:
采用第一种公式时， x2变化0.1ε时第一种方法求较小根的估计的条件数，大于理论条件数：
取x1=1, x2=0.5ε, 用改进方法计算较小根：
取x1=1, x2=0.6ε, 用改进方法计算较小根：
x2变化0.1ε时用改进方法得到的结果估计条件数，等于理论值：
第二种算法得到的相对变化与条件数一致， 而第一种算法的结果则有过大的变化。 Julia的计算与R的结果类似但不相同。
※※※※※
稳定性研究经常用倒推法： 设精确结果应该为\(y=f(x)\), \(x\)是输入， \(y\)是输出。 有误差的结果记作\(y^*=f^*(x)\), 设\(y^*\)等于某个输入\(x^*\)的精确结果\(f(x^*)\)， 则\(\| x - x^* \|\)的大小代表了算法的稳定程度, \(\| x-x^* \|\)称为倒推误差， \(\| x - x^* \| / \| x \|\)称为倒推相对误差。 如果\(\| x-x^* \| / \| x \|\)很小则称算法是倒推稳定的。
例2.13 (浮点数加法的倒推误差) 考虑两个大小相近的浮点数加法的倒推稳定性。
令\(f(\boldsymbol x) = x_1 + x_2\), 则浮点运算结果为 \[ f^*(\boldsymbol x) = \mbox{fl}( \mbox{fl}(x_1) + \mbox{fl}(x_2) ) \] (这里不考虑加法运算本身的误差)， 于是 \[\begin{align*} f^*(x) =& x_1 + \delta_1 x_1 + x_2 + x_2 + \delta_2 x_2 + \delta_3 (x_1 + x_2) \\ =& \left[ x_1 + (\delta_1 + \delta_3) x_1 \right] + \left[ x_2 + (\delta_2 + \delta_3) x_2 \right] \\ =& f\left( (x_1 + (\delta_1 + \delta_3) x_1, \; x_2 + (\delta_2 + \delta_3) x_2 )^T \right) \end{align*}\] 其中\(|\delta_i| \leq U\), \(i=1,2,3\)。 倒推相对误差为 \[\begin{align*} \frac{\| \boldsymbol x^* - \boldsymbol x \|}{\| \boldsymbol x \|} =& \frac{\| ( (\delta_1 + \delta_3) x_1, \; (\delta_2 + \delta_3) x_2 )^T \|}{\| \boldsymbol x \|} \leq 2U = \varepsilon_{\text{m}} . \end{align*}\]
※※※※※
习题
习题2
设某正数\(a\)四舍五入后保留了\(p\)位有效数字， 表示成\(a^* = (0. a_1a_2 a_3 \dots a_p)_{10} \times 10^m\) (\(a_1 \neq 0\))。 估计其绝对误差和相对误差的范围， 并分析误差与有效数字位数\(p\)的关系。
习题3
如下级数 \[\begin{aligned} S_n = 1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \dots + (-1)^{n+1} \frac{1}{n} \end{aligned}\] 收敛到\(\ln 2\)。 用下列四种不同方法计算\(n=10^6\)时级数值并比较精度：
- 按正常次序相加；
- 按相反次序相加；
- 成对组合后相加：\((1 - \frac{1}{2}) + (\frac{1}{3} - \frac{1}{4}) + \dots\)；
- 按成对组合的通分后结果相加：\(\frac{1}{1\times 2} + \frac{1}{3\times 4} + \dots\)。
习题4
对函数\(e^{-x}\)，可以用两种公式计算： \[\begin{aligned} (1)&\quad e^{-x} = 1 - x + x^2/2 - x^3/3! + \dots + (-1)^k x^k / k! + \dots;\\ (2)&\quad e^{-x} = 1 / (1 + x + x^2/2 + x^3/3 + \dots + x^k/k! + \dots) \end{aligned}\] 对\(x=1,2,3,4\)，计算到\(k=10\)的级数并比较两种算法的计算精度。
习题5
有如下一组观测数据： \[\begin{aligned} & 249,\ 254,\ 243,\ 268,\ 253,\ 269,\ 287,\ 241,\ 273,\ 306 \\ & 303,\ 280,\ 260,\ 256,\ 278,\ 344,\ 304,\ 283,\ 310 \end{aligned}\] 用两种方法计算样本方差：
- 公式 \[\begin{aligned} S^2 = \frac{1}{n-1}\left( \sum_{i=1}^n x_i^2 - n \bar x^2 \right) \end{aligned}\]
- 公式 \[\begin{aligned} S^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar x)^2 \end{aligned}\]
其中每一个中间步骤保留6位有效数字(R中用signif函数把数据四舍五入到指定有效位数)。
把结果与使用高精度计算的结果比较。