无偏估计、方差与 Monte Carlo 平均

前面一章讨论 QR algorithm、Hessenberg 形式和 SVD 时，主线是用结构化的正交变换保留矩阵的核心信息

这一章换成另一种思路：如果矩阵太大，完整计算或者完整分解都太贵，那么能不能只随机看一部分信息，再用统计意义上的平均来恢复整体？

这就是 Randomized Linear Algebra 的基本原理

当然，单次随机结果并不一定正确，我们要做的是构造一个 随机估计量，使它在 平均 意义上指向正确目标，同时通过增加 样本数 和合理选择 抽样概率 来控制 方差

从概率统计的角度看，这属于 Monte Carlo estimation 的思想：目标量本身可以是确定的，但完整计算它很贵，于是人为引入随机抽样，用随机样本的平均来估计这个目标量

在这一章里，目标量是确定的矩阵乘法

$$\mathbf{A}\mathbf{B}$$

我们把它写成很多外积项的总和，然后随机抽其中一部分来估计完整总和。这可以理解成把 Monte Carlo method 用在矩阵乘法上

为了避免后面的内容让人误解，先把几个统计学概念按依赖关系讲清楚

从真实目标到估计量

统计估计的出发点并非只局限于 “求一个平均值”，而是有一个我们想知道、但不想或不能直接完整计算的对象

这个对象叫 真实目标，记作

$$\theta$$

$\theta$ 可以是一个数，也可以是一个向量、一个矩阵，甚至是某个优化问题的最优参数。它不一定只能是均值等于我们想要的结果

比如在这一章后面的随机矩阵乘法里，真实目标就是

$$\theta=\mathbf{A}\mathbf{B}$$

因为完整计算 $\mathbf{A}\mathbf{B}$ 可能很贵，所以我们通过随机抽样构造一个替代品。这个替代品叫 估计量 estimator

如果样本是

$$X_1,X_2,\dots,X_s$$

那么估计量一般写作

$$\hat\theta = f(X_1,X_2,\dots,X_s)$$

也就是说，估计量是样本的函数

它可以是样本均值

$$\bar X=\frac{1}{s}\sum_{i=1}^s X_i$$

也可以是样本方差、样本中位数、某个优化问题的解，或者一个矩阵

例如

$$\hat\sigma^2 = \frac{1}{s-1}\sum_{i=1}^s (X_i-\bar X)^2$$

或者

$$\hat\theta = \arg\min_\theta \sum_{i=1}^s L(X_i,\theta)$$

都可以是估计量

所以 “估计” 这个过程的基本结构是

$$\boxed{ \text{真实目标 } \theta \quad \longleftarrow \quad \text{用样本构造出的估计量 } \hat\theta }$$

估计量不是 “检测抽样分布的期望”。更准确地说，估计量是样本函数；而判断它是否无偏时，才去看它的抽样分布的期望

重点是：样本只是随机看到的一部分信息，估计量才是我们准备用来估计目标的最终表达式。它不一定要复杂地修改样本，但一定是在说明 “怎样用这些样本去代表目标”。

三种容易混淆的分布

第一种是 probability distribution。它描述抽样机制本身，也就是随机变量理论上怎样取值、各自以多大概率出现

如果随机变量 $X$ 可能取

$$x_1,x_2,\dots,x_n$$

并且概率为

$$p_1,p_2,\dots,p_n$$

其中

$$\sum_{j=1}^n p_j=1$$

那么这组 $p_j$ 就是 $X$ 的概率分布。它描述的是“还没有抽样之前”的规则，而不是某一次抽样之后看到的结果

第二种是 sample distribution，也常叫 empirical distribution。它描述某一次实际抽到的样本里，各种结果出现了多少次

比如抽了 $s$ 次，结果 $x_j$ 出现了 $c_j$ 次，那么经验频率是

$$\hat p_j=\frac{c_j}{s}$$

当 $s$ 很大时，如果样本真的来自概率分布 $p_j$，经验频率通常会接近理论概率

$$\hat p_j\approx p_j$$

第三种是 sampling distribution。它不是原始随机变量 $X$ 的分布，而是估计量 $\hat\theta$ 的分布

因为样本是随机的，所以

$$\hat\theta=f(X_1,X_2,\dots,X_s)$$

本身也是随机的。每重新抽一批样本，$\hat\theta$ 就会变一次。重复很多次之后，$\hat\theta$ 形成的分布，就是抽样分布

因此常见的 bell curve 图，很多时候画的不是原始数据 $X$ 的分布，而是估计量 $\hat\theta$ 在重复抽样下如何波动

可以把三者的区别记成

概率分布描述抽样规则；经验分布描述一次样本；抽样分布描述估计量的重复实验结果。

有偏估计和无偏估计

现在才轮到 bias 和 unbiasedness

因为 $\hat\theta$ 是随机的，所以它有自己的抽样分布。我们可以对这个抽样分布求期望

$$\mathbb{E}[\hat\theta]$$

如果这个期望刚好等于真实目标

$$\boxed{ \mathbb{E}[\hat\theta]=\theta }$$

那么 $\hat\theta$ 是 $\theta$ 的 unbiased estimator

换句话说，无偏估计并不是说每一次抽样都正好能得到正确答案，而是要求在反复抽样、反复计算估计量后，它的长期平均中心落在真实目标上

如果

$$\mathbb{E}[\hat\theta]\neq\theta$$

那么 $\hat\theta$ 是 biased estimator

偏差定义为

$$\operatorname{Bias}(\hat\theta) = \mathbb{E}[\hat\theta]-\theta$$

也就是估计量的抽样分布中心偏离真实目标的量

因此这几个概念的顺序应该是

$$\boxed{ \text{真实目标 } \theta \rightarrow \text{估计量 } \hat\theta \rightarrow \text{估计量的抽样分布} \rightarrow \mathbb{E}[\hat\theta] \rightarrow \text{有偏或无偏} }$$

这里的关键是：无偏不是说每一次估计都等于真实值，也不是说每个可能结果等概率，更不是说目标 $\theta$ 必须是某个均值。它只说估计量的抽样分布期望等于目标

不要把抽样概率和样本平均混在一起

很容易混淆的是下面两件事：

首先是某个结果本来有多容易被抽到。这个由 抽样概率 决定，即

$$p_j$$

它表示第 $j$ 个可能结果被抽到的概率。所有 $p_j$ 合在一起，描述的是抽样之前的 概率分布 probability distribution

其次是已经抽到了 $s$ 次以后，怎样把这 $s$ 条样本记录平均起来。这个由样本平均里的

$$\frac{1}{s}$$

决定。它表示每一条已经抽到的样本记录，在这次平均里占相同权重

所以 $p_j$ 和 $1/s$ 说的不是同一件事：

$$\boxed{ p_j \text{ 决定某个可能结果出现得多不多，} \qquad \frac{1}{s} \text{ 决定已经抽到的某条记录有多重要} }$$

例如样本平均写成

$$\bar X=\frac{1}{s}\sum_{i=1}^s X_i$$

这里的 $1/s$ 并不表示 $X$ 的所有可能取值等概率。它的意思是：现在手里有 $s$ 次抽样结果，要把这 $s$ 条记录平均起来

如果某个结果本来更容易出现，它会更频繁地出现在样本里。比如 $X$ 只可能取 $x_1$ 和 $x_2$，并且 $x_1$ 更容易被抽到。抽 $100$ 次以后，可能得到 $90$ 次 $x_1$ 和 $10$ 次 $x_2$

这时样本平均是

$$\bar X = \frac{1}{100} \left( 90x_1+10x_2 \right) = \frac{90}{100}x_1+\frac{10}{100}x_2$$

注意，这里的 $\frac{1}{100}$ 是每一条样本记录的权重；而 $\frac{90}{100}$ 和 $\frac{10}{100}$ 是这一次样本里 $x_1$ 和 $x_2$ 出现的经验频率

如果样本平均里不是每条记录都乘同一个 $\frac{1}{100}$，那就表示我们在平均时人为给某些样本记录更大的权重。普通样本均值不是这么做的；它只是让已经抽到的每一条记录平等参与平均。

一般地，样本均值可以写成

$$\bar X = \sum_j \hat p_j x_j$$

其中

$$\hat p_j = \frac{\text{样本中 }x_j\text{ 出现的次数}}{s}$$

是这一次样本形成的经验频率。它不是抽样规则本身，而是这一次抽样之后观察到的比例

当样本数足够大时，经验频率会接近理论抽样概率

$$\hat p_j\approx p_j$$

于是样本均值才会接近理论期望

$$\bar X \approx \sum_j p_j x_j = \mathbb{E}[X]$$

所以真正表示等概率抽样的是

$$p_1=p_2=\cdots=p_n$$

而不是样本平均公式里的 $1/s$

这个区别在后面的随机矩阵乘法里尤其重要。后面会出现类似

$$\widehat{\mathbf{A}\mathbf{B}} = \frac{1}{s} \sum_{\ell=1}^s \frac{1}{p_{j_\ell}}\mathbf{M}_{j_\ell}$$

的形式。其中 $1/s$ 是对 $s$ 次抽样结果取平均，$p_{j_\ell}$ 是第 $\ell$ 次抽到的那一项原本的抽样概率，而 $1/p_{j_\ell}$ 是为了修正抽样概率带来的偏向

为什么随机估计量要多次平均

现在回到估计量本身。即使某个随机估计量是无偏的，它的单次结果也通常会波动

也就是说，单次估计量 $Y$ 可能满足

$$\mathbb{E}[Y]=\theta$$

但这只说明它的长期平均中心是 $\theta$，不说明每一次 $Y$ 都很接近 $\theta$

为了让估计更稳定，通常会独立重复抽样，得到

$$Y_1,Y_2,\dots,Y_s$$

每个 $Y_i$ 都是对同一个目标 $\theta$ 的随机估计，并且满足

$$\mathbb{E}[Y_i]=\theta$$

然后取平均

$$\bar Y=\frac{1}{s}\sum_{i=1}^s Y_i$$

多次独立抽样并取平均，会把单次估计量 $Y$ 的抽样分布变成平均估计量 $\bar Y$ 的抽样分布。接下来要分别看两件事：

• 这个新估计量的期望中心是否仍然是 $\theta$，也就是它是否仍然是无偏估计
• 这个新估计量的方差是否比单次估计更小，也就是估计结果是否更稳定、更接近真实目标

对期望的影响

先回忆一下 $\mathbb{E}$ 本身的计算公式。对于离散随机变量 $Z$，如果它可能取值为

$$z_1,z_2,\dots,z_m$$

并且

$$\mathbb{P}(Z=z_r)=q_r$$

那么它的期望是

$$\boxed{ \mathbb{E}[Z] = \sum_{r=1}^m q_r z_r }$$

也就是按概率加权求和。因为概率权重满足 $\sum_{r=1}^m q_r=1$，所以它也常被称为按概率加权的平均，但这里更应该先把它理解成加权和

在这里，$\bar Y$ 也是一个随机变量。每一次重新抽样，都会得到一个可能不同的 $\bar Y$。所以严格来说，

$$\mathbb{E}[\bar Y] = \sum_{\text{所有可能的样本组合}} \mathbb{P}(\text{该样本组合}) \cdot \bar Y(\text{该样本组合})$$

只是这个写法太长，所以通常利用期望的线性性质来简化计算

这个平均不会破坏无偏性。因为期望是线性的，

$$\begin{align*} \mathbb{E}[\bar Y] &= \mathbb{E} \left[ \frac{1}{s}\sum_{i=1}^s Y_i \right] \\ &= \frac{1}{s}\sum_{i=1}^s \mathbb{E}[Y_i] \\ &= \theta \end{align*}$$

这里的关键是：

$$\mathbb{E}[cZ]=c\mathbb{E}[Z]$$

所以外面的 $1/s$ 只是以一次方拿出来。再加上每个 $Y_i$ 的期望都等于 $\theta$，于是

$$\begin{align*} \mathbb{E}[\bar Y] &= \mathbb{E} \left[ \frac{Y_1+Y_2+\cdots+Y_s}{s} \right] \\ &= \frac{1}{s} \left( \mathbb{E}[Y_1]+\mathbb{E}[Y_2]+\cdots+\mathbb{E}[Y_s] \right) \\ &= \frac{1}{s} \left( \theta+\theta+\cdots+\theta \right) \\ &= \frac{s\theta}{s} \\ &= \theta \end{align*}$$

因此，多次独立抽样并取平均不会改变估计量的期望中心：如果单次估计量 $Y_i$ 是无偏的，那么平均估计量 $\bar Y$ 仍然是无偏的。

对方差的影响

期望中心不变，并不代表波动不变。平均估计量的好处在于：它会让抽样分布更集中，也就是方差更小

如果每次估计相互独立，并且它们由同一个抽样机制产生，那么可以设

$$\operatorname{Var}(Y_i)=\sigma^2$$

现在从方差定义开始算。因为 $\mathbb{E}[\bar Y]=\theta$，所以

$$\operatorname{Var}(\bar Y) = \mathbb{E} \left[ (\bar Y-\theta)^2 \right]$$

而

$$\bar Y-\theta = \frac{1}{s} \sum_{i=1}^s (Y_i-\theta)$$

这一步的意思是：平均估计量偏离目标的量，等于每一次估计偏离目标的量的平均

具体来说，$\theta$ 也可以写成 $s$ 个 $\theta$ 的平均：

$$\theta = \frac{1}{s}\sum_{i=1}^s \theta$$

所以，由 $\bar Y$ 的定义开始展开：

$$\begin{align*} \bar Y-\theta &= \left( \frac{1}{s}\sum_{i=1}^s Y_i \right) - \theta \\ &= \frac{1}{s}\sum_{i=1}^s Y_i - \frac{1}{s}\sum_{i=1}^s \theta \\ &= \frac{1}{s} \left( \sum_{i=1}^s Y_i - \sum_{i=1}^s \theta \right) \\ &= \frac{1}{s} \sum_{i=1}^s (Y_i-\theta) \end{align*}$$

也就是说，平均估计量偏离目标的量，等于每一次估计偏离目标的量的平均

于是从方差定义继续：

$$\begin{align*} \operatorname{Var}(\bar Y) &= \mathbb{E} \left[ \left( \frac{1}{s} \sum_{i=1}^s (Y_i-\theta) \right)^2 \right] \\ &= \frac{1}{s^2} \mathbb{E} \left[ \left( \sum_{i=1}^s (Y_i-\theta) \right)^2 \right] \\ &= \frac{1}{s^2} \mathbb{E} \left[ \sum_{i=1}^s (Y_i-\theta)^2 + 2\sum_{1\leq i<k\leq s}(Y_i-\theta)(Y_k-\theta) \right] \\ &= \frac{1}{s^2} \left( \sum_{i=1}^s \mathbb{E}[(Y_i-\theta)^2] + 2\sum_{1\leq i<k\leq s} \mathbb{E}[(Y_i-\theta)(Y_k-\theta)] \right) \end{align*}$$

上面第三步是在展开一个平方和。这里每一项都是一个偏差，展开以后会出现两类项：平方项和交叉项

比如 $s=3$ 时，

$$\begin{align*} \left[ (Y_1-\theta)+(Y_2-\theta)+(Y_3-\theta) \right]^2 &= (Y_1-\theta)^2 + (Y_2-\theta)^2 + (Y_3-\theta)^2 \\ &\quad + 2(Y_1-\theta)(Y_2-\theta) \\ &\quad + 2(Y_1-\theta)(Y_3-\theta) + 2(Y_2-\theta)(Y_3-\theta) \end{align*}$$

所以一般情形可以写成

$$\left( \sum_{i=1}^s (Y_i-\theta) \right)^2 = \sum_{i=1}^s (Y_i-\theta)^2 + 2\sum_{1\leq i<k\leq s} (Y_i-\theta)(Y_k-\theta)$$

这里的

$$1\leq i<k\leq s$$

表示在 $1,2,\dots,s$ 这些下标里，取所有不同的下标对，并且每一对只算一次。比如 $s=4$ 时，它包括

$$(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)$$

不用 $i\neq k$ 是因为那样会把 $(1,2)$ 和 $(2,1)$ 都算进去，重复计算。这里用 $i<k$，再在前面乘 $2$，正好对应平方展开中成对出现的交叉项

接下来把 $\mathbb{E}$ 分配进去，是因为期望是线性的：

$$\mathbb{E}[A+B]=\mathbb{E}[A]+\mathbb{E}[B]$$

因此平方项和交叉项可以分别求期望，也就是

$$\frac{1}{s^2} \left( \sum_{i=1}^s \mathbb{E}[(Y_i-\theta)^2] + 2\sum_{1\leq i<k\leq s} \mathbb{E}[(Y_i-\theta)(Y_k-\theta)] \right)$$

现在分别处理这两类项。这里最好把两个条件分开看：

• 无偏性 只说明单个估计量的偏差平均为 $0$
• 独立性 只说明不同估计量之间不会一起系统性偏大或偏小

先看无偏性。因为每个 $Y_i$ 都是对 $\theta$ 的无偏估计，

$$\mathbb{E}[Y_i]=\theta$$

所以

$$\mathbb{E}[Y_i-\theta] = \mathbb{E}[Y_i]-\theta = 0$$

同理，

$$\mathbb{E}[Y_k-\theta] = 0$$

这一步只说明单个偏差项的长期平均为 $0$。它还不能单独推出交叉项

$$\mathbb{E}[(Y_i-\theta)(Y_k-\theta)]$$

为 $0$

接着才用独立性。当 $i\neq k$ 时，$Y_i$ 和 $Y_k$ 独立，所以 $Y_i-\theta$ 和 $Y_k-\theta$ 也独立。于是乘积的期望可以拆开：

$$\begin{align*} \mathbb{E}[(Y_i-\theta)(Y_k-\theta)] &= \mathbb{E}[Y_i-\theta]\mathbb{E}[Y_k-\theta] \\ &= 0\cdot 0 \\ &= 0 \end{align*}$$

所以交叉项会消失，是两个条件共同作用的结果：

无偏性给出 $\mathbb{E}[Y_i-\theta]=0$；独立性允许把乘积期望拆开。

直观上说，无偏性表示单个估计的正负偏差长期会抵消；独立性表示一个估计偏大或偏小，不会系统性地带动另一个也偏大或偏小。因此不同估计之间的偏差乘积，长期平均会抵消掉

交叉项消失以后，还剩平方项。注意，平方项不会因为无偏性而消失。无偏性说的是

$$\mathbb{E}[Y_i-\theta]=0$$

也就是偏差本身正负会互相抵消；但平方以后，

$$(Y_i-\theta)^2\geq 0$$

正负号被消掉了，所以它衡量的是“偏离目标的大小”，不会因为正负抵消而变成 $0$

一个最简单的例子是

$$Y_i-\theta = \begin{cases} 1, & \text{概率 } \frac12 \\ -1, & \text{概率 } \frac12 \end{cases}$$

此时

$$\mathbb{E}[Y_i-\theta] = \frac12\cdot 1+\frac12\cdot(-1) = 0$$

这表示正偏差和负偏差在平均意义下抵消了

但是平方以后，

$$(Y_i-\theta)^2 = \begin{cases} 1, & \text{概率 } \frac12 \\ 1, & \text{概率 } \frac12 \end{cases}$$

所以

$$\mathbb{E}[(Y_i-\theta)^2] = 1$$

这说明无偏不等于没有波动。无偏只说 偏差的平均 是 $0$，方差衡量的是偏差大小的平均，也就是波动幅度

因此

$$\mathbb{E}[(Y_i-\theta)^2] = \operatorname{Var}(Y_i) = \sigma^2$$

这里没有给 $\sigma$ 加下标，是因为默认每次抽样使用同一个机制，所以 $Y_1,Y_2,\dots,Y_s$ 同分布。也就是说，

$$\operatorname{Var}(Y_1) = \operatorname{Var}(Y_2) = \cdots = \operatorname{Var}(Y_s) = \sigma^2$$

把交叉项为 $0$ 和平方项为 $\sigma^2$ 代回去：

$$\begin{align*} \operatorname{Var}(\bar Y) &= \frac{1}{s^2} \left( \sum_{i=1}^s \mathbb{E}[(Y_i-\theta)^2] + 2\sum_{1\leq i<k\leq s} \mathbb{E}[(Y_i-\theta)(Y_k-\theta)] \right) \\ &= \frac{1}{s^2} \left( \sum_{i=1}^s \mathbb{E}[(Y_i-\theta)^2] + 0 \right)\\ &= \frac{1}{s^2} \sum_{i=1}^s \sigma^2 \\ &= \frac{1}{s^2}s\sigma^2 \\ &= \frac{\sigma^2}{s} \end{align*}$$

这里有两个层次的求和，不要混在一起：

• 外面的 $\sum_{i=1}^s$ 是在把 $Y_1,\dots,Y_s$ 这 $s$ 次独立估计各自贡献的方差加起来
• 每个 $\mathbb{E}[(Y_i-\theta)^2]$ 如果按定义展开，里面还会有一个对 $Y_i$ 所有可能取值的概率加权求和

所以

$$\sum_{i=1}^s \mathbb{E}[(Y_i-\theta)^2]$$

是在说：先对每个随机估计量 $Y_i$ 求自己的方差，再把这 $s$ 个方差加起来

如果每次估计的方差不同，就应该写成

$$\operatorname{Var}(Y_i)=\sigma_i^2$$

这时结论会变成

$$\operatorname{Var}(\bar Y) = \frac{1}{s^2} \sum_{i=1}^s \sigma_i^2$$

只有当所有 $\sigma_i^2$ 都相同，才可以化简成 $\sigma^2/s$

这里也可以从一个更短的规则理解。方差和期望不同的地方是：

$$\operatorname{Var}(cZ)=c^2\operatorname{Var}(Z)$$

所以外面的 $1/s$ 进入方差以后会变成 $1/s^2$

另外，因为 $Y_i$ 之间相互独立，

$$\operatorname{Var} \left( \sum_{i=1}^s Y_i \right) = \sum_{i=1}^s \operatorname{Var}(Y_i)$$

如果不独立，中间还会多出协方差项，方差就不能这样简单相加

所以方差按 $1/s$ 下降，而标准差按 $1/\sqrt{s}$ 下降

$$\operatorname{Std}(\bar Y) = \frac{\sigma}{\sqrt{s}}$$

这就是随机算法常常“抽很多次再平均”的原因：在保持无偏性的同时降低波动。

总结

注意在进行以上分析的时候，我们使用的都是广义的估计目标，所以得出的结论并不局限于一种估计目标，而是普适的。也就是说，对于这类通过多次随机抽样并取平均来构造输出的 Monte Carlo estimator，如果希望最终估计量是无偏的，就要保证 多次抽样并取平均后的估计量仍然满足

$$\mathbb{E}[\bar Y]=\theta$$

也就是长期平均中心等于真实目标。这里的“无偏”不是说单次运行结果一定等于 $\theta$；单次结果仍然可能偏离 $\theta$，只是多次抽样并取平均以后这种偏离的波动会变得越来越小。抽样次数越多，$\bar Y$ 的抽样分布越集中，算法输出通常越靠近目标值 $\theta$

• 单次抽样 分析告诉我们这个随机估计量是否“方向正确”和单次波动有多大

• 多次抽样取平均 后的分析告诉我们实际算法输出有多稳定。

随机算法的性能通常看最终输出的估计量，而最终输出往往是多次随机估计的平均。单次分析是基础，多次平均分析是最终可靠性的来源。

期望（expectation）和方差（variance）是针对随机变量在其“真实概率分布”下所有可能结果来定义的，不是针对某一次抽样出来的 n 个对象直接定义的