利用matlab实现POD分解(在一维信号或二维流场矢量中的应用）

（2020年7月更新，增加POD算法中svd函数的’econ‘’选项，目的是加快svd方法的计算速度。）
（2020年7月更新2，增加前言相关的正交性的简单说明）
（2020年8月更新，在程序中对旧版本矩阵与向量.*不兼容进行注释）
（2021年3月更新，更改一些细节，感谢VTgrm大佬的评论建议）

0 前言

POD是一种常用的数据降维方法，全程为Proper orthogonal decomposition，翻译为中文叫做本征正交分解。本文主要是尝试利用matlab对POD方法进行实现。
本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。
由于POD可能涉及到了包括且不限于流体、机械、信号、控制、深度学习等各领域，所以请各位大佬们如果发现错误，或者有什么补充，欢迎在评论区指出，感谢。

如果想看有关DMD分解文章，可以点击下面链接
【利用matlab实现DMD分解(在一维信号或二维流场矢量中的应用）】

文章里用到的参考列出在下面：
[1]一文让你通俗理解奇异值分解
https://www.jianshu.com/p/bcd196497d94
[2]有限长平板分离再附流动非定常特性的 PIV 实验研究（张青山）
[3]Modal Analysis of Fluid Flows: An Overview（Kunihiko Taira、Steven L. Brunton等人）

0.1 matlab中特征值计算

matlab中特征值和特征向量都可以用eig()函数进行计算。

[V,D] = eig(A)

其中V为特征向量组成的矩阵，D为特征根所组成的矩阵。A、V和D三个矩阵之间的关系为：A* V = V* D。
其中V的每一列都为一个特征向量，特征值对应的D矩阵中对角线上的相应列。

在矩阵乘法中，特征向量可以表示为变化的方向，特征根表示为变化的大小。如下图所示，矩阵A=[0.95,-0.35;-0.35,0.95]，其特征值为0.6和1.3，特征向量的方向一个为45度，一个-45°。矩阵A对于随机点阵的效果如下所示。可以看到在经过矩阵A相乘后，随机点阵发生了变形，变形的方向就是特征向量的方向，变形的大小就是特征根的大小。大于1的方向拉长，小于1的方向缩短。

0.2 matlab中SVD分解计算

SVD是一种常用的数据压缩算法，可以将一个较大的矩阵压缩为几个小矩阵相乘的形式。当然我这么说不直观，接下来先用示例解释一下。

matlab中有自带的svd函数，svd计算分解如下：

[U,S,V] = svd(A)

其中A = U * S * V '。

这时有人会质疑，哪里数据压缩了，原来的A矩阵变成了3个，明显数据量增大了呀。这时就体现出SVD的奇异值分解的作用了。SVD中的S矩阵为奇异值矩阵，也是一个对角矩阵，对角线上分布着奇异值，由大到小进行排列。奇异值越大，对矩阵的贡献越大，所以根据这个原理，我们可以删除掉一些小的奇异值，只保留大的奇异值，也可以让矩阵近似成立。

比如上图，我们只保留了前3个最大的奇异值，并将其余不参与计算的矩阵删除（用灰色表示）。如果后面的舍弃的奇异值很小，则计算出的结果与原来的A矩阵会十分接近。

接下来用一个例子来解释：

%SVD的实验
%演示SVD数据的压缩
Z = peaks(100);
[U,S,V] = svd(Z);%SVD分解

K=1;%前1个奇异值
U1=U(:,1:K);S1=S(1:K,1:K);V1=V(:,1:K);
Z1=U1*S1*V1';

K=2;%前2个奇异值
U2=U(:,1:K);S2=S(1:K,1:K);V2=V(:,1:K);
Z2=U2*S2*V2';

K=5;%前5个奇异值
U5=U(:,1:K);S5=S(1:K,1:K);V5=V(:,1:K);
Z5=U5*S5*V5';

figure(1);colormap(jet);
subplot(2,2,1)
pcolor(Z);shading interp;axis off
subplot(2,2,2)
pcolor(Z1);shading interp;axis off
subplot(2,2,3)
pcolor(Z2);shading interp;axis off
subplot(2,2,4)
pcolor(Z5);shading interp;axis off

我们看到，虽然原矩阵总共有100个奇异值，但是当我们只取到前5个奇异值时，新矩阵和原来的矩阵就之间已经看不出任何差别。之后的95个奇异值的舍弃几乎不损失任何信息，我们只保存前5个奇异值和相应的矩阵U*、V*即可。

之后也简单介绍一下SVD和特征值之间的关系。
当矩阵A可以分解为下面的形式
$$A=U\cdot S\cdot V^T$$
则矩阵V可以计算为：
$$\left(A^{T}A\right)\cdot V=D\cdot V$$
上面两式的V是相同的。
特征值和奇异值之间也存在关系，上式的特征值D开根号就等于奇异值S。
$${\sqrt{D}}_i=S_i$$

知道SVD一些基本知识后，会有助于之后POD方法的理解。

0.3 信号的正交性

POD叫做本征正交分解，分解出的信号都是正交的。
正交在几何上可以描述为，向量互相的垂直，或者可以说是一个向量在另一个向量上的投影为0。

对于信号之间的正交，我们定义当两个信号x和y的内积为0时，则正交：
$$\int x\left(t\right)y\left(t\right)dt=0$$

说完了正交，我们再提一下分解。信号的分解可以理解为，我们把信号分解为很多信号叠加的形式。比如下面公式中，把信号x分解为很多${\phi }_i$叠加的形式。
$$x\left(t\right)=\sum a_i\cdot {\phi }_i\left(t\right)$$
我们把${\phi }_i$称为基，将很多‘基’作为元素乘以常系数a，再进行叠加，我们就可以还原得到信号x。

如果改变常系数$a_i$，我们则可以得到不同的信号x，这些不同的信号x组成的一个集合，定义为空间X。通过类比几何坐标系的坐标表示方式，‘基’可以看做空间X中的坐标轴，常数$a_i$可以视为坐标值。（其实POD分解就相当于给出了一种空间X的表示方式，这个学完POD之后可以对比一下）

如果我们分解得到的基之间可以两两正交，我们就称这个分解是正交分解。最常见的正交分解就是傅里叶分解，得到的正弦函数之间互相正交。

比如以下面的信号为例：

构建了一个Perlin噪声，并进行正交分解，得到3个正交基${\phi }_1,{\phi }_2,{\phi }_3$。这三个正交基可以通过下面的方式得到原信号。
$$x=-0.954{\phi }_1+0.006{\phi }_2-0.677{\phi }_3$$

正交分解有以下的特点：
1.正交性
$$⟨{\phi }_i,{\phi }_j⟩=\{\begin{array}{l}1\left(i=j\right)\\ 0\left(i\ne j\right)\end{array}$$
自身的内积等于1，和其它正交基的内积等于0。这个性质在matlab里可以利用sum(PhiU_1 .* PhiU_2)来验证。
2 信号的能量分解前后不变（Parseval定理）
由于基是正交的，所以分解后正交基本身的能量等于1，每一个模态的能量就等于系数a的平方。所以对于POD分解，每一帧时间步i，都有

sum((U_tx(i,:)-U0x).^2) - sum(An(1,:).^2) = 0

这个可以在之后POD程序里进行验证。

3 最小平方近似性质
已知信号可以被分解为如下的形式
$$x\left(t\right)=\sum _{i=1}^n{a}_n{\phi }_n\left(t\right)$$
如果我们取前L个基（L<n），作为信号x的近似，则定义这个近似信号为y
$$y\left(t\right)=\sum _{i=1}^L{b}_n{\phi }_n\left(t\right)$$
定义最小平方近似，来描述两个信号的相似度。
$$D=sum((x-y{)}^2)$$
我们希望两个信号尽可能近似，也就是令两个信号差的平方和D最小。那么，只需要令bn=an即可。也就是说，分解得到的an不需要改动，天然的就满足于最佳近似。（POD分解中，删除掉小能量模态，得到的结果也是最佳近似的）

1 一维信号POD分解

1.1 原始的POD分解

对于随时间变化的一维信号，由于添加了时间维，所以记录在一个二维矩阵Utx内。信号的空间长度为m的离散点，时间长度为N个离散点。二维矩阵具有N行m列，每一行都代表信号在某一个时间点上的形状。
U t x = [ u ( t 1 , x 1 ) u ( t 1 , x 2 ) ⋯ u ( t 1 , x m ) u ( t 2 , x 1 ) u ( t 2 , x 2 ) ⋯ u ( t 2 , x m ) ⋮ ⋮ ⋱ ⋮ u ( t N , x 1 ) u ( t N , x 2 ) ⋯ u ( t N , x m ) ] U_{tx}=\left[ \begin{matrix} u\left( t_1,x_1 \right)& u\left( t_1,x_2 \right)& \cdots& u\left( t_1,x_m \right)\\ u\left( t_2,x_1 \right)& u\left( t_2,x_2 \right)& \cdots& u\left( t_2,x_m \right)\\ \vdots& \vdots& \ddots& \vdots\\ u\left( t_N,x_1 \right)& u\left( t_N,x_2 \right)& \cdots& u\left( t_N,x_m \right)\\ \end{matrix} \right] Utx​= ​u(t1​,x1​)u(t2​,x1​)⋮u(tN​,x1​)​u(t1​,x2​)u(t