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$最小。

任何其他的估計,只要誤差向量$\tilde{\mathbf{e}}$不是垂直於$span(\mathbf{U})$,則其平方和$\|\tilde{\mathbf{e}}\|^2$將等於上述$\|\mathbf{e}\|^2$再加上誤差向量$\tilde{\mathbf{e}}$在$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}$分量與$span(\mathbf{U})$垂直,左方乘上$span(\mathbf{U}^T)$之後就消掉$\mathbf{e}$的影響了。如此,
$$\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} $$ 原本的$\mathbf{U}^T \mathbf{U} \mathbf{a} = \mathbf{U}^T \hat{\mathbf{d}}$ 也就是 $$\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 \\ 25 \\ 60 \\ 55 \end{bmatrix}+ \begin{bmatrix} 2 \\ -1 \\ 0 \\ -1 \end{bmatrix} \Bigg) $$ 其中,後面的誤差項 $ e= \begin{bmatrix} 2 \\ -1 \\ 0 \\ -1 \end{bmatrix} $ 垂直於 $\begin{bmatrix} 3 \\ 1 \\ 4 \\ 5 \end{bmatrix} $ 與 $\begin{bmatrix} 2 \\ 3 \\ 4 \\ 1 \end{bmatrix} $ 可以看到:$\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) 科普 (13) 紀の善 (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)