RESEARCH / PRACTICE / NOTES

最小二乘法

从超定矩阵方程出发,讨论普通最小二乘的求解,并以回声路径估计说明如何构造观测矩阵、检验残差及判断适用条件。

#凸优化
关键词最小二乘法线性代数
散落的观测点与穿行其中的直线投影
本文目录 · 9 个章节
  1. 矩阵方程问题描述
  2. 基本问题描述
  3. 最小二乘问题
  4. 普通最小二乘求解
  5. 一个具体问题:估计回声路径
  6. 从采样值到矩阵方程
  7. 一个可以复算的小例子
  8. 普通最小二乘与最大似然
  9. 最小二乘与梯度下降

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以估计未知参数,并使模型预测与观测数据之间的误差平方和最小。最小二乘法还可用于曲线拟合。

矩阵方程问题描述

基本问题描述

许多问题可以建模为矩阵方程形式:

Ax=b\mathbf{A} \mathbf{x}=\mathbf{b}

其中根据向量 b∈Rm×1\mathbf{b} \in \mathbb{R}^{m \times 1} 和矩阵 A∈Rm×n\mathbf{A} \in \mathbb{R}^{m \times n} 的已知条件及维度,常见问题可以分为以下几类:

1. 超定矩阵方程
m > n,b和A均已知,其中之一或二者可能存在观测误差、干扰;
2. 盲矩阵方程
仅向量b已知,矩阵A未知;
3. 欠定矩阵方程
m < n,b和A均已知;若希望得到稀疏解,还需附加稀疏性条件。

每一类问题,都有对应的方法:如对于超定矩阵方程,方程个数大于未知数个数,若矩阵列满秩,可以利用最小二乘得到唯一解;若列不满秩,最小二乘解仍可能不唯一。对于欠定矩阵方程,方程个数小于未知数个数,通常需要添加约束,例如稀疏性;但即使要求最稀疏,也不能仅凭这一要求保证解唯一。问题与对应求解方法的示意图如下:

矩阵方程
求解

超定矩阵方程

  • 普通最小二乘
  • 数据最小二乘
  • 总体最小二乘
  • 约束总体最小二乘

盲矩阵方程

  • 子空间方法
  • 非负矩阵分解

欠定稀疏矩阵方程

  • 正交匹配追踪
  • LASSO 与 LARS
  • 增广拉格朗日乘子法
矩阵方程的三类常见求解路径

对于欠定稀疏矩阵求解,核心问题为:

min⁡x∥x∥0 s.t. Ax=b\min _{x}\|\mathbf{x}\|_{0} \quad \text { s.t. } \quad \mathbf{A x}=\mathbf{b}

这也是稀疏表示和压缩感知的核心问题。在观测带有噪声时,可以把严格等式改为残差约束:

min⁡x∥x∥0 s.t. ∥Ax−b∥2≤ε\min _{x}\|\mathbf{x}\|_{0} \quad \text { s.t. } \quad\|\mathbf{A} \mathbf{x}-\mathbf{b}\|_{2} \leq \varepsilon

最小二乘问题

最小二乘根据噪声类型的不同,可以分为:普通最小二乘、数据最小二乘、总体最小二乘。

  • 普通最小二乘(Ordinary least squares, OLS) 此时,观测向量 b\mathbf b 带有误差,而 A\mathbf A 视为已知,对应的问题为:
b=Ax+ε,min⁡x∥Ax−b∥22.\mathbf b=\mathbf A\mathbf x+\boldsymbol\varepsilon, \qquad \min_{\mathbf x}\|\mathbf A\mathbf x-\mathbf b\|_2^2.
  • 数据最小二乘(Data least squares, DLS) 此时,数据矩阵 A\mathbf A 带有误差,而 b\mathbf b 视为已知,对应的问题是:
(A+ΔA)x=b,min⁡x,ΔA∥ΔA∥F2.(\mathbf A+\Delta\mathbf A)\mathbf x=\mathbf b, \qquad \min_{\mathbf x,\Delta\mathbf A}\|\Delta\mathbf A\|_F^2.
  • 总体最小二乘(Total least squares, TLS) 此时,数据矩阵 A\mathbf A 和观测向量 b\mathbf b 都带有误差,对应的问题是:
(A+ΔA)x=b+Δb,min⁡∥[ΔA,Δb]∥F2.(\mathbf A+\Delta\mathbf A)\mathbf x=\mathbf b+\Delta\mathbf b, \qquad \min\|[\Delta\mathbf A,\Delta\mathbf b]\|_F^2.

本文仅分析超定方程情况,且只讨论普通最小二乘问题。

普通最小二乘求解

问题描述:

min⁡x∥Ax−b∥22=min⁡x(Ax−b)T(Ax−b)\min_{\mathbf x}\|\mathbf A\mathbf x-\mathbf b\|_2^2 =\min_{\mathbf x}(\mathbf A\mathbf x-\mathbf b)^{T}(\mathbf A\mathbf x-\mathbf b)

令梯度为零,得到正规方程:

ATAx=ATb\mathbf{A}^{T} \mathbf{A} \mathbf{x}=\mathbf{A}^{T} \mathbf{b}

case1:列满秩,rank⁡(A)=n\operatorname{rank}(\mathbf A)=n。 此时(ATA)\left(\mathbf{A}^{T} \mathbf{A}\right)可逆,对应的解唯一,从而有:

xOLS=(ATA)−1ATb\mathbf{x}_{O L S}=\left(\mathbf{A}^{T} \mathbf{A}\right)^{-1} \mathbf{A}^{T} \mathbf{b}

case2:秩亏缺,rank⁡(A)<n\operatorname{rank}(\mathbf A)<n。 这种情况下,需要借助Moore-Penrose进行广义逆求解,从而有:

xOLS=A†b\mathbf{x}_{\mathrm{OLS}}=\mathbf A^{\dagger}\mathbf b

**秩亏缺情况下,最小二乘解一般不唯一,但由 Moore–Penrose 伪逆选出的解是唯一的。**这就不免有一个问题——它依据什么条件选解?这里直接给出结论:它是所有最小二乘解中欧氏范数最小的一个,即最小范数最小二乘解(minimum-norm least-squares solution)。数值求解时,直接使用 QR、SVD 或相应的最小二乘例程,通常比显式计算 (ATA)−1(\mathbf A^T\mathbf A)^{-1} 更稳妥。

一个具体问题:估计回声路径

在回声消除中,远端播放信号可以作为参考信号,麦克风则记录了扬声器经过房间、设备外壳等路径传来的回声。如果希望估计这条路径,可以先把它近似为一个长度为 LL 的有限冲激响应(FIR)滤波器:

d[n]=∑k=0L−1hkx[n−k]+s[n]+v[n]d[n]=\sum_{k=0}^{L-1}h_kx[n-k]+s[n]+v[n]

其中,x[n]x[n] 为播放参考,d[n]d[n] 为麦克风观测,hkh_k 为待估的路径系数,s[n]s[n] 为近端说话声,v[n]v[n] 为其他噪声。估计时先选取远端单讲、没有明显近端说话声的片段,使 s[n]s[n] 尽量小;还要使参考信号与麦克风录音在时间上对齐,并假设这段时间内回声路径没有显著变化。

远端播放信号经过声学路径形成回声,与近端声音和噪声叠加后成为麦克风观测;最小二乘利用播放参考和麦克风观测估计声学路径
用播放参考和麦克风观测估计回声路径的信号关系

从采样值到矩阵方程

把待估的滤波器系数写成 h=[h0,h1,…,hL−1]T\mathbf h=[h_0,h_1,\ldots,h_{L-1}]^T。每一个采样时刻对应一行播放参考的历史采样值:

X=[x[0]x[−1]⋯x[1−L]x[1]x[0]⋯x[2−L]⋮⋮⋱⋮x[N−1]x[N−2]⋯x[N−L]],d=[d[0]d[1]⋮d[N−1]].\mathbf X=\begin{bmatrix} x[0]&x[-1]&\cdots&x[1-L]\\ x[1]&x[0]&\cdots&x[2-L]\\ \vdots&\vdots&\ddots&\vdots\\ x[N-1]&x[N-2]&\cdots&x[N-L] \end{bmatrix},\qquad \mathbf d=\begin{bmatrix}d[0]\\d[1]\\\vdots\\d[N-1]\end{bmatrix}.

在片段开始之前没有可用采样值时,可以像下面的简化算例一样用零填充;处理连续录音时,则应尽量保留真实的历史采样值。这样,路径估计就成为上一节的普通最小二乘问题:

h^=arg⁡min⁡h  ∥Xh−d∥22.\hat{\mathbf h}=\underset{\mathbf h}{\arg\min}\;\|\mathbf X\mathbf h-\mathbf d\|_2^2.

这里的 X\mathbf X 不是抽象给定的矩阵,它的每一列就是参考信号的一个延迟版本;d\mathbf d 则是同一时刻的麦克风采样。只有在参考信号提供足够变化、X\mathbf X 的列能够被区分时,各个路径系数才有机会被稳定估计。采样点数多于系数个数只是必要的维度条件,并不能保证矩阵列满秩。

一个可以复算的小例子

下面人为设定一段八点的播放参考和麦克风观测,用两个抽头估计回声路径。数据只为展示计算步骤,并非真实录音:

x=[1,0,1,1,0,1,0,1]T,d=[0.62,0.19,0.60,0.83,0.18,0.61,0.19,0.62]T.\begin{aligned} \mathbf x&=[1,0,1,1,0,1,0,1]^T,\\ \mathbf d&=[0.62,0.19,0.60,0.83,0.18,0.61,0.19,0.62]^T. \end{aligned}

前六点用于估计,后两点留作检验。设 x[−1]=0x[-1]=0,则前六行观测矩阵及对应观测为:

Xfit=[100110110110],dfit=[0.620.190.600.830.180.61].\mathbf X_{\mathrm{fit}}= \begin{bmatrix} 1&0\\0&1\\1&0\\1&1\\0&1\\1&0 \end{bmatrix},\qquad \mathbf d_{\mathrm{fit}}= \begin{bmatrix} 0.62\\0.19\\0.60\\0.83\\0.18\\0.61 \end{bmatrix}.

代入正规方程,得到:

[4113]h^=[2.661.20],h^≈[0.6160.195].\begin{bmatrix}4&1\\1&3\end{bmatrix}\hat{\mathbf h} =\begin{bmatrix}2.66\\1.20\end{bmatrix},\qquad \hat{\mathbf h}\approx\begin{bmatrix}0.616\\0.195\end{bmatrix}.

再用第七、八点对应的参考采样预测麦克风观测:预测值分别约为 0.1950.195 和 0.6160.616,实际观测分别为 0.190.19 和 0.620.62。留出片段的残差较小,说明这组人为设定的数据与两抽头线性模型相符;但两个检验点远远不足以证明模型在真实房间或长时间运行时仍然可靠。

实际计算时,不必按上式显式求 (XTX)−1(\mathbf X^T\mathbf X)^{-1}。例如,可用最小二乘数值例程直接处理同一组数据:

import numpy as np

x = np.array([1, 0, 1, 1, 0, 1, 0, 1], dtype=float)
d = np.array([0.62, 0.19, 0.60, 0.83, 0.18, 0.61, 0.19, 0.62])
X = np.column_stack((x, np.r_[0, x[:-1]]))

h, *_ = np.linalg.lstsq(X[:6], d[:6], rcond=None)
prediction = X[6:] @ h
residual = d[6:] - prediction
print(h, prediction, residual)

对于真实录音,应该用时间上独立的较长片段检查预测残差,并分别看远端单讲、双讲和回声路径变化时的表现。双讲时,近端说话声也进入 d\mathbf d,它不属于回声路径所能解释的部分;这时若仍用普通最小二乘直接更新路径,结果可能偏移。参考与麦克风不同步、滤波器长度不足,或参考信号长时间变化很少,也都会使残差增大或估计不稳定。最小二乘给出的是所设模型与所选数据下的平方误差最优解,并不自动保证模型本身适合现场条件。

普通最小二乘与最大似然

给出数学模型(以多项式拟合为例,N次拟合,共M组样本点):

{a0+a1x1+a2x12+…aNx1N−y1=ε1a0+a1x2+a2x22+…aNx2N−y2=ε2⋯a0+a1xM+a2xM2+…aNxMN−yM=εM\left\{\begin{array}{l} a_{0}+a_{1} x_{1}+a_{2} x_{1}^{2}+\ldots a_{N} x_{1}^{N}-y_{1}=\varepsilon_{1} \\ a_{0}+a_{1} x_{2}+a_{2} x_{2}^{2}+\ldots a_{N} x_{2}^{N}-y_{2}=\varepsilon_{2} \\ \cdots \\ a_{0}+a_{1} x_{M}+a_{2} x_{M}^{2}+\ldots a_{N} x_{M}^{N}-y_{M}=\varepsilon_{M} \end{array}\right.

普通最小二乘准则函数:

J1=∑i=1M∣f(xi)−yi∣2=∑i=1M∣εi∣2J_{1}=\sum_{i=1}^{M}\left|f\left(x_{i}\right)-y_{i}\right|^{2}=\sum_{i=1}^{M}\left|\varepsilon_{i}\right|^{2}

最大似然准则:

假设误差独立同分布,且服从均值为零、方差为 σ2\sigma^2 的正态分布,则有似然函数:

L=∏i=1M12πσexp⁡ ⁣(−εi22σ2)L=\prod_{i=1}^{M}\frac{1}{\sqrt{2\pi}\sigma} \exp\!\left(-\frac{\varepsilon_i^2}{2\sigma^2}\right)

当 σ2\sigma^2 固定时,取负对数并去掉与待估参数无关的常数,最大似然准则等价于:

J2=∑i=1Mεi2.J_2=\sum_{i=1}^{M}\varepsilon_i^2.

二者等价。

最小二乘与梯度下降

上文一个大前提就是方程可以转化为线性变换:

Ax=b\mathbf{A} \mathbf{x}=\mathbf{b}

但很多时候模型不能写成对未知参数的线性形式,非线性最小二乘通常没有类似上面的闭式解,这时可以借助梯度下降等迭代方法求解。梯度下降是一种求解方法,最小二乘是一种以残差平方和为目标的建模准则;线性最小二乘也可以用迭代算法求解,而非线性最小二乘的迭代结果则可能只是局部最优。

交流与讨论

评论加载中…

随便读一篇 →

站内搜索

输入关键词开始搜索。

搜索内容不会上传 Esc 关闭

订阅明月入怀

将下面的网址添加到你的 RSS 阅读器。