【動手學MVG】矩陣分解與線性方程組的關系，求解線性方程組實戰代碼

文章目錄

- 本文解決的問題
- QR分解
- - 投影矩陣分解得到内外參
- SVD分解
- 最小二乘問題
- - 最小二乘問題的求解方法
  - - (列)滿秩最小二乘問題
    - - 正規化方法
      - QR分解方法
      - SVD 分解方法
    - 虧秩最小二乘問題
    - 齊次最小二乘問題
- 數值穩定性問題
- Ax=0與Ax=b的選擇
- 代碼實作
- - 利用SVD求解非齊次線性方程組Ax=b
  - 利用SVD求解齊次線性方程組Ax=0
- 附錄——一些概念
- 參考資料

在MVG(多視圖幾何)和機器學習領域，求解線性方程組幾乎是所有算法的根本，本文旨在幫助讀者搞懂矩陣分解與線性方程組的關系，并給出利用SVD求解線性方程組的實戰代碼。

本文解決的問題

看完本文後，應解決以下問題：

投影矩陣怎麼通過QR分解得到相機内外參？
求解非齊次線性方程組Ax=b的方法有哪些？
- 其他一些方法通用性不大，比如克萊姆法則隻用于方陣，這裡不考慮。
- 高斯消元法 | 正常方法。主要思路是講矩陣通過初等行變換，轉變為行階矩形，然後後向代入(或者稱為回代算法)可以很容易求解。
  - 根據增系數矩陣A和廣矩陣B = (A, b) 的秩可判斷解的個數。
  - 對于齊次線性方程組Ax=0，用初等行變換将A變為行最簡形矩陣A’，即可得到同解方程組A’x=0，然後計算
  - 對于非齊次線性方程組Ax=0，用初等行變換将增廣矩陣B = (A, b)變為行最簡形矩陣A’，即可得到同解方程組A’x=0，然後計算
- 方程Ax=b ,其中 A 為mxn。
  - 當m >n時，對于這樣的方程不能直接利用高斯或者其他的分解法進行求解 , 往往需要把方程轉換為另外一種更容易求解的模式，也就是我們常常說的各種分解。
  - 當m<n時，有無窮解，高斯消元等正常方法能解，但是如果我們要求使得二範數最小的最小二乘解，還是需要用SVD分解。
  - 當m=n時，一般可以用高斯消元等方法解。但是如果秩小于n，此時有無窮解，如果我們要求使得二範數最小的最小二乘解，還是需要用SVD分解。
- 通過分解，将矩陣分解為特殊的形式，這些形式求逆或者對方程求解比較簡單，比如上三角矩陣、下三角矩陣、對角矩陣、正交矩陣(求逆簡單，就是其轉職)。分解完後，方程可以轉換為比較好求解的形式，用正常方法求解即可。常用方法有：LU分解、QR分解、Cholesky分解
- 通過分解，求出矩陣廣義逆矩陣 A + A^+ A+ ，則解為 x = A + b x = A^+b x=A+b 。常用方法：正規方程( x = ( A T A ) − 1 A T b x = (A^TA)^{-1}A^Tb x=(ATA)−1ATb)、SVD分解
總結：

參考：知呼
- LU需要矩陣A可逆，Cholesky是其特殊情況，Cholesky分解的時間複雜度是 m n 2 + ( 1 / 3 ) n 3 mn^2+(1/3)n^3 mn2+(1/3)n3
- QR分解的時間複雜度是 2 m n 2 2 mn^2 2mn2 ，因為m>n，是以QR往往比Cholesky慢
- 在實際應用中，因為數值穩定性的要求，稠密矩陣往往用QR求解，對于大型的稀疏矩陣則多用Cholesky分解。
- SVD分解，比QR分解慢，但是比QR分解穩定。
非齊次線性方程組Ax=b，什麼時候有精确解，什麼時候隻有最小二乘解？
- 若矩陣A為列滿秩矩陣，
  - 當A為方陣時，有精确解
  - 當A為mxn的矩陣(m>n)時，方程個數大于未知數個數；這時候多出來的方程可能會與前面的方程沖突，這時候沒有精确解，隻有最小二乘解，是以不能直接用高斯消元等方法，隻能用分解的方法找其主要特征，然後再正常求解。若多出來的方程與前面方程不沖突，即多出來的方程都可以由前面的方程線性組合得到，這時候還是有精确解的。
- 若矩陣A為虧秩矩陣 | 即矩陣A獨立方程個數小于未知變量個數，這時候有無窮多個解。
  - 此時也隻有最小二乘解
TODO：為什麼要用QR分解來求解非齊次線性方程組Ax=b？
- A x = b ⇒ Q R x = b ⇒ R x = Q T b Ax=b \Rightarrow QRx=b \Rightarrow Rx=Q^Tb Ax=b⇒QRx=b⇒Rx=QTb 因為R是上三角矩陣，是以很容易求得方程組的解。 Q T Q^T QT表示Q的轉置。由于R是三角形，是以方程組通過右側的簡單矩陣矢量乘法和向後替換來求解。
- QR分解數值穩定性較好。
TODO：為什麼要用SVD分解來求解非齊次線性方程組Ax=b？
為什麼大多數解線性方程組用LU分解，而不用QR分解？因為QR分解時間複雜度較高，且LU分解容易并行。注意：SVD分解不容易并行看這裡
TODO：求解線性方程組Ax=b時，什麼時候用QR分解，什麼時候用SVD分解？
齊次線性方程組Ax=0的最小二乘解怎麼求？
- 對 A 進行SVD分解後， A 的最小奇異值的右奇異向量就是方程組的一個解。

QR分解

注意：除特殊說明外，讨論的矩陣均為實矩陣。

QR分解是一種正交三角分解。

定義 8.1.1 如果非奇異實矩陣 A 能夠表示為正交矩陣 Q 與上三角矩陣 R 的積，即

A ＝ Q R (8.1.1) A＝QR \tag{8.1.1} A＝QR(8.1.1)

與 QR 分解類似，還有 RQ、QL， LQ 分解，其中 L 表示下三角矩陣，R表示上三角矩陣。矩陣的 QR 分解、 RQ 分解、 QL 分解、LQ 分解統稱為矩陣的正交三角分解。其它類型分解可通過類似 QR 分解的方法獲得。
定理 8.1.1 對任意非奇異實矩陣 A 總可以分解為正交矩陣 Q 與上三角矩陣 R 的積，如果要求上三角陣 R 的對角元素均為正數，則分解是唯一的。

非奇異矩陣是行列式不為 0 的矩陣，也就是可逆矩陣。意思是n 階方陣 A 是非奇異方陣的充要條件是 A 為可逆矩陣，也即A的行列式不為零。即矩陣(方陣)A可逆與矩陣A非奇異是等價的概念。——參考360百科

反之則是非奇異矩陣。

注意：隻有方陣才有奇異矩陣和非奇異矩陣的說法。
定理 8.1.2 定理 8.1.1 可以推廣到非方陣的情形。對任意列滿秩的實矩陣 A 總可以分解為列正交的矩陣 Q 與上三角矩陣 R 的積，如果

要求上三角陣 R 的對角元素均為正數，則這種分解是唯一的。
上面用 Schmidt 正交化方法構造性地證明了非奇異矩陣總可以進行 QR 分解，但是在實踐中并不使用 Schmidt 正交化方法對矩陣作 QR 分解，而是利用實用的方法： Givens 方法或 Houesholder方法。
- Givens 旋轉方法，對于 n 階矩陣需要作 n(n-1)/2 個 Givens 旋轉矩陣的積，計算量較大。對于高階矩陣而言，用下述 Householder（反射）矩陣變換進行 QR 分解更有效。它隻需要作(n-1)個 Householder 矩陣的積，其計算量大約是 Givens 旋轉方法的一半。

投影矩陣分解得到内外參

參考：《計算機視覺中的數學方法》——吳朝福 P182、知乎、Camera Calibration

RQ 分解是實作錄影機内參數、外參數求解的最為簡便的方法。令錄影機内參數矩陣為 K，外參數矩陣為 (R, t) ，由于此時的錄影機矩陣是歐氏的，是以可寫成下述形式：

P = α K [ R , t ] = [ α K R , α K t ] P = αK[R, t] = [αKR, αKt] P=αK[R,t]=[αKR,αKt]

将 P 表示成 P = [ H , p 4 ] P = [H, p_4] P=[H,p4] ，則有

H = α K R p 4 = α K t H = αKR\\ p_4 = αKt H=αKRp4=αKt

其中， p 4 p_4 p4 表示P矩陣的第4列。對 H 作 RQ 分解： H = K ^ R ^ H = \hat{K}\hat{R} H=K^R^ ，其中 K ^ \hat{K} K^ 是對角元均為正數的上三角矩陣， R ^ \hat{R} R^ 為正交矩陣，并且種分解是唯一的。三角矩陣 K ^ \hat{K} K^ 與錄影機參數矩陣相差一個正常數倍，由于内參數矩陣最後一個元素為1，是以内參數矩陣必為： K = K ^ 33 − 1 K ^ K = \hat{K}_{33}^{-1}\hat{K} K=K^33−1K^ ，其中 K ^ 33 \hat{K}_{33} K^33 是矩陣 K ^ \hat{K} K^ 的第(3, 3)元素。最後得到的K有如下形式：

K = [ f x k c x f y c y 1 ] K = \begin{bmatrix} f_x & k & c_x \\ & f_y & c_y \\ & & 1 \end{bmatrix} K=⎣⎡fxkfycxcy1⎦⎤

其中 k = t a n θ k = tan\theta k=tanθ ， θ \theta θ 是圖像軸之間的夾角。

如果正交矩陣 R ^ \hat{R} R^ 是一個旋轉矩陣，則它就是錄影機關于世界坐标系的姿态。如果 R ^ \hat{R} R^ 是不是旋轉矩陣，則表明所估計的錄影機矩陣與實際的(歐氏)錄影機矩陣反号，即 P 中的齊次子的符号 s g n α = − 1 sgnα = −1 sgnα=−1 ，于是 R = − R ^ R = −\hat{R} R=−R^ 就是所要求解旋轉矩陣。對于平移參數 t，可通過下式來确定：

{ t = K − 1 P 4 , 若 R ^ 是旋轉矩陣 t = − K − 1 P 4 , 若 R ^ 不是旋轉矩陣 \left\{\begin{array}{l} t = K^{-1}P_4,若\hat{R}是旋轉矩陣\\ t = -K^{-1}P_4,若\hat{R}不是旋轉矩陣 \end{array}\right. {t=K−1P4,若R^是旋轉矩陣t=−K−1P4,若R^不是旋轉矩陣

這樣，就得到了錄影機内、外參數。

SVD分解

定理 8.3.2(奇異值分解， SVD) 設 A ∈ R r m × n A \in R_{r}^{m \times n} A∈Rrm×n ，則存在 m 階正交矩陣 U 和 n 階正交矩陣 V 使

A = U [ Σ 0 0 0 ] V T A = U\begin{bmatrix} \Sigma & \mathbf{0} \\ \mathbf{0} & 0 \end{bmatrix}V^T A=U[Σ000]VT

其中 Σ = d i a g ( σ 1 , σ 2 , . . . , σ r ) Σ = diag(σ_1,σ_2 ,...,σ_r) Σ=diag(σ1,σ2,...,σr) ，即 σ 1 , σ 2 , . . . , σ r σ_1,σ_2 ,...,σ_r σ1,σ2,...,σr 是 A 的非零奇異值。式(8.3.2)通常記為

A = U D V T A = UDV^T A=UDVT

注意：奇異值分解既不要求矩陣 A 是可逆的方陣，也不要求它是方陣。

定理 8.3.4 設 A ∈ R r m × r A ∈ R_{r}^{m \times r} A∈Rrm×r 的奇異值為 σ 1 ≥ σ 2 ≥ . . . ≥ σ r > σ r + 1 = . . . = σ n = 0 σ_1 ≥ σ_2 ≥ ... ≥ σ_r > σ_{r+1} = ... = σ_n = 0 σ1≥σ2≥...≥σr>σr+1=...=σn=0 ， ( A + Q ) ∈ R r m × r (A+Q) ∈ R_{r}^{m \times r} (A+Q)∈Rrm×r 的奇異值為 τ 1 ≥ τ 2 ≥ . . . ≥ τ r > τ r + 1 = . . . = τ n = 0 τ_1 ≥ τ_2 ≥ ... ≥ τ_r > τ_{r+1} = ... = τ_n = 0 τ1≥τ2≥...≥τr>τr+1=...=τn=0 ，則必有

∣ σ j − τ j ∣ ≤ ∣ ∣ Q ∣ ∣ 2 ( j = 1 , 2 , . . . , n ) |σ_j − τ_j| ≤ || Q ||_2 ( j = 1,2,...,n) ∣σj−τj∣≤∣∣Q∣∣2(j=1,2,...,n)

定理表明，矩陣 A 在攝動 Q 下，奇異值的變化量不超過 ∣ ∣ Q ∣ ∣ 2 || Q ||_2 ∣∣Q∣∣2 。是以，矩陣奇異值的計算具有良好的數值穩定性質。
奇異值分解還可以表示為

A = σ 1 u 1 v 1 T + σ 2 u 2 v 2 T + . . . + σ r u r v r T (8.3.8) A = σ_1u_1v_1^T + σ_2u_2v_2^T + ... + σ_ru_rv_r^T \tag{8.3.8} A=σ1u1v1T+σ2u2v2T+...+σrurvrT(8.3.8)

并且稱 u i u_i ui 為奇異值 σ i σ_i σi 的左奇異向量， v i v_i vi 為奇異值 $σ_i $的右奇異向量。

順便一提：SVD的代碼實作中，可以選擇隻計算U，也可以U和V一起計算，他們的時間複雜度是不一樣的。是以你可以根據自己的需求，選擇合适的方法。

最小二乘問題

考慮線性系統：

A x = b ( A ∈ R m × n ( m > n ) , b ∈ R m , x ∈ R n ) (8.4.1) Ax = b (A ∈ R^{m \times n} (m > n),b ∈ R^m, x ∈ R^n ) \tag{8.4.1} Ax=b(A∈Rm×n(m>n),b∈Rm,x∈Rn)(8.4.1)

正常來看，這個方程是沒有解的的，但在數值計算領域，我們通常的是其最小二乘解。即求 x ∈ R n x ∈ R^n x∈Rn 使得

∣ ∣ A x − b ∣ ∣ 2 = m i n { ∣ ∣ A v − b ∣ ∣ 2 : v ∈ R n } (8.4.2) || Ax-b ||_2 = min\{|| Av-b ||^2 : v \in R^n\} \tag{8.4.2} ∣∣Ax−b∣∣2=min{∣∣Av−b∣∣2:v∈Rn}(8.4.2)

記

X L S = { x ∈ R n : x 是 ( 8.4.1 ) 的解 } X_{LS} = \{ x \in R^n : x 是(8.4.1)的解\} XLS={x∈Rn:x是(8.4.1)的解}

則稱 X L S X_{LS} XLS 是最小二乘問題(8.4.1)的解集； X L S X_{LS} XLS中範數最小者稱為最小範數解，并記作 X L S X_{LS} XLS ，即

∣ ∣ x L S ∣ ∣ 2 = m i n { ∣ ∣ x ∣ ∣ 2 : x ∈ X L S } || x_{LS} ||_2 = min\{|| x ||_2 : x \in X_{LS}\} ∣∣xLS∣∣2=min{∣∣x∣∣2:x∈XLS}

命題 8.4.1 x ∈ X L S ⇔ A T ( A x − b ) = 0 x \in X_{LS} ⇔ A^T (Ax − b) = 0 x∈XLS⇔AT(Ax−b)=0 。其中，方程 A T ( A x − b ) = 0 A^T (Ax − b) = 0 AT(Ax−b)=0 稱為 A x − b = 0 Ax − b = 0 Ax−b=0 的正規方程。

根據命題 8.4.1，有下述推論：

推論 8.4.1： (1) X L S X_{LS} XLS 是凸集； (2) x L S x_{LS} xLS 是唯一的； (3) X L S = { x L S } X_{LS} = \{x_{LS}\} XLS={xLS} 的充分必要條件是 r a n k ( A ) = n rank(A) = n rank(A)=n 。

最小二乘問題的求解方法

(列)滿秩最小二乘問題

如果(8.4.1)中的矩陣 A 是列滿秩的，即 rank(A) = n ，則稱它為滿秩最小二乘問題。下面考慮滿秩最小二乘問題的數值算法。

正規化方法

将求解問題(8.4.1)轉化為求解正規化方程組

A T A x = A T b (8.4.8) A^TAx = A^Tb \tag{8.4.8} ATAx=ATb(8.4.8)

由于 rank(A) = n ，是以 A T A A^TA ATA 是對稱正定矩陣，正定矩陣可逆，是以 x L S = ( A T A ) − 1 A T b x_{LS} = (A^TA)^{-1}A^Tb xLS=(ATA)−1ATb 但由于逆矩陣 ( A T A ) − 1 (A^TA)^{-1} (ATA)−1 不是特殊的矩陣，求解逆矩陣複雜度會很高，是以一般用其他方法。

因為 A T A A^TA ATA 是對稱正定矩陣，是以(8.4.8)的唯一解 x L S x_{LS} xLS 還可以用 **Cholesky 分解法(它是LU分解的特殊情況)**求得。基本步驟為：

計算 C = A T A , d = A T b C = A^TA, d = A^Tb C=ATA,d=ATb ；
對 C 進行 Cholesky 分解 C = G G T C = GG^T C=GGT ；其中 G 為對角元素均大于零的上三角矩陣；
求解三角方程 G y = d Gy = d Gy=d 以及 G T x = y G^Tx = y GTx=y 。由于是上三角矩陣，後向疊代的高斯消元即可。

Cholesky 分解是LU分解的特殊情況。這種算法的缺點很明顯：首先， A T A A^TA ATA有時候不可逆，不能得出結果，其次，即便可逆， A T A A^TA ATA數值穩定性不好，會造成誤差。

QR分解方法

正規化方法通常比較低效，下面介紹QR分解法。

設 A 有 QR 分解：

A = Q m × m [ R n × n 0 ] m × n = Q 1 R A = Q_{m \times m}\begin{bmatrix} R_{n \times n} \\ \mathbf{0} \end{bmatrix}_{m \times n} = Q_1R A=Qm×m[Rn×n0]m×n=Q1R

其中 Q ∈ R m × m Q \in R^{m×m} Q∈Rm×m 是正交矩陣， [ R 0 ] ∈ R m × n \begin{bmatrix} R \\ \mathbf{0}\end{bmatrix} \in R^{m×n} [R0]∈Rm×n， Q1 是 Q 的前 n 列組成的矩陣，即 Q = ( Q 1 n ∣ Q 2 m − n ) Q = ( \underset{n}{Q_1} | \underset{m-n}{Q_2} ) Q=(nQ1∣m−nQ2) ， R ∈ R n × n R \in R^{n×n} R∈Rn×n 是對角線上元素均為正數的上三角矩陣。

由于正交矩陣保持範數不變，是以問題 (8.4.1) 等價于

∣ ∣ Q T ( A x − b ) ∣ ∣ 2 = m i n { ∣ ∣ Q T ( A v − b ) ∣ ∣ 2 : v ∈ R n } ||Q^T(Ax-b)||_2 = min\{||Q^T(Av-b)||_2:v\in R^n\} ∣∣QT(Ax−b)∣∣2=min{∣∣QT(Av−b)∣∣2:v∈Rn}

記

d = Q T b = [ Q 1 T Q 2 T ] b = [ d 1 d 2 ] d = Q^Tb = \begin{bmatrix} Q_1^T \\ Q_2^T \end{bmatrix}b = \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} d=QTb=[Q1TQ2T]b=[d1d2]

則有

∣ ∣ Q T ( A x − b ) ∣ ∣ 2 2 = ∣ ∣ [ R 0 ] x − [ d 1 d 2 ] ∣ ∣ 2 2 = ∣ ∣ R x − d 1 ∣ ∣ 2 2 + ∣ ∣ d 2 ∣ ∣ 2 2 ||Q^T(Ax-b)||_2^2 = || \begin{bmatrix} R \\ 0 \end{bmatrix} x- \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} ||_2^2 = ||Rx-d_1||_2^2 + ||d_2||_2^2 ∣∣QT(Ax−b)∣∣22=∣∣[R0]x−[d1d2]∣∣22=∣∣Rx−d1∣∣22+∣∣d2∣∣22

是以， x 是 (8.4.1) 的解當且僅當 x 是方程 R x = d 1 Rx=d_1 Rx=d1 的解。這樣 (8.4.1) 的解可由上三角方程組 R x = d 1 Rx=d_1 Rx=d1 求得。

而上三角方程組用後向疊代的高斯消元是很好求的。

QR 分解方法的基本步驟如下：

求 A 的 QR 分解；
計算 d 1 = Q 1 T b d_1 = Q_{1}^{T}b d1=Q1Tb ；
解方程組 R x = d 1 Rx = d_1 Rx=d1 。

注意： QR 分解方法比正規化方法有較好的數值穩定性，并且計算結果比正規化方法要精确。當然， QR 方法比正規化方法會付出更大的計算代價。

SVD 分解方法

由于 rank(A) = n ，是以 A 必有下述形式的 SVD 分解

A = U [ Σ n 0 ] V T A = U\begin{bmatrix} \Sigma_n \\ \mathbf{0} \end{bmatrix}V^T A=U[Σn0]VT

于是，A的廣義逆矩陣 A + = V [ Σ n − 1 , 0 ] U T A^+ = V[\Sigma_n^{-1}, \mathbf{0}]U^T A+=V[Σn−1,0]UT 。是以，問題(8.4.1)的解為

x = A + b = V ( Σ n − 1 , 0 ) U T b = u 1 T b σ 1 v 1 + u 2 T b σ 2 v 2 + ⋯ + u n T b σ n v n = ∑ j = 1 n u j T b σ j v j (8.4.9) x = A^+b = V(\Sigma_n^{-1}, 0)U^Tb = \frac{u_1^Tb}{\sigma_1}v_1 + \frac{u_2^Tb}{\sigma_2}v_2 + \dots + \frac{u_n^Tb}{\sigma_n}v_n = \sum_{j=1}^{n}\frac{u_j^Tb}{\sigma_j}v_j \tag{8.4.9} x=A+b=V(Σn−1,0)UTb=σ1u1Tbv1+σ2u2Tbv2+⋯+σnunTbvn=j=1∑nσjujTbvj(8.4.9)

上式給出了 SVD 分解方法關于最小二乘解的計算公式。

虧秩最小二乘問題

如果在最小二乘問題(8.4.1)中，矩陣 A 是虧秩的，即 r = rank(A) < n 。此時(8.4.1)有無窮多解，在上面介紹的處理滿秩問題(8.4.1)的正規化方法， QR 分解方法都将失敗，不能給出最小範數解 x L S x_{LS} xLS 。但是 SVD 分解方法仍然有效。具體地說，若 rank(A) = r ，則 A 有 SVD 分解：

A = U [ Σ r 0 0 0 ( m − r ) × ( n − r ) ] V T A = U\begin{bmatrix} \Sigma_r & \mathbf{0} \\ \mathbf{0} & \mathbf{0}_{(m-r)\times(n-r)} \end{bmatrix}V^T A=U[Σr000(m−r)×(n−r)]VT

是以，

x = A + b = V [ Σ r − 1 0 0 0 ] U T b = ∑ j = 1 r u j T b σ j v j (8.4.10) x = A^+b = V\begin{bmatrix} \Sigma_r^{-1} & \mathbf{0} \\ \mathbf{0} & 0 \end{bmatrix}U^Tb = \sum_{j=1}^{r}\frac{u_j^Tb}{\sigma_j}v_j \tag{8.4.10} x=A+b=V[Σr−1000]UTb=j=1∑rσjujTbvj(8.4.10)

在這裡，可以看到 SVD 分解在數值計算中的作用，不論是滿秩的還是虧秩的最小二乘問題， SVD分解方法總能給出它們的求解計算公式，并且有統一的形式。

齊次最小二乘問題

對于齊次線性方程組 A x = 0 Ax=0 Ax=0 的最小二乘解，有一點不同。

考慮齊次線性方程組：

A x = 0 ( A ∈ R m × n ( m > n ) , x ∈ R n , m > n ) Ax = 0 (A ∈ R^{m×n} (m > n), x ∈ R^n ,m > n) Ax=0(A∈Rm×n(m>n),x∈Rn,m>n)

對應的最小二乘問題是

∣ ∣ A x ∣ ∣ 2 = m i n { ∣ ∣ A v ∣ ∣ 2 : v ∈ R n } (8.4.13) || Ax ||_2 =min\{||Av||_2 : v \in R^n\} \tag{8.4.13} ∣∣Ax∣∣2=min{∣∣Av∣∣2:v∈Rn}(8.4.13)

顯然， x = 0 總是上述最小二乘問題的最小範數解。在實際中，人們所關心的并非是齊次最小二乘問題的零解，而是它的非零解。是以，我們總是考慮相應的限制最小二乘問題：(之是以取機關特征向量v，是因為若v是一個解，那乘以任意一個尺度s後還是方程的解)

∣ ∣ A x ∣ ∣ 2 = m i n { ∣ ∣ A v ∣ ∣ 2 : v ∈ R n , ∣ ∣ v ∣ ∣ 2 = 1 } (8.4.13) || Ax ||_2 =min\{||Av||_2 : v \in R^n, ||v||_2=1\} \tag{8.4.13} ∣∣Ax∣∣2=min{∣∣Av∣∣2:v∈Rn,∣∣v∣∣2=1}(8.4.13)

或者等價地寫成：

{ m i n ∣ ∣ A x ∣ ∣ 2 2 s u b j e c t t o ∣ ∣ x ∣ ∣ 2 = 1 (8.4.15) \left\{ \begin{aligned} & min \; ||Ax||_2^2 \\ & subject \; to \; ||x||_2 = 1 \end{aligned} \right. \tag{8.4.15} {min∣∣Ax∣∣22subjectto∣∣x∣∣2=1(8.4.15)

推論 8.4.6 若 rank(A) = n −1 時， (8.4.15)有唯一解（rank(A)=n時齊次線性方程組隻有零解），這個唯一解是 A T A A^TA ATA 的零特征值的機關特征向量，或者說是 A 的零奇異值的右奇異向量。當資料有誤差時， A T A A^TA ATA 的最小特征值的機關特征向量(或者說是 A 的最小奇異值的右奇異向量)是(8.4.15)的一個解(當資料無誤差時，最小特征值就是零)。

個人證明一下，可能不嚴謹：

矩陣 A ∈ R r m × n ( m > n ) , x ∈ R n , m > n A ∈ R_r^{m×n} (m > n), x ∈ R^n ,m > n A∈Rrm×n(m>n),x∈Rn,m>n .

一、若矩陣A的秩rank(A) = r ( r<n)，對矩陣A進行SVD分解

A = U [ Σ r 0 0 0 ( m − r ) × ( n − r ) ] V T A = U\begin{bmatrix} \Sigma_r & \mathbf{0} \\ \mathbf{0} & \mathbf{0}_{(m-r)\times(n-r)} \end{bmatrix}V^T A=U[Σr000(m−r)×(n−r)]VT

其中正交矩陣 U m × m U_{m \times m} Um×m 和 V n × n V_{n \times n} Vn×n 的列向量分别為 u 1 , u 2 , . . . , u m u_1,u_2,...,u_m u1,u2,...,um 和 v 1 , v 2 , . . . , v m v_1,v_2,...,v_m v1,v2,...,vm ，則将SVD分解寫成向量的形式：

A = σ 1 u 1 v 1 T + σ 2 u 2 v 2 T + . . . + σ r u r v r T A = \sigma_1u_1v_1^T + \sigma_2u_2v_2^T + ... + \sigma_ru_rv_r^T A=σ1u1v1T+σ2u2v2T+...+σrurvrT

其中 σ i \sigma_i σi 是矩陣A的奇異值，也是 Σ r \Sigma_r Σr 對角線上元素，并且 σ 1 ≥ σ 2 ≥ . . . ≥ σ r ≥ 0 \sigma_1 \ge \sigma_2 \ge ... \ge \sigma_r \ge 0 σ1≥σ2≥...≥σr≥0 。則

∣ ∣ A x ∣ ∣ 2 2 = ∣ ∣ σ 1 u 1 v 1 T x + σ 2 u 2 v 2 T x + . . . + σ r u r v r T x ∣ ∣ 2 2 ||Ax||_2^2 = ||\sigma_1u_1v_1^Tx + \sigma_2u_2v_2^Tx + ... + \sigma_ru_rv_r^Tx||_2^2 ∣∣Ax∣∣22=∣∣σ1u1v1Tx+σ2u2v2Tx+...+σrurvrTx∣∣22

由于 V n × n V_{n \times n} Vn×n 是正交矩陣，它的列向量都為機關向量并且互相正交（向量積為0），是以當取 x = v i , i > r , 滿足 ∣ ∣ x ∣ ∣ 2 = 1 x = v_i, i \gt r, 滿足||x||_2 = 1 x=vi,i>r,滿足∣∣x∣∣2=1 時， v i v_i vi 與 v 1 , v 2 , . . . , v r v_1,v_2,...,v_r v1,v2,...,vr 都正交，向量積為0，是以 ∣ ∣ A x ∣ ∣ 2 2 = 0 ||Ax||_2^2 = 0 ∣∣Ax∣∣22=0 ， x = v i , i > r x = v_i, i \gt r x=vi,i>r 是(8.4.15)的解（若rank(A)=n-1則這個解是唯一解，否則若rank(A)<n-1，則這個解隻是其中的一個解，事實上，在帶限制的條件下，它所有的解就是所有0奇異值對應的右奇異向量）。

二、若矩陣A的秩rank(A) = n (列滿秩)，對矩陣A進行SVD分解

A = U [ Σ n 0 ( m − n ) × 1 ] V T A = U\begin{bmatrix} \Sigma_n \\ \mathbf{0}_{(m-n)\times1} \end{bmatrix}V^T A=U[Σn0(m−n)×1]VT

寫成向量的形式：

A = σ 1 u 1 v 1 T + σ 2 u 2 v 2 T + . . . + σ n u n v n T A = \sigma_1u_1v_1^T + \sigma_2u_2v_2^T + ... + \sigma_nu_nv_n^T A=σ1u1v1T+σ2u2v2T+...+σnunvnT

當取 x = v n , 滿足 ∣ ∣ x ∣ ∣ 2 = 1 x = v_n, 滿足||x||_2 = 1 x=vn,滿足∣∣x∣∣2=1 時， v n v_n vn 與 v 1 , v 2 , . . . , v n − 1 v_1,v_2,...,v_{n-1} v1,v2,...,vn−1 都正交，向量積為0，而 v n , u n v_n,u_n vn,un 都是機關正交向量，有 v n T v n = 1 , ∣ ∣ u n ∣ ∣ 2 2 = 1 v_n^Tv_n = 1, ||u_n||_2^2 = 1 vnTvn=1,∣∣un∣∣22=1 ，是以 ∣ ∣ A x ∣ ∣ 2 2 = ∣ ∣ σ n u n ∣ ∣ 2 2 = σ n 2 ∣ ∣ u n ∣ ∣ 2 2 = σ n 2 ||Ax||_2^2 = ||\sigma_nu_n||_2^2 = \sigma_n^2||u_n||_2^2 = \sigma_n^2 ∣∣Ax∣∣22=∣∣σnun∣∣22=σn2∣∣un∣∣22=σn2 ，而 σ n \sigma_n σn 是所有奇異值中最小的。

注意：這裡沒有證明為什麼 x = v n x = v_n x=vn 是此時最小二乘解，隻是說明了當取這個值時，似乎它能使得(8.4.15)在限制條件下最小。

第二種情況還可以參考《計算機視覺中的數學方法》：

【動手學MVG】矩陣分解與線性方程組的關系，求解線性方程組實戰代碼

另外，看一下《計算機視覺中的多視圖幾何第一版》，講得特别好，言簡意赅且容易了解：

【動手學MVG】矩陣分解與線性方程組的關系，求解線性方程組實戰代碼

數值穩定性問題

參考：視覺SLAM常見的QR分解SVD分解等矩陣分解方式求解滿秩和虧秩最小二乘問題（最全的方法分析總結）

為什麼需要矩陣分解?

一方面，由于一些超定方程或欠定方程，沒有精确解，隻能通過分解矩陣的方式求其最小二乘解；
另一方面，當資料量很大時，将一個矩陣分解為若幹個矩陣的乘積可以大大降低存儲空間；
其次，可以減少真正進行問題處理時的計算量，畢竟算法掃描的元素越少完成任務的速度越快，這個時候矩陣的分解是對資料的一個預處理；
再次，矩陣分解可以高效和有效的解決某些問題；
最後，矩陣分解可以提高算法數值穩定性，關于這一點下面有進一步的說明。

當資料不存在誤差時，用一般的方法求解線性方程組即可；但是當我們得到的資料有誤差，或者由于計算機精度問題導緻誤差時，一般的方法求解線性方程組往往會産生很大誤差，這裡舉例說明。例如方程：

{ 5 x 1 + 7 x 2 = 0.7 7 x 1 + 10 x 2 = 1 ⇒ A x = b A = [ 5 7 7 10 ] , x = [ x 1 x 2 ] , b = [ 0.7 1 ] . \left\{ \begin{aligned} & 5x_1 + 7x_2 = 0.7 \\ & 7x_1 + 10x_2 = 1 \end{aligned} \right. \Rightarrow Ax=b \\ A = \begin{bmatrix} 5 & 7 \\ 7 & 10 \end{bmatrix}, x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, b = \begin{bmatrix} 0.7 \\ 1 \end{bmatrix}. {5x1+7x2=0.77x1+10x2=1⇒Ax=bA=[57710],x=[x1x2],b=[0.71].

直接對方程組求解可以得到：

x = [ 0.0 0.1 ] x = \begin{bmatrix} 0.0 \\ 0.1 \end{bmatrix} x=[0.00.1]

現在對方程中的 b 進行一個微小擾動：

b = [ 0.69 1.01 ] , 其中擾動項為 [ − 0.01 0.01 ] b = \begin{bmatrix} 0.69 \\ 1.01 \end{bmatrix}, 其中擾動項為 \begin{bmatrix} -0.01 \\ 0.01 \end{bmatrix} b=[0.691.01],其中擾動項為[−0.010.01]

擾動後，直接對方程組求解可以得到：

x = [ − 0.17 0.22 ] x = \begin{bmatrix} -0.17 \\ 0.22 \end{bmatrix} x=[−0.170.22]

結果變成了上式，可以看出當方程組中的常數矩陣發生微小的擾動時，會導緻最終的結果發生較大的變化，這種結果的不穩定不是因為方程求解的方法，而是方程組矩陣本身的問題。這回給我們帶來很大的危害，例如像上面的情況，計算機求解出現舍入誤差，矩陣本身不好的話會直接導緻結果失敗。

當對矩陣A或者b進行小擾動的時候，最後影響結果的是 ∣ ∣ A ∣ ∣ ⋅ ∣ ∣ A − 1 ∣ ∣ ||A|| \cdot ||A^{-1}|| ∣∣A∣∣⋅∣∣A−1∣∣，與矩陣病态性成正相關。定義矩陣的條件數 c o n d ( A ) = ∣ ∣ A ∣ ∣ ⋅ ∣ ∣ A − 1 ∣ ∣ cond(A)= ||A|| \cdot ||A^{-1}|| cond(A)=∣∣A∣∣⋅∣∣A−1∣∣ 來描述矩陣的病态程度，一般認為小于100的為良态，條件數在100到1000之間為中等程度的病态，條件數超過1000存在嚴重病态。以上面的矩陣A為例，采用2範數的條件數 c o n d ( A ) = 222.9955 cond(A)=222.9955 cond(A)=222.9955 矩陣處于中等病态程度。

矩陣其實就是一個給定的線性變換，特征向量描述了這個線性變換的主要方向，而特征值描述了一個特征向量的長度在該線性變換下縮放的比例，對開篇的例子進一步觀察發現，A是個對稱正定矩陣，A的特征值分别為 λ 1 = 14.93303437 λ1=14.93303437 λ1=14.93303437 和 λ 2 = 0.06696563 λ2=0.06696563 λ2=0.06696563 ，兩個特征值在數量級上相差很大，這意味着b發生擾動時，向量x在這兩個特征向量方向上的移動距離是相差很大的——對于λ1對應的特征向量隻需要微小的移動就可到達b的新值，而對于λ2，由于它比起λ1太小了，是以需要x做大幅度移動才能到達b的新值，于是悲劇就發生了.

Ax=0與Ax=b的選擇

對于有些問題比較靈活，既可以建構出Ax=0的方程組，也可以建構出Ax=b的方程組，那麼該如何選擇呢？

根據我看到的，好像大都選擇建構成Ax=0的形式去解，我也不知道為什麼，這裡說一下我個人的看法。

解Ax=b需要求解A的SVD分解 A = U D V T A=UDV^T A=UDVT ，解Ax=0的最小二乘解，隻需要用V的最後一列向量作為解即可，而SVD分解是可以講U、D和V分開求的，可以隻求部分，隻求部分時，時間複雜度會大大降低，是以Ax=0的一個優勢就是：不用求解完整的SVD，進而減少時間複雜度。（具體SVD的求法，可以自行搜尋）

代碼實作

為了更好地了解本章的内容，下面給出兩個求解線性方程組的實戰代碼。

如果你不知道怎麼運作下面代碼，可以點選這裡下載下傳完整的工程，工程采用VS2015建構。

利用SVD求解非齊次線性方程組Ax=b

下面代碼僅依賴于eigen庫，在VS2015上調試通過。

#include <iostream>
#include <ctime>
#include "Eigen/Dense"

using namespace std;
using namespace Eigen;

const int MOD = 100;
int main() {

	// 構造輸入資料
	uint32_t seed = time(0);
	srand(seed); // 随機初始化種子
	MatrixXd A(3, 3);
	// 當 A 是一個正常的可逆矩陣時，x的解應該是唯一的
	A << rand() % MOD, rand() % MOD, rand() % MOD,
		rand() % MOD, rand() % MOD, rand() % MOD, 
		rand() % MOD, rand() % MOD, rand() % MOD;

	 當 A 的秩為2時，Ax=b是一個欠定方程組，x的解有無數個，求解出來的可能與真值不同，但是Ax的結果仍然為b
	 注釋這一行去檢視一下結果
	//A << rand() % MOD, rand() % MOD, rand() % MOD,
	//	rand() % MOD, rand() % MOD, rand() % MOD,
	//	0, 0, 0;

	// TODO:當 A 為4x3矩陣時，Ax=b是一個超定方程組，隻能求其最小二乘解

	cout << "A: \n" << A << "\n\n";

	MatrixXd x(3, 1);
	x << rand() % MOD, rand() % MOD, rand() % MOD;

	MatrixXd b = A*x;

	// 求解方程 Ax=b，未知量為x
	// 對系數矩陣A進行SVD，A = U*D*V^t
	JacobiSVD<MatrixXd> svd(A, ComputeFullU | ComputeFullV); // SVD時不僅計算奇異值，還計算完整的U和V

	MatrixXd U_inv = svd.matrixU().transpose(); //正交矩陣的逆就是它的轉置
	MatrixXd V = svd.matrixV(); 
	cout << "U_inv: \n" << U_inv << "\n\n";
	cout << "V: \n" << V << "\n\n";

	// 求對角矩陣的逆，注意由于有0奇異值，對角陣隻對左上角非0子矩陣求逆，其他部分為0
	MatrixXd d = svd.singularValues();
	MatrixXd D_inv(3, 3);
	for (int i = 0; i < d.size(); ++i) {
		// 若奇異值大于0
		if (d(i) > 1e-6)
			D_inv(i, i) = 1.0/d(i);
		else
			break;
	}
	cout << "D_inv: \n" << D_inv << "\n\n";

	// 求A的僞逆。A_inv = V*D_inv*U^t
	MatrixXd A_inv = V * D_inv * U_inv;

	// 求解x，并對比與真值是否相同。
	// 當 A 的秩為2時，Ax=b是一個欠定方程組，x的解有無數個，求解出來的可能與真值不同，但是Ax的結果仍然為b
	MatrixXd x_solve = A_inv * b;
	cout << "x: \n" << x << "\n\n";
	cout << "x_solve: \n" << x_solve << "\n\n";

	// 對比b，b應該總是相同的
	cout << "b: \n" << b << "\n\n";
	cout << "A*x_solve: \n" << A*x_solve << "\n\n";

	return 0;
}

利用SVD求解齊次線性方程組Ax=0

下面代碼僅依賴于eigen庫，在VS2015上調試通過。

#include <iostream>
#include <ctime>
#include "Eigen/Dense"

using namespace std;
using namespace Eigen;

const int MOD = 100;
int main() {

	// 構造輸入資料
	uint32_t seed = time(0);
	srand(seed); // 随機初始化種子
	MatrixXd A(3, 3);
	 當 A 是滿秩矩陣時，隻有零解
	//A << rand() % MOD, rand() % MOD, rand() % MOD,
	//	rand() % MOD, rand() % MOD, rand() % MOD, 
	//	rand() % MOD, rand() % MOD, rand() % MOD;

	// 當 A 是虧秩矩陣時，有非零解
	// 注釋這一行去檢視一下結果
	A << rand() % MOD, rand() % MOD, rand() % MOD,
		rand() % MOD, rand() % MOD, rand() % MOD,
		0, 0, 0;


	cout << "A: \n" << A << "\n\n";

	// 求解方程 Ax=b，未知量為x
	// 對系數矩陣A進行SVD，A = U*D*V^t
	JacobiSVD<MatrixXd> svd(A, ComputeFullV); // SVD時不僅計算奇異值，還計算完整的V，不需要計算U，可以減少SVD的計算量
	MatrixXd V = svd.matrixV(); 
	cout << "V: \n" << V << "\n\n";

	// 求解x。Ax=0方程組有無窮多解，V的最後一列列向量是一個非零解，且|x|=1
	MatrixXd x_solve = V.col(V.cols()-1);
	cout << "x_solve: \n" << x_solve << "\n\n";

	// 當 A 是滿秩矩陣時，隻有零解，沒有非零解，則A*x_solve 與 0 相差比較大
	// 當 A 是虧秩矩陣時，有非零解，則 A*x_solve = 0
	cout << "A*x_solve: \n" << A*x_solve << "\n\n";

	return 0;
}

附錄——一些概念

齊次線性方程組：Ax=0
非齊次線性方程組：Ax=b
超定方程 | 超定方程組是指方程個數大于未知量個數的方程組。對于方程組Ra=y，R為n×m矩陣，如果R列滿秩，且n>m。則方程組沒有精确解，此時稱方程組為超定方程組。超定方程一般是不存在解的沖突方程，解超定方程組非常普遍。比較常用的方法是最小二乘法。
欠定方程 | 方程個數小于未知量個數的方程組。
正交矩陣一定是方陣

參考資料

《計算機視覺中的數學方法》——吳朝福 | 本文主要參考資料
線性方程組的幾種求解方法
求解線性方程組解的方法？
Solving linear equation systems | 列舉了用各種方法解線性方程組
百度百科——超定方程組
線性方程組求解 | 介紹了哪些線性方程組好解，後續的各種方法都是想辦法把方程組轉化為這些形式
SVD分解為什麼是最好的？QR分解和SVD比較？LU呢？SVD并行算法可行麼？ - 知乎 | QR分解為什麼能拿來求最小二乘，LU分解、QR分解、SVD分解的數值穩定性和時間複雜度

【動手學MVG】矩陣分解與線性方程組的關系，求解線性方程組實戰代碼

文章目錄

本文解決的問題

QR分解

投影矩陣分解得到内外參

SVD分解

最小二乘問題

最小二乘問題的求解方法

(列)滿秩最小二乘問題

正規化方法

QR分解方法

SVD 分解方法

虧秩最小二乘問題

齊次最小二乘問題

數值穩定性問題

Ax=0與Ax=b的選擇

代碼實作

利用SVD求解非齊次線性方程組Ax=b

利用SVD求解齊次線性方程組Ax=0

附錄——一些概念

參考資料

繼續閱讀

[矩陣的三角分解系列五] 三角分解中的行列變換簡介行變換分解例子引用

28個考研數學易錯點

【線性代數公開課MIT Linear Algebra】第十七課正交基和正交矩陣1. 标準正交基與正交矩陣2. 格拉姆-施密特正交化 Graham-Schmidt3. QR分解

機器學習筆記——基于奇異值分解（SVD）的圖像壓縮（PIL）

線性代數筆記16：了解PCA基本思想直覺了解PCA是最佳的仿射變換拟合PCA保留最大方差參考資料

線性代數02 線性方程組的解的情況（矩陣的秩）

矩陣分解(Matrix Factorization)與推薦系統

解線性方程組的直接方法--低階稠密矩陣三角分解法--Courant分解學習matlab實作

向量空間及相關定義

奇異值分解(SVD)的了解SVD:Sigular Value Decomposition

矩陣論（三）：矩陣分解—從Schur分解、特征值分解EVD到奇異值分解SVD（上）Schur分解（任意方陣）

矩陣的特征值分解與奇異值分解的幾何意義

線性代數 | (5) 線性方程組

拓端tecdat|python機器學習：推薦系統實作（以矩陣分解來協同過濾）1. 使用者和産品的潛在特征2. 編寫推薦系統3. 矩陣分解工作原理

SVD了解和其在PCA,LSI的應用

深度學習推薦系統-筆記08：傳統推薦算法發展彙總