笔记7-《程序员的数学3：线性代数》幂矩阵与失控判断

目录

1. 幂矩阵与失控判断
2. 为什么对角矩阵是突破口
3. 用 Ruby 验证幂矩阵
4. 数值计算上的一个提醒

幂矩阵与失控判断

• 4.2 幂矩阵：问题的提出

  • 考虑一个离散时间系统 $\mathbf{r}(t) = A\mathbf{r}(t-1)$
  • 递推展开：$\mathbf{r}(t) = A^t \mathbf{r}(0)$
  • 问题变成：$A^t$ 会不会爆炸？
  • 这是控制系统、递推算法、种群模型都会遇到的问题

• 4.3 对角矩阵的情况

  • 一维情形 $r(t) = ar(t-1)$，等价于 $r(t) = a^t r(0)$

  • 结论：

    • $|a| > 1$ → 失控
    • $|a| \le 1$ → 不失控

  • 特殊情况 $a = -1$：$r(0), -r(0), r(0), -r(0), \dots$ 数值反复跳动

    • 但不会跳到无穷大，属于「不会失控」
    • 如果不考虑这种情况而要求所有数值趋向某个目标点，术语叫渐进稳定
      • 「点 $c$ 是渐进稳定的」指：只要初始点 $\mathbf{r}(0)$ 和观测点 $c$ 之间的距离在常数 $\epsilon > 0$ 以内，当 $t \to \infty$ 时 $\mathbf{r}(t)$ 必收敛于 $c$

  • 高维情形看着吓人，实际是纸老虎

    $\mathbf{r}(t) = \begin{pmatrix} 1.3 & 0 & 0 \\ 0 & 0.2 & 0 \\ 0 & 0 & 0.8 \end{pmatrix} \mathbf{r}(t-1)$

    • 一计算就变成三个独立的标量方程：

      • $r_1(t) = 1.3 r_1(t-1)$
      • $r_2(t) = 0.2 r_2(t-1)$
      • $r_3(t) = 0.8 r_3(t-1)$

    • 每个分量各自判断是否失控即可

• 4.6 连续时间系统

  • 连续时间系统的「剧情发展」和离散系统一样，只是判断失控的条件不同
  • 背景是微分方程模型
  • 图 4.9 的 LCR 串联电路（电感 $L$、电容 $C$、电阻 $R$）
    • 无论 $L, C, R$ 取何值都不会失控
    • 真正需要担心的是加入了放大器的电路（从系统外部获得能量来源）
  • 连续时间系统对应的是 $e^{At}$ 而不是 $A^t$
    • 判据从「$|a| > 1$」变成「实部大于 0」
    • 离散系统的特征值 $\lambda$ 对应连续系统的 $\lambda$，但连续系统看的是 $\mathrm{Re}(\lambda)$

为什么对角矩阵是突破口

对角矩阵的幂很好算：$A^t$ 的对角线就是各元素自己的 $t$ 次方。所以判断「会不会爆炸」只需看对角线元素的绝对值。

但一般的 $A$ 不是对角矩阵。能不能把它变成对角的？

这正是第 4 章后半段要回答的问题：

• 4.4 可对角化的情况：$A = P D P^{-1}$，其中 $D$ 是对角矩阵

  • 于是 $A^t = P D^t P^{-1}$
  • 失控与否归结为 $D$ 的对角元素

• 4.7 不可对角化的情况：矩阵不能对角化，就没法用上面的办法

  • 这时候需要 Jordan 标准型

所以这一章的逻辑链条是：

失控问题
  → 对角矩阵能判断（简单情形）
  → 一般矩阵怎么办？
  → 能对角化的话，转化为对角问题（4.4）
  → 不能对角化的话，需要 Jordan 形（4.7）

用 Ruby 验证幂矩阵

# 矩阵乘法
def mat_mul(a, b)
  n = a.size
  m = b[0].size
  k = b.size
  n.times.map { |i| m.times.map { |j| k.times.sum { |t| a[i][t] * b[t][j] } } }
end

# 矩阵快速幂
def mat_pow(a, e)
  r = a.size
  res = Array.new(r) { |i| Array.new(r) { |j| i == j ? 1.0 : 0.0 } }
  base = a
  while e > 0
    res = mat_mul(res, base) if e & 1 == 1
    base = mat_mul(base, base)
    e >>= 1
  end
  res
end

# 系统 r(t) = A^t r(0)
a = [[1.3, 0.0], [0.0, 0.3]]
r0 = [1.0, 1.0]

[1, 2, 5, 10].each do |t|
  p = mat_pow(a, t)
  r = p.each_index.map { |i| (0...2).sum { |j| p[i][j] * r0[j] } }
  puts "t=#{t}  r = [#{r[0].round(6)}, #{r[1].round(6)}]"
end

输出：

t=1  r = [1.3, 0.3]
t=2  r = [1.69, 0.09]
t=5  r = [3.71293, 0.00243]
t=10  r = [13.785849, 6.0e-06]

两个分量的命运完全不同：第一个分量每次乘 1.3，$1.3^{10} \approx 13.79$，指数增长；第二个每次乘 0.3，$0.3^{10} \approx 6 \times 10^{-6}$，指数衰减。

整个系统失控与否，取决于最坏的那个分量。 这就是对角情形下判据「$|a| > 1$ 则失控」的直观含义。

数值计算上的一个提醒

书里 3.7 节提到，可以用消元求出的解作为初始值反复迭代，得到精度更高的结果。这叫「迭代改进」。

对幂矩阵来说，同样可以用 $A^t = (A^{t/2})^2$ 的思路。上面代码里的快速幂就是这个道理：把指数按二进制拆开，$A^{10} = A^8 \cdot A^2$，只需要 4 次矩阵乘法而不是 10 次。

矩阵乘法的复杂度是 $O(n^3)$，所以快速幂把 $O(n^3 t)$ 降到了 $O(n^3 \log t)$。当 $t$ 很大时（比如 $t = 10^6$）差别巨大。