【机器学习】基于稀疏识别方法的洛伦兹混沌系统预测

1. 引言

1.1. DNN模型的来由

从数据中识别非线性动态学意味着什么?
假设我们有时间序列数据,这些数据来自一个(非线性)动态学系统。

  • 识别一个系统意味着基于数据推断该系统的控制方程。换句话说,就是找到动态系统方程 x ˙ = f ( x ) \mathbf{\dot x} = f(\mathbf{x}) x˙=f(x)中的 f f f(其中 x \mathbf{x} x 可能是向量值)。
    在这里插入图片描述
    例如,对于Lorenz系统,我们希望从时间序列数据中学习到方程右边的部分。

1.2. 研究稀疏性的意义

为什么我们需要稀疏性?
在这里,稀疏性意味着控制方程中的项数很少。稀疏性是有益的,因为它更具:

  • 可解释性。在需要理解变量及其相互作用的应用中至关重要,例如在需要安全关键保证的应用中。
  • 泛化能力。如果正确,方程将准确描述训练数据所填充的状态空间区域之外的动态。
    通常,人们可以将SINDy识别的模型视为与物理方程相对的模型,而不是大型、不透明的深度神经网络。

2. SINDy算法

SINDy试图找到适合数据 X ˙ = f ( X ) \mathrm{\dot X} = f(\mathrm{X}) X˙=f(X)的动态系统 f f f。这个函数逼近问题被表述为线性回归 X ˙ = Θ ( X ) Ξ \mathrm{\dot X} = \Theta(\mathrm{X}) \Xi X˙=Θ(X)Ξ,其中系数为 Ξ \Xi Ξ 和回归项库 Θ ( X ) \Theta(X) Θ(X)。算法分为三个步骤:

  • 从动态系统生成数据 X X X 并计算导数 X ˙ \dot X X˙
  • 建立候选项库 Θ ( X ) \Theta(X) Θ(X) 作为 X X X上的函数。
  • 稀疏回归系数 Ξ \Xi Ξ,以最好地描述数据。
  1. SINDy假设是测量了 n n n 维数据点的时间序列 x = ( x 1 , … x n ) \mathbf{x}=(x_1, \ldots x_n) x=(x1,xn) m m m 个时间步 t 1 , … , t m t_1, \ldots, t_m t1,,tm,我们定义数据矩阵 X X X 和导数矩阵 X ˙ \dot X X˙ 为:
    X = [ x 1 ( t 1 ) x 2 ( t 1 ) ⋯ x n ( t 1 ) x 1 ( t 2 ) x 2 ( t 2 ) ⋯ x n ( t 2 ) x 1 ( t 3 ) x 2 ( t 3 ) ⋯ x n ( t 3 ) ⋮ ⋮ ⋱ ⋮ x 1 ( t m ) x 2 ( t m ) ⋯ x n ( t m ) ] X=\begin{bmatrix} x_1(t_1)&x_2(t_1)&\cdots&x_n(t_1)\\ x_1(t_2)&x_2(t_2)&\cdots&x_n(t_2)\\ x_1(t_3)&x_2(t_3)&\cdots&x_n(t_3)\\ \vdots&\vdots&\ddots&\vdots\\ x_1(t_m)&x_2(t_m)&\cdots&x_n(t_m)\\ \end{bmatrix} X= x1(t1)x1(t2)x1(t3)x1(tm)x2(t1)x2(t2)x2(t3)x2(tm)xn(t1)xn(t2)xn(t3)xn(tm)
    X ˙ = [ x ˙ 1 ( t 1 ) x ˙ 2 ( t 1 ) ⋯ x ˙ n ( t 1 ) x ˙ 1 ( t 2 ) x ˙ 2 ( t 2 ) ⋯ x ˙ n ( t 2 ) x ˙ 1 ( t 3 ) x ˙ 2 ( t 3 ) ⋯ x ˙ n ( t 3 ) ⋮ ⋮ ⋱ ⋮ x ˙ 1 ( t m ) x ˙ 2 ( t m ) ⋯ x ˙ n ( t m ) ] \dot{X}=\begin{bmatrix} \dot{x}_1(t_1)&\dot{x}_2(t_1)&\cdots&\dot{x}_n(t_1)\\ \dot{x}_1(t_2)&\dot{x}_2(t_2)&\cdots&\dot{x}_n(t_2)\\ \dot{x}_1(t_3)&\dot{x}_2(t_3)&\cdots&\dot{x}_n(t_3)\\ \vdots&\vdots&\ddots&\vdots\\ \dot{x}_1(t_m)&\dot{x}_2(t_m)&\cdots&\dot{x}_n(t_m)\\ \end{bmatrix} X˙= x˙1(t1)x˙1(t2)x˙1(t3)x˙1(tm)x˙2(t1)x˙2(t2)x˙2(t3)x˙2(tm)x˙n(t1)x˙n(t2)x˙n(t3)x˙n(tm)
  2. 接下来,我们定义库矩阵 Θ ( X ) \Theta(X) Θ(X),其列是应用于数据的一组基函数 { θ l } l = 1 , … , L \{\theta_l\}_{l=1,\ldots, L} {θl}l=1,,L
    Θ ( X ) = [ ∣ ∣ ∣ Θ 1 ( X ) Θ 2 ( X ) ⋯ Θ n ( X ) ∣ ∣ ∣ ] \Theta(X)=\begin{bmatrix} |&|&&|\\ \Theta_1(X)&\Theta_2(X)&\cdots&\Theta_n(X)\\ |&|&&|\\ \end{bmatrix} Θ(X)= Θ1(X)Θ2(X)Θn(X)
    以下是对上述文本的改写:
    一些简单的例子包括多项式的基础,比如在泰勒级数展开中使用的多项式,或者三角函数,如 sin(x_1), cos(x_1), sin(2x_1), 等等,这些在傅里叶级数展开中常见。然而,根据不同的问题,可能需要使用更复杂的基础函数,例如贝塞尔函数。
  • 最后,我们使用稀疏线性回归算法(例如LASSO)找到系数 Ξ \Xi Ξ,使得:
    Ξ = [ ∣ ∣ ∣ ε 1 ( X ) ε 2 ( X ) ⋯ ε n ( X ) ∣ ∣ ∣ ] \Xi=\begin{bmatrix} |&|&&|\\ \varepsilon_1(X)&\varepsilon_2(X)&\cdots&\varepsilon_n(X)\\ |&|&&|\\ \end{bmatrix} Ξ= ε1(X)ε2(X)εn(X)
    因此
    X ˙ = θ ( X ) Ξ \dot{X}=\theta(X)\Xi X˙=θ(X)Ξ
    在这里插入图片描述

3.线性动态系统示例

假设我们测量了一个由以下动力学系统控制的粒子轨迹:

d d t ( x y ) = ( − 2 x y ) = ( − 2 0 0 1 ) ( x y ) \frac{d}{dt} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} -2x \\ y \end{pmatrix} = \begin{pmatrix} -2 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} dtd(xy)=(2xy)=(2001)(xy)

3.1.构建数据矩阵

以初始条件 x 0 = 3 x_0=3 x0=3 y 0 = 1 3 y_0=\frac{1}{3} y0=31我们构建数据矩阵 X X X

import numpy as np
import pysindy as pst = np.linspace(0, 1, 100)
x = 3 * np.exp(-2 * t)
y = 0.5 * np.exp(t)
X = np.stack((x, y), axis=-1)

3.2.拟合模型

在构建SINDy模型时,我们有几个关键选项需要确定,包括微分方法、基函数库以及优化器的选择。现在,让我们具体探讨一下以傅里叶基作为基函数的情况。

model = ps.SINDy(differentiation_method=ps.FiniteDifference(order=2),feature_library=ps.FourierLibrary(),optimizer=ps.STLSQ(threshold=0.2),feature_names=["x", "y"]
)
model.fit(X, t=t)

3.3.测试模型

拟合模型后,可以使用 print 成员函数检查模型。

model.print()
(x)' = 0.772 sin(1 x) + 2.097 cos(1 x) + -2.298 sin(1 y) + -3.115 cos(1 y)
(y)' = 1.362 sin(1 y) + -0.222 cos(1 y)

然后使用测试数据验证

import matplotlib.pyplot as pltdef plot_simulation(model, x0, y0):t_test = np.linspace(0, 1, 100)x_test = x0 * np.exp(-2 * t_test)y_test = y0 * np.exp(t_test)sim = model.simulate([x0, y0], t=t_test)plt.figure(figsize=(6, 4))plt.plot(x_test, y_test, label="Ground truth", linewidth=4)plt.plot(sim[:, 0], sim[:, 1], "--", label="SINDy estimate", linewidth=3)plt.plot(x0, y0, "ko", label="Initial condition", markersize=8)plt.xlabel("x")plt.ylabel("y")plt.legend()plt.show()x0 = 6
y0 = -0.1
plot_simulation(model, x0, y0)
x0 = 6
y0 = -0.1
plot_simulation(model, x0,y0)

在这里插入图片描述
正如我们预期的那样,对于这个问题,傅里叶基并不是最佳选择。现在,让我们尝试一个不同的基函数!

3.4.自定义基函数

实验1:研究当我们选择一个不同的基函数时,会有什么样的结果。

为了实现这个目标,我们将编写一个Sindy算法类,并使用一个自定义的基函数库,而不是默认的ps.PolynomialLibrary()

我们将这样初始化ps.SINDy():将feature_library属性设置为我们的自定义傅里叶库(如果可用)。然后,我们将使用.fit(X, t=t)方法来拟合这个实例化的算法。

# 导入所需的库(假设 ps 是 PySINDy 的别名)  
# 注意:PySINDy 库的实际导入可能需要根据你的环境进行调整  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
from pysindy.utils import plot_simulation  # 初始化一个 SINDy 模型,使用有限差分法(二阶)进行微分  
# 这里使用了多项式库(一阶)作为特征库  
# 使用了阈值为 0.2 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = SINDy(  differentiation_method=FiniteDifference(order=2),  # 使用二阶有限差分法进行微分  feature_library=PolynomialLibrary(degree=1),       # 使用一阶多项式库作为特征库  optimizer=STLSQ(threshold=0.2),                    # 使用阈值为 0.2 的 STLSQ 优化器  feature_names=["x", "y"]                           # 特征名称设置为 'x' 和 'y'  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # 打印模型结果  
model_1.print()  # 假设 plot_simulation 是一个自定义函数,用于绘制模拟结果  
# 这里设置初始条件 x0=6, y0=-0.1  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述
在此简单示例中,模型利用多项式基函数准确地识别了背后的数学方程。

我们面对的是一个一阶多项式微分方程,因此,使用与这些模型假设完全一致的回归方法来识别它是合情合理的。

实验1.1:尝试逐步提高多项式的度数,观察在哪个度数时SINDy无法正确识别方程,以及在哪个度数时测试数据集上的预测开始出现偏差。这将帮助我们了解SINDy在不同复杂度下的效能表现。

# 导入 PySINDy 相关的库(假设 ps 是 PySINDy 的别名)  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
from pysindy.utils import plot_simulation  # 假设 plot_simulation 函数已经定义或可用  # 初始化一个 SINDy 模型  
# 使用二阶有限差分法进行微分  
# 使用四阶多项式库作为特征库  
# 使用阈值为 0.2 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = ps.SINDy(  differentiation_method=ps.FiniteDifference(order=2),  # 使用二阶有限差分法进行微分  feature_library=ps.PolynomialLibrary(degree=4),       # 使用四阶多项式库作为特征库  optimizer=ps.STLSQ(threshold=0.2),                    # 使用阈值为 0.2 的 STLSQ 优化器  feature_names=["x", "y"]                              # 特征名称设置为 'x' 和 'y'  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # (注意:以下部分是错误的,因为不应该再次初始化 SINDy,并且没有缩进)  
# 正确的方式是上面已经初始化和拟合了 model_1,接下来直接使用它  
# SINDy(differentiation_method=FiniteDifference(),  
#       feature_library=PolynomialLibrary(degree=4), feature_names=['x', 'y'],  
#       optimizer=STLSQ(threshold=0.2))  # 打印模型结果  
model_1.print()  # 设置初始条件 x0=6, y0=-0.1  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述
实验1.2:探讨阈值设定对模型性能的影响。如果阈值设定过低,将导致哪些问题?相反,如果阈值设定过高,又将引发哪些后果?

STLSQ(序列阈值最小二乘法):这种方法通过迭代进行最小二乘优化,并对于权重向量中小于特定阈值的元素进行置零处理,以此来最小化目标函数 ∣ ∣ y − X w ∣ ∣ 2 2 + α ∣ ∣ w ∣ ∣ 2 2 ||y-X_w||^2_2+\alpha||w||^2_2 ∣∣yXw22+α∣∣w22

阈值:这是权重向量中系数的最小允许幅度。任何幅度小于这个阈值的系数都将被清零,从而有助于实现模型的稀疏性。

# 导入 PySINDy 库(假设 ps 是 PySINDy 的别名)  
from pysindy import SINDy, FiniteDifference, PolynomialLibrary, STLSQ  
# 假设 plot_simulation 函数已经定义或可用  
from pysindy.utils import plot_simulation  # 初始化 SINDy 模型  
# 使用二阶有限差分法进行微分  
# 使用四阶多项式库作为特征库  
# 使用阈值为 0.1 的 STLSQ 优化器  
# 特征名称设置为 'x' 和 'y'  
model_1 = ps.SINDy(  differentiation_method=ps.FiniteDifference(order=2),  # 微分方法:二阶有限差分法  feature_library=ps.PolynomialLibrary(degree=4),       # 特征库:四阶多项式库  optimizer=ps.STLSQ(threshold=0.1),                    # 优化器:STLSQ,阈值设置为 0.1  feature_names=["x", "y"]                              # 特征名称:['x', 'y']  
)  # 使用数据 X 和时间 t 拟合模型  
model_1.fit(X, t=t)  # 注意:以下部分是错误的,它试图再次初始化一个 SINDy 模型,但并未赋值给任何变量  
# 这里我们忽略它,因为 model_1 已经正确初始化和拟合  
# SINDy(differentiation_method=FiniteDifference(),  
#       feature_library=PolynomialLibrary(degree=4), feature_names=['x', 'y'],  
#       optimizer=STLSQ())  # 打印模型结果  
model_1.print()  # 设置初始条件  
x0 = 6  
y0 = -0.1  # 调用 plot_simulation 函数绘制模拟结果  
plot_simulation(model_1, x0, y0)

在这里插入图片描述

4. 洛伦兹吸引子实验

当使用 SINDy 对 洛伦兹吸引子(Lorenz Attractor)进行测试时,我们首先通过其真实的动力学方程来模拟一条轨迹。这样,我们可以使用模拟的数据来验证 SINDy 方法是否能够有效地识别出 Lorenz 系统的动态结构。

import numpy as np  
import matplotlib.pyplot as plt  
from scipy.integrate import odeint  
from mpl_toolkits.mplot3d import Axes3D  
# 假设 ps 是 PySINDy 的别名  
from pysindy import SINDy, STLSQ, PolynomialLibrary  # Lorenz 吸引子的参数  
rho = 28.0  
sigma = 10.0  
beta = 8.0 / 3.0  
dt = 0.01  # Lorenz 系统的动力学方程  
def f(state, t):  x, y, z = state  return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]  # 初始状态  
state0 = [1.0, 1.0, 1.0]  
# 时间步长数组  
time_steps = np.arange(0.0, 40.0, dt)  # 模拟 Lorenz 吸引子的轨迹  
x_train = odeint(f, state0, time_steps)  # 初始化 SINDy 模型  
model = SINDy(  optimizer=STLSQ(threshold=0.05),  # 使用阈值为 0.05 的 STLSQ 优化器  feature_library=PolynomialLibrary(degree=2),  # 使用二阶多项式库  
)  # 使用模拟的数据拟合 SINDy 模型  
model.fit(x_train, t=dt)  # 使用 SINDy 模型进行模拟  
x_sim = model.simulate(x_train[0], time_steps)  # 打印 SINDy 模型的识别结果  
model.print()  # 绘制 Lorenz 吸引子的轨迹  
plt.figure(figsize=(6, 4))  
# 绘制真实轨迹  
plt.plot(x_train[:, 0], x_train[:, 2], label='真实轨迹')  
# 绘制 SINDy 估计的轨迹  
plt.plot(x_sim[:, 0], x_sim[:, 2], '--', label='SINDy 估计轨迹')  
# 绘制初始条件  
plt.plot(x_train[0, 0], x_train[0, 2], "ko", label="初始条件", markersize=8)  
plt.legend()  
plt.xlabel('x')  
plt.ylabel('z')  
plt.title('Lorenz 吸引子轨迹')  
plt.draw()  
plt.show()

在这里插入图片描述

import matplotlib.pyplot as plt# 定义一个函数,用于绘制特定维度的时间序列数据
def plot_dimension(dim, name):# 创建一个图形对象,设置图形大小为9x2英寸fig = plt.figure(figsize=(9, 2))# 获取当前的坐标轴对象ax = fig.gca()# 在坐标轴上绘制训练数据的时间序列ax.plot(time_steps, x_train[:, dim])# 在相同的坐标轴上绘制模拟数据的时间序列,使用虚线表示ax.plot(time_steps, x_sim[:, dim], "--")# 设置x轴标签为"时间"plt.xlabel("时间")# 设置y轴标签为维度名称plt.ylabel(name)# 调用函数绘制x、y、z三个维度的时间序列数据
plot_dimension(0, 'x')  # 绘制x维度
plot_dimension(1, 'y')  # 绘制y维度
plot_dimension(2, 'z')  # 绘制z维度

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

5. 挑战

将SINDy应用于真实世界数据更具挑战性,主要体现在以下几个方面:

  1. 数据:需要多少数据以及数据的质量如何?SINDy默认依赖于干净、快速采样且噪声小的数据。数据需要足够干净,以便计算导数。但这些要求可以相互权衡:例如,如果数据干净,少量即可。反之,对于大量但嘈杂的数据,可以通过积分来平均噪声。
  2. 坐标:应该测量哪些变量?如果可能,使用专家领域知识。或者应用奇异值分解来找到最相关的变量。或者将深度自编码器与SINDy结合训练。
  3. :哪些库项最能捕捉动态?实际上,从简单的开始,例如线性项,然后是二次项等。也可以使用领域知识。注意,不要过度扩展库,以避免列的线性依赖导致 (\Theta(X)) 矩阵病态。
  4. 优化:使用哪种算法?如何实际找到稀疏项?最简单的是LASSO,实践中,顺序阈值最小二乘(迭代硬阈值系数)效果更好。在这里,也可以纳入物理信息约束:例如,可以直接在优化中实施能量守恒或全局稳定性。

在过去的几年中,这方面取得了很多进展,SINDy已成功应用于各种真实世界问题,包括发现由于问题复杂性而无法用解析方法表示的等离子体流体动力学方程。

附录 A:水库计算实例

水库计算是一种简单的方法,用于在没有反向传播通过时间的情况下训练递归神经网络,以及与之相关的臭名昭著的梯度消失和爆炸问题。基本步骤如下:

  1. 随机初始化递归神经网络权重。
  2. 固定隐藏连接权重。
  3. 使用线性回归训练线性输出层。

更正式地说,递归神经网络的状态由其隐藏激活 h ∈ R M \mathbf{h} \in \mathbb{R}^{M} hRM给出,这些激活通过随机初始化并固定的矩阵 W h \mathbf{W}_{\mathrm{h}} Wh 连接。输入序列 X 1 : τ \mathbf{X}_{1:\tau} X1:τ 通过线性映射 W i \mathbf{W}_{\mathrm{i}} Wi嵌入到状态 h τ \mathbf{h}_{\tau} hτ
h 1 = θ ( W h h 0 + W i X 1 ) ⋮ h τ = θ ( W h h τ − 1 + W i X τ ) h_1=\theta(W_hh_0+W_iX_1)\\ \vdots\\ h_{\tau}=\theta(W_hh_{\tau-1}+W_iX_{\tau}) h1=θ(Whh0+WiX1)hτ=θ(Whhτ1+WiXτ)

然后,水库可以基于隐藏状态 h τ \mathbf{h}_{\tau} hτ线性预测下一个输入 x τ + 1 \mathbf{x}_{\tau +1} xτ+1
x ^ τ + 1 = W o h τ \mathbf{\hat {x}}_{\tau +1} = \mathbf{W}_{\mathrm{o}} \mathbf{h}_\tau x^τ+1=Wohτ
其中 W o \mathbf{W}_{\mathrm{o}} Wo 将隐藏激活映射到输出。只有这个输出矩阵 W o \mathbf{W}_{\mathrm{o}} Wo的参数被训练。因此,可以通过岭回归(正则化参数 λ \lambda λ)以解析方式获得最优值:
W o = ( H T H + λ I ) − 1 H T X W_{\mathrm{o}}=(H^TH+\lambda\mathrm{I})^{-1}H^TX Wo=(HTH+λI)1HTX
其中 H \mathbf{H} H 是隐藏状态, I \mathbf{I} I是单位矩阵, X \mathbf{X} X是回归的目标。

A.1. 设置

# 导入scipy库中的sparse模块  
from scipy import sparse  # 定义一些参数  
# 半径  
radius = 0.6  
# 稀疏度  
sparsity = 0.01  
# 输入的维度  
input_dim = 3  
# 水库(或称为储备池)的大小  
reservoir_size = 1000  
# 预热步数  
n_steps_prerun = 10  
# 正则化参数  
regularization = 1e-2   

A.2.随机初始化网络权重

代码是用于初始化水库计算中的网络权重的Python代码。

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import eigs# 定义水库大小和输入维度
reservoir_size = 1000  # 假设的水库大小
input_dim = 3         # 输入数据的维度
sparsity = 0.01      # 稀疏性参数
radius = 0.6          # 权重的半径# 随机初始化隐藏层权重矩阵,大小为reservoir_size x reservoir_size
weights_hidden = sparse.random(reservoir_size, reservoir_size, density=sparsity)# 计算隐藏层权重矩阵的特征值
eigenvalues, _ = eigs(weights_hidden)# 标准化权重以确保最大的特征值的绝对值不超过radius
weights_hidden = weights_hidden / np.max(np.abs(eigenvalues)) * radius# 初始化输入到隐藏层的权重矩阵,大小为reservoir_size x input_dim
weights_input = np.zeros((reservoir_size, input_dim))# 为每个输入维度分配一个子空间
q = int(reservoir_size / input_dim)
for i in range(input_dim):# 为每个输入维度随机生成-1到1之间的值weights_input[i * q:(i + 1) * q, i] = 2 * np.random.rand(q) - 1# 初始化输出层权重矩阵,大小为input_dim x reservoir_size
weights_output = np.zeros((input_dim, reservoir_size))

代码功能说明:

  • 首先,使用scipy.sparse.random函数随机生成一个稀疏的weights_hidden矩阵,该矩阵表示隐藏层神经元之间的连接权重。
  • 然后,使用scipy.sparse.linalg.eigs函数计算weights_hidden矩阵的特征值,以确保权重矩阵的条件数在合理范围内。
  • 接着,根据最大的绝对特征值调整权重矩阵的规模,以满足指定的radius参数。
  • 为输入层到隐藏层的连接初始化一个全零矩阵weights_input,然后为每个输入维度随机分配权重。
  • 最后,初始化输出层的权重矩阵weights_output为一个全零矩阵,这个矩阵将在后续的训练过程中被优化。

这些权重矩阵将用于水库计算模型中,该模型通过固定隐藏层权重并仅训练输出层来简化训练过程。

B.3 嵌入序列

将序列嵌入到网络的隐藏状态中。
代码定义了几个函数,用于初始化和处理水库计算中的隐藏状态。

import numpy as np# 初始化隐藏状态函数
def initialize_hidden(reservoir_size, n_steps_prerun, sequence):# 创建一个全零的隐藏状态矩阵,大小为(reservoir_size, 1)hidden = np.zeros((reservoir_size, 1))# 预热阶段,不更新隐藏状态,只用于填充初始隐藏状态for t in range(n_steps_prerun):# 获取输入序列并重塑形状input = sequence[t].reshape(-1, 1)# 更新隐藏状态,使用tanh激活函数hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)return hidden# 增强隐藏状态的函数,例如通过平方操作
def augment_hidden(hidden):# 复制隐藏状态矩阵h_aug = hidden.copy()# 对每隔一个元素进行平方操作,增强状态表示h_aug[::2] = pow(h_aug[::2], 2.0)return h_aug# 使用序列数据和预热步数调用初始化隐藏状态函数
hidden = initialize_hidden(reservoir_size, n_steps_prerun, sequence)# 初始化存储隐藏状态和目标状态的列表
hidden_states = []
targets = []# 对序列中每个时间步进行循环
for t in range(n_steps_prerun, len(sequence) - 1):# 重塑输入和目标向量形状input = np.reshape(sequence[t], (-1, 1))target = np.reshape(sequence[t + 1], (-1, 1))# 更新隐藏状态hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)# 增强隐藏状态hidden = augment_hidden(hidden)# 将当前隐藏状态添加到列表中hidden_states.append(hidden)# 将目标状态添加到列表中targets.append(target)# 将列表转换为数组,并去除单维度条目
targets = np.squeeze(np.array(targets))
hidden_states = np.squeeze(np.array(hidden_states))

代码功能说明:

  • initialize_hidden 函数用于初始化隐藏状态,通过预热步骤填充初始状态,但不更新状态。
  • augment_hidden 函数用于增强隐藏状态,例如通过平方操作来增加状态的表达能力。
  • 在主循环中,对于序列中的每个时间步,更新隐藏状态,并通过augment_hidden函数增强状态,然后将状态和目标添加到各自的列表中。
  • 最后,将隐藏状态和目标状态的列表转换为NumPy数组,并使用np.squeeze去除单维度条目,以便用于后续的线性回归或其他处理。

B.4.获取线性输出层权重

使用岭回归来获取线性输出层的权重。
代码展示了如何使用Python实现水库计算中的输出权重计算和预测函数,以及如何绘制预测结果与实际数据的对比图。

import numpy as np
import matplotlib.pyplot as plt# 计算输出权重矩阵
weights_output = (np.linalg.inv(hidden_states.T @ hidden_states + regularization * np.eye(reservoir_size)) @hidden_states.T @ targets).T# 定义预测函数
def predict(sequence, n_steps_predict):# 初始化隐藏状态hidden = initialize_hidden(reservoir_size, n_steps_prerun, sequence)# 重塑输入向量input = sequence[n_steps_prerun].reshape((-1, 1))# 初始化输出列表outputs = []for t in range(n_steps_prerun, n_steps_prerun + n_steps_predict):# 更新隐藏状态hidden = np.tanh(weights_hidden @ hidden + weights_input @ input)# 增强隐藏状态hidden = augment_hidden(hidden)# 计算输出output = weights_output @ hidden# 更新输入为当前输出input = output# 将输出添加到列表中outputs.append(output)# 返回预测结果作为数组return np.array(outputs)# 使用predict函数进行预测
x_sim = predict(sequence, 4000)# 绘制真实数据和预测结果的对比图
plt.figure(figsize=(6, 4))
plt.plot(x_train[:4000, 0], x_train[:4000, 2], label="Ground truth")
plt.plot(x_sim[:, 0], x_sim[:, 2], '--', label="Reservoir computing estimate")
plt.plot(x_train[0, 0], x_train[0, 2], "ko", label="Initial condition", markersize=8)
plt.legend()
plt.show()

代码功能说明:

  • weights_output 的计算使用了岭回归,其中 hidden_states 是隐藏状态的集合,targets 是目标值,regularization 是正则化参数,用于防止过拟合。
  • predict 函数接收一个序列和预测步数 n_steps_predict,初始化隐藏状态,并循环进行预测。在每一步,使用当前隐藏状态和输入来计算下一个时间步的输出。
  • initialize_hidden 函数用于初始化隐藏状态,augment_hidden 函数用于增强隐藏状态,例如通过平方操作。
  • 最后,使用 matplotlib 库绘制真实数据和水库计算估计结果的对比图,以及初始条件的标记。

附录B. 洛伦兹吸引子简介

洛伦兹吸引子(Lorenz Attractor)是一个在混沌理论中非常著名的动态系统模型,由气象学家爱德华·洛伦兹(Edward Lorenz)在1963年提出。这个模型描述了一个三变量的自治常微分方程组,通常用于展示混沌现象,特别是在大气对流的研究中。

洛伦兹吸引子的数学模型由以下三个方程组成:

d x d t = σ ( y − x ) \frac{dx}{dt} = \sigma(y - x) dtdx=σ(yx)
d y d t = x ( ρ − z ) − y \frac{dy}{dt} = x(\rho - z) - y dtdy=x(ρz)y
d z d t = x y − β z \frac{dz}{dt} = xy - \beta z dtdz=xyβz

其中,(x), (y), 和 (z) 是状态变量,而 (\sigma), (\rho), 和 (\beta) 是系统的参数。这组方程描述了流体对流过程中的三个关键变量:对流强度((x)),水平温度梯度((y)),以及垂直温度梯度((z))。

当参数 (\sigma), (\rho), 和 (\beta) 被设置为某些值时(例如,(\sigma = 10), (\rho = 28), (\beta = \frac{8}{3})),系统将会展现出混沌行为。此时,洛伦兹吸引子将表现为一个双螺旋状的吸引子,系统状态((x), (y), (z))的轨迹会在这个吸引子上无限循环,但永远不会重复。

洛伦兹吸引子的特性包括:

  1. 敏感性:系统对初始条件非常敏感,即使初始条件只有微小的差异,长期下来也会导致截然不同的结果。这就是所谓的“蝴蝶效应”。
  2. 混沌性:系统具有混沌特性,即其动态行为看起来是随机的,但实际上是由确定性方程产生的。
  3. 吸引性:尽管系统表现出混沌行为,但其轨迹仍然被限制在一个有限的区域内,即洛伦兹吸引子。

洛伦兹吸引子不仅在气象学和流体力学中有着广泛的应用,还成为了研究混沌现象和复杂系统动力学的重要工具。它揭示了自然界中许多复杂系统可能存在的普遍行为,对于理解这些系统的行为和控制它们具有重要的指导意义。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/856464.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

XXL-Job实战(一)

​需求介绍:构建一个分布式短信发送系统,应对双十一活动需向1000万用户快速推送营销短信的挑战,每条数据的业务处理逻辑为0.1s。对于普通任务来说,只有一个线程来处理 可能需要10万秒才能处理完,业务则严重受影响。 常…

5款堪称变态的AI神器,焊死在电脑上永不删除!

一 、AI视频合成工具——Runway: 第一款RunWay,你只需要轻轻一抹,视频中的元素就会被擦除,再来轻轻一抹,直接擦除,不喜欢这个人直接擦除,一点痕迹都看不出来。 除了视频擦除功能外,…

【AI大模型】Transformers大模型库(十):repetition_penalty惩罚系数

目录​​​​​​​ 一、引言 二、惩罚系数repetition_penalty 2.1 概述 2.2 使用说明 2.3 使用示例 三、总结 一、引言 这里的Transformers指的是huggingface开发的大模型库,为huggingface上数以万计的预训练大模型提供预测、训练等服务。 🤗 T…

韩顺平0基础学Java——第27天

p548-568 明天开始坦克大战 Entry 昨天没搞明白的Map、Entry、EntrySet://GPT教的 Map 和 Entry 的关系 1.Map 接口:它定义了一些方法来操作键值对集合。常用的实现类有 HashMap、TreeMap 等。 2. Entry接口:Entry 是 Map 接口的一个嵌…

vivado TILE

TILE是包含一个或多个SITE对象的设备对象。可编程逻辑TILE 包括各种各样的对象,如SLICE/CLB、BRAM、DSP、I/O块、时钟资源,以及 GT块。从结构上讲,每个瓦片都有许多输入和输出,并且可编程 互连以将瓦片的输入和输出连接到任何其他…

实现一个简易动态线程池

项目完整代码:https://github.com/YYYUUU42/Yu-dynamic-thread-pool 如果该项目对你有帮助,可以在 github 上点个 ⭐ 喔 🥰🥰 1. 线程池概念 2. ThreadPoolExecutor 介绍 2.1. ThreadPoolExecutor是如何运行,如何同时…

elementUI的el-table自定义表头

<el-table-column label"昨日仪表里程(KM)" align"left" min-width"190" :render-header"(h, obj) > renderHeader(h, obj, 参数)" > <template slot-scope"scope"> <span>{{ scope.row.firstStartMil…

流程图工具评测:十大热门软件对比

流程图是一种用图形符号和箭头表示工作流程的图形表示方法。它展示了一系列相互关联的步骤&#xff0c;以显示过程中数据或物质的流动、决策点和操作步骤。流程图广泛用于各种领域&#xff0c;包括业务流程、软件开发、工程等&#xff0c;以帮助人们更好地理解和分析工作流程。…

大模型应用开发实践:RAG与Agent

RAG planning是任务拆解的一些方法。 Agent RAG现在基本上推荐LangChain开发框架。而Agent目前没有一个通用的好的开发框架/范式。 学习路径

达梦8 兼容MySQL语法支持非分组项作为查询列

MySQL 数据库迁移到达梦后&#xff0c;部分GROUP BY语句执行失败&#xff0c;报错如下&#xff1a; 问题原因&#xff1a; 对于Oracle数据库&#xff0c;使用GROUP BY时&#xff0c;SELECT中的非聚合列必须出现在GROUP BY后面&#xff0c;否则就会报上面的错误&#xff0c;达梦…

使用宝塔面板搭建Flask项目保姆级喂饭教程

目录 零.前言 一.准备工作 1.1创建requirements.txt文件 1.2将项目打包为压缩文件 1.3租一台服务器 1.4部署宝塔面板 二.宝塔面板(服务器)上的操作 2.1将本地Flask项目上传到服务器 2.2添加Python项目 2.3配置Python项目 2.4配置Nginx 2.5宝塔面板放行端口 2.6在服…

首个AI高考评测结果出炉,GPT-4o排名第二

近日&#xff0c;上海人工智能实验室利用其自主研发的“司南”评测体系OpenCompass&#xff0c;对国内外多个知名大模型进行了一场特殊的“高考”。这些来自阿里巴巴、智谱AI、Mistral等机构&#xff0c;以及OpenAI的GPT-4o等“考生”&#xff0c;接受了新课标I卷“语数外”的全…

百万级 QPS 接入层网关架构方案演进

文章目录 前言1、单机架构2、DNS 轮询3、Nginx 单机4、Nginx 主备 Keepalived5、LVS 主备 Keepalived Nginx 集群6、LVS 主备 Keepalived Nginx 集群 DNS 轮询 前言 随着PC、移动互联网的快速发展&#xff0c;越来越多的人通过手机、电脑、平板等设备访问各种各样APP、网…

找不到com.fasterxml.jackson.core.exc.StreamWriteException的类文件

1. 前言: 使用springboot搭建的项目, 需要使用 jackson 更改json文件的内容; maven管理jar包, 导入jar包版本信息如下: <!-- 读写json文件所需依赖 --> <dependency><groupId>com.fasterxml.jackson.core</groupId><artifactId>jackson-databin…

C++语法06 格式化输出及保留小数点后指定位数

格式化输出 格式化输出所用的函数为 printf&#xff0c;它可以输出任意位数的小数。 使用格式&#xff1a;printf(“%.nf”,a)。这句话的作用是将变量a保留n位小数输出。 注意事项&#xff1a; 1、这里的n&#xff0c;需要具体化为一个数字&#xff0c;保留几位小数&#x…

【ARMv8/v9 GIC 系列 3 -- GIC 的 类型寄存器 GICD_TYPER】

文章目录 GIC 类型寄存器 GICD_TYPERESPI_Range, 位[31:27]RSS, 位[26]No1N, 位[25]A3V, 位[24]IDBits, 位[23:19]DVIS, 位[18]LPIs, 位[17]MBIS, 位[16]NUM_LPIs, 位[15:11]SecurityExtn, 位[10]NMI, 位[9]ESPI, 位[8]CPUNumber, 位[7:5]ITLinesNumber, 位[4:0]GIC 类型寄存器…

朗科HD10M2Pr震撼上市,自带风扇极速降温,匹敌私有云

近日,存储领域的领军企业朗科旗下全资子公司朗科创新宣布,其最新款磁吸硬盘盒HD10M2Pr正式上市。这款产品凭借超薄设计、极速降温、高速传输等多项优势,迅速成为了行业内的讨论焦点。 随着移动设备使用的普及和短视频内容的日益丰富,对于存储空间不断增长的需求逐渐成为日常生活…

导入导出带下拉框模版(EasyExcel)

前言 项目进行到新的一个迭代了&#xff0c;赶了1周需求&#xff0c;接口终于处理完了。分享记录下迭代中处理导入、导出、下载模版功能的细节吧。 一、场景 EasyExcel&#xff08;阿里&#xff09;实现Excel数据处理三层表头&#xff0c;第二、三层表头动态数据根据第二、三层…

RabbitMQ(六)仲裁队列、流式队列、异地容灾(联邦队列Federation Queue)

文章目录 仲裁队列1、创建交换机2、创建仲裁队列3、验证主节点宕机不影响消息发送和接收 流式队列&#xff08;不推荐&#xff0c;Kafka主场&#xff09;概念 异地容灾一、Federation插件概述 二、Federation交换机1、总体说明2、准备工作3、启用联邦插件4、添加上游连接端点5、…

NetSuite Inventory Transfer Export Saved Search

用户之前有提出一个实际的需求&#xff0c;大致意思是想要导出Inventory Transfer的相关明细行信息&#xff0c;且要包含From Location&#xff0c;To Location&#xff0c;Quantity等信息。 我们知道From Location和To Location在IT Form中应该是在Main的部分&#xff0c;在D…