YSYSIMONCODE. ART. SCIENCE. VISION.

Woodbury 公式、RLS 与 Kalman Filter

从 rank-one inverse update 推到 Sherman-Morrison 和 Woodbury，并连接递归最小二乘与 Kalman Filter。

May 16, 2026约 32 分钟MATRIXysysimon
MathMatrix MethodsData AnalysisNumerical Linear AlgebraWoodburyKalman Filter

前面一章讨论随机矩阵乘法时，我们把一个矩阵乘法拆成很多 rank-one terms 的总和：

AB=∑j=1najbj⊤\mathbf{A}\mathbf{B} = \sum_{j=1}^n \mathbf{a}_j\mathbf{b}_j^\top

这种写法的重点是：大矩阵可以拆分为单独处理的小结构。随机算法利用这些小结构做抽样估计，而这一章要讨论的是另一个方向：如果一个大矩阵只发生了 rank-one 或 low-rank 的变化，那么能不能不重新求完整的逆矩阵？

答案就是 Sherman-Morrison formula 和 Woodbury formula。

它们背后的直觉是：

大矩阵的低秩更新⟶小矩阵的求逆或少量线性方程求解\boxed{ \text{大矩阵的低秩更新} \quad \longrightarrow \quad \text{小矩阵的求逆或少量线性方程求解} }

这种思想会自然连接到 Recursive Least Squares，也就是 递归最小二乘；如果再往前走一步，则能看到它和 Kalman filter 的关系。

Rank-one 更新的逆

先从最简单的形式开始。设

u,v∈Rn\mathbf{u},\mathbf{v}\in\mathbb{R}^n

那么

uv⊤\mathbf{u}\mathbf{v}^\top

是一个 n×nn\times n 矩阵，但它的 rank 最多只有 11。因为它把任意向量 x\mathbf{x} 变成

uv⊤x=u(v⊤x)\mathbf{u}\mathbf{v}^\top\mathbf{x} = \mathbf{u}(\mathbf{v}^\top\mathbf{x})

也就是说，输出结果永远落在 u\mathbf{u} 张成的一维方向上。

现在考虑

I−uv⊤\mathbf{I}-\mathbf{u}\mathbf{v}^\top

它看起来是一个完整的 n×nn\times n 矩阵，但实际上只是单位矩阵被一个 rank-one 矩阵修正了一下。

Sherman-Morrison formula 给出它的逆可以直接由下列表达式计算：

(I−uv⊤)−1=I+uv⊤1−v⊤u\boxed{ (\mathbf{I}-\mathbf{u}\mathbf{v}^\top)^{-1} = \mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}} }

前提是

1−v⊤u≠01-\mathbf{v}^\top\mathbf{u}\neq 0

这里最容易让人奇怪的地方是分母，为什么会出现

1−v⊤u1-\mathbf{v}^\top\mathbf{u}

而不是别的东西？

关键在于 v⊤u\mathbf{v}^\top\mathbf{u} 是一个标量：

v⊤u∈R\mathbf{v}^\top\mathbf{u}\in\mathbb{R}

所以当两个 rank-one 矩阵连续相乘时，中间那一段可以被压缩成一个普通数：

uv⊤uv⊤=u(v⊤u)v⊤=(v⊤u)uv⊤\begin{align*} \mathbf{u}\mathbf{v}^\top\mathbf{u}\mathbf{v}^\top &= \mathbf{u}(\mathbf{v}^\top\mathbf{u})\mathbf{v}^\top \\ &= (\mathbf{v}^\top\mathbf{u})\mathbf{u}\mathbf{v}^\top \end{align*}

这就是整个公式能化简的核心，我们直接验证一下，将 I−uv⊤\mathbf{I}-\mathbf{u}\mathbf{v}^\top 和它的逆 I+uv⊤1−v⊤u\mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}} 相乘，看看能否得到单位矩阵：

(I−uv⊤)(I+uv⊤1−v⊤u)=I+uv⊤1−v⊤u−uv⊤−uv⊤uv⊤1−v⊤u=I+uv⊤1−v⊤u−uv⊤−(v⊤u)uv⊤1−v⊤u\begin{align*} &(\mathbf{I}-\mathbf{u}\mathbf{v}^\top) \left( \mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}} \right) \\ &= \mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}} - \mathbf{u}\mathbf{v}^\top - \frac{ \mathbf{u}\mathbf{v}^\top\mathbf{u}\mathbf{v}^\top }{ 1-\mathbf{v}^\top\mathbf{u} } \\ &= \mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}} - \mathbf{u}\mathbf{v}^\top - \frac{ (\mathbf{v}^\top\mathbf{u})\mathbf{u}\mathbf{v}^\top }{ 1-\mathbf{v}^\top\mathbf{u} } \end{align*}

为了把最后三项看得更清楚，先记

a=v⊤ua=\mathbf{v}^\top\mathbf{u}

于是上面的式子可以写成

I+uv⊤1−a−uv⊤−auv⊤1−a=I+(11−a−1−a1−a)uv⊤\begin{align*} \mathbf{I} &+ \frac{\mathbf{u}\mathbf{v}^\top}{1-a} - \mathbf{u}\mathbf{v}^\top - \frac{a\mathbf{u}\mathbf{v}^\top}{1-a} \\ &= \mathbf{I} + \left( \frac{1}{1-a} -1 - \frac{a}{1-a} \right) \mathbf{u}\mathbf{v}^\top \end{align*}

现在只需要检查括号里的标量系数：

11−a−1−a1−a=1−a1−a−1=1−1=0\begin{align*} \frac{1}{1-a} -1 - \frac{a}{1-a} &= \frac{1-a}{1-a} -1 \\ &= 1-1 \\ &= 0 \end{align*}

也就是说，所有带有 uv⊤\mathbf{u}\mathbf{v}^\top 的项正好相互抵消，所以最后只剩下

I\mathbf{I}

这说明上面的式子确实是 (I−uv⊤)−1(\mathbf{I}-\mathbf{u}\mathbf{v}^\top)^{-1}。

从这个验证过程可以看出，分母来自 rank-one 矩阵反复相乘以后，唯一留下来的那个标量反馈：

v⊤u\mathbf{v}^\top\mathbf{u}

从 Sherman-Morrison 到 Woodbury

如果不是 rank-one 更新，而是 rank-kk 更新，就把两个向量换成两个矩阵：

U,V∈Rn×k\mathbf{U},\mathbf{V}\in\mathbb{R}^{n\times k}

于是

UV⊤∈Rn×n\mathbf{U}\mathbf{V}^\top\in\mathbb{R}^{n\times n}

但它的 rank 最多只有 kk。当 k≪nk\ll n 时，它仍然是一个低秩扰动。

对应的公式是

(I−UV⊤)−1=I+U(Ik−V⊤U)−1V⊤\boxed{ (\mathbf{I}-\mathbf{U}\mathbf{V}^\top)^{-1} = \mathbf{I} + \mathbf{U} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{U})^{-1} \mathbf{V}^\top }

这里默认小矩阵 Ik−V⊤U\mathbf{I}_k-\mathbf{V}^\top\mathbf{U} 可逆，否则右侧的逆不存在。

这可以看成 Sherman-Morrison 的 rank-kk 版本。原来 rank-one 公式里的标量

1−v⊤u1-\mathbf{v}^\top\mathbf{u}

现在变成了一个 k×kk\times k 小矩阵：

Ik−V⊤U\mathbf{I}_k-\mathbf{V}^\top\mathbf{U}

这就是低秩更新真正节省的地方。左边要求的是一个 n×nn\times n 矩阵的逆，而右边真正新增的求逆只是

(Ik−V⊤U)−1(\mathbf{I}_k-\mathbf{V}^\top\mathbf{U})^{-1}

如果 k≪nk\ll n，这比重新处理一个大矩阵便宜得多。

更一般地，如果被更新的不是单位矩阵，而是一个可逆矩阵 A\mathbf{A}，则有

(A−UV⊤)−1=A−1+A−1U(Ik−V⊤A−1U)−1V⊤A−1\boxed{ (\mathbf{A}-\mathbf{U}\mathbf{V}^\top)^{-1} = \mathbf{A}^{-1} + \mathbf{A}^{-1}\mathbf{U} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U})^{-1} \mathbf{V}^\top\mathbf{A}^{-1} }

这里同样默认 Ik−V⊤A−1U\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U} 可逆。

这就是 Woodbury formula 的一种常见形式。

这里要特别注意：Woodbury formula 不是说我们完全不需要知道 A−1\mathbf{A}^{-1}。它更准确的意思是：如果 A\mathbf{A} 的分解已经可用，或者至少我们已经有办法高效求解

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

那么当 A\mathbf{A} 被低秩扰动以后，就不必从头处理新的大矩阵。

这里的重点在于 矩阵分解可以复用。

例如如果已经做过 LU decomposition

A=LU\mathbf{A}=\mathbf{L}\mathbf{U}

那么求解

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

会变成两步：

Ly=b\mathbf{L}\mathbf{y}=\mathbf{b} Ux=y\mathbf{U}\mathbf{x}=\mathbf{y}

第一步是 forward substitution，第二步是 back substitution。真正昂贵的是得到 L\mathbf{L} 和 U\mathbf{U} 的分解过程，但是分解只需要进行一次；一旦分解已经有了，后面如果只是右端项从 b\mathbf{b} 换成另一个向量 c\mathbf{c}，就不需要重新分解 A\mathbf{A}，只需要再做一次代入求解：

Ly=c\mathbf{L}\mathbf{y}=\mathbf{c} Ux=y\mathbf{U}\mathbf{x}=\mathbf{y}

如果 A\mathbf{A} 是 symmetric positive definite matrix，也可以用 Cholesky decomposition

A=LL⊤\mathbf{A}=\mathbf{L}\mathbf{L}^\top

道理一样：分解一次，之后面对不同右端项，只做三角系统的代入求解。

三角系统之所以容易解，是因为它的未知量之间有天然顺序。比如上三角系统

Ux=y\mathbf{U}\mathbf{x}=\mathbf{y}

可以写成

[u11u12⋯u1n0u22⋯u2n⋮⋮⋱⋮00⋯unn][x1x2⋮xn]=[y1y2⋮yn]\begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}

展开后最后一行只含有最后一个未知量：

unnxn=ynu_{nn}x_n=y_n

所以可以先得到

xn=ynunnx_n=\frac{y_n}{u_{nn}}

倒数第二行只含有 xn−1x_{n-1} 和已经求出的 xnx_n，于是可以继续求出 xn−1x_{n-1}。一直往上代回去，就能依次得到

xn−2,xn−3,…,x1x_{n-2},x_{n-3},\dots,x_1

这就是 back substitution。

下三角系统

Ly=b\mathbf{L}\mathbf{y}=\mathbf{b}

可以写成

[ℓ110⋯0ℓ21ℓ22⋯0⋮⋮⋱⋮ℓn1ℓn2⋯ℓnn][y1y2⋮yn]=[b1b2⋮bn]\begin{bmatrix} \ell_{11} & 0 & \cdots & 0 \\ \ell_{21} & \ell_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ \ell_{n1} & \ell_{n2} & \cdots & \ell_{nn} \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}

则刚好相反，第一行只含有第一个未知量，所以可以从 y1y_1 开始一路向下求到 yny_n，这叫 forward substitution。

因此三角系统不需要重新做 Gaussian elimination。对于一个 n×nn\times n 三角矩阵，代入求解的代价大约是

O(n2)O(n^2)

而完整分解一个一般稠密矩阵通常需要大约

O(n3)O(n^3)

所以当同一个矩阵 A\mathbf{A} 要面对很多不同右端项时，先分解一次，再反复做三角回代，会比每次从头求解便宜很多。

所以在数值计算里，表达式

A−1b\mathbf{A}^{-1}\mathbf{b}

通常不应该理解成“先显式形成 A−1\mathbf{A}^{-1}，再乘 b\mathbf{b}”，而应该理解成：

用已经分解好的 A 去解 Ax=b\boxed{ \text{用已经分解好的 } \mathbf{A} \text{ 去解 } \mathbf{A}\mathbf{x}=\mathbf{b} }

有了这个结论以后，我们再来看 Woodbury。原公式是

(A−UV⊤)−1=A−1+A−1U(Ik−V⊤A−1U)−1V⊤A−1(\mathbf{A}-\mathbf{U}\mathbf{V}^\top)^{-1} = \mathbf{A}^{-1} + \mathbf{A}^{-1}\mathbf{U} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U})^{-1} \mathbf{V}^\top\mathbf{A}^{-1}

如果我们真正关心的是解线性方程

(A−UV⊤)x=b(\mathbf{A}-\mathbf{U}\mathbf{V}^\top)\mathbf{x} = \mathbf{b}

也就是计算

x=(A−UV⊤)−1b\mathbf{x} = (\mathbf{A}-\mathbf{U}\mathbf{V}^\top)^{-1}\mathbf{b}

那么只要把 Woodbury 公式右边整体乘上 b\mathbf{b}：

x=(A−UV⊤)−1b=A−1b+A−1U(Ik−V⊤A−1U)−1V⊤A−1b\begin{align*} \mathbf{x} &= (\mathbf{A}-\mathbf{U}\mathbf{V}^\top)^{-1}\mathbf{b} \\ &= \mathbf{A}^{-1}\mathbf{b} + \mathbf{A}^{-1}\mathbf{U} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U})^{-1} \mathbf{V}^\top\mathbf{A}^{-1}\mathbf{b} \end{align*}

这一步只是公式代入，但它暴露出实际需要计算的对象：

A−1b,A−1U,(Ik−V⊤A−1U)−1\mathbf{A}^{-1}\mathbf{b}, \qquad \mathbf{A}^{-1}\mathbf{U}, \qquad (\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U})^{-1}

它们看起来写着 A−1\mathbf{A}^{-1}，但根据我们前面介绍矩阵分解的结论，实际计算时不应该显式形成 A−1\mathbf{A}^{-1}。按照前面矩阵分解的观点，A−1b\mathbf{A}^{-1}\mathbf{b} 的意思就是解

Az=b\mathbf{A}\mathbf{z}=\mathbf{b}

并把结果记作

z=A−1b\mathbf{z}=\mathbf{A}^{-1}\mathbf{b}

同理，A−1U\mathbf{A}^{-1}\mathbf{U} 的意思就是解

AW=U\mathbf{A}\mathbf{W}=\mathbf{U}

并把结果记作

W=A−1U\mathbf{W}=\mathbf{A}^{-1}\mathbf{U}

这里 U\mathbf{U} 有 kk 列，所以 AW=U\mathbf{A}\mathbf{W}=\mathbf{U} 可以理解成同时解 kk 个右端项。由于 A\mathbf{A} 没变，这些求解都可以 复用 同一个 矩阵分解，只是更换右端项并做回代。

把

z=A−1b,W=A−1U\mathbf{z}=\mathbf{A}^{-1}\mathbf{b}, \qquad \mathbf{W}=\mathbf{A}^{-1}\mathbf{U}

代回上面的式子，得到

x=z+W(Ik−V⊤W)−1V⊤z\mathbf{x} = \mathbf{z} + \mathbf{W} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{W})^{-1} \mathbf{V}^\top\mathbf{z}

现在再把小矩阵和小向量单独记出来：

S=Ik−V⊤W,y=V⊤z\mathbf{S} = \mathbf{I}_k-\mathbf{V}^\top\mathbf{W}, \qquad \mathbf{y} = \mathbf{V}^\top\mathbf{z}

那么上面的公式可以写成

x=z+WS−1y\mathbf{x} = \mathbf{z} + \mathbf{W}\mathbf{S}^{-1}\mathbf{y}

这里的 S−1y\mathbf{S}^{-1}\mathbf{y} 仍然不应该理解成先显式求出 S−1\mathbf{S}^{-1}。它和前面的 A−1b\mathbf{A}^{-1}\mathbf{b} 是同一种思想：要求的是逆矩阵作用在一个右端项上的结果。

所以令

t=S−1y\mathbf{t} = \mathbf{S}^{-1}\mathbf{y}

等价于解一个 k×kk\times k 小系统

St=y\mathbf{S}\mathbf{t}=\mathbf{y}

这个小系统同样可以通过矩阵分解来解。只是这里的矩阵大小是 k×kk\times k，而不是 n×nn\times n；当 k≪nk\ll n 时，这个分解和回代都比重新处理大矩阵便宜得多。

最后得到

x=z+Wt\boxed{ \mathbf{x} = \mathbf{z} + \mathbf{W}\mathbf{t} }

这套流程里，大矩阵 A\mathbf{A} 的部分只是反复更换右端项，然后用已有分解回代。真正新出现的求解问题只有一个 k×kk\times k 的小系统。

所以我们并不

显式更新整个逆矩阵\text{显式更新整个逆矩阵}

而是

复用旧矩阵的分解，快速得到低秩更新后的新解\boxed{ \text{复用旧矩阵的分解，快速得到低秩更新后的新解} }

Rank-one 情况：Woodbury 退化为 Sherman-Morrison

上面讨论的是 rank-kk 更新：

A−UV⊤\mathbf{A}-\mathbf{U}\mathbf{V}^\top

如果 k=1k=1，那么 U\mathbf{U} 和 V\mathbf{V} 就各自只有一列。把它们写成向量

U=u,V=v\mathbf{U}=\mathbf{u}, \qquad \mathbf{V}=\mathbf{v}

于是 low-rank update 就变成 rank-one update：

A−uv⊤\mathbf{A}-\mathbf{u}\mathbf{v}^\top

这时 Woodbury formula 会退化成 Sherman-Morrison formula。换句话说，Sherman-Morrison 是 Woodbury 在 k=1k=1 时的特殊情况。

从这个角度也能解释为什么有时不直接写逆矩阵，而是直接讨论线性方程，rank-one 情况下，新的方程是

(A−uv⊤)x=b(\mathbf{A}-\mathbf{u}\mathbf{v}^\top)\mathbf{x} = \mathbf{b}

这和 Woodbury 是同一件事

假设旧问题已经通过 A\mathbf{A} 的分解解过：

Aw=b\mathbf{A}\mathbf{w}=\mathbf{b}

得到

w=A−1b\mathbf{w}=\mathbf{A}^{-1}\mathbf{b}

这里依旧不会显式形成 A−1\mathbf{A}^{-1}，而是用已有分解回代得到 w\mathbf{w}。再额外更换右端项，解一个辅助问题：

Az=u\mathbf{A}\mathbf{z}=\mathbf{u}

也就是

z=A−1u\mathbf{z}=\mathbf{A}^{-1}\mathbf{u}

这里同样不需要显式形成 A−1\mathbf{A}^{-1}，只是复用 A\mathbf{A} 的分解，把右端项换成 u\mathbf{u}。

现在把新方程改写。因为

u=Az\mathbf{u}=\mathbf{A}\mathbf{z}

并且

b=Aw\mathbf{b}=\mathbf{A}\mathbf{w}

所以

(A−uv⊤)x=b(A−Azv⊤)x=AwA(I−zv⊤)x=Aw\begin{align*} (\mathbf{A}-\mathbf{u}\mathbf{v}^\top)\mathbf{x} &= \mathbf{b} \\ (\mathbf{A}-\mathbf{A}\mathbf{z}\mathbf{v}^\top)\mathbf{x} &= \mathbf{A}\mathbf{w} \\ \mathbf{A}(\mathbf{I}-\mathbf{z}\mathbf{v}^\top)\mathbf{x} &= \mathbf{A}\mathbf{w} \end{align*}

如果 A\mathbf{A} 可逆，两边消去 A\mathbf{A}，得到

(I−zv⊤)x=w(\mathbf{I}-\mathbf{z}\mathbf{v}^\top)\mathbf{x} = \mathbf{w}

这就回到了最开始的 Sherman-Morrison 形式。也就是

(I−zv⊤)−1=I+zv⊤1−v⊤z(\mathbf{I}-\mathbf{z}\mathbf{v}^\top)^{-1} = \mathbf{I} + \frac{\mathbf{z}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{z}}

因此

x=(I−zv⊤)−1w\mathbf{x} = (\mathbf{I}-\mathbf{z}\mathbf{v}^\top)^{-1}\mathbf{w}

代入 rank-one 公式：

x=(I+zv⊤1−v⊤z)w=w+zv⊤w1−v⊤z=w+zv⊤w1−v⊤z\begin{align*} \mathbf{x} &= \left( \mathbf{I} + \frac{\mathbf{z}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{z}} \right)\mathbf{w} \\ &= \mathbf{w} + \frac{\mathbf{z}\mathbf{v}^\top\mathbf{w}}{1-\mathbf{v}^\top\mathbf{z}} \\ &= \mathbf{w} + \mathbf{z} \frac{\mathbf{v}^\top\mathbf{w}}{1-\mathbf{v}^\top\mathbf{z}} \end{align*}

所以得到

x=w+zv⊤w1−v⊤z\boxed{ \mathbf{x} = \mathbf{w} + \mathbf{z} \frac{\mathbf{v}^\top\mathbf{w}}{1-\mathbf{v}^\top\mathbf{z}} }

这说明两种说法完全一致：

更新逆矩阵⟺更新线性方程的解\boxed{ \text{更新逆矩阵} \quad \Longleftrightarrow \quad \text{更新线性方程的解} }

但在实际计算中，后者通常更重要。因为我们很多时候并不真的需要整个

(A−uv⊤)−1(\mathbf{A}-\mathbf{u}\mathbf{v}^\top)^{-1}

而只需要它作用在 b\mathbf{b} 上的结果。

递归最小二乘

现在把这个思想应用到 least squares。

假设旧的 least squares 问题已经求解出来了：

Ax≈b\mathbf{A}\mathbf{x}\approx\mathbf{b}

它的 normal equation 是

A⊤Ax^=A⊤b\mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}} = \mathbf{A}^\top\mathbf{b}

也就是说，我们已经有了旧估计 x^old\hat{\mathbf{x}}_{\text{old}}，并且旧的 normal matrix

A⊤A\mathbf{A}^\top\mathbf{A}

的逆，或者它的分解，也已经可以复用。

这里默认 A\mathbf{A} 满列秩，因此 normal matrix A⊤A\mathbf{A}^\top\mathbf{A} 可逆；实际 RLS 里也常通过初始协方差或正则化项保证这个逆稳定存在。

现在问题变成：如果此时新来了一条观测

v⊤x≈bnew\mathbf{v}^\top\mathbf{x} \approx b_{\text{new}}

我们当然可以把它加入原数据，再从头求一次 least squares。但如果数据是一条一条来的，这样每次都重新构造 normal equation、重新分解矩阵，会很浪费。

RLS 想做的是：不从头重算，而是在旧解的基础上，用这条新观测做一次递归更新。

把新观测放进原数据，相当于在原来的数据矩阵下面增加一行：

Anew=[Av⊤],bnew all=[bbnew]\mathbf{A}_{\text{new}} = \begin{bmatrix} \mathbf{A} \\ \mathbf{v}^\top \end{bmatrix}, \qquad \mathbf{b}_{\text{new all}} = \begin{bmatrix} \mathbf{b} \\ b_{\text{new}} \end{bmatrix}

新的 normal matrix 是

Anew⊤Anew=[A⊤v][Av⊤]=A⊤A+vv⊤\begin{align*} \mathbf{A}_{\text{new}}^\top\mathbf{A}_{\text{new}} &= \begin{bmatrix} \mathbf{A}^\top & \mathbf{v} \end{bmatrix} \begin{bmatrix} \mathbf{A} \\ \mathbf{v}^\top \end{bmatrix} \\ &= \mathbf{A}^\top\mathbf{A} + \mathbf{v}\mathbf{v}^\top \end{align*}

右端项也变成

Anew⊤bnew all=[A⊤v][bbnew]=A⊤b+vbnew\begin{align*} \mathbf{A}_{\text{new}}^\top\mathbf{b}_{\text{new all}} &= \begin{bmatrix} \mathbf{A}^\top & \mathbf{v} \end{bmatrix} \begin{bmatrix} \mathbf{b} \\ b_{\text{new}} \end{bmatrix} \\ &= \mathbf{A}^\top\mathbf{b} + \mathbf{v}b_{\text{new}} \end{align*}

所以新问题满足

(A⊤A+vv⊤)x^new=A⊤b+vbnew(\mathbf{A}^\top\mathbf{A} + \mathbf{v}\mathbf{v}^\top) \hat{\mathbf{x}}_{\text{new}} = \mathbf{A}^\top\mathbf{b} + \mathbf{v}b_{\text{new}}

这里的关键是

A⊤A⟶A⊤A+vv⊤\mathbf{A}^\top\mathbf{A} \quad \longrightarrow \quad \mathbf{A}^\top\mathbf{A} + \mathbf{v}\mathbf{v}^\top

只是一个 rank-one update。

如果令

P=(A⊤A)−1\mathbf{P} = (\mathbf{A}^\top\mathbf{A})^{-1}

那么新增一条观测后，新的逆 normal matrix 是

Pnew=(A⊤A+vv⊤)−1\mathbf{P}_{\text{new}} = (\mathbf{A}^\top\mathbf{A}+\mathbf{v}\mathbf{v}^\top)^{-1}

对

M=A⊤A\mathbf{M}=\mathbf{A}^\top\mathbf{A}

使用 Sherman-Morrison 的加号版本：

(M+vv⊤)−1=M−1−M−1vv⊤M−11+v⊤M−1v(\mathbf{M}+\mathbf{v}\mathbf{v}^\top)^{-1} = \mathbf{M}^{-1} - \frac{ \mathbf{M}^{-1}\mathbf{v}\mathbf{v}^\top\mathbf{M}^{-1} }{ 1+\mathbf{v}^\top\mathbf{M}^{-1}\mathbf{v} }

代入 P=M−1\mathbf{P}=\mathbf{M}^{-1}，得到

Pnew=P−Pvv⊤P1+v⊤Pv\boxed{ \mathbf{P}_{\text{new}} = \mathbf{P} - \frac{ \mathbf{P}\mathbf{v}\mathbf{v}^\top\mathbf{P} }{ 1+\mathbf{v}^\top\mathbf{P}\mathbf{v} } }

这说明新逆矩阵可以从旧的 P\mathbf{P} 直接更新出来，而不是重新求

(Anew⊤Anew)−1(\mathbf{A}_{\text{new}}^\top\mathbf{A}_{\text{new}})^{-1}

和前面讨论线性方程时一样，实际计算里也可以把这里的“乘以逆矩阵”理解成：复用旧 normal matrix 的分解，或者更新已有分解，然后通过解线性系统得到需要的结果，而不是显式形成完整逆矩阵。

为了把这个公式和估计量更新联系起来，先定义

K=Pv1+v⊤Pv\boxed{ \mathbf{K} = \frac{\mathbf{P}\mathbf{v}}{1+\mathbf{v}^\top\mathbf{P}\mathbf{v}} }

这个 K\mathbf{K} 叫 gain。它来自 Sherman-Morrison 公式里的低秩修正项。

因为

K=Pv1+v⊤Pv\mathbf{K} = \frac{\mathbf{P}\mathbf{v}}{1+\mathbf{v}^\top\mathbf{P}\mathbf{v}}

所以 Pnew\mathbf{P}_{\text{new}} 的更新式也可以写成

Pnew=P−Kv⊤P\mathbf{P}_{\text{new}} = \mathbf{P} - \mathbf{K}\mathbf{v}^\top\mathbf{P}

旧估计满足

x^old=PA⊤b\hat{\mathbf{x}}_{\text{old}} = \mathbf{P}\mathbf{A}^\top\mathbf{b}

原因是旧 normal equation 为

A⊤Ax^old=A⊤b\mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}}_{\text{old}} = \mathbf{A}^\top\mathbf{b}

两边左乘

P=(A⊤A)−1\mathbf{P}=(\mathbf{A}^\top\mathbf{A})^{-1}

就得到上面的式子。

而新估计是

x^new=Pnew(A⊤b+vbnew)\hat{\mathbf{x}}_{\text{new}} = \mathbf{P}_{\text{new}} (\mathbf{A}^\top\mathbf{b}+\mathbf{v}b_{\text{new}})

原因也一样：新 normal equation 为

(A⊤A+vv⊤)x^new=A⊤b+vbnew(\mathbf{A}^\top\mathbf{A}+\mathbf{v}\mathbf{v}^\top) \hat{\mathbf{x}}_{\text{new}} = \mathbf{A}^\top\mathbf{b}+\mathbf{v}b_{\text{new}}

两边左乘新的逆 normal matrix

Pnew=(A⊤A+vv⊤)−1\mathbf{P}_{\text{new}} = (\mathbf{A}^\top\mathbf{A}+\mathbf{v}\mathbf{v}^\top)^{-1}

就得到新估计的表达式。

现在把 Pnew=P−Kv⊤P\mathbf{P}_{\text{new}}=\mathbf{P}-\mathbf{K}\mathbf{v}^\top\mathbf{P} 代入：

x^new=(P−Kv⊤P)(A⊤b+vbnew)=PA⊤b+Pvbnew−Kv⊤PA⊤b−Kv⊤Pvbnew\begin{align*} \hat{\mathbf{x}}_{\text{new}} &= (\mathbf{P}-\mathbf{K}\mathbf{v}^\top\mathbf{P}) (\mathbf{A}^\top\mathbf{b}+\mathbf{v}b_{\text{new}}) \\ &= \mathbf{P}\mathbf{A}^\top\mathbf{b} + \mathbf{P}\mathbf{v}b_{\text{new}} - \mathbf{K}\mathbf{v}^\top\mathbf{P}\mathbf{A}^\top\mathbf{b} - \mathbf{K}\mathbf{v}^\top\mathbf{P}\mathbf{v}b_{\text{new}} \end{align*}

现在先把能认出来的旧估计替换掉。因为

PA⊤b=x^old\mathbf{P}\mathbf{A}^\top\mathbf{b} = \hat{\mathbf{x}}_{\text{old}}

所以

v⊤PA⊤b=v⊤x^old\mathbf{v}^\top\mathbf{P}\mathbf{A}^\top\mathbf{b} = \mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}}

代回展开式：

x^new=x^old−Kv⊤x^old+(Pv−Kv⊤Pv)bnew\begin{align*} \hat{\mathbf{x}}_{\text{new}} &= \hat{\mathbf{x}}_{\text{old}} - \mathbf{K}\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}} + \left( \mathbf{P}\mathbf{v} - \mathbf{K}\mathbf{v}^\top\mathbf{P}\mathbf{v} \right)b_{\text{new}} \end{align*}

现在只剩下括号里的系数需要处理。令

α=v⊤Pv\alpha = \mathbf{v}^\top\mathbf{P}\mathbf{v}

因为

K=Pv1+α\mathbf{K} = \frac{\mathbf{P}\mathbf{v}}{1+\alpha}

所以有

Pv=(1+α)K\mathbf{P}\mathbf{v} = (1+\alpha)\mathbf{K}

因此括号里的系数变成

Pv−Kv⊤Pv=(1+α)K−αK=K\begin{align*} \mathbf{P}\mathbf{v} - \mathbf{K}\mathbf{v}^\top\mathbf{P}\mathbf{v} &= (1+\alpha)\mathbf{K} - \alpha\mathbf{K} \\ &= \mathbf{K} \end{align*}

于是

x^new=x^old−Kv⊤x^old+Kbnew=x^old+K(bnew−v⊤x^old)\begin{align*} \hat{\mathbf{x}}_{\text{new}} &= \hat{\mathbf{x}}_{\text{old}} - \mathbf{K}\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}} + \mathbf{K}b_{\text{new}} \\ &= \hat{\mathbf{x}}_{\text{old}} + \mathbf{K} (b_{\text{new}}-\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}}) \end{align*}

这就得到 RLS 的核心形式：

x^new=x^old+K(bnew−v⊤x^old)\boxed{ \hat{\mathbf{x}}_{\text{new}} = \hat{\mathbf{x}}_{\text{old}} + \mathbf{K} (b_{\text{new}}-\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}}) }

这里

bnew−v⊤x^oldb_{\text{new}}-\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}}

叫 residual，也常叫 innovation。它表示“新观测”和“旧估计对这条观测的预测”之间差了多少。

而 K\mathbf{K} 控制这次 residual 应该怎样修正当前估计。

所以 RLS 的直觉形式是

新估计=旧估计+gain×新观测误差\boxed{ \text{新估计} = \text{旧估计} + \text{gain} \times \text{新观测误差} }

这个形式非常重要。它说明 recursive least squares 并不是每来一条新数据就从头解一次 least squares，而是保存旧的信息，再根据新观测做一次低秩修正。

和 Kalman Filter 的关系

Kalman filter 可以看成 RLS 思想的进一步推广，但它估计的对象不再是固定参数，而是随时间变化的状态。

RLS 的典型目标是一个固定的未知向量：

x\mathbf{x}

每来一条观测，就把对 x\mathbf{x} 的估计更新一下。

Kalman filter 的目标是一个动态状态：

xt\mathbf{x}_t

例如位置、速度、姿态，或者某个系统在时刻 tt 的隐藏状态。它会先根据模型预测当前状态，再用新观测修正这个预测。

只看更新部分，Kalman filter 的核心形式是

xt+=xt−+Kt(yt−Htxt−)\boxed{ \mathbf{x}_t^+ = \mathbf{x}_t^- + \mathbf{K}_t (\mathbf{y}_t-\mathbf{H}_t\mathbf{x}_t^-) }

这里 xt−\mathbf{x}_t^- 是观测到来前的预测，xt+\mathbf{x}_t^+ 是观测修正后的估计。

观测模型通常写成

yt≈Htxt\mathbf{y}_t \approx \mathbf{H}_t\mathbf{x}_t

所以

yt−Htxt−\mathbf{y}_t-\mathbf{H}_t\mathbf{x}_t^-

就是新观测和旧预测之间的差，也就是 Kalman filter 里的 innovation。

对应的 gain 是

Kt=Pt−Ht⊤(HtPt−Ht⊤+R)−1\mathbf{K}_t = \mathbf{P}_t^-\mathbf{H}_t^\top (\mathbf{H}_t\mathbf{P}_t^-\mathbf{H}_t^\top+\mathbf{R})^{-1}

这和 RLS 的 gain 在结构上很像。它们都在做同一件事情：

根据不确定性决定这次新信息应该被相信多少\boxed{ \text{根据不确定性决定这次新信息应该被相信多少} }

如果旧估计很不确定，或者新观测很可靠，那么 gain 会让新观测产生更大的修正；如果旧估计已经很可靠，或者新观测噪声很大，那么修正就会比较小。

更深一层看，Kalman filter 的测量更新也可以理解成一个加权 least squares 问题：

min⁡x∥x−xt−∥(Pt−)−12+∥yt−Htx∥R−12\min_{\mathbf{x}} \|\mathbf{x}-\mathbf{x}_t^-\|_{(\mathbf{P}_t^-)^{-1}}^2 + \|\mathbf{y}_t-\mathbf{H}_t\mathbf{x}\|_{\mathbf{R}^{-1}}^2

第一项表示不要离原来的预测太远，第二项表示要解释新的观测。两项的权重由协方差控制。

对应的 normal matrix 会出现

(Pt−)−1+Ht⊤R−1Ht(\mathbf{P}_t^-)^{-1} + \mathbf{H}_t^\top\mathbf{R}^{-1}\mathbf{H}_t

如果观测维度远小于状态维度，那么

Ht⊤R−1Ht\mathbf{H}_t^\top\mathbf{R}^{-1}\mathbf{H}_t

就是一种低维观测带来的低秩或低维更新。Woodbury formula 可以把状态维度上的大矩阵求逆，转化成观测维度上的小矩阵求逆：

(HtPt−Ht⊤+R)−1(\mathbf{H}_t\mathbf{P}_t^-\mathbf{H}_t^\top+\mathbf{R})^{-1}

这就是 Kalman gain 里那个小矩阵逆出现的原因之一。

因此这条线可以总结为：

Woodbury→RLS→Kalman filter\boxed{ \text{Woodbury} \rightarrow \text{RLS} \rightarrow \text{Kalman filter} }

它们共同的结构是一种反复出现的模式：

旧信息+低秩新信息⟹用小系统更新结果\boxed{ \text{旧信息} + \text{低秩新信息} \quad \Longrightarrow \quad \text{用小系统更新结果} }

总结

这一章的主线是：当矩阵只发生低秩变化时，不要急着重新求完整逆矩阵。

rank-one 情况下，Sherman-Morrison formula 把

(I−uv⊤)−1(\mathbf{I}-\mathbf{u}\mathbf{v}^\top)^{-1}

写成一个单位矩阵加 rank-one 修正：

(I−uv⊤)−1=I+uv⊤1−v⊤u(\mathbf{I}-\mathbf{u}\mathbf{v}^\top)^{-1} = \mathbf{I} + \frac{\mathbf{u}\mathbf{v}^\top}{1-\mathbf{v}^\top\mathbf{u}}

rank-kk 情况下，Woodbury formula 把 n×nn\times n 的大矩阵逆，转化成 k×kk\times k 的小矩阵逆：

(A−UV⊤)−1=A−1+A−1U(Ik−V⊤A−1U)−1V⊤A−1(\mathbf{A}-\mathbf{U}\mathbf{V}^\top)^{-1} = \mathbf{A}^{-1} + \mathbf{A}^{-1}\mathbf{U} (\mathbf{I}_k-\mathbf{V}^\top\mathbf{A}^{-1}\mathbf{U})^{-1} \mathbf{V}^\top\mathbf{A}^{-1}

但在数值计算中，更应该把它理解成复用 A\mathbf{A} 的矩阵分解：当右端项改变时，只需要做 forward substitution 和 back substitution，而不是显式构造 A−1\mathbf{A}^{-1}。

在 least squares 里，每新增一条观测，normal matrix 发生 rank-one update：

A⊤A⟶A⊤A+vv⊤\mathbf{A}^\top\mathbf{A} \quad \longrightarrow \quad \mathbf{A}^\top\mathbf{A} + \mathbf{v}\mathbf{v}^\top

于是 RLS 可以递归更新解：

x^new=x^old+K(bnew−v⊤x^old)\hat{\mathbf{x}}_{\text{new}} = \hat{\mathbf{x}}_{\text{old}} + \mathbf{K} (b_{\text{new}}-\mathbf{v}^\top\hat{\mathbf{x}}_{\text{old}})

Kalman filter 则把这个结构推广到动态状态估计：

xt+=xt−+Kt(yt−Htxt−)\mathbf{x}_t^+ = \mathbf{x}_t^- + \mathbf{K}_t (\mathbf{y}_t-\mathbf{H}_t\mathbf{x}_t^-)

所以从 Woodbury 到 RLS，再到 Kalman filter，核心思想一直没有变：

不要从头重算；利用旧信息，用低秩更新吸收新信息\boxed{ \text{不要从头重算；利用旧信息，用低秩更新吸收新信息} }