POD 方法的介绍如下所示:
https://blog.csdn.net/lusongno1/article/details/125944587
经典的 POD 方法无法处理非线性项,所以发展了很多的方法,比如 DEIM:
https://blog.csdn.net/lusongno1/article/details/125955245
如果我们要考虑的方程里面有一组参数,不同的参数都能获得一个 Matrix,那么我们应该如何把多组参数下的得到的数据放在一块作为先验,来给新到来的一组参数做模型降阶呢?这就是这篇文章要考虑的问题,即高阶的张量分解,和参数维度上的插值。简单地说,就是如果已经有了一组参数下的快照,那么新来一个参数,它之下的 SVD 分解不需要重新找快照,只需要由原来的参数插值得到。
这篇博客有点难理解,可以大概看看,再去看一些文章。我们略去一些约定俗成的符号的含义的解释,不理解的符号含义可以看一下参考文献。我们只介绍方法,证明分析和实验,可以阅读论文。
我们先介绍一些基础的定义。一个代数张量空间中的张量可以写为 elementary 张量的线性组合:
v
=
∑
i
=
1
m
v
i
(
1
)
⊗
⋯
⊗
v
i
(
d
)
v=\sum_{i=1}^{m} v_{i}^{(1)} \otimes \cdots \otimes v_{i}^{(d)}
v=i=1∑mvi(1)⊗⋯⊗vi(d)
注意到这里其实是没有组合系数的。
一个已知张量的最小子空间,它是一组空间。它是 v ∈ ⨂ v = 1 d U v v \in \bigotimes_{v=1}^{d} U_{v} v∈⨂v=1dUv 中, U ν U_{\nu} Uν 的极小化。也就说,能够包住 v v v 的最最小的空间分量。
一个二阶张量的秩,是最小的
r
r
r 使得满足:
u
=
∑
i
=
1
r
s
i
⊗
v
i
u=\sum_{i=1}^{r} s_{i} \otimes v_{i}
u=i=1∑rsi⊗vi
对于二阶张量,奇异值分解一般是这么写的:
u
=
∑
i
=
1
n
σ
i
s
i
⊗
v
i
u=\sum_{i=1}^{n} \sigma_{i} s_{i} \otimes v_{i}
u=i=1∑nσisi⊗vi
对于高阶张量,它的 Canonical 秩就是二阶张量秩的推广,即最小的
r
r
r 使得,
v
=
∑
i
=
1
r
v
i
(
1
)
⊗
⋯
⊗
v
i
(
d
)
v=\sum_{i=1}^{r} v_{i}^{(1)} \otimes \cdots \otimes v_{i}^{(d)}
v=i=1∑rvi(1)⊗⋯⊗vi(d)
所谓的
α
\alpha
α 值得是最小子空间的维数,
α
\alpha
α 表示索引,
rank
α
(
v
)
=
dim
(
U
α
min
(
v
)
)
.
\operatorname{rank}_{\alpha}(v)=\operatorname{dim}\left(U_{\alpha}^{\min }(v)\right) .
rankα(v)=dim(Uαmin(v)).
Tucker Rank 其实是一个向量,它是用最小子空间的维数来定义的。即一个具有 Tucker Rank
r
=
(
r
ν
)
v
∈
D
r=\left(r_{\nu}\right)_{v \in D}
r=(rν)v∈D 的张量集合,它满足在所有的指标集合里,对应的最小子空间的维数是限定的,即:
T
r
=
{
v
∈
X
:
rank
v
(
v
)
=
dim
(
U
v
min
(
v
)
)
≤
r
ν
,
v
∈
D
}
,
\mathscr{T}_{r}=\left\{v \in X: \operatorname{rank}_{v}(v)=\operatorname{dim}\left(U_{v}{ }^{\min }(v)\right) \leq r_{\nu}, v \in D\right\},
Tr={v∈X:rankv(v)=dim(Uvmin(v))≤rν,v∈D},
那么 Tucker 张量里面的一个元素就可以写成:
v
=
∑
i
1
=
1
r
1
…
∑
i
d
=
1
r
d
C
i
1
,
…
,
i
d
v
i
1
(
1
)
⊗
⋯
⊗
v
i
d
(
d
)
v=\sum_{i_{1}=1}^{r_{1}} \ldots \sum_{i_{d}=1}^{r_{d}} C_{i_{1}, \ldots, i_{d}} v_{i_{1}}^{(1)} \otimes \cdots \otimes v_{i_{d}}^{(d)}
v=i1=1∑r1…id=1∑rdCi1,…,idvi1(1)⊗⋯⊗vid(d)
Tree-Based Rank 指的是,我们现在不限定单个的指标对应的分量最小子空间的维数,我们限定给定指标集的对应的张量分解空间的张成的空间的维数,即
B
T
r
=
{
v
∈
X
:
rank
α
(
v
)
=
dim
(
U
α
min
(
v
)
)
≤
r
α
,
α
∈
T
D
}
.
\mathscr{B} \mathscr{T}_{r}=\left\{v \in X: \operatorname{rank}_{\alpha}(v)=\operatorname{dim}\left(U_{\alpha}^{\min }(v)\right) \leq r_{\alpha}, \alpha \in T_{D}\right\} .
BTr={v∈X:rankα(v)=dim(Uαmin(v))≤rα,α∈TD}.
注意到这里,一个指标集合,就给一个限定,有多少个指标集合,就给多少个限定。Tucker Rank 是它去每个指标集合单个指标的时候的特例。
TT(Tensor-train) Rank 也是 Tree-Based Rank 的一种特例。它取得指标集是有序的,从尾巴开始增长的一组指标集合,即
T
T
r
=
{
v
∈
X
:
rank
{
k
+
1
,
…
,
d
}
(
v
)
≤
r
k
}
\mathscr{T} \mathscr{T}_{r}=\left\{v \in X: \operatorname{rank}_{\{k+1, \ldots, d\}}(v) \leq r_{k}\right\}
TTr={v∈X:rank{k+1,…,d}(v)≤rk}
下面简单介绍一下高阶奇异值分解(HOSVD)。给定一个指标集合
α
\alpha
α,可以把 vector 空间分为两部分。一般的截断 SVD 可以写为:
u
α
,
r
α
=
∑
i
=
1
r
α
σ
i
(
α
)
u
i
(
α
)
⊗
u
i
(
α
c
)
u_{\alpha, r_{\alpha}}=\sum_{i=1}^{r_{\alpha}} \sigma_{i}^{(\alpha)} u_{i}^{(\alpha)} \otimes u_{i}^{\left(\alpha^{c}\right)}
uα,rα=∑i=1rασi(α)ui(α)⊗ui(αc)
HOSVD 在 Tucker 格式和 tree-based Tucker 格式,都是基于投影来实现的,具体就不说了。
假设我们求解的问题如下,
u
t
=
F
(
t
,
u
,
α
)
,
t
∈
(
0
,
T
)
,
and
u
∣
t
=
0
=
u
0
\mathbf{u}_{t}=F(t, \mathbf{u}, \boldsymbol{\alpha}), \quad t \in(0, T), \quad \text { and }\left.\mathbf{u}\right|_{t=0}=\mathbf{u}_{0}
ut=F(t,u,α),t∈(0,T), and u∣t=0=u0
我们把这个问题的解一列一列地排成一个矩阵,叫做 snapshots,
Φ
pod
(
α
)
=
[
ϕ
1
(
α
)
,
…
,
ϕ
N
(
α
)
]
∈
R
M
×
N
\Phi_{\text {pod }}(\boldsymbol{\alpha})=\left[\boldsymbol{\phi}_{1}(\boldsymbol{\alpha}), \ldots, \boldsymbol{\phi}_{N}(\boldsymbol{\alpha})\right] \in \mathbb{R}^{M \times N}
Φpod (α)=[ϕ1(α),…,ϕN(α)]∈RM×N
对 snapshots 矩阵做奇异值分解,
Φ
p
o
d
(
α
)
=
U
Σ
V
T
\Phi_{\mathrm{pod}}(\boldsymbol{\alpha})=\mathrm{U \Sigma V}^{T}
Φpod(α)=UΣVT
取前几个左奇异向量,作为子空间的一组基,我们叫做 POD 基。
现在我们来处理,假设上述的参数不是一个,而是有很多很多个,我们又应该如何做 SVD 分解,又应该如何得到 POD 基呢?
我们不妨假设,上述的参数是位于一个高维的盒子里面,维数是
D
D
D,
A
=
⨂
i
=
1
D
[
α
i
min
,
α
i
max
]
\mathcal{A}=\bigotimes_{i=1}^{D}\left[\alpha_{i}^{\min }, \alpha_{i}^{\max }\right]
A=i=1⨂D[αimin,αimax]
我们可以把这个高维的盒子在各个方向上都切上几刀,就得到了一些节点集合,
A
^
=
{
α
^
=
(
α
^
1
,
…
,
α
^
D
)
T
:
α
^
i
∈
{
α
^
i
j
}
j
=
1
,
…
,
n
i
,
i
=
1
,
…
,
D
}
\widehat{\mathcal{A}}=\left\{\widehat{\boldsymbol{\alpha}}=\left(\widehat{\alpha}_{1}, \ldots, \widehat{\alpha}_{D}\right)^{T}: \widehat{\alpha}_{i} \in\left\{\widehat{\alpha}_{i}^{j}\right\}_{j=1, \ldots, n_{i}}, i=1, \ldots, D\right\}
A
={α
=(α
1,…,α
D)T:α
i∈{α
ij}j=1,…,ni,i=1,…,D}
总的节点数目为,
K
=
∏
i
=
1
D
n
i
K=\prod_{i=1}^{D} n_{i}
K=i=1∏Dni
每一个点,就是一个参数,每一个参数,都可以做 snapshots 得到一个矩阵,这些所有矩阵想办法组织在一起,就形成了高阶张量,
(
Φ
)
:
,
i
1
,
…
,
i
D
,
k
=
ϕ
k
(
α
^
1
i
1
,
…
,
α
^
D
i
D
)
(\boldsymbol{\Phi})_{:, i_{1}, \ldots, i_{D}, k}=\boldsymbol{\phi}_{k}\left(\widehat{\alpha}_{1}^{i_{1}}, \ldots, \widehat{\alpha}_{D}^{i_{D}}\right)
(Φ):,i1,…,iD,k=ϕk(α
1i1,…,α
DiD)
张量的 Frobenius 范数的定义和矩阵是类似的定义方式,
∥
Φ
∥
F
:
=
(
∑
j
=
1
M
∑
i
1
=
1
n
1
…
∑
i
D
=
1
n
D
∑
k
=
1
N
Φ
j
,
i
1
,
…
,
i
D
,
k
2
)
1
/
2
\|\boldsymbol{\Phi}\|_{F}:=\left(\sum_{j=1}^{M} \sum_{i_{1}=1}^{n_{1}} \ldots \sum_{i_{D}=1}^{n_{D}} \sum_{k=1}^{N} \boldsymbol{\Phi}_{j, i_{1}, \ldots, i_{D}, k}^{2}\right)^{1 / 2}
∥Φ∥F:=(j=1∑Mi1=1∑n1…iD=1∑nDk=1∑NΦj,i1,…,iD,k2)1/2
在这个定义之下,我们要假设张量的低阶逼近
Φ
~
\widetilde{\boldsymbol{\Phi}}
Φ
要满足,
∥
Φ
~
−
Φ
∥
F
≤
ε
~
∥
Φ
∥
F
\|\tilde{\boldsymbol{\Phi}}-\boldsymbol{\Phi}\|_{F} \leq \widetilde{\varepsilon}\|\boldsymbol{\Phi}\|_{F}
∥Φ~−Φ∥F≤ε
∥Φ∥F
对于一个张量而言,我们可以定义某一个方向上的矩阵乘向量操作如下:
(
Ψ
×
k
a
)
j
1
,
…
,
j
k
−
1
,
j
k
+
1
,
…
,
j
m
=
∑
j
k
=
1
N
k
Ψ
j
1
,
…
,
j
m
a
j
k
\left(\boldsymbol{\Psi} \times_{k} \mathbf{a}\right)_{j_{1}, \ldots, j_{k-1}, j_{k+1}, \ldots, j_{m}}=\sum_{j_{k}=1}^{N_{k}} \boldsymbol{\Psi}_{j_{1}, \ldots, j_{m}} a_{j_{k}}
(Ψ×ka)j1,…,jk−1,jk+1,…,jm=jk=1∑NkΨj1,…,jmajk
$ \times_{k}$ 表示在第 k k k 个方向上做线性组合压缩。
我们可以定义类似于 one-hot 编码的算子
e
i
(
α
^
)
=
(
e
1
i
(
α
^
)
,
…
,
e
n
i
i
(
α
^
)
)
T
∈
R
n
i
,
i
=
1
,
…
,
D
\mathbf{e}^{i}(\widehat{\boldsymbol{\alpha}})=\left(e_{1}^{i}(\widehat{\boldsymbol{\alpha}}), \ldots, e_{n_{i}}^{i}(\widehat{\boldsymbol{\alpha}})\right)^{T} \in \mathbb{R}^{n_{i}}, i=1, \ldots, D
ei(α
)=(e1i(α
),…,enii(α
))T∈Rni,i=1,…,D 为,
e
j
i
(
α
^
)
=
{
1
if
α
^
i
=
α
^
i
j
0
otherwise
,
j
=
1
,
…
,
n
i
e_{j}^{i}(\widehat{\boldsymbol{\alpha}})=\left\{1 if ˆαi=ˆαji0 otherwise
事实上,
e
i
\mathbf{e}^{i}
ei 表示的是,对于
α
^
=
(
α
^
1
,
…
,
α
^
D
)
T
\widehat{\boldsymbol{\alpha}}=\left(\widehat{\alpha}_{1}, \ldots, \widehat{\alpha}_{D}\right)^{T}
α
=(α
1,…,α
D)T 的第
i
i
i 个维度,它的取值是第几个,那么
e
i
\mathbf{e}^{i}
ei 的第几个位置就为 1 ,其他为零。
有了上面的符号定义,我们就可以提取每个参数节点上的 snapshot 如下所示:
Φ
e
(
α
^
)
=
Φ
×
2
e
1
(
α
^
)
×
3
e
2
(
α
^
)
⋯
×
D
+
1
e
D
(
α
^
)
∈
R
M
×
N
\Phi_{e}(\widehat{\boldsymbol{\alpha}})=\boldsymbol{\Phi} \times_{2} \mathbf{e}^{1}(\widehat{\boldsymbol{\alpha}}) \times_{3} \mathbf{e}^{2}(\widehat{\boldsymbol{\alpha}}) \cdots \times_{D+1} \mathbf{e}^{D}(\widehat{\boldsymbol{\alpha}}) \in \mathbb{R}^{M \times N}
Φe(α
)=Φ×2e1(α
)×3e2(α
)⋯×D+1eD(α
)∈RM×N
同样地,对于逼近张量,我们也可以用相同的方式进行提取,
Φ
~
e
(
α
^
)
=
Φ
~
×
2
e
1
(
α
^
)
×
3
e
2
(
α
^
)
⋯
×
D
+
1
e
D
(
α
^
)
∈
R
M
×
N
\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})=\widetilde{\boldsymbol{\Phi}} \times_{2} \mathbf{e}^{1}(\widehat{\boldsymbol{\alpha}}) \times_{3} \mathbf{e}^{2}(\widehat{\boldsymbol{\alpha}}) \cdots \times_{D+1} \mathbf{e}^{D}(\widehat{\boldsymbol{\alpha}}) \in \mathbb{R}^{M \times N}
Φ
e(α
)=Φ
×2e1(α
)×3e2(α
)⋯×D+1eD(α
)∈RM×N
若
N
~
=
rank
(
Φ
~
e
(
α
^
)
)
\tilde{N}=\operatorname{rank}\left(\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\right)
N~=rank(Φ
e(α
)),我们有如下估计式成立,
∑
i
=
1
N
∥
ϕ
i
−
∑
j
=
1
N
~
⟨
ϕ
i
,
z
j
⟩
z
j
∥
ℓ
2
2
≤
ε
~
2
∥
Φ
∥
F
2
\sum_{i=1}^{N}\left\|\phi_{i}-\sum_{j=1}^{\tilde{N}}\left\langle\phi_{i}, \mathbf{z}_{j}\right\rangle \mathbf{z}_{j}\right\|_{\ell^{2}}^{2} \leq \widetilde{\varepsilon}^{2}\|\boldsymbol{\Phi}\|_{F}^{2}
i=1∑N∥
∥ϕi−j=1∑N~⟨ϕi,zj⟩zj∥
∥ℓ22≤ε
2∥Φ∥F2
类似于上面
e
i
\mathbf{e}^{i}
ei 的定义方式,我们也可以将其定义为参数盒子第
i
i
i 个维度上的拉格朗日插值节点基函数,
e
i
:
α
→
R
n
i
,
i
=
1
,
…
,
D
\mathbf{e}^{i}: \boldsymbol{\alpha} \rightarrow \mathbb{R}^{n_{i}}, \quad i=1, \ldots, D
ei:α→Rni,i=1,…,D
e
j
i
(
α
)
=
{
∏
m
=
1
,
m
≠
k
p
(
α
^
i
i
m
−
α
i
)
/
∏
m
=
1
,
m
≠
k
p
(
α
^
i
i
m
−
α
^
i
j
)
,
if
j
=
i
k
∈
{
i
1
,
…
,
i
p
}
,
0
,
otherwise
,
e_{j}^{i}(\boldsymbol{\alpha})= {∏pm=1,m≠k(ˆαimi−αi)/∏pm=1,m≠k(ˆαimi−ˆαji), if j=ik∈{i1,…,ip},0, otherwise ,
由多项式插值的定义,我们就可以讲一个函数在一个插值节点上进行插值如下,
g
(
α
i
)
≈
∑
j
=
1
n
i
e
j
i
(
α
)
g
(
α
^
i
j
)
,
i
=
1
,
…
,
D
g\left(\alpha_{i}\right) \approx \sum_{j=1}^{n_{i}} e_{j}^{i}(\boldsymbol{\alpha}) g\left(\widehat{\alpha}_{i}^{j}\right), \quad i=1, \ldots, D
g(αi)≈j=1∑nieji(α)g(α
ij),i=1,…,D
在这些定义之下,对于高维张量,我们就可以通过插值的方式,把参数对应的那些维度抹成一个量,这个量是个函数,
Φ
~
e
(
α
)
=
Φ
~
×
2
e
1
(
α
)
×
3
e
2
(
α
)
⋯
×
D
+
1
e
D
(
α
)
∈
R
M
×
N
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\tilde{\boldsymbol{\Phi}} \times_{2} \mathbf{e}^{1}(\boldsymbol{\alpha}) \times_{3} \mathbf{e}^{2}(\boldsymbol{\alpha}) \cdots \times_{D+1} \mathbf{e}^{D}(\boldsymbol{\alpha}) \in \mathbb{R}^{M \times N}
Φ
e(α)=Φ~×2e1(α)×3e2(α)⋯×D+1eD(α)∈RM×N
当然,如果我们对参数盒子,只是在每个维度上进行部分地踩点,譬如,
A
^
p
:
=
{
α
^
=
(
α
^
1
,
…
,
α
^
D
)
T
:
α
^
i
∈
{
α
^
i
j
}
j
∈
{
i
1
,
…
,
i
p
}
,
i
=
1
,
…
,
D
}
⊂
A
^
\widehat{\mathcal{A}}_{p}:=\left\{\widehat{\boldsymbol{\alpha}}=\left(\widehat{\alpha}_{1}, \ldots, \widehat{\alpha}_{D}\right)^{T}: \widehat{\alpha}_{i} \in\left\{\widehat{\alpha}_{i}^{j}\right\}_{j \in\left\{i_{1}, \ldots, i_{p}\right\}}, i=1, \ldots, D\right\} \subset \widehat{\mathcal{A}}
A
p:={α
=(α
1,…,α
D)T:α
i∈{α
ij}j∈{i1,…,ip},i=1,…,D}⊂A
也可以有类似的结果。具有这些准备之后,下面我来介绍一下三种格式之下的 POD 基的选取做法。不再赘述,只力图把算法清楚,以能够编程为目标。事实上,这三种方式的核心思想都是一样的,即我们先对张量快照进行某些方式下的张量分解,然后对参数部分的元张量做插值,得到中心矩阵,最后对中心矩阵做奇异值分解。选取分解后左奇异向量的前
n
n
n 个。
canonical polyadic (CP)分解是最基本的一个张量分解方式。这种分解表述下的张量逼近,简单来说是如下一个过程。
离线阶段:
首先,我们肯定是通过一些传统的数值方法求解问题,得到解的一个张量式的 snapshots(为了中化,下面不妨就叫做快照)。对于这个快照张量,给定目标 canonical 秩
R
R
R,我们使用使用交替最小二乘法算法(ALS)去得到它的低秩 CP 分解,
Φ
≈
Φ
~
=
∑
r
=
1
R
u
r
∘
σ
1
r
∘
⋯
∘
σ
D
r
∘
v
r
\boldsymbol{\Phi} \approx \widetilde{\boldsymbol{\Phi}}=\sum_{r=1}^{R} \mathbf{u}^{r} \circ \boldsymbol{\sigma}_{1}^{r} \circ \cdots \circ \boldsymbol{\sigma}_{D}^{r} \circ \mathbf{v}^{r}
Φ≈Φ
=r=1∑Rur∘σ1r∘⋯∘σDr∘vr
ALS 算法参考:Harshman R A. Foundations of the PARAFAC procedure: Models and conditions for an" explanatory" multimodal factor analysis[J]. 1970.
其次,我们把得到的分解的
u
r
\mathbf{u}^{r}
ur 和
V
r
\mathbf{V}^{r}
Vr 按列拼接成两个矩阵,
U
^
=
[
u
1
,
…
,
u
R
]
∈
R
M
×
R
,
V
^
=
[
v
1
,
…
,
v
R
]
∈
R
N
×
R
\widehat{\mathrm{U}}=\left[\mathbf{u}^{1}, \ldots, \mathbf{u}^{R}\right] \in \mathbb{R}^{M \times R}, \quad \widehat{\mathrm{V}}=\left[\mathbf{v}^{1}, \ldots, \mathbf{v}^{R}\right] \in \mathbb{R}^{N \times R}
U
=[u1,…,uR]∈RM×R,V
=[v1,…,vR]∈RN×R
最后,我们把得到的这两个矩阵分别做 QR 分解,
U
^
=
U
R
U
,
V
^
=
V
R
V
\widehat{\mathrm{U}}=\mathrm{UR}_{U}, \quad \widehat{\mathrm{V}}=\mathrm{VR}_{V}
U
=URU,V
=VRV
在线阶段:
指定一个约化空间维数
n
≤
R
n \leq R
n≤R,以及给定一组参数向量
α
∈
A
\alpha \in \mathcal{A}
α∈A。
首先,我们可以把张量在各个维度上在点
α
∈
A
\alpha \in \mathcal{A}
α∈A 处做插值,对于 CP 分解来说,事实上就是各个维度上的节点基函数在
α
∈
A
\alpha \in \mathcal{A}
α∈A 取值,并和 elementary 张量做内积,变成一些数并做积,即
Φ
~
e
(
α
)
=
∑
r
=
1
R
s
r
u
r
∘
v
r
∈
R
M
×
N
,
with
s
r
=
∏
i
=
1
D
⟨
σ
i
r
,
e
i
(
α
)
⟩
∈
R
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\sum_{r=1}^{R} s_{r} \mathbf{u}^{r} \circ \mathbf{v}^{r} \in \mathbb{R}^{M \times N}, \text { with } s_{r}=\prod_{i=1}^{D}\left\langle\boldsymbol{\sigma}_{i}^{r}, \mathbf{e}^{i}(\boldsymbol{\alpha})\right\rangle \in \mathbb{R}
Φ
e(α)=r=1∑Rsrur∘vr∈RM×N, with sr=i=1∏D⟨σir,ei(α)⟩∈R
其次,把这些
s
r
s_r
sr 张成对角阵
S
(
α
)
=
diag
(
s
1
,
…
,
s
R
)
\mathrm{S}(\boldsymbol{\alpha})=\operatorname{diag}\left(s_{1}, \ldots, s_{R}\right)
S(α)=diag(s1,…,sR),并构造 core Matrix 如下,
C
(
α
)
=
R
U
S
(
α
)
R
V
T
\mathrm{C}(\boldsymbol{\alpha})=\mathrm{R}_{U} \mathrm{~S}(\boldsymbol{\alpha}) \mathrm{R}_{V}^{T}
C(α)=RU S(α)RVT
接着,对 core Matrix 做奇异值分解,
C
(
α
)
=
U
c
Σ
c
V
c
T
\mathrm{C}(\boldsymbol{\alpha})=\mathrm{U}_{c} \Sigma_{c} \mathrm{~V}_{c}^{T}
C(α)=UcΣc VcT
最后,我们就得到了这个在线参数下的奇异值分解,
Φ
~
e
(
α
)
=
U
C
(
α
)
V
T
=
(
U
U
c
)
Σ
c
(
V
V
c
)
T
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\mathrm{UC}(\boldsymbol{\alpha}) \mathrm{V}^{T}=\left(\mathrm{UU}_{c}\right) \Sigma_{c}\left(\mathrm{VV}_{c}\right)^{T}
Φ
e(α)=UC(α)VT=(UUc)Σc(VVc)T
那么,
{
β
i
(
α
)
}
i
=
1
n
⊂
R
R
\left\{\boldsymbol{\beta}_{i}(\boldsymbol{\alpha})\right\}_{i=1}^{n} \subset \mathbb{R}^{R}
{βi(α)}i=1n⊂RR 就可以直接取为
U
c
\mathrm{U}_{c}
Uc 的前
n
n
n 列。
高阶奇异值分解(higher order SVD ,HOSVD)比起 CP 分解,能有一个保证极小的性质。即满足,
∥
Φ
−
Φ
~
∥
≤
D
+
2
∥
Φ
−
Φ
o
p
t
∥
and
∥
Φ
−
Φ
~
∥
≤
(
∑
i
=
1
D
+
1
Δ
i
2
)
1
2
\|\boldsymbol{\Phi}-\tilde{\boldsymbol{\Phi}}\| \leq \sqrt{D+2}\left\|\boldsymbol{\Phi}-\boldsymbol{\Phi}^{\mathrm{opt}}\right\| \quad \text { and } \quad\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\| \leq\left(\sum_{i=1}^{D+1} \Delta_{i}^{2}\right)^{\frac{1}{2}}
∥Φ−Φ~∥≤D+2∥
∥Φ−Φopt∥
∥ and ∥Φ−Φ
∥≤(i=1∑D+1Δi2)21
离线阶段:
我们给定张量的快照和目标的精度
ε
~
\widetilde{\varepsilon}
ε
。
离线阶段就是计算 HOSVD 分解,得到分解矩阵和系数(core tensor),可以用的算法有很多,比如:
De Lathauwer L, De Moor B, Vandewalle J. A multilinear singular value decomposition[J]. SIAM journal on Matrix Analysis and Applications, 2000, 21(4): 1253-1278.
把得到的分量都拼一拼,
U
=
[
u
1
,
…
,
u
M
~
]
∈
R
M
×
M
~
,
V
=
[
v
1
,
…
,
v
N
~
]
∈
R
N
×
N
~
,
S
i
=
[
σ
i
1
,
…
,
σ
i
n
~
i
]
T
∈
R
n
~
i
×
n
i
,
i
=
1
,
…
,
D
.
U=[u1,…,u˜M]∈RMטM, V=[v1,…,v˜N]∈RNטN, Si=[σ1i,…,σ˜nii]T∈R˜ni×ni,i=1,…,D.
注意,这里的 S i \mathrm{S}_{i} Si 表达有个转置。
在线阶段:
给定约化后的空间维数
n
≤
min
{
M
~
,
N
~
}
n \leq \min \{\widetilde{M}, \widetilde{N}\}
n≤min{M
,N
},和在线参数
α
∈
A
\boldsymbol{\alpha} \in \mathcal{A}
α∈A。
首先,和 CP 一样,依然对参数分量进行某种插值表达,得到 core Matrix(下称中心矩阵),
C
e
(
α
)
=
C
×
2
(
S
1
e
1
(
α
)
)
×
3
(
S
2
e
2
(
α
)
)
⋯
×
D
+
1
(
S
D
e
D
(
α
)
)
∈
R
M
~
×
N
~
\mathrm{C}_{e}(\boldsymbol{\alpha})=\mathbf{C} \times_{2}\left(\mathrm{~S}_{1} \mathrm{e}^{1}(\boldsymbol{\alpha})\right) \times_{3}\left(\mathrm{~S}_{2} \mathrm{e}^{2}(\boldsymbol{\alpha})\right) \cdots \times_{D+1}\left(\mathrm{~S}_{D} \mathbf{e}^{D}(\boldsymbol{\alpha})\right) \in \mathbb{R}^{\widetilde{M} \times \widetilde{N}}
Ce(α)=C×2( S1e1(α))×3( S2e2(α))⋯×D+1( SDeD(α))∈RM
×N
那么事实上,
Φ
~
e
(
α
)
=
U
C
e
(
α
)
V
T
\widetilde{\Phi}_{e}(\boldsymbol{\alpha}) = \mathrm{UC}_{e}(\boldsymbol{\alpha}) \mathrm{V}^{T}
Φ
e(α)=UCe(α)VT
其次,对中心矩阵做 SVD 分解,
C
e
(
α
)
=
U
c
Σ
c
V
c
T
\mathrm{C}_{e}(\boldsymbol{\alpha})=\mathrm{U}_{c} \Sigma_{c} \mathrm{~V}_{c}^{T}
Ce(α)=UcΣc VcT
最后,就得到参数维度插值逼近后矩阵的正交分解(
U
\mathrm{U}
U 和
V
\mathrm{V}
V 一开始就是正交的,不用管),
Φ
~
e
(
α
)
=
(
U
U
c
)
Σ
c
(
V
V
c
)
T
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\left(\mathrm{UU}_{c}\right) \Sigma_{c}\left(\mathrm{VV}_{c}\right)^{T}
Φ
e(α)=(UUc)Σc(VVc)T
这时,令
{
β
i
(
α
)
}
i
=
1
n
⊂
R
M
~
\left\{\boldsymbol{\beta}_{i}(\boldsymbol{\alpha})\right\}_{i=1}^{n} \subset\mathbb{R}^{\widetilde{M}}
{βi(α)}i=1n⊂RM
是
U
c
\mathrm{U}_{c}
Uc 的前面
n
n
n 列就可以了。
Tensor Train 分解是一种特殊的 tree-based 的分解。
离线阶段:
给定快照张量和目标精度。
使用论文
Oseledets I V. Tensor-train decomposition[J]. SIAM Journal on Scientific Computing, 2011, 33(5): 2295-2317.
提到的算法,对于给定的目标精度,做张量快照的 TT 分解,
Φ
≈
Φ
~
=
∑
j
1
=
1
r
~
1
⋯
∑
j
D
+
1
=
1
r
~
D
+
1
u
j
1
∘
σ
1
j
1
,
j
2
∘
⋯
∘
σ
D
j
D
,
j
D
+
1
∘
v
j
D
+
1
\boldsymbol{\Phi} \approx \widetilde{\boldsymbol{\Phi}}=\sum_{j_{1}=1}^{\tilde{r}_{1}} \cdots \sum_{j_{D+1}=1}^{\widetilde{r}_{D+1}} \mathbf{u}^{j_{1}} \circ \boldsymbol{\sigma}_{1}^{j_{1}, j_{2}} \circ \cdots \circ \boldsymbol{\sigma}_{D}^{j_{D}, j_{D+1}} \circ \mathbf{v}^{j_{D+1}}
Φ≈Φ
=j1=1∑r~1⋯jD+1=1∑r
D+1uj1∘σ1j1,j2∘⋯∘σDjD,jD+1∘vjD+1
依然对得到的量做一些拼接,
U
=
[
u
1
,
…
,
u
r
~
1
]
∈
R
M
×
r
~
1
,
V
=
[
v
1
,
…
,
v
r
~
D
+
1
]
∈
R
N
×
r
~
D
+
1
\mathrm{U}=\left[\mathbf{u}^{1}, \ldots, \mathbf{u}^{\widetilde{r}_{1}}\right] \in \mathbb{R}^{M \times \widetilde{r}_{1}}, \quad \mathrm{~V}=\left[\mathbf{v}^{1}, \ldots, \mathbf{v}^{\widetilde{r}_{D+1}}\right] \in \mathbb{R}^{N \times \widetilde{r}_{D+1}}
U=[u1,…,ur
1]∈RM×r
1, V=[v1,…,vr
D+1]∈RN×r
D+1
(
S
i
)
j
k
q
=
(
σ
i
j
q
)
k
,
j
=
1
,
…
,
r
~
i
,
k
=
1
,
…
,
n
i
,
q
=
1
,
…
,
r
~
i
+
1
\left(\mathbf{S}_{i}\right)_{j k q}=\left(\boldsymbol{\sigma}_{i}^{j q}\right)_{k}, \quad j=1, \ldots, \widetilde{r}_{i}, \quad k=1, \ldots, n_{i}, \quad q=1, \ldots, \widetilde{r}_{i+1}
(Si)jkq=(σijq)k,j=1,…,r
i,k=1,…,ni,q=1,…,r
i+1
注意这里的每一个参数分量,都可以拼接成一个三维的 tensor。
另外,对
V
\mathrm{V}
V 取模放对角线上,构造一个量,
W
c
=
diag
(
∥
v
1
∥
,
…
,
∥
v
r
~
D
+
1
∥
)
∈
R
r
~
D
+
1
×
r
~
D
+
1
\mathrm{W}_{c}=\operatorname{diag}\left(\left\|\mathbf{v}^{1}\right\|, \ldots,\left\|\mathbf{v}^{\widetilde{r}_{D+1}}\right\|\right) \in \mathbb{R}^{\widetilde{r}_{D+1} \times \widetilde{r}_{D+1}}
Wc=diag(∥
∥v1∥
∥,…,∥
∥vr
D+1∥
∥)∈Rr
D+1×r
D+1
在线阶段:
给定约化空间的维数和在线的一组参数。
第一步,对参数张量进行插值,并做连乘得到得到中心矩阵:
C
e
(
α
)
=
∏
i
=
1
D
(
S
i
×
2
e
i
(
α
)
)
\mathrm{C}_{e}(\boldsymbol{\alpha})=\prod_{i=1}^{D}\left(\mathbf{S}_{i} \times_{2} \mathbf{e}^{i}(\boldsymbol{\alpha})\right)
Ce(α)=i=1∏D(Si×2ei(α))
则
Φ
~
e
(
α
)
=
U
C
e
(
α
)
V
T
\widetilde{\Phi}_{e}(\boldsymbol{\alpha}) = \mathrm{UC}_{e}(\boldsymbol{\alpha}) \mathrm{V}^{T}
Φ
e(α)=UCe(α)VT
第二步,对缩放过的中心矩阵做奇异值分解:
C
e
(
α
)
W
c
=
U
c
Σ
c
V
c
T
\mathrm{C}_{e}(\boldsymbol{\alpha}) \mathrm{W}_{c}=\mathrm{U}_{c} \Sigma_{c} \mathrm{~V}_{c}^{T}
Ce(α)Wc=UcΣc VcT
第三步,得到插值后矩阵的奇异值分解:
Φ
~
e
(
α
)
=
U
C
e
(
α
)
W
c
W
c
−
1
V
T
=
(
U
U
c
)
Σ
c
(
V
W
c
−
1
V
c
)
T
.
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\mathrm{UC}_{e}(\boldsymbol{\alpha}) \mathrm{W}_{c} \mathrm{~W}_{c}^{-1} \mathrm{~V}^{T}=\left(\mathrm{UU}_{c}\right) \Sigma_{c}\left(\mathrm{VW}_{c}^{-1} \mathrm{~V}_{c}\right)^{T} .
Φ
e(α)=UCe(α)Wc Wc−1 VT=(UUc)Σc(VWc−1 Vc)T.
第四步,POD 即便取决于
{
β
i
(
α
)
}
i
=
1
n
⊂
R
r
~
1
\left\{\boldsymbol{\beta}_{i}(\boldsymbol{\alpha})\right\}_{i=1}^{n} \subset \mathbb{R}^{\tilde{r}_{1}}
{βi(α)}i=1n⊂Rr~1,它为
U
c
\mathrm{U}_{c}
Uc 的前
n
n
n 列。
上述的方法都是在给定了网格 based 的参数空间快照,并不十分经济,并且当参数区域不是盒子的时候,会有麻烦。这种情况下,我们可以在线参数表示为已有参数的线性组合,
α
=
∑
j
=
1
K
a
j
α
^
j
,
e
(
α
)
=
(
a
1
,
…
,
a
K
)
T
\boldsymbol{\alpha}=\sum_{j=1}^{K} a_{j} \widehat{\boldsymbol{\alpha}}_{j}, \quad \mathbf{e}(\boldsymbol{\alpha})=\left(a_{1}, \ldots, a_{K}\right)^{T}
α=j=1∑Kajα
j,e(α)=(a1,…,aK)T
e
\mathbf{e}
e 是提取系数的,
e
:
α
→
R
K
\mathbf{e}: \boldsymbol{\alpha} \rightarrow \mathbb{R}^{K}
e:α→RK
我们假定对光滑的函数满足某种线性性,
g
(
α
)
≈
∑
j
=
1
K
a
j
g
(
α
^
j
)
g(\boldsymbol{\alpha}) \approx \sum_{j=1}^{K} a_{j} g\left(\widehat{\boldsymbol{\alpha}}_{j}\right)
g(α)≈j=1∑Kajg(α
j)
那么,对于每个离线参数,我们可以构造三阶张量,
(
Φ
)
i
j
k
=
u
i
(
t
k
,
α
^
j
)
,
i
=
1
,
…
,
M
,
j
=
1
,
…
,
K
,
k
=
1
,
…
,
N
(\boldsymbol{\Phi})_{i j k}=u_{i}\left(t_{k}, \widehat{\boldsymbol{\alpha}}_{j}\right), \quad i=1, \ldots, M, \quad j=1, \ldots, K, \quad k=1, \ldots, N
(Φ)ijk=ui(tk,α
j),i=1,…,M,j=1,…,K,k=1,…,N
那么,新来一个参数,它可以表示为已有参数的线性组合。对应快照刚好是已有快照的线性组合:
Φ
~
e
(
α
)
=
Φ
~
×
2
e
(
α
)
\widetilde{\Phi}_{e}(\boldsymbol{\alpha})=\widetilde{\boldsymbol{\Phi}} \times_{2} \mathbf{e}(\boldsymbol{\alpha})
Φ
e(α)=Φ
×2e(α)
这合理吗?
简单地说, e ( α ) \mathbf{e}(\boldsymbol{\alpha}) e(α) 从原来的插值函数变成了线性组合函数,别的 CP、HOSVD、TT ,是一样地做。
1、Model reduction and approximation: theory and algorithms[M]. Society for Industrial and Applied Mathematics, 2017.
2、Mamonov A V, Olshanskii M A. Interpolatory tensorial reduced order models for parametric dynamical systems[J]. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115122.