最小二乘問題詳解5：非線性最小二乘求解實例

通過擬合非線性函數的實例，詳細講解了非線性最小二乘問題中Gauss-Newton方法的實現過程，并給出了基于Eigen庫的C++代碼示例及收斂分析。

1. 引言

在上一篇文章《最小二乘問題詳解4：非線性最小二乘》中，介紹了非線性最小二乘問題的基本定義、求解思路及其核心算法Gauss-Newton方法，強調通過局部線性化將非線性問題轉化為迭代的線性最小二乘子問題來求解。由于非線性最小二乘問題起來比線性最小二乘復雜多了，這里就通過一個擬合曲線$y = \exp(a x^2 + b x + c)$的實例來加深對非線性最小二乘問題的理解。

2. 雅可比矩陣

最麻煩的還是計算雅可比矩陣。設參數向量為：

\[\boldsymbol{\theta} = \begin{bmatrix} a \\ b \\ c \end{bmatrix} \]

模型函數為：

\[f(x; \boldsymbol{\theta}) = \exp(a x^2 + b x + c) \]

令中間變量：

\[u = a x^2 + b x + c \quad \Rightarrow \quad f = e^u \]

根據定義，雅可比矩陣是模型輸出對參數的偏導形成的矩陣，即：

\[\mathbf{J}_{i,j} = \frac{\partial f(x_i; \boldsymbol{\theta})}{\partial \theta_j} \]

每一行是模型在第$i$個樣本處對參數的梯度。計算每個偏導數的解析形式：

對 $ a $ 求偏導：

\[\frac{\partial f}{\partial a} = \frac{\partial}{\partial a} e^{a x^2 + b x + c} = e^{a x^2 + b x + c} \cdot \frac{\partial}{\partial a}(a x^2 + b x + c) = e^{u} \cdot x^2 = f(x; \boldsymbol{\theta}) \cdot x^2 \]
對 $ b $ 求偏導：

\[\frac{\partial f}{\partial b} = e^{u} \cdot \frac{\partial}{\partial b}(a x^2 + b x + c) = e^{u} \cdot x = f(x; \boldsymbol{\theta}) \cdot x \]
對 $ c $ 求偏導：

\[\frac{\partial f}{\partial c} = e^{u} \cdot \frac{\partial}{\partial c}(a x^2 + b x + c) = e^{u} \cdot 1 = f(x; \boldsymbol{\theta}) \]

因此，對于 $ N $ 個數據點 $ x_1, x_2, \dots, x_N $，雅可比矩陣為：

\[\mathbf{J} = \begin{bmatrix} \frac{\partial f(x_1)}{\partial a} & \frac{\partial f(x_1)}{\partial b} & \frac{\partial f(x_1)}{\partial c} \\ \frac{\partial f(x_2)}{\partial a} & \frac{\partial f(x_2)}{\partial b} & \frac{\partial f(x_2)}{\partial c} \\ \vdots & \vdots & \vdots \\ \frac{\partial f(x_N)}{\partial a} & \frac{\partial f(x_N)}{\partial b} & \frac{\partial f(x_N)}{\partial c} \end{bmatrix}= \begin{bmatrix} f(x_1) \cdot x_1^2 & f(x_1) \cdot x_1 & f(x_1) \\ f(x_2) \cdot x_2^2 & f(x_2) \cdot x_2 & f(x_2) \\ \vdots & \vdots & \vdots \\ f(x_N) \cdot x_N^2 & f(x_N) \cdot x_N & f(x_N) \end{bmatrix} \]

或者寫成：

\[\boxed{ \mathbf{J} = \begin{bmatrix} e^{a x_1^2 + b x_1 + c} \cdot x_1^2 & e^{a x_1^2 + b x_1 + c} \cdot x_1 & e^{a x_1^2 + b x_1 + c} \\ e^{a x_2^2 + b x_2 + c} \cdot x_2^2 & e^{a x_2^2 + b x_2 + c} \cdot x_2 & e^{a x_2^2 + b x_2 + c} \\ \vdots & \vdots & \vdots \\ e^{a x_N^2 + b x_N + c} \cdot x_N^2 & e^{a x_N^2 + b x_N + c} \cdot x_N & e^{a x_N^2 + b x_N + c} \end{bmatrix} } \]

3. 實例

其實要求解非線性最小二乘問題可以使用現成的庫（比如Ceres Solver），不過本文主要為了理解非線性最小二乘的求解過程，尤其是Gauss-Newton方法。因此這個實例還是使用基礎的矩陣庫Eigen來實現，具體代碼如下：

#include <Eigen/Dense>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>

using namespace std;
using namespace Eigen;

// 模型函數: y = exp(a*x^2 + b*x + c)
double model(double x, const Vector3d& theta) {
  double a = theta(0);
  double b = theta(1);
  double c = theta(2);
  double exponent = a * x * x + b * x + c;
  // 防止溢出
  if (exponent > 300) return exp(300);
  if (exponent < -300) return exp(-300);
  return exp(exponent);
}

// 計算 Jacobian 矩陣（數值導數或解析導數）
MatrixXd computeJacobian(const vector<double>& x_data,
                         const vector<double>& y_data, const Vector3d& theta) {
  int N = x_data.size();
  MatrixXd J(N, 3);  // 每行對應一個點，三列：?f/?a, ?f/?b, ?f/?c

  for (int i = 0; i < N; ++i) {
    double x = x_data[i];
    double exponent = theta(0) * x * x + theta(1) * x + theta(2);
    double f = model(x, theta);  // 當前預測值

    // 解析導數（鏈式法則）
    double df_de = f;  // d(exp(u))/du = exp(u)
    double de_da = x * x;
    double de_db = x;
    double de_dc = 1.0;

    J(i, 0) = df_de * de_da;  // ?f/?a
    J(i, 1) = df_de * de_db;  // ?f/?b
    J(i, 2) = df_de * de_dc;  // ?f/?c
  }

  return J;
}

// 計算殘差向量 r_i = y_i - f(x_i; theta)
VectorXd computeResiduals(const vector<double>& x_data,
                          const vector<double>& y_data, const Vector3d& theta) {
  int N = x_data.size();
  VectorXd residuals(N);
  for (int i = 0; i < N; ++i) {
    double pred = model(x_data[i], theta);
    residuals(i) = y_data[i] - pred;
  }
  return residuals;
}

int main() {
  // ========================
  // 1. 設置真實參數
  // ========================
  Vector3d true_params;
  true_params << 0.05, -0.4, 1.0;  // a, b, c
  double a_true = true_params(0), b_true = true_params(1),
         c_true = true_params(2);

  cout << "真實參數: a=" << a_true << ", b=" << b_true << ", c=" << c_true
       << endl;

  // ========================
  // 2. 生成觀測數據（帶高斯噪聲）
  // ========================
  int N = 50;
  vector<double> x_data(N), y_data(N);

  random_device rd;
  mt19937 gen(rd());
  uniform_real_distribution<double> x_dis(-5.0, 5.0);  // x 在 [-5, 5]
  normal_distribution<double> noise_gen(0.0, 0.1);     // 噪聲 ~ N(0, 0.1)

  for (int i = 0; i < N; ++i) {
    x_data[i] = x_dis(gen);
    double y_true = model(x_data[i], true_params);
    y_data[i] = y_true + noise_gen(gen);  // 添加噪聲
  }

  // ========================
  // 3. 初始化參數（隨便猜）
  // ========================
  Vector3d theta = Vector3d::Zero();  // 初始猜測: a=0, b=0, c=0
  cout << "初始猜測: a=" << theta(0) << ", b=" << theta(1) << ", c=" << theta(2)
       << endl;

  // ========================
  // 4. Gauss-Newton 迭代
  // ========================
  int max_iter = 50;
  double tol = 1e-6;
  double prev_cost = 0.0;

  cout << "\n開始 Gauss-Newton 迭代...\n";

  for (int iter = 0; iter < max_iter; ++iter) {
    // 計算殘差 r
    VectorXd r = computeResiduals(x_data, y_data, theta);

    // 計算代價函數 ||r||^2
    double cost = r.squaredNorm();
    cout << "迭代 " << iter << ": 殘差平方和 = " << cost << endl;

    // 收斂判斷
    if (iter > 0 && abs(cost - prev_cost) < tol) {
      cout << "收斂！" << endl;
      break;
    }
    prev_cost = cost;

    // 計算 Jacobian 矩陣
    MatrixXd J = computeJacobian(x_data, y_data, theta);

    // 求解 Gauss-Newton 方程: (J^T J) Δ = J^T r
    // 使用 QR 分解求解（更穩定）
    VectorXd delta = J.colPivHouseholderQr().solve(r);

    // 更新參數: theta = theta + delta
    theta += delta;

    // 可選：限制更新幅度防止發散
    // if (delta.norm() > 1.0) delta *= 1.0 / delta.norm();

    // 檢查參數是否合理（防止溢出）
    if (!theta.allFinite()) {
      cout << "警告：參數發散！停止迭代。" << endl;
      break;
    }
  }

  // ========================
  // 5. 輸出結果
  // ========================
  cout << "\n--- 擬合完成 ---" << endl;
  cout << "估計參數: a=" << theta(0) << ", b=" << theta(1) << ", c=" << theta(2)
       << endl;
  cout << "真實參數: a=" << a_true << ", b=" << b_true << ", c=" << c_true
       << endl;

  // 最終殘差
  VectorXd final_r = computeResiduals(x_data, y_data, theta);
  cout << "最終殘差平方和: " << final_r.squaredNorm() << endl;

  // ========================
  // 6. （可選）計算參數協方差與標準差
  // ========================
  MatrixXd J_final = computeJacobian(x_data, y_data, theta);
  int n = N, p = 3;
  double sigma_squared = final_r.squaredNorm() / (n - p);  // 估計噪聲方差
  MatrixXd cov_theta =
      sigma_squared * (J_final.transpose() * J_final).inverse();

  Vector3d std_error;
  std_error << sqrt(cov_theta(0, 0)), sqrt(cov_theta(1, 1)),
      sqrt(cov_theta(2, 2));

  cout << "\n參數標準差 (近似):" << endl;
  cout << "a: ±" << std_error(0) << endl;
  cout << "b: ±" << std_error(1) << endl;
  cout << "c: ±" << std_error(2) << endl;

  return 0;
}

運行的結果如下：

真實參數: a=0.05, b=-0.4, c=1
初始猜測: a=0, b=0, c=0

開始 Gauss-Newton 迭代...
迭代 0: 殘差平方和 = 15595.7
迭代 1: 殘差平方和 = 9.83388e+37
迭代 2: 殘差平方和 = 1.33087e+37
迭代 3: 殘差平方和 = 1.80114e+36
迭代 4: 殘差平方和 = 2.43757e+35
迭代 5: 殘差平方和 = 3.2989e+34
迭代 6: 殘差平方和 = 4.46457e+33
迭代 7: 殘差平方和 = 6.04214e+32
迭代 8: 殘差平方和 = 8.17715e+31
迭代 9: 殘差平方和 = 1.10666e+31
迭代 10: 殘差平方和 = 1.4977e+30
迭代 11: 殘差平方和 = 2.02691e+29
迭代 12: 殘差平方和 = 2.74313e+28
迭代 13: 殘差平方和 = 3.71242e+27
迭代 14: 殘差平方和 = 5.02421e+26
迭代 15: 殘差平方和 = 6.79954e+25
迭代 16: 殘差平方和 = 9.20217e+24
迭代 17: 殘差平方和 = 1.24538e+24
迭代 18: 殘差平方和 = 1.68544e+23
迭代 19: 殘差平方和 = 2.28099e+22
迭代 20: 殘差平方和 = 3.08698e+21
迭代 21: 殘差平方和 = 4.17778e+20
迭代 22: 殘差平方和 = 5.65401e+19
迭代 23: 殘差平方和 = 7.65187e+18
迭代 24: 殘差平方和 = 1.03557e+18
迭代 25: 殘差平方和 = 1.40149e+17
迭代 26: 殘差平方和 = 1.89671e+16
迭代 27: 殘差平方和 = 2.56691e+15
迭代 28: 殘差平方和 = 3.47393e+14
迭代 29: 殘差平方和 = 4.70143e+13
迭代 30: 殘差平方和 = 6.36261e+12
迭代 31: 殘差平方和 = 8.61052e+11
迭代 32: 殘差平方和 = 1.16541e+11
迭代 33: 殘差平方和 = 1.57677e+10
迭代 34: 殘差平方和 = 2.13228e+09
迭代 35: 殘差平方和 = 2.87975e+08
迭代 36: 殘差平方和 = 3.87593e+07
迭代 37: 殘差平方和 = 5.17234e+06
迭代 38: 殘差平方和 = 677744
迭代 39: 殘差平方和 = 86534.7
迭代 40: 殘差平方和 = 10405.1
迭代 41: 殘差平方和 = 543.309
迭代 42: 殘差平方和 = 4.93922
迭代 43: 殘差平方和 = 0.579241
迭代 44: 殘差平方和 = 0.577034
迭代 45: 殘差平方和 = 0.577034
收斂！

--- 擬合完成 ---
估計參數: a=0.0498221, b=-0.400562, c=1.00071
真實參數: a=0.05, b=-0.4, c=1
最終殘差平方和: 0.577034

參數標準差 (近似):
a: ±0.00041604
b: ±0.00273643
c: ±0.00568421

需要注意的就是一下幾點：

Gauss-Newton方法所使用的迭代逼近過程，在代碼中通常使用一個while/for循環來實現。因此初值選擇和停止條件特別重要，否則容易在循環過程中發散和震蕩，導致不能實現收斂合適的求解值：

for (int iter = 0; iter < max_iter; ++iter) {
    //...

    // 收斂判斷
    if (iter > 0 && abs(cost - prev_cost) < tol) {
    cout << "收斂！" << endl;
        break;
    }

    //...

    // 更新參數: theta = theta + delta
    theta += delta;
}

初值選擇不太容易，需要對求解問題的領域知識有一定的先驗經驗，或者通過使用近似的線性最小二乘問題的解作為初值。這里初值選擇(0,0,0)因為是示例沒有考慮太多的因素，實際應用中具體的問題需要具體分析。
常見的停止條件有三種：
1. 參數變化很小：
\[\boxed{\| \Delta \theta_k \| < \epsilon_1} \]
1. 目標函數$S(\theta) = \| \mathbf{r}(\theta) \|^2$變化很小：
\[\boxed{| S(\theta_{k+1}) - S(\theta_k) | < \epsilon_2} \]
1. 梯度足夠小。
\[\boxed{\| J_k^T \mathbf{r}_k \| < \epsilon_3} \]
目標函數的梯度是$\nabla S(\theta) = 2 J(\theta)^T \mathbf{r}(\theta)$，$J^T \mathbf{r}$ 正是正規方程的右邊項。

關鍵在于雅可比矩陣的計算：

// 計算 Jacobian 矩陣（數值導數或解析導數）
MatrixXd computeJacobian(const vector<double>& x_data,
                        const vector<double>& y_data, const Vector3d& theta) {
    int N = x_data.size();
    MatrixXd J(N, 3);  // 每行對應一個點，三列：?f/?a, ?f/?b, ?f/?c

    for (int i = 0; i < N; ++i) {
        double x = x_data[i];
        double exponent = theta(0) * x * x + theta(1) * x + theta(2);
        double f = model(x, theta);  // 當前預測值

        // 解析導數（鏈式法則）
        double df_de = f;  // d(exp(u))/du = exp(u)
        double de_da = x * x;
        double de_db = x;
        double de_dc = 1.0;

        J(i, 0) = df_de * de_da;  // ?f/?a
        J(i, 1) = df_de * de_db;  // ?f/?b
        J(i, 2) = df_de * de_dc;  // ?f/?c
    }

    return J;
}

可以將這段代碼與推導的雅可比矩陣進行對照。

這一段代碼生成的隨機值不一定總是能夠收斂，因為有可能會遇到雅可比矩陣病態導致更新方向錯誤，初始猜測太差導致發散等問題。Gauss-Newton也理論上易于理解的方法，更加工程化的實踐需要使用Levenberg-Marquardt算法。

posted @ 2025-10-16 20:07 charlee44 閱讀(145) 評論(0) 收藏舉報

刷新頁面返回頂部

人生海海山山而川

分說，不分說，不由分說！

最小二乘問題詳解5：非線性最小二乘求解實例

1. 引言

2. 雅可比矩陣

3. 實例

公告

人生海海 山山而川

分說，不分說，不由分說！

最小二乘問題詳解5：非線性最小二乘求解實例

1. 引言

2. 雅可比矩陣

3. 實例

公告

人生海海山山而川