【统计学习】递归最小二乘算法与奇异值分解

这篇文章讲关于自适应算法的又一个认识：递归最小二乘（Recursive Least Square）

关于自适应算法的认识（LMS）

假设有一组数据（在不同时刻）：$x(1),x(2),...,x(n)$

在每一时刻有不同的目标：$d(1),d(2),...,d(n)$

我们打算用数据对目标进行估计/预测/逼近。也就是找到一种变换，使得数据经过变换之后能够逼近目标：$g(x(n))\to d(n)$

我们对这个变换要求很简单，线性： $g(x(n))={\omega }^{T}x(n)$

我们只需把精力放在如何寻找这个系数 $\omega$上

接下来我们假定一个损失，均方误差：$\underset{\omega }{\min }E|d(n)-{\omega }^{T}x(n){|}^2$

我们知道这样做出来的线性变换，是维纳滤波。

维纳滤波可以在每一个时间点上做：$E|d(n)-{\omega }^T(n)x(n){|}^2$

问题本身有一个时间轴，也就是数据发生变化的时间轴。而在每一个时间点上，求解系数也需要时间，这也对应一个时间轴。这两个时间轴很难统一起来，因为你在求解上一时刻数据的系数过程当中，数据本身已经发生了变化。因此我们才考虑到自适应（Adaptive）。

自适应的一种理解是把两个时间轴耦合在一起。通常的做法是：在每一个时间点上，无需完全求解系数，只需做出一步调整，随着时间的发展，数据所对应的问题实质地不断变化，系数也在相应地不断做调整，大量的实验表明，这样的做法最后能够得到一个有效的结果。

LMS算法的损失函数把前面的期望去掉了： $\underset{\omega }{\min }|d(n)-{\omega }^{T}x(n){|}^2$，这样做的考虑是基于每一时刻只能取一个样本，用这一个样本代表期望。

LMS的这个操作所导致的就是参数的如下迭代：${\hat{\omega }}_{n+1}={\hat{\omega }}_n+\mu x(n)e(n)$，其中$e(n)=d(n)-{\hat{\omega }}_n^{T}x(n)$是误差项。

典型的最小二乘解（LS）

现在我们要重新考量Loss的定义：
$$\underset{\omega }{\min }E|d(n)-{\omega }^{T}x(n){|}^2$$上述损失函数有什么问题？——只考虑了一个时刻。包括LMS也是，调整所用的数据信息都只用了一个时刻。如果我们对一段历史已经有完全把握，当我们试图做n+1时刻的某些事情，该如何把整段历史利用起来？我们可以这么调整损失函数：
$$\underset{\omega }{\min }\sum _{k=1}^{n}E|d(k)-{\omega }^{T}x(k){|}^2$$这么做虽然可行，但是有点违背马尔可夫假设，因为马尔可夫性质告诉我们，历史上的东西重要性没那么大。如果我既想利用历史信息，又不想太过“负重前行”，那么我可以添加遗忘因子：
$$\underset{\omega }{\min }\sum _{k=1}^n\underset{\begin{array}{rl}& \text{Forgetting}\\ & \text{Factor}\end{array}}{\underbrace{\lambda (n,k)}}E|d(k)-{\omega }^{T}x(k){|}^2$$最常见的遗忘因子是指数遗忘：
$$\underset{\omega }{\min }\sum _{k=1}^n{\lambda }^{n-k}E|d(k)-{\omega }^{T}x(k){|}^2$$由于我们已经是在多个样本数据上做某种意义上的期望了，现在完全可以抛开期望符号，就得到了最终的目标函数：
$$\underset{\omega }{\min }g(\omega )=\underset{\omega }{\min }\sum _{k=1}^n{\lambda }^{n-k}|d(k)-{\omega }^{T}x(k){|}^2$$下面我们来做这个优化：
$$\begin{array}{rl}g(\omega )& =\sum _{k=1}^n{\lambda }^{n-k}{d}^2(k)-2\sum _{k=1}^{n}d(k){\omega }^{T}x(k){\lambda }^{n-k}+\sum _{k=1}^n{\lambda }^{n-k}({\omega }^{T}x(k){)}^2\\ & =\sum _{k=1}^n{\lambda }^{n-k}{d}^2(k)-2{\omega }^T(\sum _{k=1}^n{\lambda }^{n-k}x(k)d(k))+{\omega }^T(\sum _{k=1}^n{\lambda }^{n-k}x(k)x^T(k))\omega \end{array}$$我们用 $Z(n)={\sum }_{k=1}^n{\lambda }^{n-k}x(k)d(k)$ 和 $\Phi (n)={\sum }_{k=1}^n{\lambda }^{n-k}x(k)x^T(k)$ 来简化符号：
$$g(\omega )=C-2{\omega }^{T}Z(n)+{\omega }^T\Phi (n)\omega$$这样的优化目标我们很熟悉，是关于这个线性系数$\omega$​的二次型，我们只需对它求梯度：
$${\nabla }_{\omega }g(\omega )=-2Z(n)+2\Phi (n)\omega =0\Rightarrow {\omega }_n={\Phi }^{-1}(n)Z(n)$$这是典型的最小二乘解，反应了线性估计的精神实质：分母是自相关，分子是互相关。

递归最小二乘（RLS）

现在我们想要让最小二乘解也递推自适应起来：
$${\omega }_{n+1}\leftarrow {\omega }_n$$最小二乘求解和计算的关键在于这个逆${\Phi }^{-1}(n)$，那么自然要从它入手寻求递推关系：
$$\begin{array}{rl}\Phi (n)& =\sum _{k=1}^n{\lambda }^{n-k}x(k)x^T(k)=\sum _{k=1}^{n-1}{\lambda }^{n-k}x(k)x^T(k)+x(n)x^T(n)\\ & =\lambda \sum _{k=1}^{n-1}{\lambda }^{n-1-k}x(k)x^T(k)+x(n)x^T(n)\\ & =\lambda \Phi (n-1)+x(n)x^T(n)\end{array}$$下面考察它的逆：
$${\Phi }^{-1}(n)=(\lambda \Phi (n-1)+x(n)x^T(n){)}^{-1}$$

矩阵求逆公式（Matrix Inversion）

往往一些大的突破，都是在一些小技巧上硬核展开的结果。勿以善小而不为

$$(A+BCD{)}^{-1}=?$$我们想对分块矩阵$\left[\begin{array}{cc}A& B\\ D& -C^{-1}\end{array}\right]$​做一个对角化：
$$\left[\begin{array}{cc}I& BC\\ 0& I\end{array}\right]\left[\begin{array}{cc}A& B\\ D& -C^{-1}\end{array}\right]\left[\begin{array}{cc}I& 0\\ CD& I\end{array}\right]=\left[\begin{array}{cc}A+BCD& 0\\ 0& -C^{-1}\end{array}\right]$$$$\left[\begin{array}{cc}I& 0\\ -D{A}^{-1}& I\end{array}\right]\left[\begin{array}{cc}A& B\\ D& -C^{-1}\end{array}\right]\left[\begin{array}{cc}I& -A^{-1}B\\ 0& I\end{array}\right]=\left[\begin{array}{cc}A& 0\\ 0& -(D{A}^{-1}B+C^{-1})\end{array}\right]$$$$\begin{array}{rl}{\left[\begin{array}{cc}A+BCD& 0\\ 0& -C^{-1}\end{array}\right]}^{-1}& =\left[\begin{array}{cc}I& 0\\ -CD& I\end{array}\right]\left[\begin{array}{cc}I& -A^{-1}B\\ 0& I\end{array}\right]\left[\begin{array}{cc}{A}^{-1}& 0\\ 0& -(D{A}^{-1}B+C^{-1}{)}^{-1}\end{array}\right]\left[\begin{array}{cc}I& 0\\ -D{A}^{-1}& I\end{array}\right]\left[\begin{array}{cc}I& -BC\\ 0& I\end{array}\right]\\ & =\left[\begin{array}{cc}I& -A^{-1}B\\ -CD& \ast \end{array}\right]\left[\begin{array}{cc}{A}^{-1}& 0\\ 0& -(D{A}^{-1}B+C^{-1}{)}^{-1}\end{array}\right]\left[\begin{array}{cc}I& -BC\\ -D{A}^{-1}& \ast \end{array}\right]\\ & =\left[\begin{array}{cc}{A}^{-1}& A^{-1}B(D{A}^{-1}B+C^{-1}{)}^{-1}\\ \ast & \ast \end{array}\right]\left[\begin{array}{cc}I& -BC\\ -D{A}^{-1}& \ast \end{array}\right]\end{array}$$$$(A+BCD{)}^{-1}=A^{-1}-A^{-1}B(D{A}^{-1}B+C^{-1}{)}^{-1}D{A}^{-1}$$

---

回到我们的问题，我们要对 ${\Phi }^{-1}(n)=(\lambda \Phi (n-1)+x(n)x^T(n){)}^{-1}$ 求逆，那么用矩阵求逆公式，我们给定 $A=\lambda \Phi (n-1)$、$B=x(n)$、$C=1$、$D=x^T(n)$：
$${\Phi }^{-1}(n)={\lambda }^{-1}{\Phi }^{-1}(n-1)-\frac{{\lambda }^{-1}{\Phi }^{-1}(n-1)x(n)x^T(n){\Phi }^{-1}(n-1){\lambda }^{-1}}{1+{\lambda }^{-1}{x}^T(n){\Phi }^{-1}(n-1)x(n)}$$从而，我们从本质上通过递推的方式绕过了矩阵求逆。

下一步，我们用 $K(n)=\frac{{\lambda }^{-1}{\Phi }^{-1}(n-1)x(n)}{1+{\lambda }^{-1}{x}^T(n){\Phi }^{-1}(n-1)x(n)}$ 简化递推公式：
$${\Phi }^{-1}(n)={\lambda }^{-1}{\Phi }^{-1}(n-1)-{\lambda }^{-1}K(n)x^T(n){\Phi }^{-1}(n-1)$$我们稍微研究一下这个$K(n)$：
$$K(n)+{\lambda }^{-1}K(n)x^T(n){\Phi }^{-1}(n-1)x(n)={\lambda }^{-1}{\Phi }^{-1}(n-1)x(n)$$$$K(n)=({\lambda }^{-1}{\Phi }^{-1}(n-1)-{\lambda }^{-1}K(n)x^T(n){\Phi }^{-1}(n-1))x(n)={\Phi }^{-1}(n)x(n)$$好了，我们再对 $Z(n)={\sum }_{k=1}^n{\lambda }^{n-k}x(k)d(k)$ 进行递推的转换：
$$\begin{array}{rl}Z(n)& =\sum _{k=1}^n{\lambda }^{n-k}x(k)d(k)=\lambda \sum _{k=1}^{n-1}{\lambda }^{n-1-k}x(k)d(k)+x(n)d(n)\\ & =\lambda Z(n-1)+x(n)d(n)\end{array}$$从而我们就可以改写最小二乘解为递推的形式：
$$\begin{array}{rl}{\omega }_n& ={\Phi }^{-1}(n)Z(n)=\\ & =({\lambda }^{-1}{\Phi }^{-1}(n-1)-{\lambda }^{-1}K(n)x^T(n){\Phi }^{-1}(n-1))(\lambda Z(n-1)+x(n)d(n))\\ & ={\Phi }^{-1}(n-1)Z(n-1)-K(n)x^T(n){\Phi }^{-1}(n-1)Z(n-1)+{\Phi }^{-1}(n)x(n)d(n)\\ & ={\omega }_{n-1}+K(n)(d(n)-{\omega }_{n-1}^{T}x(n))\end{array}$$递归最小二乘的这个结果，从本质上来说，是一种预测矫正过程，换句话说，就是卡尔曼滤波。或者说，我们可以对卡尔曼滤波有一个新认识：卡尔曼滤波所做的优化是跟踪整个历史数据，拿到一起，来做最小二乘优化。

所谓预测矫正，通俗地说就是：先用老结果预测，老结果在新数据上的误差用一个增益加权，把它融汇到老结果里面来，对老结果形成矫正。

奇异值分解（Singular Value Decomposition, SVD）

假设我要分解一个长矩阵 $A=\left(\begin{array}{c}{a}_1^T\\ ⋮\\ a_n^T\end{array}\right)\in {\mathbb{R}}^{n\times m,(n>m)}$

那么 $a_i^T$ 全部是 $m$ 维空间的点，我现在想要找一个超平面，使得这$n$个点到这个超平面的距离的平方和最小。这个问题跟最小二乘有微妙的差异，最小二乘是在做点与超平面垂直方向上的距离优化，而不是直接点与超平面的距离优化。所以，这本质上是对最小二乘理论的前进。

点 $a$ 到超平面 $\mathbf{v}$的距离：
$$d(a,\mathbf{v})=\frac{|a^T\mathbf{v}|}{\Vert \mathbf{v}\Vert }$$假定 $\Vert \mathbf{v}\Vert =1$：
$$d^2(a,\mathbf{v})=|a^T\mathbf{v}{|}^2$$优化目标：
$$\sum _{k=1}^n{d}^2(a_k,\mathbf{v})=\sum _{k=1}^n|a_k^T\mathbf{v}{|}^2=\Vert A\mathbf{v}{\Vert }_2^2$$即，我们要做如下优化：
$$\begin{gathered} \underset{\mathbf{v}}{\min }\Vert A\mathbf{v}{\Vert }_2^2 \\ s.t.\Vert \mathbf{v}\Vert =1 \end{gathered}$$这个解的形式：
$${\mathbf{v}}_1={\mathrm{argmin}}_{\Vert \mathbf{v}\Vert =1}\Vert A\mathbf{v}{\Vert }_2^2$$我们得到了第一个矢量，下面我们要找第二个矢量，使得这 $n$个点到这两个矢量张成的子空间的距离最小。
d(a,{v1,v2})=d(a,v1)+d(a,v2) d(a,\{\mathbf v_1,\mathbf v_2\})=d(a,\mathbf v_1)+d(a,\mathbf v_2) d(a,{v1​,v2​})=d(a,v1​)+d(a,v2​)
我们自然希望上式能够成立，这就要求我们得在${\mathbf{v}}_1$的正交补空间里找${\mathbf{v}}_2$​。
$${\mathbf{v}}_2={\mathrm{argmin}}_{\mathbf{v}\perp {\mathbf{v}}_1,\Vert \mathbf{v}\Vert =1}\Vert A\mathbf{v}{\Vert }_2^2$$
类似的，剩余的矢量我们也是这么找下去：
v3=arg⁡min⁡v⊥{v1，v2},∥v∥=1∥Av∥22 \mathbf v_3=\mathop{\arg\min}_{\mathbf v\perp \{\mathbf v_1，\mathbf v_2\},\|\mathbf v\|=1}\|A\mathbf v\|_2^2 v3​=argminv⊥{v1​，v2​},∥v∥=1​∥Av∥22​
以此类推。假设 $A$ 矩阵的秩为 $\text{rank}(A)=r$，那么找到的矢量就有 $r$ 个 $({\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_r)$ 标准正交。

我们可以进一步扩充为标准正交基 $({\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_r,{\mathbf{v}}_{r+1},...,{\mathbf{v}}_m)$

我们令 ${\mathbf{u}}_k=A{\mathbf{v}}_k/\Vert A{\mathbf{v}}_k\Vert$、${\sigma }_k=\Vert A{\mathbf{v}}_k{\Vert }_2,(k=1,2,...,m)$，于是有
$$A({\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_m)=({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_m)\left[\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_m\end{array}\right]$$

$$A=({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_m)\left[\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_m\end{array}\right]\left(\begin{array}{c}{\mathbf{v}}_1^T\\ ⋮\\ {\mathbf{v}}_m^T\end{array}\right)$$

其中 $({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_m)$ 彼此正交，因为有如下对称矩阵分解 $A{A}^T=U\Lambda U^T$成立：
$$\begin{array}{rl}A{A}^T& =({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_m)\left[\begin{array}{ccc}{\sigma }_1^2& & \\ & \ddots & \\ & & {\sigma }_m^2\end{array}\right]\left(\begin{array}{c}{\mathbf{u}}_1^T\\ ⋮\\ {\mathbf{u}}_m^T\end{array}\right)\\ & =({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_m,{\mathbf{u}}_{m+1},...,{\mathbf{u}}_n)\left[\begin{array}{cccccc}{\sigma }_1^2& & & & & \\ & \ddots & & & & \\ & & {\sigma }_m^2& & & \\ & & & 0& & \\ & & & & \ddots & \\ & & & & & 0\end{array}\right]\left(\begin{array}{c}{\mathbf{u}}_1^T\\ ⋮\\ {\mathbf{u}}_n^T\end{array}\right)\end{array}$$

$\text{rank}(A{A}^T)=\text{rank}(A)=\text{rank}(A^T)$

于是我们便证明了奇异值分解：$A=U\Sigma V^T$

其中 $U\in {\mathbb{R}}^{n\times n},U{U}^T=I_n$、$V\in {\mathbb{R}}^{m\times m},V{V}^T=I_m$

---

奇异值分解下的最小二乘（伪逆）

下面我们打算用奇异值分解对最小二乘重新做一个透视，原问题：
$$\underset{\omega }{\min }\Vert \left(\begin{array}{c}d(1)\\ ⋮\\ d(n)\end{array}\right)-\left(\begin{array}{c}{x}^T(1)\\ ⋮\\ x^T(n)\end{array}\right)\omega {\Vert }_2^2$$也可以写为：
$$\underset{\omega }{\min }\Vert \mathbf{d}-\mathbf{x}\omega {\Vert }_2^2=(\mathbf{d}-\mathbf{x}\omega {)}^T(\mathbf{d}-\mathbf{x}\omega )$$加入指数遗忘因子 $\Lambda =\text{diag}({\lambda }^{n-1},{\lambda }^{n-2},...,1)$：
$$\underset{\omega }{\min }(\mathbf{d}-\mathbf{x}\omega {)}^T\Lambda (\mathbf{d}-\mathbf{x}\omega )$$我们知道最小二乘解为 ${\omega }_{LS}=({\mathbf{x}}^T\mathbf{x}{)}^{-1}{\mathbf{x}}^T\mathbf{d}$

对$\mathbf{x}$做奇异值分解： $\mathbf{x}=U\Lambda V^T$，$\Lambda =\left[\begin{array}{ccc}{\sigma }_1^2& & \\ & \ddots & \\ & & {\sigma }_m^2\\ & 0& \end{array}\right]=\left[\begin{array}{c}\tilde{\Lambda }\\ 0\end{array}\right]$​，就得到如下形式：
$$\begin{array}{rl}{\omega }_{LS}& =(V{\Lambda }^T{U}^{T}U\Lambda V^T{)}^{-1}V{\Lambda }^T{U}^T\mathbf{d}\\ & =(V{\Lambda }^T\Lambda V^T{)}^{-1}V{\Lambda }^T{U}^T\mathbf{d}\\ & =(V{\tilde{\Lambda }}^2{V}^T{)}^{-1}V{\Lambda }^T{U}^T\mathbf{d}\\ & =V\left[\begin{array}{ccc}{\sigma }_1^{-2}& & \\ & \ddots & \\ & & {\sigma }_m^{-2}\end{array}\right]V^{T}V\left[\begin{array}{c}\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_m\end{array}0\end{array}\right]U^T\mathbf{d}\\ & =V\left[\begin{array}{c}\begin{array}{ccc}{\sigma }_1^{-1}& & \\ & \ddots & \\ & & {\sigma }_m^{-1}\end{array}0\end{array}\right]U^T\mathbf{d}\\ & ={\mathbf{x}}^{+}\mathbf{d}\end{array}$$其中，${\mathbf{x}}^{+}=V\left[\begin{array}{c}\begin{array}{ccc}{\sigma }_1^{-1}& & \\ & \ddots & \\ & & {\sigma }_m^{-1}\end{array}0\end{array}\right]U^T$ 称为伪逆（Pseudo Inversion Generalized）

之所以称为伪逆，是因为原先按照线性方程组求解： $\mathbf{x}\omega =\mathbf{d}\Rightarrow \omega ={\mathbf{x}}^{-1}\mathbf{d}$

又因为 $n>m$ 是一个超定方程，无法求精确解，所以最小二乘是在做 $\text{Fitting}$ 。

而现在通过奇异值分解的角度进一步对最小二乘进行透视，我们发现实际上是在求伪逆：$\omega ={\mathbf{x}}^{+}\mathbf{d}$

伪逆求法两部曲：

• 先把 $\Lambda$ 转秩
• 能求逆的部分就求，求不了的那些个零就扔在那。

$${\mathbf{x}}^{+}=V{\Lambda }^{+}{U}^T$$伪逆的性质：

• $\mathbf{x}{\mathbf{x}}^{+}\mathbf{x}=\mathbf{x}$
• ${\mathbf{x}}^{+}\mathbf{x}{\mathbf{x}}^{+}={\mathbf{x}}^{+}$

矩阵求逆，用奇异值分解，它的数值稳定性和效率都是最好的。用奇异值分解伪逆的思想求最小二乘，是大部分软件的做法。原因是求逆这个操作本身，很怕矩阵的刚性比（=|最大特征值|/|最小特征值|）很大。从优化角度来说，这样会导致优化曲面很扁，从而梯度下降时Zigzag的现象会比较严重。