病態方程的例子

從零開始幾何處理：函數擬合 - 知乎 (zhihu.com)

前言

之前學完了 GAMES102，現在打算以 101 類似的形式，以作業為單位整理一下涉及的知識。

不過明顯感覺到 102 涉及的知識更雜，包括函數擬合、離散微分幾何，以及后面的三維重建等知識，后面有機會的話可以單獨寫個總結（挖個大坑）捋一遍。

同樣盡我所能寫的簡單易懂點，如有錯誤歡迎指出……

引入

函數擬合這個詞並不陌生，高中時我們就學過最小二乘法線性擬合二維點，這就是一個函數擬合的例子。

函數擬合，顧名思義，就是有一對數據點，我們要根據這些數據點找到一個函數，來描述數據集。

在實際生活中，許多應用也滿足這一條件，比如三維重建就是用激光雷達等掃描儀采集數據后，再找到一個表面（函數）來擬合這些數據點。

函數擬合通常分為兩大類：插值問題、逼近問題。插值問題限定函數必須經過那些數據點，逼近問題則不一定。

下面來分別介紹一下。

插值問題

首先來定義一下插值問題：已知離散的數據點集，構建（估計、尋找）一個新的數據點集。

換句話說，給出 m 個點  ，求一個函數  作為點集的估計，使得滿足  。我們可以通過函數 f 來獲取新的點集。

函數空間

函數 f 是什么？我們並不知道。它有可能是多項式函數，有可能是三角函數的組合，有可能根本不存在解析式。我們的目標是找到一個比較好的估計。

一般來說，為了更方便地尋找函數 f，我們需要規定一個函數空間，在這個空間內尋找 f。和線性代數里的基向量類似，函數空間的若干基底稱之為基函數，該函數空間內的任意函數都是基函數的線性組合。

一個栗子：多項式

一個最常用的函數空間就是多項式函數空間，基函數為 ​​​

假定函數 f 是一個 n-1 次多項式，那么 f 可以寫為

我們的目的是求得所有的  。換句話說，我們求解的是如下的線性方程組

其中

注意到 G 是一個 ​​ 的矩陣。假設方程之間線性無關，由線性代數知識可知，求解線性方程組需要比較未知數個數 n 和方程個數 m 的大小：

1. 若 m > n，即插值問題中給定點數比未知數要多，則無解。
2. 若 m = n，則有唯一解。
3. 若 m < n，則有無窮解。

值得一提的是，對於多項式函數的插值，當 m = n 時，存在一種快速求解 f 的方式，稱為拉格朗日插值法，如下

即構造了  ，使得  。所以多項式插值可以在  下完成。

任意函數空間

事實上，對於任意函數空間，假設基函數為 ​ ， 可以寫作

類似於上述，插值問題即求解 ​​ 中的 ​，其中  與上述類似，而

解的情況與上述討論類似，不同的是，非多項式函數空間可能並不存在像拉格朗日插值一樣方便的求解方法。

另外，在實際操作中，如果出現 m<n 的情況，即方程個數（給定數據點個數）小於未知量個數時，可以使用某些方法添加約束（添加方程/數據點）使得滿足 m=n。例如取平均等。

多項式函數空間是萬能的嗎？

首先看一個問題：既然我們有了拉格朗日插值法，可以快速插出任意的多項式，那么其他函數空間還有存在的必要嗎？

答案當然是有必要的。事實上，多項式函數空間存在一些問題。

注：本小節中假設 n=m，即插值點數（方程數）等於基函數個數

我們已經知道，插值問題等價於一個線性方程組求解問題，而線性方程組的求解問題又等價於矩陣求逆問題。

然而在實際操作中，逆矩陣可能遭遇數值精度等問題，導致難以求解。例如下圖所示的病態矩陣，對其做很小的擾動便導致了解的劇變

可以把每個方程看做一個超平面，求解方程組變成若干超平面的求交問題，如果超平面接近於平行，那么很小的擾動便會導致解發生巨大的變化。

衡量矩陣的病態程度可以使用矩陣的條件數，即最大與最小的特征值之比，越大表明某一個基底的支配程度越高，則對擾動越敏感。

在插值問題中，系數矩陣只與基函數和給定數據點有關。如果基函數選的不好，那么可能會導致病態問題。

多項式插值的系數矩陣是病態的。因為多項式插值的系數矩陣，即范德蒙德矩陣，的條件數呈指數增加

另外，我們知道任意函數都可以由無窮多的多項式函數任意逼近，但是在實際應用時，我們不可能使用無窮多個基函數。我們希望，在使用有限基函數時，隨着基函數（插值點數）的增加，插值函數逐漸逼近於目標函數。

然而，多項式插值具有強烈的振盪現象，逼近的過程並不穩定。這使得我們對插值函數的可信度有所質疑

綜上所述，多項式函數空間存在其本身的缺點：病態問題、振盪問題。所以我們需要更好的函數空間。

總結

插值問題，即給定若干離散的數據點，構建出新的數據點的問題。

插值問題可以認為是一種函數擬合問題，我們需要在特定的函數空間中尋找，以找到一個函數來描述數據。

在特定函數空間中尋找函數，本質上是一個線性方程組求解問題。

在實際應用中，一個好的函數空間可以讓問題更容易解決。

逼近問題

在實驗學科中，為了得到離散的數據點，常常使用采樣等方式。通常情況下，這不可避免地帶來噪聲、異常值等問題，也就是說不是所有的數據點都 100% 可信。另一方面，如果必須完美擬合所有數據點，這帶來的代價可能是巨大的（例如過高的參數量、計算量）。

舍棄完美擬合所有的數據點，而是尋找一個逼近（近似）函數，就是逼近問題。

更准確的，逼近問題是，給出 m 個點 ​​​​ ，求一個函數 ​​​​ 使其為目標函數的一種最佳逼近（近似）。

類似於插值問題，我們仍然需要一個函數空間。

為了求解此問題，我們需要先定義“什么是最佳逼近（近似）”。事實上這並不唯一。

一個栗子：最小二乘逼近

一個廣為人知的逼近算法是最小二乘法，它定義最佳逼近為：

通過將 f 寫成某一組基函數的線性組合，我們可以通過求導來求解 f。沿用插值問題的符號，即求

則

事實上，在深度學習中，如果將神經網絡看做一個函數空間，那么 loss 函數就是最佳逼近的一種定義。

實驗

1. 實現多項式函數空間的的拉格朗日插值法
2. 實現高斯基函數空間的插值法
3. 實現最小二乘法
4. 實現最小二乘法增加正則項，即嶺回歸

拉格朗日插值法，套公式，兩個 for 即可

void drawPoly(ImDrawList* draw_list, std::vector<pointf2>& points, ImVec2 origin) {
 if (points.size() <= 1) return;

 auto f = [&](float x) {
  int n = points.size(); // a0 + a_1*x + ... + a_{n-1}*x^{n-1}
  float ans = 0;
  for (int i = 0; i < n; i++) {
   float tmp_fz = 1, tmp_fm = 1;
   for (int j = 0; j < n; j++) {
    if (i == j) continue;
    tmp_fz *= points[j][0] - x;
    tmp_fm *= points[j][0] - points[i][0];
   }
   ans += tmp_fz / tmp_fm * points[i][1];
  }
  return ans;
 };

 drawFunc(draw_list, points, origin, f, IM_COL32(255, 0, 255, 255));
}

高斯基函數，直接求解線性方程組

void drawGauss(ImDrawList* draw_list, std::vector<pointf2>& points, ImVec2 origin, float sigma) {
 if (points.size() <= 1) return;
 
 auto g = [&](int i, float x) {
  if (i == 0) return 1.0f;
  float xi = points[i - 1][0];
  return std::exp(-(x - xi) * (x - xi) / (2 * sigma * sigma));
 };

 int n = points.size();

 float nx = 0, ny = 0; // new point
 for (int i = 0; i < n; i++) nx += points[i][0], ny += points[i][1];
 nx /= n, ny /= n;

 // Y=GA
 Eigen::MatrixXf G(n + 1, n + 1);
 for (int i = 0; i <= n; i++) { 
  if (i < n) { // row in [0, n-1], use origin points
   for (int j = 0; j <= n; j++) {
    G(i, j) = g(j, points[i][0]);
   }
  }
  else { // row=n
   for (int j = 0; j <= n; j++) {
    G(i, j) = g(j, nx);
   }
  }
 }

 Eigen::VectorXf Y(n + 1);
 for (int i = 0; i <= n; i++) {
  Y(i) = i < n ? points[i][1] : ny;
 }
 Eigen::VectorXf A = G.colPivHouseholderQr().solve(Y);

 // f(x) = a0 + a_1*g_1(x) + ... + a_n*g_n(x)
 auto f = [&](float x) {
  int n = points.size();  
  float ans = 0;
  for (int i = 0; i <= n; i++) {
   ans += A[i] * g(i, x);
  }
  return ans;
 };

 drawFunc(draw_list, points, origin, f, IM_COL32(255, 255, 0, 255));
}

最小二乘法，對求導的式子直接求解

void drawLeastSquares(ImDrawList* draw_list, std::vector<pointf2>& points, ImVec2 origin, int n) {
 if (points.size() <= 1) return;
 int m = points.size();
 if (n >= m) return;
 //return;
 
 // G^T G A = G^T Y

 Eigen::MatrixXf G(m, n);
 for (int i = 0; i < m; i++) {
  for (int j = 0; j < n; j++) {   
   G(i, j) = j == 0 ? 1.0f : G(i, j - 1) * points[i][0];
  }
 }

 Eigen::VectorXf Y(m);
 for (int i = 0; i < m; ++i)
  Y(i) = points[i][1];

 Eigen::VectorXf A = (G.transpose() * G).inverse() * G.transpose() * Y;

 // f(x) = a0 + a_1*g_1(x) + ... + a_n*g_n(x)
 auto f = [&](float x) {
  float ans = 0, xx = 1.0f;
  for (int i = 0; i < n; i++) {
   ans += A(i) * xx; xx *= x;
  }
  return ans;
 };

 drawFunc(draw_list, points, origin, f, IM_COL32(0, 0, 255, 255));
}

嶺回歸，和上述差不多

void drawRidge(ImDrawList* draw_list, std::vector<pointf2>& points, ImVec2 origin, int n, float lambda) {
 if (points.size() <= 1) return;
 int m = points.size();
 if (n >= m) return;


 // (G^T G + lambda E) A = G^T Y
 Eigen::MatrixXf G(m, n);
 for (int i = 0; i < m; i++) {
  for (int j = 0; j < n; j++) {
   G(i, j) = j == 0 ? 1.0f : G(i, j - 1) * points[i][0];
  }
 }

 Eigen::VectorXf Y(m);
 for (int i = 0; i < m; ++i)
  Y(i) = points[i][1];

 Eigen::MatrixXf E(n, n);
 E.setIdentity();
 Eigen::VectorXf A = (G.transpose() * G + lambda * E).inverse() * G.transpose() * Y;

 // f(x) = a0 + a_1*g_1(x) + ... + a_n*g_n(x)
 auto f = [&](float x) {
  float ans = 0, xx = 1.0f;
  for (int i = 0; i < n; i++) {
   ans += A(i) * xx; xx *= x;
  }
  return ans;
 };

 drawFunc(draw_list, points, origin, f, IM_COL32(0, 255, 0, 255));
}