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

news2024/10/16 15:26:22

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 ps

t = 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 plt

def 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.coloradmin.cn/o/1843570.html

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

相关文章

【etcd】etcd单机安装及简单操作

https://blog.csdn.net/Mr_XiMu/article/details/125026635 https://blog.csdn.net/m0_73192864/article/details/136509244 etcd在生产环境中一般为集群方式部署 etcd使用的2个默认端口号:2379和2380 2379:用于客户端通信(类似于sqlserver的1433&#x…

视频融合共享平台LntonCVS视频监控安防系统运用多视频协议建设智慧园区方案

智慧园区,作为现代化城市发展的重要组成部分,不仅推动了产业的升级转型,也成为了智慧城市建设的核心力量。随着产业园区之间的竞争日益激烈,如何打造一个功能完善、智能化程度高的智慧园区,已经成为了业界广泛关注的焦…

五十、openlayers官网示例JSTS Integration解析——使用JSTS 库来处理几何缓冲区并在地图上显示结果

官网demo地址: JSTS Integration 这篇讲了如何在地图上添加缓冲图形 什么叫做缓冲几何? 几何缓冲(Geometric Buffering)是指在 GIS(地理信息系统)和计算几何中,围绕一个几何对象创建一个具有…

时空预测 | 基于深度学习的碳排放时空预测模型

时空预测 模型描述 数据收集和准备:收集与碳排放相关的数据,包括历史碳排放数据、气象数据、人口密度数据等。确保数据的质量和完整性,并进行必要的数据清洗和预处理。 特征工程:根据问题的需求和领域知识,对数据进行…

Walrus:去中心化存储和DA协议,可以基于Sui构建L2和大型存储

Walrus是为区块链应用和自主代理提供的创新去中心化存储网络。Walrus存储系统今天以开发者预览版的形式发布,面向Sui开发者征求反馈意见,并预计很快会向其他Web3社区广泛推广。 通过采用纠删编码创新技术,Walrus能够快速且稳健地将非结构化数…

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

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

第一个Neety程序

&#x1f4dd;个人主页&#xff1a;五敷有你 &#x1f525;系列专栏&#xff1a;Netty ⛺️稳中求进&#xff0c;晒太阳 加入依赖 <dependency><groupId>io.netty</groupId><artifactId>netty-all</artifactId><version>4.1.39.F…

【Spine学习10】之 创建新骨骼时,自动绑定图片和插槽的快捷方式

两天没更新了。 遇到一些难解的难题 用的版本是破解版 不知道为啥现在的教程非常地快 明明有些细节很重要还略过讲 所以创建骨骼这里 基本创建是都会 可是骨骼一多 实际工作中的重命名也太麻烦了 。 这就需要学习快捷创建方式&#xff1a; <将对应图片自动绑定到新骨骼上并…

哔哩哔哩视频URL解析原理

哔哩哔哩视频URL解析原理 视频网址解析视频的原理通常涉及以下几个步骤&#xff1a; 1、获取视频页面源代码&#xff1a;通过HTTP请求获取视频所在网页的HTML源代码。这一步通常需要处理反爬虫机制&#xff0c;如验证码或用户登录。 2、解析页面源代码&#xff1a;分析HTML源代…

【上海交大】博士生年度进展报告模板

上海交通大学 博士生年度进展报告模板 比较不好找&#xff0c;在交我办中发起申请流程后才能看到链接&#xff0c;链接如下&#xff1a; https://www.gs.sjtu.edu.cn/xzzx/pygl/15

byte[]转MultipartFile、byte[]转File一次看个够

目录 需求背景 当你需要将byte[]、MultipartFile、File实现互转时&#xff0c;无外乎以下场景&#xff1a; 保存第三方接口返回二进制流前/后端文件流上传微服务间调用文件格式转换 正如你所需要的&#xff0c;通过搜索引擎筛选到我的本篇文章是因为你在开发中需要将byte[]转…

剑指offer 算法题(数组中重复的数据)

剑指offer 第一题 去力扣里测试算法 思路一&#xff1a; 排序后&#xff0c;前一个与后一个相比是否相同。 class Solution { public:vector<int> findDuplicates(vector<int>& nums) {sort(nums.begin(), nums.end());int n 0;int len nums.size();vect…

如何设计一个点赞系统

首先我们定义出一个点赞系统需要对外提供哪些接口&#xff1a; 1.用户对特定的消息进行点赞&#xff1b; 2.用户查看自己发布的某条消息点赞数量以及被哪些人赞过&#xff1b; 3.用户查看自己给哪些消息点赞过&#xff1b; 这里假设每条消息都有一个message_id, 每一个用户都…

文件防篡改监控工具 - WGCLOUD全面介绍

WGCLOUD是一款优秀的运维监控软件&#xff0c;免费、轻量、高效&#xff0c;部署容易&#xff0c;上手简单&#xff0c;对新手非常友好 WGCLOUD部署完成后&#xff0c;点击菜单【文件防篡改】&#xff0c;可以看到如下页面 我们点击【添加】按钮&#xff0c;输入监控文件的信息…

python之Bible快速检索器

内容将会持续更新&#xff0c;有错误的地方欢迎指正&#xff0c;谢谢! python之Bible快速检索器 TechX 坚持将创新的科技带给世界&#xff01; 拥有更好的学习体验 —— 不断努力&#xff0c;不断进步&#xff0c;不断探索 TechX —— 心探索、心进取&#xff01; 助力快…

【网络安全的神秘世界】已解决Failed to start proxy service on 127.0.0.1:8080

&#x1f31d;博客主页&#xff1a;泥菩萨 &#x1f496;专栏&#xff1a;Linux探索之旅 | 网络安全的神秘世界 | 专接本 | 每天学会一个渗透测试工具 解决burpsuite无法在 127.0.0.1&#xff1a;8080 上启动代理服务端口被占用以及抓不到本地包的问题 Burpsuite无法启动proxy…

ISSN和CN到底有什么不同呢?

在学术出版领域&#xff0c;我们常常会遇到各种标识码&#xff0c;其中ISSN和CN无疑是两个最为常见的。尽管它们都是用于标识出版物的重要信息&#xff0c;但两者在定义、功能和应用上却有着显著的区别。那么&#xff0c;ISSN和CN到底有什么不同呢&#xff1f;接下来&#xff0…

深度解析RocketMq源码-持久化组件(四) CommitLog

1.绪论 commitLog是rocketmq存储的核心&#xff0c;前面我们介绍了mappedfile、mappedfilequeue、刷盘策略&#xff0c;其实commitlog的核心组件我们基本上已经介绍完成。 2.commitLog的组成 commitLog的核心其实就是MqppedFilequeue&#xff0c;它本质上就是多个mappedFile…

图像数字化基础

一、像素 1、获取图像指定位置的像素 import cv2 image cv2.imread("E:\\images\\2.png") px image[291,218] print("坐标(291,218)上的像素的BGR值是&#xff1a;",px) &#xff08;1&#xff09;RGB色彩空间 R通道&#xff1a;红色通道 G通道&…

LVS(Linux Virtual Server)集群,(1)NAT模式

Cluster&#xff1a;集群&#xff0c;为了解决某个特定问题将多台计算机组合起来形成的单个系统。 集群分为三种类型&#xff1a; LB(Load Balancing)&#xff0c;负载均衡&#xff0c;多个主机组成&#xff0c;每个主机只承担一部分访问请求 HA(High Availiablity)&#xf…