MATLAB奇异值分解(SVD)从入门到精通 – wiki基地

SVD (Singular Value Decomposition) is a powerful matrix factorization technique with widespread applications in various fields such as data analysis, image processing, and machine learning. This article provides a comprehensive guide to understanding and implementing SVD in MATLAB, covering concepts from basic introduction to advanced applications.

1. 奇异值分解 (SVD) 基础

1.1 什么是SVD？

SVD 是一种将任何 $m \times n$ 矩阵 $A$ 分解为三个矩阵乘积的形式：
$A = U \Sigma V^T$

其中：
* $U$ 是一个 $m \times m$ 的正交矩阵，其列向量是 $AA^T$ 的特征向量，称为左奇异向量。
* $\Sigma$ 是一个 $m \times n$ 的对角矩阵，其对角线上的元素是非负实数，称为奇异值。奇异值通常按降序排列。除了对角线元素外，所有元素都为零。
* $V$ 是一个 $n \times n$ 的正交矩阵，其列向量是 $A^TA$ 的特征向量，称为右奇异向量。$V^T$ 是 $V$ 的转置。

奇异值是矩阵 $A$ 的“重要性”的度量，较大的奇异值对应着更重要的信息。

1.2 为什么SVD很重要？

SVD 具有以下重要特性：
* 普适性： 任何实数或复数矩阵都可以进行SVD分解，无论其是否为方阵，也无论其是否可逆。
* 数据压缩： 可以通过保留最大的几个奇异值来近似原始矩阵，实现数据降维和压缩。
* 噪声去除： 较小的奇异值通常与噪声或不重要的信息相关联，去除这些奇异值有助于降噪。
* 基础分析： SVD 提供了矩阵的低秩近似，揭示了数据的主要结构和关系。

2. MATLAB 中的 SVD

MATLAB 提供了内置函数 svd 来执行奇异值分解。

2.1 基本用法

最简单的用法是直接对矩阵 $A$ 调用 svd 函数。

“`matlab
% 创建一个示例矩阵
A = [1 2 3; 4 5 6; 7 8 9; 10 11 12];

% 执行 SVD
[U, S, V] = svd(A);

% 显示结果
disp(‘矩阵 U:’);
disp(U);
disp(‘奇异值矩阵 S:’);
disp(S);
disp(‘矩阵 V:’);
disp(V);

% 验证分解: A 应该约等于 U * S * V’
A_reconstructed = U * S * V’;
disp(‘重构矩阵 A:’);
disp(A_reconstructed);
“`

输出解释：
* U 是一个 $4 \times 4$ 的正交矩阵。
* S 是一个 $4 \times 3$ 的对角矩阵。注意，MATLAB 默认会生成一个与 $A$ 大小相同的对角矩阵 S，其对角线上的元素是奇异值。如果 $A$ 是 $m \times n$，则 S 也是 $m \times n$。
* V 是一个 $3 \times 3$ 的正交矩阵。

2.2 经济型 SVD (Economy SVD)

当 $m > n$ 时 (行数多于列数)，矩阵 $U$ 将包含许多零向量。为了节省计算和存储资源，可以使用经济型 SVD。

“`matlab
% 经济型 SVD
[U_e, S_e, V_e] = svd(A, ‘econ’);

% 显示结果
disp(‘经济型 SVD – 矩阵 U_e:’);
disp(U_e); % U_e 是 4×3 矩阵
disp(‘经济型 SVD – 奇异值矩阵 S_e:’);
disp(S_e); % S_e 是 3×3 矩阵
disp(‘经济型 SVD – 矩阵 V_e:’);
disp(V_e); % V_e 是 3×3 矩阵

% 验证分解
A_reconstructed_econ = U_e * S_e * V_e’;
disp(‘经济型 SVD 重构矩阵 A:’);
disp(A_reconstructed_econ);
“`
在经济型 SVD 中，如果 $m > n$，则 $U$ 的前 $n$ 列和 $S$ 的前 $n$ 行 $n$ 列会被保留，而 $V$ 保持不变。如果 $n > m$，则 $U$ 保持不变，$S$ 的前 $m$ 行 $m$ 列和 $V$ 的前 $m$ 列会被保留。

2.3 仅获取奇异值

如果你只需要奇异值，而不需要 $U$ 和 $V$ 矩阵，可以直接调用 svd 返回一个向量：

matlab
s = svd(A);
disp('奇异值向量 s:');
disp(s);
s 将是一个列向量，包含按降序排列的奇异值。

3. SVD 的实际应用

SVD 在 MATLAB 中的应用非常广泛，以下是一些常见的例子。

3.1 图像压缩与重建

SVD 可以用于图像压缩，通过保留最重要的奇异值来近似图像。

“`matlab
% 读取图像 (确保你有 ‘peppers.png’ 文件在当前目录或路径中)
try
img = imread(‘peppers.png’);
catch
disp(‘peppers.png not found. Downloading a sample image…’);
img_url = ‘https://i.stack.imgur.com/nJ5QJ.png’; % A sample image URL
websave(‘sample_image.png’, img_url);
img = imread(‘sample_image.png’);
end

% 如果是彩色图像，转换为灰度图像 (SVD通常在灰度图像上操作)
if size(img, 3) == 3
img_gray = rgb2gray(img);
else
img_gray = img;
end

% 转换为双精度浮点类型进行计算
img_double = im2double(img_gray);

% 执行 SVD
[U, S, V] = svd(img_double);

% 选择不同数量的奇异值进行图像重建
num_singular_values = [10, 50, 100]; % 可以尝试不同值

figure;
subplot(2,2,1);
imshow(img_double);
title(‘原始图像’);

for i = 1:length(num_singular_values)
k = num_singular_values(i);

% 构建近似奇异值矩阵
Sk = S;
Sk(k+1:end, :) = 0; % 将 k 之后的奇异值设为 0
Sk(:, k+1:end) = 0;

% 重建图像
img_reconstructed = U * Sk * V';

subplot(2,2,i+1);
imshow(img_reconstructed);
title(sprintf('保留 %d 个奇异值', k));

end
“`
通过观察不同数量奇异值重建的图像，你可以看到随着保留奇异值数量的增加，重建图像的质量也随之提高。

3.2 矩阵的秩和近似

SVD 可以帮助确定矩阵的数值秩，并找到最佳的低秩近似。矩阵的秩等于非零奇异值的数量。

“`matlab
% 创建一个秩亏缺的矩阵
A_rank_deficient = [1 2 3; 2 4 6; 3 6 9; 1 1 1];

% 计算奇异值
s_rd = svd(A_rank_deficient);
disp(‘秩亏缺矩阵的奇异值:’);
disp(s_rd);

% 设定一个小的容差来判断非零奇异值 (由于浮点数误差)
tolerance = max(size(A_rank_deficient)) * eps(max(s_rd));
numerical_rank = sum(s_rd > tolerance);
fprintf(‘矩阵的数值秩为: %d\n’, numerical_rank);

% 低秩近似 (例如，秩为 1 的近似)
k_approx = 1;
[U_rd, S_rd, V_rd] = svd(A_rank_deficient);

% 保留 k_approx 个奇异值进行近似
S_approx = zeros(size(S_rd));
S_approx(1:k_approx, 1:k_approx) = S_rd(1:k_approx, 1:k_approx);
A_approx = U_rd * S_approx * V_rd’;
disp(‘秩为 1 的近似矩阵:’);
disp(A_approx);
“`
SVD 允许我们系统地找到给定矩阵的最佳低秩近似（在 Frobenius 范数意义下）。

3.3 求解线性方程组

当线性方程组 $Ax = b$ 中的矩阵 $A$ 是奇异的、病态的或非方阵时，SVD 可以提供一个稳定的最小二乘解。

“`matlab
% 创建一个欠定系统 (m < n)
A_under = [1 1 1; 2 1 3];
b_under = [6; 10];

% 使用 SVD 求解最小二乘解
[U_under, S_under, V_under] = svd(A_under);

% 奇异值矩阵的伪逆
S_plus = zeros(size(S_under’)); % S_plus 的大小是 S 的转置
non_zero_singular_values = diag(S_under);
for i = 1:length(non_zero_singular_values)
if non_zero_singular_values(i) > eps * max(size(A_under)) * max(non_zero_singular_values)
S_plus(i, i) = 1 / non_zero_singular_values(i);
end
end

x_svd = V_under * S_plus * U_under’ * b_under;
disp(‘SVD 得到的最小二乘解 x:’);
disp(x_svd);

% 验证 A * x_svd 应该接近 b
disp(‘A * x_svd:’);
disp(A_under * x_svd);
“`

3.4 主成分分析 (PCA)

SVD 是执行主成分分析 (PCA) 的基础。PCA 用于数据降维和特征提取。

“`matlab
% 生成一些示例数据
data = [1 2; 2 3; 3 4; 4 5; 5 6] + randn(5,2)*0.1;

% 数据中心化 (减去均值)
mean_data = mean(data);
centered_data = data – mean_data;

% 计算协方差矩阵 (可选，SVD可以直接作用于中心化数据)
% covariance_matrix = cov(centered_data);

% 对中心化数据执行 SVD
[U_pca, S_pca, V_pca] = svd(centered_data);

% V_pca 的列就是主成分 (特征向量)
principal_components = V_pca;
disp(‘主成分 (V 矩阵的列):’);
disp(principal_components);

% 投影数据到主成分空间 (降维)
% 假设我们想降到 1 维 (保留第一个主成分)
k_pca = 1;
reduced_data = centered_data * principal_components(:, 1:k_pca);
disp(sprintf(‘降维后的数据 (保留 %d 个主成分):’, k_pca));
disp(reduced_data);

% 可视化
figure;
scatter(data(:,1), data(:,2), ‘b’, ‘filled’);
hold on;
quiver(mean_data(1), mean_data(2), principal_components(1,1), principal_components(2,1), ‘r’, ‘LineWidth’, 2);
quiver(mean_data(1), mean_data(2), principal_components(1,2), principal_components(2,2), ‘g’, ‘LineWidth’, 2);
title(‘原始数据和主成分’);
legend(‘数据点’, ‘第一主成分’, ‘第二主成分’);
grid on;
hold off;
``
在这个例子中，V_pca矩阵的列即为主成分的方向。我们可以通过将中心化数据乘以V_pca` 的前 $k$ 列来实现降维。

4. 深入理解和高级主题

4.1 奇异值与特征值

奇异值是与矩阵 $A$ 和其转置 $A^T$ 相关的两个对称矩阵的特征值的平方根：
* $A^TA$ 的特征值的平方根是奇异值。
* $AA^T$ 的特征值的平方根也是奇异值。

左奇异向量 $U$ 是 $AA^T$ 的特征向量，右奇异向量 $V$ 是 $A^TA$ 的特征向量。

4.2 几何解释

从几何角度看，SVD 可以解释为：任何线性变换都可以分解为一次旋转 (由 $V^T$ 表示)、一次缩放 (由 $\Sigma$ 表示) 和另一次旋转 (由 $U$ 表示)。对于一个单位球体，SVD 描述了该球体在经过线性变换后如何变为一个椭球体，其中奇异值是椭球体的半轴长度，奇异向量是半轴的方向。

4.3 截断 SVD (Truncated SVD)

截断 SVD 指的是只保留前 $k$ 个最大的奇异值来近似原始矩阵。这正是图像压缩和 PCA 降维所采用的核心思想。它提供了一种最优的低秩近似，这意味着在给定秩 $k$ 的情况下，通过截断 SVD 得到的近似矩阵与原始矩阵之间的差异（Frobenius 范数）是最小的。

4.4 稀疏 SVD

对于非常大的稀疏矩阵，计算完整的 SVD 效率低下。MATLAB 提供了 svds 函数来计算稀疏矩阵的部分奇异值和奇异向量。

“`matlab
% 创建一个稀疏矩阵
A_sparse = sprand(1000, 1000, 0.01); % 1000×1000, 1% 非零元素

% 计算前 k 个最大的奇异值和奇异向量
k_sparse = 5;
[U_s, S_s, V_s] = svds(A_sparse, k_sparse);

disp(sprintf(‘稀疏矩阵的前 %d 个奇异值:’, k_sparse));
disp(diag(S_s));
``svds` 在处理大规模稀疏数据时非常有用，例如在推荐系统和自然语言处理中。

5. 总结

SVD 是线性代数中一个极其重要的工具，它的应用遍布科学和工程的各个领域。在 MATLAB 中，svd 和 svds 函数使得执行奇异值分解变得简单而高效。从基本的矩阵分解到复杂的图像处理、数据降维和数值分析，掌握 SVD 及其在 MATLAB 中的实现，将极大地增强你的数据处理和分析能力。