数据压缩（九）——线性预测器推导及最小二乘总结

本周任务：推导线性预测器及总结最小二乘——梯度下降、牛顿法、高斯牛顿法

线性预测器

最小二乘法

最小二乘法是数学优化技术。

原理：通过最小误差的平方和寻求数据的最佳函数匹配。

用于求得未知的数据且求得数据与实际数据之间误差的平方和最小；曲线拟合。

公式

∑ j = 1 n X i j β j = y i , ( i = 1 , 2 , 3... m m > n ) \sum_{j=1}^n X_{ij}β_j=y_i,(i=1,2,3...m m>n) ∑j=1nXijβj=yi,(i=1,2,3...mm>n)，量化后

X β = y Xβ=y Xβ=y

X = [ X 11 X 12 . . . X 1 n X 21 X 22 . . . X 2 n ⋮ ⋮ ⋱ ⋮ X m 1 X m 2 . . . X m n ] X=\begin{bmatrix}X_{11} & X_{12} & ... & X_{1n}\\ X_{21}&X_{22}&...&X_{2n}\\ \vdots&\vdots&\ddots&\vdots\\X_{m1}&X_{m2}&...&X_{mn}\end{bmatrix} X=⎣⎢⎢⎢⎡X11X21⋮Xm1X12X22⋮Xm2......⋱...X1nX2n⋮Xmn⎦⎥⎥⎥⎤， β = [ β 1 β 2 ⋮ β n ] β=\begin{bmatrix}β_{1} \\ β_{2}\\ \vdots\\β_{n}\end{bmatrix} β=⎣⎢⎢⎢⎡β1β2⋮βn⎦⎥⎥⎥⎤， y = [ y 1 y 2 ⋮ y m ] y=\begin{bmatrix}y_{1} \\ y_{2}\\ \vdots\\y_{m}\end{bmatrix} y=⎣⎢⎢⎢⎡y1y2⋮ym⎦⎥⎥⎥⎤

引入残差平方函数和 S ( β ) = ∣ ∣ X β − y ∣ ∣ 2 S(β)=||Xβ-y||^2 S(β)=∣∣Xβ−y∣∣2，当 β = β ^ β=\hat{β} β=β^时， S ( β ) S(β) S(β)取最小值

对 S ( β ) S(β) S(β)进行微分求最值，得 X T X β ^ = X T y X^TX\hat{β}=X^Ty XTXβ^=XTy

当矩阵X^TX非奇异，则β有唯一解， β ^ = ( X T X ) − 1 X T y \hat{β}=(X^TX)^{-1}X^Ty β^=(XTX)−1XTy

拟合

对给定数据点集合 ( X i , y i ) ( i = 0 , 1 , 2 , … , m ) {(X_i,y_i)}(i=0,1,2,…,m) (Xi,yi)(i=0,1,2,…,m)，在取定得函数类 φ \varphi φ中，求 p ( X ) ∈ φ p(X)\in\varphi p(X)∈φ，使误差得平方和 E 2 E^2 E2最小， E 2 = ∑ [ p ( X i ) − y i ] 2 E^2=\sum[p(X_i)-y_i]^2 E2=∑[p(Xi)−yi]2。

从几何意义上讲，就是寻求与给定点集 ( X i , y i ) ( i = 0 , 1 , 2 , . . . , m ) {(X_i,y_i)}(i=0,1,2,...,m) (Xi,yi)(i=0,1,2,...,m)的距离平方和为最小的曲线 y = p ( x ) y=p(x) y=p(x)。

函数 p ( x ) p(x) p(x)为拟合函数（最小二乘解）。求拟合函数的方法为曲线拟合的最小二乘法。
写成矩阵形式：
A x = b Ax=b Ax=b
A = [ A 11 A 12 . . . A 1 k A 21 A 22 . . . A 2 k ⋮ ⋮ ⋱ ⋮ A n 1 A n 2 . . . X n k ] A=\begin{bmatrix}A_{11} & A_{12} & ... & A_{1k}\\ A_{21}&A_{22}&...&A_{2k}\\ \vdots&\vdots&\ddots&\vdots\\A_{n1}&A_{n2}&...&X_{nk}\end{bmatrix} A=⎣⎢⎢⎢⎡A11A21⋮An1A12A22⋮An2......⋱...A1kA2k⋮Xnk⎦⎥⎥⎥⎤， x = [ x 1 x 2 ⋮ x k ] x=\begin{bmatrix}x_{1} \\ x_{2}\\ \vdots\\x_{k}\end{bmatrix} x=⎣⎢⎢⎢⎡x1x2⋮xk⎦⎥⎥⎥⎤， b = [ b 1 b 2 ⋮ b n ] b=\begin{bmatrix}b_{1} \\ b_{2}\\ \vdots\\b_{n}\end{bmatrix} b=⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤
分 n > k n>k n>k和 n < k n<k n<k两种情况。

法一：计算 m i n ∣ ∣ A x − b ∣ ∣ min||Ax-b|| min∣∣Ax−b∣∣，即 A T A x = A T b A^TAx=A^Tb ATAx=ATb，得 x x x。低效

法二：对 A A A进行QR分解（ A = Q R A=QR A=QR），

m i n ∣ ∣ A x − b ∣ ∣ = m i n ∣ ∣ Q R x − b ∣ ∣ = m i n ∣ ∣ R x − Q T b ∣ ∣ min||Ax-b||=min||QRx-b||=min||Rx-Q^Tb|| min∣∣Ax−b∣∣=min∣∣QRx−b∣∣=min∣∣Rx−QTb∣∣

matlab实现：

x=[0.5,1,1.5,2,2.5,3];
y=[1.75,2.45,3.81,4.8,7,8.6];
p=polyfit(x,y,2);
x1=0.5:0.5:3;
y1=polyval(p,x1);
plot(x,y,'*r',x1,y1,'-b');

梯度下降

沿梯度下降得方向求解极小值或沿上升方向求最大值。

迭代公式： a k + 1 = a k + ρ k s ⃗ ( k ) a_{k+1}=a_k+\rho_k\vec{s}^{(k)} ak+1=ak+ρks (k)， s ⃗ ( k ) \vec{s}^{(k)} s (k)表示梯度方向， ρ k \rho_k ρk表示梯度方向上得搜索步长。

步长的确定：线性搜索算法，把下一个点的坐标看作 a k + 1 a_{k+1} ak+1的函数，求满足 f ( a k + 1 ) f(a_{k+1}) f(ak+1)的最小值{a_k+1}。

%求f(x)=x^2的最小值
sym x;f=x^2;
step=0.1;
x=2;
k=0;
f_change=x^2;
f_current=x^2;
ezplot(@(x,f)f-x.^2)
axis([-2,2,-0.2,3])
hold on;
while f_change>0.000000001x=x-step*2*x;f_change=f_current-x^2;f_current=x^2;plot(x,f_current,'ro','markrsize',7)drawnow;plause(0.2);k=k+1;
end
hold off

牛顿法

求解方程：
对于需要求解优化函数 f ( x ) f(x) f(x)，利用泰勒公式，在 x 0 x_0 x0处展开到一阶， f ( x ) = f ( x 0 ) + ( x − x 0 ) f ′ ( x 0 ) f(x)=f(x_0)+(x-x_0)f'(x_0) f(x)=f(x0)+(x−x0)f′(x0)，求解方程 f ( x ) = 0 f(x)=0 f(x)=0，即 f ( x 0 ) + ( x − x 0 ) f ′ ( x 0 ) = 0 f(x_0)+(x-x_0)f'(x_0)=0 f(x0)+(x−x0)f′(x0)=0，求解 x = x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x=x_1=x_0-\frac{f(x_0)}{f'(x_0)} x=x1=x0−f′(x0)f(x0)。

利用泰勒公式的一阶展开， f ( x ) = f ( x 0 ) + ( x − x 0 ) f ′ ( x 0 ) f(x)=f(x_0)+(x-x_0)f'(x_0) f(x)=f(x0)+(x−x0)f′(x0)处不完全相等，为近似相等， x 1 x_1 x1不能让 f ( x ) = 0 f(x)=0 f(x)=0，但 f ( x 1 ) f(x_1) f(x1)的值比 f ( x 0 ) f(x_0) f(x0)更接近于 f ( x ) = 0 f(x)=0 f(x)=0。

通过迭代求解，可以推得 x ( n + 1 ) = x ( n ) − f [ x ( n ) ] f ′ [ x ( n ) ] x(n+1)=x(n)-\frac{f[x(n)]}{f'[x(n)]} x(n+1)=x(n)−f′[x(n)]f[x(n)]，通过迭代，该式在 f ( x ∗ ) = 0 f(x*)=0 f(x∗)=0时收敛。过程用图表示：

最优化：
牛顿法通常用于非线性优化问题。假设任务为优化一个目标函数f，求函数的极大极小问题，则可转化为求解函数f的导数 f ′ = 0 f'=0 f′=0的问题。

二维情况：
求解 f ′ ( x ) = 0 f'(x)=0 f′(x)=0，把 f ( x ) f(x) f(x)展开到二阶形式：
f ( x + δ x ) = f ( x ) + f ′ ( x ) Δ x + 1 2 f ′ ′ ( x ) Δ x 2 f(x+\delta x)=f(x)+f'(x)\Delta x+\frac{1}{2}f''(x)\Delta x^2 f(x+δx)=f(x)+f′(x)Δx+21f′′(x)Δx2

当且仅当Δx无限趋近于0时，上式等价为 f ′ ( x ) + f ′ ′ ( x ) Δ x = 0 f'(x)+f''(x)\Delta x=0 f′(x)+f′′(x)Δx=0
解得： Δ x = − f ′ ( x n ) f ′ ′ ( x n ) \Delta x=-\frac{f'(x_n)}{f''(x_n)} Δx=−f′′(xn)f′(xn)

迭代公式： x n + 1 = x n − f ′ ( x n ) f ′ ′ ( x n ) , n = 0 , 1 , . . . x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)},n=0,1,... xn+1=xn−f′′(xn)f′(xn),n=0,1,...

拓展到n维：
x n + 1 = x n − [ H f ( x n ) ] − 1 ∇ f ( x n ) , n ≥ 0 x_{n+1}=x_{n}-[Hf(x_n)]^{-1}\nabla f(x_n),n\geq0 xn+1=xn−[Hf(xn)]−1∇f(xn),n≥0

牛顿法和梯度下降法对比

牛顿法可以利用到曲线本身的信息，比梯度下降法更容易收敛（迭代更少次数），如下图是一个最小化一个目标方程的例子，红色曲线是利用牛顿法迭代求解，绿色曲线是利用梯度下降法求解。

高斯牛顿法

仅用于解决非线性最小二乘问题。

和牛顿法比较：优点是不需要计算二阶导数矩阵H(f)，缺点是高斯牛顿法仅适用于与最小二乘的问题。

假设 γ ( X i ) \gamma(X_i) γ(Xi)是最小二乘的残差，其目标应使残差和 f ( θ ) f(\theta) f(θ)最小。
即 f ( θ ) = 1 2 ∑ i = 1 m [ r ( X i ) ] 2 f(\theta)=\frac{1}{2}\sum_{i=1}^{m}{[r(X_i)]^2} f(θ)=21∑i=1m[r(Xi)]2

利用牛顿法求解 f ( θ ) f(\theta) f(θ)的最小值，计算其梯度向量和Hessian矩阵。
梯度向量： ∇ γ ( θ ) = α f α θ = ∑ i = 1 m γ i α γ i α θ = [ ∇ θ γ ( x 1 ) T ∇ θ γ ( x 2 ) T ⋮ ∇ θ γ ( x m ) T ] T [ γ ( x 1 ) γ ( x 2 ) ⋮ γ ( x m ) ] \nabla_\gamma(\theta)=\frac{\alpha f}{\alpha\theta}=\sum_{i=1}^{m}\gamma_i\frac{\alpha\gamma_i}{\alpha\theta}=\begin{bmatrix}\nabla_\theta\gamma(x_1)^T \\ \nabla_\theta\gamma(x_2)^T \\ \vdots\\\nabla_\theta\gamma(x_m)^T \end{bmatrix}^T\begin{bmatrix}\gamma(x_1) \\ \gamma(x_2) \\ \vdots\\ \gamma(x_m) \end{bmatrix} ∇γ(θ)=αθαf=∑i=1mγiαθαγi=⎣⎢⎢⎢⎡∇θγ(x1)T∇θγ(x2)T⋮∇θγ(xm)T⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎡γ(x1)γ(x2)⋮γ(xm)⎦⎥⎥⎥⎤，
假设 J γ ( θ ) J_\gamma(\theta) Jγ(θ)是 γ \gamma γ对 θ \theta θ求导的Jacobian矩阵，得，
∇ θ f = J θ ( γ ) ( n × m ) T γ = ∑ i = 1 m γ i ∇ θ γ i \nabla_\theta f=J_\theta(\gamma)_{(n×m)}^T\gamma=\sum_{i=1}{m}\gamma_i\nabla_\theta\gamma_i ∇θf=Jθ(γ)(n×m)Tγ=∑i=1mγi∇θγi

得到Hessian矩阵 H = ∑ i = 1 m ( ∇ θ γ i ) T ( ∇ θ γ i ) + γ i ( H θ γ i ) H=\sum_{i=1}^{m}(\nabla_\theta\gamma_i)^T(\nabla_\theta\gamma_i)+\gamma_i(H_\theta\gamma_i) H=∑i=1m(∇θγi)T(∇θγi)+γi(Hθγi)

假设残差 γ i \gamma_i γi很小，可忽略残差项，即
H = ∑ i = 1 m ( ∇ θ γ i ) ( ∇ θ γ i ) T = J θ T ( γ ) J θ ( γ ) H=\sum_{i=1}^{m}(\nabla_\theta\gamma_i)(\nabla_\theta\gamma_i)^T=J_\theta^T(\gamma)J_\theta(\gamma) H=∑i=1m(∇θγi)(∇θγi)T=JθT(γ)Jθ(γ)

将梯度向量和HESSIAN矩阵带入后，得到高斯牛顿法得迭代公式：
θ i = θ i − 1 − ( J γ T J γ ) ( n × n ) − 1 J γ T γ ( m × 1 ) \theta_i=\theta_{i-1}-(J_\gamma^TJ_\gamma)_{(n×n)}^{-1}J_\gamma^T\gamma_{(m×1)} θi=θi−1−(JγTJγ)(n×n)−1JγTγ(m×1)