RNN 与状态空间模型

17.1 节的 MLP 和 CNN 都是“一次看一张”：同样的输入永远给出同样的输出。可机器人关心的很多量不在任何一帧读数里。关节编码器只给角度，速度要靠前后几帧之差；腿式机器人的本体感觉（proprioception，关节角、角速度、IMU）里没有“地面有多滑”这一项，只有当腿踩下去、关节怎样响应时才看得出来；夹爪闭合后物体被手挡住，相机这一帧看不到它，但上一秒看到过。这些问题的共同点是：答案藏在历史里，用 20.1 节的话说，观测只是部分可观测的。

经典的答案是滤波：4.8 节的卡尔曼滤波维护一个状态估计，每来一个测量就更新一次，更新规则来自人写的运动模型和噪声模型。学习的答案是同一个结构，但让数据来写规则：给网络一个隐藏状态 \(\mathbf h_t\)，每一步用当前输入和上一步的状态算出新的状态。这就是循环神经网络（RNN）。它有两个大家族：带非线性和“门”的 RNN（LSTM、GRU），以及把递推保持为线性的状态空间模型（S4、Mamba），后者正是 13.2 节的 \(\mathbf x_{k+1}=A_d\mathbf x_k+B_d\mathbf u_k\)，矩阵由训练得到。

为什么不干脆把最近 \(n\) 帧拼起来喂给 MLP，或者用 17.3 节的 Transformer 看整段历史？也可以，而且很多机器人系统就是这样做的（第 6 节）。循环模型的卖点是：每个控制周期只花固定的计算和固定的内存，历史再长也一样；这在几十到几百赫兹的控制回路、跑在机载计算机上的策略里很重要。本节从“状态 = 记忆”出发，现场训练一个 RNN，看清它训练时的主要困难（梯度随时间消失或爆炸）和各种解法。

RNN 的全部定义就是一个递推式：

\(\mathbf x_t\) 是第 \(t\) 步的输入（比如关节角读数），\(\mathbf h_t\) 是隐藏状态，\(\hat{\mathbf y}_t\) 是输出，\(\phi\) 通常是 \(\tanh\)。所有时刻共用同一组参数 \(\theta=(W_h,W_x,\mathbf b,W_y,\mathbf b_y)\)：这和 CNN 在空间上共用卷积核是同一个想法，只是换到了时间上。注意记号：13.2 节里 \(\mathbf x\) 是状态、\(\mathbf u\) 是输入；机器学习的习惯是 \(\mathbf x\) 为输入、\(\mathbf h\) 为状态，本节沿用后者。把 \(\phi\) 换成恒等映射，RNN 就是 13.2 节的离散线性系统，\(W_h\) 对应 \(A_d\)，\(W_x\) 对应 \(B_d\)。

线性的情形最好理解记忆。在第 5 步给一个脉冲 \(\mathbf x_5\)，之后输入为零，则 \(\mathbf h_{5+k}=W_h^k W_x\mathbf x_5\)：脉冲的痕迹每一步乘一次 \(W_h\)。把 \(W_h\) 分解到特征向量上，每个分量按特征值 \(\hl{lam}{\lambda}\) 的幂变化：\(|\lambda|<1\) 时按 \(|\lambda|^k\) 衰减，半衰期是 \(\ln\tfrac12/\ln|\lambda|\) 步；\(|\lambda|>1\) 时指数增长；复特征值 \(\lambda=\rho e^{\pm i\omega}\) 让状态每步转过 \(\omega\) 弧度，在平面上画出螺旋。这正是 13.2 节“离散时间稳定 \(\iff\) 特征值在单位圆内”。下图是一个 2 维隐藏状态，\(W_h=\rho R(\omega)\)（旋转乘缩放），它的特征值就是 \(\rho e^{\pm i\omega}\)；拖动特征值，看同一个脉冲被记住多久。

任务：一个关节的角度传感器有噪声（标准差 \(\sigma\)，单位 rad；比如廉价电位器或从图像估计的角度），以 50 Hz 采样（\(\Delta t=0.02\) s）。我们要从读数序列 \(z_1,z_2,\dots\) 实时估计关节角速度 \(\dot q_t\)。关节的运动是两个随机正弦之和（频率 0.3 到 1.5 Hz，幅值 0.15 到 0.6 rad），每段 80 步（1.6 s）。三种估计器同台比较：

• 差分：\((z_t-z_{t-1})/\Delta t\)。几乎没有偏差（只滞后半步），但噪声被放大 \(\sqrt2/\Delta t\) 倍，\(\sigma=0.01\) rad 时就是约 0.7 rad/s 的抖动。
• 卡尔曼滤波（4.8 节）：常速度模型 \([q,\dot q]\)，加速度当白噪声，过程噪声强度在另一批轨迹上网格搜索调到最好。
• RNN：12 个 tanh 隐藏单元，输入 \(x_t=2z_t\)，输出 \(\hat y_t\approx\dot q_t/3\)（缩放只为让数值在 1 附近）。每步训练抽 8 段新轨迹，损失是前 10 步（0.2 s 的热身）之后的均方误差，用 Adam（学习率 0.01 并逐渐减小，梯度范数裁剪到 1）更新。

怎么求梯度？把递推按时间展开：80 步的 RNN 就是一个 80 层、每层权重相同的深网络，再用 16.3 节的反向传播从最后一步往回算。这就叫沿时间反向传播（BPTT）。因为每一层是同一个 \(W_h\)，它的梯度是各步贡献之和。长序列上常用截断 BPTT：每 \(K\) 步把往回传的梯度切断（状态照样往前传，只是梯度不再往更早传），省内存，也让每次更新更快。下面的 \(K\) 就是它。

BPTT 里，第 \(T\) 步的误差要影响第 \(T-k\) 步的参数，必须穿过 \(k\) 个雅可比的乘积 \(\partial\mathbf h_T/\partial\mathbf h_{T-k}=J_T J_{T-1}\cdots J_{T-k+1}\)，每个 \(J_t=\operatorname{diag}(\phi'(\mathbf a_t))\,W_h\)。16.3 节在深的前馈网络里已经见过这种乘积；时间上的版本更糟，因为每一层都是同一个 \(W_h\)，乘积就像矩阵幂：典型的增益略小于 1，几十步后就小到没有意义（消失）；略大于 1，就大到数值溢出（爆炸）。下图对随机初始化的网络直接把这个乘积算出来（隐藏维数 16，输入是随机噪声），画它的最大奇异值随 \(k\) 的变化。三种单元：普通 tanh RNN、GRU、LSTM（门的定义见下一小节和第 5 节）。“增益” \(g\) 是循环权重的尺度（谱半径约为 \(g\)），“门偏置” \(b\) 让门偏向“保留”：LSTM 的遗忘门偏置取 \(+b\)，GRU 的更新门偏置取 \(-b\)。

爆炸好办：训练时把梯度范数裁剪到一个上限（FIG 2 就裁剪到 1），方向不变只缩短步长（Pascanu 等，2013）。消失难办得多：梯度小到不出错也不报警，网络只是学不会长程依赖。门控单元给出的解法是：给记忆开一条加法的通路。

LSTM（Hochreiter 与 Schmidhuber，1997；遗忘门由 Gers 等在 2000 年加入）另设一个“细胞状态” \(\mathbf c_t\)，更新方式是 \(\mathbf c_t=\mathbf f_t\odot\mathbf c_{t-1}+\mathbf i_t\odot\mathbf g_t\)：遗忘门 \(\mathbf f_t\) 决定保留多少旧内容，输入门 \(\mathbf i_t\) 决定写入多少新内容，门都是 0 到 1 之间的 sigmoid，由当前输入和上一步输出算出。沿 \(\mathbf c\) 这条直通路径，\(\partial\mathbf c_t/\partial\mathbf c_{t-1}=\operatorname{diag}(\mathbf f_t)\)：没有矩阵相乘、没有 tanh 的压缩，只要 \(\mathbf f_t\approx1\)，梯度就几乎原样传回去。GRU（Cho 等，2014）更省：只有一个状态，用更新门 \(\hl{gate}{\mathbf z_t}\) 在“保留旧状态”和“换成候选值”之间插值，\(\mathbf h_t=(1-\mathbf z_t)\odot\mathbf h_{t-1}+\mathbf z_t\odot\tilde{\mathbf h}_t\)，\(\mathbf z_t\approx0\) 就是原样保留。

下面用一个需要长记忆的任务对比普通 RNN 和 GRU。机器人在 \(t=0\) 看到一个提示（比如“放进左边还是右边的箱子”，编码为 \(+1\) 或 \(-1\)），接下来 \(D\) 步同一个通道上只有干扰噪声（标准差 1，和提示一样大），最后一步来一个“执行”信号，网络要输出当初的提示。两个网络都是 8 个隐藏单元，用同样的种子、同样的数据、同样的优化器同时训练。GRU 的更新门偏置初始化为 \(-3\)（\(z\approx0.05\)，每步保留 95%），这和给 LSTM 的遗忘门加正偏置是同一个技巧。

还可以把 GRU 的门偏置改成 0 再训练 \(D=30\) 或 \(D=60\)：默认种子下 2000 步内 GRU 也学不会。门本身只提供了“可以保留”的通路，初始化得让网络一开始就倾向于保留，梯度才过得去：结构和初始化缺一不可。

另一条路是反过来：把递推里的非线性去掉，非线性放到层与层之间。一层就是一个连续时间线性系统 \(\dot{\mathbf h}=A\mathbf h+B x,\ y=C\mathbf h\)（每个输入通道一份），按步长 \(\Delta\) 用 13.2 节的零阶保持离散化：\(\bar A=e^{\Delta A}\)，\(\bar B=A^{-1}(e^{\Delta A}-I)B\)，于是 \(\mathbf h_t=\bar A\mathbf h_{t-1}+\bar Bx_t\)。线性时不变带来一个大礼：把递推展开，

整层就是输入与一个固定卷积核 \(\bar{\mathbf K}\) 的卷积。训练时整段序列已知，先算出卷积核，再用 FFT 一次卷完，\(O(L\log L)\) 且在时间上完全并行；推理时又可以退回递推，每步只做 \(O(N)\) 的运算、只存 \(N\) 维状态。S4（Gu、Goel 与 Ré，2022）用 HiPPO 理论设计 \(A\) 的结构与初值，让状态能记很长的历史；之后的 S4D 把 \(A\) 取成对角复矩阵，例如 S4D-Lin 的 \(A_n=-\tfrac12+i\pi n\)，每个模态就是一个独立的复数标量递推。步长 \(\Delta\) 是可学习的“时间尺度”：卷积核的长度与 \(1/\Delta\) 成正比。

递推本身也能并行。\(h_t=a_th_{t-1}+b_t\) 是仿射映射的复合，而复合满足结合律：\((a_2,b_2)\hl{scan}{\bullet}(a_1,b_1)=(a_2a_1,\ a_2b_1+b_2)\)。所以所有前缀 \(h_1,\dots,h_L\) 可以用并行前缀扫描（parallel scan）在 \(\lceil\log_2L\rceil\) 轮里算完，每一轮所有位置同时计算。扫描的另一个好处是：它不要求 \(a_t\) 与 \(t\) 无关；卷积则要求。这一点对下一小节至关重要。

拖动 \(\Delta\) 看卷积核：S4D-Lin、8 个模态时，\(\Delta=0.05\) 的核在约 190 步后衰减到峰值的 1% 以下，\(\Delta=0.01\) 时约 960 步：\(\Delta\) 缩小 5 倍，记忆拉长 5 倍。直观地说，\(\Delta\) 决定了“一步”在连续时间里有多长，模态衰减率 \(\operatorname{Re}A_n=-\tfrac12\) 是按连续时间计的。

线性时不变的层有个根本局限：它对每个输入一视同仁，同一个卷积核作用在所有位置上，没办法“这一帧重要就记住，那一帧是垃圾就跳过”。Mamba（Gu 与 Dao，2023）的做法是让 \(\Delta,B,C\) 都成为当前输入的函数，例如 \(\hl{delta}{\Delta_t}=\operatorname{softplus}(\mathbf w^\top\mathbf x_t+\beta)\)。取最简单的标量 \(A=-1\)、\(B=1\)，零阶保持给出

\(\Delta_t\to0\) 时 \(h_t\approx h_{t-1}\)：保留旧状态，忽略这一步的输入；\(\Delta_t\gg1\) 时 \(h_t\approx x_t\)：用这一步的输入覆盖状态。这就是 GRU 的更新门 \(z_t=1-e^{-\Delta_t}\)！区别在于：门只由输入决定、不依赖 \(h_{t-1}\)，状态更新仍然是线性的，所以依旧可以并行扫描（但不能再用卷积）。下面的任务需要这种能力：一串传感器读数 \(v_t\) 里只有带标记 \(m_t=1\) 的才有效（比如只有在足端着地时力传感器的读数才有意义，或者在夹爪闭合那一刻要锁存一个读数），目标是输出最近一次被标记的读数。一个是 4 个模态、所有参数都可学的线性时不变层；一个是只有 1 个模态、5 个参数的选择性层。

记预激活 \(\mathbf a_t=W_h\mathbf h_{t-1}+W_x\mathbf x_t+\mathbf b\)，\(\mathbf h_t=\phi(\mathbf a_t)\)，总损失 \(L=\sum_t\ell_t(\hat{\mathbf y}_t,\mathbf y_t)\)。令 \(\bar{\mathbf g}_t=\partial L/\partial\mathbf h_t\) 为对 \(\mathbf h_t\) 的全导数（既经过本步输出，也经过所有未来步）。链式法则给出从后往前的递推：

递推从 \(\boldsymbol\delta_{T+1}=\mathbf 0\) 开始。这就是 16.3 节反向传播在展开图上的样子：前向时要存下全部 \(\mathbf h_t\)（内存 \(O(TH)\)），反向的计算量是前向的一个小倍数。截断 BPTT 只改一处：若 \(t+1\) 是一个长度为 \(K\) 的块的起点，就把 \(W_h^\top\boldsymbol\delta_{t+1}\) 这一项置零。FIG 2 的训练代码就是这几行，页面里的实现用中心差分做过梯度检查（误差小于 \(10^{-8}\)）。

把递推展开，第 \(T\) 步对第 \(t\) 步状态的敏感度是

tanh 的 \(\gamma=1\)。所以 \(W_h\) 的最大奇异值小于 1 时梯度必然按指数消失；反过来，最大奇异值大于 \(1/\gamma\) 是爆炸的必要条件（Pascanu 等，2013）。实际网络还要乘上 \(\phi'(\mathbf a_s)\le1\)：输入越大、tanh 越饱和，有效增益越小，所以 FIG 3 里 \(g=1\) 时 RNN 仍在消失，要到 \(g\approx1.8\) 左右才不再衰减，\(g\approx2\) 以上才开始爆炸。爆炸用梯度裁剪 \(\mathbf g\leftarrow\mathbf g\cdot\min(1,c/\|\mathbf g\|)\) 处理；消失要靠结构（门、线性递推）和初始化。

LSTM 沿细胞状态的直通雅可比是 \(\partial\mathbf c_t/\partial\mathbf c_{t-1}=\operatorname{diag}(\mathbf f_t)\)，GRU 的雅可比里有一项 \(\operatorname{diag}(1-\mathbf z_t)\)，再加上经过门和候选值的其余项；门偏向“保留”时，这一项接近单位阵，乘很多次也不衰减。所以常把 LSTM 的 \(\mathbf b_f\) 初始化为 1 或更大（Jozefowicz 等，2015 推荐 1），GRU 的 \(\mathbf b_z\) 取负值（FIG 4 取 \(-3\)）。两点约定上的差别要留意：PyTorch 的 GRU 写成 \(\mathbf h_t=(1-\mathbf z_t)\odot\tilde{\mathbf h}_t+\mathbf z_t\odot\mathbf h_{t-1}\)，它的 \(\mathbf z\) 是这里的 \(1-\mathbf z\)；它还把 \(\mathbf r_t\) 乘在 \(U_n\mathbf h_{t-1}\) 之后而不是之前，是一个小变体。

沿用 4.8 节的记号（那里 \(\mathbf x\) 是状态、\(\mathbf z\) 是测量）：预测 \(\hat{\mathbf x}^-=F\hat{\mathbf x}_{k-1}\)，更新 \(\hat{\mathbf x}_k=\hat{\mathbf x}^-+K_k(\mathbf z_k-H\hat{\mathbf x}^-)\)。代入得

增益 \(K_k\) 很快收敛到稳态值，此后这就是一个 \(\phi\) 为恒等映射、参数固定的线性 RNN，状态维数 2（FIG 2 的卡尔曼基线就是它，输出取速度分量）。差别只在参数从哪里来：卡尔曼滤波从模型 \((F,H,Q,R)\) 推出来，RNN 从数据里学出来。这也解释了 FIG 2：当模型（常速度、白噪声加速度）与真实运动不符时，从数据学参数的那一方有机会赢。

对角化之后每个模态 \(n\) 独立：

复模态成对出现（\(A_n\) 与 \(\bar A_n\) 共轭），只存一半再取 \(2\operatorname{Re}\) 就得到实数输出。只要 \(\operatorname{Re}A_n<0\)，任何步长下都有 \(|\bar A_n|=e^{\Delta\operatorname{Re}A_n}<1\)：零阶保持把左半平面映进单位圆（13.2 节），所以这样参数化的层在训练中不会变得不稳定，这正是 FIG 1 那种“越过单位圆就失控”的问题的结构性解法。实际实现里常把 \(\operatorname{Re}A_n\) 写成 \(-e^{p_n}\) 之类的形式，保证它恒为负（FIG 6 的时不变层用 \(\bar a_n=e^{-\operatorname{softplus}(p_n)}\in(0,1)\)，同一个道理）。

把第 \(t\) 步写成元素 \(e_t=(a_t,b_t)\)，表示映射 \(h\mapsto a_th+b_t\)。两个映射先后作用的结果仍是这种形式：

有结合律就可以随意加括号，于是“所有前缀积”能分轮并行计算。FIG 5 用的 Hillis–Steele 扫描在第 \(\ell\) 轮让每个 \(t\ge2^{\ell-1}\) 同时做 \(e_t\leftarrow e_t\bullet e_{t-2^{\ell-1}}\)，\(\lceil\log_2L\rceil\) 轮后完成，总工作量 \(O(L\log L)\)；Blelloch（1990）的工作量最优版本只需 \(O(L)\) 次合并、约 \(2\log_2L\) 轮。注意整个推导只用到“\(a_t,b_t\) 已知”，从没用到它们与 \(t\) 无关，所以选择性的 \(a_t=e^{-\Delta_t}\) 照样适用；Mamba 正是用一个针对 GPU 内存层级优化过的扫描来训练的。

Mamba 让 \(\hl{delta}{\Delta_t}=\operatorname{softplus}(\mathbf w^\top\mathbf x_t+\beta)\)，\(B_t\)、\(C_t\) 也是 \(\mathbf x_t\) 的线性函数，\(A\) 本身不依赖输入；FIG 6 的选择性层是它最小的版本：\(\Delta_t=\operatorname{softplus}(w_vv_t+w_mm_t+\beta)\)，\(A=-1\)，\(B=1\)，\(y_t=c\,h_t+d\)。默认种子训练 1500 步后 \(w_m\approx9.8\)、\(\beta\approx-5.9\)：有标记的步 \(\Delta_t\approx3.9\)（写入 98%），没有标记的步 \(\Delta_t\approx0.003\)（每步保留 99.7%），读数里能看到这几个数随训练变化。

为什么时不变的层做不到？读数 \(v_t\) 在 \([-1,1]\) 上独立均匀分布（方差 \(\tfrac13\)），标记以概率 \(p=0.12\) 独立出现。任何时不变线性层的输出都形如 \(\hat y_t=\sum_k\alpha_kv_{t-k}+\sum_k\beta_km_{t-k}+c_0\)（\(c_0\) 为常数）。目标 \(y_t=v_{t-\tau}\)，\(\tau\) 是距上一个标记的步数，\(P(\tau=k)=p(1-p)^k\)，所以 \(\operatorname{Cov}(y_t,v_{t-k})=\tfrac13p(1-p)^k\)，而 \(y_t\) 与所有 \(m_s\) 不相关（\(v\) 均值为零且与标记独立）。输入在时间上互不相关，最优系数就是 \(\alpha_k=p(1-p)^k\)、\(\beta_k=0\)，能解释的方差比例为

（这是长序列的极限，此时 \(\operatorname{Var}(y_t)\to\tfrac13\)。）最好的时不变层也只能比“输出 0”好约 6%；FIG 6 里训练出的时不变层停在约 8%（64 步的短序列开头有一段目标恒为 0，比长序列极限稍容易一点），也就在这附近；它学到的衰减率是什么都无所谓，问题在结构。选择性层把乘积 \(m_t\cdot v_t\) 这种输入之间的乘法交互放进了递推，线性时不变系统天生没有这种能力。

FIG 2 是一个缩小的典型用法：机器人测得到 \(q\)，测不到 \(\dot q\)，就从历史里推。放大之后，同样的思路用来推断那些根本没有传感器的量：地面摩擦和软硬、负载质量、电机实际输出的力矩、物体在手里的姿态。做法分两派：一派用固定长度的历史窗口（几十到几百步的本体感觉），交给 MLP、一维卷积（TCN）或 Transformer；另一派用循环状态，每个控制周期更新一次。前者训练简单、没有“状态带到哪里”的问题；后者每步代价固定、能记住窗口之外的东西。选哪个往往取决于控制频率、机载算力和任务需要的记忆长度。