矩陣的數值秩

本文的閱讀等級:高級

1879年,德國數學家弗羅貝尼烏斯 (Ferdinand Georg Frobenius) 提出秩 (rank) 作為矩陣所攜帶訊息量的一種測度。一矩陣 A 的秩,記作 \hbox{rank}A,定義為最大可逆子陣的階數,即最大非零餘子式 (minor) 的行列式階數,因此也稱為行列式秩。近代線性代數已捨棄行列式定義,改用較富幾何意義的向量空間定義:秩是矩陣的行空間 (column space) 維數,同樣也是列空間 (row space) 維數,換句話說,秩是最大的線性獨立行 (列) 向量數 (見“你不能不知道的矩陣秩”)。一般基礎線性代數教科書多以高斯消去法揭示矩陣秩 (見“利用行列式計算矩陣秩”),但實際數值計算卻不採用高斯消去法,原因在於 \text{rank}A 並非 A 的連續函數。舉例來說,A=\begin{bmatrix}  1&0\\  0&\epsilon  \end{bmatrix}。若 \epsilon\neq 0,則 \text{rank}A=2。若 \epsilon=0,則 \text{rank}A=1。這個例子顯示極微小的擾動量便足以造成矩陣秩的跳躍。由於計算機的浮點運算難免引入誤差,我們必須審慎處理矩陣秩的數值計算問題。本文介紹奇異值分解在矩陣秩的擾動分析以及數值計算的應用,原因無他,奇異值分解是現今最可靠的數值秩算法 (另一個有效的算法是計算成本較低的 QR 分解[1])。

 
以下考慮實矩陣。令 m\times n 階矩陣 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\\  0&0  \end{bmatrix} 是一 m\times n 階矩陣,分塊 D=\text{diag}(\sigma_1,\ldots,\sigma_r) 由非零奇異值 \sigma_1\ge\cdots\ge\sigma_r>0 構成 (見“奇異值分解 (SVD)”)。因為可逆矩陣 UV^T 不改變 A 的秩,\text{rank}A 等於非零奇異值總數,即 \text{rank}A=r

 
秩是矩陣的真實尺寸的測度,矩陣範數則是矩陣的數值大小的度量。對於一 m\times n 階矩陣 A,2-範數定義為

\displaystyle  \Vert A\Vert_2=\max_{\Vert\mathbf{x}\Vert=1}\Vert A\mathbf{x}\Vert

我們可以證明矩陣的2-範數等於最大奇異值,即 \Vert A\Vert_2=\max_{1\le i\le r}\sigma_i=\sigma_1 (見“矩陣範數”)。下面的矩陣近似定理是必備的分析工具 (見“SVD 於矩陣近似的應用”)。

 
定理一:令 ABm\times n 階矩陣。若 r=\text{rank}A\sigma_1\ge\cdots\ge\sigma_r>0A 的非零奇異值,則每一 k<r

\displaystyle  \sigma_{k+1}=\min_{\text{rank}B=k}\Vert A-B\Vert_2

 
我們利用奇異值分解來證明。令 A 的奇異值分解為 A=U\begin{bmatrix}  D&0\\  0&0  \end{bmatrix}V^T,其中 D=\text{diag}(\sigma_1,\ldots,\sigma_r)。定義 D_{k+1}=\text{diag}(\sigma_1,\ldots,\sigma_{k+1}),且 V=\begin{bmatrix}  V_1&V_2  \end{bmatrix},其中 V_1n\times(k+1) 階分塊,V_2n\times(n-k-1) 階分塊。假設 \text{rank}B=k,矩陣乘法不增大秩,故 \text{rank}(BV_1)\le\text{rank}B=k。根據秩─零度定理,\dim N(BV_1)=k+1-\text{rank}(BV_1)\ge 1,即知存在一 (k+1) 維向量 \mathbf{x} 滿足 BV_1\mathbf{x}=\mathbf{0}\Vert\mathbf{x}\Vert=1。利用正交矩陣性質

\displaystyle  V^TV=\begin{bmatrix}  V_1^T\\  V_2^T  \end{bmatrix}\begin{bmatrix}  V_1&V_2  \end{bmatrix}=\begin{bmatrix}  V_1^TV_1&V_1^TV_2\\  V_2^TV_1&V_2^TV_2  \end{bmatrix}=I

可得

\displaystyle\begin{aligned}  AV_1\mathbf{x}&=U\begin{bmatrix}  D&0\\  0&0  \end{bmatrix}\begin{bmatrix}  V_1^T\\  V_2^T  \end{bmatrix}V_1\mathbf{x}=U\begin{bmatrix}  D&0\\  0&0  \end{bmatrix}\begin{bmatrix}  V_1^TV_1\mathbf{x}\\  V_2^TV_1\mathbf{x}  \end{bmatrix}\\  &=U\begin{bmatrix}  D&0\\  0&0  \end{bmatrix}\begin{bmatrix}  I_{k+1}\mathbf{x}\\  \mathbf{0}  \end{bmatrix}=U\begin{bmatrix}  D_{k+1}\mathbf{x}\\  \mathbf{0}  \end{bmatrix}.\end{aligned}

因為 \Vert V_1\mathbf{x}\Vert^2=\mathbf{x}^TV_1^TV_1\mathbf{x}=\Vert\mathbf{x}\Vert^2=1,使用2-範數定義 \Vert A-B\Vert_2=\max_{\Vert\mathbf{z}\Vert=1}\Vert (A-B)\mathbf{z}\Vert,並代入上式,可得

\displaystyle\begin{aligned}  \Vert A-B\Vert_2^2&\ge\Vert(A-B)V_1\mathbf{x}\Vert^2=\Vert AV_1\mathbf{x}\Vert^2=\left \|U\begin{bmatrix}  D_{k+1}\mathbf{x}\\  \mathbf{0}  \end{bmatrix} \right \|^2\\  &=\Vert{D}_{k+1}\mathbf{x}\Vert^2=\sum_{i=1}^{k+1}\sigma_i^2x_i^2\ge\sigma_{k+1}^2\sum_{i=1}^{k+1}x_i^2=\sigma_{k+1}^2.  \end{aligned}

上面使用了 \Vert U\mathbf{y}\Vert=\Vert\mathbf{y}\Vert。若 B=U\begin{bmatrix}  D_k&0\\  0&0  \end{bmatrix}V^T,其中 D_k=\text{diag}(\sigma_1,\ldots,\sigma_k),則 \Vert A-B\Vert_2 等於 A-B 的最大奇異值 \sigma_{k+1},即證明所求。

 
接著討論矩陣秩的擾動分析。若一矩陣摻雜微小擾動量,矩陣秩將會發生甚麼變化?定理二闡明一個奇特的現象:如果擾動量足夠小,矩陣秩只可能增大或不變,但絕不會減小。

 
定理二:令 AEm\times n 階矩陣。若 r=\text{rank}A\sigma_1\ge\cdots\ge\sigma_r>0A 的非零奇異值,且 \Vert E\Vert_2<\sigma_r,則 \text{rank}(A+E)\ge\text{rank}A

 
使用反證法。假設 k=\text{rank}(A+E)<r,使用定理一,可得

\displaystyle  \Vert E\Vert_2=\Vert A-(A+E)\Vert_2\ge\min_{\text{rank}B=k}\Vert A-B\Vert_2=\sigma_{k+1}\ge\sigma_r

上式與奇異值的設定排序矛盾,證得 \text{rank}(A+E)\ge r=\text{rank}A

 
我們進一步分析這個情況。令 X 為一 m\times m 階可逆矩陣且 Y 為一 n\times n 階可逆矩陣使得 A=X\begin{bmatrix}  I_r&0\\  0&0  \end{bmatrix}Y,稱為等價標準型 (見“每週問題 April 20, 2009”)。設 E=X\begin{bmatrix}  E_{11}&E_{12}\\  E_{21}&E_{22}  \end{bmatrix}Y,其中 E_{11}r\times r 階分塊,則

\displaystyle  A+E=X\begin{bmatrix}  I_r+E_{11}&E_{12}\\  E_{21}&E_{22}  \end{bmatrix}Y

\Vert E_{11}\Vert_2<1,則 I_r+E_{11} 是可逆矩陣 (見“Neumann 無窮級數”)。在此情形下,上式的分塊矩陣可分解為

\displaystyle  \begin{bmatrix}  I&0\\  -E_{21}(I+E_{11})^{-1}&I  \end{bmatrix}\begin{bmatrix}  I+E_{11}&E_{12}\\  E_{21}&E_{22}  \end{bmatrix}\begin{bmatrix}  I&-(I+E_{11})^{-1}E_{12}\\  0&I  \end{bmatrix}=\begin{bmatrix}  I+E_{11}&0\\  0&F  \end{bmatrix}

其中 F=E_{22}-E_{21}(I+E_{11})^{-1}E_{12}。可逆矩陣乘法運算不改變矩陣秩,可知

\displaystyle \text{rank}(A+E)=\text{rank}\begin{bmatrix}  I+E_{11}&0\\  0&F  \end{bmatrix}=\text{rank}(I+E_{11})+\text{rank}F=r+\text{rank}F

如果擾動矩陣 E 具有隨機性,當 r<mr<nE_{22}=E_{21}(I+E_{11})^{-1}E_{12} (即 F=0) 發生的機率非常地小,也就是說,我們幾乎可以肯定 \text{rank}(A+E)>\text{rank}A。如果 A 為一不可逆矩陣,微小誤差 E 很容易製造出可逆矩陣 A+E,使用高斯消去法演算矩陣秩或解齊次方程 A\mathbf{x}=\mathbf{0} 注定要以失敗收場。奇異值分解能夠扭轉劣勢嗎?定理三說明擾動 E 的矩陣範數與 AA+E 的奇異值變動量之間的關係,提供了奇異值於計算矩陣秩的理論基礎。

 
定理三:對於 m\times n 階實矩陣 A 和擾動矩陣 Ep=\min\{m,n\},若 \sigma_1\ge\cdots\ge\sigma_p\ge 0A 的奇異值 (包含正數以及零),\varsigma_1\ge\cdots\ge\varsigma_p\ge 0A+E 的奇異值,則

\displaystyle \vert\sigma_k-\varsigma_k\vert\le \Vert E\Vert_2,~~1\le k\le p

 
證明於下。寫出 A 的奇異值分解的向量表達式

\displaystyle  A=\sum_{i=1}^p\sigma_i\mathbf{u}_i\mathbf{v}_i^T

並令

\displaystyle  A_{k}=\sum_{i=1}^k\sigma_i\mathbf{u}_i\mathbf{v}_i^T

使用兩次定理一,以及三角不等式,

\displaystyle\begin{aligned}  \sigma_k&=\min_{\text{rank}B=k-1}\Vert A-B\Vert_2=\min_{\text{rank}B=k-1}\Vert A+E-B-E\Vert_2\\  &\le\min_{\text{rank}B=k-1}\Vert A+E-B\Vert_2+\Vert E\Vert_2=\varsigma_k+\Vert E\Vert_2.  \end{aligned}

另一方面,矩陣2-範數等於最大的奇異值,使用反向三角不等式 \Vert X-Y\Vert\ge\Vert X\Vert-\Vert Y\Vert (見“每週問題 May 17, 2010”),可得

\displaystyle\begin{aligned}  \sigma_k&=\Vert A-A_{k-1}\Vert_2=\Vert A+E-A_{k-1}-E\Vert_2\\  &\ge \Vert A+E-A_{k-1}\Vert_2-\Vert E\Vert_2\\  &\ge \min_{\text{rank}B=k-1}\Vert A+E-B\Vert_2-\Vert E\Vert_2=\varsigma_k-\Vert E\Vert_2.  \end{aligned}

合併以上結果,即證明 \vert \sigma_k-\varsigma_k\vert\le\Vert E\Vert_2

 
假設 \hbox{rank}A=r,則 \sigma_{r+1}=\cdots=\sigma_p=0。對於擾動矩陣 E,定理三表明 \vert\varsigma_k\vert\le\Vert E\Vert_2r+1\le k\le p。基於這個事實,如果 A+E 的奇異值滿足

\displaystyle  \varsigma_1\ge\cdots\ge\varsigma_{\tilde{r}}>\Vert E\Vert_2\ge\varsigma_{\tilde{r}+1}\ge\cdots\varsigma_p

\tilde{r}A 的一個合理數值秩。在多數的情況下,\Vert E\Vert_2 是一未知門檻,但仍可由所採用的算法推導出估計值。一般而言,好的奇異值分解算法有 \Vert E\Vert_2\approx 5\times 10^{-u}\Vert A\Vert_2,其中 u 是浮點運算的小數位數[2]。奇異值分解不僅揭示矩陣的數值秩,同時也提供齊次方程 A\mathbf{x}=\mathbf{0} 的完整解。一旦確定了數值秩 \tilde{r},使用 A+E 的奇異值分解 A+E=\tilde{U}\begin{bmatrix}  \tilde{D}&0\\  0&0  \end{bmatrix}\tilde{V}^T 或表示為 (A+E)\tilde{V}=U\begin{bmatrix}  \tilde{D}&0\\  0&0  \end{bmatrix},其中 \tilde{D}=\text{diag}(\varsigma_1,\ldots,\varsigma_p),捨去小於 \Vert E\Vert_2 的奇異值 \varsigma_{\tilde{r}+1},\ldots,\varsigma_pA\mathbf{x}=\mathbf{0} 的解空間 (A 的零空間) 即由 \tilde{V}=\begin{bmatrix}  \tilde{V}_1&\tilde{V}_2  \end{bmatrix}n\times (n-\tilde{r}) 階分塊 \tilde{V}_2 的行向量擴張而成。

 
參考來源:
[1] 維基百科:RRQR factorization
[2] Carl Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, 2000, pp 422.

廣告
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s