模型降阶方法之 POD

模型降阶方法之 POD

POD 方法是 PCA 在更一般的函数空间的推广，POD 是 PCA 的一般化，PCA 是 POD 的特例。

从简单的例子说起

例子一

我们考虑一个向量值常微分方程组（ ODEs ）：
$$\{\begin{array}{l}\frac{d{u}_1}{dt}+u_1=f_1\\ \frac{d{u}_2}{u_2}+u_2=f_2\\ \frac{d{u}_3}{dt}+u_3=f_1\end{array}$$
我们认真观察可以发现，这个 ODEs 是具有冗余的，也就是说，$u_3=u_1$，我们只要求解前两个方程（约化方程）：
$$\{\begin{array}{l}\frac{d{u}_1}{dt}+u_1=f_1\\ \frac{d{u}_2}{u_2}+u_2=f_2\end{array}$$
再令 $u_3=u_1$ 就可以了，削减了计算的复杂度。那么，如何从第一个方程组到第二个方程组呢？事实上，增加一个辅助矩阵，就可以把第一个方程表达组表达为：
$$\left(\begin{array}{ll}1& 0\\ 0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{l}\frac{d{u}_1}{dt}\\ \frac{d{u}_2}{dt}\end{array}\right)+\left(\begin{array}{ll}1& 0\\ 0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{l}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{ll}1& 0\\ 0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{l}{f}_1\\ f_2\end{array}\right)$$
左右两边再左乘矩阵
$$V^T:={\left(\begin{array}{ll}1& 0\\ 0& 1\\ 1& 0\end{array}\right)}^T$$
便可得到方程：
$$\{\begin{array}{l}2\frac{d{u}_1}{dt}+2u_1=2f_1\\ \frac{d{u}_2}{u_2}+u_2=f_2\end{array}$$
这就和上述约化方程等价类。从这个例子中，我们看出，一个三阶的方程，我们经过简单的处理，就变成了二阶的方程。我们看到了模型降阶的可能性。

例子二

我们再看一个例子：
$$\left(\begin{array}{lll}1& 2& 0\\ 2& 1& 0\\ 0& 0& 3\end{array}\right)\left(\begin{array}{c}\frac{d{u}_1}{dt}\\ \frac{d{u}_2}{dt}\\ \frac{d{u}_3}{dt}\end{array}\right)+\left(\begin{array}{ccc}1& 2& 0\\ 2& 1& 0\\ 0& 0& 3\end{array}\right)\left(\begin{array}{l}{u}_1\\ u_2\\ u_3\end{array}\right)=\left(\begin{array}{c}{f}_1\\ f_2\\ f_1+f_2\end{array}\right)$$
假设，我们令 $u_3=u_1+u_2$，那么，只要前两个方程满足，第三个方程自然满足。令，
$$\left(\begin{array}{l}{u}_1\\ u_2\\ u_3\end{array}\right)=\left(\begin{array}{ll}1& 0\\ 0& 1\\ 1& 1\end{array}\right)\left(\begin{array}{l}{u}_1\\ u_2\end{array}\right):=V\left(\begin{array}{l}{u}_1\\ u_2\end{array}\right)$$
代入上式，并左乘 $V^T$ 可得：
$$\left(\begin{array}{ll}4& 5\\ 5& 4\end{array}\right)\left(\begin{array}{l}\frac{d{u}_1}{dt}\\ \frac{d{u}_2}{dt}\end{array}\right)+\left(\begin{array}{ll}4& 5\\ 5& 4\end{array}\right)\left(\begin{array}{l}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}2f_1+f_2\\ 2f_2+f_1\end{array}\right)$$
容易知道，这个方程组的求解其实就等价于原方程租的前两个方程的求解。如此，我们便把一个三阶的问题，降成了二阶的问题。

POD 的核心思想

上面的两个例子，分别满足 $u_3=u_1$ 和 $u_3=u_1+u_2$。这揭示一种共性，即虽然是解空间是一个三维的空间，但是他们的解都位于一个二维的超平面里面，也就是说，用两个变量，就可以把三维的向量值解表达出来了。从几何直观上看， $u_3=u_1$ 和 $u_3=u_1+u_2$ 反应的是 $z=x$ 和 $z=x+y$ ，分别表示的是两个三维空间中的平面，而上述方程组的解就位于这两个平面内，就给了我们足够充分的降阶理由。也就是说，我们可以用低维的超平面的一个基组去表达高维中间中的解表达进行降维，这不就是主成分分析的思想嘛。

在更高维的空间中，一般的 ODEs 的解的分布区域也不是充满整个空间的，它是位于一个特定的平面或者流形上的。由上述例子得到的启发，我们似乎可以把 ODEs 进行降阶。我们只要找到解所在的低维超平面即可。

那么问题来了，上述的两个例子很简单，我们用肉眼都能观察到方程之间的关系，从而进行降阶，对于更一般的问题，我们应该怎么处理呢？这里，我们就可以用所谓的 POD（Proper Orthogonal Decomposition）方法。

我们想要找到解所在的超平面，那么事实上，我们只要找到这个超平面上足够多个“典型”的点，那么事实上，我们就获得了平面的信息。这些离散的典型点，我们可以通过 “snapshot”（快照）的方式获得，即通过一般的数值解法，获得不同时刻的解的数值表征。

这些数值算法算出来的典型点，因为数值误差和方程固有的冗余两方面的原因，不一定都位于同一个超平面内，但是我们可以令这种 “差距” 尽可能地小，主成分分析就能很好地做到这点。事实上，这种做法在降低维度的同时刚好舍弃数值误差带来的“噪声”，所以效果很好。天下没有免费的午餐，这是一种基于已有数值解的降阶方法。

POD 的一般做法

我们来够了一下 POD 方法的一般做法。首先我们考虑一个含参的模型：

$$\mathbf{E}(\mathbf{\mu })\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}+\mathbf{A}(\mathbf{\mu })\mathbf{u}=\mathbf{f}(\mathbf{\mu },\mathbf{u})$$
其中，$\mathbf{u}\in {\mathbb{R}}^n$，$\mathbf{A}(\mathbf{\mu })\in {\mathbb{R}}^{n\times n}$，$\mathbf{E}(\mathbf{\mu })\in {\mathbb{R}}^{n\times n}$，$\mathbf{f}(\mathbf{\mu },\mathbf{u})\in {\mathbb{R}}^n$。我们先考虑右端项是线性的不含参的模型：
$$\mathbf{E}\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}+\mathbf{A}\mathbf{u}=\mathbf{b}$$

我们取 $n_s$ 个时刻的数值方法得到解（ snapshots ）每个时刻一列，把它排成一个矩阵：
$$\mathbf{U}\in {\mathbb{R}}^{n\times n_s},\mathbf{U}=\left[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_{n_s}\right]$$

我们对 $\mathbf{U}$ 做（满秩）奇异值分解如下：

$$\mathbf{U}=\mathbf{V}\mathbf{\Sigma }{\mathbf{W}}^{\mathrm{T}}$$

取 $\mathbf{X}$ 的前几列，构成矩阵
$${\mathbf{V}}_{\mathbf{r}}=\left[{\mathbf{v}}_1,\dots ,{\mathbf{v}}_r\right]$$
一般我们称之为 POD 基。很容易理解，这其实就是 PCA 的过程，参考：
https://blog.csdn.net/lusongno1/article/details/80517808
https://blog.csdn.net/lusongno1/article/details/79683461

PCA 告诉我们，POD 基表示对于投影距离平方和的极小化：
$$\begin{array}{rl}& \underset{V\in {\mathbb{R}}^{n\times r}}{\min }{\Vert \mathbf{U}-\mathbf{V}{\mathbf{V}}^{\mathrm{T}}\mathbf{U}\Vert }_F^2\\ & =\underset{V\in {\mathbb{R}}^{n\times r}}{\min }\sum _{i=1}^{n_s}{\Vert {\mathbf{u}}_i-\mathbf{V}{\mathbf{V}}^{\mathrm{T}}{\mathbf{u}}_i\Vert }_2^2=\sum _{i=r+1}^{n_s}{\sigma }_i^2\end{array}$$

有了 POD 基之后，我们令 ${\mathbf{u}}_r={{\mathbf{V}}_{\mathbf{r}}}^T\mathbf{u}$，即将 $$\mathbf{u}={\mathbf{V}}_{\mathbf{r}}{\mathbf{u}}_r$$ 代入方程，并在最后的表达式上再左乘一个 ${{\mathbf{V}}_{\mathbf{r}}}^T$，最后可以得到一个低阶的方程，
$${\mathbf{E}}_r\frac{\mathrm{d}{\mathbf{u}}_r}{\mathrm{d}t}+{\mathbf{A}}_r{\mathbf{u}}_r={\mathbf{b}}_r$$
其中，${\mathbf{E}}_r={{\mathbf{V}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{E}{\mathbf{V}}_{\mathbf{r}}$，${\mathbf{A}}_r={{\mathbf{V}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{A}{\mathbf{V}}_{\mathbf{r}}$，${\mathbf{b}}_r={{\mathbf{V}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{b}$

这样，就把一个 $n$ 阶的模型降维成了 $r$ 阶的模型，怎么样，是不是很简单呢？

关于 POD 方法的讨论

注意到，上述的过程是对右端项线性的特殊情况，如果右端是非线性的，那么， POD 方法就无法在数量级上降低计算了。譬如对于一般的非线性系统：
$$\mathbf{E}(\mathbf{\mu })\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}+\mathbf{A}(\mathbf{\mu })\mathbf{u}=\mathbf{f}(\mathbf{\mu },\mathbf{u})$$

用上述的方式，我们得到：
$${\mathbf{E}}_r(\mathbf{\mu })\frac{\mathrm{d}{\mathbf{u}}_r}{\mathrm{d}t}+{\mathbf{A}}_r(\mathbf{\mu }){\mathbf{u}}_r={\mathbf{f}}_r\left(\mathbf{\mu },{\mathbf{u}}_r\right)={{\mathbf{V}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{f}\left(\mathbf{\mu },{\mathbf{V}}_{\mathbf{r}}{\mathbf{u}}_r\right)$$
假设我们使用迭代法求解这个问题，那么我们可以看到，在每一次迭代当中，我们必须计算和 $n$ 行的量 ${\mathbf{V}}_{\mathbf{r}}$ 的乘法，以及 $n$ 维非线性函数 $\mathbf{f}$ 的计算。这个拌住了整个方程求解的计算量。为了解决这个问题，我们必须想办法搞掉 $\mathbf{f}$ 对 ${\mathbf{u}}_r$ 的直接依赖，否则就避不开非线性项的复杂性问题。

下一个主题，会介绍 DEIM ，就是为了解决非线性项的问题，持续关注。