利用 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}$ 。

按同樣方式， $C$ 、 $D$ 和 $E$ 也沿用方陣 $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 k$ 階 $B$ 分塊的 $(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}$ 整除等號右邊的行列式。

剩下來的問題是設計一個算法求出 $A$ 的 $k\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=2$ 和 $k=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.

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

利用 Bareiss 算法判別正定矩陣

Leave a Reply Cancel reply

搜尋(繁體中文或英文)

訊息看板

近期文章

線性代數專欄

其他主題專欄

每週問題

數據充分性問題

其他分類

Recent Comments

近期最多人點閱

分類

Archives

標籤雲

線代線上影音課程

線代學習網站

線代電子書

矩陣計算器

LaTeX

Blogroll

訂閱

閱讀導引

學習資源

專題探究

急救查詢

其他分頁

Meta

網路狀態

Blog Stats

利用 Bareiss 算法判別正定矩陣

Share this:

Like this:

Leave a Reply Cancel reply

搜尋(繁體中文或英文)

訊息看板

近期文章

線性代數專欄

其他主題專欄

每週問題

數據充分性問題

其他分類

Recent Comments

近期最多人點閱

分類

Archives

標籤雲

線代線上影音課程

線代學習網站

線代電子書

矩陣計算器

LaTeX

Blogroll

訂閱