稀疏促进动态模态分解（SPDMD）结合DMD与稀疏优化技术，深入探讨其在多领域的应用。
原文标题：稀疏促进动态模态分解（SPDMD）详细介绍以及应用
原文作者：数据派THU
冷月清谈：
怜星夜思：
2、如何在金融领域利用SPDMD进行风险分析？
3、SPDMD如何应用于图像处理领域？
原文内容
来源：Deephub Imba本文约5800字，建议阅读5分钟本文介绍了稀疏促进动态模态分解。
SPDMD的应用
SPDMD的优势
SPDMD的核心思想
SPDMD的数学基础
-
X 是包含快照数据的矩阵
-
Φ 是DMD模态矩阵
-
b 是模态幅度向量
-
λ 是控制重建精度和稀疏性平衡的稀疏参数
实例演示：简单的模型示例
import matplotlib.colors import matplotlib.pyplot as plt import numpy as np import pandas as pd import scipy as sp from matplotlib.colors import ListedColormap import os import io %matplotlib inline plt.rcParams.update({'font.size' : 18, 'font.family' : 'Times New Roman', "text.usetex": True})
savePath = ‘.’
x = np.linspace(-10, 10, 100) t = np.linspace(0, 20, 200) Xm, Tm = np.meshgrid(x, t)M1 = 5.0 / np.cosh(0.5Xm) * np.tanh(0.5Xm) * np.exp(1.5jTm) # 主要模态
M2 = 0.5 * np.sin(2.0Xm) * np.exp(2.0jTm) # 次要模态
M3 = 0.5 * np.sin(3.2Xm) * np.cosh(3.2j*Tm) # 三级模态
M4 = 0.5 * np.random.randn(Xm.shape[0], Tm.shape[1]) # 噪声
data_matrix = M1 + M2 + M3 + M4
d = 0.1 Ub = 0.015fig, axs = plt.subplots(2, 2, figsize=(15, 15))
ax = axs[0,0]
p = ax.contourf(t, x, M1.T, levels = 501, cmap = ‘RdBu’)
ax.xaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.yaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.xaxis.set_ticks_position(‘both’)
ax.yaxis.set_ticks_position(‘both’)ax.set_ylabel(r’$x$‘)
ax.set_xlabel(r’$t$')
ax.set_title(‘Mode 1’)ax = axs[0,1]
p = ax.contourf(t, x, M2.T, levels = 501, cmap = ‘RdBu’)
ax.xaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.yaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.xaxis.set_ticks_position(‘both’)
ax.yaxis.set_ticks_position(‘both’)ax.set_ylabel(r’$x$‘)
ax.set_xlabel(r’$t$')
ax.set_title(‘Mode 2’)ax = axs[1,0]
p = ax.contourf(t, x, M3.T, levels = 501, cmap = ‘RdBu’)
ax.xaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.yaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.xaxis.set_ticks_position(‘both’)
ax.yaxis.set_ticks_position(‘both’)ax.set_ylabel(r’$x$‘)
ax.set_xlabel(r’$t$')
ax.set_title(‘Mode 3’)ax = axs[1,1]
p = ax.contourf(t, x, M4.T, levels = 501, cmap = ‘RdBu’)
ax.xaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.yaxis.set_tick_params(direction=‘in’, which=‘both’)
ax.xaxis.set_ticks_position(‘both’)
ax.yaxis.set_ticks_position(‘both’)ax.set_ylabel(r’$x$‘)
ax.set_xlabel(r’$t$')
ax.set_title(‘Noise’)plt.show()
DMD算法实现
def DMD(X1, X2, r, dt): U, s, Vh = np.linalg.svd(X1, full_matrices=False) # 截断SVD矩阵 Ur = U[:, :r] Sr = np.diag(s[:r]) Vr = Vh.conj().T[:, :r]构建Atilde并计算其特征值和特征向量
Atilde = Ur.conj().T @ X2 @ Vr @ np.linalg.inv(Sr)
Lambda, W = np.linalg.eig(Atilde)计算DMD模态
Phi = X2 @ Vr @ np.linalg.inv(Sr) @ W
omega = np.log(Lambda)/dt # 连续时间特征值
计算幅度
alpha1 = np.linalg.lstsq(Phi, X1[:, 0], rcond=None)[0] # DMD模态幅度
b = np.linalg.lstsq(Phi, X2[:, 0], rcond=None)[0] # DMD模态幅度
return Phi, omega, Lambda, alpha1, b
DMD重构
SPDMD优化问题
ADMM算法实现
from numpy import dot, multiply, diag from numpy.linalg import inv, eig, pinv, norm, solve, cholesky from scipy.linalg import svd, svdvals from scipy.sparse import csc_matrix as sparse from scipy.sparse import vstack as spvstack from scipy.sparse import hstack as sphstack from scipy.sparse.linalg import spsolvedef admm_for_dmd(P, q, s, gamma_vec, rho=1, maxiter=1000, eps_abs=1e-6, eps_rel=1e-4):
初始化返回值对象
answer = type(‘ADMMAnswer’, (object,), {})()
输入验证
P = np.squeeze(P)
q = np.squeeze(q)[:,np.newaxis]
gamma_vec = np.squeeze(gamma_vec)
if P.ndim != 2:
raise ValueError(‘invalid P’)
if q.ndim != 2:
raise ValueError(‘invalid q’)
if gamma_vec.ndim != 1:
raise ValueError(‘invalid gamma_vec’)优化变量的数量
n = len(q)
单位矩阵
I = np.eye(n)
为依赖于gamma的输出变量分配内存
answer.gamma = gamma_vec
answer.Nz = np.zeros([len(gamma_vec),]) # 非零幅度的数量
answer.Jsp = np.zeros([len(gamma_vec),]) # Frobenius范数的平方（抛光前）
answer.Jpol = np.zeros([len(gamma_vec),]) # Frobenius范数的平方（抛光后）
answer.Ploss = np.zeros([len(gamma_vec),]) # 最优性能损失（抛光后）
answer.xsp = np.zeros([n, len(gamma_vec)], dtype=‘complex’) # 幅度向量（抛光前）
answer.xpol = np.zeros([n, len(gamma_vec)], dtype=‘complex’) # 幅度向量（抛光后）矩阵P + (rho/2)*I的Cholesky分解
Prho = P + (rho/2) * I
Plow = cholesky(Prho)
Plow_star = Plow.conj().T稀疏P（用于KKT系统）
Psparse = sparse(P)
for i, gamma in enumerate(gamma_vec):
初始条件
y = np.zeros([n, 1], dtype=‘complex’) # 拉格朗日乘子
z = np.zeros([n, 1], dtype=‘complex’) # x的副本使用ADMM求解gamma参数化问题
for step in range(maxiter):
x-最小化步骤
u = z - (1/rho) * y
xnew = solve(Plow_star, solve(Plow, q + (rho/2) * u))z-最小化步骤
a = (gamma/rho) * np.ones([n, 1])
v = xnew + (1/rho) * yv的软阈值处理
znew = multiply(multiply(np.divide(1 - a, np.abs(v)), v), (np.abs(v) > a))
原始和对偶残差
res_prim = norm(xnew - znew, 2)
res_dual = rho * norm(znew - z, 2)拉格朗日乘子更新步骤
y += rho * (xnew - znew)
停止准则
eps_prim = np.sqrt(n) * eps_abs + eps_rel * np.max([norm(xnew, 2), norm(znew, 2)])
eps_dual = np.sqrt(n) * eps_abs + eps_rel * norm(y, 2)if (res_prim < eps_prim) and (res_dual < eps_dual):
break
else:
z = znew记录输出数据
answer.xsp[:,i] = z.squeeze() # 幅度向量
answer.Nz[i] = np.count_nonzero(answer.xsp[:,i]) # 非零幅度的数量
answer.Jsp[i] = (
np.real(dot(dot(z.conj().T, P), z))
- 2 * np.real(dot(q.conj().T, z))
- s) # Frobenius范数（抛光前）
非零幅度的抛光
形成约束矩阵E，使得E^T x = 0
ind_zero = np.flatnonzero(np.abs(z) < 1e-12) # 找到z的零元素的索引
m = len(ind_zero) # 零元素的数量
if m > 0: # 形成最优条件的KKT系统 E = I[:,ind_zero] E = sparse(E, dtype='complex') KKT = spvstack([ sphstack([Psparse, E], format='csc'), sphstack([E.conj().T, sparse((m, m), dtype='complex')], format='csc'), ], format='csc') rhs = np.vstack([q, np.zeros([m, 1], dtype='complex')]) # 垂直堆叠求解KKT系统
sol = spsolve(KKT, rhs)
else:
sol = solve(P, q)抛光后（最优）的幅度向量
xpol = sol[:n]
记录输出数据
answer.xpol[:,i] = xpol.squeeze()
抛光后（最优）的最小二乘残差
answer.Jpol[i] = (
np.real(dot(dot(xpol.conj().T, P), xpol))
- 2 * np.real(dot(q.conj().T, xpol))
- s)
抛光后（最优）的性能损失
answer.Ploss[i] = 100 * np.sqrt(answer.Jpol[i]/s)
print(i)
return answer
SPDMD应用
U, s, Vh = np.linalg.svd(data_matrix.T, full_matrices=False)
Vand = np.vander(Lambda, len(t), True)
P = multiply(dot(Phi.conj().T, Phi), np.conj(dot(Vand, Vand.conj().T)))
q = np.conj(diag(dot(dot(Vand, (dot(dot(U, diag(s)), Vh)).conj().T), Phi)))
s = norm(diag(s), ord=‘fro’)**2
gamma_vec = np.logspace(np.log10(0.05), np.log10(200), 150)
bNew = admm_for_dmd(P, q, s, gamma_vec, rho=1, maxiter=1000, eps_abs=1e-6, eps_rel=1e-4)
结果分析
-
Nz表示解中非零模态的数量
-
PLoss表示性能损失 - 一个归一化的度量，反映了使用最优b导致的整体模型误差
-
xpol包含对应于特定gamma值的最优b的值
