State Space Models

参考文献：Fundamentals of Recurrent Neural Network (RNN) and Long Short-Term Memory (LSTM) Network

理论起源：微分方程视角

微分方程视角：RNN 的动力学方程可由一阶非齐次 ODE 推导而来。设 $\vec{s}(t)\in {\mathbb{R}}^d$ 为 $d$ 维状态信号向量，其随时间 $t$ 的演化可表示为：

$$\frac{d\vec{s}(t)}{dt}=\vec{f}(t)+\vec{\varphi }$$

• $\vec{x}(t)$ 是 $d$ 维输入信号向量，$s$ 为 $d$ 为状态向量；
• $\vec{f}(t)$ 是 $d$ 维向量函数，通常依赖于 $\vec{s}(t)$ 和 $\vec{x}(t)$，即 $\vec{f}(t)=\vec{h}(\vec{s}(t),\vec{x}(t))$；
• $\vec{\varphi }$ 是常数 $d$ 维向量（偏差项）。

加性模型：这里采用加性模型 (Additive Model)，将 $\vec{f}(t)$ 分解为三个独立项的线性组合：

$$\vec{f}(t)=\vec{a}(t)+\vec{b}(t)+\vec{c}(t)$$

记 ${\tau }_s,{\tau }_r,{\tau }_x$ 分别表示状态、读出和输入的延迟时间常数。

• 状态反馈项 $\vec{a}(t)$：表示状态信号自身的延迟影响（"模拟"分量）

$$\vec{a}(t)=\sum _{k=0}^{K_s-1}{\vec{a}}_k(\vec{s}(t-{\tau }_s(k)))$$

• 读出反馈项 $\vec{b}(t)$：表示经过非线性变换后的读出信号的延迟影响：

$$\vec{b}(t)=\sum _{k=0}^{K_r-1}{\vec{b}}_k(\vec{r}(t-{\tau }_r(k)))$$

• 外部输入项 $\vec{c}(t)$：表示外部输入信号的延迟影响

$$\vec{c}(t)=\sum _{k=0}^{K_x-1}{\vec{c}}_k(\vec{x}(t-{\tau }_x(k)))$$

延迟微分方程 (DDE) 系统：

$$\frac{d\vec{s}(t)}{dt}=\underset{\text{状态反馈}}{\underbrace{\sum _{k=0}^{K_s-1}{\vec{a}}_k(\vec{s}(t-{\tau }_s(k)))}}+\underset{\text{读出反馈}}{\underbrace{\sum _{k=0}^{K_r-1}{\vec{b}}_k(\vec{r}(t-{\tau }_r(k)))}}+\underset{\text{外部输入}}{\underbrace{\sum _{k=0}^{K_x-1}{\vec{c}}_k(\vec{x}(t-{\tau }_x(k)))}}+\vec{\varphi }$$

线性系数 DDE：当上述函数为线性时，可以转化为具有线性矩阵系数的非线性 DDE：

$$\frac{d\vec{s}(t)}{dt}=\sum _{k=0}^{K_s-1}{A}_k\vec{s}(t-{\tau }_s(k))+\sum _{k=0}^{K_r-1}{B}_k\vec{r}(t-{\tau }_r(k))+\sum _{k=0}^{K_x-1}{C}_k\vec{x}(t-{\tau }_x(k))+\vec{\varphi }$$

其中 $A_k,B_k,C_k$ 分别为对应项的系数矩阵。

简化模型 (The Simplified Model)：通过以下特定约束，可将通用模型简化为包含 Continuous Hopfield Network 和 Cellular Neural Network 的特例形式：

• 单项约束：设 $K_s=K_r=K_x=1$（每种反馈仅有一项）。
• 延迟约束：设 ${\tau }_s(0)=0,{\tau }_x(0)=0$（状态和输入无延迟），仅保留读出信号的单一延迟 ${\tau }_r(0)={\tau }_0$。
• 矩阵重命名：$A_0=A,B_0=B,C_0=C$。

$$\frac{d\vec{s}(t)}{dt}=A\vec{s}(t)+B\vec{r}(t-{\tau }_0)+C\vec{x}(t)+\vec{\varphi }$$

时间离散与 RNN 公式

命题(离散化公式)：使用向后欧拉离散上述方程可得到

$$\vec{s}[n]=W_s\vec{s}[n-1]+\underset{\text{循环/记忆项}}{\underbrace{((\Delta T)W_{s}B)\vec{r}[n-1]}}+\underset{\text{输入项}}{\underbrace{((\Delta T)W_{s}C)\vec{x}[n]}}+\underset{\text{偏差项}}{\underbrace{((\Delta T)W_s\vec{\varphi })}}$$

证明：设采样时间步长为 $\Delta T$，时间 $t=n\Delta T$。使用后向欧拉法近似导数：

$$\frac{d\vec{s}(t)}{dt}\approx \frac{\vec{s}(n\Delta T+\Delta T)-\vec{s}(n\Delta T)}{\Delta T}$$

令延迟时间 ${\tau }_0$ 等于单步采样时间 $\Delta T$（即 ${\tau }_0=\Delta T$）：

$$\frac{\vec{s}[n+1]-\vec{s}[n]}{\Delta T}=A\vec{s}[n+1]+B\vec{r}[n]+C\vec{x}[n+1]+\vec{\varphi }$$

这里采用了简化记号 $\vec{s}[n]\equiv \vec{s}(n\Delta T)$。注意方程右侧使用了 $n+1$ 时刻的 $\vec{s}$ 和 $\vec{x}$（后向法的特性），但使用了 $n$ 时刻的 $\vec{r}$（由延迟 ${\tau }_0$ 决定）。接着

$$(I-(\Delta T)A)\vec{s}[n+1]=\vec{s}[n]+((\Delta T)B)\vec{r}[n]+((\Delta T)C)\vec{x}[n+1]+(\Delta T)\vec{\varphi }$$

定义 $W_s=(I-(\Delta T)A{)}^{-1}$，并在方程两边左乘 $W_s$，得到：

$$\vec{s}[n+1]=W_s\vec{s}[n]+((\Delta T)W_{s}B)\vec{r}[n]+((\Delta T)W_{s}C)\vec{x}[n+1]+((\Delta T)W_s\vec{\varphi })$$

将索引向前平移一步 ($n\to n-1$) 可得结论。

标准 RNN 形式：为了进一步简化表示，定义新的权重矩阵和偏差向量：

$$\begin{array}{rl}\vec{s}[n]& =W_s\vec{s}[n-1]+W_r\vec{r}[n-1]+W_x\vec{x}[n]+{\vec{\theta }}_s\\ \vec{r}[n]& =G(\vec{s}[n])\end{array}$$

其中

$$W_r=(\Delta T)W_{s}B,W_x=(\Delta T)W_{s}C,{\vec{\theta }}_s=(\Delta T)W_s\vec{\varphi }$$

• $W_s$: 状态自身的递归权重（State-to-State）。
• $W_r$: 读出信号（即上一时刻的激活值）的反馈权重（Readout-to-State）。
• $W_x$: 当前输入的权重（Input-to-State）。

稳定性分析

稳定性条件：上述系统稳定的条件是矩阵 $\hat{W}=W_s+W_r$ 的所有特征值必须位于复平面的单位圆内。

标准 RNN 定义：在实际应用中，常作进一步简化以获得最简形式：

• 单位时间步长：设 $\Delta T=1$。
• 快速状态衰减：假设矩阵 $A$ 为对角阵且对角元为很大的负数（$a_{ii}\ll 0$），这意味着状态的衰减非常快。
• 忽略状态记忆：由此导致 $W_s\approx -A^{-1}$ 为对角阵且元素为很小的正数。在这种情况下，状态信号 $\vec{s}[n-1]$ 对当前轨迹的显式影响（第一项）可以忽略不计（尽管通过 $\vec{r}[n-1]$ 的隐式影响依然存在）。
• 忽略 $W_s$ 项（即设 $W_s\approx 0$）

我们得到最常见的标准 RNN 定义，其中 $G(s)$ 通常是 $\tanh$ 函数

$$\begin{array}{rl}\vec{s}[n]& =W_r\vec{r}[n-1]+W_x\vec{x}[n]+{\vec{\theta }}_s\\ \vec{r}[n]& =G(\vec{s}[n])\end{array}$$

在此简化下，系统的稳定性完全取决于 $W_r$ 的特征值 ${\mu }_i$。在“小信号”区域（$\Vert \vec{s}[n]\Vert \ll 1$），稳定性的充要条件是 $0<{\mu }_i<1$。

Diagonal SSM 系统

连续 Diagonal SSM 系统：一般考虑如下 Diagonal SSM LTI 系统

$$\begin{array}{rl}\dot{x}(t)& =Ax(t)+Bu(t),x(0)=x_0\\ y(t)& =Cx(t)\end{array}$$

其中 $u(t)\in \mathbb{R}$，$x(t)\in {\mathbb{C}}^N$，$y(t)\in \mathbb{R}$ （一般使用取实部操作），$A\in {\mathbb{C}}^{N\times N}$，$B\in {\mathbb{C}}^{N\times 1}$，$C\in {\mathbb{C}}^{1\times N}$。

ZOH 离散化公式：连续矩阵转化为离散矩阵：

$$\bar{A}=\exp (A\Delta )=\text{diag}(e^{A_1\Delta },\dots ,e^{A_N\Delta })\in {\mathbb{C}}^{N\times N}$$ $$\bar{B}=(\bar{A}-I)A^{-1}B\in {\mathbb{C}}^{N\times 1}$$

其中 $\Delta \in {\mathbb{R}}^{+}$ 是可学习的参数，$C$ 不参与离散。

离散递归方程 (RNN Mode): $$x_k=\bar{A}{x}_{k-1}+\bar{B}{u}_k,y_k=\text{Re}(C{x}_k)+D{u}_k$$

能控性与能观性分析

Kalman 结论保证如果连续系统是能控/能观的，那么离散化系统也是能控的，除非采样周期 $\Delta$ 满足特定的病态条件（未调研）。下面默认使用连续系统矩阵。

能控性分析

能控性矩阵：能控性矩阵为：

$$\mathcal{C}=\left[\begin{array}{ccccc}B& AB& A^{2}B& \dots & A^{n-1}B\end{array}\right]$$

能控性分析：Diagonal SSM 能控当且仅当 $A$ 的特征值互不相同且 $B$ 所有元素非零。

证明：设 $A=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$，且 $B=[b_1,b_2,\dots ,b_n{]}^T$。计算能控性矩阵 $\mathcal{C}$ 的各项：

• 第 1 列 ($B$): $[b_1,b_2,\dots ,b_n{]}^T$
• 第 2 列 ($AB$): $[{\lambda }_1{b}_1,{\lambda }_2{b}_2,\dots ,{\lambda }_n{b}_n{]}^T$
• 第 $k$ 列 ($A^{k-1}B$): $[{\lambda }_1^{k-1}{b}_1,{\lambda }_2^{k-1}{b}_2,\dots ,{\lambda }_n^{k-1}{b}_n{]}^T$

此时，能控性矩阵 $\mathcal{C}$ 可以分解为两个矩阵的乘积：

$$\mathcal{C}=\underset{\text{对角阵 }{\Sigma }_B}{\underbrace{\left[\begin{array}{cccc}{b}_1& 0& \dots & 0\\ 0& b_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_n\end{array}\right]}}\cdot \underset{\text{范德蒙德矩阵 (Vandermonde Matrix) }V}{\underbrace{\left[\begin{array}{ccccc}1& {\lambda }_1& {\lambda }_1^2& \dots & {\lambda }_1^{n-1}\\ 1& {\lambda }_2& {\lambda }_2^2& \dots & {\lambda }_2^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\lambda }_n& {\lambda }_n^2& \dots & {\lambda }_n^{n-1}\end{array}\right]}}$$

即 $\mathcal{C}=\text{diag}(B)\cdot V$。为了使 $\mathcal{C}$ 满秩 ($\det (\mathcal{C})\ne 0$)，必须满足两个条件：

• 输入必须能直接影响每一个状态分量：$\det (\text{diag}(B))\ne 0$，即 $b_i\ne 0,\forall i$。
• 范德蒙德矩阵满秩：$\det (V)\ne 0$，而 $\det (V)={\prod }_{1\le i<j\le n}({\lambda }_j-{\lambda }_i)$。这要求 $A$ 的所有特征值 ${\lambda }_i$ 必须互不相同 (Distinct Eigenvalues)。

能观性分析

能观性矩阵：

$$\mathcal{O}=\left[\begin{array}{c}C\\ CA\\ ⋮\\ C{A}^{n-1}\end{array}\right]$$

能观性分析：Diagonal SSM 能观当且仅当 $A$ 所有特征值不相同，$C$ 所有元素非零。

证明：与能控性类似。

背景与预备知识

连续系统与卷积格式

系统方程：模型由以下线性微分方程定义

$$\dot{x}(t)=Ax(t)+Bu(t),y(t)=Cx(t)$$

$A\in {\mathbb{C}}^{N\times N}$ 为系统矩阵（State Matrix），$B\in {\mathbb{R}}^{C\times 1}$ 为输入矩阵，$C\in {\mathbb{C}}^{1\times N}$ 为输出矩阵。

卷积形式 (Convolution View)：若 $x(0)=0\in {\mathbb{C}}^N$，则上述系统可写为

$$\begin{array}{r}K(t)=C{e}^{At}B,y(t)=(K\ast u)(t)\end{array}$$

其中卷积定义为 $(f\ast g)(t)={\int }_{-\infty }^{\infty }f(\tau )g(t-\tau )d\tau$。

证明：根据常数变异法可解出方程的解

$$x(t)={\int }_0^t{e}^{A(t-\tau )}Bu(\tau )\mathrm{d}\tau \Rightarrow y(t)={\int }_0^{t}C{e}^{A(t-\tau )}Bu(\tau )\mathrm{d}\tau .$$

令 $K(t)=C{e}^{At}B$，则得到

$$y(t)={\int }_0^{t}K(t-\tau )u(\tau )\mathrm{d}\tau \Rightarrow y(t)=(K\ast u)(t).$$

在对角 SSM 的情形下，$A$ 为对角矩阵。此时符号 $A_n,B_n,C_n$ 表示对应参数的元素分量。

Structured State Spaces Model (S4)

HiPPO 矩阵：为了有效地捕捉长期依赖，$A$ 矩阵被初始化为特定的 HiPPO-LegS 矩阵。其元素定义为：

$$A_{nk}=\{\begin{array}{ll}-(2n+1{)}^{1/2}(2k+1{)}^{1/2}& \text{if }n>k\\ -(n+1)& \text{if }n=k\\ 0& \text{if }n<k\end{array}$$

DPLR 结构：为了加速矩阵指数 $e^{At}$ 的计算，S4 将矩阵 $A$ 约束为对角阵加低秩矩阵 (Diagonal Plus Low-Rank, DPLR) 结构：

$$A=A^{(N)}-P{P}^T$$

• 低秩项 $P\in {\mathbb{R}}^{N\times 1}$： 其元素定义为 $P_n=(n+1/2{)}^{1/2}$。
• 正规矩阵 $A^{(N)}\in {\mathbb{R}}^{N\times N}$： 这是一个斜对称矩阵（Skew-symmetric）加上对角阵。其元素满足：

$$A_{nk}^{(N)}=\{\begin{array}{ll}-(n+1/2{)}^{1/2}(k+1/2{)}^{1/2}& \text{if }n>k\\ -1/2& \text{if }n=k\\ (n+1/2{)}^{1/2}(k+1/2{)}^{1/2}& \text{if }n<k\end{array}$$

Diagonal State Spaces (DSS)

核心结论：近期研究（由 Gupta 等人提出）发现，S4 中的低秩项部分（即 $P{P}^T$）在某些情况下可以省略，从而进一步简化模型。DSS 证明了通过对 S4 矩阵进行对角化处理，依然能够保持捕捉长程记忆的能力，同时极大地简化了代码实现和计算流程。

S4D 模型

连续系统的离散：给定步长 $\Delta t$，连续参数 $(A,B)$ 可通过以下两种常用方法转换为离散参数 $(\overline{A},\overline{B})$：

• 双线性变换 (Bilinear):

$$\overline{A}=(I-\Delta t/2A{)}^{-1}(I+\Delta t/2A),\overline{B}=(I-\Delta t/2A{)}^{-1}\cdot \Delta tB$$

• 零阶保持 (ZOH):

$$\overline{A}=\exp (\Delta tA),\overline{B}=(\Delta tA{)}^{-1}(\exp (\Delta tA)-I)\cdot \Delta tB$$

离散卷积核：离散时间 SSM 的输出为 $y=u\ast \overline{K}$，其中离散卷积核 $\overline{K}$ 具有 Krylov 子空间结构：

$$\overline{K}=(C\overline{B},C\overline{A}\overline{B},\dots ,C{\overline{A}}^{L-1}\overline{B})$$

卷积核计算：

A 参数化：为保证稳定性（防止 $t\to \infty$ 时核爆炸），通常强制 $A$ 的实部为负。$A=-\exp (A_{Re})+i\cdot A_{Im}$ （或使用 ReLU 等激活函数）。

对角状态矩阵的初始化

对角系数矩阵的稠密性

命题：任意状态空间模型 $(A,B,C)$ 等价于 $(V^{-1}AV,V^{-1}B,CV)$ 模型，其中 $V\in {\mathbb{C}}^{N\times N}$ 为可逆矩阵。

证明：考虑如下线性系统

$$\dot{x}(t)=Ax(t)+Bu(t),y(t)=Cx(t).$$

令 $x=V\ell (t)$，则上述系统等价于 $V\dot{\ell }(t)=AV\ell (t)+Bu(t)$ 以及 $y(t)=CV\ell (t)$，因此我们得到新系统

$$\dot{\ell }(t)=V^{-1}AV\ell (t)+V^{-1}Bu(t),y(t)=CV\ell (t).$$

命题：所有 ${\mathbb{C}}^{N\times N}$ 中的可对角化矩阵集合 $\mathcal{D}\subset {\mathbb{C}}^{N\times N}$ 在 ${\mathbb{C}}^{N\times N}$ 中稠密。

• 不可对角化矩阵之所以“稀有”，是因为它要求特征方程 $\det (\lambda I-A)=0$ 必须有重根，并且这些重根对应的几何重数小于代数重数。
• 上述命题证明了纯对角矩阵在数学上可以表示任何线性系统，但这并不意味着能通过梯度下降轻易地“练”出这个效果。事实上 Albert Gu 等人证明了如果使用随机初始化，稠密实矩阵和对角复矩阵的表达能力都较差。
• 数值稳定性至关重要。具有相同频谱（即等价于相同对角阵）的两个初始化可能有截然不同的性能表现。

S4D 格式推导

S4D-LegS：令 $A=A^{(N)}-P{P}^{\top }$ 为 HiPPO-LegS 矩阵。当状态维度 $N\to \infty$ 时，由正规部分 $A^{(N)}$ (即 S4D-LegS) 生成的 kernel bases $K_{A^{(N)},B/2}(t)$ 收敛于原始 S4 的 kernel bases $K_{A,B}(t)$。

S4D-Inv：为了进一步简化 S4D-LegS，我们分析 $A^{(D)}$ 的结构。实部恒定为 $\mathrm{Re}(A)=-\frac{1}{2}$，虚部遵循逆缩放律。猜测当 $N\to \infty$ 时，虚部 $\mathrm{Im}(A{)}_n\approx \Theta (n^{-1})$。

$$A_n=-\frac{1}{2}+i\frac{N}{\pi }\left(\frac{N}{2n+1}-1\right)$$

S4D-Lin：更简单的初始化 S4D-Lin 近似于傅里叶级数频率（对应 S4-FouT 变体）。其基函数 $e^{tA}B$ 表现为阻尼傅里叶基函数。

$$A_n=-\frac{1}{2}+i\pi n.$$

卷积核格式推导

递推格式展开：考虑如下离散后的 S4D 系统，其中 $t\in [0,L]$，

$$\begin{array}{rl}{h}_t& =\bar{A}{h}_{t-1}+\bar{B}{x}_t\\ y_t& =C{h}_t\end{array}$$

其中 $x_t\in \mathbb{R}$，$h_t\in {\mathbb{C}}^N$，$y_t\in {\mathbb{R}}^1$，$\bar{B}\in {\mathbb{C}}^{N\times 1}$，$\bar{A}=\text{diag}({\bar{A}}_1,\dots ,{\bar{A}}_N)\in {\mathbb{C}}^{N\times N}$，$C\in {\mathbb{C}}^{1\times N}$。我们有

$$h_t=\sum _{\tau =0}^t{\bar{A}}^{t-\tau }\bar{B}{x}_{\tau },y_t=\sum _{\tau =0}^t\bar{C}{\bar{A}}^{t-\tau }\bar{B}{x}_{\tau }.$$

卷积核定义：卷积核可以定义为如下

$$K_t=C{\bar{A}}^t\bar{B}.$$

由于 $\bar{A}$ 是对角矩阵，$\bar{B}$ 是向量，$C$ 是向量，因此 $K_t$ 是一个数，可以写为

$$K_t=\sum _{n=1}^N{C}_n{\bar{A}}_n^t{\bar{B}}_n.$$

从而

$$\mathbf{y}=\mathbf{x}\ast K\Rightarrow y_t=(\mathbf{x}\ast K{)}_t=\sum _{\tau =0}^t{x}_{\tau }\cdot K_{t-\tau }$$

其中 $K=[K_0,K_1,\cdots ,K_{L-1}]\in {\mathbb{C}}^L$ 满足如上式子。

Pytorch 实现：先根据 A, B, C, L 计算出 Kernel K (L,)。然后再转换到频域计算。

import torch
import torch.fft

def compute_ssm_kernel(A, B, C, L):
    """
    步骤 1: 生成卷积核 K
    对应公式: K_t = sum_{n=1}^N C_n * (A_n)^t * B_n
    
    参数:
        A: (N,) 复数对角矩阵的对角元素
        B: (N,) 复数输入投影
        C: (N,) 复数输出投影
        L: 序列长度
    返回:
        K: (L,) 复数卷积核
    """
    # 1. 构造时间步向量 t = [0, 1, ..., L-1]
    # shape: (L)
    t = torch.arange(L, device=A.device)
    
    # 2. 计算 A 的幂次 (A_n)^t
    # 利用广播机制: (N, 1) ** (L) -> (N, L)
    # 这一步计算了所有状态 n 在所有时刻 t 的衰减项
    A_powers = torch.pow(A.unsqueeze(-1), t)
    
    # 3. 计算各项乘积 C_n * (A_n)^t * B_n
    # term shape: (N, L)
    term = (C * B).unsqueeze(-1) * A_powers
    
    # 4. 对状态维度 N 求和 (sum_{n=1}^N)
    # 这一步将 N 个独立状态的响应混合成一个系统的脉冲响应
    K = torch.sum(term, dim=0) 
    
    return K # shape: (L,)

def fft_convolution(x, K):
    """
    步骤 2: FFT 卷积加速
    对应公式: y = x * K
    
    参数:
        x: (Batch, L) 实数输入序列
        K: (L,) 复数卷积核
    返回:
        y: (Batch, L) 实数输出序列
    """
    L = x.shape[-1]
    
    # 1. 确定 FFT 长度 (通常设为 2*L 以避免循环卷积混叠)
    fft_len = 2 * L
    
    # 2. 输入 x (实数) -> 频域
    # 使用 rfft 因为 x 是实数，只计算一半频谱
    x_f = torch.fft.rfft(x, n=fft_len)
    
    # 3. 卷积核 K (复数) -> 频域
    # 使用 fft 因为 K 本身是复数
    k_f = torch.fft.fft(K, n=fft_len)
    
    # 4. 对齐频谱
    # rfft 得到的频谱长度为 fft_len//2 + 1
    # 我们需要截取 k_f 对应的前半部分正频率
    k_f = k_f[..., :x_f.shape[-1]]
    
    # 5. 频域乘法 (对应时域卷积)
    y_f = x_f * k_f
    
    # 6. 逆变换回时域 + 截断
    y = torch.fft.irfft(y_f, n=fft_len)
    y = y[..., :L] # 去掉补零产生的部分
    
    return y

• 原论文：Arxiv
• 2025.12.10 完成笔记
• 2025.12.11 补充 SSM 经典模型的能观性和能控性，具体见 Diagonal SSM

预备知识

离散线性时不变系统 (Discrete LTI Systems)

定义(离散 LTI 系统)：设 $\mathcal{G}$ 为一个离散线性时不变 (Linear Time-Invariant, LTI) 系统，其状态方程描述如下：

$$\begin{array}{rl}h(k+1)& =Ah(k)+Bx(k),h(0)=h_0\\ y(k)& =Ch(k)+Dx(k)\end{array}$$

• $h\in {\mathbb{R}}^n$ 为状态向量 (state)；
• $x\in {\mathbb{R}}^p$ 为输入向量 (input)；
• $y\in {\mathbb{R}}^q$ 为输出向量 (output)；
• 系统矩阵维度为 $A\in {\mathbb{R}}^{n\times n},B\in {\mathbb{R}}^{n\times p},C\in {\mathbb{R}}^{q\times n},D\in {\mathbb{R}}^{q\times p}$。

假设(稳定性 Stability)：系统是稳定的，即矩阵 $A$ 的所有特征值模长均小于 1：

$$|{\lambda }_i(A)|<1,\forall i$$

假设(可控性 Controllability)：系统对 $(A,B)$ 是可控的。即存在有限时间，使得状态 $h$ 可以从任意初始状态被控制（驱动）到任意最终状态。可控性矩阵 $\mathcal{C}$ 必须满秩：

$$\mathcal{C}=\left[\begin{array}{ccccc}B& AB& A^{2}B& \dots & A^{n-1}B\end{array}\right],\text{rank}(\mathcal{C})=n$$

如果系统不能控，则状态空间中存在某些维度（子空间），无论输入 $x(k)$ 是什么，都无法影响这些维度的状态。

假设(可观性 Observability)：系统对 $(A,C)$ 是可观的。即通过观测有限时间内的输出 $y$ 和输入 $x$，足以确定初始状态 $h_0$。能观性矩阵 $\mathcal{O}$ 必须满秩：

$$\mathcal{O}=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^{n-1}\end{array}\right],\text{rank}(\mathcal{O})=n$$

如果系统不可观，意味着状态空间中存在某些维度的变化，根本不会反映在输出 $y(k)$ 中。

离散 Lyapunov 方程理论 (Discrete Lyapunov Equation)

Lyapunov 方程起源：设系统的输入 $x(k)$ 不是确定的控制量，而是零均值的高斯白噪声 $w(k)$，其协方差为 $M$（即 $E[w{w}^T]=M$）。

$$h(k+1)=Ah(k)+w(k)$$

我们想知道当时间趋于无穷时，系统状态 $h(k)$ 的分布是什么样的？定义状态的协方差矩阵为 $X(k)=E[h(k)h(k{)}^T]$。我们可以推导 $X(k)$ 随时间的变化：

$$\begin{array}{rl}X(k+1)& =E[h(k+1)h(k+1{)}^T]\\ & =E[(Ah(k)+w(k))(Ah(k)+w(k){)}^T]\\ & =E[Ah(k)h(k{)}^T{A}^T+Ah(k)w(k{)}^T+w(k)h(k{)}^T{A}^T+w(k)w(k{)}^T]\end{array}$$

由于 $h(k)$ 是过去的状态，与当前的噪声 $w(k)$ 不相关（独立），中间两项期望为 0。于是得到： $$X(k+1)=A\underset{X(k)}{\underbrace{E[h(k)h(k{)}^T]}}{A}^T+\underset{M}{\underbrace{E[w(k)w(k{)}^T]}}$$

$$X(k+1)=AX(k)A^T+M$$

如果系统是稳定的（$A$ 使得能量衰减），且噪声源源不断地输入（$M$），那么最终状态的“云团”大小（协方差）会达到一个动态平衡。此时 $X(k+1)=X(k)=X$。

定义(离散 Lyapunov 方程)：给定矩阵 $A\in {\mathbb{R}}^{n\times n}$ 和对称矩阵 $M\in {\mathbb{R}}^{n\times n}$，寻找对称矩阵 $X\in {\mathbb{R}}^{n\times n}$ 满足：

$$AX{A}^T-X+M=0$$

或者写作等价形式 $X=AX{A}^T+M$。

Lyapunov 方程有两种理解：

• 输入噪声：方程解描述了给定持续不断的随机误差扰动 $M$，系统最终会在多大的范围（方差）内波动。
• 输入控制：方程解描述了给定输入控制 $M$，对系统造成的能量大小/状态变化。

推导(级数解的构造)：我们可以通过迭代法直观地构造解。假设 $X$ 满足 $X=AX{A}^T+M$，我们将 $X$ 自身反复代入方程右边：

$$\begin{array}{rl}X& =M+AX{A}^T\\ & =M+A(M+AX{A}^T)A^T\\ & =M+AM{A}^T+A^{2}X(A^T{)}^2\\ & =M+AM{A}^T+A^2(M+AX{A}^T)(A^T{)}^2\\ & =M+AM{A}^T+A^{2}M(A^T{)}^2+A^{3}X(A^T{)}^3\\ & ⋮\\ & =\sum _{k=0}^{N-1}{A}^{k}M(A^T{)}^k+A^{N}X(A^T{)}^N\end{array}$$

当 $N\to \infty$ 时，如果余项 $A^{N}X(A^T{)}^N\to 0$，则我们得到级数解：

$$X=\sum _{k=0}^{\infty }{A}^{k}M(A^T{)}^k$$

定理(存在性与唯一性)：若矩阵 $A$ 是稳定的（即谱半径 $\rho (A)<1$），则对于任意对称矩阵 $M$：

• 级数 ${\sum }_{k=0}^{\infty }{A}^{k}M(A^T{)}^k$ 绝对收敛，且该级数是离散 Lyapunov 方程的唯一解。
• 若 $M$ 是正定 (或半正定) 的，则解 $X$ 也是正定 (或半正定) 的。

这个方程本质上是在计算一个线性系统在无限时间内的“能量累积”。只有当系统是稳定的（能量会耗散，而不是无限增长），这个累积和（级数）才是一个有限值。

可控性与可观性格拉姆矩阵 (Gramians)

Gramians 提供了比秩判据更丰富的信息：它们不仅告诉我们状态是否可控/可观，还量化了控制或观测的难易程度（能量代价）。

定义(离散可控性格拉姆矩阵)：假设上述系统是稳定且可控的，则存在唯一的对称正定矩阵 $P\in {\mathbb{R}}^{n\times n}$ 满足离散 Lyapunov 方程：

$$AP{A}^T-P+B{B}^T=0$$

$$P=\sum _{i=0}^{\infty }{A}^{i}B{B}^T(A^T{)}^i$$

物理意义：我们要寻找一个能量最小的输入序列 $x(0),x(1),\dots ,x(N)$，把系统从零状态 $h(0)=0$ 驱动到目标状态 $h_{\text{target}}$。根据系统方程 $h(k+1)=Ah(k)+Bx(k)$，经过 $N$ 步后，最终状态 $h(N)$ 可以写成输入的线性组合：

$$h(N)=\sum _{i=0}^{N-1}{A}^{N-1-i}Bx(i)$$

$$h_{\text{target}}=\underset{{\mathcal{C}}_N}{\underbrace{\left[\begin{array}{cccc}B& AB& \dots & A^{N-1}B\end{array}\right]}}\underset{\mathbf{x}}{\underbrace{\left[\begin{array}{c}x(N-1)\\ ⋮\\ x(0)\end{array}\right]}}$$

即 $h_{\text{target}}={\mathcal{C}}_N\mathbf{x}$。我们需要最小化输入能量 $E=\Vert \mathbf{x}{\Vert }^2={\mathbf{x}}^T\mathbf{x}$，同时满足约束 $h_{\text{target}}={\mathcal{C}}_N\mathbf{x}$。这是一个标准的最小范数解问题（求伪逆）。根据拉格朗日乘子法，最优解 ${\mathbf{x}}_{opt}$ 为：

$${\mathbf{x}}_{opt}={\mathcal{C}}_N^T({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-1}{h}_{\text{target}}$$

现在，我们计算这个最优输入的能量 $E_{\min }$：

$$\begin{array}{rl}{E}_{\min }& ={\mathbf{x}}_{opt}^T{\mathbf{x}}_{opt}\\ & =\left[h_{\text{target}}^T({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-T}{\mathcal{C}}_N\right]\left[{\mathcal{C}}_N^T({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-1}{h}_{\text{target}}\right]\\ & =h_{\text{target}}^T({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-1}\underset{I}{\underbrace{({\mathcal{C}}_N{\mathcal{C}}_N^T)({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-1}}}{h}_{\text{target}}\\ & =h_{\text{target}}^T({\mathcal{C}}_N{\mathcal{C}}_N^T{)}^{-1}{h}_{\text{target}}\end{array}$$

当 $N\to \infty$ 时，矩阵 ${\mathcal{C}}_N{\mathcal{C}}_N^T={\sum }_{i=0}^{\infty }{A}^{i}B{B}^T(A^T{)}^i$ 正是我们的 可控性 Gramian $P$。所以：

$$E_{\min }=h_{\text{target}}^T{P}^{-1}{h}_{\text{target}}$$

物理意义：$P$ 衡量了到达某个状态所需的最小输入能量，准确来说是 $E_{\min }=h_{\text{target}}^T{P}^{-1}{h}_{\text{target}}$，因此 $P$ 的特征值越大，所需的能量越小。

定义(离散可观性格拉姆矩阵)：假设上述系统是稳定且可观的，则存在唯一的对称正定矩阵 $Q\in {\mathbb{R}}^{n\times n}$ 满足离散 Lyapunov 方程：

$$A^{T}QA-Q+C^{T}C=0$$

其显式解形式为：

$$Q=\sum _{i=0}^{\infty }(A^T{)}^i{C}^{T}C{A}^i$$

物理意义：$Q$ 衡量了 $h_0$ 能产生多少输出能量，$J=h_0^{T}Q{h}_0$。

平衡实现 (Balanced Realizations)

定义(状态空间实现)：离散时间线性系统 $\mathcal{G}$ 由其输入-输出映射完全刻画：

$$\mathcal{G}:\{x(k){\}}_{k\ge 0}\mapsto \{y(k){\}}_{k\ge 0}$$

若四元组 $(A,B,C,D)$ 及状态 $h(k)$ 能够实现上述映射，则称其为 $\mathcal{G}$ 的一个实现 (Realization)。

注：实现不唯一。若 $(A,B,C,D)$ 是一个实现，则对于任意可逆矩阵 $T\in {\mathbb{R}}^{n\times n}$，$(T^{-1}AT,T^{-1}B,CT,D)$ 也是同一个系统的实现。

定义(最小实现 Mimimal)：若一个实现既是可控的又是可观的，则称其为最小实现。其状态维度 $n$ 称为实现的阶数 (Order)。

定义(平衡实现 Balanced)：若一个实现的格拉姆矩阵满足 $P=Q$，则称该实现是平衡的。此时记公共矩阵为 $W$，称为平衡系统的格拉姆矩阵。

Balanced realization 的意义是给系统找一把“公平的尺子”，让所有状态在“输入端”和“输出端”的重要性被同等看待。即对于某个状态 $h$，对于输入端 $P$ 和输出端 $Q$ 是同等重要的。

定理(Antoulas, 2005)：任意稳定、最小的离散 LTI 系统都存在一个平衡实现，其可控性和可观性格拉姆矩阵相等且为对角阵：

$$W=\text{diag}(\sigma )=\text{diag}({\sigma }_1,\dots ,{\sigma }_n)$$

其中 ${\sigma }_1\ge \cdots \ge {\sigma }_n>0$ 称为 Hankel 奇异值 (Hankel Singular Values, HSV)。

计算方法：HSV 也可以通过非平衡实现的格拉姆矩阵计算得出：

$$\sigma ={\text{sort}}_{\downarrow }(\sqrt{\text{spec}(PQ)})$$

Hankel 奇异值量化了每个状态的联合可控性和可观性。大的 ${\sigma }_i$ 对应既容易被控制又容易被观测的状态（即对系统动力学贡献大的状态）。

SRBT 算法：找到变换矩阵 $T$ 同时对角化 $P$ 和 $Q$。核心思想为利用 $P$ 和 $Q$ 的“平方根”（Cholesky 因子）以及 SVD 分解来构造 $T$。

• Cholesky 分解 (求平方根)：由于 $P,Q$ 是对称正定矩阵，对其进行 Cholesky 分解：

$$P=L_p{L}_p^T,Q=L_q{L}_q^T$$

其中 $L_p,L_q$ 为下三角矩阵。

• 奇异值分解 (SVD)：计算交叉乘积矩阵 $L_q^T{L}_p$ 并进行 SVD 分解

$$L_q^T{L}_p=U\Sigma V^T$$

其中 $\Sigma =\text{diag}({\sigma }_1,\dots ,{\sigma }_n)$ 即为包含 Hankel 奇异值的对角阵，$U,V$ 为正交矩阵。

• 构造变换矩阵：平衡变换矩阵 $T$ 及其逆矩阵 $T^{-1}$ 可由下式给出：

$$T=L_{p}V{\Sigma }^{-1/2}$$

$$T^{-1}={\Sigma }^{-1/2}{U}^T{L}_q^T$$

直观理解：$T$ 融合了可控性信息 ($L_p$) 和可观性信息 ($U,\Sigma$)，将系统投影到一个两者“势均力敌”的坐标系中。

平衡截断 (Balanced Truncation)

平衡截断是一种模型降阶 (Model Order Reduction, MOR) 方法，利用 Hankel 奇异值的大小对状态进行取舍。

定义(平衡截断步骤)：考虑一个稳定、最小的平衡实现 $(A,B,C,D)$，其格拉姆矩阵 $W=\text{diag}({\Sigma }_1,{\Sigma }_2)$，其中 ${\Sigma }_1$ 包含较大的 $r$ 个奇异值，${\Sigma }_2$ 包含剩余的 $n-r$ 个较小奇异值。将系统矩阵分块：

$$A=\left[\begin{array}{cc}{A}_{1,1}& A_{1,2}\\ A_{2,1}& A_{2,2}\end{array}\right],B=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right],C=\left[\begin{array}{cc}{C}_1& C_2\end{array}\right]$$

其中 $A_{1,1}\in {\mathbb{R}}^{r\times r},B_1\in {\mathbb{R}}^{r\times p},C_1\in {\mathbb{R}}^{q\times r}$。

性质(降阶系统)：由 $(A_{1,1},B_1,C_1,D)$ 定义的降阶系统 $\hat{\mathcal{G}}$ 是稳定的，且满足以下误差界：

$$||\mathcal{G}-\hat{\mathcal{G}}|{|}_{\infty }\le 2\sum _{i=r+1}^n{\sigma }_i$$

这表明截断误差由被丢弃的 Hankel 奇异值之和控制。

Hermitian 矩阵的谱稳定性 (Spectral Stability)

在训练状态空间模型 (SSMs) 时，梯度下降会逐步修改状态矩阵。了解这种扰动如何影响 Hankel 奇异值至关重要。

定理(Weyl, 1912)：设 $W$ 和 $W^{\prime }$ 为 $n$ 阶 Hermitian 矩阵（实对称矩阵），令 $\delta W=W^{\prime }-W$。设 ${\lambda }_i(W)$ 表示 $W$ 的第 $i$ 大特征值。则 ${\lambda }_i(W)$ 在 Hermitian 矩阵空间上关于算子范数是 Lipschitz 连续的：

$$|{\lambda }_i(W^{\prime })-{\lambda }_i(W)|\le \underset{i=1,\dots ,n}{\max }(|{\lambda }_i(\delta W)|)=\max (|{\lambda }_1(\delta W)|,|{\lambda }_n(\delta W)|)$$

物理意义：$W$ 的每个特征值的波动幅度不会超过扰动矩阵 $\delta W$ 的最大绝对特征值。这为训练过程中的谱变化提供了理论上界

CompreSSM 算法

CompreSSM 算法的核心思想是在训练过程中动态地对 SSM 层的线性系统进行平衡截断。该算法针对每一个 SSM 块（Block）独立进行。

算法输入：模型权重中的离散线性系统矩阵 $A,B,C$。当前系统阶数（秩）$n$。能量阈值 $\tau \in [0,1]$（预设超参数）。

Reduction Analysis 部分：求出系统中哪些状态是重要的

• 提取矩阵：从模型权重中提取 $A,B,C$。
• 求解 Gramians：解离散 Lyapunov 方程获得 $P$ 和 $Q$。$P$ 显示出哪些状态容易被输入激发（存能量），$Q$ 显示出哪些状态容易影响输出（放能量）。当 $A$ 是对角矩阵时，$P,Q$ 的无穷级数可以用等比数列求和计算。
• 计算 HSV：计算 Hankel 奇异值 $\sigma$。${\sigma }_i$ 是第 $i$ 个状态的【综合重要性得分】。

$$\sigma ={\text{sort}}_{\downarrow }(\sqrt{\text{spec}(PQ)})$$

• 确定截断阶数 $r$：找到满足总能量比例 $\tau$ 的最小阶数 $r$：

$$r=\min \left\{k\in \{1,\dots ,n\}:\sum _{i=1}^k{\sigma }_i\ge \tau \sum _{i=1}^n{\sigma }_i\right\}$$

Balanced Truncation 部分：

• 判断是否截断：若 $r<n$（即存在冗余），计算平衡变换矩阵 $T$。否则，保持系统不变。
• 执行平衡变换：将系统转换到对角平衡实现：

$$(A_b,B_b,C_b)=(T^{-1}AT,T^{-1}B,CT)$$

• 执行截断：保留前 $r$ 个维度（使用张量切片表示）：

$$(A_r,B_r,C_r)=(A_b[:r,:r],B_b[:r,:],C_b[:,:r])$$

System Replacement 部分：将模型中的原始矩阵替换为截断后的矩阵：

$$(A,B,C)\leftarrow (A_r,B_r,C_r)$$

训练中的动态降阶 (In-Training Reduction)

CompreSSM 主张在训练的早期（如学习率预热阶段）进行降阶。这依赖于以下理论和实验观察，证明了训练动力学有利于早期截断。

Hankel 奇异值的连续性

在梯度下降过程中，模型参数 $A,B,C$ 发生微小变化 $(\delta A,\delta B,\delta C)$，导致系统变为新的动力系统 $(A^{\prime },B^{\prime },C^{\prime })$。我们需要保证这种微小变化不会导致 Hankel 奇异值（即状态的重要性）发生剧烈跳变。

定义：令 $H=\sqrt{P^{1/2}Q{P}^{1/2}}$，其特征值即为 Hankel 奇异值。 $H$ 对于系统扰动是连续的，记 $H^{\prime }=H+\delta H$。

引理 (训练更新下 HSV 的连续性)： 根据 Weyl 定理，在梯度步之间，每个 Hankel 奇异值的变化幅度有上界：

$$|{\sigma }_i^{\prime }-{\sigma }_i|\le \max (|{\lambda }_{\max }(\delta H)|,|{\lambda }_{\min }(\delta H)|)=\Vert \delta H{\Vert }_2$$

这意味着状态的重要性评分 $\sigma$ 是关于模型权重的连续函数，不会突变。

相对次序的稳定性 (Stability of Relative Ordering)

仅有连续性是不够的。如果 HSV 的相对大小频繁交叉（即一个不重要的状态突然变得非常重要），那么早期的截断就是危险的。

实验观察：

• 图1 次序保持：实验表明 Hankel 奇异值的相对次序在训练初期迅速稳定。大奇异值保持较大，小奇异值保持较小。
• 图2 连续性：每步的 $\delta H$ 最大为 $O(1)$ 的量级，与 Singular Value 本身 $O(10-100)$ 的量级比起来可以算是连续的。
• 图3 次序保持：对每个维度进行追踪，次序保持稳定。
• 图4 能量贡献：底部 $r$ 个（被截断的）维度的累积能量贡献在训练过程中始终维持在低位，很少获得实质性的能量增长。

结论： 在训练早期被识别为“可忽略”的维度，通常在整个训练过程中都保持可忽略状态。因此，早期截断决策与最终的重要性排名很少冲突，这使得 In-Training Reduction 既有效又鲁棒。

实验结果

• CompreSSM 模型：从大模型开始（如 CIFAR10 初始维度 384），训练中途自动变小。
• Baseline 模型：直接从一开始就训练小模型。