The article describing MATLAB SVD: principles, algorithms, and code implementation is complete.
Here is the full article:
MATLAB SVD：原理、算法与代码实现
奇异值分解（Singular Value Decomposition，简称 SVD）是线性代数中一个极其重要且应用广泛的矩阵分解技术。它能够将任意实数或复数矩阵分解为三个特殊矩阵的乘积，揭示了矩阵深层的几何和代数结构。SVD 不仅在理论上具有优雅的数学性质，在实际应用中更是无处不在，例如数据压缩、降维（如主成分分析 PCA）、图像处理、推荐系统、自然语言处理以及解决最小二乘问题等。
本文将深入探讨 SVD 的核心原理，介绍其计算算法，并通过 MATLAB 代码演示如何利用内置函数以及如何从头实现 SVD，旨在帮助读者全面理解并掌握 SVD 这一强大工具。
一、 SVD 的数学原理
奇异值分解将任意一个 m × n 的矩阵 A 分解为以下形式：
A = U Σ Vᵀ
其中：
* U 是一个 m × m 的正交矩阵（Orthogonal Matrix），其列向量被称为 A 的左奇异向量。U 满足 UᵀU = I。
* Σ (Sigma) 是一个 m × n 的对角矩阵（Diagonal Matrix），其对角线上的元素 σ_i 被称为 A 的奇异值 (Singular Values)。奇异值是非负的，通常按降序排列：σ_1 ≥ σ_2 ≥ ... ≥ σ_r > 0，其中 r 是矩阵 A 的秩。对角线以外的元素均为零。
* V 是一个 n × n 的正交矩阵，其列向量被称为 A 的右奇异向量。V 满足 VᵀV = I。Vᵀ 是 V 的转置。
SVD 与特征值分解的关系
SVD 与特征值分解（Eigenvalue Decomposition）紧密相关。对于一个方阵 A，其特征值分解为 A = P Λ P⁻¹。然而，特征值分解只能应用于方阵，且并非所有方阵都可以被对角化。SVD 则适用于任何形状的矩阵。
这种联系体现在以下几个方面：
-
右奇异向量 (
V) 与AᵀA的特征向量：
将A = U Σ Vᵀ代入AᵀA，我们得到：
AᵀA = (U Σ Vᵀ)ᵀ (U Σ Vᵀ) = V Σᵀ Uᵀ U Σ Vᵀ = V Σᵀ Σ Vᵀ由于
Σ是对角矩阵，ΣᵀΣ也是对角矩阵，其对角线元素为σ_i²。
因此，AᵀA的特征值是A的奇异值的平方 (σ_i²)，而AᵀA的特征向量就是A的右奇异向量 (V的列向量)。 -
左奇异向量 (
U) 与AAᵀ的特征向量：
类似地，将A = U Σ Vᵀ代入AAᵀ，我们得到：
AAᵀ = (U Σ Vᵀ) (U Σ Vᵀ)ᵀ = U Σ Vᵀ V Σᵀ Uᵀ = U Σ Σᵀ UᵀΣΣᵀ也是对角矩阵，其对角线元素同样是σ_i²。
因此，AAᵀ的特征值也是A的奇异值的平方 (σ_i²)，而AAᵀ的特征向量就是A的左奇异向量 (U的列向量)。
值得注意的是，尽管 AᵀA 和 AAᵀ 的非零特征值相同（均为奇异值的平方），但它们的特征向量（即 V 和 U）通常是不同的。
二、 SVD 的计算算法
虽然实际应用中我们通常会使用高度优化的库函数来计算 SVD，但理解其背后的算法原理对于深入掌握 SVD 至关重要。这里介绍一种基于特征值分解的 SVD 求解方法，这也是许多从零开始实现 SVD 的基础。
假设我们有一个 m × n 的矩阵 A。
算法步骤：
-
计算右奇异向量 (
V) 和奇异值 (Σ)：- 首先，计算
A的转置与其自身的乘积：C = AᵀA。这是一个n × n的对称矩阵。 - 对矩阵
C进行特征值分解。C的特征向量构成了矩阵V的列向量（即右奇异向量）。 C的特征值λ_i将是非负的。取这些特征值的平方根σ_i = √λ_i，它们就是矩阵A的奇异值。- 将奇异值
σ_i降序排列，并相应地重新排列V的列，使得V的第i列对应于第i大的奇异值。 - 构建对角矩阵
Σ，将降序排列的奇异值放在其主对角线上。Σ的大小为m × n。
- 首先，计算
-
计算左奇异向量 (
U)：- 一旦有了奇异值
σ_i和右奇异向量v_i，我们可以利用关系Av_i = σ_i u_i来求解左奇异向量u_i。 - 对于每个非零奇异值
σ_i，对应的左奇异向量u_i可以通过u_i = Av_i / σ_i计算得到。 - 如果
m > n或者A的秩小于min(m, n)，则通过上述方法得到的u_i不足以构成完整的U矩阵（m × m的正交矩阵）。 - 为了得到完整的
U矩阵，我们可以计算AAᵀ的特征向量。AAᵀ的特征值也是A的奇异值的平方，其特征向量就是U的列向量（左奇异向量）。需要注意的是，通过AAᵀ得到的特征向量的符号可能与通过Av_i = σ_i u_i关系计算出的u_i的符号不一致。因此，在实践中，通常会使用Av_i = σ_i u_i来确定前rank(A)个u_i的方向，并对AAᵀ得到的U的相应列进行符号校正以确保A = U Σ Vᵀ成立。
- 一旦有了奇异值
数值稳定性考虑：
在实际数值计算中，直接使用 AᵀA 或 AAᵀ 来计算特征值可能会导致数值精度问题，特别是在处理病态矩阵（ill-conditioned matrices）时。这是因为计算 AᵀA 会使得条件数平方，放大误差。因此，专业的 SVD 算法（例如 Golub-Kahan 算法或 QR 迭代法）通常会避免直接计算 AᵀA 或 AAᵀ，而是采用更稳定的迭代方法。然而，对于教学和理解原理而言，基于特征值分解的方法是一个很好的起点。
三、 MATLAB 中的 SVD 实现
MATLAB 作为一个强大的数值计算环境，提供了高效且易于使用的内置函数来执行奇异值分解。
1. 内置 svd 函数
MATLAB 的 svd 函数是计算矩阵奇异值分解的首选工具。它基于高度优化的数值算法（例如，通常采用 QR 迭代算法或相关的变体），确保了计算的准确性和效率。
svd 函数有几种调用形式：
s = svd(A): 返回一个包含矩阵A的奇异值向量s。这些奇异值以列向量形式返回，并按降序排列。[U, S, V] = svd(A): 返回三个矩阵：U：一个m × m的正交矩阵，其列是A的左奇异向量。S：一个m × n的对角矩阵，其对角线元素是A的奇异值，按降序排列。V：一个n × n的正交矩阵，其列是A的右奇异向量。- 此形式使得
A = U * S * V'成立。
[U, S, V] = svd(A, 'econ')(Economy Size SVD):- 当
m > n时，U将是一个m × n的矩阵，S是一个n × n的方阵，V是一个n × n的矩阵。此时U的列只包含前n个左奇异向量，S只包含前n个奇异值。 - 当
m <= n时，结果与[U, S, V] = svd(A)相同。 - 这种”经济型”分解在许多应用中非常有用，因为它减少了计算量和存储空间，尤其是在
m远大于n的情况下。它仍然满足A = U * S * V'。
- 当
示例：
“`matlab
% 定义一个示例矩阵
A = [1 2 3; 4 5 6; 7 8 9; 10 11 12]; % 4×3 矩阵
% 计算奇异值
s = svd(A);
disp(‘奇异值 (s):’);
disp(s);
% 计算完整的 SVD
[U, S, V] = svd(A);
disp(‘完整的 SVD (U, S, V):’);
disp(‘U = ‘); disp(U);
disp(‘S = ‘); disp(S);
disp(‘V = ‘); disp(V);
% 计算经济型 SVD
[U_econ, S_econ, V_econ] = svd(A, ‘econ’);
disp(‘经济型 SVD (U_econ, S_econ, V_econ):’);
disp(‘U_econ = ‘); disp(U_econ);
disp(‘S_econ = ‘); disp(S_econ);
disp(‘V_econ = ‘); disp(V_econ);
% 验证分解结果
A_reconstructed = U * S * V’;
fprintf(‘原始矩阵 A 与完整 SVD 重构矩阵的差的范数: %e\n’, norm(A – A_reconstructed));
A_reconstructed_econ = U_econ * S_econ * V_econ’;
fprintf(‘原始矩阵 A 与经济型 SVD 重构矩阵的差的范数: %e\n’, norm(A – A_reconstructed_econ));
“`
运行上述代码，你会发现 norm(A - A_reconstructed) 和 norm(A - A_reconstructed_econ) 的结果都非常接近于零，这验证了 SVD 分解的正确性。
2. 从零实现 mySVD 函数
为了更好地理解 SVD 的底层原理，我们可以尝试基于特征值分解手动实现一个简化的 SVD 函数。以下是 mySVD 函数的代码，它遵循了前面介绍的基于 AᵀA 和 AAᵀ 特征值分解的算法：
“`matlab
function [U, S, V] = mySVD(A)
% mySVD Custom implementation of Singular Value Decomposition (SVD)
% [U, S, V] = mySVD(A) computes the SVD of matrix A such that A = USV’.
% This implementation derives V and S from A’*A and then U from A, V, and S
% (by aligning with AA’ eigenvectors for a full orthonormal basis).
[m, n] = size(A);
% --- 步骤 1: 从 A'*A 计算 V (右奇异向量) 和奇异值 (S 的对角线元素) ---
% A'*A 是一个 n x n 的对称矩阵。它的特征向量构成了 V，其特征值是奇异值的平方。
ATA = A' * A;
% 对 A'*A 进行特征值分解
% D_V 包含特征值 (对角线元素)，V_raw 包含对应的特征向量 (列向量)
[V_raw, D_V] = eig(ATA);
% 提取特征值并按降序排序
eigenvalues_ATA = diag(D_V);
[sorted_eigenvalues, idx] = sort(eigenvalues_ATA, 'descend');
% 奇异值是特征值的非负平方根。
% 使用 max(0, ...) 来处理由于数值精度可能导致的微小负特征值。
singular_values = sqrt(max(0, sorted_eigenvalues));
% 根据排序后的特征值重新排列 V_raw 的列，得到最终的 V
V = V_raw(:, idx);
% --- 步骤 2: 构建对角矩阵 S ---
% S 是一个 m x n 的矩阵，其主对角线上是奇异值。
S = zeros(m, n);
min_dim = min(m, n);
% 将计算出的奇异值放置在 S 的对角线上
S(1:min_dim, 1:min_dim) = diag(singular_values(1:min_dim));
% --- 步骤 3: 计算 U (左奇异向量) 并确保符号一致性 ---
% U 是一个 m x m 的正交矩阵。
% 前 'rank(A)' 列的 U 可以通过 U_i = (A * V_i) / sigma_i 计算。
% 对于剩余的列（如果 m > rank(A)），它们构成了 A' 的零空间的正交基。
% 为了得到完整的 U 矩阵并确保正交性，我们可以使用 AA' 的特征向量。
AAT = A * A';
[U_raw, D_U] = eig(AAT);
eigenvalues_AAT = diag(D_U);
[~, idx_U] = sort(eigenvalues_AAT, 'descend');
U_full = U_raw(:, idx_U); % 这是 U 的一个候选，但其列的符号可能是任意的。
% 根据非零奇异值确定矩阵的有效秩。
% 使用一个小的容差来处理浮点数精度问题。
tolerance = max(m, n) * eps(max(singular_values));
rank_A = sum(singular_values > tolerance);
% 校正 U_full 的列的符号，使其与通过 A*V/sigma 关系得到的方向一致。
% 这确保了 A = U*S*V' 的正确性。
for i = 1:rank_A
% 计算 U(:,i) 的预期方向，基于 A*V(:,i) / sigma_i
expected_U_dir = (A * V(:, i)) / singular_values(i);
% 如果点积为负，表示方向相反，则翻转 U_full 对应列的符号。
if dot(U_full(:, i), expected_U_dir) < 0
U_full(:, i) = -U_full(:, i);
end
end
U = U_full; % 最终的符号对齐的 U 矩阵。
end
“`
代码解释：
-
初始化与预处理：
[m, n] = size(A);获取输入矩阵A的维度。
-
计算
V和奇异值：ATA = A' * A;计算AᵀA。[V_raw, D_V] = eig(ATA);对AᵀA进行特征值分解。V_raw包含特征向量（列），D_V是特征值对角矩阵。eigenvalues_ATA = diag(D_V);提取特征值。[sorted_eigenvalues, idx] = sort(eigenvalues_ATA, 'descend');将特征值降序排序，并记录排序索引。singular_values = sqrt(max(0, sorted_eigenvalues));计算奇异值。max(0, ...)用于避免因浮点误差导致的负数平方根问题。V = V_raw(:, idx);根据排序索引重新排列V_raw的列，得到正确的V矩阵。
-
构建
S矩阵：S = zeros(m, n);初始化一个m × n的零矩阵。min_dim = min(m, n);确定奇异值的数量，即min(m, n)。S(1:min_dim, 1:min_dim) = diag(singular_values(1:min_dim));将计算出的奇异值填充到S的主对角线上。
-
计算
U矩阵并校正符号：AAT = A * A';计算AAᵀ。[U_raw, D_U] = eig(AAT);对AAᵀ进行特征值分解，得到U_raw。[~, idx_U] = sort(diag(D_U), 'descend');排序AAᵀ的特征值（它们的平方根也是奇异值），并获取索引。U_full = U_raw(:, idx_U);根据排序重新排列U_raw的列，得到一个初步的U矩阵。- 符号校正循环： 由于特征向量的符号是任意的，
U_full的列可能与实际需要的左奇异向量方向相反。为了确保A = U Σ Vᵀ关系成立，我们遍历前rank_A个奇异值。expected_U_dir = (A * V(:, i)) / singular_values(i);根据Av = σu的关系计算第i个左奇异向量的期望方向。if dot(U_full(:, i), expected_U_dir) < 0：检查U_full中对应列与期望方向的点积。如果点积为负，说明它们方向相反。U_full(:, i) = -U_full(:, i);：翻转U_full对应列的符号，使其与expected_U_dir方向一致。
U = U_full;将校正后的U_full赋值给U。
使用 mySVD 函数的示例：
“`matlab
% 定义一个示例矩阵
A = [1 2; 3 4; 5 6]; % 3×2 矩阵
% 使用自定义的 mySVD 函数
[U_my, S_my, V_my] = mySVD(A);
disp(‘Custom SVD U:’);
disp(U_my);
disp(‘Custom SVD S:’);
disp(S_my);
disp(‘Custom SVD V:’);
disp(V_my);
% 验证分解结果
A_reconstructed_my = U_my * S_my * V_my’;
fprintf(‘原始矩阵 A 与 mySVD 重构矩阵的差的范数: %e\n’, norm(A – A_reconstructed_my));
% 与 MATLAB 内置 svd 结果进行比较
[U_builtin, S_builtin, V_builtin] = svd(A);
disp(‘Built-in SVD U:’);
disp(U_builtin);
disp(‘Built-in SVD S:’);
disp(S_builtin);
disp(‘Built-in SVD V:’);
disp(V_builtin);
% 注意：U 和 V 的列的符号可能有所不同，但其基是相同的，因为奇异值分解是唯一的，但 U 和 V 的列的符号可以反转。
% 只要 USV’ 结果一致即可。
“`
通过这个 mySVD 函数，我们可以直观地看到 SVD 是如何通过特征值分解一步步构建出来的。
四、 总结与应用
奇异值分解（SVD）是矩阵理论和应用中一个极其重要且强大的工具。它以一种优雅的方式将任意矩阵分解为三个具有特定结构和物理意义的矩阵乘积：左奇异向量矩阵 U、奇异值对角矩阵 Σ 和右奇异向量矩阵 V 的转置。
本文详细阐述了 SVD 的核心原理，揭示了其与特征值分解的深层联系，并通过 AᵀA 和 AAᵀ 的特征值和特征向量来推导奇异值和奇异向量。我们还介绍了 SVD 的基本计算算法，并特别关注了在 MATLAB 中的实现：
- MATLAB 内置
svd函数：这是在实际工程和科学计算中最常用和推荐的方法，它提供了高度优化和数值稳定的实现，支持多种调用形式，包括经济型 SVD，以适应不同需求。 - 自定义
mySVD函数：通过从零实现一个基于特征值分解的mySVD函数，我们深入理解了 SVD 的数学机制和计算步骤，包括了如何处理奇异值的排序和奇异向量的符号校正。
SVD 的应用范围极其广泛，包括但不限于：
- 数据降维：如主成分分析（PCA），通过保留最大的奇异值及其对应的奇异向量来捕捉数据的主要方差，从而减少数据的维度。
- 图像处理：用于图像压缩、去噪和特征提取。
- 推荐系统：通过对用户-物品评分矩阵进行 SVD，可以发现潜在的关联性，从而进行个性化推荐。
- 自然语言处理 (NLP)：如潜在语义分析 (LSA)，用于发现文本中词语和文档之间的隐藏语义关系。
- 机器学习：在各种算法中作为预处理步骤，提高模型性能和效率。
- 数值稳定性：解决最小二乘问题、计算矩阵的秩、伪逆等。
通过掌握 SVD 的原理和在 MATLAB 中的实现，读者将能够有效地运用这一强大工具解决复杂的科学计算和数据分析问题。理解其背后的数学逻辑，将有助于更好地选择和优化 SVD 在特定场景下的使用策略。