搬砖笔记

SourceLizi's Blog

0%

矩阵数值算法总结

发表于 分类于

LU分解

  • 线性方程组求解稳定性判断：条件数
    通常情况下，所有求解形如的方程组算法的数值稳定性与矩阵d本身特性是密切相关的，稳定性的衡量可以用矩阵的条件数（通常实际计算中会计算奇异值的最大值和最小值之比表示条件数）

    矩阵条件数越大代表矩阵越接近病态，而条件数越接近1则矩阵越接近良态。病态矩阵对数值变化十分敏感，当病态矩阵的系数发生微小变化时，其求解结果会发生巨大变化。

LU分解主要用于线性方程组的求解，它是把高斯消元法写成了矩阵的形式表示计算，首先朴素的高斯消元法是用增广矩阵第行的方程依次消去第列对角线下方的元素，将矩阵变换为上三角阵后，再从最后一行开始依次解出的每个元素。首先线性方程组的增广矩阵可以表示为：

消去对角线下方第一列的元素可以得到

以此类推消去其他列对角线下方的元素，得到上三角矩阵

这个时候可以很容易计算出最后一个元素，再将代入到上一条方程中可以计算出，按照这个顺序依次操作即可解出所有元素的值
基于上面的过程可以用矩阵表示，将矩阵A分解为下三角矩阵L和上三角矩阵U，首先将上面的变换过程（即行操作）写为矩阵，变换结果的上三角矩阵写为，表示为。其次将作为为未知数，得到新的线性方程组，先求解再重新求解即可。
上述LU分解存在除数为0或者由于过于接近0导致数值误差较大，为了解决这个问题，我们可以交换行来尽可能令行中最大的元素出现在对角线上，该改进被称为部分选主元(partial pivoting)。这一改进通过增加一个矩阵交换矩阵中的行实现，即改进后的LU矩阵分解表示为。

LU分解的Python实现如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def gaussian_elimination(A,b):
    n = len(A)
    assert len(b) == n, "A and b size not match"
    for l in A:
        assert len(l) == n, "Require square matrix"
    for i in range(n-1):
        ##搜索最大值主元
        for j in range(i+1,n):
            if abs(A[j][i]) > abs(A[i][i]):
                A[j], A[i] = A[i], A[j]
                b[j], b[i] = b[i], b[j]
        ##消元步骤
        for j in range(i+1,n):
            r = A[j][i]/A[i][i]
            A[j][i] = 0.0
            for k in range(i+1,n):
                A[j][k] -= r*A[i][k]
            b[j] -= r*b[i]
    ##回代方程求出x
    x = [0.0 for _ in range(n)]
    for i in range(n-1,-1,-1):
        for j in range(i+1,n):
            b[i] -= A[i][j]*x[j]
        x[i] = b[i]/A[i][i]
    return x

另外若为实对称正定矩阵，上述分解可以用Cholesky分解代替，将正定对称矩阵分解为下三角矩阵与其转置的乘积即

迭代法（解线性方程组）

  • Jacobi法
    Jacobi法是一种矩阵形式的不动点迭代法，用于求解线性方程组。简单来讲Jacobi方法先将第k个未知数用除了这个未知数以外的其他未知数表示，以此作为不定点不动点迭代法的迭代式进行迭代。相比于矩阵分解等方法求解，其在求解大规模的稀疏矩阵中具有明显的优势。针对求解方程组，假设矩阵的对角元素矩阵，下三角元素矩阵和上三角元素矩阵分别为，则Jacobi法的迭代式可以表示为

    与标量版本的不动点迭代法一样，这里迭代也有收敛条件。Jacobi法收敛的充分必要条件是矩阵A的谱半径（特征值绝对值的最大值）小于1，另外一种常用的判断依据（充分条件）是矩阵是严格对角占优(strictly diagonally dominant)，即对于
  • Gauss–Seidel法
    与Jacobi法基本相同，但是每计算其中一个变量的值会用上本轮已经更新后的其他变量的值，而不是像Jacobi法那样每轮迭代完之后更新的值。从代码上来讲这种方法可以直接使用同一块内存储存，实现起来更加简单方便。
  • SOR(Successive Over-Relaxation)法
    类似动量法，在每一步计算中叠加上一次计算的值，迭代式为

    其中 为SOR法在第k步计算出的下一步更新值，设置会让迭代过程略微过冲，达到加速收敛的目的

Jacobi法和Gauss–Seidel法的算法Python实现如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def jacobi_iter(A,b, max_iter:int = 100, tol: float = 1e-5):
    ##D^-1(b-(L+U)x) -> D^-1(b-Ax)+x
    n = len(A)
    assert len(b) == n, "A and b size not match"
    for l in A:
        assert len(l) == n, "Require square matrix"
    x = [0.0 for _ in range(n)]
    for _ in range(max_iter):
        xk = x.copy()
        max_d = float("inf")
        for i in range(n):
            s = 0.0
            for j in range(n):
                s += A[i][j]*x[j]
            delta = (b[i] - s)/A[i][i]
            max_d = max(max_d, abs(delta))
            xk[i] += delta
        x = xk
        if max_d < tol:
            break
    return x
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def gauss_seidel_iter(A,b, max_iter:int = 100, tol: float = 1e-5):
    ##D^-1(b-(L+U)x) -> D^-1(b-Ax)+x
    n = len(A)
    assert len(b) == n, "A and b size not match"
    for l in A:
        assert len(l) == n, "Require square matrix"
    x = [0.0 for _ in range(n)]
    for _ in range(max_iter):
        max_d = float("inf")
        for i in range(n):
            s = 0.0
            for j in range(n):
                s += A[i][j]*x[j]
            delta = (b[i] - s)/A[i][i]
            max_d = max(max_d, abs(delta))
            x[i] += delta
        if max_d < tol:
            break
    return x

QR 分解

QR分解指的是将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积，其实现方法可以通过GS(Gram–Schmidt)正交化或Householder变换

GS(Gram–Schmidt)正交化

Gram–Schmidt正交化用于对一组矢量进行正交化。给定一组向量作为输入，其目标是这组矢量所张成的子空间找到一组正交基表示。或者说说，给定n个线性无关的输入向量，Gram–Schmidt正交化可以计算出n个互相垂直的单位向量，这些向量张成的子空间与输入向量所张成的子空间相同。
Gram–Schmidt正交化的原理是在输入的线性无关矢量中，不断减去已经选出的正交矢量在新的矢量上的投影并归一化后获得下一个矢量。其主要过程如下：
假设输入n个线性无关的列矢量为，首先选择第一个矢量并归一化，获得第一个正交基

在选择第二个矢量中减去第一个正交基在上面的投影，并且归一化结果获得第一个正交基

可以发现和是正交的，因此和也是正交的。依此类推，计算出的第个正交基为

定义，上述正交化结果可以重新表示每个矢量，写成矩阵形式为

矩阵为正交矩阵，而组成的矩阵为上三角矩阵，这样就完成了原矩阵的QR分解。
考虑到存在矢量元素个数的情况，我们可以在矩阵末尾补上个矢量，而在矩阵末尾补0，这样能够令分解出来的为与原矩阵行数相同的方阵而大小与原矩阵相同。

另外一点，在实现QR分解的时候为了避免不正交的带来的数值误差，在计算投影时，不使用原来的矢量，而是使用之前已经移除掉部分矢量投影的临时值，即每一步减去投影的矢量会立即影响到下一步投影的计算。基于改进的GS(Modified Gram–Schmidt)正交化的QR分解Python实现如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import math
def qr_factorization_gs(A):
    m = len(A)
    n = len(A[0]) if m > 0 else 0
    for l in A:
        assert len(l) == n, "Not a matrix"
    R = [[0.0 for _i in range(n)] for _j in range(n)]
    Q = [[0.0 for _i in range(n)] for _j in range(m)]
    for j in range(n):
        y = [A[k][j] for k in range(m)]
        for i in range(j):
            for k in range(m):
                R[i][j] += Q[k][i]*y[k]
            for k in range(m):
                y[k] -= R[i][j]*Q[k][i]
        R[j][j] = math.sqrt(sum([yk*yk for yk in y]))
        for k in range(m):
            Q[k][j] = y[k] / R[j][j]
    return Q, R

Householder变换

Householder变换又称初等反射变换，该变换是一种线性变换，可以用初等反射矩阵表示。这个矩阵可以想象为把某个向量照着某个平面进行反射翻转，将出射/入射进行变换。假设出射和入射矢量为，且，则平面法向量为，Householder变换矩阵可以表示为，该矩阵满足：

  • 对称性和正交性：
  • 镜像变换：

基于这个矩阵的性质，我们可以设计一种Householder变换将矩阵中的某一列后行元素置零，这样既能够以Householder变换矩阵的正交性获得正交矩阵，又可以置零获得上三角矩阵。我们可以将列的第一个元素置为矢量模长而其余元素置为0，这样的变换依次对原矩阵的右下角分块进行，即可获得一个上三角矩阵。以4*3大小的矩阵为例，可以设计三个Householder变换将矩阵变换为上三角矩阵，矩阵

基于HouseHolder变换的QR分解Python实现如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import math
##矩阵乘法
def mat_mul(a,b):
    m_a,n_a = len(a),len(a[0])
    n_b,k_b = len(b),len(b[0])
    assert n_a == n_b, "Matrix size must match"
    c = [[0.0 for _j in range(k_b)] for _i in range(m_a)]
    for i in range(m_a):
        for j in range(k_b):
            for k in range(n_a):
                c[i][j] += a[i][k]*b[k][j]
    return c

def qr_factorization_householder(A):
    m = len(A)
    n = len(A[0]) if m > 0 else 0
    for l in A:
        assert len(l) == n, "Not a matrix"
    Q = [[1.0 if i == j else 0.0 for i in range(m)]
         for j in range(m)]
    for i in range(min(m-1,n)):
        v0 = math.sqrt(sum([A[k][i]*A[k][i] for k in range(i,m)]))
        v = [-A[k][i] for k in range(i,m)]
        v[0] += -math.copysign(v0,v[0])
        vn2 = sum([p*p for p in v])
        H = [[1.0 if i == j else 0.0 for i in range(m)] for j in range(m)]
        for ib in range(i,m):
            for jb in range(i,m):
                H[ib][jb] -= 2*v[ib-i]*v[jb-i]/vn2
        Q = mat_mul(Q,H)
        A = mat_mul(H,A)
    return Q, A

应用

  • 求解最小二乘问题
    最小二乘问题可以表述为

    一种解析解的方法是计算伪逆，即，但是求逆也可能出现病态的问题，这里可以采用QR分解作为数值更加稳定的方法。对进行完整QR分解，原式可以表示为，由于为正交矩阵，因此原式等价为。令，则

    观察的组成可以发现，后面到之间的元素均恒等于，同时当前n个元素等于零时最小。故以上面n*n的方阵作为系数矩阵，前个元素作为方程右边的矢量，求解获得解

特征值求解

幂迭代法

幂迭代法(Power Iteration)是用于求解最大特征值（如果存在）及其对应特征向量的一种数值方法，它的思路非常简单：对于给定的一个任意矢量，用待求解矩阵的n次方左乘，当n足够大时其结果即为最大特征值对应的特征矢量，求解对应特征值可以通过瑞利商

与其对应的还有反幂迭代法(Inverse Power Iteration)，即对求逆的矩阵进行幂迭代计算其最小特征值。

基于QR分解

在上面幂迭代法的基础上，用QR分解的方法同时对m个正交矢量进行幂迭代，可以获得m*m矩阵的特征向量矩阵。这个方法要求矩阵为实对称矩阵（对称矩阵的特征矢量都是正交的）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def qr_eig(A, max_iter: int = 100):
    m = len(A)
    for l in A:
        assert len(l) == m, "Require square matrix"
    Q = [[1.0 if i == j else 0.0 for i in range(m)]
         for j in range(m)]
    Q_res = Q
    R = A
    for _ in range(max_iter):
        Q, R = qr_factorization_householder(mat_mul(R,Q))
        Q_res = mat_mul(Q_res,Q)
    mat_eig = mat_mul(R,Q)
    eig = [mat_eig[i][i] for i in range(m)]
    return eig, Q_res

在这个基础上可以改进得到求解任意方阵特征值的Shifted QR Algorithm，在QR迭代前为矩阵特征值增加一个偏置，使得调整后的矩阵特征值接近0，在迭代结束后去掉这个偏置，这样可以加快迭代收敛。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def qr_shifted_eig(A, tol : float = 1e-8):
    m = len(A)
    for l in A:
        assert len(l) == m, "Require square matrix"
    eig = [0.0 for _ in range(m)]
    for i in range(m-1,0,-1):
        while max([abs(A[i][j]) for j in range(i)]) > tol:
            mu = A[i][i]
            M1 = [[A[j][k] - (mu if j == k else 0.0)
                    for k in range(i+1)] for j in range(i+1)]
            Q, R = qr_factorization_householder(M1)
            A = mat_mul(R,Q)
            for k in range(i+1):
                A[k][k] += mu
        eig[i] = A[i][i]
        A = [[A[j][k] for k in range(i)] for j in range(i)]
    eig[0] = A[0][0]
    return eig

上面这个算法主要的问题时没有处理复数特征值。迭代求解的过程实际上最好的结果是Schur form矩阵（上三角矩阵上可能有2*2的块，代表复数特征值），在算法迭代过程中矩阵A可能会出现无法收敛的问题，这个时候要检查是否存在2*2块并做出处理。另外为了处理多重复数特征值（特征方程存在重根）还可能需要使用HouseHolder变换将原矩阵先转换为为Hessenberg形式的矩阵。

经典Jacobi方法

Jacobi方法是一种迭用于求解实对称矩阵的全部特征值和特征向量的数值迭代。其核心思想是通过不断构造正交的旋转矩阵消去非对角元素，使矩阵逐步趋近于对角阵，则对角线元素即为特征值，对应的正交矩阵列向量为特征向量。该过程表示为

其中矩阵是一个Givens旋转矩阵，用于调整原矩阵两行两列的对应四个位置的值，它可以表示为

其中四个对角的位置分别为，假设调整后的四个元素分别为，则

令，即，解得。这样迭代一次两个对角上的元素即被置为0.
算法设计上，每次搜索绝对值最大的非对角元素消去。由于每次迭代两行两列元素均会受到影响，因此可能出现之前迭代为0的元素在后面的迭代里面又变为非零，故算法可能存在某个非对角位置被重复迭代的情况。基于Jacobi方法的特征值和特征向量求解算法Python实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def jacobi_eig(A, max_iter: int = 100,tol : float = 1e-8):
    m = len(A)
    for l in A:
        assert len(l) == m, "Require square matrix"
    V = [[1.0 if i == j else 0.0 for i in range(m)]
         for j in range(m)]
    for _ in range(max_iter):
        p, q, max_val = 0, 1, abs(A[0][1])
        for i in range(m):
            for j in range(i+1,m):
                if abs(A[i][j]) > max_val:
                    p, q, max_val = i, j, abs(A[p][q])
        if max_val < tol:
            break
        theta = 0.5 * math.atan2(2*A[p][q], A[p][p] - A[q][q])
        c, s = math.cos(theta), math.sin(theta)

        ## 矩阵J比较稀疏，如果要优化效率乘法可以拆开来直接计算
        J = [[1.0 if i == j else 0.0 for i in range(m)]
         for j in range(m)]
        J[p][p], J[p][q] = c, -s
        J[q][p], J[q][q] = s, c

        J_T = [[1.0 if i == j else 0.0 for i in range(m)]
         for j in range(m)]
        J_T[p][p], J_T[p][q] = c, s
        J_T[q][p], J_T[q][q] = -s, c
        A = mat_mul(mat_mul(J_T,A),J)
        V = mat_mul(V,J)
    eig = [A[i][i] for i in range(m)]
    return eig, V

奇异值分解

假设矩阵是一个的矩阵，则矩阵A的SVD定义为：

其中分别为和的酉矩阵，是一个的对角阵，对角线上的元素为奇异值（矩阵特征值的算数平方根）。相比于特征值分解，奇异值分解能够处理任意大小的矩阵而不仅仅是方阵，因此更加通用。

单边Jacobi

与奇异值分解里面的Jacobi方法几乎一样，只不过对于矩阵通过右乘若干个Givens旋转矩阵，使得最终结果的矩阵每一列均正交，

则最终矩阵的列的模为奇异值，归一化后的列为左奇异矩阵，累乘的结果为右奇异矩阵的转置。假设要正交化的原矩阵中的两列矢量分别，只考虑输出的两列

要令两列正交则

可以得到类似特征值分解求解过程的推导结果，由于每次右乘都只能正交化两列向量，后面的计算可能会影响到前面计算两列的正交性，因此算法同样需要迭代。用Python编写的基于单边Jacobi的SVD算法实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import math
def jacobi_svd(A, max_iter: int = 100,tol : float = 1e-8):
    m = len(A)
    n = len(A[0]) if m > 0 else 0
    for l in A:
        assert len(l) == n, "Not a matrix"
    V = [[1.0 if i == j else 0.0 for i in range(n)]
         for j in range(n)]
    a_col2 = [0.0 for _ in range(n)]
    for _ in range(max_iter):
        converged = True
        for i in range(n-1):
            for j in range(i+1, n):
                a_col2[i] = sum([A[k][i]*A[k][i] for k in range(m)])
                a_col2[j] = sum([A[k][j]*A[k][j] for k in range(m)])
                ai_aj = sum([A[k][j]*A[k][i] for k in range(m)])
                if abs(ai_aj) > tol * math.sqrt(a_col2[i] * a_col2[j]):
                    converged = False
                    theta = 0.5 * math.atan2(2*ai_aj, a_col2[j] - a_col2[i])
                    c, s = math.cos(theta), math.sin(theta)
                    for k in range(m):
                        A[k][i], A[k][j] = c*A[k][i] - s*A[k][j], s*A[k][i] + c*A[k][j]
                        V[k][i], V[k][j] = c*V[k][i] - s*V[k][j], s*V[k][i] + c*V[k][j]
        if converged:
            break
    s = [math.sqrt(a_col2[i]) for i in range(n)]
    U = [[A[i][j]/s[j] for i in range(n)]
         for j in range(m)]
    return U, s, V

Golub-Kahan方法

Golub-Kahan方法类似上面用HouseHolder变换完成QR分解，它将原矩阵分解为两个正交矩阵和一个双对角矩阵，再用其他方法分解双对角阵为对角阵和两个正交矩阵，组合这些正交矩阵得到SVD分解结果

应用

  • 求解最小二乘问题
    对于经典的最小二乘，SVD可以很方便地表示出结果
    另外一种特殊情况是，原问题可以表述为

    其中是为了取基础解系里唯一一个解而增加的约束。将SVD结果代入得

    令，由可知的模长同样等于1，原优化的目标函数可以写为

    当且仅当时目标函数取最小值，此时，即取正交矩阵最小特征值对应的列向量。
  • 求矩阵的秩
    矩阵的秩定义为非零奇异值的数量，因此求解出非零奇异值的数量即可
  • 低秩近似
    根据Eckart-Young-Mirsky定理，求秩为矩阵使其最佳近似原矩阵（范数最小化），只需要将奇异值分解之后取前大的奇异值重构该矩阵即可
  • 正交矩阵近似
    机器视觉里面求解出的旋转矩阵不一定正交，我们可以采用矩阵近似的方法求一个与结果最佳逼近的正交矩阵，对进行原矩阵SVD分解，则该矩阵的最佳近似正交矩阵为

    推导参考这里

参考

  • Timothy Sauer - Numerical Analysis(Third Edition)