高斯有研究矩阵分解？

你可能会说: 我知道，高斯消元法可以用矩阵乘积的形式来表示，相当于对方程组的系数矩阵作分解。

是可以这么说，但这毕竟是后人的观点。本文并不是指分解，而是想说说由另一个主题引出的矩阵分解。

高斯消元法是为了求解线性方程组而生的，但是，它也可以拿来计算二次型的标准型，即对称矩阵对角化。

关于二次型与标准型的知识，可以阅读下面这篇。

二次型和矩阵合同原来是这么一回事

高斯那会儿还没有矩阵的概念，更没有矩阵的运算。但为了让过程简洁清晰，我们还是使用矩阵的语言来描述和证明高斯的二次型标准化方法（对角化对称矩阵）。

本文就是为了将这个过程清晰地展现出来，让大家分分钟就能整明白。对，我们的目标就是: 只花几分钟，知识、方法带回家。

1高斯与二次型

矩阵分解法是随着行列式、线性方程组，尤其是双线性形式和二次型等问题的研究而逐渐显现出来的。

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

大概是为了计算二次型的极值，拉格朗日在 1759 年提出了所谓的配方法，用现在的话说就是构造以三角矩阵表示的线性变换来实现对二次型的标准化处理。

例如，可以用如下方式来将一个二元二次型转化为标准型，

\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}

其中，为系数矩阵的阶顺序主子式，以及

\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].

好家伙，又发现行列式的身影了。

拉格朗日并没有处理元一般情况，后来雅可比在 1834 年左右推广了上面方法，有些资料将其称为雅可比公式。

本文主要看高斯的工作，他在 1823 年撰写的手稿中使用他早在 1801 年就使用的消元法来实现将一般二次型转化为标准型的任务（也有一说认为，高斯是在拉格朗日方法的启发下实现最小二乘法，这确实也是一个二次型的极值问题）。

这里我们用高斯的符号来书写。具体而言，可以将函数（关于的一个二次型）简化为以下形式：

\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

其中，除数等是常数，而等是等的线性函数。但是，第二个函数独立于；第三个独立于和；第四个独立于和，依此类推。最后一个函数仅取决于最后一个未知数。此外，系数在中分别乘以。

从中我们可以看到高斯这里的方法将二次型的矩阵分解为乘积，其中是对角矩阵，而是与具有相同对角元素的上三角矩阵。

高斯这里的等函数是向量对应的元素。

值得注意的是，说它们是矩阵分解，肯定是以后人的眼光来看待这些工作。

2推导

接下来我们用矩阵运算那一套来证明高斯的对角化方法。

首先，我们知道通过左乘一系列初等矩阵的乘积，可以将矩阵转化为一个上三角矩阵，

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

由于矩阵是对称的，再右乘矩阵的转置，相当于对矩阵作了列的初等变换，最终可以将矩阵转化为对角矩阵。

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

于是，我们就得到了矩阵的如下分解，

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

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

验证

由式得，由式可得，将它们代入式，得

\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}

上式中红色框框里的矩阵刚好消去了，最后再根据矩阵是对称的，验证通过。

当然，这种方法依赖高斯消元，那么消元法的弊端也可能会遇到，即在消元过程中碰到对角元素为零的情况。

下面，我们对式作进一步说明。

我们知道，高斯消元的过程相当于对矩阵不停地左乘初等矩阵。

先构造如下初等矩阵，

\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]

对系数矩阵左乘该初等矩阵，对第 2 行首个元素进行消元，

\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}

由于矩阵是对称的，即。此时如果对上面结果再右乘，相当于对列作初等变换，可得对称矩阵，

\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 行继续消元，得

\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}

接下来除去首行和首列，是不是变成一个对称矩阵啦。

两边继续左右开弓，最后得到对角矩阵，

\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}

二次型的标准化

高斯这里的方法本来就是为了对二次型作标准化处理。现在已经有了如下分解，

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

则对应的二次型为，

\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}

于是找到了变量代换，即，实现了对二次型的标准化处理。对，这就是高斯论文里做的事情。

程序小实验

下面我用 NumPy 做个实验来实际感受一番。

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.   ]])

A@P1.T

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

用 P1 左右开弓，是不是消元的同时将 A 转化为另一个对称矩阵啦。

P1@A@P1.T

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]])

A@P1.T@P2.T

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

A1 = P2@P1@A@P1.T@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@A1@P3.T
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]])

与其他分解的关系

我们知道，对称矩阵可以通过它的特征值和特征向量来对角化，与本文介绍的这种分解方法是不同的。它们的结果不同，需要的计算量也是不同的。

另外，这里的分解跟后来所谓的 LDL 分解十分相似，仅仅是上三角矩阵的对角元素和中间的对角矩阵取法有一点区别。看来这个方法又可以归功于高斯啊。但实际看貌似并没有这么命名，不清楚这里面的原因，或许高斯这篇论文并没有引起多少人关注吧。

下面的代码展示了对上面矩阵作 LDL 分解后的结果。

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小结

矩阵分解至少有两个方面的意义，

理论分析。矩阵分解可以用于进一步分析其他问题。这里更注重分解的意义以及在理论推导中的应用。
数值计算。一次分解，结果可以重复利用，从而大大降低计算量。这里更注重分解的效率与数值稳定性等问题。

除了批量解线性方程组的例子，我们再来看一个其他例子。

很多时候需要计算一个矩阵的次幂，如

\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 次幂，计算量大大降低了。当然，本文中的分解结果中，是，而不是，因此并不适合上面的应用。

本文介绍了高斯利用消元法实现了二次型的标准化处理，从矩阵的角度看，恰恰就是对称矩阵的对角化以及分解。高斯这个分解方法的意义或许并不是很大，但可以说是非常早期的分解方法，对后续方法有一定启发作用。

矩阵分解对于矩阵的应用十分重要，本篇只是引子，更多精彩在后头。

相关阅读

矩阵前传 - 克莱姆没能证明的法则被他两行搞定

矩阵前传 - 矩阵之父 Sylvester 为什么提出 Matrix

矩阵前传 - 柯西-比内公式及其用初等矩阵的证明

二次型和矩阵合同原来是这么一回事

拉格朗日乘子法的来历与直观解释

矩阵特征值是这么来的，以及有趣的盖尔圆

万能的 SVD 分解是哪位牛人提出来的？

矩阵分解术，不得不从高斯说起

1高斯与二次型

2推导

验证

二次型的标准化

程序小实验

与其他分解的关系

3小结

矩阵和线性代数原来是这么来的

概率论原来可以这样优雅地入门

机器学习的数学基础之向量范数

机器学习的数学基础之矩阵范数

拉格朗日乘子法的来历与直观解释

矩阵分解术，不得不从高斯说起

1高斯与二次型

2推导

验证

二次型的标准化

程序小实验

与其他分解的关系

3小结

矩阵和线性代数原来是这么来的

概率论原来可以这样优雅地入门

机器学习的数学基础 之 向量范数

机器学习的数学基础 之 矩阵范数

拉格朗日乘子法的来历与直观解释

机器学习的数学基础之向量范数

机器学习的数学基础之矩阵范数