State Space Models

• 原论文：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 模型：直接从一开始就训练小模型。