2019/07/26

Least Squares Approximation 最小平方近似

我的技術顧問服務,是提供各種應用於混合模式積體電路設計(mixed mode IC design)的線性代數理論與演算法,以及數位信號處理和數位積體電路設計的相關技術指導。

最近,不同的客戶們的幾個的問題,雖然以各種不同的形式展現,但它們其實都是least square approximation的問題。高斯在二百多年前就知道它的解法。因此我把它寫在底下,給大家做為參考。


一個 over-determined system $$\mathbf{U} \cdot \mathbf{a}= \mathbf{d}$$ 的least squares (LS) approximation是
$$\hat{\mathbf{a}} = (\mathbf{U}^{T} \cdot \mathbf{U})^{-1}\mathbf{U}^T \cdot \mathbf{d}$$例:10元硬幣重x克,5元硬幣重y克,以下四個等式是四組測量數據,請問10元硬幣和5元硬幣重量的『合理的』估計值分別是多少?
$$3x+2y=42$$ $$x+3y=24$$ $$4x+4y=60$$ $$5x+y=54$$ 例題中的量測數據,有四個等式,但只有兩個變數,並沒有一組解可以同時滿足四個等式。我們最多只能找出硬幣重量的估計值$x$與$y$,使得四個等式的誤差盡量小一點,無法讓它完全相等。

四個等式用矩陣寫成 $$ \begin{bmatrix} 3 & 2 \\ 1 & 3 \\ 4 & 4 \\ 5 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 42 \\ 24 \\ 60 \\ 54 \end{bmatrix} $$ 或是一般的形式 $$\mathbf{U} \cdot \mathbf{a} = \mathbf{d}$$ $m \times n$的矩陣在$m > n$時,多數的情況之下,等式往往是無解的。此時只能找到$\mathbf{U} \cdot \mathbf{a}=\hat{\mathbf{d}} \approx \mathbf{d}$。我們能做的只有讓誤差$\mathbf{e}=\mathbf{d}-\hat{\mathbf{d}}$盡量小一點。

在各種可能的近似估計$\mathbf{U} \cdot \mathbf{a}=\hat{\mathbf{d}} \approx \mathbf{d}$當中,其中一個合理的估計方式,是找到『誤差的平方和』最小的那一組當做最佳近似(least squares approximation)。亦即

$min \ \|\mathbf{e}\|^2$
$s.t. \ \ \mathbf{e}=\mathbf{d}-\hat{\mathbf{d}} = \mathbf{d} -\mathbf{U}\cdot \mathbf{a}$


幾何的觀察:
image/svg+xmle d span( U ) d ˆ

$\mathbf{d}=\hat{\mathbf{d}}+\mathbf{e}$

要讓誤差向量的平方和$ \|\mathbf{e}\|^2$最小,則:

$\hat{\mathbf{d}}$是$\mathbf{d}$在$span(\mathbf{U})$上的投影,或者說,$\hat{\mathbf{d}}$是$\mathbf{d}$在子空間$span(\mathbf{U})$的分量。

$\mathbf{e}$是$\mathbf{d}$在$span(\mathbf{U})$投影為0的分量,也就是說,$\mathbf{e}$與子空間$span(\mathbf{U})$垂直,或者寫成$\mathbf{U}^T \cdot \mathbf{e} = 0 $

此時的$\|\mathbf{e}\|^2$最小。

其他的估計,由於誤差向量不是垂直於$span(\mathbf{U})$,其平方和將等於上述$\|\mathbf{e}\|^2$再加上誤差向量在的分量的平方和。

因此,欲求得$\mathbf{U} \cdot \mathbf{a} \approx \mathbf{d}$的 LS approximation,可以同時在$\mathbf{U} \cdot \mathbf{a}$以及$\mathbf{d}$的左方乘上$\mathbf{U}^T$以消除$\mathbf{e}$分量。倘若為full rank則
$$\mathbf{U}^T (\mathbf{U} \mathbf{a}) = \mathbf{U}^T \mathbf{d} = \mathbf{U}^T (\hat{\mathbf{d}}+\mathbf{e}) = \mathbf{U}^T \hat{\mathbf{d}}+\mathbf{U}^T \mathbf{e} = \mathbf{U}^T \hat{\mathbf{d}}$$ 就有唯一解 
$$\hat{\mathbf{a}} = (\mathbf{U}^{T} \cdot \mathbf{U})^{-1}\mathbf{U}^T \cdot \mathbf{d}$$
在例題中$\begin{bmatrix} 3 \\ 1\\ 4 \\ 5 \end{bmatrix}$ 與 $\begin{bmatrix} 2 \\ 3\\ 4 \\ 1\end{bmatrix}$ 的所有線性組合,構成了二維子平面$span(\mathbf{U})$。然而$\mathbf{d}=\begin{bmatrix} 42 \\ 24 \\ 60 \\ 54 \end{bmatrix}$並不落在這個子平面上,因此無解。套用上述的做法:
$$ \begin{bmatrix} 3 & 1 & 4 & 5 \\ 2 & 3 & 4 &1 \end{bmatrix} \begin{bmatrix} 3 & 2 \\ 1 & 3 \\ 4 & 4 \\ 5 & 1 \end{bmatrix} \begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = \begin{bmatrix} 3 & 1 & 4 & 5 \\ 2 & 3 & 4 &1 \end{bmatrix} \begin{bmatrix} 42 \\ 24 \\ 60 \\ 54 \end{bmatrix} $$ 也就是
$$ \begin{bmatrix} 51 & 30 \\ 30 & 30 \end{bmatrix} \begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = \begin{bmatrix} 660 \\ 450 \end{bmatrix}$$ 亦即 $$ \begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = \begin{bmatrix} 10 \\ 5 \end{bmatrix} $$
將這一組解代回原等式得到$$ \hat{\mathbf{d}} = \begin{bmatrix} 3 & 2 \\ 1 & 3 \\ 4 & 4 \\ 5 & 1 \end{bmatrix} \begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = \begin{bmatrix} 40 \\ 25 \\ 60 \\ 55 \end{bmatrix} $$ 以及 $$ e= \begin{bmatrix} 42 \\ 24 \\ 60 \\ 54 \end{bmatrix} - \begin{bmatrix} 40 \\ 25 \\ 60 \\ 55 \end{bmatrix} = \begin{bmatrix} 2 \\ -1 \\ 0 \\ -1 \end{bmatrix} $$ 以及 $$\begin{bmatrix} 3 & 1 & 4 & 5 \\ 2 & 3 & 4 &1 \end{bmatrix} \begin{bmatrix} 3 & 2 \\ 1 & 3 \\ 4 & 4 \\ 5 & 1 \end{bmatrix} \begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = \begin{bmatrix} 3 & 1 & 4 & 5 \\ 2 & 3 & 4 &1 \end{bmatrix} \Bigg( \begin{bmatrix} 42 \\ 24 \\ 60 \\ 54 \end{bmatrix}+ \begin{bmatrix} 2 \\ -1 \\ 0 \\ -1 \end{bmatrix} \Bigg) $$ 可以看到:因為$\mathbf{U}^T \cdot \mathbf{e} = 0$ ,消去了$\mathbf{d}$當中沒有辦法從$\begin{bmatrix} 3 \\ 1\\ 4 \\ 5 \end{bmatrix}$與$\begin{bmatrix} 2 \\ 3\\ 4 \\ 1 \end{bmatrix}$組合而成的分量$\mathbf{e}=\begin{bmatrix} 2 \\ -1 \\ 0 \\ -1 \end{bmatrix}$。因此,兩個變數四個等式的無解的情況,變為兩個變數兩個等式。

1 則留言:

匿名 提到...

古兄能否發篇比較好啃的文章??? 哈

Labels

七條 (1) 大かまど芝 (2) 小狗教我的事 (5) 巴黎印象 (7) 日本大阪 (22) 日本北海道 (40) 日本印象 (27) 日本東京 (40) 日航關西機場飯店 (1) 水井茶堂 (4) 火星在線 (4) 北工房 (1) 台南遠東香格里拉 (4) 台灣旅遊 (46) 本湖月 (4) 甘味處 花 (1) 甲殼公司 (30) 吉崎食堂 (3) 竹北 (10) 牡蠣料理 開 (1) 咖啡 (26) 東京麵通團 (1) 南半球 (5) 叛徒心情日記 (120) 科普 (12) 紀の善 (1) 美瑛選果 (1) 音樂 (15) 香港旅遊 (11) 香港機場富豪酒店 (1) 倫敦旅遊 (7) 宮本 (1) 馬雅烘焙 (3) 強羅花壇 (2) 涵碧樓 (1) 野球 (6) 森の時計 (2) 菡萏咖啡 (3) 新竹喜來登飯店 (1) 溫味 (2) 銀座 古川 (1) 線性代數與信號處理 (1) 閱讀與影劇 (100) 韓國印象 (17) 攝影習作 (40) ANA Gate Tower Hotel Osaka (1) Auberge de Tefutefu (7) Cafe de la Paix (1) COEX (4) Conrad Hotel Singapore (1) Dommy's Dining and Bar (1) Gender (16) Haneda Excel Hotel Tokyu (1) Hilton KL (1) Hotel Monterey Sanno (1) Hotel Nikko Narita (1) Hotel Okura Sapporo (3) HOTEL日航福岡 (1) JAL修行 (36) KAMADO MARUYAMA (3) KLCC Trader Hotel (1) La Sante (2) Mitsui Garden Hotel (1) Moliere (1) Pres Vert (2) Relais Saint Germain (2) Rihga Royal Hotel (4) Sapporo Grand Hotel (3) Sapporo Park Hotel (2) Swissotel Nankai Osaka (1) うさぎや(兔屋) (1) おしどり (2) おばんざい とらや (2) きじ (1) サンパレス球陽館 (2) ねぎ焼 やまもと (2) ホテルコムズ新千歳空港 (2)