最小二乘問題詳解3：線性最小二乘實例

本文通過雙線性變換擬合案例，詳細講解了線性最小二乘的求解過程，并深入探討了參數估計的精度評估方法，重點介紹了協方差矩陣的計算與意義。

1. 引言

在上一篇文章《最小二乘問題詳解2：線性最小二乘求解》中筆者詳細介紹了如何求解線性最小二乘問題，一般使用QR分解或者SVD分解法，這里筆者就實現一個具體的案例來驗證一下。

2. 案例

總是舉擬合直線的例子實在太簡單了，這里就使用一個更加復雜一點問題模型：雙線性變換。具體來說，假設存在兩幅地圖需要配置，并且找到了各自地圖上的同名點，可以使用雙線性變換模型來進行快速、初步的校正。也就是說一組同名點滿足如下關系：

\[\begin{cases} x_1 = a_0 x_0 + b_0 y_0 + c_0 x_0 y_0 + d_0 \\ y_1 = a_1 x_0 + b_1 y_0 + c_1 x_0 y_0 + d_1 \end{cases} \]

其中，\((x_0,y_0)\)和\((x_1,y_1)\)是對應的同名點，也就是觀測值。而\(a_0,b_0,c_0,d_0\)和\(a_1,b_1,c_1,d_1\)則是X和Y方向上的雙線性變換系數，也就是待定值。對于觀測值來說，這個函數是非線性的，但是對于待定值來說這個問題模型則是線性的。這也是筆者在《最小二乘問題詳解1：線性最小二乘》中強調的一點：最小二乘問題是線性還是非線性，需要通過待定值來判斷。

要實現這個案例，那么就需要準備一組同名點，\((x_0,y_0)\)可以通過隨機數生成，\((x_1,y_1)\)則可以通過雙線性變換公式加上一點噪聲值來得到。另外，也不需要從頭實現矩陣運算的，一定規模的矩陣運算庫就可以實現線性方程組的求解，比如這里筆者使用的Eigen。具體的案例代碼實現如下：

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

using namespace std;
using namespace Eigen;

int main() {
  // ========================
  // 1. 設置真實參數（我們想要求解的目標）
  // ========================
  Vector4d params_x_true, params_y_true;
  params_x_true << 1.5, -0.3, 0.1, 2.0;    // a0, b0, c0, d0
  params_y_true << 0.4, 1.8, -0.05, -1.0;  // a1, b1, c1, d1

  cout << "真實參數 x: a0=" << params_x_true(0) << ", b0=" << params_x_true(1)
       << ", c0=" << params_x_true(2) << ", d0=" << params_x_true(3) << endl;

  cout << "真實參數 y: a1=" << params_y_true(0) << ", b1=" << params_y_true(1)
       << ", c1=" << params_y_true(2) << ", d1=" << params_y_true(3) << endl;

  // ========================
  // 2. 生成 N 個隨機點 (x0, y0) 并計算 (x1, y1)
  // ========================
  int N = 100;  // 點的數量
  vector<double> x0(N), y0(N), x1(N), y1(N);

  // 隨機數生成器
  random_device rd;  //真隨機數生成器，生成不可預測的隨機種子
  mt19937 gen(rd());                                   //偽隨機數生成器
  uniform_real_distribution<double> dis(-10.0, 10.0);  // 均勻分布，模擬觀測值
  normal_distribution<double> noise(0.0, 0.1);  // 正態分布，模擬噪聲

  for (int i = 0; i < N; ++i) {
    x0[i] = dis(gen);
    y0[i] = dis(gen);

    // 應用雙線性變換
    double x1_true = params_x_true(0) * x0[i] + params_x_true(1) * y0[i] +
                     params_x_true(2) * x0[i] * y0[i] + params_x_true(3);
    double y1_true = params_y_true(0) * x0[i] + params_y_true(1) * y0[i] +
                     params_y_true(2) * x0[i] * y0[i] + params_y_true(3);

    // 添加噪聲
    x1[i] = x1_true + noise(gen);
    y1[i] = y1_true + noise(gen);
  }

  // ========================
  // 3. 構造設計矩陣 A 和觀測向量 b
  // ========================
  // 對于 x 方向：x1 = a0*x0 + b0*y0 + c0*x0*y0 + d0
  MatrixXd A_x(N, 4);  // N 行，4 列（a0, b0, c0, d0）
  VectorXd b_x(N);

  // 對于 y 方向：y1 = a1*x0 + b1*y0 + c1*x0*y0 + d1
  MatrixXd A_y(N, 4);  // N 行，4 列（a1, b1, c1, d1）
  VectorXd b_y(N);

  for (int i = 0; i < N; ++i) {
    A_x(i, 0) = x0[i];          // a0
    A_x(i, 1) = y0[i];          // b0
    A_x(i, 2) = x0[i] * y0[i];  // c0
    A_x(i, 3) = 1.0;            // d0
    b_x(i) = x1[i];

    A_y(i, 0) = x0[i];          // a1
    A_y(i, 1) = y0[i];          // b1
    A_y(i, 2) = x0[i] * y0[i];  // c1
    A_y(i, 3) = 1.0;            // d1
    b_y(i) = y1[i];
  }

  // ========================
  // 4. 使用 Eigen 求解最小二乘
  // ========================
  Vector4d theta_x = A_x.colPivHouseholderQr().solve(b_x);
  Vector4d theta_y = A_y.colPivHouseholderQr().solve(b_y);

  // ========================
  // 5. 輸出結果
  // ========================
  cout << "\n--- 擬合結果 ---" << endl;
  cout << "估計參數 x: a0=" << theta_x(0) << ", b0=" << theta_x(1)
       << ", c0=" << theta_x(2) << ", d0=" << theta_x(3) << endl;

  cout << "估計參數 y: a1=" << theta_y(0) << ", b1=" << theta_y(1)
       << ", c1=" << theta_y(2) << ", d1=" << theta_y(3) << endl;

  // ========================
  // 6. 驗證：計算殘差
  // ========================
  VectorXd residual_x = b_x - A_x * theta_x;
  VectorXd residual_y = b_y - A_y * theta_y;

  cout << "\n殘差平方和 (x): " << residual_x.squaredNorm() << endl;
  cout << "殘差平方和 (y): " << residual_y.squaredNorm() << endl;
  cout << "總殘差平方和: "
       << (residual_x.squaredNorm() + residual_y.squaredNorm()) << endl;

  return 0;
}

有以下幾點需要注意：

隨機生成觀測值數據的實現

// 隨機數生成器
random_device rd;  //真隨機數生成器，生成不可預測的隨機種子
mt19937 gen(rd());                                   //偽隨機數生成器
uniform_real_distribution<double> dis(-10.0, 10.0);  // 均勻分布，模擬觀測值
normal_distribution<double> noise(0.0, 0.1);  // 正態分布，模擬噪聲

中噪聲的隨機值不是隨意生成的，而要符合正態分布（高斯分布）。這是因為偶然誤差（隨機誤差）通常符合正態分布。

這里分成X方向和Y方向求解了兩組線性方程組，因為有兩組不相關的待定系數\((a_0,b_0,c_0,d_0)\)和\((a_1,b_1,c_1,d_1)\)。

本例使用的QR分解法求解的線性最小二乘問題，如果想使用SVD也很簡單，可以將colPivHouseholderQr替換成如下接口：

Vector4d theta_x = A_x.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b_x);
Vector4d theta_y = A_y.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b_y);

最終運行結果如下：

真實參數 x: a0=1.5, b0=-0.3, c0=0.1, d0=2
真實參數 y: a1=0.4, b1=1.8, c1=-0.05, d1=-1

--- 擬合結果 ---
估計參數 x: a0=1.5006, b0=-0.299571, c0=0.100173, d0=1.98956
估計參數 y: a1=0.398688, b1=1.80047, c1=-0.0501635, d1=-0.994459

殘差平方和 (x): 0.745402
殘差平方和 (y): 0.945377
總殘差平方和: 1.69078

3. 精度

3.1 引出

雖然把最小二乘解求出來了，不過筆者更加關心一個問題，那就是求解的精度是多少？如果我只需要求解一個大致的解，那么隨便取四組點求解出來就可以了，反正不能精確求解，得到結果也大差不差——其實這就是最小二乘的意義：我不僅僅求解出來了，還可以明確計算出求解的精度誤差，使得觀測值與求解的符合度始終在這個誤差范圍之內，從而實現科學上從定性到定量的質變。

以這個例子來說，由于是先設置的真實參數再進行仿真求解，那么可以直接計算估計值與真實值的差值來確定誤差。但是在真實場景中，肯定是不知道真實值的，那么就只能估計精度，具體來說就是計算待定參數的協方差矩陣。

不過要說清楚協方差矩陣不是一件容易的事情，因為這個概念的背景知識是《概率論與數理統計》和《誤差理論》和這兩門課程的內容。大概論述一下就是，協方差矩陣描述了待定參數的估值的不確定性有多大，它的對角線元素是每個參數的方差，開根號就是標準差，而標準差越小，說明估計越精確，所以標準差就是“精度”的量化。

3.2 協方差

從頭來講，可能我們大學學習的《概率論與數理統計》都忘光了，但是高中數學學習的概率與統計知識應該還是記得的。假設我們使用尺子測量一根棍子的長度，一共測了5次，得到：

\[x = [10.1, 9.9, 10.0, 10.2, 9.8] \text{ cm} \]

那么平均值（估計值）是：

\[\bar{x} = \frac{1}{5} \sum x_i = 10.0 \text{ cm} \]

這個值就是棍子“真實長度”的最佳估計，那么方差（不確定性）是：

\[s^2 = \frac{1}{n-1} \sum (x_i - \bar{x})^2 = 0.025 \]

標準差（精度）則是：

\[s = \sqrt{0.025} \approx 0.158 \text{ cm} \]

那么可以這樣評估精度：“長度是 \(10.0 \pm 0.16\) cm”。

那么現在進行升級，不是估計一個數，而是在估多個參數，例如本例中的雙線性變換中，需要估計：

\[\theta = [a_0, b_0, c_0, d_0] \]

這就好比需要同時測量四根棍子的長度，并且由于數據的相關性和估計過程的耦合，這些參數并不獨立。為了更好的描述多個參數的估計，就引入了協方差矩陣（Covariance Matrix）：當有多個參數 \(\theta = [\theta_1, \theta_2, \dots, \theta_p]\)，它們的不確定性用一個矩陣表示：

\[\mathrm{Cov}(\theta) = \begin{bmatrix} \mathrm{Var}(\theta_1) & \mathrm{Cov}(\theta_1,\theta_2) & \cdots \\ \mathrm{Cov}(\theta_2,\theta_1) & \mathrm{Var}(\theta_2) & \cdots \\ \vdots & \vdots & \ddots \end{bmatrix} \]

對角線：每個參數的方差 → 表示它自己的不確定性
非對角線：兩個參數之間的協方差 → 表示它們的估計是否相互影響

方差我們還好理解，那么協方差是什么的呢？這個概念也有點抽象，假設有兩個隨機變量 \(X\) 和 \(Y\)（比如身高和體重）：

如果 \(X\) 大時 \(Y\) 也大 → 正協方差
如果 \(X\) 大時 \(Y\) 小 → 負協方差
如果無關 → 協方差接近 0

也就是說，協方差衡量的是兩個變量一起變化的趨勢。如果還要說的更加具體一點，那么可以這樣定義：對于兩個參數 \(\hat{a}_0\) 和 \(\hat_0\)，它們的協方差是：

\[\mathrm{Cov}(\hat{a}_0, \hat_0) = \frac{1}{N-1} \sum_{i=1}^N (\hat{a}_0^{(i)} - \bar{a}_0) (\hat_0^{(i)} - \bar_0) \]

其中：

\(N\)：重復實驗次數
\(\hat{a}_0^{(i)}\)：第 \(i\) 次實驗估計的 \(a_0\)
\(\bar{a}_0\)：所有 \(\hat{a}_0^{(i)}\) 的平均值
\(\bar_0\)：所有 \(\hat_0^{(i)}\) 的平均值

也即是兩個參數的協方差就是兩者誤差乘積的平均，如果誤差經常同號，乘積為正，那么協方差就大于0；如果誤差經常異號，乘積為負，那么協方差就小于0。

3.3 計算

那么，如何得到這個協方差矩陣呢？這里進行推導說明有點復雜，先直接給出結論，以后有機會再進行詳細說明：

\[\mathrm{Cov}(\hat{\theta}) = \sigma^2 (A^T A)^{-1} \]

其中 \(\sigma^2\) 是噪聲方差的估計：

\[\hat{\sigma}^2 = \frac{\|r\|^2}{n - p} \]

\(r\): 殘差
\(n\): 數據點數
\(p\): 參數個數

回到問題，要評價最小二乘解的精度，就可以取協方差矩陣對角線的值，開平方得到每個待定參數的標準差。

具體的代碼實現如下：

// 假設求解了 theta_x
VectorXd residual = b_x - A_x * theta_x;
int n = A_x.rows();  // 數據點數
int p = A_x.cols();  // 參數數 = 4
double sigma_squared = residual.squaredNorm() / (n - p);

// 計算 (A^T A)^{-1}
MatrixXd AtA_inv = (A_x.transpose() * A_x).inverse();

// 協方差矩陣
MatrixXd cov_theta = sigma_squared * AtA_inv;

// 每個參數的標準差（即“精度”的量化）
VectorXd std_error = cov_theta.diagonal().cwiseSqrt();

cout << "參數標準差 (std error):" << endl;
cout << "a0: ±" << std_error(0) << endl;
cout << "b0: ±" << std_error(1) << endl;
cout << "c0: ±" << std_error(2) << endl;
cout << "d0: ±" << std_error(3) << endl;

4. 其他

其實本文關于精度的內容沒有完全展開說清楚，比如精度的指標不僅僅包含標準差，還包含置信區間等概念。不過寫的太多就感覺有點偏題了，在后續的文章中再進一步論述清楚吧。

posted @ 2025-10-05 23:12 charlee44 閱讀(213) 評論(0) 收藏舉報

刷新頁面返回頂部

人生海海山山而川

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

最小二乘問題詳解3：線性最小二乘實例

1. 引言

2. 案例

3. 精度

3.1 引出

3.2 協方差

3.3 計算

4. 其他

公告

人生海海 山山而川

分說，不分說，不由分說！

最小二乘問題詳解3：線性最小二乘實例

1. 引言

2. 案例

3. 精度

3.1 引出

3.2 協方差

3.3 計算

4. 其他

公告

人生海海山山而川

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