模型降阶方法之 DEIM

POD 方法可以参考：
https://blog.csdn.net/lusongno1/article/details/125944587

POD 方法固有的局限性，使得无法很好地处理非线性项，所以必须有一个合适的方法处理非线性项。

模型问题

我们考虑的问题依然是，
$$\frac{d}{dt}\mathbf{y}(t)=\mathbf{Ay}(t)+\mathbf{F}(\mathbf{y}(t))$$
它对应的稳态问题是，
$$\mathbf{Ay}(\mu )+\mathbf{F}(\mathbf{y}(\mu ))=0$$
一些迭代法，比如牛顿法，需要雅克比矩阵。把 $(\mathbf{y}$ 看做变量，稳态问题对应的雅克比是，
$$\mathbf{J}(\mathbf{y}(\mu )):=\mathbf{A}+{\mathbf{J}}_{\mathbf{F}}(\mathbf{y}(\mu ))$$
其中，
J F ( y ( μ ) ) = diag ⁡ { F ′ ( y 1 ( μ ) ) , … , F ′ ( y n ( μ ) ) } \mathbf{J}_{\mathbf{F}}(\mathbf{y}(\mu))=\operatorname{diag}\left\{F^{\prime}\left(\mathbf{y}_{1}(\mu)\right), \ldots, F^{\prime}\left(\mathbf{y}_{n}(\mu)\right)\right\} JF​(y(μ))=diag{F′(y1​(μ)),…,F′(yn​(μ))}

POD 方法回顾

根据之前提的 POD 方法，我们取 snapshot 矩阵 $\mathbf{Y}=\left[{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{n_s}\right]$ 的奇异值分解，
$$\mathbf{Y}=\mathbf{V}\Sigma {\mathbf{W}}^T$$

取左奇异向量前 $k$ 列，

$${\mathbf{V}}_{\mathbf{k}}=\left[{\mathbf{v}}_1,\dots ,{\mathbf{v}}_k\right]$$

那么，上述的稳态问题和非稳态问题就约为了，

$$\frac{d}{dt}\tilde{\mathbf{y}}(t)=\underset{\tilde{\mathbf{A}}}{\underbrace{{\mathbf{V}}_k^T\mathbf{A}{\mathbf{V}}_k}}\tilde{\mathbf{y}}(t)+{\mathbf{V}}_k^T\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)$$

$$\underset{\tilde{\mathbf{A}}}{\underbrace{{\mathbf{V}}_k^T\mathbf{A}{\mathbf{V}}_k}}\tilde{\mathbf{y}}(\mu )+{\mathbf{V}}_k^T\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(\mu )\right)=0$$

稳态问题的雅克比为，
$$\tilde{\mathbf{J}}(\tilde{\mathbf{y}}(\mu )):=\tilde{\mathbf{A}}+{\mathbf{V}}_k^T{\mathbf{J}}_{\mathbf{F}}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(\mu )\right){\mathbf{V}}_k.$$

离散经验插值方法

离散经验插值方法（DISCRETE EMPIRICAL INTERPOLATION METHOD），简称 DEIM，是处理非线性项的一种方法。下面我们来介绍一下，它的计算流程，最后再来考虑它的原理。

如前所述，我们只要考虑非线性项，

$$\tilde{\mathbf{N}}(\tilde{\mathbf{y}}):=\underset{k\times n}{\underbrace{{\mathbf{V}}_k^T}}\underset{n\times 1}{\underbrace{\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)}}.$$

为了方便，我们统一用记号 $\mathbf{f}(\tau )$ 来表示 $\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(\mu )\right)$ 和 $\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)$。

我们企图获得非线性项的线性表示，即
$$\mathbf{f}(\tau )\approx \mathbf{Uc}(\tau )$$
其中，$\mathbf{U}=\left[{\mathbf{u}}_1,\dots ,{\mathbf{u}}_m\right]$ 是对应于 $\mathbf{f}$ 的 POD 基。为了确定 $\mathbf{c}(\tau )$ 的表达，我们不妨求解矛盾方程组，
$$\mathbf{f}(\tau )=\mathbf{Uc}(\tau )$$
求解这个矛盾方程组，我们可以玄奇其中的 $m$ 行来求解，即
$${\mathbf{P}}^T\mathbf{f}(\tau )=\left({\mathbf{P}}^T\mathbf{U}\right)\mathbf{c}(\tau )$$
其中，
$$\mathbf{P}=\left[{\mathbf{e}}_{{\wp }_1},\dots ,{\mathbf{e}}_{{\wp }_m}\right]$$

${\mathbf{e}}_j$ 表示只有第 $j$ 个位置是 1，其他都是 0 的列向量。那么，
$$\mathbf{f}(\tau )\approx \mathbf{Uc}(\tau )=\mathbf{U}{\left({\mathbf{P}}^T\mathbf{U}\right)}^{-1}{\mathbf{P}}^T\mathbf{f}(\tau )$$

这里看起来，似乎还需要把非线性项施加于每个分量算出来，其实并不然，这里的 ${\mathbf{P}}^T$ 其实就是个选行操作。我们只要计算选出的哪些行即可。

如何选定这个 $\mathbf{U}$ 的若干行去逼近这个 $\mathbf{U}$，我们可以用 greedy 的策略，也就是所谓的 DEIM 方法。

这里的 $\mathbf{P}$ （选哪些行）怎么确定？这边是 DEIM 的精髓。

如此一来，DEIM 是怎么操作的就很清楚了。为什么叫插值，为什么叫 greedy，可以看得出来，我们选的就是投影之后，距离最大的一个点，然后把它固定。固定点不就是插值么，选最大不就是贪婪么。

DEIM 误差界

容易证明，非线性项的逼近误差其实是满足，
$$\Vert \mathbf{f}(\tau )-\hat{\mathbf{f}}(\tau ){\Vert }_2\le \mathbf{C}{\mathcal{E}}_{\ast }(\tau )$$
其中，
$${\mathcal{E}}_{\ast }(\tau )={\Vert \left(\mathbf{I}-{\mathbf{UU}}^T\right)\mathbf{f}(\tau )\Vert }_2$$
$$\mathbf{C}=\Vert {\left({\mathbf{P}}^T\mathbf{U}\right)}^{-1}\Vert \le (1+\sqrt{2n}{)}^m$$

DEIM 应用细节

上面提到的 $\mathbf{P}$ 其实只是一些选行的操作，那么，容易想到
$$\begin{array}{rl}\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)& \approx \mathbf{U}{\left({\mathbf{P}}^T\mathbf{U}\right)}^{-1}{\mathbf{P}}^T\mathbf{F}\left({\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)\\ & =\mathbf{U}{\left({\mathbf{P}}^T\mathbf{U}\right)}^{-1}\mathbf{F}\left({\mathbf{P}}^T{\mathbf{V}}_k\tilde{\mathbf{y}}(t)\right)\end{array}$$

也就是说，我们对于 $ \mathbf{F}$ 选行，我们只要计算里面必要分量作用即可，也就是可以把选行操作放到 $ \mathbf{F}$ 里面去。

这么操作之后，前置项可以提前算好。那么，在每次迭代中，就没有跟 $n$ 次相关的操作了，大大缩减了计算量。

同理，对于稳态相关的问题，我们需要计算其雅克比，非线性项的雅克比矩阵为，

$${\tilde{\mathbf{J}}}_{\mathbf{F}}(\tilde{\mathbf{y}}(\mu ))\approx \underset{\text{precomputed:k }\times m}{\underbrace{{\mathbf{V}}_k^T\mathbf{U}{\left({\mathbf{P}}^T\mathbf{U}\right)}^{-1}}}\underset{m\times m}{\underbrace{{\mathbf{J}}_{\mathbf{F}}\left({\mathbf{P}}^T{\mathbf{V}}_k\tilde{\mathbf{y}}(\mu )\right)}}\underset{m\times k}{\underbrace{{\mathbf{P}}^T{\mathbf{V}}_k}},$$

其中，
J F ( P T V k y ~ ( μ ) ) = diag ⁡ { F ′ ( y 1 r ( μ ) ) , … , F ′ ( y m r ( μ ) ) } \mathbf{J}_{\mathbf{F}}\left(\mathbf{P}^{T} \mathbf{V}_{k} \tilde{\mathbf{y}}(\mu)\right)=\operatorname{diag}\left\{F^{\prime}\left(\mathbf{y}_{1}^{r}(\mu)\right), \ldots, F^{\prime}\left(\mathbf{y}_{m}^{r}(\mu)\right)\right\} JF​(PTVk​y~​(μ))=diag{F′(y1r​(μ)),…,F′(ymr​(μ))}
$${\mathbf{y}}^r(\mu )={\mathbf{P}}^T{\mathbf{V}}_k\tilde{\mathbf{y}}(\mu )$$
同样可知，它的计算量也不大。

参考

Chaturantabut S, Sorensen D C. Discrete empirical interpolation for nonlinear model reduction[C]//Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference. IEEE, 2009: 4316-4321.

Model reduction and approximation: theory and algorithms[M]. Society for Industrial and Applied Mathematics, 2017.