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

2020-11-01 406

高斯有研究矩阵分解？

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

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

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

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

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

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

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

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 分解是哪位牛人提出来的？

本文分享自微信公众号 - 机器学习与数学（Mathinside2016）。
如有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一起分享。

微信关注我们

原文链接：https://my.oschina.net/mathinside/blog/4699711

转载内容版权归作者及来源网站所有！

低调大师中文资讯倾力打造互联网数据资讯、行业资源、电子商务、移动互联网、网络营销平台。持续更新报道IT业界、互联网、市场资讯、驱动更新,是最及时权威的产业资讯及硬件资讯报道平台。

持续坚持原创输出，点击蓝字关注我吧作者：小傅哥博客：https://bugstack.cn 沉淀、分享、成长，让自己和他人都能有所收获！😜 目录一、前言二、2万元私活收入 1. 企业门户网站(5000元) 2. 卖家具宣传网站(2000元) 3. Netty通信框架(2000元) 4. 毕设、讲课、数据采集(11000元) 三、心得体会四、系列推荐一、前言 5:30下班让我有更多的时间！ 🤔13年~15年，我还从来没想过一天上班，还能干到6点以后！因为我上班的第一家公司是偏传统的外包公司，与互联网公司不同。不知道现在如何了，但当时： 9点到公司，6点下班，但基本大家5:30也就走了。每天下班都能看见日落☀️ 一个C#写的项目代码能运行十年。一个Java程序员写了两年C# 一年会出那么几次差，全国各个地方都几乎会去。一年去了四次上海！也正是因为这样，5:30 下班回家的我，吃完饭也就7点。剩下了大把的时间让我不停的折腾。也就是那段时间接了一些私活，共计赚了2万块。那段期间租住在距离卢沟桥很近的大瓦窑，村落里的平房小院子住起来还是蛮舒服的，开开心心的住了两年。 ...

2020-11-02

546

【专题】多角度深入解析开放原子开源基金会

2020 年 9 月 9 日，开放原子开源基金会正式对外发声，同天，宣布百度超级链正式成为基金会首个捐赠项目。 9 月 10日，华为 OpenHarmony 操作系统开源，开放原子开源基金会获捐 OpenHarmony。 9 月 28 日，开放原子开源基金会组建基金会孵化项目导师团队——“伯乐招募计划”，以为更好地推动开放原子开源基金会项目的孵化毕业。 10 月 15 日，开放原子开源基金会和深圳市开放原子开源技术服务中心主办的“有你有我开源节”开源技术研讨会暨深圳·国际开源谷揭牌仪式在深圳启动。正式全面向外界介绍基金会组织和项目，传播开源文化。至此，开放原子开源基金会的组织架构、性质、能力等已经基本呈现，早期的项目捐赠、孵化工作正在推进，这些信息共同勾勒出一个较为完整的基金会形象。本专题围绕开源节论坛所讨论的话题展开，以期较为全面地介绍我国首个开源基金会。重新认识基金会基金会是为公益目的或者其他非营利目的成立，不向出资人、设立人或者会员分配所取得利润的法人，为非营利法人。基于此，开放原子开源基金会宗旨目标是推动科学领域的相关公益事业的发展，不分配所取得利润，其所取得利润属于...

2020-11-02

574

发表评论

资源下载

更多资源

优质分享App

近一个月的开发和优化，本站点的第一个app全新上线。该app采用极致压缩，本体才4.36MB。系统里面做了大量数据访问、缓存优化。方便用户在手机上查看文章。后续会推出HarmonyOS的适配版本。

Mario

马里奥是站在游戏界顶峰的超人气多面角色。马里奥靠吃蘑菇成长，特征是大鼻子、头戴帽子、身穿背带裤，还留着胡子。与他的双胞胎兄弟路易基一起，长年担任任天堂的招牌角色。

腾讯云软件源

为解决软件依赖安装时官方源访问速度慢的问题，腾讯云为一些软件搭建了缓存服务。您可以通过使用腾讯云软件源站来提升依赖包的安装速度。为了方便用户自由搭建服务架构，目前腾讯云软件源站支持公网访问和内网访问。

Spring

Spring框架（Spring Framework）是由Rod Johnson于2002年提出的开源Java企业级应用框架，旨在通过使用JavaBean替代传统EJB实现方式降低企业级编程开发的复杂性。该框架基于简单性、可测试性和松耦合性设计理念，提供核心容器、应用上下文、数据访问集成等模块，支持整合Hibernate、Struts等第三方框架，其适用范围不仅限于服务器端开发，绝大多数Java应用均可从中受益。

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

1高斯与二次型

2推导

验证

二次型的标准化

程序小实验

与其他分解的关系

3小结

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

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

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

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

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

刚毕业不久，接私活赚了2万块！

【专题】多角度深入解析开放原子开源基金会

相关文章

发表评论

资源下载

优质分享App

Mario

腾讯云软件源

Spring

欢迎您来访！

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

1高斯与二次型

2推导

验证

二次型的标准化

程序小实验

与其他分解的关系

3小结

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

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

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

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

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

刚毕业不久，接私活赚了2万块！

【专题】多角度深入解析开放原子开源基金会

相关文章

发表评论

资源下载

优质分享App

Mario

腾讯云软件源

Spring

欢迎您来访！

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

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