30、数据处理与神经网络近似方法解析

Mars5

于 2025-08-30 16:02:33 发布

阅读量33

点赞数

CC 4.0 BY-SA版权

分类专栏：数据科学的艺术与实践文章标签：迭代最小二乘法梯度下降法函数预计算

本文链接：https://blog.youkuaiyun.com/mars5/article/details/151202876

数据科学的艺术与实践专栏收录该内容

36 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

数据处理与神经网络近似方法解析

1. 迭代最小二乘法

迭代最小二乘法是一种用于更新试验解的重要方法。曲率矩阵 $B(0)=2[G(0)]^TG$ ，更新试验解 $m(k)$ 的步骤如下：
1. 计算数据残差：$\Delta d = d_{obs} - d_{pre}(m(k))$
2. 计算数据核矩阵 $G(k)$ 的元素：$G_{ij}(k)=\frac{\partial d_{pre}^i}{\partial m^j}\big| {m(k)}$
3. 计算模型参数的更新量：$\Delta m = [G(k)]^T[G(k)]^{-1}[G(k)]^Td {obs}$
4. 更新试验解：$m(k + 1)=m(k)+\Delta m$

对于具有均匀协方差 $\sigma_d^2$ 的不相关数据，后验协方差近似为：$C_m=\sigma_d^2[G(K)]^T[G(K)]^{-1}$ ，其中 $K$ 表示最终迭代。

迭代停止的条件有多种可能：
- 当拟合误差（以后验方差量化）下降到指定值以下，表示拟合效果可接受。
- 当迭代之间的误差不再显著减小，即 $\Delta E = E(m(k)) - E(m(k - 1))$ 下降到指定值以下。
- 当迭代之间的解不再显著变化，即 $[\Delta m(k)^T\Delta m(k)]^{\frac{1}{2}}$ 下降到指定值以下。

此外，可通过线性方程 $Hm = h_{pri}$ 加入先验信息，将 $m(k + 1)=m(k)+\Delta m$ 代入先验信息方程可得：$H\Delta m = h_{pri}-Hm(k)$

2. 拟合未知频率的正弦曲线

假设数据由叠加在背景水平上的单个正弦曲线组成，正弦曲线的振幅、相位、频率以及背景水平均未知，该场景可用非线性方程描述：$d_{pre}^i = m_1d + m_2A\sin(m_4\omega_at_i)+m_3A\cos(m_4\omega_at_i)$

为了对模型参数进行归一化，引入了常数 $d$、$A$ 和 $\omega_a$ 。矩阵 $G(k)$ 的元素为：
- $G_{i1}(k)=d$
- $G_{i2}(k)=A\sin(m_4(k)\omega_at_i)$
- $G_{i3}(k)=A\cos(m_4(k)\omega_at_i)$
- $G_{i4}(k)=m_2(k)A\omega_at_i\cos(m_4(k)\omega_at_i)-m_3(k)A\omega_at_i\sin(m_4(k)A\omega_at_i)$

以黑岩森林温度数据为例，可将 $d$ 设置为数据的样本均值，$A$ 设置为数据范围的一半，$\omega_a$ 设置为 $2\pi/365$ rad/day 。合理的试验解为 $m(0)=[1.0, 0.1, -0.9, 1.0]^T$ 。使用迭代过程（上述迭代最小二乘法步骤）来细化该估计，当模型参数从一次迭代到下一次的变化低于预定值（如 $10^{-4}$ ）时，迭代过程终止。

该方法得出的估计周期为 $365.61\pm0.06$ (95%) 天，但其非线性拟合是否比标准振幅谱密度分析能更好地估计周期存在争议。不过，振幅谱密度的年度峰值中心与非线性分析确定的值非常接近。

3. 梯度下降法

梯度下降法用于寻找函数的最小值。在第 11.6 节讨论的线性化问题中，梯度向量 $b_i(k)=\frac{\partial E}{\partial m^i}\big|_{m(k)}$ 表示误差 $E(m)$ 在试验解 $m(k)$ 处的斜率，它指向误差最快速增加的方向，与减少误差时 $m(k)$ 需扰动的方向相反。

更新试验解 $m(k)$ 的规则可重写为：$\Delta m=-[B(k)]^{-1}b(k)=M(k)v(k)$ ，其中 $M(k)=\frac{b(k)}{\vert\vert b(k)\vert\vert}[B(k)]^{-1}$ ，$v(k)=-\frac{b(k)}{\vert\vert b(k)\vert\vert}$ 。向量 $v(k)$ 是下坡方向，即 $m(k)$ 应扰动的方向以最大程度减少误差；矩阵 $M(k)$ 编码了达到误差 $E(m)$ 最小值所需的扰动大小信息。

梯度向量 $b(k)$ 可从线性化数据核 $G(k)$ 轻松计算：$b(k)=-2[G(k)]^T(d_{obs}-d_{pre})$

计算矩阵 $M(k)$ 较为困难，因为它依赖于曲率矩阵 $B(k)$ 且涉及矩阵求逆。梯度法将 $M(k)$ 近似为 $\alpha I$ ，其中 $\alpha$ 是未知步长。在误差最小化过程开始时对 $\alpha$ 进行初始猜测，并根据规则 $m(k + 1)=m(k)+\alpha v(k)$ 暂时更新试验解。如果误差 $E(k + 1)$ 减小，则迭代过程继续；如果误差增加，说明步长太大，新解超过了最小值。一种更好的拒绝临时结果的准则是 $E(k + 1)>E(k)+c\alpha[v(k)]^Tb(k)$ （即 Armijo 规则），其中 $c$ 通常设置为 $10^{-4}$ 。

与迭代最小二乘法相比，梯度法的主要限制是需要更多的迭代才能达到可接受的解。但对于具有大量数据和模型参数的问题，梯度下降迭代本质上比迭代最小二乘法更简单、更快。

此外，当最小二乘误差 $E(m)=\sum_{i = 1}^{N}e_i^2$ （其中 $e_i = d_{obs}^i - d_{pre}^i(m)$ ）中的个体误差数量非常大时，可采用随机梯度法。该方法在每次迭代中使用不同的随机个体误差子集，只要 $v(k)$ 不太不准确，每个步骤仍会使解向最小值移动。在黑岩森林数据集的情况下，将个体误差数量减少 1000 倍仍可得到可接受的解。

4. 函数预计算与表查找

许多数据分析方法需要多次评估同一函数。如果函数复杂但只有一两个已知范围的参数，预计算函数并将其值存储在表中可能在计算上更高效。只要表查找所需的时间小于评估函数所需的时间，就能减少计算时间。虽然表查找只能提供函数的近似值，但在许多情况下，这种不准确性会被其他误差源（尤其是测量误差）所掩盖。

一维函数 $d(x)$ 的查找表是一个时间序列 $d_k = d(x_k)$ ，其中 $x_k=(k - 1)\Delta x + x_{min}$ ，$1 < k < K$ 。$x_k$ 的范围必须涵盖所有需要估计的 $x$ 值，且采样间隔必须足够小以忠实地表示函数的波动。函数 $d(x)$ 的简单估计是相应表项的值 $d_n$ ，其中 $n = \lfloor\frac{x - x_{min}}{\Delta x}\rfloor + 1$ 。也可通过在两个相邻值 $d_n$ 和 $d_{n + 1}$ 之间进行线性插值来估计 $d(x)$ ，但插值过程会增加计算时间。

以使用网格搜索确定两个模型参数 $(m_1, m_2)$ 为例，数据 $d_i = d(t_i)$ 满足 $d_i=\sin(\pi\sin(\pi x_i)\sin(\pi(1 - x_i^2)))$ ，其中 $x_i = m_1 + m_2t_i$ 。假设 $0 < m_1 < 1$ ，$0 < m_2 < 1$ ，$0 < t < 1$ ，可计算出 $0 < x < 2$ 。查找表为该范围的 $d(x)$ 制表，$d_i$ 的值通过先计算 $x_i = m_1 + m_2t_i$ ，然后在表中查找 $d(x_i)$ 来计算。使用查找表的网格搜索版本的执行时间仅约为直接每次计算 $d_i$ 的版本的一半。

查找表适用于两个或更多变量的函数，但超过二维的表很少有用，因为构建它们所需的时间非常长。

5. 人工神经网络

当函数 $d(x)$ 不准确已知时，可先进行初步的粗略估计，然后随着新信息的出现进行改进。使用查找表表示函数的一个优点是易于更新，但它也具有一定的局限性，如元素宽度 $\Delta x$ 一旦定义就难以更改，且 $d(x)$ 在区间边界处会突然跳跃，可能违反函数平滑性的先验信息。

人工神经网络（简称神经网络）是一种新的近似方法，与查找表有一些相似之处，但更具适应性和光滑性。任意函数 $d(x)$ 可示意性地表示为一行三个框，信息从左到右流动。左框表示输入 $x$ ，右框表示输出 $d(x)$ ，中间框表示函数的运算过程，可能非常复杂，因此被称为隐藏层。表近似将中间框的运算替换为简单的表查找，而表查找过程可视为对箱形或塔形函数的求和。

每个表项的左右边缘可通过将塔形函数视为两个阶跃函数的和来分离：$d_cH(w_Lx + b_L)-d_cH(w_Rx + b_R)$ ，其中 $H(z)$ 是 Heaviside 阶跃函数。权重 $w$ 和偏置 $b$ 参数化了边缘的位置，常数 $d_c$ 指定了塔形函数的高度。在查找表中，相邻塔形函数的边缘相互匹配；而在神经网络中，这种要求被放宽，阶跃函数可以任意重叠。

为了解决平滑性问题，用 Sigmoid 函数 $\sigma(z)=\frac{1}{1 + \exp(-z)}$ （其中 $z = wx + b$ ）代替阶跃函数。Sigmoid 函数在 $z \to \pm\infty$ 时分别渐近于 0 和 1 ，在 $z = 0$ （或 $x_0 = -\frac{b}{w}$ ）处有过渡，且是平滑的。权重 $w$ 控制过渡的尖锐程度，最大斜率 $s = \frac{d\sigma}{dx}=\frac{w}{4}$ 发生在 $x_0$ 处。

6. 神经网络中的信息流

神经网络由 $L$ 列（或层）的框组成，第 $k$ 层包含 $N(k)$ 个框，这些框被称为神经元。信息从左到右沿着连接相邻层框的线段流动。第 $k$ 层的神经元 $i$ 有一个偏置 $b_i(k)$ ，神经元 $i$ 与第 $k - 1$ 层的神经元 $j$ 之间的连接有一个权重 $w_{ij}(k)$ ，权重为 0 表示没有连接。神经元 $i$ 还有一个输入值 $z_i(k)$ 和一个输出值（或活动） $a_i(k)$ 。第一层神经元的活动表示函数的输入，第 $L$ 层神经元的活动表示函数的输出。

信息在层与层之间传播的规则如下：
- 第 $k$ 层神经元 $i$ 的输入：$z_i(k)=\sum_{j = 1}^{N(k - 1)}w_{ij}(k)a_j(k - 1)+b_i(k)$
- 第 $k$ 层神经元 $i$ 的活动：$a_i(k)=\begin{cases}\sigma(z_i(k)) & k < L\z_i(k) & k = L\end{cases}$

以下是 MATLAB 和 Python 中实现神经网络信息传播的代码：

function [ a ] = eda_net_eval(N,w,b,a)
% evaluates network with properites N, w, b, and
% initial activities a, returns final activities, a
Nmax = max(N);
Lmax = length(N);
for i=[2:Lmax]
    z(1:N(i),i)=w(1:N(i),1:N(i-1),i)*a(1:N(i-1),i-1)+b(1:N(i),i);
    if i < Lmax
        a(1:N(i),i) = 1.0./(1.0+exp(-z(1:N(i),i)));
    else
        a(1:N(i),i) = z(1:N(i),i);
    end
end
end

import numpy as np

def eda_net_eval(N,w,b,a):
    # eda_eval_net evaluates a network with properties N, w and b
    # the "input" activities a[0:N[0,1],0] must be set
    # the output are the final activities and the final z0s
    Nmax = np.max(N);
    Lmax, tmp = np.shape(N);
    z = np.zeros( (Nmax, Lmax) );
    for i in range(1,Lmax):
        M = w[0:N[i,0],0:N[i-1,0],i];
        v = a[0:N[i-1,0],i-1];
        u = b[0:N[i,0],i];
        zz = np.zeros((N[i,0],1));
        zz[:,0] = np.matmul( M, v ) + u;
        z[0:N[i,0],i] = zz[:,0];
        if i < (Lmax-1):
            a[0:N[i,0],i] = eda_sigma(1.0,zz).ravel();
        else:
            a[0:N[i,0],i] = z[0:N[i,0],i];
    return a, z;

平滑阶跃函数和塔形函数的神经网络示例见图 11.12 。由于输出层省略了 Sigmoid 函数，该层的活动只是其输入的加权和，这使得网络易于合并。例如，几个单塔网络可以合并以产生塔形函数的和。以下是初始化网络为塔形函数和的 MATLAB 和 Python 代码：

function [N,w,b] = eda_net_1dtower(Nc,slope,Xc,Dx,h,N,w,b)
% eda_net_linear
% defines a network that computes a set of 1D towers
% input Nc: number of towers
%     slope: slope of all the towers
%     Xc: centers of the towers
%     Dx: width of the towers
%     h: heights of the towaers
%     N: layer array from eda_net_init()
%     w: weights array from eda_net_init()
%     b: bias array from from eda_net_init()
% output: updated N, w, b arrays
...

def eda_net_1dtower( Nc, slope, Xc, Dx, h, N, w, b ):
    # eda_net_linear
    # defines a network that computes a set of 1D towers
    # input Nc: number of towers
    #     slope: slope of all the towers
    #     Xc: centers of the towers
    #     Dx: width of the towers
    #     h: heights of the towaers
    #     N: layer array from eda_net_init()
    #     w: weights array from eda_net_init()
    #     b: bias array from from eda_net_init()
    # output: updated N, w, b arrays
    ...

二维塔形函数（如输入为 $(x, y)$ ）可通过合并两个一维塔形网络（$x$ 方向和 $y$ 方向的塔形网络）创建。但它们的和不是完全的二维塔形函数，而是由两个相交的脊组成，需要一个第四层来删除这些脊并形成塔形函数。以下是初始化网络为二维塔形函数和的 MATLAB 和 Python 代码：

function [N,w,b] = ...
eda_net_2dtower(Nc,slope,Xc,Yc,Dx,Dy,h,N,w,b)
% eda_net_linear
% defines a network that superimposes several towers
% input:
%     Nc: number of towers
%     slope: slope of all the towers
%     Xc, Yc: vector of the centers of the towers in x and y
%     Dx, Dy: vector of the widths of the towers in x and y
%     h: vector of heights of the towers
%     N: layer array from eda_net_init()
%     w: weights array from eda_net_init()
%     b: bias array from from eda_net_init()
% output: updated N, w, b arrays

def eda_net_2dtower( Nc, slope, Xc, Yc, Dx, Dy, h, N, w, b ):
    # defines a network that superimposes several towers
    # input:
    #     Nc: number of towers
    #     slope: slope of all the towers
    #     Xc, Yc: vector of the centers of the towers in x and y
    #     Dx, Dy: vector of the widths of the towers in x and y
    #     h: vector of heights of the towers
    #     N: layer array from eda_net_init()
    #     w: weights array from eda_net_init()
    #     b: bias array from from eda_net_init()
    # output: updated N, w, b arrays
    ...

通过叠加足够多的塔形函数，神经网络能够以任意精度近似任何连续函数。但以这种方式构建的网络可能不是能够充分近似该函数的最小网络。后续可开发技术来训练任意网络以产生其所能达到的函数的最佳近似，这是一种简单的机器学习形式。

总结

本文介绍了多种数据处理和近似方法，包括迭代最小二乘法、梯度下降法、函数预计算与表查找以及人工神经网络。这些方法在不同场景下各有优劣，可根据具体问题选择合适的方法。迭代最小二乘法适用于需要精确拟合的问题，但计算量较大；梯度下降法计算简单，适用于大规模数据问题；函数预计算与表查找可提高计算效率；人工神经网络则具有强大的函数近似能力。

流程图：迭代最小二乘法流程

graph TD;
    A[开始] --> B[初始化试验解m(0)];
    B --> C[计算数据残差Δd];
    C --> D[计算数据核矩阵G(k)];
    D --> E[计算模型参数更新量Δm];
    E --> F[更新试验解m(k+1)];
    F --> G{是否满足停止条件};
    G -- 是 --> H[结束];
    G -- 否 --> C;

表格：不同方法对比

方法	优点	缺点	适用场景
迭代最小二乘法	拟合精度高	计算量大，需要矩阵求逆	数据量较小，对拟合精度要求高的问题
梯度下降法	计算简单，无需矩阵求逆	需要更多迭代，收敛速度慢	大规模数据和模型参数的问题
函数预计算与表查找	计算效率高	只能提供近似值	多次评估同一复杂函数的问题
人工神经网络	强大的函数近似能力，适应性强	训练过程复杂，计算资源需求大	函数形式复杂且未知的问题

数据处理与神经网络近似方法解析（续）

7. 人工神经网络的训练与应用

人工神经网络虽然能够近似任何连续函数，但构建的网络可能并非是能充分近似函数的最小网络。因此，需要对网络进行训练，使其产生函数的最佳近似，这也是一种简单的机器学习形式。

在训练神经网络时，一般会采用误差最小化的策略。通过不断调整网络中的权重 $w$ 和偏置 $b$ ，使得网络输出与真实数据之间的误差尽可能小。常见的误差度量方式有均方误差（MSE）等，即 $E = \frac{1}{n}\sum_{i = 1}^{n}(y_{i}^{pred} - y_{i}^{true})^2$ ，其中 $y_{i}^{pred}$ 是网络的预测输出，$y_{i}^{true}$ 是真实数据。

训练的具体步骤如下：
1. 初始化参数 ：随机初始化网络中的权重 $w$ 和偏置 $b$ 。
2. 前向传播 ：根据输入数据和当前的权重、偏置，按照信息传播规则计算网络的输出。
3. 计算误差 ：使用选定的误差度量方式计算网络输出与真实数据之间的误差。
4. 反向传播 ：根据误差计算权重和偏置的梯度，通过链式法则将误差从输出层反向传播到输入层。
5. 更新参数 ：根据计算得到的梯度，使用优化算法（如梯度下降法）更新权重和偏置。
6. 重复步骤 2 - 5 ：不断迭代，直到误差满足停止条件（如误差小于某个阈值，或者迭代次数达到上限）。

以下是一个简单的 Python 示例，展示如何使用梯度下降法训练一个简单的神经网络：

import numpy as np

# 定义 Sigmoid 函数
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

# 定义 Sigmoid 函数的导数
def sigmoid_derivative(x):
    return x * (1 - x)

# 输入数据
X = np.array([[0, 0, 1],
              [0, 1, 1],
              [1, 0, 1],
              [1, 1, 1]])

# 真实输出数据
y = np.array([[0], [1], [1], [0]])

# 随机初始化权重
np.random.seed(1)
weights0 = 2 * np.random.random((3, 4)) - 1
weights1 = 2 * np.random.random((4, 1)) - 1

# 训练参数
epochs = 60000
learning_rate = 0.1

# 训练过程
for epoch in range(epochs):
    # 前向传播
    layer0 = X
    layer1 = sigmoid(np.dot(layer0, weights0))
    layer2 = sigmoid(np.dot(layer1, weights1))

    # 计算误差
    layer2_error = y - layer2

    # 反向传播
    layer2_delta = layer2_error * sigmoid_derivative(layer2)
    layer1_error = layer2_delta.dot(weights1.T)
    layer1_delta = layer1_error * sigmoid_derivative(layer1)

    # 更新权重
    weights1 += layer1.T.dot(layer2_delta) * learning_rate
    weights0 += layer0.T.dot(layer1_delta) * learning_rate

    if epoch % 10000 == 0:
        print(f'Epoch {epoch}: Error = {np.mean(np.abs(layer2_error))}')

print("Training finished!")

8. 不同方法在实际问题中的应用案例

为了更直观地了解各种方法的应用效果，下面以一个实际问题为例，对比不同方法的表现。

假设我们要对一个复杂的时间序列数据进行拟合和预测，该数据具有一定的周期性和噪声。

迭代最小二乘法 ：
- 操作步骤 ：
  1. 对数据进行预处理，去除噪声和异常值。
  2. 选择合适的模型函数，例如多项式函数或三角函数的组合。
  3. 初始化试验解 $m(0)$ 。
  4. 按照迭代最小二乘法的步骤，不断更新试验解，直到满足停止条件。
- 应用效果 ：能够较好地拟合数据的主要特征，但对于噪声的处理能力较弱，可能会出现过拟合的情况。
梯度下降法 ：
- 操作步骤 ：
  1. 定义误差函数，如均方误差。
  2. 初始化权重和偏置。
  3. 计算梯度向量。
  4. 根据梯度向量更新权重和偏置，不断迭代。
- 应用效果 ：计算简单，能够在一定程度上处理噪声，但收敛速度较慢，需要更多的迭代次数。
函数预计算与表查找 ：
- 操作步骤 ：
  1. 分析数据的特征，确定需要预计算的函数。
  2. 确定函数的参数范围，构建查找表。
  3. 在计算过程中，通过查找表获取函数值。
- 应用效果 ：显著提高了计算效率，但由于是近似计算，可能会引入一定的误差。
人工神经网络 ：
- 操作步骤 ：
  1. 设计网络结构，包括层数、每层的神经元数量。
  2. 初始化权重和偏置。
  3. 使用训练数据对网络进行训练，调整权重和偏置。
  4. 使用训练好的网络进行预测。
- 应用效果 ：能够很好地捕捉数据的复杂特征和模式，但训练过程复杂，需要大量的计算资源和时间。

流程图：神经网络训练流程

graph TD;
    A[开始] --> B[初始化权重和偏置];
    B --> C[前向传播计算输出];
    C --> D[计算误差];
    D --> E[反向传播计算梯度];
    E --> F[更新权重和偏置];
    F --> G{是否满足停止条件};
    G -- 是 --> H[结束];
    G -- 否 --> C;

表格：不同方法在实际问题中的表现对比

方法	拟合效果	计算效率	抗噪声能力	复杂度
迭代最小二乘法	较好，可能过拟合	低	弱	高
梯度下降法	一般，收敛慢	中	中	中
函数预计算与表查找	有一定误差	高	中	低
人工神经网络	好，能捕捉复杂模式	低	强	高

总结

本文详细介绍了多种数据处理和近似方法，包括迭代最小二乘法、梯度下降法、函数预计算与表查找以及人工神经网络。通过对这些方法的原理、操作步骤和应用案例的分析，我们可以看到它们在不同场景下各有优劣。在实际应用中，需要根据具体问题的特点，如数据规模、复杂度、噪声水平等，选择合适的方法。同时，也可以考虑将多种方法结合使用，以达到更好的效果。未来，随着数据科学和机器学习的不断发展，这些方法也将不断改进和完善，为解决更复杂的问题提供更强大的工具。