曲線配適

本文的閱讀等級:初級

曲線配適 (curve fitting,或稱曲線擬合) 是一種廣泛使用於科學和工程學的模型分析法。假設我們從領域知識得知變數 xy 的關係為一數學函數 y=f(x),但是我們不知道 f 的內部參數值,例如,f(x)=ax^2+bx+c,其中參數 a,b,c 尚待決定。曲線配適是指利用實驗或調查收集得來的一組數據 (x_1,y_1),\ldots,(x_n,y_n) 找出函數的最佳參數,此處最佳的意義是所有數據的誤差平方和最小,之所以選擇誤差平方和作為評斷標準是因為它具有容易計算與便利理論推導的優點。下文我先介紹形式較為簡單的直線配適,接著再推廣至一般曲線配適。

 
考慮以下樣本大小為 n=4 的數據:

\begin{matrix}    \hline    x_i&\vert&1&2&3&5\\    \hline    y_i&\vert&2&3&5&8 \\    \hline    \end{matrix}

設想變數 xy 的關係由一直線描述,即 y=b_0+b_1x。假若 4 個點都同在一直線上,參數 (b_0,b_1) 就是下列線性方程組的解:

\begin{aligned}  b_0+b_1&=2\\  b_0+2b_1&=3\\  b_0+3b_1&=5\\  b_0+5b_1&=8\end{aligned}

事實是這些點並不在同一條直線上 (見下圖),所以此方程式不存在解。自然地,我們希望找到最「靠近」所有資料點的直線。

資料散布與直線配適

令第 i 個資料點的誤差 e_i 為等號兩邊的差額,e_i=y_i-b_0-b_1x_i,我們的目標是求出參數 (b_0,b_1) 使得誤差平方和最小。定義目標函數

\begin{aligned}  E(b_0,b_1)\displaystyle&\overset{\underset{\mathrm{def}}{}}{=}\sum_{i=1}^ne_i^2=\sum_{i=1}^n(y_i-{b_0}-{b_1}x_i)^2\end{aligned}

微分學給出 E 有極小值的產生條件為

\begin{aligned}  \displaystyle\frac{\partial E}{\partial{b_0}}&=-2\sum_{i=1}^n(y_i-{b_0}-{b_1}x_i)=0\\    \frac{\partial E}{\partial{b_1}}&=-2\sum_{i=1}^n(y_i-{b_0}-{b_1}x_i)x_i=0\end{aligned}

將參數 b_0b_1 提出並整理成下列方程式:

\begin{aligned}  nb_0+\displaystyle\left(\sum_{i=1}^nx_i\right)b_1&=\sum_{i=1}^ny_i\\  \left(\sum_{i=1}^nx_i\right)b_0+\left(\sum_{i=1}^nx_i^2\right)b_1&=\sum_{i=1}^nx_iy_i\end{aligned}

或寫為矩陣形式,

\begin{bmatrix}    n&\sum_{i}x_i\\    \sum_{i}x_i&\sum_{i}x_i^2    \end{bmatrix}\begin{bmatrix}    b_0\\  b_1    \end{bmatrix}=\begin{bmatrix}    \sum_{i}y_i\\    \sum_{i}x_iy_i    \end{bmatrix}

注意我們以簡寫符號 \sum_{i} 替代 \sum_{i=1}^n。若 2\times 2 階係數矩陣可逆,由上式即可解得最佳參數,以 \hat{b}_0\hat{b}_1 表示,因此根據給定資料所配適出的直線即為 \hat{y}=\hat{b}_0+\hat{b}_1x

 
直線配適問題也可以用線性代數方法解決。第一個步驟是寫出誤差的矩陣算式

\begin{bmatrix}    e_1\\    \vdots\\    e_n    \end{bmatrix}=\begin{bmatrix}    y_1\\    \vdots\\    y_n    \end{bmatrix}-\begin{bmatrix}    1&x_1\\    \vdots&\vdots\\    1&x_n    \end{bmatrix}\begin{bmatrix}    b_0\\  b_1    \end{bmatrix}

簡記為 \mathbf{e}=\mathbf{y}-X\mathbf{b},誤差平方和即誤差向量長度平方,

E(\mathbf{b})=\Vert\mathbf{y}-X\mathbf{b}\Vert^2

這是一個最小平方問題 (見“從線性變換解釋最小平方近似”)。令 \hat{\mathbf{y}} 代表 n 維向量 \mathbf{y} 至行空間 (column space) C(X) 的正交投影,正交原則指出 \hat{\mathbf{y}} 是行空間 C(X) 中與 \mathbf{y} 距離最小的向量,也就是對於任意 \mathbf{v}\in C(X)\Vert\mathbf{y}-\hat{\mathbf{y}}\Vert\le\Vert\mathbf{y}-\mathbf{v}\Vert,所以最佳解 \hat{\mathbf{b}} 滿足 \hat{\mathbf{y}}=X\hat{\mathbf{b}}。另一種說法是具有最小 (長度) 的誤差向量 \mathbf{y}-X\hat{\mathbf{b}} 正交於行空間 C(X),也就屬於 C(X) 的正交補餘 N(X^T) (見“線性代數基本定理(二)”),故 X^T(\mathbf{y}-X\hat{\mathbf{b}})=\mathbf{0}。我們得到了一個可供計算的條件式,稱為正規方程 (normal equation):

X^TX\hat{\mathbf{b}}=X^T\mathbf{y}

將係數矩陣乘開,

\begin{aligned}  X^TX&=\displaystyle\begin{bmatrix}    1&\cdots&1\\    x_1&\cdots&x_n    \end{bmatrix}\begin{bmatrix}    1&x_1\\    \vdots&\vdots\\    1&x_n    \end{bmatrix}=\begin{bmatrix}    n&\sum_{i}x_i\\    \sum_{i}x_i&\sum_{i}x_i^2    \end{bmatrix}\end{aligned}

再將常數向量乘開,

\begin{aligned}  X^T\mathbf{y}&=\begin{bmatrix}    1&\cdots&1\\    x_1&\cdots&x_n    \end{bmatrix}\begin{bmatrix}    y_1\\    \vdots\\    y_n    \end{bmatrix}=\begin{bmatrix}    \sum_{i}y_i\\    \sum_{i}x_iy_i    \end{bmatrix}\end{aligned}

確認上式和微分法得到的方程式完全相同。

 
回到前面的例子,令

X=\begin{bmatrix}    1&1\\    1&2\\    1&3\\    1&5    \end{bmatrix},~\hat{\mathbf{b}}=\begin{bmatrix}    \hat{b}_0\\    \hat{b}_1    \end{bmatrix},~\mathbf{y}=\begin{bmatrix}    2\\    3\\    5\\    8    \end{bmatrix}

將數值代入正規方程式 X^TX\hat{\mathbf{b}}=X^T\mathbf{y},計算得到

\begin{bmatrix}    4&11\\    11&39    \end{bmatrix}\begin{bmatrix}    \hat{b}_0\\    \hat{b}_1    \end{bmatrix}=\begin{bmatrix}    18\\    63    \end{bmatrix}

解出 \hat{b}_0=9/35\hat{b}_1=54/35,故最佳配適直線是

\displaystyle y=\frac{9}{35}+\frac{54}{35}x

 
接下來我們將直線配適推廣至一般曲線配適。給定 n 個資料點 (x_1,y_1),\ldots,(x_n,y_n),一般曲線配適具有下列形式:

y=b_1\phi_1(x)+\cdots+b_k\phi_k(x)

其中 \phi_1(x),\ldots,\phi_k(x) 是已知函數,b_1,\ldots,b_k 是未知參數。舉兩個例子,第一個例子包含四個函數 \{1,x,x^2,x^3\}y=b_0+b_1x+b_2x^2+b_3x^3;第二個例子包含兩個函數 \{\sin x,\cos x\}y=b_1\sin x+b_2\cos x

 
直線配適方法依然適用於一般曲線配適嗎?是的。一開始我們還是假設所有資料點都在同一曲線上,就有

\begin{aligned}  b_1\phi_1(x_1)+b_2\phi_2(x_1)+\cdots+b_k\phi_k(x_1)&=y_1\\  b_1\phi_1(x_2)+b_2\phi_2(x_2)+\cdots+b_k\phi_k(x_2)&=y_2\\  &\vdots\\  b_1\phi_1(x_n)+b_2\phi_2(x_n)+\cdots+b_k\phi_k(x_n)&=y_n\end{aligned}

上述方程組雖然並不是變數 x 的線性方程,但卻是函數 \phi_j 的線性方程,矩陣表達式為 X\mathbf{b}=\mathbf{y},其中

X=\begin{bmatrix}    \phi_1(x_1)&\phi_2(x_1)&\cdots&\phi_k(x_1)\\    \vdots&\vdots&\ddots&\vdots\\    \phi_1(x_n)&\phi_2(x_n)&\cdots&\phi_k(x_n)    \end{bmatrix},~\mathbf{b}=\begin{bmatrix}    b_1\\  \vdots\\  b_k    \end{bmatrix},~\mathbf{y}=\begin{bmatrix}    y_1\\  \vdots\\  y_n    \end{bmatrix}

同樣地,最佳參數 \hat{\mathbf{b}} 可由正規方程式 X^TX\hat{\mathbf{b}}=X^T\mathbf{y} 解出,其中 X^TXk\times k 階 Gramian 矩陣 (見“特殊矩陣(14):Gramian 矩陣”),

X^TX=\begin{bmatrix}    \sum_{i}\phi_1(x_i)\phi_1(x_i)& \sum_{i}\phi_1(x_i)\phi_2(x_i)&\cdots&\sum_{i}\phi_1(x_i)\phi_k(x_i)\\    \vdots&\vdots&\ddots&\vdots\\    \sum_{i}\phi_k(x_i)\phi_1(x_i)& \sum_{i}\phi_k(x_i)\phi_2(x_i)&\cdots&\sum_{i}\phi_k(x_i)\phi_k(x_i)    \end{bmatrix}

X^T\mathbf{y}k 維向量,

X^T\mathbf{y}=\begin{bmatrix}    \sum_{i}\phi_1(x_i)y_i\\    \vdots\\    \sum_{i}\phi_k(x_i)y_i    \end{bmatrix}

X 的行向量線性獨立,\mathrm{rank}X=k,正規方程式便有唯一解

\hat{\mathbf{b}}=(X^TX)^{-1}X^T\mathbf{y}

如果 X 的 Gram-Schmidt 正交化產生 QR 分解 X=QR,從等價方程式 R\hat{\mathbf{b}}=Q^T\mathbf{y} 很容易解得最佳參數 \hat{\mathbf{b}},相關討論請見“Gram-Schmidt 正交化與 QR 分解”。

Advertisements
This entry was posted in 線性代數專欄, 應用之道 and tagged , , , , . Bookmark the permalink.

4 則回應給 曲線配適

  1. Liang Dai 說:

    周老师,用QR分解来计算那个正规方程的解比直接代入X与b计算有什么优势呢?

  2. ccjou 說:

    https://ccjou.wordpress.com/2010/04/22/gram-schmidt-%E6%AD%A3%E4%BA%A4%E5%8C%96%E8%88%87-qr-%E5%88%86%E8%A7%A3/
    該文最後提到…
    讀者或許納悶:以高斯消去法解線性方程式比 QR 分解簡單容易,何須如此費事呢?從表面上看,似乎如此。然而,在一些特定的問題中,QR 分解的精確度比高斯消去法高;另外,QR 分解有它自己的典型應用——求解最小平方法。

    不過該文沒有提到…
    Xm\times n 階,當 m>>n,正規方程僅使用約 QR 分解的一半運算量,儲存量也較少。

    換句話說,選擇"正確"的方法是很困難的,沒有一定的答案,這件事到現在還是爭論不休。

  3. ming 說:

    老師,方程式的顯示壞掉了…

    • ccjou 說:

      出現黃色盒子的latex path not specified?WordPress的LaTeX rendering錯誤,但我用Chrome,IE11都可以顯示本文。

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s