针对循环神经网络RNN梯度消失的缺点，造成早期的序列数据信息在迭代过程中丢失，产生了长短记忆网络（long short time memory，LSTM）的 RNN。

基本思想

与 RNN 相比，LSTM 多了个遗忘与记忆功能：

• $h_t$ 表示短期记忆
• $c_t$ 表示长期记忆
• $g_t$ 表示候选的新记忆
• $f_t$ 遗忘门，表示决定保留多少旧记忆，靠近0或1的一个数
• $i_t$ 输入门，表示决定要接收多少新记忆，靠近0或1的一个数
• $o_t$ 输出门，表示决定要输出多少当前的所有记忆，靠近0或1的一个数

LSTM 的基本思想由下面这个公式体现：
$$c_t=f_t\odot c_{t-1}+i_t\odot g_t$$

其中：

• $\odot$ 表示哈达玛乘积 (Hadamard product)，其输入为两个相同形状的矩阵或向量，输出是具有同样形状的、各个位置的元素等于两个输入矩阵相同位置元素的乘积的矩阵或向量.
• 哈达玛乘积与向量的点积不同，点积是对两个向量逐元素相乘再求和
• 这个公式表示新的长期记忆为保留的旧记忆与输入的新记忆的和

最终的输出信息为：
$$h_t=o_t\odot \text{tanh}(c_t)$$

其中：

• tanh 为激活函数 tanh (常见的激活函数介绍见这篇博客)
• 在 lstm 网络中，每个时间步的输出信息一般为 $h_t$，等同于 RNN 或传统神经网络中的 ${\hat{y}}_t$.

下图是一个经典的 LSTM cell 的内部结构：

在该图中：

• $\odot$ 表示哈达玛乘积

• $+$ 表示相加

• 输入 $x_t$, $h_{t-1}$, $c_{t-1}$ (若不是第一个隐含层，则输入只有 $h_{t-1},c_{t-1}$)，输出 $h_t,c_t$

• $\sigma$ 为 sigmoid 函数

• tanh 为 tanh 激活函数

遗忘门 $f_t$ 的表达式为：
$$f_t=\sigma (W_f[h_{t-1},x_t]+b_f)$$

• $[h_{t-1},x_t]$ 为列向量 $h_{t-1}$ 与 $x_t$ 的上下拼接

输入门 $i_t$ 的表达式为：
$$i_t=\sigma (W_i[h_{t-1},x_t]+b_i)$$

输出门 $o_t$ 的表达式为：
$$o_t=\sigma (W_o[h_{t-1},x_t]+b_o)$$

候选的新记忆 $g_t$ 的表达式为：
$$g_t=\sigma (W_c[h_{t-1},x_t]+b_c)$$

下图是一个时间步长为 3，单层的 LSTM 结构图：

下面是一个双层的 LSTM 结构：

其中，上标表示层数。

• 同一时间步内，也可以出现多个 LSTM cell 并行的情况，此时它们共享输入

步骤

下面，我们推导出广义的权重更新公式，考虑 LSTM 可以有多层 LSTM cell，用上标表示层数，假设一共有 $L$ 层。

1 符号约定

• 层数：$l=1,2,\dots ,L$

• 时间步：$t=1,2,\dots ,T$

• 第 $l$ 层的隐藏状态与细胞状态：
  $$h_t^{(l)}\in {\mathbb{R}}^n,c_t^{(l)}\in {\mathbb{R}}^n$$

• 第 $l$ 层输入：
  $$\boxed{x_t^{(l)}=\{\begin{array}{ll}{x}_t& l=1\\ h_t^{(l-1)}& l>1\end{array}}$$

• 参数：
  $$W^{(l)}\in {\mathbb{R}}^{4n\times m_l},U^{(l)}\in {\mathbb{R}}^{4n\times n},b^{(l)}\in {\mathbb{R}}^{4n},W_y\in {\mathbb{R}}^{n\times m_o}$$
  其中 $m_l$ 是第 $l$ 层输入的维度，$m_o$ 为最后输出层的维度。

2 前向传播

将传播公式用矩阵形式写出来：

$$\boxed{\begin{array}{rl}{a}_t^{(l)}& =\left[\begin{array}{c}{a}_{i,t}^{(l)}\\ a_{f,t}^{(l)}\\ a_{o,t}^{(l)}\\ a_{g,t}^{(l)}\end{array}\right]=W^{(l)}{x}_t^{(l)}+U^{(l)}{h}_{t-1}^{(l)}+b^{(l)},\\ i_t^{(l)}& =\sigma (a_{i,t}^{(l)}),f_t^{(l)}=\sigma (a_{f,t}^{(l)}),\\ o_t^{(l)}& =\sigma (a_{o,t}^{(l)}),g_t^{(l)}=\tanh (a_{g,t}^{(l)}),\\ c_t^{(l)}& =f_t^{(l)}\odot c_{t-1}^{(l)}+i_t^{(l)}\odot g_t^{(l)},\\ h_t^{(l)}& =o_t^{(l)}\odot \tanh (c_t^{(l)})\\ a_t^y& =W_t^y{h}_t^{(l)}+b_t^y\end{array}}$$

3 反向传播（BPTT through layers and time）

定义反向梯度：
$${\bar{h}}_t^{(l)}=\frac{\partial \mathcal{L}}{\partial h_t^{(l)}},{\bar{c}}_t^{(l)}=\frac{\partial \mathcal{L}}{\partial c_t^{(l)}},{\delta }_t^{(l)}=\frac{\partial \mathcal{L}}{\partial a_t^{(l)}}=\left[\begin{array}{c}{\delta }_{i,t}^{(l)}\\ {\delta }_{f,t}^{(l)}\\ {\delta }_{o,t}^{(l)}\\ {\delta }_{g,t}^{(l)}\end{array}\right],{\delta }_t^{(L+1)}=\frac{\partial \mathcal{L}}{\partial a_t^y}$$

上划线表示对中间状态（输出 $h$, $c$）的梯度

3.1 梯度递推关系

• 对隐藏状态的反向传递：
  $$\boxed{{\bar{h}}_t^{(l)}=\frac{\partial \mathcal{L}}{\partial h_t^{(l)}}{|}_{direct}+(U^{(l)}{)}^{\top }{\delta }_{t+1}^{(l)}+(W^{(l+1)}{)}^{\top }{\delta }_t^{(l+1)}}$$

• 第一项为当前时间步，当前层的直接损失；若损失函数不直接依赖于该层输出，则此项为 0
• 第二项为来自同层下一时间步的反向传播
• 第三项为来自下一层同一时间步的反向传播
• 第二项或第三项其实是 $f,g,i,o$ 四个门相应矩阵乘积的加和

• 对细胞状态的反向传递：
  $$\boxed{{\bar{c}}_t^{(l)}=\frac{\partial \mathcal{L}}{\partial c_t^{(l)}}{|}_{direct}+{\bar{h}}_t^{(l)}\odot o_t^{(l)}\odot (1-{\tanh }^2(c_t^{(l)}))+{\bar{c}}_{t+1}^{(l)}\odot f_{t+1}^{(l)}}$$

• 第一项通常为 0， 因为 $c_t$ 不直接影响当前的损失
• 第二项是公式 $h_t^{(l)}=o_t^{(l)}\odot \tanh (c_t^{(l)})$ 对 $c_t^{(l)}$ 的求导
• 第三项是公式 $c_{t+1}^{(l)}=f_{t+1}^{(l)}\odot c_t^{(l)}+i_{t+1}^{(l)}\odot g_{t+1}^{(l)}$ 对 $c_t^{(l)}$ 的求导

• 对输出层的 pre-activation 梯度以及临界值：
  $$\boxed{\begin{array}{rl}{\delta }_t^{(L+1)}& ={\hat{y}}_t-y_t\\ {\delta }_{T+1}^{(l)}& =0\end{array}}$$
• 对各门的 pre-activation 梯度：
  $$\boxed{\begin{array}{rl}{\delta }_{i,t}^{(l)}& ={\bar{c}}_t^{(l)}\odot g_t^{(l)}\odot i_t^{(l)}\odot (1-i_t^{(l)}),\\ {\delta }_{f,t}^{(l)}& ={\bar{c}}_t^{(l)}\odot c_{t-1}^{(l)}\odot f_t^{(l)}\odot (1-f_t^{(l)}),\\ {\delta }_{o,t}^{(l)}& ={\bar{h}}_t^{(l)}\odot \tanh (c_t^{(l)})\odot o_t^{(l)}\odot (1-o_t^{(l)}),\\ {\delta }_{g,t}^{(l)}& ={\bar{c}}_t^{(l)}\odot i_t^{(l)}\odot (1-(g_t^{(l)}{)}^2).\end{array}}$$

3.2 权重梯度累积公式

对第 $l$ 层：
$$\boxed{\begin{array}{rl}\frac{\partial \mathcal{L}}{\partial W^{(l)}}& =\sum _{t=1}^T{\delta }_t^{(l)}(x_t^{(l)}{)}^{\top },\\ \frac{\partial \mathcal{L}}{\partial U^{(l)}}& =\sum _{t=1}^T{\delta }_t^{(l)}(h_{t-1}^{(l)}{)}^{\top },\\ \frac{\partial \mathcal{L}}{\partial b^{(l)}}& =\sum _{t=1}^T{\delta }_t^{(l)},\\ \frac{\partial \mathcal{L}}{\partial W^y}& =\sum _{t=1}^T{\delta }_t^{(L+1)}(h_t^{(L)}{)}^{\top },\\ \frac{\partial \mathcal{L}}{\partial b^y}& =\sum _{t=1}^T{\delta }_t^{(l)}\end{array}}$$

编程实现时，上面几个圈中的公式比较重要。

3.3 参数更新

以 SGD 为例

$$\begin{array}{rl}{W}^{(l)}& \leftarrow W^{(l)}-\eta \frac{\partial \mathcal{L}}{\partial W^{(l)}},\\ U^{(l)}& \leftarrow U^{(l)}-\eta \frac{\partial \mathcal{L}}{\partial U^{(l)}},\\ b^{(l)}& \leftarrow b^{(l)}-\eta \frac{\partial \mathcal{L}}{\partial b^{(l)}}.\end{array}$$

ADAM 的更新情况可以类似推出。

4 bidirection 双向

lSTM 还可以双向传播，即每个时间步输出拼接了前向 + 后向 hidden：

• 双向 LSTM 的传播是两个方向同时计算，然后将输出拼接，即

每个时间步包含前向 + 后向信息

• h_t_forward: 从序列头到 t 的信息

• h_t_backward: 从序列尾到 t 的信息

拼接后，h_t = [h_t_forward ; h_t_backward]

具体实现顺序（伪算法）

1. 前向传播（Forward Pass）

   • 按时间 $t=1,\dots ,T$ 和层 $l=1,\dots ,L$ 遍历。
   • 对第 1 层输入为 $x_t$，对第 $l$ 层（$l>1$）输入为上一层的隐藏状态 $h_t^{(l-1)}$
   • 计算每个 LSTM 单元的输入门、遗忘门、输出门和候选状态门的线性组合及非线性激活。
   • 更新细胞状态 $c_t^{(l)}$ 和隐藏状态 $h_t^{(l)}$。
   • 保存所有中间量（门激活、线性组合、细胞状态、隐藏状态等），以便反向传播使用。

2. 初始化梯度容器

   • 对每层的权重矩阵 $W^{(l)},U^{(l)},b^{(l)}$ 初始化梯度累加器为 0。
   • 对时间 $T+1$ 的隐藏状态梯度和细胞状态梯度初始化为 0。

3. 反向传播（Backward Pass）

   • 按层从 L→1，按时间从 T→1 遍历。
   • 计算每个时间步每层的隐藏状态和细胞状态的累积梯度。
   • 根据链式法则计算各门的梯度（输入门、遗忘门、输出门、候选门）。
   • 将各门梯度合并形成 LSTM 单元的总梯度。

4. 累积梯度

   • 使用门的梯度和对应的输入或上一时刻隐藏状态累积到权重梯度 $d{W}^{(l)},d{U}^{(l)},d{b}^{(l)}$。

5. 参数更新

   • 根据 SGD 或者 ADAM 更新参数。

一个 C++ 例子

用 C++ 代码编写 LSTM，预测 airpassengers 数据集，用了 Eigen 库进行一些矩阵运算。

预测效果图：

误差

Total absolute error 1174.41 RMSE: 11.3162, MAE(mean absolute error):
8.89707

/**
 * Author: Zhen Chen
 * Email: chen.zhen5526@gmail.com
 * Created on: 01/11/2025, 20:14
 * Description:
 *
 */
#include <Eigen/Dense>
#include <cmath>
#include <iostream>
#include <vector>

using namespace Eigen;

#include "../utils/matplotlibcpp.h"
namespace plt = matplotlibcpp;

void draw_pic(const std::vector<double> &prediction, const std::vector<double> &true_value) {
  // 强制使用 TkAgg 后端，避免 Qt6Agg 卡死
  plt::backend("TkAgg"); // 等价于 python 中的命令 matplotlib.use("TkAgg")
  std::vector<double> x(prediction.size());
  for (int i = 0; i < prediction.size(); ++i) {
    x[i] = static_cast<double>(i);
  }

  plt::plot(x, true_value, {{"label", "true value"}});
  plt::plot(x, prediction, {{"label", "prediction value"}});
  plt::grid(true);
  plt::legend(); // 显示图例
  plt::show();
}

// Sigmoid
MatrixXd sigmoid(const MatrixXd &x) { return 1.0 / (1.0 + (-x.array()).exp()); }
MatrixXd sigmoid_derivative(const MatrixXd &x) { return x.array() * (1.0 - x.array()); }
MatrixXd tanh_derivative(const MatrixXd &x) { return 1.0 - x.array().tanh().square(); }

class LSTM {
public:
  int input_size, hidden_size;
  double lr;
  double beta1 = 0.9, beta2 = 0.999, epsilon = 1e-8;
  int t_adam = 0;

  // LSTM 参数
  MatrixXd W_i, W_f, W_o, W_g;
  MatrixXd U_i, U_f, U_o, U_g;
  VectorXd b_i, b_f, b_o, b_g;

  // Adam 一阶、二阶矩
  MatrixXd m_W_i, m_W_f, m_W_o, m_W_g;
  MatrixXd m_U_i, m_U_f, m_U_o, m_U_g;
  VectorXd m_b_i, m_b_f, m_b_o, m_b_g;
  MatrixXd v_W_i, v_W_f, v_W_o, v_W_g;
  MatrixXd v_U_i, v_U_f, v_U_o, v_U_g;
  VectorXd v_b_i, v_b_f, v_b_o, v_b_g;

  // 用于存储前向传播中各变量数据
  std::vector<MatrixXd> h;
  std::vector<MatrixXd> c;
  std::vector<MatrixXd> i;
  std::vector<MatrixXd> f;
  std::vector<MatrixXd> o;
  std::vector<MatrixXd> g;

  LSTM(const int input_size_, const int hidden_size_, const double lr_ = 0.001)
      : input_size(input_size_), hidden_size(hidden_size_), lr(lr_) {

    W_i = MatrixXd::Random(hidden_size, input_size) * 0.1;
    W_f = MatrixXd::Random(hidden_size, input_size) * 0.1;
    W_o = MatrixXd::Random(hidden_size, input_size) * 0.1;
    W_g = MatrixXd::Random(hidden_size, input_size) * 0.1;

    U_i = MatrixXd::Random(hidden_size, hidden_size) * 0.1;
    U_f = MatrixXd::Random(hidden_size, hidden_size) * 0.1;
    U_o = MatrixXd::Random(hidden_size, hidden_size) * 0.1;
    U_g = MatrixXd::Random(hidden_size, hidden_size) * 0.1;

    b_i = VectorXd::Zero(hidden_size);
    b_f = VectorXd::Zero(hidden_size);
    b_o = VectorXd::Zero(hidden_size);
    b_g = VectorXd::Zero(hidden_size);

    // Adam
    m_W_i = MatrixXd::Zero(hidden_size, input_size);
    m_W_f = m_W_i;
    m_W_o = m_W_i;
    m_W_g = m_W_i;
    m_U_i = MatrixXd::Zero(hidden_size, hidden_size);
    m_U_f = m_U_i;
    m_U_o = m_U_i;
    m_U_g = m_U_i;
    m_b_i = VectorXd::Zero(hidden_size);
    m_b_f = m_b_i;
    m_b_o = m_b_i;
    m_b_g = m_b_i;
    v_W_i = m_W_i;
    v_W_f = m_W_i;
    v_W_o = m_W_i;
    v_W_g = m_W_i;
    v_U_i = m_U_i;
    v_U_f = m_U_i;
    v_U_o = m_U_i;
    v_U_g = m_U_i;
    v_b_i = m_b_i;
    v_b_f = m_b_i;
    v_b_o = m_b_i;
    v_b_g = m_b_i;
  }

  void forward(const MatrixXd &x) {
    const size_t T = x.rows(); // x 的每一行是一个时间步
    if (h.size() != T) {
      h.resize(T);
      c.resize(T);
      i.resize(T);
      f.resize(T);
      o.resize(T);
      g.resize(T);
    }
    for (size_t t = 0; t < T; t++) {
      h[t] = MatrixXd::Zero(hidden_size, 1);
      c[t] = MatrixXd::Zero(hidden_size, 1);
      i[t] = MatrixXd::Zero(hidden_size, 1);
      f[t] = MatrixXd::Zero(hidden_size, 1);
      o[t] = MatrixXd::Zero(hidden_size, 1);
      g[t] = MatrixXd::Zero(hidden_size, 1);
    }

    // std::cout << std::endl; // 换行并刷新
    auto last_c = c[0];
    auto last_h = h[0];
    for (size_t t = 0; t < T; ++t) {
      auto x_t = x.row(static_cast<Index>(t)).transpose();
      i[t] = sigmoid(W_i * x_t + U_i * last_h + b_i);
      f[t] = sigmoid(W_f * x_t + U_f * last_h + b_f);
      o[t] = sigmoid(W_o * x_t + U_o * last_h + b_o);
      // .array() 的作用是将一个矩阵或向量 从线性代数（matrix）对象转换为数组（Array）对象，
      // 从而可以进行 逐元素（element-wise）操作
      // .tanh() 是 Eigen 数组对象自带的逐元素双曲正切函数
      g[t] = (W_g * x_t + U_g * last_h + b_g).array().tanh();
      c[t] = i[t].array() * g[t].array() + f[t].array() * last_c.array();
      h[t] = o[t].array() * c[t].array().tanh();

      last_c = c[t];
      last_h = h[t];
    }
  }

  void backward_with_dh(const MatrixXd &x, const MatrixXd &dh_external) {
    const size_t T = x.rows();
    MatrixXd dW_i = MatrixXd::Zero(hidden_size, input_size);
    MatrixXd dW_f = MatrixXd::Zero(hidden_size, input_size);
    MatrixXd dW_o = MatrixXd::Zero(hidden_size, input_size);
    MatrixXd dW_g = MatrixXd::Zero(hidden_size, input_size);

    MatrixXd dU_i = MatrixXd::Zero(hidden_size, hidden_size);
    MatrixXd dU_f = MatrixXd::Zero(hidden_size, hidden_size);
    MatrixXd dU_o = MatrixXd::Zero(hidden_size, hidden_size);
    MatrixXd dU_g = MatrixXd::Zero(hidden_size, hidden_size);

    VectorXd db_i = VectorXd::Zero(hidden_size);
    VectorXd db_f = VectorXd::Zero(hidden_size);
    VectorXd db_o = VectorXd::Zero(hidden_size);
    VectorXd db_g = VectorXd::Zero(hidden_size);

    MatrixXd delta_h = MatrixXd::Zero(hidden_size, 1);
    MatrixXd delta_c = MatrixXd::Zero(hidden_size, 1);
    MatrixXd delta_f = MatrixXd::Zero(hidden_size, 1);
    MatrixXd delta_g = MatrixXd::Zero(hidden_size, 1);
    MatrixXd delta_i = MatrixXd::Zero(hidden_size, 1);
    MatrixXd delta_o = MatrixXd::Zero(hidden_size, 1);
    for (int t = static_cast<int>(T - 1); t >= 0; --t) {
      if (t == T - 1) {
        delta_h = dh_external;
      } else {
        delta_h = U_f.transpose() * delta_f + U_g.transpose() * delta_g +
                  U_o.transpose() * delta_o + U_i.transpose() * delta_i;
      }
      if (t == T - 1)
        delta_c = delta_h.array() * o[t].array() * (1 - c[t].array().tanh().square());
      else {
        delta_c = delta_h.array() * o[t].array() * (1 - c[t].array().tanh().square()) +
                  delta_c.array() * f[t + 1].array();
      }

      delta_i = delta_c.array() * g[t].array() * i[t].array() * (1 - i[t].array());
      delta_f = delta_c.array() * c[t].array() * f[t].array() * (1 - f[t].array());
      delta_o = delta_h.array() * c[t].array().tanh() * o[t].array() * (1 - o[t].array());
      delta_g = delta_c.array() * i[t].array() * (1 - g[t].array().square());

      dW_i += delta_i * x.row(t).transpose(); // 这个 x[t] 不需要转置，因为x[t]是一个行向量
      if (t > 0)
        dU_i += delta_i * h[t - 1].transpose();
      db_i += delta_i.col(0);
      dW_f += delta_f * x.row(t).transpose();
      if (t > 0)
        dU_f += delta_f * h[t - 1].transpose();
      db_f += delta_f.col(0);
      dW_o += delta_o * x.row(t).transpose();
      if (t > 0)
        dU_o += delta_o * h[t - 1].transpose();
      db_o += delta_o.col(0);
      dW_g += delta_g * x.row(t).transpose();
      if (t > 0)
        dU_g += delta_g * h[t - 1].transpose();
      db_g += delta_g.col(0);
    }

    // // SDG updating
    // W_i -= lr * dW_i;
    // U_i -= lr * dU_i;
    // b_i -= lr * db_i;
    // W_f -= lr * dW_f;
    // U_f -= lr * dU_f;
    // b_f -= lr * db_f;
    // W_g -= lr * dW_g;
    // U_g -= lr * dU_g;
    // b_g -= lr * db_g;
    // W_o -= lr * dW_o;
    // U_o -= lr * dU_o;
    // b_o -= lr * db_o;

    // ADAM updating
    t_adam++;
    auto adam_update = [this](MatrixXd &param, MatrixXd &m, MatrixXd &v, const MatrixXd &grad) {
      m = beta1 * m + (1 - beta1) * grad;
      v = beta2 * v + (1 - beta2) * grad.array().square().matrix();
      MatrixXd m_hat = m / (1 - pow(beta1, t_adam));
      MatrixXd v_hat = v / (1 - pow(beta2, t_adam));
      MatrixXd update = lr * m_hat.array() / (v_hat.array().sqrt() + epsilon).array();
      param = param - update;
    };
    auto adam_update_vec = [this](VectorXd &param, VectorXd &m, VectorXd &v, const VectorXd &grad) {
      m = beta1 * m + (1 - beta1) * grad;
      v = beta2 * v + (1 - beta2) * grad.array().square().matrix();
      VectorXd m_hat = m / (1 - pow(beta1, t_adam));
      VectorXd v_hat = v / (1 - pow(beta2, t_adam));
      VectorXd update = lr * m_hat.array() / (v_hat.array().sqrt() + epsilon).array();
      param = param - update;
    };
    adam_update(W_i, m_W_i, v_W_i, dW_i);
    adam_update(U_i, m_U_i, v_U_i, dU_i);
    adam_update_vec(b_i, m_b_i, v_b_i, db_i);
    adam_update(W_f, m_W_f, v_W_f, dW_f);
    adam_update(U_f, m_U_f, v_U_f, dU_f);
    adam_update_vec(b_f, m_b_f, v_b_f, db_f);
    adam_update(W_o, m_W_o, v_W_o, dW_o);
    adam_update(U_o, m_U_o, v_U_o, dU_o);
    adam_update_vec(b_o, m_b_o, v_b_o, db_o);
    adam_update(W_g, m_W_g, v_W_g, dW_g);
    adam_update(U_g, m_U_g, v_U_g, dU_g);
    adam_update_vec(b_g, m_b_g, v_b_g, db_g);
  }
};

int main() {
  // AirPassengers 数据
  std::vector<double> data = {
      112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118, 115, 126, 141, 135, 125, 149,
      170, 170, 158, 133, 114, 140, 145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166,
      171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194, 196, 196, 236, 235, 229, 243,
      264, 272, 237, 211, 180, 201, 204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229,
      242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278, 284, 277, 317, 313, 318, 374,
      413, 405, 355, 306, 271, 306, 315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336,
      340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337, 360, 342, 406, 396, 420, 472,
      548, 559, 463, 407, 362, 405, 417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432};

  const double max_val = *std::ranges::max_element(data);
  const double min_val = *std::ranges::min_element(data);
  for (auto &v : data)
    v = (v - min_val) / (max_val - min_val);

  constexpr int input_size = 1;
  constexpr int seq_size = 12; // 用过去12月作为输入
  constexpr int hidden_size = 10;
  constexpr double lr = 0.01;
  LSTM lstm(input_size, hidden_size, lr);

  constexpr int epochs = 300;
  const int train_len = static_cast<int>(data.size()) - 12;

  MatrixXd W_y = MatrixXd::Random(1, hidden_size) * 0.1;
  MatrixXd b_y = MatrixXd::Zero(1, 1);

  for (int e = 0; e < epochs; e++) {
    double loss = 0.0;
    for (int t = 0; t < train_len; t++) {
      // 构建输入 x = 过去12个月
      MatrixXd x(seq_size, 1);
      for (int k = 0; k < seq_size; k++) {
        x(k, 0) = data[t + k];
      }
      MatrixXd y(1, 1);
      y(0, 0) = data[t + seq_size]; // 预测下一个月

      // 前向传播
      lstm.forward(x);

      // 输出层前向
      MatrixXd y_pred = W_y * lstm.h[seq_size - 1] + b_y;
      MatrixXd dy = y_pred - y;
      loss += dy(0, 0) * dy(0, 0);

      // 反向传播到隐藏层
      MatrixXd dWy = dy * lstm.h[seq_size - 1].transpose();
      MatrixXd delta_h = W_y.transpose() * dy;

      lstm.backward_with_dh(x, delta_h);

      // 输出层梯度更新
      W_y -= lr * dWy;
      b_y -= lr * dy;
    }
    if (e % 50 == 0)
      std::cout << "Epoch " << e << ", Loss: " << loss / train_len << std::endl;
  }

  // 使用滑动窗口预测整个序列
  std::vector<double> predictions(data.size() - seq_size);
  std::vector<double> true_values(data.size() - seq_size);
  for (int t = 0; t < train_len; t++) {
    MatrixXd x(seq_size, 1);
    for (int k = 0; k < seq_size; k++)
      x(k, 0) = data[t + k];

    lstm.forward(x);

    const double y_pred = (W_y * lstm.h[seq_size - 1] + b_y)(0, 0);
    predictions[t] = y_pred * (max_val - min_val) + min_val; // 反归一化
  }

  // 计算 RMSE 与 MAD
  double rmse = 0.0, mad = 0.0;
  for (int t = 0; t < train_len; t++) {
    const double true_val = data[t + seq_size] * (max_val - min_val) + min_val;
    const double error = predictions[t] - true_val;
    true_values[t] = true_val;
    rmse += error * error;
    mad += fabs(error);
    std::cout << "idx=" << t << " true = " << true_val << " pred = " << predictions[t] << std::endl;
  }
  rmse = sqrt(rmse / train_len);
  std::cout << std::endl;
  std::cout << "Total absolute error " << mad << std::endl;
  mad /= train_len;

  std::cout << "RMSE: " << rmse << ", MAE(mean absolute error): " << mad << std::endl;

  // 画图
  draw_pic(predictions, true_values);

  return 0;
}

一个 Python

上面的例子，用 python 编程：

"""
Python version: 3.12.7
Author: Zhen Chen, chen.zhen5526@gmail.com
Date: 2025/10/20 20:34
Description: lstm, criteria L1Loss() is better than MSELoss()


"""

import pandas as pd
import torch
import torch.nn as nn
import numpy as np
from statsmodels.datasets import get_rdataset


# -------------------------------
# 加载数据
# -------------------------------
dataset = get_rdataset("AirPassengers").data
raw_data = dataset["value"].values.astype(float)

# 归一化
data_min = raw_data.min()
data_max = raw_data.max()
data = (raw_data - data_min) / (data_max - data_min)


# -------------------------------
# 创建序列
# -------------------------------
def create_sequences(data, seq_length):
    xs, ys = [], []
    for i in range(len(data) - seq_length):
        x = data[i : i + seq_length]
        y = data[i + seq_length]
        xs.append(x)
        ys.append(y)
    return np.array(xs), np.array(ys)


seq_length = 12  # 用过去一年数据预测下一月
X, y = create_sequences(data, seq_length)
X = torch.tensor(X, dtype=torch.float32).unsqueeze(-1)  # 对于一维数据，增加维度
y = torch.tensor(y, dtype=torch.float32).unsqueeze(-1)


# -------------------------------
# 定义 LSTM
# -------------------------------
class LSTMModel(nn.Module):
    def __init__(self, input_size=1, hidden_size=10, num_layers=1, output_size=1):
        super().__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.linear = nn.Linear(hidden_size, output_size)  # 线性输出

    def forward(self, x):
        # 类中有 __call__() 函数，所以类可以直接调用
        out, _ = self.lstm(x)
        # out 的维度 (batch_size, seq_length, hidden_size)
        out = out[:, -1, :]
        out = self.linear(out)
        return out


model = LSTMModel()
criterion = nn.L1Loss()  # nn.MSELoss()  #
# model.parameters() 包含模型训练的各项权重与偏置
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

# -------------------------------
# 训练
# -------------------------------
epochs = 300
for epoch in range(epochs):
    optimizer.zero_grad()  # 每次循环时梯度清零，不累加；模拟大 batch 时累加
    output = model(X)

    loss = criterion(output, y)
    loss.backward()
    optimizer.step()  # 更新参数
    if (epoch + 1) % 50 == 0:
        print(f"Epoch {epoch+1}/{epochs}, Loss: {loss.item():.6f}")

# -------------------------------
# 预测
# -------------------------------
model.eval()  # 切换到评估模式
pred = model(X)  # 得到模型输出（带梯度）
pred = pred.detach()  # 返回一个共享相同数据但不需要梯度的新张量
pred = pred.numpy()  # 转化为 numpy

# -------------------------------
# 计算预测误差
# -------------------------------
y_true = y.detach().numpy() * (data_max - data_min) + data_min
y_pred = pred * (data_max - data_min) + data_min

# 平均绝对误差 MAE
mae = np.mean(np.abs(y_pred - y_true))

# 总绝对误差 SAE（sum of absolute errors）
sae = np.sum(np.abs(y_pred - y_true))

# （可选）均方根误差 RMSE
rmse = np.sqrt(np.mean((y_pred - y_true) ** 2))

print(f"平均绝对误差 (MAE): {mae:.6f}")
print(f"总绝对误差 (SAE): {sae:.6f}")
print(f"均方跟误差 (RMSE): {rmse:.6f}")

import matplotlib

matplotlib.use("TkAgg")  # 或者 "Qt5Agg"，具体取决于环境中装了哪个
import matplotlib.pyplot as plt

plt.plot(range(len(data)), data, label="True")
plt.plot(range(seq_length, len(data)), pred, label="Predicted")
plt.legend()
plt.show()

参考文献¹²³

---

1. https://www.geeksforgeeks.org/deep-learning/deep-learning-introduction-to-long-short-term-memory/ ↩︎

2. https://en.wikipedia.org/wiki/Long_short-term_memory ↩︎

3. https://medium.com/@ottaviocalzone/an-intuitive-explanation-of-lstm-a035eb6ab42c ↩︎