約束最小平方問題

本文的閱讀等級：中級

令 $A$ 為一個 $m\times n$ 階實矩陣， $\mathbf{b}\in\mathbb{R}^m$ 。如果線性方程 $A\mathbf{x}=\mathbf{b}$ 是不一致的 (即不存在解)，實務的作法是將線性方程問題改為最小平方近似問題：

$\displaystyle \min_{\mathbf{x}}\Vert A\mathbf{x}-\mathbf{b}\Vert^2$ ，

其中 $\Vert A\mathbf{x}-\mathbf{b}\Vert$ 是2-範數 (見“向量範數”)，即 $A\mathbf{x}$ 與 $\mathbf{b}$ 的歐幾里得距離。根據正交原則，最小平方解 $\hat{\mathbf{x}}$ 滿足正規方程 (normal equation) $A^TA\hat{\mathbf{x}}=A^T\mathbf{b}$ (見“從線性變換解釋最小平方近似”)。若 $\hbox{rank}A=n$ ，也就是說 $A$ 的行向量 (column vector) 構成一個線性獨立集合，則存在唯一的最小平方解

$\displaystyle \hat{\mathbf{x}}=(A^TA)^{-1}A^T\mathbf{b}$ 。

如果最小平方解必須滿足某些束縛條件，則稱為約束最小平方問題 (constrained least-squares problem)。本文討論兩種常出現在多種應用場合的約束形式。線性約束最小平方問題是指束縛條件為線性方程^[1]：

$\displaystyle \min_{C\mathbf{x}=\mathbf{d}}\Vert A\mathbf{x}-\mathbf{b}\Vert^2$ ，

其中 $C$ 是一個 $p\times n$ 階實矩陣， $\mathbf{d}\in\mathbb{R}^p$ 。正則 (regularized) 最小平方問題限制未知向量的長度必須固定：

$\displaystyle \min_{\Vert\mathbf{x}\Vert=d}\Vert A\mathbf{x}-\mathbf{b}\Vert^2$ 。

線性約束最小平方問題

考慮下面的問題：

$\displaystyle \min_{C\mathbf{x}=\mathbf{d}}\Vert A\mathbf{x}-\mathbf{b}\Vert^2$ ，

其中 $A$ 是一個 $m\times n$ 階實矩陣， $C$ 是一個 $p\times n$ 階實矩陣， $\mathbf{b}\in\mathbb{R}^m$ ， $\mathbf{d}\in\mathbb{R}^p$ 。我們要求限制條件數不大於未知向量的維數，即 $p\le n$ 。極小範數問題 (minimum norm problem) 是一種特殊的線性約束最小平方問題：

$\displaystyle \min_{C\mathbf{x}=\mathbf{d}}\Vert \mathbf{x}\Vert^2$ 。

若 $\hbox{rank}C=p$ ，則存在唯一的極小範數解 (見“極小範數解”)

$\displaystyle \hat{\mathbf{x}}=C^T(CC^T)^{-1}\mathbf{d}$ 。

底下介紹一般的線性約束最小平方問題的解法。

Lagrange 乘數法是最常被使用的約束優化問題解法 (見“Lagrange 乘數法”)。定義 Lagrangian 函數

$L(\mathbf{x},\boldsymbol{\alpha})=\Vert A\mathbf{x}-\mathbf{b}\Vert^2+2\boldsymbol{\alpha}^T(C\mathbf{x}-\mathbf{d})$ ，

其中 $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_p)^T$ 是 Lagrange 乘數組成的向量。展開 Lagrangian 函數，

$\displaystyle\begin{aligned} L(\mathbf{x},\boldsymbol{\alpha})&=(A\mathbf{x}-\mathbf{b})^T(A\mathbf{x}-\mathbf{b})+2\boldsymbol{\alpha}^T(C\mathbf{x}-\mathbf{d})\\ &=\mathbf{x}^TA^TA\mathbf{x}-2\mathbf{x}^TA^T\mathbf{b}+\mathbf{b}^T\mathbf{b}+2\boldsymbol{\alpha}^T(C\mathbf{x}-\mathbf{d}),\end{aligned}$

接著計算偏導數 (見“矩陣導數”)，

$\displaystyle \begin{aligned} \frac{\partial L}{\partial \mathbf{x}}&=2A^TA\mathbf{x}-2A^T\mathbf{b}+2C^T\boldsymbol{\alpha}\\ \frac{\partial L}{\partial\boldsymbol{\alpha}}&=2(C\mathbf{x}-\mathbf{d}). \end{aligned}$

令上面兩式等於零可得優化條件式：

$\begin{bmatrix} A^TA&C^T\\ C&0 \end{bmatrix}\begin{bmatrix} \mathbf{x}\\ \boldsymbol{\alpha} \end{bmatrix}=\begin{bmatrix}A^T\mathbf{b}\\ \mathbf{d}\end{bmatrix}$ 。

我們可以證明 $(n+p)\times (n+p)$ 階係數矩陣 $\begin{bmatrix} A^TA&C^T\\ C&0 \end{bmatrix}$ 是可逆的一個充分條件為 $\hbox{rank}A=n$ 且 $\hbox{rank}C=p$ 。考慮齊次方程

$\displaystyle \begin{bmatrix} A^TA&C^T\\ C&0 \end{bmatrix}\begin{bmatrix} \mathbf{x}\\ \mathbf{y} \end{bmatrix}=\mathbf{0}$ ，

即 $A^TA\mathbf{x}+C^T\mathbf{y}=\mathbf{0}$ 且 $C\mathbf{x}=\mathbf{0}$ 。第一式左乘 $\mathbf{x}^T$ 並使用第二式，可得

$0=\mathbf{x}^TA^TA\mathbf{x}+\mathbf{x}^TC^T\mathbf{y}=(A\mathbf{x})^T(A\mathbf{x})+(C\mathbf{x})^T\mathbf{y}=\Vert A\mathbf{x}\Vert^2$ ，

故 $A\mathbf{x}=\mathbf{0}$ 並有 $C^T\mathbf{y}=-A^TA\mathbf{x}=\mathbf{0}$ 。若 $\hbox{rank}A=n$ 和 $\hbox{rank}C=\hbox{rank}C^T=p$ ，推論 $\mathbf{x}=\mathbf{0}$ 和 $\mathbf{y}=\mathbf{0}$ ，即證明 $\begin{bmatrix} A^TA&C^T\\ C&0 \end{bmatrix}$ 是一個可逆矩陣。在此情況下，解開優化條件式可得唯一的約束最小平方解 $\hat{\mathbf{x}}$ ，不過這個方法的主要缺點在於 $A^TA$ 容易引入誤差。

除了以 Lagrangian 函數推導，另一個作法使用 QR 分解 (見“線代膠囊──QR 分解”)。為清楚說明，我們假設 $\hbox{rank}A=n$ 且 $\hbox{rank}C=p$ ， $n\ge p$ 。設 $n\times p$ 階矩陣 $C^T$ 的 QR 分解為

$C^T=Q\begin{bmatrix} R\\ 0 \end{bmatrix}$ ，

其中 $Q$ 是一個 $n\times n$ 階實正交矩陣， $Q^T=Q^{-1}$ ， $R$ 是 $p\times p$ 階上三角可逆矩陣。令

$AQ=\begin{bmatrix} A_1&A_2\end{bmatrix},~~Q^T\mathbf{x}=\begin{bmatrix} \mathbf{x}_1\\ \mathbf{x}_2\end{bmatrix}$ ，

其中 $A_1$ 是 $m\times p$ 階， $A_2$ 是 $m\times (n-p)$ 階， $\mathbf{x}_1\in\mathbb{R}^p$ 且 $\mathbf{x}_2\in\mathbb{R}^{n-p}$ 。寫出

$A\mathbf{x}-\mathbf{b}=AQQ^T\mathbf{x}-\mathbf{b}=\begin{bmatrix} A_1&A_2 \end{bmatrix}\begin{bmatrix} \mathbf{x}_1\\ \mathbf{x}_2 \end{bmatrix}-\mathbf{b}=A_1\mathbf{x}_1+A_2\mathbf{x}_2-\mathbf{b}$

以及

$C\mathbf{x}=\begin{bmatrix} R^T&0 \end{bmatrix}Q^T\mathbf{x}=\begin{bmatrix} R^T&0 \end{bmatrix}\begin{bmatrix} \mathbf{x}_1\\ \mathbf{x}_2 \end{bmatrix}=R^T\mathbf{x}_1$ 。

原線性約束最小平方問題等價於

$\displaystyle \min_{R^T\mathbf{x}_1=\mathbf{d}}\left\Vert A_1\mathbf{x}_1+A_2\mathbf{x}_2-\mathbf{b}\right\Vert^2$ 。

因此， $\mathbf{x}_1$ 由一致的線性方程 $R^T\mathbf{x}_1=\mathbf{d}$ 決定，而 $\mathbf{x}_2$ 則為下列無約束最小平方解：

$\displaystyle \min_{\mathbf{x}_2}\left\Vert A_2\mathbf{x}_2-(\mathbf{b}-A_1\mathbf{x}_1)\right\Vert^2$ 。

合併以上結果可得 $\hat{\mathbf{x}}=Q\begin{bmatrix} \mathbf{x}_1\\ \mathbf{x}_2 \end{bmatrix}$ 。

附帶一提，約束優化問題還可用無約束優化問題近似。考慮

$\displaystyle \min_{\mathbf{x}}\Vert A\mathbf{x}-\mathbf{b}\Vert^2+\beta^2\Vert C\mathbf{x}-\mathbf{d}\Vert^2$ ，

其中 $\Vert C\mathbf{x}-\mathbf{d}\Vert^2$ 稱為懲罰項， $\beta$ 是一個很大的數，稱為懲罰係數。上式可改寫為標準的最小平方問題

$\displaystyle \min_{\mathbf{x}}\left\Vert \begin{bmatrix} A\\ \beta C \end{bmatrix}\mathbf{x}-\begin{bmatrix} \mathbf{b}\\ \beta\mathbf{d}\end{bmatrix}\right\Vert^2$ 。

當 $\beta\to\infty$ ， $\Vert C\mathbf{x}-\mathbf{b}\Vert\to 0$ ，故最小平方解趨於受限的最小平方解。看一個例子^[2]，線性約束最小平方問題

$\displaystyle \min_{x_1=x_2}\left\Vert\begin{bmatrix} 1&2\\ 3&4\\ 5&6\end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}-\begin{bmatrix} 7\\ 1\\ 3\end{bmatrix}\right\Vert^2$

的最佳解為 $\hat{\mathbf{x}}=(0.3407821,0.3407821)^T$ 。如果以下列最小平方問題近似

$\displaystyle \min_{\mathbf{x}}\left\Vert\left[\!\!\begin{array}{rr} 1&2\\ 3&4\\ 5&6\\ 1000&-1000 \end{array}\!\!\right]\begin{bmatrix} x_1\\ x_2\end{bmatrix}-\begin{bmatrix}7\\1\\3\\0\end{bmatrix}\right\Vert^2$ ，

最佳解則為 $\hat{\mathbf{x}}'=(0.3407810,0.3407829)^T$ 。

正則最小平方問題

給定一個 $m\times n$ 階實矩陣 $A$ ， $\hbox{rank}A=n$ ，以及 $\mathbf{b}\in\mathbb{R}^m$ ，考慮這個問題：

$\displaystyle \min_{\Vert\mathbf{x}\Vert=d}\Vert A\mathbf{x}-\mathbf{b}\Vert^2$ 。

下面分開兩種情況討論： $\mathbf{b}=\mathbf{0}$ 和 $\mathbf{b}\neq\mathbf{0}$ 。若 $\mathbf{b}=\mathbf{0}$ ，在不失一般性的原則下，正則最小平方問題可表示為

$\displaystyle \min_{\Vert\mathbf{x}\Vert=1}\Vert A\mathbf{x}\Vert^2$ 。

設 $A$ 的奇異值分解為 $A=U\Sigma V^T$ ，其中 $U$ 是 $m\times m$ 階實正交矩陣， $U^T=U^{-1}$ ， $V$ 是 $n\times n$ 階實正交矩陣， $V^T=V^{-1}$ ， $\Sigma=\begin{bmatrix}D\\0\end{bmatrix}$ 是 $m\times n$ 階矩陣， $D=\hbox{diag}(\sigma_1,\ldots,\sigma_n)$ ， $\sigma_1\ge\cdots\ge\sigma_n>0$ (見“奇異值分解的幾何意義”)。令 $\mathbf{y}=V^T\mathbf{x}$ 。因此， $\Vert\mathbf{y}\Vert^2=\mathbf{x}^TVV^T\mathbf{x}=\Vert\mathbf{x}\Vert^2=1$ 且 $\Vert A\mathbf{x}\Vert^2=\Vert U\Sigma V^T\mathbf{x}\Vert^2=\Vert \Sigma \mathbf{y}\Vert^2$ 。我們得到約化的等價問題：

$\displaystyle \min_{\Vert\mathbf{y}\Vert=1}\Vert \Sigma\mathbf{y}\Vert^2=\min_{y_1^2+\cdots+y_n^2=1}(\sigma_1y_1)^2+\cdots+(\sigma_ny_n)^2$ 。

明顯地， $\mathbf{y}=\mathbf{e}_n=(0,\ldots,0,1)^T$ 最小化目標函數。寫出正交矩陣 $V$ 的行向量表示式 $V=\begin{bmatrix} \mathbf{v}_1&\cdots&\mathbf{v}_n \end{bmatrix}$ ，正則最小平方解為 $V$ 的最末一個行向量， $\hat{\mathbf{x}}=V\mathbf{e}_n=\mathbf{v}_n$ 。在此情形下，最小的目標函數值為 $\Vert A\hat{\mathbf{x}}\Vert^2=\Vert\Sigma \mathbf{e}_n\Vert^2=\sigma_n^2$ 。

假設 $\mathbf{b}\neq\mathbf{0}$ 。定義 Lagrangian 函數

$L(\mathbf{x},\lambda)=\Vert A\mathbf{x}-\mathbf{b}\Vert^2+\lambda(\Vert\mathbf{x}\Vert^2-d^2)$ ，

其中 $\lambda$ 是 Lagrange 乘數。將上式寫為

$\displaystyle\begin{aligned} L(\mathbf{x},\lambda)&=(A\mathbf{x}-\mathbf{b})^T(A\mathbf{x}-\mathbf{b})+\lambda(\mathbf{x}^T\mathbf{x}-d^2)\\ &=\mathbf{x}^TA^TA\mathbf{x}-2\mathbf{x}^TA^T\mathbf{b}+\mathbf{b}^T\mathbf{b}^2+\lambda(\mathbf{x}^T\mathbf{x}-d^2). \end{aligned}$

求出偏導數並設為零，

$\displaystyle \frac{\partial L}{\partial\mathbf{x}}=2A^TA\mathbf{x}-2A^T\mathbf{b}+2\lambda\mathbf{x}=\mathbf{0}$ ，

可得優化條件式^[3]

$(A^TA+\lambda I)\mathbf{x}=A^T\mathbf{b}$ 。

如前，代入 $A$ 的奇異值分解 $A=U\Sigma V^T$ ，則

$A^TA+\lambda I=V\Sigma^T U^TU\Sigma V^T+\lambda I=V(\Sigma^T\Sigma +\lambda I)V^T=V(D^2+\lambda I)V^T$

且 $A^T\mathbf{b}=V\begin{bmatrix}D&0 \end{bmatrix}U^T\mathbf{b}$ 。設 $\mathbf{c}=(c_1,\ldots,c_m)^T=U^T\mathbf{b}$ ，優化條件式左乘 $V^T$ ，就得到簡化條件式

$\displaystyle (D^2+\lambda I)V^T\mathbf{x}=\begin{bmatrix}D&0\end{bmatrix}\mathbf{c}=\begin{bmatrix}\sigma_1c_1\\ \vdots\\ \sigma_nc_n \end{bmatrix}$ 。

若係數矩陣 $D^2+\lambda I$ 是可逆的，則存在唯一的正則最小平方解

$\displaystyle \hat{\mathbf{x}}=V(D^2+\lambda I)^{-1}\begin{bmatrix}\sigma_1c_1\\ \vdots\\ \sigma_nc_n \end{bmatrix}=V\begin{bmatrix} \frac{\sigma_1c_1}{\sigma_1^2+\hat\lambda}\\ \vdots\\ \frac{\sigma_nc_n}{\sigma_n^2+\hat\lambda} \end{bmatrix}=\sum_{i=1}^n\frac{\sigma_ic_i}{\sigma_i^2+\lambda}\mathbf{v}_i$ 。

接下來要找出 $\hat\lambda$ 使滿足 $\Vert\hat{\mathbf{x}}\Vert^2=d^2$ 。使用上面結果，寫出

$\displaystyle d^2=\Vert\hat{\mathbf{x}}\Vert^2=\Vert V^T\hat{\mathbf{x}}\Vert^2=\left\Vert\begin{bmatrix} \frac{\sigma_1c_1}{\sigma_1^2+\lambda}\\ \vdots\\ \frac{\sigma_nc_n}{\sigma_n^2+\lambda} \end{bmatrix}\right\Vert^2=\sum_{i=1}^n\left(\frac{\sigma_ic_i}{\sigma_i^2+\lambda}\right)^2$ ，

稱為特徵方程 (secular equation)。定義

$\displaystyle \phi(\lambda)=\sum_{i=1}^n\left(\frac{\sigma_ic_i}{\sigma_i^2+\lambda}\right)^2$ 。

當 $\lambda>0$ ， $\phi(\lambda)$ 是一個單調遞減函數。若 $\phi(0)>d^2$ ，則存在唯一的 $\hat\lambda>0$ 滿足 $\phi(\hat\lambda)=d^2$ ，因此保證 $D^2+\hat\lambda I$ 是可逆的。特徵方程 $\phi(\lambda)=d^2$ 是一個非線性方程，使用數值方法可求出 $\hat{\lambda}$ (例如牛頓法，見“牛頓法──非線性方程的求根方法”)，代回即得到正則最小平方解 $\hat{\mathbf{x}}$ 。

舉一個例子^[4]，

$\displaystyle \min_{\Vert\mathbf{x}\Vert=1}\left\Vert\begin{bmatrix} 2&0\\ 0&1\\ 0&0 \end{bmatrix}\mathbf{x}-\begin{bmatrix} 4\\ 2\\ 3 \end{bmatrix}\right\Vert^2$

的特徵方程為

$\displaystyle \left(\frac{8}{4+\lambda}\right)^2+\left(\frac{2}{1+\lambda}\right)^2=1$ 。

這個問題的解是 $\hat{\lambda}=4.57132$ 與 $\hat{\mathbf{x}}=(0.93334,0.35898)^T$ 。

註解
[1] 這裡使用較簡潔的寫法，約束優化問題的傳統表達式為

$\displaystyle \begin{array}{ll} \hbox{minimize}&\Vert A\mathbf{x}-\mathbf{b}\Vert^2\\ \hbox{subject to}&C\mathbf{x}=\mathbf{d}.\end{array}$

[2] 取自 G.H. Golub and C.F. Van Loan, Matrix Computations, 2nd edition, 1989, pp 568.
[3] 對於 $\lambda>0$ ，Tikhonov 正則化問題

$\displaystyle \min_{\mathbf{x}}\Vert A\mathbf{x}-\mathbf{b}\Vert^2+\lambda\Vert\mathbf{x}\Vert^2$

的最佳解滿足 $(A^TA+\lambda I)\mathbf{x}=A^T\mathbf{b}$ (見“每週問題 September 5, 2016”)。
[4] 同 [1], pp 565.

4 Responses to 約束最小平方問題

Meiyue Shao says:

12/20/2016 at 4:27 pm

那个n+p阶矩阵可逆的充要条件应当改成充分条件，必要性没有验证（事实上必要性不成立，比如n=p, A=0, C=I）。
或者把rank(A)=n作为大前提，这样那个n+p阶矩阵可逆的充要条件就是rank(C)=p。

- ccjou says:
  
  12/20/2016 at 4:39 pm
  
  謝謝指正，我疏忽了，確實只是充分條件。
  
kaerfuka says:

12/20/2016 at 10:59 pm

Reblogged this on 站点标题.

Lin says:

12/27/2016 at 7:20 pm

rank(A)=n这个条件太关键了。如果这个条件不成立，很多推导都不成立了。

	陳倍恩 on 線性代數的第一堂課──矩陣乘法的定義
	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨