1 LU分解

1.1 直观理解

LU分解（lower-upper decomposition）的原理是将一个方阵分解为一个下三角矩阵和一个上三角矩阵的乘积。即满足：LU分解的原理是高斯消元法，具体来说就是通过行变换将化为上三角矩阵。

举例，对于矩阵：进行初等行变换，将下三角部分消去，首先将第一行的1.5倍加到第二行，：第一行的1倍加到第三行，：第二行的-4倍加到第三行，：此时已经化为上三角矩阵，记。

上述行变换操作，可以视为一系列初等矩阵依次左乘原始矩阵：其中每个标注的元素称作消元因子，其位置和大小记录了此次消去操作的信息。

此时，高斯消去过程可以写作：将一系列合并：在上面的例子中，可以计算出为满秩矩阵，因此可以求得：对（1）式两边左乘，得：，这个就是下三角矩阵。

1.2 存在性和唯一性

存在性：

并非所有矩阵存在LU分解。对于方阵存在LU分解的前提是可逆，且的前个顺序主子式不为零（在LU分解中，我们限制了对系数矩阵中行的互换）。

唯一性：

对于可逆阵有，其中为下三角矩阵，为上三角矩阵，且，其中，则LU分解唯一。

：

设可逆阵有两个LU分解：，显然可逆。

因此有：注意到为下三角阵，为上三角阵，其逆矩阵仍然为上/下三角阵。

因此左式为下三角阵，右式为上三角阵，因此：因为和对角元，故对角元仍全部为1，所以：即：

1.3 数值稳定性问题

对于主元为零的矩阵，例如：实际上只需要在这之前添加行交换的步骤即可。

对于矩阵中出现很小的数，比如：分解后得到的和会出现数值很大的数，为了避免以上两种情况，利用换行的思想，可以对上文的高斯消元法进行改进，也就是PLU分解。

2 PLU分解

2.1 直观理解

PLU分解就是在LU分解的基础上添加了行交换的操作。

行交换的准则为：处理每一行时，找到该行主元下方绝对值最大的元素，把这个最大元素所处的行和当前行交换。

以矩阵为例：第一行的主元为零，我们找到主元为零所在列中最大元素，将其所在行（第二行）与当前行（第一行）互换：然后进行消元，得到：然后继续处理第二行，找到所在列最大的元素：交换消元：

2.2 形式表达

由于初等矩阵左乘一个矩阵可以看做初等行变换，比如交换单位矩阵前两行：我们将这种表示行交换的矩阵称为置换矩阵，对此行交换可以看做多个置换矩阵依次左乘。

实际上所有的都可以写作的形式，当没有行互换时，。

置换矩阵满足，因此可以转换为：。

相比LU来说，PLU具有很高的数值稳定性，在scipy中，scipy.linalg.lu默认采用的就是PLU分解。

3 Cholesky分解

3.1 定义

假设矩阵为正定矩阵，根据正定的充要条件可知，存在可逆阵，使得成立。这里的不唯一，因此cholesky分解希望寻找一个唯一的矩阵，也就是上三角矩阵。

但是仅仅上三角矩阵还不能满足唯一性，考虑对角阵，其中对角元素为1或-1的情况，此时有，因此满足，有无数种，因此这里限制为上三角矩阵，且对角线元素均大于零。

同理，可以得到cholesky分解的另一种表示形式，其中为下三角矩阵且对角元均大于零。

下面通过一个例子演示cholesky的具体过程。假设矩阵为：由于我们寻找的满足：因此，观察一阶顺序子式，可以得到，得到正数解。然后观察二阶顺序子式，将代入，可以得到，所以，将前面的结果代入，解得。

依次类推，这是一个递归的求解过程，每次可以求得矩阵的一列，最终我们可以解出所有的，此时即可确定。

3.2 代码实现cholesky分解

该方法可以用代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np


def cholesky(X):
    L = np.zeros((X.shape[0], X.shape[0]))
    for i in range(X.shape[0]):
        square_sum = 0
        for j in range(i):
            sum_ = 0
            for k in range(j):
                sum_ += L[k][j] * L[k][i]
            L[j][i] = (X[j][i] - sum_) / L[j][j]
            square_sum += pow(L[j][i], 2)
        L[i][i] = np.sqrt(X[i][i] - square_sum)
    return L

时间复杂度为级别。要求矩阵的正定非常重要，通过这个代码可以看出，如果矩阵非正定，则可能会出现对负数开平方、除数为零的错误。

3.3 cholesky分解的应用

下面介绍一些常用的cholesky分解的应用。

首先是解方程组，若方程组中，为对称正定矩阵，则可以通过cholesky分解求解，步骤如下：

  1. 对进行cholesky分解，得到
  2. 原方程组变为，令，对原式换元得到方程组，求解这个方程组得到
  3. 将解出的代入，解这个方程组得到

虽然看起来“多了一道手续”，但实际上对于上三角/下三角矩阵，方程组求解是非常容易的。

另外一个应用是矩阵求逆，对于一个对称正定矩阵，要求它的逆实际上就是求，步骤如下：

  1. 对进行cholesky分解，得到
  2. 原方程组变为，令，对原式换元得到方程组，求解这个方程组得到
  3. 将解出的代入，解这个方程组得到，它就是

在numpy中，numpy.linalg.cholesky可以对矩阵进行cholesky分解，分解出来的是下三角矩阵，对于scipy，scipy.linalg.cholesky同样可以进行分解，得到的是上三角矩阵。

4 LDL分解

3.1 定义

在cholesky分解的基础上，我们可以进一步限制，也就是在下三角且对角元为正的基础上，限制对角线元素均为1。此时，对于正定矩阵，存在唯一的对角元素均为1的下三角矩阵和对角元素均为正的对角矩阵，使得成立。

LDL与cholesky分解是等价的，也就是说二者可以相互转换。

：是对角阵，有，因此：令，显然它是下三角矩阵，因此有：

根据这个结论，我们可以先求出cholesky分解，然后求解；也可以先求出LDL分解，再解出对应的cholesky分解。

3.2 scipy中使用LDL分解

在scipy中，可以使用scipy.linalg.ldl进行LDL分解：

1
2
3
X = np.array([[2, -1, 3], [0, 2, 0], [0, 0, 1]])
l, d, perm = scipy.linalg.ldl(X)
l, d, perm = scipy.linalg.ldl(X, lower=False)

其中l和d分别是下三角矩阵和对角矩阵，perm是l的行排列索引数组，取决于参数lower，该参数会在因式分解的下三角和上三角分解之间切换。

5 特征值分解

5.1 特征值与特征向量

矩阵的本质是一种线性映射，如果是在同一个空间中的映射就是线性变换，对应的矩阵就是一个方阵。因此对于矩阵和列向量，就是经过的由到的线性映射。

在几何上，线性映射可以分为伸缩变换和旋转变换，为了方便，以二维矩阵为例。假设，为单位正交基向量，则有：可以发现，经过变换后，的坐标分量变为原先的4倍，的坐标分量变为原先的1倍，在图形上来看就是分别对两个基向量在各自方向上进行不同比例的伸缩变换，4和1分别是它们的伸缩因子。

image-20230730160703930

可以发现，在矩阵的作用下，变换前后的向量并没有发生旋转，而仅仅是原方向上的伸缩变化，因此和被称为矩阵的特征向量。不难发现，在坐标轴上任意取两个正交向量组，它们与和共线，显然这些向量同样具有这样的性质，因此可以使用和表示它们。对应地，伸缩因子4和1被称为矩阵的特征值，这些向量被称为特征值对应的特征向量。

下面考虑平面上不和和共线的向量，它经过映射后，得到。

image-20230730161936776

显然映射后与映射前的向量不是共线的向量，其中，横坐标被伸长为原来的4倍，纵坐标是原先的1倍。这个变换发生了旋转，但可以看出，以特征向量线性组合的任意向量经过对角阵的变换，各个坐标会独立地成比例变换，比例因子就是对应的特征值。

显然，上面的矩阵是比较特殊的对角矩阵，对于一般的矩阵，同样存在这样性质的特征值和特征向量。比如，令矩阵为普通的二阶方阵，经过它映射后的单位正交基向量如下所示：发现映射后得到的向量和不仅发生了伸缩，还发生了旋转变换：

image-20230730163131997

这说明，无法使得以和为基只进行伸缩变换，换句话说，在对和变换时，它们已经离开了与它们共线的直线上。因此，我们要寻找的特征向量不是和。

简单分析可知，我们需要寻找的向量需要具有如下的性质：

  1. 不能是零向量，否则线性变换全部变成零向量，没有意义
  2. 矩阵作用于向量时，只会在与这个向量的同一直线上进行伸缩变换
  3. 寻找到两个线性无关的向量，那么就可以在平面内以它们作为一组基，使得用这组基表达的向量在经过的映射后，使用这组基表达，仅仅是各个坐标独立地成比例变换

这个结论可以扩展到维，因此就有了如下定义，如果设比例因子（也就是特征值）为，经过阶方阵线性变换作用在我们要找的向量上，结果仅仅是将原来的伸缩为原来的倍，则就是的一个特征值，对应的和与之共线的非零向量称为特征值对应的特征向量。用公式表达就是教科书上的：那么，找到的这些特征向量能否构成维空间中的一组基？显然，这需要找到方阵的个线性无关的特征向量。下面的问题就是如何求出矩阵的特征值和特征向量。

显然，如果是简单的矩阵，特征值和特征向量可以直接观察得出，但是对于一般的矩阵，就不容易了。这里从特征值的定义入手，可以得到：这是一个以为未知量的齐次线性方程组，由于特征向量不为零向量，方程一定有非零解，所以方程组的秩，由于和为方阵，因此可知行列式，利用这个性质我们可以求出特征值，然后将特征值代入到原方程中，解得对应的特征向量。另外，，使用哪一种都可以。

对于上面的矩阵，，解得两个特征值。分别代入原方程组，得到：解得对应的特征向量为，对应的特征向量为，其中均不为零。发现这两个特征向量线性无关。

可以证明，不同特征值对应的特征向量是线性无关的，另外的特征值之和等于的迹（主对角元素之和）；的行列式值等于的特征值之积。

5.2 重根特征值，代数重数与几何重数

情况A：设矩阵，根据可以求出其特征值为，这里每个特征值的代数重数都是1，也就是分解出来的指数为1，所以每个特征值都对应一个线性无关的特征向量。

情况B：设矩阵，根据可以求出其特征值为，这里特征值为1的代数重数是2，也就是分解出来的指数为2，特征值为2的代数重数为1。解得对应的特征向量为，对应的特征向量为，此时找不到的三个线性无关的特征向量作为三维空间的一组基。

情况C：设矩阵，根据可以求出其特征值为，这里特征值为1的代数重数是2，也就是分解出来的指数为2，特征值为2的代数重数为1。解得对应的特征向量为，对应的特征向量为，此时重根特征值对应的特征向量是两个线性无关的向量的线性组合，可以找到的三个线性无关的特征向量作为三维空间的一组基。

可以发现，重根的特征值对应的线性无关的特征向量个数不确定，我们称能够找到的线性无关的特征向量的个数叫做几何重数，几何重数一定小于等于代数重数（构成的空间维数），这是因为这些特征向量不能表示超过代数重数的维度的空间。因此可以推出，单根特征值只能对应一个线性无关的特征向量，它的几何重数与代数重数相等，都为1。

5.3 相似矩阵

矩阵，，对的映射结果是，设这个结果为。我们一般默认以单位正交基向量作为基底描述向量的坐标，如图所示：

image-20230730190527349

所以，，。

下面用另一组基向量来描述的坐标，设新坐标为，只需要求解即可，这可以看做自变量是的非齐次线性方程组，求解这个方程组，可以得到，。

以为基向量，我们可以表示出这个变换过程为：这里；

同样的变换，如果以为基向量，变换过程为：，我们设变换矩阵为，由于基和坐标都发生了变化，显然这里的，下面求。

由于前面求出了由两组基下的坐标表示，也就是如下等式成立：可以很容易地求出由到的基坐标变换矩阵：设这个基变换矩阵为，为可逆矩阵。因此，以表示的坐标可以由乘以由为基下的坐标：并且可以反过来用来表示：此时，以为基，为变换矩阵的的坐标，代入（a）和（b）式得：观察发现，显然上面式子中中所求的就等于，就此推导出了和的关系。

可以发现，和都可以表示同一个线性变换，只不过由于和选取的基底不同，因此和也不同，但是二者之间可以通过基变换矩阵进行联系。

同样的结论可以扩展到维，假设是维向量，在一组基下的坐标为，在另一组基下的坐标为，则有：假设在下经过的线性变换后，向量变为，对应的坐标就是；同样的变换在下需要经过的线性变换，对应的坐标为，得到如下等式：如何求得与的关系呢？答案就是两组基向量的基变换矩阵。设为从到的基变换矩阵，即，则根据坐标变换公式可以得到：因此，代入（c）式中可以得到：所以，，我们称这个变换为到的相似变换，为的相似矩阵，记做。

可见，和表示的线性空间中相同的线性变换，只不过二者使用的基底不同，因此变换矩阵和也不同。

简单介绍相似矩阵的性质，如果：

  • 的特征值、特征多项式相同，即
  • 的秩相同，
  • 的迹相同
  • 特征值乘积相同，行列式相同，即
  • 若一方可逆，则另一方也可逆，且，同理可得
  • ，这个结论可以推广，即

5.4 相似对角化

相似矩阵的一个重要的应用，就是当如果很复杂，而很简单时，则对应的线性变换可以使用来表达，特别是当为对角矩阵时，可以得到一些非常好的性质。求对应的相似对角矩阵的过程就称为矩阵的相似对角化。

并不是所有的矩阵都可以相似对角化。假设一个维方阵存在相似对角矩阵，则存在可逆矩阵，使得：成立。将上式变形，得到：，

其中：则得到：相当于的每个列向量乘以对角阵的各个元素，对比可以发现每个元素都有如下等式：成立，这就是特征值与特征向量的定义。因此得出结论，如果通过可逆矩阵的相似变换得到了对角阵，则的各个元素是的特征值，的各列就是特征值对应的特征向量。现在问题就转换为求的特征值和特征向量。

显然，阶方阵一定有个特征值，每个特征值都对应各自的特征向量，所以一定存在对角阵和方阵，使得成立。

但是，由于特征值可能有重根，也就不一定具有个线性无关的特征向量，因此不一定可逆，无法得出成立。这说明，必须是可逆矩阵，也就是必须具有个线性无关的特征向量时，才成立，这就是可以相似对角化的充要条件。

5.5 特征值分解

前面介绍了所需的前置知识，下面就非常简单了。特征值分解就是基于相似对角化而来的。假设可以进行相似对角化，那么一定存在对角矩阵和可逆矩阵，使得成立。

将上面的式子变形，即可得到：也就是说，如果一个方阵可以相似对角化，那么可以被分解为三个矩阵相乘的形式，其中为可逆阵，为对角阵，由前面的推导可知，U的每个元素都是的特征值，因此也被称为特征值矩阵，被称为特征向量矩阵。换句话说，可以由的特征值矩阵和的特征向量矩阵的乘积共同表示。

在numpy中，可以通过numpy.linalg.eig计算阶方阵的特征值和特征向量。它返回两个值，这里用w和v表示。前者是特征值，后者是特征值对应的特征向量组成的矩阵。

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np

X = np.array([
    [1, 6, 2],
    [4, 5, 1],
    [7, 1, 1],
])

w, v = np.linalg.eig(X)
# 根据公式可以计算出原矩阵X
origin_X = v @ np.diag(w) @ np.linalg.inv(v)
# 或者计算出与X相似的对角矩阵
diag = np.linalg.inv(v) @ X @ v

需要注意的是，如果输入的矩阵特征向量之间线性相关，则可能无法求逆。如果矩阵特征值都不同，那么特征向量是线性无关的，就可以通过上述公式求得相似对角矩阵。

5.6 应用

根据相似矩阵的性质，我们可以利用特征值分解简化计算。比如，若要计算方阵的幂次，可以将其分解为对角矩阵，计算的幂次，由于对角矩阵的幂次就等于对角元素各自的幂次，计算变得非常简单，然后将计算后的幂次通过可逆阵转换得到的幂次。

下面是一个矩阵求平方根的例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
import scipy

X = np.array([
    [1.25, -0.75],
    [-0.75, 1.25],
])

w, v = np.linalg.eig(X)
# 计算出与X相似的对角矩阵
diag = np.linalg.inv(v) @ X @ v
sqrt_diag = scipy.linalg.sqrtm(diag)
# 计算X的平方根
result = v @ sqrt_diag @ np.linalg.inv(v)
# 直接计算平方根
sqrt = scipy.linalg.sqrtm(X)

这里的scipy.linalg.sqrtm用于计算矩阵的平方根。

6 谱分解

在特征值分解的基础上，如果为实对称矩阵，则它一定可以相似对角化，且其不同特征值对应的特征向量相互正交。因此有如下定理：假设为阶实对称矩阵，则一定存在正交阵和对角矩阵，使得：成立。

将上面的式子变形，即可得到：也就是说，一个实对称矩阵一定可以分解为正交阵和与之相似的对角阵的表示形式，这就是谱分解。

此时程序中就可以用转置代替求逆：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import numpy as np
import scipy

X = np.array([
    [1, 2, 3],
    [2, 1, 4],
    [3, 4, 1],
])

w, v = np.linalg.eig(X)
# 根据公式可以计算出原矩阵X
origin_X = v @ np.diag(w) @ v.T
# 或者计算出与X相似的对角矩阵
diag = v.T @ X @ v