本征正交分解

本征正交分解（proper orthogonal decomposition，POD）从一组状态快照中寻找低维正交子空间，使快照投影后的平方重建误差最小。它把多维状态写成少数空间方向与相应系数的组合，是建立降阶模型的一种基底选择方法。这里的最优性针对给定快照、内积和误差标准，并不自动保证模型预测准确。

快照的投影与丢失部分

设每个快照 $x_{j} \in ℝ^{n}$ ，共 $m$ 个，按列组成 $X = [x_{1} \dots x_{m}]$ 。选择 $r$ 个标准正交方向，组成 $U_{r} \in ℝ^{n \times r}$ ，满足 $U_{r}^{𝖳} U_{r} = I_{r}$ 。快照的系数、重建与残差分别为 $a_{j} = U_{r}^{𝖳} x_{j}, {\hat{x}}_{j} = U_{r} a_{j}, e_{j} = x_{j} - {\hat{x}}_{j} .$ POD所解的问题是 $min_{U_{r}^{𝖳} U_{r} = I_{r}} \sum_{j = 1}^{m} ‖ x_{j} - U_{r} U_{r}^{𝖳} x_{j} ‖_{2}^{2} .$ 其中“正交”指所选内积下的正交，“能量”若被使用，指这里的平方范数；只有模型赋予它相应单位和权重后，才可解释成特定物理能量。^[1]

取四个人工快照 $X = (\begin{matrix} 2 & 0 & - 2 & 0 \\ 0 & 1 & 0 & - 1 \end{matrix}) .$ 选横轴方向 $u = (1, 0)^{𝖳}$ 后，四个系数为 $2, 0, - 2, 0$ ，竖直两个快照都投影到原点。总残差平方为 $1 + 1 = 2$ 。若改选纵轴，总残差平方变成 $4 + 4 = 8$ 。

四个快照在两条坐标轴上，横轴方向的POD投影把两个竖直快照送到原点 — 虚线为两个非零投影残差。横向快照已经在所选子空间中，残差为零。

奇异值与最优误差

作快照矩阵的奇异值分解 $X = U Σ V^{𝖳}$ ，奇异值依次为 $σ_{1} \geq σ_{2} \geq \dots \geq 0$ 。前 $r$ 个左奇异向量给出最优子空间，最小残差为 $\sum_{k > r} σ_{k}^{2} .$

证明。 设 $C = X X^{𝖳}$ ，对称半正定，其特征值为 $λ_{k} = σ_{k}^{2}$ ，相应标准正交特征向量为 $u_{k}$ 。对任意候选 $W \in ℝ^{n \times r}$ ，投影的正交分解给出 $‖ X - W W^{𝖳} X ‖_{F}^{2} = ‖ X ‖_{F}^{2} - tr (W^{𝖳} C W) .$ 令 $b_{k} = ‖ W^{𝖳} u_{k} ‖_{2}^{2}$ ，则 $0 \leq b_{k} \leq 1$ 且 $\sum_{k} b_{k} = r$ 。保留部分为 $\sum_{k} λ_{k} b_{k}$ 。当 $1 \leq r < n$ 时， $\begin{matrix} \sum_{k \leq r} λ_{k} - \sum_{k} λ_{k} b_{k} \\ \geq λ_{r} (\sum_{k \leq r} (1 - b_{k}) - \sum_{k > r} b_{k}) \\ = 0. \end{matrix}$ 不等式使用前 $r$ 个特征值不小于 $λ_{r}$ ，其余不大于它；括号因权重总和而为零。取 $W = (u_{1}, \dots, u_{r})$ 达到等号，故误差等于剩余特征值之和。 $r = 0$ 与 $r = n$ 分别是全部丢弃和完整保留，直接验证即可。

在本例， $X X^{𝖳} = diag (8, 2)$ ，奇异值为 $2 \sqrt{2}, \sqrt{2}$ 。保留一个方向时的平方范数比例为 $8 / (8 + 2) = 80 %$ ，误差平方为二。

两个POD方向的奇异值平方分别为8和2 — 保留第一方向留下80%的快照平方范数。该比例衡量当前数据的重建，不是未来预测的准确率。

若截断处 $σ_{r} = σ_{r + 1}$ ，最优子空间可能不唯一。若数据矩阵为 $I_{2}$ 而只保留一维，任何过原点的直线都给出相同误差一，这提供了一个可直接复算的变式。

权重、均值与物理内积

时间采样不均匀时，可以给快照正权重 $w_{j}$ ，最小化 $\sum_{j} w_{j} ‖ e_{j} ‖^{2}$ ，等价于对列被乘以 $\sqrt{w_{j}}$ 的矩阵做普通POD。本例若权重依次为 $1, 4, 1, 4$ ，加权相关矩阵变成 $diag (8, 8)$ ，不再有唯一的首方向。改变权重确实改变了所教给模型的数据分布。

若离散物理场的范数为 $‖ x ‖_{M}^{2} = x^{𝖳} M x$ ，其中 $M$ 是给定对称正定质量矩阵，应先对 $M^{1 / 2} X$ （以及相应快照权重）做SVD，再把方向乘 $M^{- 1 / 2}$ 变回原坐标。所得基满足 $U_{r}^{𝖳} M U_{r} = I$ 。直接忽略网格体积或不同单位，会使“主要方向”的含义改变。

主成分分析通常先减去样本均值。POD也可以针对均值之外的波动做分解，但必须说明是否中心化；若中心化，重建应写成 $\bar{x} + U_{r} a_{j}$ 。未中心化POD近似经过原点的子空间，中心化后则近似经过均值的仿射空间，两者处理的是不同问题。

Python重建

下面按欧氏内积处理未中心化的有限实快照矩阵，返回基、投影系数、全部奇异值及残差平方，不改写输入。 $r = 0$ 被允许，表示重建为零。

Python 3

import numpy as np

def pod(X, r):
    X = np.asarray(X, dtype=float)
    if X.ndim != 2 or min(X.shape) == 0 or not np.isfinite(X).all():
        raise ValueError("require a finite nonempty matrix")
    if type(r) is not int or not 0 <= r <= min(X.shape):
        raise ValueError("rank outside allowed range")
    U, s, _ = np.linalg.svd(X, full_matrices=False)
    basis = U[:, :r]
    coefficients = basis.T @ X
    residual = X - basis @ coefficients
    return basis, coefficients, s, float(np.sum(residual ** 2))

X = np.array([[2, 0, -2, 0], [0, 1, 0, -1]])
basis, coefficients, singular_values, error = pod(X, 1)
print(singular_values ** 2)  # [8., 2.]
print(error)               # 2.0
assert np.allclose(basis.T @ basis, np.eye(1))
assert np.isclose(error, np.sum(singular_values[1:] ** 2))
assert np.isclose(pod(X, 0)[3], 10)
assert np.isclose(pod(X, 2)[3], 0)

稠密SVD时间为 $O (n m min (n, m))$ ，还需 $O (n m r)$ 的投影与重建运算；存储规模至少包括 $n m$ 个快照元素和SVD工作数组。实际大规模实现可用快照相关矩阵或迭代SVD，具体数值稳定性另需检查。

POD仅选出重建方向，尚未给出系数怎样随时间变化。要预测，仍需投影原方程或拟合系数演化模型。快照没有覆盖的参数、外力和罕见状态可能被丢失；小的训练重建误差也不能单独证明降阶系统稳定。

参考文献

↑ Carmen Gräßle、Michael Hinze、Stefan Volkwein：Model Order Reduction by Proper Orthogonal Decomposition，§2，有限维快照与POD最优化问题。

[1] Carmen Gräßle、Michael Hinze、Stefan Volkwein：Model Order Reduction by Proper Orthogonal Decomposition，§2，有限维快照与POD最优化问题。

[1]