雅可比SVD算法：高精度矩阵分解的经典方法

原创于 2025-10-28 22:57:04 发布 · 629 阅读

23 ·

CC 4.0 BY-SA版权

文章标签：

#矩阵 #线性代数 #SVD #雅可比SVD #Jacobi SVD #人工智能 #矩阵分解

人工智能同时被 3 个专栏收录

290 篇文章

订阅专栏

OTHER

261 篇文章

订阅专栏

Python

258 篇文章

订阅专栏

本文由「大千AI助手」原创发布，专注用真话讲AI，回归技术本质。拒绝神话或妖魔化。搜索「大千AI助手」关注我，一起撕掉过度包装，学习真实的AI技术！

雅可比SVD算法是计算矩阵奇异值分解的一种经典数值方法，以其高精度和天然并行性在科学计算和工程应用中占有重要地位。虽然相比现代SVD算法（如QR-based方法）速度较慢，但雅可比算法在需要高精度的场景和并行计算环境中仍具有不可替代的价值。

1 雅可比SVD算法的历史与背景

雅可比SVD算法源于19世纪卡尔·古斯塔夫·雅可比（Carl Gustav Jacob Jacobi）提出的特征值求解方法，后来被推广到奇异值分解领域。这种算法本质上是通过一系列正交变换（雅可比旋转）逐步将矩阵对角化，从而得到奇异值和奇异向量。

传统的雅可比SVD算法主要有两种变体：

双边雅可比SVD：通过左右两侧交替施加正交变换，直接对原矩阵进行对角化
单边雅可比SVD：通过正交变换将矩阵的列向量逐步正交化，更适合并行计算

本文由「大千AI助手」原创发布，专注用真话讲AI，回归技术本质。拒绝神话或妖魔化。搜索「大千AI助手」关注我，一起撕掉过度包装，学习真实的AI技术！

往期文章推荐:

2 雅可比SVD的数学原理

2.1 核心思想

雅可比SVD算法的核心思想是通过一系列平面旋转矩阵（雅可比旋转矩阵）逐步将矩阵的非对角元素零化。对于任意矩阵A ∈ ℝᵐˣⁿ，算法寻求正交矩阵U和V，使得：

UᵀAV = Σ

其中Σ是对角矩阵，其对角线元素即为A的奇异值。

2.2 雅可比旋转矩阵

雅可比旋转矩阵J(p, q, θ)是一个在(p, q)坐标平面内旋转θ角度的正交矩阵，形式为：

J = [ I   0   0 ]
    [ 0   c  -s ]
    [ 0   s   c ]

其中c = cos(θ)，s = sin(θ)，I是单位矩阵。

2.3 双边雅可比算法步骤

双边雅可比SVD算法的基本迭代步骤如下：

选择枢轴对(p, q)，其中p < q
构造雅可比旋转矩阵J，使得(AJ)ᵀ(AJ)在(p, q)和(q, p)位置变为零
更新矩阵：A ← AJ
累积变换：V ← VJ
对左侧执行类似过程：A ← JᵀA，U ← UJ
重复迭代直到非对角元素足够小

算法收敛后，我们可以这样恢复奇异值分解M=USVᵀ：矩阵V是雅可比旋转矩阵的累积，将转换后矩阵M的列归一化得到矩阵U，而奇异值则由转换后矩阵M列的范数给出。

3 雅可比SVD的算法实现

3.1 基本实现框架

以下是雅可比SVD算法的简化Python实现，展示了算法的核心逻辑：

import numpy as np
import matplotlib.pyplot as plt

def jacobi_svd(A, max_iter=100, tol=1e-8):
    """
    简化的雅可比SVD算法实现
    :param A: 输入矩阵 (m x n)
    :param max_iter: 最大迭代次数
    :param tol: 收敛容差
    :return: U, S, Vt
    """
    m, n = A.shape
    U = np.eye(m)
    V = np.eye(n)
    B = A.copy()

    for iteration in range(max_iter):
        max_off_diag = 0
        p, q = 0, 0

        # 找出最大的非对角元素
        for i in range(min(m, n)):
            for j in range(i + 1, min(m, n)):
                off_diag = abs(B[i, j]) + abs(B[j, i])
                if off_diag > max_off_diag:
                    max_off_diag = off_diag
                    p, q = i, j

        # 收敛检查
        if max_off_diag < tol:
            print(f"算法在{iteration}次迭代后收敛")
            break

        # 计算旋转角度
        if B[p, p] == B[q, q]:
            theta = np.pi / 4
        else:
            theta = 0.5 * np.arctan(2 * B[p, q] / (B[p, p] - B[q, q]))

        c, s = np.cos(theta), np.sin(theta)

        # 构造旋转矩阵
        J = np.eye(n)
        J[p, p] = c
        J[p, q] = -s
        J[q, p] = s
        J[q, q] = c

        # 更新矩阵
        B = B @ J
        V = V @ J

        # 左侧更新
        J_left = np.eye(m)
        J_left[p, p] = c
        J_left[p, q] = -s
        J_left[q, p] = s
        J_left[q, q] = c

        B = J_left.T @ B
        U = U @ J_left

    # 提取奇异值
    S = np.zeros((m, n))
    for i in range(min(m, n)):
        S[i, i] = B[i, i]

    return U, S, V.T

# 算法测试
print("🧪 雅可比SVD算法测试")
print("=" * 50)

# 创建一个测试矩阵
A_test = np.array([[3, 2, 1],
                   [1, 4, 2],
                   [2, 1, 3]], dtype=float)

print("输入矩阵A:")
print(A_test)

# 使用我们的雅可比SVD实现
U_jacobi, S_jacobi, Vt_jacobi = jacobi_svd(A_test)

print("\n雅可比SVD结果:")
print("U矩阵:")
print(U_jacobi)
print("奇异值:", np.diag(S_jacobi))
print("Vt矩阵:")
print(Vt_jacobi)

# 与NumPy官方实现对比
U_np, S_np, Vt_np = np.linalg.svd(A_test)
print("\nNumPy SVD结果对比:")
print("官方奇异值:", S_np)

# 重构误差计算
A_reconstructed = U_jacobi @ S_jacobi @ Vt_jacobi
reconstruction_error = np.linalg.norm(A_test - A_reconstructed)
print(f"\n重构误差: {reconstruction_error:.2e}")