EE7402-Statistical-Signal-Processing-Part1-3:线性模型与最佳线性无偏估计

引入-线性模型

在第一章节中，介绍MVUE，说它比较难以寻找，通常有两种方式。前面已经介绍了从CRLB观察出预测函数$g(x)$，本章节介绍另一种：线性模型。

直接寻找MVUE特别困难，有时也没办法从CRLB直接整出$g(x)$。但是许多问题可以用线性模型来表达，一旦确定了线性模型，就一定可以写成$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \mathbf{I}(\boldsymbol{\theta}) [\mathbf{g}(\boldsymbol{x}) - \boldsymbol{\theta}]$，然后拿到MVUE。

线性模型的表达

假设一个信号，它有原始参数为$\theta$，在信道中传输叠加了噪声，采集到的序列$\boldsymbol{x}$可以表达成：

$$\boldsymbol{x=H\theta+w}$$

其中：

• $\pmb x$是$N\times1$的观测向量（即，采集到的序列），$\pmb x=[x[0],x[1],x[2]…x[N-1]]^T$
• $\pmb \theta$是$P\times1$的参数矩阵，$\pmb \theta=[\theta_1,\theta_2,\theta_3…\theta_P]^T$
• $\pmb H$是$N\times P$的观测矩阵，它描述了参数 (\boldsymbol{\theta}) 与观测 (\boldsymbol{x}) 之间的线性关系（比如不同参数对观测的 “贡献” 方式）
• $\pmb w$是$N\times 1$的噪声矩阵

所谓线性模型，就是认为原始信号的参数与采集的信号存在某种线性的关系。即$\boldsymbol{H\theta}$这一部分

线性模型的VMUE

假设噪声向量服从多元正态分布$\boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0}, \sigma^2\mathbf{I})$($\mathbf{I}$ 是单位矩阵，保证各噪声分量独立）

那么根据高斯分布的PDF，似然函数$p$为：

$$p(\boldsymbol{x}) = \frac{1}{(2\pi\sigma^2)^{N/2}} \exp\left[ -\frac{1}{2\sigma^2} (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})^T (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta}) \right]$$

$\ln p$为：

$$\ln p(\boldsymbol{x}) = -\ln(2\pi\sigma^2)^{N/2} - \frac{1}{2\sigma^2} (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})^T (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})$$

对其求导，拿到一阶导数，首先处理一下$\mathbf{H}\boldsymbol{\theta})^T (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})$这个二次型：

$$(\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})^T (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta}) = \boldsymbol{x}^T\boldsymbol{x} - \boldsymbol{x}^T (\mathbf{H}\boldsymbol{\theta}) - (\mathbf{H}\boldsymbol{\theta})^T \boldsymbol{x} + (\mathbf{H}\boldsymbol{\theta})^T (\mathbf{H}\boldsymbol{\theta})$$

由于(\boldsymbol{x}^T (\mathbf{H}\boldsymbol{\theta}))是一个标量（1×1 矩阵），而标量的转置等于自身，即$\boldsymbol{x}^T (\mathbf{H}\boldsymbol{\theta}) = [\boldsymbol{x}^T (\mathbf{H}\boldsymbol{\theta})]^T = \boldsymbol{\theta}^T \mathbf{H}^T \boldsymbol{x}$；

因此：$(\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta})^T (\boldsymbol{x} - \mathbf{H}\boldsymbol{\theta}) =\boldsymbol{x}^T\boldsymbol{x} - 2\boldsymbol{\theta}^T \mathbf{H}^T \boldsymbol{x} + \boldsymbol{\theta}^T (\mathbf{H}^T \mathbf{H}) \boldsymbol{\theta}$

矩阵求导公式：

1. 若(f(\boldsymbol{\theta}) = \boldsymbol{a}^T \boldsymbol{\theta})（(\boldsymbol{a})是与(\boldsymbol{\theta})无关的向量），则(\frac{\partial f}{\partial \boldsymbol{\theta}} = \boldsymbol{a})；
2. 若(f(\boldsymbol{\theta}) = \boldsymbol{\theta}^T \mathbf{A} \boldsymbol{\theta})（(\mathbf{A})是与(\boldsymbol{\theta})无关的对称矩阵），则(\frac{\partial f}{\partial \boldsymbol{\theta}} = 2\mathbf{A} \boldsymbol{\theta})。

$$\frac{\partial \ln p(\boldsymbol{x})}{\partial \boldsymbol{\theta}}=-\frac{1}{2\sigma^2}[0-2\mathbf{H}^T \boldsymbol{x}+2(\mathbf{H}^T \mathbf{H}) \boldsymbol{\theta}]$$

整理一下，尝试写成$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \mathbf{I}(\boldsymbol{\theta}) [\mathbf{g}(\boldsymbol{x}) - \boldsymbol{\theta}]$：

$$\frac{\partial \ln p(\boldsymbol{x})}{\partial \boldsymbol{\theta}} = \frac{1}{\sigma^2} (\mathbf{H}^T\mathbf{H}) \left[ (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\boldsymbol{x} - \boldsymbol{\theta} \right]$$

此时就可以拿到MVUE的估计函数：

$$\hat{\boldsymbol{\theta}} = (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\boldsymbol{x}$$

同时Fisher信息：

$$\mathbf{I}(\boldsymbol{\theta}) = \frac{1}{\sigma^2} (\mathbf{H}^T\mathbf{H})$$

协方差矩阵最小为：

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}} = \sigma^2 (\mathbf{H}^T\mathbf{H})^{-1}$$

这是任一线性模型的通式。所以说只要模型能够线性表达，就可以很轻易地找到VMUE。

例子

直线拟合

考虑在CRLB中的直线拟合的例子：假设有一条直线，它的采样收到了WGN的影响，有N个观测量：$x[n]=A+Bn+w[n]$。对AB为已知常数，根据观测量对这条曲线的参数A和B进行拟合。

在线性模型中，它可以被表达成$\boldsymbol{x=H\theta+w}$，其中：

$$\boldsymbol{x}=[x[0],x[1],..,x[N-1]]^T;\quad \boldsymbol{w}[w[0],w[1],...,w[N-1]]^T$$ $$\theta=[A,B]^T$$ $$H=\begin{bmatrix} 1 & 0\\ 1 & 1\\ 1 & 2\\ ... & ...\\ 1 & N-1 \end{bmatrix}$$ $$H^TH= \begin{bmatrix} N & \frac{1}{2}N(N-1)\\ \frac{1}{2}N(N-1)& \frac{1}{6}N(N-1)(2N-1) \end{bmatrix}$$

把$(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$乘上观测矩阵$\boldsymbol{x}$就可以拿到最后的观测结果。

曲线拟合

假设有曲线$x(t_n)=\theta_1+\theta_2t+\theta_3t^2+…+\theta_pt^{p-1}+w(t_n)$。$n=0,1,…,N-1$

观测矩阵$\boldsymbol{x}=[x(t_0),x(t_1),x(t_2)…x(t_{N-1})]$。

参数矩阵$\boldsymbol{\theta}=[\theta_1,\theta_2,\theta_3…\theta_p]^T$。

观测矩阵$
\mathbf{H} =
\begin{bmatrix}
1 & t_0 & t_0^2 & \cdots & t_0^{p-1} \\
1 & t_1 & t_1^2 & \cdots & t_1^{p-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & t_{N-1} & t_{N-1}^2 & \cdots & t_{N-1}^{p-1}
\end{bmatrix}
\quad (N \times P)
$

把$(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$乘上观测矩阵$\boldsymbol{x}$就可以拿到最后的观测结果。

最佳线性无偏估计(Best Linear Unbiased Estimators）

引入

实际应用中，即便最小方差无偏估计量（MVUE）存在，也常难以找到。比如数据的PDF未知时，基于CRLB和完全充分统计量的方法就无法使用。

这种情况下，可改用次优估计量。但代价是，我们不知道与 MVUE 相比，性能损失了多少。不过，若次优估计量的方差能确定，且满足当前问题的要求，那使用它也是合理的。

最佳线性无偏估计量（BLUE）是一种常见的次优估计量。这里，估计量的形式被限制为对数据的线性组合，然后我们去寻找无偏且方差最小的那个。

BLUE限定估计量的形式为：

$$\hat{\theta}=\sum_{n = 0}^{N - 1}a_nx[n]$$

其中，常数(a_n)需通过施加 “无偏（零偏差）” 和 “最小方差” 这两个条件来确定。这样就可以不使用噪声的PDF。

然而BLUE存在如下限制：

• BLUE 只有在 MVUE 恰好是线性形式时，才是最优的。
• BLUE 可能是次优的，其性能也可能很差。
• 对于某些问题，使用 BLUE 可能完全不合适。

求解BLUE（BLUE估计的通式）

观测到的数据向量是$\mathbf{x} = [x[0], x[1],\cdots, x[N - 1]]^T$，估计的系数向量是$\mathbf{a} = [a_0, a_1,\cdots, a_{N - 1}]^T$

估计量$\hat{\theta}=\sum_{n = 0}^{N - 1}a_nx[n]=\mathbf{a}^T\mathbf{x}$

估计量的方差可以被表示为：

$$\text{var}(\hat{\theta})=E\left\{\left(\mathbf{a}^T\mathbf{x}-E\{\mathbf{a}^T\mathbf{x}\}\right)^2\right\}=E\left\{\left(\mathbf{a}^T(\mathbf{x}-E\{\mathbf{x}\})\right)^2\right\}$$ $$=E\left\{\mathbf{a}^T(\mathbf{x}-E\{\mathbf{x}\})(\mathbf{x}-E\{\mathbf{x}\})^T\mathbf{a}\right\}=\mathbf{a}^T\mathbf{C}_x\mathbf{a}$$

其中，$\mathbf{C}_x = E\left\{(\mathbf{x}-E\{\mathbf{x}\})(\mathbf{x}-E\{\mathbf{x}\})^T\right\}$是数据的协方差矩阵。

估计量的均值是：

$$E\{\hat{\theta}\}=\sum_{n = 0}^{N - 1}a_nE\{x[n]\}=\mathbf{a}^TE\{\mathbf{x}\}$$

为使估计量无偏，数据的均值必须是(\theta)的线性形式，将其表达成：(E\{\mathbf{x}\}=\mathbf{s}\theta)（其中(\mathbf{s} = [s_0, s_1,\cdots, s_{N - 1}]^T)，用于限定其线性）。那么估计量的均值可以改写为：$E\{\hat{\theta}\}=\mathbf{a}^T\mathbf{s}\theta$。

1. 写出拉格朗日函数（代价函数）

要让(E\{\hat{\theta}\}=\theta)，则需满足(\mathbf{a}^T\mathbf{s}=1)​，这被称为约束条件。这里使用拉格朗日乘数法。待优化的拉格朗日函数（或代价函数）为：

$$L(\mathbf{a},\lambda)=\mathbf{a}^T\mathbf{C}_x\mathbf{a}+\lambda(\mathbf{a}^T\mathbf{s}-1)$$

这里的目标函数（要优化的函数）是$\text{var}(\hat{\theta})=\mathbf{a}^T\mathbf{C}_x\mathbf{a}$。约束条件是(\mathbf{a}^T\mathbf{s} = 1)，将其调整为约束条件=0的构型是项$\mathbf{a}^T\mathbf{s}-1$

拉格朗日函数是在有约束的优化问题中，将约束条件与目标函数结合起来的一种函数，用于找到在满足约束条件下目标函数的极值。

拉格朗日函数的构造是:

$目标函数(要优化的函数) + 拉格朗日乘数 × 约束条件(调整后，使约束等式为 0 的形式)$

其中拉格朗日乘数$\lambda$的作用是把 “约束条件” 和 “目标函数” 融合成一个新的函数（拉格朗日函数）,可以想象成 “约束对目标的影响程度” 的调节因子：

• 如果约束很强（比如 “无偏” 是必须严格满足的硬约束），(\lambda)会调整大小，让优化过程必须优先满足约束；
• 求解时，(\lambda)会和变量一起被算出，最终得到既满足约束、又让目标函数最优的结果

2.对拉格朗日函数进行求解

拉格朗日函数对(\mathbf{a})求偏导：(\frac{\partial L}{\partial \mathbf{a}} = 2\mathbf{C}_x\mathbf{a}+\lambda\mathbf{s})

拉格朗日函数对(\lambda)求偏导：(\frac{\partial L}{\partial \lambda}=\mathbf{a}^T\mathbf{s}-1)

令这两个偏导数为零，得到：

$$\mathbf{a}=-\frac{1}{2}\lambda\mathbf{C}_x^{-1}\mathbf{s},\quad\mathbf{a}^T\mathbf{s}=1 = \mathbf{s}^T\mathbf{a}$$

将$\boldsymbol a$代入$a^Ts$：

$$\mathbf{s}^T\left(-\frac{1}{2}\lambda\mathbf{C}_x^{-1}\mathbf{s}\right)=1,\quad \lambda=-\frac{2}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}}$$

最后得到最优估计系数向量：

$$\mathbf{a}_{opt}=\frac{\mathbf{s}^T\mathbf{C}_x^{-1}}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}}$$

这个估计下，方差是：

$$\text{var}(\hat{\theta})=\frac{1}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}}$$

3.计算估计值

计算该估计量时，我们只需要知道 “缩放后的均值”(\mathbf{s})和协方差矩阵(\mathbf{C}_x)。有了这些之后，直接代入$\mathbf{a}_{opt}\boldsymbol{x}$就行。此时PDF不再需要。

数个例子

白噪声下的DC信号

有直流信号在信道传输，均值为0，方差为$\sigma^2$的白噪声$w$（不是高斯白噪声，不知道PDF），采集信号为$x[n]=A+w[n]$。采集样本有N个，估计直流信号的电平值。

1.拿到期望的线性表达

根据题目，已知：$E\{x[n]\} = A$。那么，这一组观测向量可以表达成：$E\{\mathbf{x}\} = \mathbf{s}A$，其中$ \mathbf{s} = [1,1,\dots,1]^T$是长度为N的全1向量。

2. 拿到协方差矩阵

对于随机向量 (\mathbf{x} = [x[0], x[1], \dots, x[N-1]]^T)，其协方差矩阵 (\mathbf{C}_x) 的元素 ([\mathbf{C}_x]_{i,j}) 定义为：

$$[\mathbf{C}_x]_{i,j} = \text{Cov}(x[i], x[j]) = E\left\{ (x[i] - E\{x[i]\})(x[j] - E\{x[j]\}) \right\}$$ $$x[i] - E\{x[i]\} = (A + w[i]) - A = w[i]$$

因此协方差矩阵可以化简为：

$$[\mathbf{C}_x]_{i,j} = E\left\{ w[i] \cdot w[j] \right\}$$

若两个随机变量 (w[i]) 和 (w[j]) 满足 “不相关”，则他们的协方差为0，即$E\{w[i]w[j]\} - E\{w[i]\}E\{w[j]\} = 0$。因此$E\{w[i]w[j]\} =E\{w[i]\}E\{w[j]\} $

而又由于白噪声的期望是0，所以$[\mathbf{C}_x]_{i,j} = E\left\{ w[i] \cdot w[j] \right\}$在$j\neq i$时为0。在$j=i$时为噪声自己的方差$\sigma^2$。因此，他们的协方差矩阵可以写成：(其中$\mathbf{I}$为$N\times N$的单位矩阵，对角线为1，其余为0)

$$\mathbf{C}_x = \sigma^2 \mathbf{I}$$ $$\mathbf{C}_x^{-1} = \frac{1}{\sigma^2} \mathbf{I}$$

3.对信号进行估计

最优估计是$\mathbf{\theta}_{opt}=\frac{\mathbf{s}^T\mathbf{C}_x^{-1}}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}} \mathbf{x}$，对分子分母分别求一下：

分子：

$$\mathbf{s}^T\mathbf{C}_x^{-1} \mathbf{x} = \mathbf{s}\begin{bmatrix} \frac{1}{\sigma^2} \cdot x[0] + 0 \cdot x[1] + \dots + 0 \cdot x[N-1] \\ 0 \cdot x[0] + \frac{1}{\sigma^2} \cdot x[1] + \dots + 0 \cdot x[N-1] \\ \vdots \\ 0 \cdot x[0] + 0 \cdot x[1] + \dots + \frac{1}{\sigma^2} \cdot x[N-1] \end{bmatrix} =\mathbf{s}^T \frac{1}{\sigma^2} \begin{bmatrix} x[0] \\ x[1] \\ \vdots \\ x[N-1] \end{bmatrix}$$ $$= \frac{1}{\sigma^2} \left( 1 \cdot x[0] + 1 \cdot x[1] + \dots + 1 \cdot x[N-1] \right) = \frac{1}{\sigma^2} \sum_{n=0}^{N-1} x[n]$$

分母：

$$\mathbf{s}^T \mathbf{C}_x^{-1} \mathbf{s} = \mathbf{1}^T \left( \frac{1}{\sigma^2} \mathbf{I} \right) \mathbf{1} = \frac{1}{\sigma^2} \sum_{n=0}^{N - 1} 1 = \frac{N}{\sigma^2}$$

故最估计量为：

$$\theta=\frac{\frac{1}{\sigma^2} \sum_{n=0}^{N-1} x[n]}{\frac{N}{\sigma^2}}=\overline x$$

估计的方差是：

$$\text{var}(\hat{\theta})=\frac{1}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}}=\frac{\sigma^2}{N}$$

讨论：非白噪声下的DC估计

延续之前 “白噪声中直流电平估计” 的例子，但现在噪声不再是 “白噪声”，而是均值仍为0的不相关噪声（不同时刻噪声不相关，但各时刻噪声方差可能不同）。我们仍要估计直流电平 A，观测模型为 (x[n] = A + w[n])（(n = 0,1,\dots,N - 1)），其中 (w[n]) 是不相关噪声。

由于是不相关噪声，其方差$\sigma$要用$\sigma_n$表示。因为它还是独立的，因此协方差矩阵还是对角矩阵，如下：

$$\mathbf{C}_x = \begin{bmatrix} \sigma_0^2 & 0 & \cdots & 0 \\ 0 & \sigma_1^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma_{N - 1}^2 \end{bmatrix}$$

逆矩阵是：

$$\mathbf{C}_x^{-1} = \begin{bmatrix} \frac{1}{\sigma_0^2} & 0 & \cdots & 0 \\ 0 & \frac{1}{\sigma_1^2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{1}{\sigma_{N - 1}^2} \end{bmatrix}$$

均值向量仍是 ( \mathbf{s} = [1,1,\dots,1]^T)（因为噪声期望还是0）

对其进行估计$\mathbf{\theta}_{opt}=\frac{\mathbf{s}^T\mathbf{C}_x^{-1}}{\mathbf{s}^T\mathbf{C}_x^{-1}\mathbf{s}} \mathbf{x}$,分别求分子与分母：

$$\mathbf{s}^T\mathbf{C}_x^{-1} \mathbf{x} = \mathbf{s}\begin{bmatrix} \frac{1}{\sigma_0^2} \cdot x[0] + 0 \cdot x[1] + \dots + 0 \cdot x[N-1] \\ 0 \cdot x[0] + \frac{1}{\sigma_1^2} \cdot x[1] + \dots + 0 \cdot x[N-1] \\ \vdots \\ 0 \cdot x[0] + 0 \cdot x[1] + \dots + \frac{1}{\sigma_n^2} \cdot x[N-1] \end{bmatrix} =\mathbf{s}^T \begin{bmatrix} \frac{1}{\sigma_0^2}x[0] \\ \frac{1}{\sigma_1^2}x[1] \\ \vdots \\ \frac{1}{\sigma_n^2}x[N-1] \end{bmatrix}$$ $$=\sum_{n = 0}^{N - 1} \frac{x[n]}{\sigma_n^2}$$ $$\mathbf{s}^T \mathbf{C}_x^{-1} \mathbf{s} = \mathbf{1}^T \mathbf{C}_x^{-1} \mathbf{1} = \sum_{n=0}^{N - 1} \frac{1}{\sigma_n^2}$$

因此估计值：

$$\hat \theta=\frac{\sum_{n = 0}^{N - 1} \frac{x[n]}{\sigma_n^2}}{\sum_{n=0}^{N - 1} \frac{1}{\sigma_n^2}}$$

估计的方差：

$$\text{var}(\hat{A}) = \frac{1}{\sum_{n = 0}^{N - 1} \frac{1}{\sigma_n^2}}$$

使用BLUE估计多个参数

估计参数的向量表达

与前面CRLB一样，多个参数估计需要将被估计参数$\hat \theta$写成向量，假设由M个被估计的参数：$\hat{\boldsymbol{\theta}} = [\hat{\theta}_1, \hat{\theta}_2, \dots, \hat{\theta}_M]^T$

而根据线性估计的假设，估计量要满足线性关系：$\hat{\boldsymbol{\theta}} = \mathbf{A}\mathbf{x}$。此时估计系数就由以前的$1\times N$向量转换为了$M\times N$的系数矩阵$\mathbf{A}$。被估计向量的的第m个元素值是$\hat{\theta}_m = \sum_{n=0}^{N-1} a_{mn} x[n] \quad (m = 1,2,\dots,M)$。

向量估计下的无偏约束

由于$\hat{\boldsymbol{\theta}} = \mathbf{A}\mathbf{x}$，无偏估计要求(E\{\hat{\boldsymbol{\theta}}\} = \boldsymbol{\theta})，进一步改写就是$E\{\hat{\boldsymbol{\theta}}\} = \boldsymbol{\theta}=AE\{\mathbf{x}\}$

为使估计量无偏，数据的均值必须是(\theta)的线性形式：(E\{\mathbf{x}\} = \mathbf{S}\boldsymbol{\theta})（(\mathbf{S})是(N \times M)矩阵，刻画数据均值与参数(\boldsymbol{\theta})的线性关系）。联立上面二式得到：

$$E\{\hat{\boldsymbol{\theta}}\} = \mathbf{A}E\{\mathbf{x}\} = \mathbf{A}\mathbf{S}\boldsymbol{\theta}$$

不难看出，线性约束的条件已经转化为了$\mathbf{A}\mathbf{S}=\mathbf{I}$，其中$\mathbf{I}$是$M\times M$的单位矩阵。

估计量的协方差矩阵

协方差矩阵(\mathbf{C}_{\hat{\boldsymbol{\theta}}})描述估计量各分量的方差与协方差，定义为：

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}} = E\left\{ (\hat{\boldsymbol{\theta}} - E\{\hat{\boldsymbol{\theta}}\})(\hat{\boldsymbol{\theta}} - E\{\hat{\boldsymbol{\theta}}\})^T \right\}$$

将(\hat{\boldsymbol{\theta}} = \mathbf{A}\mathbf{x})和(E\{\hat{\boldsymbol{\theta}}\} = \mathbf{A}E\{\mathbf{x}\})代入，推导：

矩阵乘后转置满足：$(\mathbf{AB})^T = \mathbf{B}^T \mathbf{A}^T$

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}} = E\left\{ (\mathbf{A}\mathbf{x} - \mathbf{A}E\{\mathbf{x}\})(\mathbf{A}\mathbf{x} - \mathbf{A}E\{\mathbf{x}\})^T \right\}$$ $$= E\left\{ \mathbf{A}(\mathbf{x} - E\{\mathbf{x}\})(\mathbf{x} - E\{\mathbf{x}\})^T \mathbf{A}^T \right\}$$

这里面的$ E\left\{ (\mathbf{x} - E\{\mathbf{x}\})(\mathbf{x} - E\{\mathbf{x}\})^T \right\}$恰好就是观测量的协方差矩阵$\mathbf{C}_x$，因此：

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}}=\mathbf{A}\mathbf{C}_x\mathbf{A}^T$$

使用拉格朗日函数来求最小方差

将(\mathbf{A})按列分块为(\mathbf{A} = [\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_M]^T)（(\mathbf{a}_m)是(\mathbf{A})的第m行，长度为N的行向量）。将(\mathbf{S})按列分块为(\mathbf{S} = [\mathbf{s}_1, \mathbf{s}_2, \dots, \mathbf{s}_M])（(\mathbf{s}_m)是(\mathbf{S})的第m列，长度为N的列向量）。

此时无偏约束(\mathbf{A}\mathbf{S} = \mathbf{I})等价于：

$$\mathbf{a}_m^T \mathbf{s}_n = \delta_{mn}$$

(\delta_{mn})是克罗内克函数，(m = n)时为 1，否则为 0

构造出的拉格朗日函数为：

$$(L_k = \mathbf{a}_k^T \mathbf{C}_x \mathbf{a}_k + \sum_{j=1}^M \lambda_{kj} \left( \mathbf{a}_k^T \mathbf{s}_j - \delta_{kj} \right)$$

对$a$求偏导：

$$\frac{\partial L_k}{\partial \mathbf{a}_k} = 2\mathbf{C}_x \mathbf{a}_k + \sum_{j=1}^M \lambda_{kj} \mathbf{s}_j = \mathbf{0}$$

令(\boldsymbol{\lambda}_k = [\lambda_{k1}, \lambda_{k2}, \dots, \lambda_{kM}]^T)，则上式可简写为：

$$2\mathbf{C}_x \mathbf{a}_k + \mathbf{S} \boldsymbol{\lambda}_k = \mathbf{0}$$

解得：

$$\mathbf{a}_k = -\frac{1}{2} \mathbf{C}_x^{-1} \mathbf{S} \boldsymbol{\lambda}_k$$

代入$\mathbf{a}_k^T \mathbf{s}_j = \delta_{kj}$:

$$\left( -\frac{1}{2} \mathbf{C}_x^{-1} \mathbf{S} \boldsymbol{\lambda}_k \right)^T \mathbf{s}_j =-\frac{1}{2} \boldsymbol{\lambda}_k^T \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{s}_j = \delta_{kj}$$

令(\mathbf{e}_k)为第k个标准基向量（只有第k位为 1，其余为 0），则(\delta_{kj} = \mathbf{e}_k^T \mathbf{e}_j)，进一步推导可得：

$$\boldsymbol{\lambda}_k = -2 \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{e}_k$$

带回$\mathbf{a}_k = -\frac{1}{2} \mathbf{C}_x^{-1} \mathbf{S} \boldsymbol{\lambda}_k$，求得最优估计矩阵：

$$\mathbf{a}_k^{opt} = \mathbf{C}_x^{-1} \mathbf{S} \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{e}_k$$

最优矩阵系数$\mathbf{A}=\left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1}$

估计量$\hat{\boldsymbol{\theta}} = \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{x}$

估计量的协方差：$\mathbf{C}_{\hat{\boldsymbol{\theta}}}=\left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1}$

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}} = E\left\{ (\hat{\boldsymbol{\theta}} - E\{\hat{\boldsymbol{\theta}}\})(\hat{\boldsymbol{\theta}} - E\{\hat{\boldsymbol{\theta}}\})^T \right\}=\mathbf{A}\mathbf{C}_x\mathbf{A}^T$$

代入$\mathbf{A}=\left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1}$：

$$= \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{C}_x \mathbf{C}_x^{-1} \mathbf{S} \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} = \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1}$$

高斯-马尔科夫定理

高斯 - 马尔可夫定理表明：在 “线性模型 + 噪声零均值 + 噪声协方差已知” 的条件下，上述形式的使用BLUE估计的$\hat{\boldsymbol{\theta}} = \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{x}$ 是最佳线性无偏估计量.

证明如下：

对于一般线性模型：$\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w}$，噪声均值为0时，其期望为：$E\{\mathbf{x}\} = \mathbf{H}\boldsymbol{\theta}$。对于噪声，已知其协方差矩阵为$E\{\mathbf{w}\mathbf{w}^T\} = \mathbf{C}$

(\mathbf{C}_x)是数据(\mathbf{x})的协方差矩阵。由于(\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w})，且(\boldsymbol{\theta})是确定参数，因此(\mathbf{x})的协方差矩阵等于噪声(\mathbf{w})的协方差矩阵(\mathbf{C})（即(\mathbf{C}_x = \mathbf{C})）

结合之前 “向量参数 BLUE” 的内容（数据均值需满足(E\{\mathbf{x}\} = \mathbf{S}\boldsymbol{\theta})），这里可将关联矩阵(\mathbf{S})识别为设计矩阵(\mathbf{H})，即(\mathbf{S} = \mathbf{H})。

代入(\mathbf{S} = \mathbf{H})和$mathbf{C}_x = \mathbf{C}$那么向量参数BLUE的通式$\hat{\boldsymbol{\theta}} = \left( \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{S} \right)^{-1} \mathbf{S}^T \mathbf{C}_x^{-1} \mathbf{x}$可以写成：

$$\hat{\boldsymbol{\theta}} = \left( \mathbf{H}^T \mathbf{C}^{-1} \mathbf{H} \right)^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{x}$$

估计量的协方差可以写成：

$$\mathbf{C}_{\hat{\boldsymbol{\theta}}} = \left( \mathbf{H}^T \mathbf{C}^{-1} \mathbf{H} \right)^{-1}$$