线性方程组

高斯消元法

高斯消元法(Gaussian Elimination Method)即加减消元法,是求解线性方程组的经典方法

通过将一个方程减掉另一个方程的倍数消掉末知数,得到阶梯型方程组,然后依次解出每一个末知数

下面用一个简单的例子进行说明,对于如下的线性方程组

先消去方程2和方程3的第一个末知数,将方程2减去方程1的3倍,将方程3加上方程1,消掉方程2和方程3中的$x_{1}$,得

然后将方程3加上方程2的3倍,消掉方程3中的$x_{2}$,得

可以解得

下面用矩阵的形式描述这一求解过程,如下所示

下面将这种消元法进行推广,对于任意的线性方程组,采用如下的初等变换对其进行变形,方程组的解不变

  1. 交换两个方程的位置
  2. 用非0的常数乘以某方程的两端
  3. 将一个方程的常数倍加到另一个方程上去

采用这种初等变换,每次消掉一个末知数,最后得到一个阶梯形方程组,即可求出方程的解

齐次方程组

齐次线性方程组(Homogeneous Linear Equations)是常数项全部为0的线性方程组,可以写成如下形式

其中$\boldsymbol{A} \in \mathbb{R}^{m \times n}, \boldsymbol{x} \in \mathbb{R}^{n}$

将系数矩阵$\boldsymbol{A}$按列分块为,齐次方程可以写成

以解向量为组合系数,向量组的线性组合为0向量

显然$\boldsymbol{x}=\mathbf{0}$是方程组的解,因此齐次方程一定有解,更重要的是,除$\boldsymbol{x}=\mathbf{0}$ 外的解,称为非0解,下面讨论这种解的存在性

根据线性相关性的定义,如果向量组线性无关,则不存在一组不全为0的系数使得其线性组合为0

如果向量组线性相关,则存在一组不全为0的系数使得其线性组合为0,这就是方程组的非0解

前者对应于矩阵$A$的秩为$n$,后者秩小于$n$,由此得到齐次方程组解的存在性判定条件,分下面两种情况

  1. 如果$r(\boldsymbol{A})=n$,方程组只有0解
  2. 如果$r(\boldsymbol{A}) \lt n$,方程组有非0解

方程组有非$\boldsymbol{0}$解的充分必要条件是$r(\boldsymbol{A})<n$

如果$\boldsymbol{A}$是方阵,$r(\boldsymbol{A})<n$等同于$\boldsymbol{A}$不可逆, 如果$m<n$,即(方程的数量)小于未知数的数量,则有

此时方程组必定有非$\boldsymbol{0}$解,对于如下的线性方程组

其系数矩阵的秩为

因此方程组只有$\mathbf{0}$解,对于如下的方程组

其系数矩阵的秩为

也是方程组的解,证明如下

假设都是方程组的解,如果这组解线性无关且方程组的任意一个解都可以由这组解线性表示,则称是方程组的一个基础解系

如果$r(\boldsymbol{A})<n$,则存在基础解系,且基础解系中包含$n-r(\boldsymbol{A})$个解

齐次线性方程组的求解方法

通常采用的是初等行变换法,对应高斯消元法,经过初等行变换将系数矩阵化为阶梯形矩阵之后

如果出现自由未知数,可以将它们设为特殊值,形成基础解系,然后得到方程组的通解(General Solution)

如果是自由末知数,通常将它们的值依次设为

这是$\mathbb{R}^{n-r}$空间一组最简单的标准正交基,然后根据它们的值解出其他的末知数,对于如下的方程组

对其系数矩阵进行初等行变换

由于$r(\boldsymbol{A})=2<4$,因此方程组有非$\mathbf{0}$解,最后丙个末知数为自由变量

,得到基础解系的第一个解

,得到基础解系的第二个解

方程组的通解为

其中为任意常数

非齐次方程组

非齐次线性方程组(Non-homogeneous Linear Equations)的常数项不全为0,可写成如下形式

这与一元一次方程$a x=0$在形式上是统一的,方程组的增广矩阵是系数矩阵和常数向量合并构成的矩阵

对于如下的线性方程组

其系数矩阵为

增广矩阵为

假设$\boldsymbol{A} \in \mathbb{R}^{m \times n}, \boldsymbol{x} \in \mathbb{R}^{n}$ ,系数矩阵$\boldsymbol{A}$扩列分块为,非齐次方程可以写成

以$x$为组合系数,向量组的线性组合为向量,如果可以由的列向量线性表示,则方程组有解,否则方程组无解

用初等行变换将增广矩阵化为阶梯矩阵

如果,则意味着出现矛盾方程,方程无解,如果,则方程组有解,对于第二种情况,增广矩阵的秩与系数矩阵的秩小于增广矩阵的秩,且

由此得到非齐次方程组解的存在性判定条件

  1. 如果$r(\boldsymbol{A})=r(\boldsymbol{B})$,那么方程组的解存在
  2. 如果$r(\boldsymbol{A})<r(\boldsymbol{B})$,那么方程组的解不存在

对于第一种情况,如果$r(\boldsymbol{A})=n$,那么方程组有唯一解,如果$r(\boldsymbol{A})<n$,那么方程组有无穷多组解

解的性质与结构

如果是非齐次方程组所对应的齐次方程的一组解,是非齐次方程的一个解,则是非齐次方程的解,显然

如果是齐次方程组的基础解系,是非齐次方程组的一个解,则非齐次方程组的解可以表示为

同样可以用初等行变换求解非齐次方程组,其解为对应的齐次方程组的通解加上它的一个特解(Particular Solution)

齐次方程组通解的求解方法在前面已经介绍,非齐次方程组的特解可以任意选取,通常令自由末知数的值全为0

用初等行变换解下面的非齐次线性方程组

对其增广矩阵进行初等行变换

这里是自由未知数,令,得到一个特解

方次方程组的基础解系为

因此方程的解为

其中为任意常数

Python中linalg的solve函数提供了求解非齐次线性方程组的功能,函数的传人参数为系数矩阵$A$,以及常数向量$\boldsymbol{b}$,返回值是方程组$\boldsymbol{A x}=\boldsymbol{b}$的解向量$x$,对于方程组

下面是求解该方程组的代码

1
2
3
4
5
import numpy as np
A = np.array([3,1], [1,2])
b = np.array([9,8])
x = np.linalg.solve(A, b)
print(x)

程序运行结果为$[2,3]$

特征值和特征向量

特征值(Eigenvalue,也称为本征值)与特征向量(Eigenvector,也称为本征向量)决定了矩阵的很多性质

从几何的角度来看,持征向量是经过矩阵的线性变换仍然处于同一条直线上的向量

特征值与特征向量

特征值与特征向量

对于$n$阶矩阵$A$,其特征向量是经过这个矩阵的线性变换之后仍然处于同一条直线上的向量,新向量的方向可能会相反,长度可能会改变,即存在一个数$\lambda$以及非$0$向量$x$,满足

则称$\lambda$为矩阵$A$的特征值,$x$为该特征值对应的特征向量

特征值是特征向量在矩阵的线性变换下的缩放代,如果持征值大于0,那么经过线性变换之后特征向量的方向不变;如果特征值小于0,那么经过线性变换之后特征向量的方向相反,如果特征值为0,则经过线性变换之后特征向量收缩回原点

特征矩阵和特征方程

上式变形后可以得到

$A-\lambda I$ 称为特征矩阵,按照线性方程组的理论,上面的齐次方程有非0解的条件是系数矩阵的行列式必须为0,即

式子称为特征方程(Eigenvalue Equation)

对于矩阵

其特征方程为

特征多项式

上面的行列式展开之后是$\lambda$的$n$次多项式,称为矩阵的特征多项式(Characteristic Polynomial),为如下形式

求解这个特征多项式对应的特征方程可以得到所有特征值,方程的根可能是复数,此时的特征值为复数,特征向量为复向量

根据对角行列式的计算公式,对角矩阵的特征为其主对角线元素,对于如下对角矩阵

其特征方程为

类似地,上三角矩阵的特征值为其主对角线元素,对于如下的上三角矩阵

其特征方程为

对于下三角矩阵有相同的结论,一种计算特征值的方法是通过相似变换将矩阵变为上三角矩阵

代数重数、谱

根据多项式分解定理,特征方程可以写成

其中称为特征值代数重数(Algebraic Multiplicity),根据代数方程的理论,有

所有$N_{\lambda}$个不同的特征值构成的集合称为矩阵的(Spectrum)

矩阵的谱半径(Spectral Radius)定义为所有特征值绝对值的最大值,记为

如果矩阵$\boldsymbol{A}$不可逆,则

因此0是它的特征值,反之,如果可逆,则0不是它的特征值

几何重数

得到每个特征值$\lambda_{i}$之后,解下面的线性方程组

即可得到其对应的特征向量,此方程组有$1 \leq m_i \leq n_i$,个线性无关的解,称$m_i$为$\lambda _i$的几何重数(Geometric Multiplicity)

这些线性无关的解构成的空间称为矩阵$\boldsymbol{A}$关于特征值的特征子空间,记为

根据齐次线性方程组解的理论,特征子空间的维数为

属于不同特征值的特征向量线性无关,矩阵所有线性无关的特征向量的数量为

显然有$N_{x} \leqslant n$

特征值与特征向量的计算

下面用一个例子来说明特征值与特征向量的计算,对于如下矩阵

其特征多项式为

特征方程$-(1-\lambda)(1+\lambda)=0$的根为$\lambda=1$与$\lambda=-1$,因此该矩阵的特征值为1与-1

将特征值1代人,可得

该齐次方程的解为

此即特征值1对应的特征向量

将另外一个特征值-1代人,可得

该齐次方程的解为

即特征值-1对应的特征向量

上三角矩阵的特征值为其主对角线元素,对于如下矩阵

根据前面的结论,其特征多项式为

其特征值为 $1 、 2 、 3$

对于不超过4阶的矩阵,可通过解特征方程得到特征值

但更高次方程的求根存在困难,阿贝尔-鲁菲尼(Abel-Ruffini)定理指出,4次以上的代数方程没有公式解,对于一般的高次方程,方程系数的有限次四则运算、开方运算的结果均不可能是方程的根,因此高阶矩阵的特征值只能求近似解

直接求解特征方程并不是一种好的选择,更有效的方法是迭代法

Python中linalg的eig函数实现了计算矩阵的特征值与特征向量的功能,函数的输入为方阵,输出为所有的特征值以及这些特征值对应的特征向量

1
2
3
4
5
import numpy as np
A = np.array([[1,0,0],[0,1,0],[0,0,5]])
eigvalues, eigvectors = np.1inalg.eig(A)
print(eigvalues)
print(eigvectors)

程序运行结果为

1
2
3
4
5
6
#  特征值
[1. 1. 5.]
# 特征向量
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]

矩阵的迹

下面介绍特征值与矩阵主对角线元素以及行列式的关系,矩阵的(Trace)定义为其主对角线元素之和

对于如下矩阵

其迹为

关于矩阵的迹,有下面的公式成立

根据韦达定理,下面的$n$次方程

所有根之和为

所有根的乘积为

下面计算$n$阶矩阵的特征多项式,首先将行列式写成下面的形式

然后按照第1列拆开、变为两个行列式之和

接下来将这两个行列式均按照第2列拆开,变为4个行列式之和

依此类推,将上一步的结果中所有行列式按照下一列拆开,最后可以得到$2^{n}$个行列式,特征值多项式是它们之和

这些行列式的展开结果中,含有$\lambda^{n}$的只有

因此特征多项式的首次项就是$(-1)^{n} \lambda^{n}$ ,含有$\lambda^{n-1}$的是下面$n$个行列式

它们之和为

因此特征多项式的$\lambda^{n-1}$项系数是,不含的只有下面一个行列式

因此特征多项式中常数项的系数为$|\boldsymbol{A}|$,由此可以得到特征多项式为

将特征多项式乘以$(-1)^{n}$可以变为

因此矩阵所有特征值的和为矩阵的迹

所有特征值的积为矩阵的行列式

特征值的若干重要性质

1️⃣如果矩阵可逆且为它的特征值,则的特征值,根据特征值与特征向量的定义有

上式两边同时左乘$A^{-1}$,可以得到

因此的特征值,为对应的特征向量

2️⃣如果$\lambda$是矩阵$A$的特征值,则$\lambda^{n}$是$\boldsymbol{A}^{n}$的特征值,根据特征值与特征向量的定义有

反复利用此式,有

因此$\lambda^{n}$是$\boldsymbol{A}^{n}$的特征值,类似地可以证明如果$\lambda$是矩阵$\boldsymbol{A}$的特征值,则$k \lambda$是$k \boldsymbol{A}$的特征值,对于如下的多项式

3️⃣如果$\lambda$是矩阵$\boldsymbol{A}$的特征值,则$f(\lambda)$是$f(\boldsymbol{A})$的特征值

4️⃣矩阵$A$与$A^{\mathrm{T}}$有相同的特征值,显然

因此

特征向量的若干重要性

1️⃣如果向量都是矩阵关于同一个特征值的特征向量,则它们的非$\mathbf{0}$线性组合

也是矩阵$\boldsymbol{A}$关于$\lambda$的特征向量,根据特征值与特征向量的定义有

因此是关于的特征向量

2️⃣矩阵属于不同特征值的特征向量线性无关

假设矩阵$A$的$l$个不同特征值为,它们对应的特征向量为,下面用归纳法进行证明

当$l=1$时结论成立,因为,如果,则必定有

假设当$l=m$时结论成立,当$l=m+1$时,有

式子两边左乘$\boldsymbol{A}$,有

由于

因此

将第一个式子乘以$\lambda_{m+1}$,然后减去上式,可得

由于线性无关,因此

,因此,将代入第一个式子可得

由于是特征向量,因此,故,因此线性无关

3️⃣实对称矩阵的特征值均为实数

首先定义矩阵的共轭运算,复数矩阵$\boldsymbol{A}$的共轭$\overline{\boldsymbol{A}}$为将其所有元素共轭后形成的矩阵,例如对于下面的矩阵

其共轭矩阵为

可以证明共轭运算满足下面的性质

显然对于实矩阵有

假设$\lambda$是实对称矩阵$\boldsymbol{A}$的特征值,$x$是对应的特征向量

由于是实对称矩阵,因此$\bar{A}^{\mathrm{T}}=A$,由于$A x=\lambda x$,因此

上式两边同时右乘$\boldsymbol{x}$可以得到

从而有

由于$x \neq 0$,因此$\overline{x^T} x>0$,可以得到 $\lambda=\bar{\lambda}$,这意呠着$\lambda$是实数

4️⃣实对称矩阵属于不同特征值的特征向量相互正交

假设$A$为实对称矩阵,是它的两个不同的特征值,分别为属于的特征向量,则有

上式的第一式两边左乘$x_{2}^{\mathrm{T}}$可以得到

因此有

由于,因此

机器学习中使用的矩阵一般为实对称矩阵,因此特征值均为实数,且不同特征值的特征向量正交

特征值和特征向量被大量用于机器学习算法,典型的包括主成分分析(PCA),线性判别分析(LDA),流形学习等降维算法

相似变换

通过相似变换可以将一个矩阵变为对角矩阵

定义

如果有两个矩阵$A$、$B$以及一个可逆矩阵$P$满足

则称矩阵$A$,$B$相似,记为$A \sim B$,上式称为相似变换,$P$为相似变换矩阵

性质

相似具有自反性,矩阵与其自身相似,即$A \sim A$,显然

相似具有对称性,如果$A \sim B$,则$B \sim A$,由于

上式两边左乘$P$,右乘$P^{-1}$,可以得到

相似具有传递性,如果$A \sim B$且$B \sim C$,则$A \sim C$,由于$A \sim B$且$B \sim C$,因此有

从而有

相似矩阵有相同的特征值,这意昩着相似变换保持矩阵的特征值不变

假设$A \sim B$,则存在可逆矩阵$\boldsymbol{P}$使得

因此

这一性质可用于求解特征值,通过相似变换将矩阵$A$变为对角矩阵或三角矩阵

特征值不变,对角矩阵或三角矩阵的主对角线元素即为$A$的特征值

如果矩阵满足一定的条件,通过相似变换可将其转为对角矩阵

假设是$n$阶矩阵$A$的$n$个特佂值,是它们对应的特征向量根据特征值与特征向量的定义有

如果令矩阵,对角矩阵

根据右乘对角矩阵的性质有

如果矩阵$P$可逆,那么上式两边同时左乘$P^{-1}$可以得到

通过这种相似变换可以将矩阵化为对角矩阵,称为矩阵的相似对角化

式子意味着可以以矩阵$A$的特征向量为列构造一个矩阵$P$,通过它将矩阵对角化,得到以$A$的特征值为主对角线的对角矩阵$\boldsymbol{\Lambda}$

这种做法成立的条件是矩阵$\boldsymbol{P}$可逆,即矩阵$A$有$n$个线性无关的特征向量

正交变换

对于实对称矩阵,我们可以构造一个正交的相似变换将其对角化

用归纳法证明实对称矩阵一定可以对角化

这意味着$n$阶实对称矩阵有$n$个线性无关的特征向量,实对称矩阵$A$属于不同特征值的特征向量是相互正交的

如果用格拉姆-施密特正交化将同一个特征值的所有特征向量正交化,然后将所有特征向量单位化,可以得到一组标准正交基

以它们为列构造相似变换矩阵$\boldsymbol{P}$,则矩阵$\boldsymbol{P}$是正交矩阵,可通过正交变换(Orthogonal Transformation)将矩阵化为对角阵

由于$P^{T}=P^{-1}$,因此这是一种更特殊的相似变换,实现时只需要对同一个特征值的不同特征向量正交化,然后将所有正交化之后的特征向量进行标准化即可

例子

将矩阵通过正交变换化为对角矩阵。对于下面的矩阵

其特征多项式为

因此其特征值为$2,-1 ,-1$,当$\lambda=2$时,有

解得

当$\lambda=-1$时,有

解得

正交单位化之后为

则有

正交变换具有一个优良的性质,它可以保持矩阵的对称性

假设$A$是对称矩阵,$P$是正交矩阵,使用下面的正交变换

$B$仍然是对称矩阵,下面给出证明,显然有

豪斯霍尔德(Householder)变换

特殊的正交变换豪斯霍尔德(Householder)变换,它在QR算法以及其他矩阵算法中有重要的应用

首先定义Householder矩阵,为如下形式

其中$\boldsymbol{w}$是$n$维非0列向量,且有$|\boldsymbol{w}|=1$,显然矩阵$\boldsymbol{P}$是对称矩阵,并且是正交矩阵

由于$P$是对称矩阵,因此有

故该矩阵是正交矩阵,通常将$P$写成如下形式

其中$u$为任意非$\mathbf{0}$向量且

这里用$H$对$u$进行了标准化

对于$n$维列向量$\boldsymbol{x}$,构造下面的向量

其中单位向量$\boldsymbol{e}_{1}=\left(\begin{array}{llll}1 & 0 & \cdots & 0\end{array}\right)^{\mathrm{T}}$

根据$\boldsymbol{u}$用上面的式子构造Householder矩阵$\boldsymbol{P}$,下面将向量$x$左乘$P$的结果

其中$x_{1}$是$x$的第1个分量

这表明将列向量$\boldsymbol{x}$左乘$\boldsymbol{P}$之后将零化$\boldsymbol{x}$除第1个元素之外的所有元素,同时保持向量的长度不变,将行向量右乘该矩阵之后有类似的效果

根据这一特性,我们可以构造以Householder矩阵为基础的正交变换,将矩阵转化为类似对角矩阵的形式,零化主对角线之外的元素

对于对称矩阵$\boldsymbol{A}$,使用它的第1列计算向量$u$,按照上面的式子构造Householder矩阵$P$

然后对$A$进行正交变换,这里的正交变换通过将矩阵$A$先左乘$P$,然后右乘$P$实现

左乘$P$实现$A$的第1列的零化,右乘$P$实现$A$的第1行的零化

下面来看矩阵$P$的构造

如果用$A$的整个第1列作为向量,按照上面的式子构造$P$,虽然可在左乘$P$之后将$A$的第1列除第1个元素之外的所有元素全部零化,但会改变$A$的第1行所有元素的值

接下来在右乘$P$的时候无法保证将$P A$的第1行零化,因此$P$需要保证将$A$的第1列的元素零化的同时确保$A$的第1行的元素不变,以便在右乘$P$的时候将这个行零化

我们可以按照下面的形式构造$P$

其中

是用$\boldsymbol{A}$的第1列的后面$n-1$个元素按照上面式子构造的

我们将上上式的矩阵作为第1次豪斯霍尔德变换的矩阵,记为,将左乘之后可以保证的第1行元素不变,同时将的第1列的后面个元素全部变为0

接下来右乘,由于是对称矩阵,因此第1列和第1行相同,右乘可以将第1行后面个元素全部变为0,并且不改变第1列所有元素的值,因此不会破坏前面的列零化结果

然后进行第2次豪斯霍尔德变换

由于正交变换可以保持矩阵的对称性,因此仍然是对称矩阵,用的第2列的后面个元素构造

其中

根据$A_{1}$的第2列的后面$n-2$个元素按照上面的式子构造

经过第2次豪斯霍尔德变换可以将$A_{1}$的第2列的后面$n-3$个元素,第2行的后面$n-3$个元素全部变为0

依此类推,经过$n-2$次豪斯霍尔德变换,可以将对称矩阵化为如下的对称三对角矩阵(Tridiagonal Matrix)

这种矩阵除主对角线、主对角线以上及以下的对角线之外,其他元素均为0

对于一般的$n$阶矩阵$\boldsymbol{A}$,用同样的方法构造豪斯霍尔德矩阵

左乘$P$之后将$A$的第1列后面$n-2$个元素零化,同时保持$\boldsymbol{A}$的第1行元素不变

由于$\boldsymbol{A}$不是对角矩阵,其行和列不相等,因此右乘$\boldsymbol{P}$ 时候无法将其第1行元素零化

同样不能用完整的列构造豪斯霍尔德变换矩 阵,否则右乘该矩阵的时候会破坏前面零化的结果

用和对称矩阵相同的方法构造变换矩阵,第一次豪斯霍尔德变换之后的结果为

第二次豪斯霍尔德变换可以将第2列的后面$n-3$个元素零化,变换之后的结果为

依次类推,通过$n-2$次豪斯霍尔德变换可以将$A$化为如下形式的上海森堡矩阵(upper-Hessenberg form)

这种矩阵除主对角线及以上,主对角线下面的对角线的元素外,其他的元素均为0

QR 算法

下面介绍求解高阶矩阵特征值的QR筫法,它被誉为20世纪十大算法之一

它依赖于$\mathrm{QR}$分解,对于一个矩阵$A$,$\mathrm{QR}$分解将其化为一个正交矩阵$Q$与一个上三角矩阵$R$的乘积

$\mathrm{QR}$算法是一种迭代法,从矩阵开始,每次构造一个相似变换,将变换为,最后收敛到一个上角矩阵,主对角线元素即为其特征值

由于矩阵相似,因此它们有相同的特征值,得到了的特征值即得到了的特征值

问题的核心是如何构造这种相似变换,这借助于$\mathrm{QR}$分解实现,每次迭代时,首先对$\boldsymbol{A}_{k}$进行$\mathrm{QR}$分解

然后用分解结果构造一个新的矩阵$\boldsymbol{A}_{k+1}$,这里将$\mathrm{QR}$分解的结果矩阵交换顺序后相乘

上面两个式子给出了根据当前矩阵构造下一个矩阵的方式,称为$\mathrm{QR}$迭代

因为是相似的,式两边同时左乘可以得到

将其代人式$A_{k+1}$可得

由于相似具有传递性,因此$A$与$A_{k}, k=1, \cdots, n$相似

如果$\boldsymbol{A}$满足一定的条件,那么$\mathrm{QR}$迭代所产生的矩阵序列${\boldsymbol{A}_{k}}$将收敛到一个上三角矩阵,主对角线元素即为$\boldsymbol{A}$的特征值

$Q R$迭代是正交变换,如果$A$是对称矩阵,这种变换将保持对称性,且收敛到上三角矩阵,因此最终会收敛到对角矩阵

对于实对称矩阵$A$,$Q R$迭代代产生的所有正交矩阵$Q_{k}$的乘积的所有列,即为$\boldsymbol{A}$的特征向量,由于

因此

如果令

因此$P$的列为$A$的特征向量
下面来看$\mathrm{QR}$算法的一个例子,对于如下矩阵

$Q R$算法第1次迭代的结果为

再经过5次迭代后的结果为

此时矩阵$A_{6}$已经接近于上三角矩阵,主对角线元素接近于$A$的特征值,事实上,矩阵$A$的特征值为$1,2,3$

为了加速收敛,通常采用带位移的$QR$算法(Shifted QR Algorithm),在第$k$次迭代时,对于设定的移位常数$s_{k}$,迭代公式为

即先将矩阵减掉后再进行分解,在构造时再将加回来

按照这种迭代公式,也是相似的,根据上式的1式有

将其代人上式的2式,可得

因此相似

移位系数$s_{k}$的计算

一种方案是选择右下角子矩阵的两个特征值中接近于元素的那个特征值,以便于使得经过此次迭代后变为0

计算此子矩阵的特征值可以通过求解特征方程实现

这里即为的一个特征值,对剩下的矩阵继续进行算法迭代, 得到所有的特征值

对于一般的矩阵,QR算法的收敛速度可能很慢,如果将矩阵变换成接近于三角阵,则能加快收敛速度

如果是对称矩阵,可先用豪斯霍尔德变换将其化为对称三对角矩阵,然后用$Q R$算法进行达代,收敛到对角矩阵

求解特征值的整个流程为: 对称矩阵$\rightarrow$对称三对角矩阵$\rightarrow$对角矩阵

对于普通矩阵,可用豪斯霍尔德变换将其化为上海森堡矩阵,然后用QR算法进行迭代,收敛到上三角矩阵

整个流程为: 普通矩阵$\rightarrow$上海森堡矩阵$\rightarrow$三角矩阵

广义特征值

广义特征值(Generalized Eigenvalue)是特征值的推广,定义于两个矩阵之上

对于方阵$A$和$B$,如果存在一个数$\lambda$及非0向量$\boldsymbol{x}$,满足

则称$\lambda$为广义特征值,$x$为广义特征向量,类似的有特征方程

如果矩阵$B$可逆,对第一个式子左乘$B^{-1}$,则式1的问题等价于下面的特征值问题

广义特征值在机器学习中被广泛使用,包括流形学习、谱聚类算法,以及线性判别分析等

瑞利商

瑞利商

方阵$A$和非0向量$x$的瑞利商(Rayleigh Quotient)定义为如下比值

根据式子的定义,对于任意的非0实数$k$,有

即对向量缩放之后其瑞利商不变,瑞利商存在冗余,证明如下

假设是矩阵的最小特征值,是最大特征值,则有

即瑞利商的最小值为矩阵的最小特征值,最大值为矩阵的最大特征值

并且当$x$分别为最小和最大的特征值对应的特征向量的时候取得这两个值,可以用拉格朗日乘数法证明

由于将向量乘以非0系数之后瑞利商不变,因此式1极值问题的解不唯一

为此增加一个约束条件以保证解的唯一性,同时简化问题的表述,限定$x$为单位向量,有下面的等式约束

增加此约束之后,瑞利商变为

构造拉格朗日乘子函数

对$\boldsymbol{x}$求样度并令栖度为$\boldsymbol{0}$可以得到

这里利用了矩阵与向量求导公式,上式等价于

此结果意味着瑞利商的所有极值在矩阵的特征值与特征向量处取得

假设的第个待征值、是其对应的特征向量,将它们代人瑞利商的定义,可以得到

因此,在最大的特征值处,瑞利商有最大值;在最小的特征值处,瑞利商有最小值

瑞利商在机器学习领域的典型应用是主成分分析

广义瑞利商

对瑞利商进行推广,得到广义瑞利商,定义为

同样,广义瑞利商存在余,将向量$x$缩放之后广义瑞利商不变

假设对任意的非0向量$x$,有$\boldsymbol{x}^{\mathrm{T}} B x>0$

如果令$B=C C^{\mathrm{T}}$,这是对矩阵$B$的楚列斯基(Cholesky)分解,同时令$x=\left(C^{\mathrm{T}}\right)^{-1} y$,则可以将广义瑞利商转化成瑞利商的形式

根据瑞利商的结论,广义瑞利商的最大值和最小值由矩阵$C^{-1} A\left(C^{\mathrm{T}}\right)^{-1}$的最大和最小特佂值决定

也可以直接通过广义特征值得到广义瑞利商的极值,与瑞利商类似,加上等式约束消掉最优解的冗余

广义瑞利商变为

可以用拉格朗日乘数法求解,构造拉格朗日乘子函数

对$\boldsymbol{x}$求梯度并令梯度为$\boldsymbol{0}$,可以得到

这等价于

这是广义特征值问题,如果矩阵$\boldsymbol{B}$可逆,那么上式两边同乘以其逆矩阵可以得到

因此广义瑞利商的所有极值在上面的广义特征值处取得,假设是第个广义特征值,是其对应的特征向量,将它们代人广义瑞利商的定义,可以得到

因此广义瑞利商的极大值在最大广义特征值处取得,极小值在最小广义特征值处取得

线性判別分析的优化目标函数即为广义瑞利商

谱范数与特征值的关系

前面定义了谱范数的概念,可以证明矩阵$\boldsymbol{W}$的谱范数等于$\boldsymbol{W}^{\mathrm{T}} \boldsymbol{W}$的最大特征值的平方根,即$\boldsymbol{W}$最大的奇异值

其中的奇异值,是的特征值的平方根,奇异值将在后面详细介绍,根据定义,谱范数的平方为

它就是瑞利商的极大值,上一节已经证明了这一最优化问题的解是矩阵$\boldsymbol{W}^{\mathrm{T}} \boldsymbol{W}$的最大特征值,因此结论成立

Python中linalg的norm函数提供了计算矩阵范数的功能,函数的输人参数为要计算的矩阵,以及范数的类型,如果类型值为2,则计算谱范数

1
2
3
4
import numpy as np

A = np.array([[0,1,2,3],[4,5,6,7],[8,9,10,11]])
print(n)

程序运行结果为22

条件数

如果矩阵$W$可逆,条件数(Condition Number)定义为它的范数与它的逆矩阵范数的乘积

这里的范数可以是任何一种范数

如果使用谱范数,则$|\boldsymbol{W}|$等于$\boldsymbol{W}$的最大奇异值,$\left|\boldsymbol{W}^{-1}\right|$等于$W$最小奇异值的逆,根据谱范数的定义有

上式第二步进行了换元,令$\boldsymbol{W}^{-1} y=x$,根据前二节的结论,$\frac{|\boldsymbol{W} \boldsymbol{x}|}{|\boldsymbol{x}|}$的最小奇异值

此时条件数等于矩阵的最大奇异值与最小奇异值的比值

其中的奇异值,显然矩阵的条件数总是大于或等于1

条件数决定了矩阵的稳定性,一个矩阵的条件数越大,则它越接近于不可逆矩阵,矩阵越”病态”

条件数在诸多算法的稳定性分析中有重要的作用,Python中linalg的cond函数实现了计算矩阵条件数的功能

1
2
3
4
5
import numpy as np

A = np.array([[1,0,0],[0,1,0],[0,0,5]])
c = np.linalg.cond(A)
print(c)

程序运行结果为$5.0$,是该矩阵最大特征值与最小特征值的比值,在这里奇异值等于特征值

应用——谱归一化与谱正则化

正则化

正则化是机器学习中减轻过拟合的一种技术,它迫使模型的参数值很小,使模型变得更简单,一般情况下,简单的模型有更好的泛化性能

正则化可以通过在目标函数中增加正则化项实现,正则化项通常为参数向量的L1范数,或L2范数的平方

谱正则化(Spectral Regularization)用谱范数构造正则化项

另外一种技术是谱归一化(Spectral Normalization),通过用谱范数对线性映射的矩阵进行谱归一化而确保映射有较小的李普希茨常数,从而保证机器学习模型对输人数据的扰动不敏感

神经网络中用权重矩阵$\boldsymbol{W}$与偏置向量$\boldsymbol{b}$对输人数据$\boldsymbol{x}$进行映射,得到输出结果$W x+b$,如果此映射满足李普希茨条件(这里将其推广到多元函数,将绝对值改为向量的范数),则有

其中$K$为李普希茨常数,将式子变形后可以得到

上式左侧部分的极大值就是权重矩阵的谱范数

如果权重矩阵的谱范数存在一个较小的上界,则神经网络该层的映射有较小的李普希茨常数,从而保证输人值的较小改变不会导致输出值的突变,映射更为平滑

假设权重矩阵的谱范数为$\sigma(\boldsymbol{W})$,如果用它对矩阵进行归一化

则能保证归一化之后的权重矩阵满足$\sigma(\boldsymbol{W})=1$,这由矩阵乘以常数之后的特征值的性质保证

前两节已经证明谱范数是矩阵$W$的最大奇异值,计算矩阵奇异值的代价太大,因此在实现时需要对谱范数$\sigma(\boldsymbol{W})$的值近似计算,可以采用幂迭代法

接下来考虑如何用谱范数为神经网络的目标函数构造正则化项(Spectral Norm Regularizer)

给定训练样本集为输人向量,为标签向量,加上谱正则化项后的目标函数为

其中为对单个样本的损失函数,为神经网络的层数,为第层的权重矩阵

上式第2项为谱正则化项,$\lambda>0$为正则化项的权重,谱正则化项由神经网络所有层权重矩阵的谱范数平方之和构成

可以防止权重矩阵出现大的谱范数,从而保证神经网络的映射有较小的李普希茨常数