数值分析：矩阵奇异值分解

1. 奇异值分解(SVD)

(1)奇异值分解

已知矩阵\(\bm{A} \in \R^{m \times n}\), 其奇异值分解为：

\[\bm{A} = \bm{U}\bm{S}\bm{V}^T \]

其中\(\bm{U} \in \R^{m \times m}\)，\(\bm{V} \in \R^{n \times n}\)是正交矩阵，\(\bm{S} \in \R^{m \times n}\)是对角线矩阵。\(\bm{S}\)的对角线元素\(\bm{s}_1, \bm{s}_2,..., \bm{s}_{\min(m,n)}\)是矩阵的奇异值。

(2) 奇异值分解的求解

而求矩阵的奇异值的算法非常简单，对于实数域下的矩阵\(A\)，我们只需要求\(A^TA\)的特征值和特征向量。其特征向量归一化后即右奇异向量\(\bm{v}_1,\bm{v}_2,...,\bm{v}_n\)，其特征值开根号即对应的奇异值\(\bm{s}_1, \bm{s}_2,..., \bm{s}_{\min(m,n)}\)。 然后由等式

\[\bm{A}\bm{v}_1=s_1\bm{u}_1, \\ \bm{A}\bm{v}_2=s_2\bm{u}_2, \\ ..., \\ A\bm{v}_{\min(m,n)} = \bm{s}_{\min(m,n)}\bm{u}_{\min(m,n)}\]

依次计算出相应的\(\bm{u}_i\)向量的值。
至于特征值的计算，采用 QR 算法，此处不予介绍，这里可以直接调用 np.linalg.eig()函数实现。以下给出奇异值计算代码实例(此处仅为知识演示，具体的工业级别的奇异值计算算法要复杂得多，参考 Golub 与 Van Loan《矩阵计算》)

import numpy as np
def svd(A):
    eigen_values, eigen_vectors = np.linalg.eig(A.T.dot(A))
    singular_values = np.sqrt(eigen_values)
    #这里奇异值要从大到小排序，特征向量也要随之从大到小排
    val_vec = [] #存储奇异值-特征向量对
    for i in range(len(eigen_values)):
        val_vec.append((singular_values[i], eigen_vectors[:, i]))
    val_vec.sort(key = lambda x:-x[0])
    singular_values = [ pair[0] for pair in val_vec]
    eigen_vectors = [ pair[1] for pair in val_vec]

    # 在计算左奇异向量之前，先要对右奇异向量也就是特征向量组成的基正交化
    # 不过linalg.eig返回的是已经正交化的，这一步可省略

    # 由等式Avi = siui(vi是右奇异向量, ui是左奇异向量)
    # 依次计算左奇异向量
    U = np.zeros((A.shape[0], A.shape[1]))
    for i in range(A.shape[1]):
        u = A.dot(eigen_vectors[i])/singular_values[i]
        U[:, i] = u
    # 给U加上标准正交基去构造R3的基
    for i in range(A.shape[1], A.shape[0]):
        basis = np.zeros((A.shape[0], 1))
        basis[i] = 1
        U = np.concatenate([U, basis], axis=1)
    eigen_vectors = [vec.reshape(-1, 1) for vec in eigen_vectors]
    eigen_vectors = np.concatenate(eigen_vectors, axis=1)
    return U, singular_values, eigen_vectors

if __name__ == '__main__':
    # 例一：普通矩阵
    A = np.array(
        [
            [0, 1],
            [0, -1]
        ]
    )
    # 例二：对称矩阵
    # A = np.array(
    #     [
    #         [0, 1],
    #         [1, 3/2]
    #     ]
    # )
    U, S, V = svd(A)
    print("我们实现的算法结果：")
    print(U, "\n", S, "\n", V)
    print("\n")
    print("调用库函数的计算结果：")
    # 调用api核对
    U2, S2, V2 = np.linalg.svd(A)
    print(U2, "\n", S2, "\n", V2)

对普通矩阵\(\left(\begin{matrix}0 & 1 \\0 & -1\end{matrix}\right)\)运行该算法的结果为：

我们实现的算法结果：
[[ 0.70710678  0.70710678]
 [-0.70710678 -0.70710678]] 
 [1.4142135623730951, 0.0] 
 [[0. 1.]
 [1. 0.]]


调用库函数的计算结果：
[[