利用 Bareiss 算法判別正定矩陣

本文的閱讀等級:中級

A 為一 n\times n 階實對稱矩陣且任一 n 維非零實向量 \mathbf{x} 滿足 \mathbf{x}^T A\mathbf{x}>0,則 A 稱為正定矩陣 (見“特殊矩陣 (6):正定矩陣”)。如果 A 為一 Hermitian 複矩陣,將上述實向量改成複向量,轉置 (\cdot)^T 替換為共軛轉置 (\cdot)^\ast 即可。正定矩陣有多種判別方式,當下列任一條件成立時,A 是正定矩陣 (見“正定矩陣的性質與判別方法”):

  • A 的特徵值皆為正數;
  • A 的軸元 (pivot) 皆為正數;
  • A 的領先主子陣 (leading principal submatrix) 的行列式皆為正數;
  • A 可分解為 A=B^\ast B,其中 B 是一可逆矩陣。

本文介紹一個基於領先主子陣的行列式的正定矩陣判別法,稱為 Bareiss 算法[1]。在開始討論前,我們先說明本文所使用的子陣表達記號。令 A_{R,C} 代表 n\times n 階矩陣 A 的一個 \vert R\vert\times\vert C\vert 階子陣,其中 R\subseteq\{1,\ldots,n\} 是列 (row) 指標集合,C\subseteq\{1,\ldots,n\} 是行 (column) 指標集合,\vert R\vert\vert C\vert 分別表示列指標集合的基數 (cardinal number) 和行指標集合的基數。見下例,

A=\left[\!\!\begin{array}{rrrr}  5&-1&3&-1\\  -1&2&-2&-1\\  3&-2&3&1\\  -1&-1&1&7  \end{array}\!\!\right]

包含的四個領先主子陣依序是

\displaystyle\begin{aligned} A_{\{1\},\{1\}}&=\begin{bmatrix} 5 \end{bmatrix}\\ A_{\{1,2\},\{1,2\}}&=\left[\!\!\begin{array}{rr}  5&-1\\  -1&2  \end{array}\!\!\right]\\ A_{\{1,2,3\},\{1,2,3\}}&=\left[\!\!\begin{array}{rrr}  5&-1&3\\  -1&2&-2\\  3&-2&3 \end{array}\!\!\right]\\ A_{\{1,2,3,4\},\{1,2,3,4\}}&=A.\end{aligned}

Bareiss 算法將 \det A_{\{1,\ldots,k\},\{1,\ldots,k\}}k=1,2,3,4,的計算整合成一套演算程序:以二階行列式運算實現高斯消去法,化簡 A 至上三角矩陣

\displaystyle \begin{bmatrix} \det A_{\{1\},\{1\}}&\ast&\ast&\ast\\ 0&\det A_{\{1,2\},\{1,2\}}&\ast&\ast\\ 0&0&\det A_{\{1,2,3\},\{1,2,3\}}&\ast\\ 0&0&0&\det A_{\{1,2,3,4\},\{1,2,3,4\}} \end{bmatrix}

其中主對角元即為 A 的領先主子陣的行列式。

 
Bareiss 算法建立在 Sylvester 恆等式之上。令 A=[a_{ij}] 為一 n\times n 階矩陣,以分塊矩陣表示為

\displaystyle A=\begin{bmatrix} B&C\\ D&E \end{bmatrix}

其中 B=A_{\{1,\ldots,k\},\{1,\ldots,k\}}k\times k 階領先主子陣。假設 B 是可逆矩陣。使用基本列運算將方陣 A 的左下分塊 D 消去,算式如下:

\displaystyle \begin{bmatrix} I_k&0\\ -DB^{-1}&I_{n-k} \end{bmatrix}\begin{bmatrix} B&C\\ D&E \end{bmatrix}=\begin{bmatrix} B&C\\ 0&F \end{bmatrix}

其中 (n-k)\times(n-k) 階分塊 F=E-DB^{-1}C 稱為 B 的 Schur 互補 (complement)。使用矩陣乘積的行列式可乘公式和分塊三角矩陣的行列式性質 (見“分塊矩陣的行列式”),可得

\displaystyle \det A=(\det B)(\det F)

稱之為 Schur 恆等式。直白地說,一方陣的行列式等於任一可逆領先主子陣與其 Schur 互補的行列式乘積。為簡化記號,以下令 Schur 互補 F=[f_{ij}] 的列指標 i 和行指標 j 為該元於 A 中所在位置,即

\displaystyle F=\begin{bmatrix} f_{k+1,k+1}&\cdots&f_{k+1,n}\\ \vdots&\ddots&\vdots\\ f_{n,k+1}&\cdots&f_{nn} \end{bmatrix}

按同樣方式,CDE 也沿用方陣 A 的列行指標。在這個規範下,Schur 互補 F 的每一元可表示為

\displaystyle f_{ij}=(E-DB^{-1}C)_{ij}=(E)_{ij}-(DB^{-1}C)_{ij},~~k+1\le i,j\le n

 
對於 k+1\le i,j\le n,考慮包含 k\times kB 分塊的 (k+1)\times(k+1) 階共邊矩陣 (bordered matrix)

\displaystyle A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}}=\begin{bmatrix} B&C'\\ D'&E' \end{bmatrix}=\begin{bmatrix} a_{11}&\cdots&a_{1k}&a_{1j}\\ \vdots&\ddots&\vdots&\vdots\\ a_{k1}&\cdots&a_{kk}&a_{kj}\\ a_{i1}&\cdots&a_{ik}&a_{ij} \end{bmatrix}

其中 C'=(a_{1j},\ldots,a_{kj})^T 是一 k 維行向量,D'=(a_{i1},\ldots,a_{ik}) 是一 k 維列向量,E'=a_{ij}。因為 \displaystyle F'=E'-D'B^{-1}C'=(E)_{ij}-(DB^{-1}C)_{ij}=f_{ij},共邊矩陣 A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}} 的 Schur 恆等式如下:

\displaystyle \det A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}}=(\det B)(\det F')=(\det B)f_{ij}

定義 (n-k)\times(n-k) 階矩陣 A^{(k)}=\begin{bmatrix}\alpha^{(k)}_{ij}\end{bmatrix}k+1\le i,j\le n,其中

\displaystyle \alpha^{(k)}_{ij}=\det A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}}=\begin{vmatrix} a_{11}&\cdots&a_{1k}&a_{1j}\\ \vdots&\ddots&\vdots&\vdots\\ a_{k1}&\cdots&a_{kk}&a_{kj}\\ a_{i1}&\cdots&a_{ik}&a_{ij} \end{vmatrix}

合併上面二式,可得 \alpha^{(k)}_{ij}=(\det B)f_{ij},矩陣表達式為 A^{(k)}=(\det B)F。使用 Schur 恆等式,推得

\displaystyle\begin{aligned} \det A^{(k)}&=\det((\det B)F)=(\det B)^{n-k}\det F\\ &=(\det B)^{n-k}\frac{\det A}{\det B}=(\det B)^{n-k-1}\det A. \end{aligned}

這個結果稱為 Sylvester 恆等式。以上推論雖假設 B 為可逆矩陣,但 Sylvester 恆等式對於不可逆矩陣 B 依然成立 (使用連續論證法即可證明,見“連續論證法”)。值得注意的是

\displaystyle \det B=\begin{vmatrix} a_{11}&\cdots&a_{1,k-1}&a_{1k}\\ \vdots&\ddots&\vdots&\vdots\\ a_{k-1,1}&\cdots&a_{k-1,k-1}&a_{k-1,k}\\ a_{k1}&\cdots&a_{k,k-1}&a_{kk} \end{vmatrix}=\alpha^{(k-1)}_{kk}

故 Sylvester 恆等式可以明確地表示為

\displaystyle \left(\alpha_{kk}^{(k-1)}\right)^{n-k-1}\det A=\begin{vmatrix} \alpha^{(k)}_{k+1,k+1}&\cdots&\alpha^{(k)}_{k+1,n}\\[0.3em] \vdots&\ddots&\vdots\\[0.3em] \alpha^{(k)}_{n,k+1}&\cdots&\alpha^{(k)}_{nn} \end{vmatrix}

由於 \alpha^{(k)}_{ij}=\det A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}} 是一 (k+1) 階行列式,套用 Sylvester 恆等式,將 A 替換為 A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}},可得下列表達式:

\displaystyle \alpha_{ij}^{(k)}=\frac{1}{\left(a_{pp}^{(p-1)}\right)^{k-p}}\begin{vmatrix} \alpha_{p+1,p+1}^{(p)}&\cdots&\alpha_{p+1,k}^{(p)}&\alpha_{p+1,j}^{(p)}\\[0.3em] \vdots&\ddots&\vdots&\vdots\\[0.3em] \alpha_{k,p+1}^{(p)}&\cdots&\alpha_{kk}^{(p)}&\alpha_{kj}^{(p)}\\[0.5em] \alpha_{i,p+1}^{(p)}&\cdots&\alpha_{ik}^{(p)}&\alpha_{ij}^{(p)} \end{vmatrix},~~1\le p\le k-1

如果 A 的所有元都是整數,子陣 A_{\{1,\ldots,k,i\},\{1,\ldots,k,j\}} 的行列式 \alpha_{ij}^{(k)} 必定是整數,也就是說,上式 \left(a_{pp}^{(p-1)}\right)^{k-p} 整除等號右邊的行列式。

 
剩下來的問題是設計一個算法求出 Ak\times k 階領先主子陣 A_{\{1,\ldots,k\},\{1,\ldots,k\}} 的行列式,即 \alpha_{kk}^{(k-1)}1\le k\le n。使用 Sylvester 恆等式導出的 \alpha_{ij}^{(k)} 表達式,設 p=k-1,可得 \alpha_{ij}^{(k)} 的二階行列式算式:

\displaystyle \alpha_{ij}^{(k)}=\frac{1}{\alpha_{k-1,k-1}^{(k-2)}}\begin{vmatrix} \alpha_{kk}^{(k-1)}&\alpha_{kj}^{(k-1)}\\[0.5em] \alpha_{ik}^{(k-1)}&\alpha_{ij}^{(k-1)} \end{vmatrix}

根據這個遞迴關係式,設定初始值,我們得到下列算法。

Bareiss 算法


初始化 \alpha_{00}^{(-1)}=1\alpha_{ij}^{(0)}=a_{ij}i,j=1,\ldots,n
對於 k=1,2,\ldots,n-1
    對於 i=k+1,\ldots,n
        對於 j=k+1,\ldots,n

\displaystyle \alpha_{ij}^{(k)}\leftarrow\frac{1}{\alpha_{k-1,k-1}^{(k-2)}}\begin{vmatrix} \alpha_{kk}^{(k-1)}&\alpha_{kj}^{(k-1)}\\[0.5em] \alpha_{ik}^{(k-1)}&\alpha_{ij}^{(k-1)} \end{vmatrix}


 
如何利用 Bareiss 算法判別正定矩陣?我們需要一個易於操作的簿記方法。因為 k 值隨演算程序遞增,且 k+1\le i,j\le n,我們可以將算得的 \alpha_{ij}^{(k)} 取代 \alpha_{ij}^{(k-1)}。上例中,當 k=1,計算過程如下:

\displaystyle \left[\!\!\begin{array}{rrrr}  5&-1&3&-1\\  -1&2&-2&-1\\  3&-2&3&1\\  -1&-1&1&7  \end{array}\!\!\right]\xrightarrow[]{k=1}\left[\!\!\begin{array}{cccc} 5&\ast&\ast&\ast\\ 0&\left|\!\!\begin{array}{rr} 5&-1\\ -1&2 \end{array}\!\!\right| &\left|\!\!\begin{array}{rr} 5&3\\ -1&-2 \end{array}\!\!\right|&\left|\!\!\begin{array}{rc} 5&-1\\ -1&-1 \end{array}\!\!\right|\\[0.8em] 0&\left|\!\!\begin{array}{cr} 5&-1\\ 3&-2 \end{array}\!\!\right| &\left|\!\!\begin{array}{rc} 5&3\\ 3&3 \end{array}\!\!\right|&\left|\!\!\begin{array}{rr} 5&-1\\ 3&1 \end{array}\!\!\right|\\[0.8em] 0&\left|\!\!\begin{array}{rc} 5&-1\\ -1&-1 \end{array}\!\!\right| &\left|\!\!\begin{array}{rc} 5&3\\ -1&1 \end{array}\!\!\right|&\left|\!\!\begin{array}{rr} 5&-1\\ -1&7 \end{array}\!\!\right| \end{array}\!\!\right]=\left[\!\!\begin{array}{crrr}  5&\ast&\ast&\ast\\ 0&9&-7&-6\\ 0&-7&6&8\\ 0&-6&8&34  \end{array}\!\!\right]

符號 \ast 表示該元不影響下一階段的計算,故可忽略。另外,軸元 5 底下的所有元都是 0,原因在於 Bareiss 算法其實是高斯消去法的一種變形 (見“高斯消去法”)。對於每一列,Bareiss 算法包含三個基本列運算步驟:伸縮、取代,並再次伸縮。以第二列為例,先令第二列通乘 5,將第二列減去 -1 通乘第一列的結果,最後第二列再通除 \alpha_{00}^{(-1)}=1,如下所示:

\displaystyle\begin{aligned} \left[\!\!\begin{array}{rrrr}  5&-1&3&-1\\  -1&2&-2&-1\\  3&-2&3&1\\  -1&-1&1&7  \end{array}\!\!\right]&\to\left[\!\!\begin{array}{cccc}  5&-1&3&-1\\  5\cdot(-1)&5\cdot 2&5\cdot(-2)&5\cdot(-1)\\  3&-2&3&1\\  -1&-1&1&7  \end{array}\!\!\right]\\ &\to\left[\!\!\begin{array}{rccc}  5&-1&3&-1\\  0&\left|\!\!\begin{array}{rr} 5&-1\\ -1&2 \end{array}\!\!\right| &\left|\!\!\begin{array}{rr} 5&3\\ -1&-2 \end{array}\!\!\right|&\left|\!\!\begin{array}{rc} 5&-1\\ -1&-1 \end{array}\!\!\right|\\  3&-2&3&1\\  -1&-1&1&7  \end{array}\!\!\right]\end{aligned}

除開引進通除運算步驟,Bareiss 算法無異於我們曾經介紹過的一個基於二階行列式的矩陣秩算法 (見“利用行列式計算矩陣秩”)。當 k=2k=3,Bareiss 算法的結果如下:

\displaystyle\begin{aligned} \left[\!\!\begin{array}{crrr}  5&\ast&\ast&\ast\\ 0&9&-7&-6\\ 0&-7&6&8\\ 0&-6&8&34  \end{array}\!\!\right]&\xrightarrow[]{k=2}\left[\!\!\begin{array}{cccc} 5&\ast&\ast&\ast\\ 0&9&\ast&\ast\\ 0&0&\frac{1}{5}\left|\!\!\begin{array}{rr} 9&-7\\ -7&6 \end{array}\!\!\right|& \frac{1}{5}\left|\!\!\begin{array}{rr} 9&-6\\ -7&8 \end{array}\!\!\right|\\[0.8em] 0&0&\frac{1}{5}\left|\!\!\begin{array}{rr} 9&-7\\ -6&8 \end{array}\!\!\right|& \frac{1}{5}\left|\!\!\begin{array}{rr} 9&-6\\ -6&34 \end{array}\!\!\right| \end{array}\!\!\right]\\ &=\left[\!\!\begin{array}{cccr}  5&\ast&\ast&\ast\\ 0&9&\ast&\ast\\ 0&0&1&6\\ 0&0&6&54  \end{array}\!\!\right]\\ &\xrightarrow[]{k=3}\left[\!\!\begin{array}{cccc} 5&\ast&\ast&\ast\\ 0&9&\ast&\ast\\ 0&0&1&\ast\\ 0&0&0&\frac{1}{9}\left|\!\!\begin{array}{cr} 1&6\\ 6&54 \end{array}\!\!\right| \end{array}\!\!\right]\\ &=\left[\!\!\begin{array}{cccr}  5&\ast&\ast&\ast\\ 0&9&\ast&\ast\\ 0&0&1&\ast\\ 0&0&0&2 \end{array}\!\!\right]\end{aligned}

從上三角矩陣的主對角元讀出

\displaystyle \alpha_{11}^{(0)}=5,~~\alpha_{22}^{(1)}=9,~~\alpha_{33}^{(2)}=1,~~\alpha_{44}^{(3)}=2

因此判定 A 是一正定矩陣。

 
最後我們討論關於 Bareiss 算法的幾個問題。第一,演算過程中如發現 \alpha_{kk}^{(k-1)}\le 0,表明 A 不是正定矩陣,即可停止計算。第二,因為 A 是實對稱矩陣,不難驗證 \alpha_{ij}^{(k)}=\alpha_{ji}^{(k)}k+1\le i,j\le n,故實際計算量幾可減半。第三,針對正定矩陣判定問題,高斯消去法計算軸元,Bariess 算法計算領先主子陣的行列式。Bariess 算法的計算量是高斯消去法的三倍,Bariess 算法具備甚麼優點?若給定矩陣的所有元都是整數,Bariess 算法必產生整數,而且數值不會無節制地變大 (因為通除步驟抑制了數值增長)。最後一個問題:Bareiss 算法是否適用於計算一般 n\times n 階矩陣的行列式?當然可行,不過前提是所有的除數,即 k\times k 階領先主子陣的行列式 \alpha_{kk}^{(k-1)}1\le k\le n-2,必須不等於零。這個問題與 Dodgson 縮合法的遭遇極類似 (見“Dodgson 縮合法──奇特的行列式運算法”),為了避免麻煩,讀者不妨選用效果相當的 Chiò 算法 (見“Chiò 演算法──另類行列式計算法”)。

 
參考來源:
[1] Erwin Bareiss, Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination, Mathematics of computation, (PDF) vol. 22, 1968, pp 565–578.

Advertisements
本篇發表於 線性代數專欄, 行列式 並標籤為 , , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s