高斯有研究矩陣分解？

你可能會說: 我知道，高斯消元法可以用矩陣乘積的形式來表示，相當(dāng)于對方程組的系數(shù)矩陣作分解。

是可以這么說，但這畢竟是后人的觀點(diǎn)。本文并不是指分解，而是想說說由另一個主題引出的矩陣分解。

高斯消元法是為了求解線性方程組而生的，但是，它也可以拿來計算二次型的標(biāo)準(zhǔn)型，即對稱矩陣對角化。

關(guān)于二次型與標(biāo)準(zhǔn)型的知識，可以閱讀下面這篇。

二次型和矩陣合同原來是這么一回事

高斯那會兒還沒有矩陣的概念，更沒有矩陣的運(yùn)算。但為了讓過程簡潔清晰，我們還是使用矩陣的語言來描述和證明高斯的二次型標(biāo)準(zhǔn)化方法（對角化對稱矩陣）。

本文就是為了將這個過程清晰地展現(xiàn)出來，讓大家分分鐘就能整明白。對，我們的目標(biāo)就是: 只花幾分鐘，知識、方法帶回家。

1高斯與二次型

矩陣分解法是隨著行列式、線性方程組，尤其是雙線性形式和二次型等問題的研究而逐漸顯現(xiàn)出來的。

拉格朗日、高斯和雅可比等可以說是做了一些早期工作。

大概是為了計算二次型的極值，拉格朗日在 1759 年提出了所謂的配方法，用現(xiàn)在的話說就是構(gòu)造以三角矩陣表示的線性變換來實(shí)現(xiàn)對二次型的標(biāo)準(zhǔn)化處理。

例如，可以用如下方式來將一個二元二次型轉(zhuǎn)化為標(biāo)準(zhǔn)型，

\begin{aligned} \mathbf{x}^{\top} \mathbf{A} \mathbf{x} &=a_{11} x_{1}^{2}+2 a_{12} x_{1} x_{2}+a_{22} x_{2}^{2} \\ &=a_{11}\left(x_{1}+\frac{a_{12}}{a_{11}} x_{2}\right)^{2}+\frac{a_{11} a_{22}-a_{12}^{2}}{a_{11}} x_{2}^{2} \\ &=D_{1} y_{1}^{2}+\frac{D_{2}}{D_{1}} y_{2}^{2} \end{aligned}

其中，為系數(shù)矩陣的階順序主子式，以及

\left[\begin{array}{l} y_{1} \\ y_{2} \end{array}\right]=\left[\begin{array}{cc} 1 & \frac{a_{12}}{a_{11}} \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right].

好家伙，又發(fā)現(xiàn)行列式的身影了。

拉格朗日并沒有處理元一般情況，后來雅可比在 1834 年左右推廣了上面方法，有些資料將其稱為雅可比公式。

本文主要看高斯的工作，他在 1823 年撰寫的手稿中使用他早在 1801 年就使用的消元法來實(shí)現(xiàn)將一般二次型轉(zhuǎn)化為標(biāo)準(zhǔn)型的任務(wù)（也有一說認(rèn)為，高斯是在拉格朗日方法的啟發(fā)下實(shí)現(xiàn)最小二乘法，這確實(shí)也是一個二次型的極值問題）。

這里我們用高斯的符號來書寫。具體而言，可以將函數(shù) （關(guān)于的一個二次型）簡化為以下形式：

\frac{u^{0} u^{0}}{\mathcal{A}^{0}}+\frac{u^{\prime} u^{\prime}}{\mathcal{B}^{\prime}}+\frac{u^{\prime \prime} u^{\prime \prime}}{\mathcal{C}^{\prime \prime}}+\frac{u^{\prime \prime \prime} u^{\prime \prime \prime}}{\mathcal{D}^{\prime \prime \prime}}+\mathrm{etc.}+M

其中，除數(shù) 等是常數(shù)，而等是等的線性函數(shù)。但是，第二個函數(shù) 獨(dú)立于；第三個獨(dú)立于和；第四個獨(dú)立于和，依此類推。最后一個函數(shù) 僅取決于最后一個未知數(shù)。此外，系數(shù) 在中分別乘以。

從中我們可以看到高斯這里的方法將二次型的矩陣分解為乘積，其中是對角矩陣，而是與具有相同對角元素的上三角矩陣。

高斯這里的等函數(shù)是向量對應(yīng)的元素。

值得注意的是，說它們是矩陣分解，肯定是以后人的眼光來看待這些工作。

2推導(dǎo)

接下來我們用矩陣運(yùn)算那一套來證明高斯的對角化方法。

首先，我們知道通過左乘一系列初等矩陣的乘積，可以將矩陣轉(zhuǎn)化為一個上三角矩陣，

\mathbf{PA} = \mathbf{U} \qquad (1)

由于矩陣是對稱的，再右乘矩陣的轉(zhuǎn)置，相當(dāng)于對矩陣作了列的初等變換，最終可以將矩陣轉(zhuǎn)化為對角矩陣。

\mathbf{PAP}^{\top} = \mathbf{UP}^{\top}=\mathbf{\Sigma} \qquad (2)

于是，我們就得到了矩陣的如下分解，

\mathbf{U}^{\top}\mathbf{\Sigma}^{-1}\mathbf{U} \qquad (3)

注意，這里的矩陣只是上三角矩陣，并不是正交矩陣。

驗(yàn)證

由式得，由式可得，將它們代入式，得

\begin{array}{l} \mathbf{U}^{\top}\mathbf{\Sigma}^{-1}\mathbf{U} &= \mathbf{(PA)}^{\top} (\mathbf{PAP}^{\top})^{-1} \mathbf{PA}\\[0.68em] &= (\mathbf{A}^{\top}\mathbf{P}^{\top})(\mathbf{P}^{\top})^{-1}\mathbf{A}^{-1}\mathbf{P}^{-1}\mathbf{PA}\\[0.68em] &= \mathbf{A}^{\top}\color{red}{\boxed{\mathbf{P}^{\top}(\mathbf{P}^{\top})^{-1}\mathbf{A}^{-1}\mathbf{P}^{-1}\mathbf{PA}}}\\[0.68em] &= \mathbf{A}^{\top}\\[0.68em] &= \mathbf{A} \end{array}

上式中紅色框框里的矩陣剛好消去了，最后再根據(jù)矩陣是對稱的，驗(yàn)證通過。

當(dāng)然，這種方法依賴高斯消元，那么消元法的弊端也可能會遇到，即在消元過程中碰到對角元素為零的情況。

下面，我們對式作進(jìn)一步說明。

我們知道，高斯消元的過程相當(dāng)于對矩陣不停地左乘初等矩陣。

先構(gòu)造如下初等矩陣，

\mathbf{P}_{1} = \left[\begin{array}{cccc} 1 & 0 & \cdots & 0 \\ -\frac{1}{a_{11}} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{array}\right]

對系數(shù)矩陣左乘該初等矩陣，對第 2 行首個元素進(jìn)行消元，

\begin{array}{l} \mathbf{P}_{1} \mathbf{A} &= \left[\begin{array}{cccc} 1 & 0 & \cdots & 0 \\ -\frac{a_{21}}{a_{11}} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{array}\right] \left[\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right]\\[1em] &= \left[\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ \color{red}{0} & \color{red}{a'_{22}} & \color{red}{\cdots} & \color{red}{a'_{2 n}} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right] \end{array}

由于矩陣是對稱的，即。此時如果對上面結(jié)果再右乘，相當(dāng)于對列作初等變換，可得對稱矩陣，

\begin{array}{l} \mathbf{P}_{1} \mathbf{A} \mathbf{P}_{1}^{\top} &= \left[\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ \color{red}{0} & \color{red}{a'_{22}} & \color{red}{\cdots} & \color{red}{a'_{2 n}} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right] \left[\begin{array}{cccc} 1 & -\frac{a_{21}}{a_{11}} & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{array}\right] \\[1em] &= \left[\begin{array}{cccc} a_{11} & \color{red}{0} & \cdots & a_{1 n} \\ \color{red}{0} & \color{red}{a'_{22}} & \color{red}{\cdots} & \color{red}{a'_{2 n}} \\ \vdots & \color{red}{\vdots} & \ddots & \vdots \\ a_{n 1} & \color{red}{a'_{n 2}} & \cdots & a_{n n} \end{array}\right] \end{array}

對第 1 列和第 1 行繼續(xù)消元，得

\begin{array}{l} \mathbf{P}_{n-1} \cdots \mathbf{P}_{1} \mathbf{A} \mathbf{P}_{1}^{\top} \cdots \mathbf{P}_{n-1}^{\top} = \left[\begin{array}{cccc} a_{11} & \color{blue}{0} & \cdots & \color{blue}{0} \\ \color{blue}{0} & \color{blue}{\hat{a}_{22}} & \cdots & \color{blue}{\hat{a}_{2 n}} \\ \vdots & \vdots & \ddots & \vdots \\ \color{blue}{0} & \color{blue}{\hat{a}_{n 2}} & \cdots & \color{blue}{\hat{a}_{n n}} \end{array}\right] \end{array}

接下來除去首行和首列，是不是變成一個對稱矩陣?yán)病?/p>

兩邊繼續(xù)左右開弓，最后得到對角矩陣，

\mathbf{P}_{\mathcal{N}-1} \cdots \mathbf{P}_{1} \mathbf{A} \mathbf{P}_{1}^{\top} \cdots \mathbf{P}_{\mathcal{N}-1}^{\top} =\mathbf{\Sigma}

把所有初等矩陣乘起來得到矩陣，于是得到式??，

\mathbf{PAP}^{\top} =\mathbf{\Sigma}

二次型的標(biāo)準(zhǔn)化

高斯這里的方法本來就是為了對二次型作標(biāo)準(zhǔn)化處理?，F(xiàn)在已經(jīng)有了如下分解，

\mathbf{A} = \mathbf{U}^{\top}\mathbf{\Sigma}^{-1}\mathbf{U}

則對應(yīng)的二次型為，

\begin{array}{l} \mathbf{x}^{\top}\mathbf{A}\mathbf{x} &= \mathbf{x}^{\top}\mathbf{U}^{\top}\mathbf{\Sigma}^{-1}\mathbf{U}\mathbf{x} \\[0.68em] &=\left(\mathbf{Ux}\right)^{\top} \mathbf{\Sigma}^{-1} \left(\mathbf{U}\mathbf{x}\right) \\[0.68em] &= \mathbf{\xi}^{\top} \mathbf{\Sigma}^{-1} \mathbf{\xi} \end{array}

于是找到了變量代換，即，實(shí)現(xiàn)了對二次型的標(biāo)準(zhǔn)化處理。對，這就是高斯論文里做的事情。

程序小實(shí)驗(yàn)

下面我用 NumPy 做個實(shí)驗(yàn)來實(shí)際感受一番。

import?numpy?as?np

#?生成一個對稱矩陣
np.random.seed(9)
a?=?np.random.randint(0,5,(3,3))
A?=?a+a.T
A

array([[8, 1, 3],
       [1, 8, 5],
       [3, 5, 2]])

P1?=?np.eye(3)
P1[1,0]=-1/8
P1

array([[ 1.   ,  0.   ,  0.   ],
       [-0.125,  1.   ,  0.   ],
       [ 0.   ,  0.   ,  1.   ]])

P1@A

array([[8.   , 1.   , 3.   ],
       [0.   , 7.875, 4.625],
       [3.   , 5.   , 2.   ]])

[email protected]

array([[8.   , 0.   , 3.   ],
       [1.   , 7.875, 5.   ],
       [3.   , 4.625, 2.   ]])

用 P1 左右開弓，是不是消元的同時將 A 轉(zhuǎn)化為另一個對稱矩陣?yán)病?/p>

P1@[email protected]

array([[8.   , 0.   , 3.   ],
       [0.   , 7.875, 4.625],
       [3.   , 4.625, 2.   ]])

P2?=?np.eye(3)
P2[2,0]?=?-3/8
P2

array([[ 1.   ,  0.   ,  0.   ],
       [ 0.   ,  1.   ,  0.   ],
       [-0.375,  0.   ,  1.   ]])

P2@A

array([[8.   , 1.   , 3.   ],
       [1.   , 8.   , 5.   ],
       [0.   , 4.625, 0.875]])

P2@P1@A

array([[8.   , 1.   , 3.   ],
       [0.   , 7.875, 4.625],
       [0.   , 4.625, 0.875]])

[email protected]@P2.T

array([[8.   , 0.   , 0.   ],
       [1.   , 7.875, 4.625],
       [3.   , 4.625, 0.875]])

A1?=?P2@P1@[email protected]@P2.T
A1

array([[8.   , 0.   , 0.   ],
       [0.   , 7.875, 4.625],
       [0.   , 4.625, 0.875]])

P3?=?np.eye(3)
P3[2,1]?=?-A1[2,1]/A1[1,1]
P3

array([[ 1.        ,  0.        ,  0.        ],
       [ 0.        ,  1.        ,  0.        ],
       [ 0.        , -0.58730159,  1.        ]])

P3@A1

array([[ 8.        ,  0.        ,  0.        ],
       [ 0.        ,  7.875     ,  4.625     ],
       [ 0.        ,  0.        , -1.84126984]])

D?=?P3@[email protected]
D

array([[ 8.        ,  0.        ,  0.        ],
       [ 0.        ,  7.875     ,  0.        ],
       [ 0.        ,  0.        , -1.84126984]])

U?=?P3@P2@P1@A
U

array([[ 8.        ,  1.        ,  3.        ],
       [ 0.        ,  7.875     ,  4.625     ],
       [ 0.        ,  0.        , -1.84126984]])

D_inv?=?np.linalg.inv(D)

U.T@D_inv@U

array([[8., 1., 3.],
       [1., 8., 5.],
       [3., 5., 2.]])

array([[8, 1, 3],
       [1, 8, 5],
       [3, 5, 2]])

注意，這里的上三角矩陣并不是正交矩陣。U.T@U

array([[64.        ,  8.        , 24.        ],
       [ 8.        , 63.015625  , 39.421875  ],
       [24.        , 39.421875  , 33.78089963]])

與其他分解的關(guān)系

我們知道，對稱矩陣可以通過它的特征值和特征向量來對角化，與本文介紹的這種分解方法是不同的。它們的結(jié)果不同，需要的計算量也是不同的。

另外，這里的分解跟后來所謂的 LDL 分解十分相似，僅僅是上三角矩陣的對角元素和中間的對角矩陣取法有一點(diǎn)區(qū)別?？磥磉@個方法又可以歸功于高斯啊。但實(shí)際看貌似并沒有這么命名，不清楚這里面的原因，或許高斯這篇論文并沒有引起多少人關(guān)注吧。

下面的代碼展示了對上面矩陣作 LDL 分解后的結(jié)果。

from?scipy?import?linalg?as?sl

LD?=?sl.ldl(A,lower=True)

LD[0]

array([[1.        , 0.        , 0.        ],
       [0.125     , 1.        , 0.        ],
       [0.375     , 0.58730159, 1.        ]])

LD[1]

array([[ 8.        ,  0.        ,  0.        ],
       [ 0.        ,  7.875     ,  0.        ],
       [ 0.        ,  0.        , -1.84126984]])

LD[0]@LD[1]@LD[0].T

array([[8., 1., 3.],
       [1., 8., 5.],
       [3., 5., 2.]])

3小結(jié)

矩陣分解至少有兩個方面的意義，

理論分析。矩陣分解可以用于進(jìn)一步分析其他問題。這里更注重分解的意義以及在理論推導(dǎo)中的應(yīng)用。
數(shù)值計算。一次分解，結(jié)果可以重復(fù)利用，從而大大降低計算量。這里更注重分解的效率與數(shù)值穩(wěn)定性等問題。

除了批量解線性方程組的例子，我們再來看一個其他例子。

很多時候需要計算一個矩陣的次冪，如

\mathbf{A}^{100}= \underbrace{\mathbf{A} \mathbf{A} \cdots \mathbf{A} \mathbf{A}}_{100 個\; \mathbf{A} \; 相乘}

直接計算顯然不劃算，而且還可能會帶來很大誤差。但如果可以對角化，那么有?。于是，

\begin{array}{l} \mathbf{A}^{100}&=\left(\mathbf{P} \mathbf{\Sigma} \mathbf{P}^{-1}\right)\left(\mathbf{P} \mathbf{\Sigma} \mathbf{P}^{-1}\right) \cdots \left(\mathbf{P} \mathbf{\Sigma} \mathbf{P}^{-1}\right)\left(\mathbf{P} \mathbf{\Sigma} \mathbf{P}^{-1}\right)\\[0.618em] &=\mathbf{P} \mathbf{\Sigma}^{100} \mathbf{P}^{-1} \end{array}

最后只需要計算對角陣的 n 次冪，計算量大大降低了。當(dāng)然，本文中的分解結(jié)果中，是，而不是，因此并不適合上面的應(yīng)用。

本文介紹了高斯利用消元法實(shí)現(xiàn)了二次型的標(biāo)準(zhǔn)化處理，從矩陣的角度看，恰恰就是對稱矩陣的對角化以及分解。高斯這個分解方法的意義或許并不是很大，但可以說是非常早期的分解方法，對后續(xù)方法有一定啟發(fā)作用。

矩陣分解對于矩陣的應(yīng)用十分重要，本篇只是引子，更多精彩在后頭。

相關(guān)閱讀

矩陣前傳 -?克萊姆沒能證明的法則被他兩行搞定

矩陣前傳 -?矩陣之父 Sylvester 為什么提出 Matrix

矩陣前傳 -?柯西-比內(nèi)公式及其用初等矩陣的證明

二次型和矩陣合同原來是這么一回事

拉格朗日乘子法的來歷與直觀解釋

矩陣特征值是這么來的，以及有趣的蓋爾圓

萬能的 SVD 分解是哪位牛人提出來的？

矩陣分解術(shù)，不得不從高斯說起

1高斯與二次型

2推導(dǎo)

驗(yàn)證

二次型的標(biāo)準(zhǔn)化

程序小實(shí)驗(yàn)

與其他分解的關(guān)系

3小結(jié)

矩陣和線性代數(shù)原來是這么來的

概率論原來可以這樣優(yōu)雅地入門

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之向量范數(shù)

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之矩陣范數(shù)

拉格朗日乘子法的來歷與直觀解釋

矩陣分解術(shù)，不得不從高斯說起

1高斯與二次型

2推導(dǎo)

驗(yàn)證

二次型的標(biāo)準(zhǔn)化

程序小實(shí)驗(yàn)

與其他分解的關(guān)系

3小結(jié)

矩陣和線性代數(shù)原來是這么來的

概率論原來可以這樣優(yōu)雅地入門

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之 向量范數(shù)

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之 矩陣范數(shù)

拉格朗日乘子法的來歷與直觀解釋

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之向量范數(shù)

機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ) 之矩陣范數(shù)