奇异值分解的应用及MATLAB实现（图像压缩、谱聚类、潜在语义分析、主成分分析、背景删除等）

相关文章

奇异值分解介绍
从几何视角看奇异值分解（MATLAB实现）

奇异值分解的应用

图像压缩

数字图像通常可以被表示为一个矩阵，其中每个元素代表图像像素的亮度或颜色值。通过奇异值分解，我们可以将图像矩阵分解为多个矩阵的乘积，这些矩阵包含了图像的主要特征信息，而图像的大部分信息通常集中在少量的奇异值上，较小的奇异值往往包含的是噪声或细微的细节。因此，保留较大的奇异值可以有效地保持图像的主要特征。 图像压缩步骤如下：

• 图像矩阵化： 首先，将图像转换为矩阵形式。例如，一幅灰度图像可以表示为一个 $m\times n$ 的矩阵，每个元素表示像素的灰度值。
• 应用奇异值分解： 对图像矩阵进行奇异值分解，得到三个矩阵 $\mathbf{U}$、$\mathbf{\Sigma }$和 ${\mathbf{V}}^{\text{T}}$。
• 选择前 k 个奇异值： 通常，大部分图像信息都集中在较大的奇异值上。通过选择前 k 个最大的奇异值及其对应的奇异向量，可以近似原图像矩阵。

${\mathbf{A}}_{\mathbf{k}}={\mathbf{U}}_{\mathbf{k}}{\mathbf{\Sigma }}_{\mathbf{k}}{\mathbf{V}}_{\mathbf{k}}^{\text{T}}$
其中 ${\mathbf{U}}_{\mathbf{k}}$ 是前 k 个左奇异向量构成的矩阵，${\mathbf{\Sigma }}_{\mathbf{k}}$ 是前 k个奇异值构成的对角矩阵，${\mathbf{V}}_{\mathbf{k}}$ 是前 k 个右奇异向量构成的矩阵。

• 图像重构： 使用选定的 k个奇异值及其对应的奇异向量重构图像矩阵，尽管这种重构的图像会丢失一些细节信息，但其主要特征仍然保留。
• 压缩效果： 通过减少保留的奇异值数量，可以显著降低数据存储需求。例如，如果原矩阵的大小是 $m\times n$，那么重构矩阵只需要存储 ${\mathbf{U}}_{\mathbf{k}}$ （ $m\times k$）、${\mathbf{\Sigma }}_{\mathbf{k}}$ （大小为 $k\times k$）和${\mathbf{V}}_{\mathbf{k}}$ （大小为 $k\times n$），这大大减少了数据量。

下图展示利用奇异值分解进行图像压缩的结果，分别取前$k=150,100,80,50,20$个最大的奇异值及其对应的奇异向量压缩后得到的图像。当$k=150,100$时，压缩后的图像对原图像的还原度较高，视觉差异较小；当$k=80$时，以及可以感受到轻微失真，随着$k$值的减少这种失真越来越明显。

图像压缩MATLAB实现：

% 读取图像
img = imread('pic/lenna_gray.jpg');

% 转换为灰度图像
% gray_img = rgb2gray(img);

% 将图像数据转换为双精度类型
A = double(img);

% 进行 SVD 分解
[U, S, V] = svd(A);

% 选择奇异值的数量
k = 200;

% 使用前 k 个奇异值重构图像
A_compressed = U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k)';

% 显示原始图像
subplot(1, 2, 1);
imshow(uint8(A));
title('原始图像');

% 显示压缩后的图像
subplot(1, 2, 2);
imshow(uint8(A_compressed));
title(['压缩后的图像，k = ' num2str(k)]);

% 保存压缩后的图像
imwrite(uint8(A_compressed), 'compressed_image.png');

% 计算压缩率
original_size = numel(U) + numel(S) + numel(V);
compressed_size = k * (size(U, 1) + size(V, 2) + 1);
compression_ratio = compressed_size / original_size;
fprintf('压缩率：%.2f%%\n', compression_ratio * 100);

谱聚类

谱聚类是一种基于图论的聚类方法，主要通过研究样本数据的相似性矩阵来进行聚类，它能够识别复杂的数据结构，因此可以用于数据聚类，其步骤如下：

• 构建相似度矩阵： 根据数据点之间的距离构建一个相似度矩阵。常见的方法是高斯核函数（RBF kernel），定义如下：

$S_{ij}=\exp \left(-\frac{\Vert x_i-x_j{\Vert }^2}{2{\sigma }^2}\right)$
其中，$x_i$ 和 $x_j$ 是样本点，$\sigma$是核函数的带宽参数。

• 构建图的拉普拉斯矩阵：从相似度矩阵出发，构建拉普拉斯矩阵，标准的无向图拉普拉斯矩阵定义为：

$L=D-S$
其中，$D$ 是对角度矩阵，其对角线元素是相似度矩阵中每行元素的和。

• 计算拉普拉斯矩阵的特征值和特征向量：计算拉普拉斯矩阵的特征值和对应的特征向量（根据构造拉普拉斯矩阵的过程可知拉普拉斯矩阵可以对角化）。选择最小的$k$个特征值对应的特征向量。
• 使用特征向量进行聚类：将特征向量作为新的特征空间，使用传统的聚类方法（K-means）在这个空间中进行聚类。

下图展示了$k$分别取$6、31$的谱聚类的结果，从图中可以看到两种聚类方法都取得了较好的分类结果：

谱聚类MATLAB实现：

% 生成示例数据
load('data/galaxy.mat','data'); % 假设数据存储在名为 'galaxy.mat' 的变量中

% 假设 data 是已加载的两列数据
% n = size(data, 1);

% 计算高斯相似性矩阵
sigma = 2;  % 高斯核的带宽参数
S = exp(-pdist2(data, data).^2 / (2*sigma^2));

% 度矩阵 D 和拉普拉斯矩阵 L
D = diag(sum(S, 2));
L = D - S;

% 计算特征值和特征向量
[V, D] = eig(L);

% 基于第二小的非零特征值的特征向量进行 k-means 聚类
k = 31;  % 假设我们预计有31/6个聚类
second_smallest_eigenvector = V(:,2:k+1);
clusters = kmeans(second_smallest_eigenvector, k);
% clusters = kmeans(data, k);

% 可视化聚类结果
gscatter(data(:,1), data(:,2), clusters);
% scatter(data(:,1), data(:,2), 50, clusters, 'filled');
title('谱聚类结果');

潜在语义索引

潜在语义索引（Latent Semantic Indexing, LSI）是一种通过对原始词-文档矩阵进行降阶近似来改进词-文档搜索和检索的方法。以下是LSI的主要内容和步骤：

• 词-文档矩阵构建： 构建一个$M\times N$的矩阵$L$，其中$[i,j]$表示词$i$在文档$j$中的出现情况，可以使用0和1表示词的出现，也可以使用权重表示词的重要性或频率。
• 查询与文档的相似度计算： 使用余弦相似度来计算查询向量$\mathbf{q}$和文档的列向量之间的相似度，公式为：$\cos ({\theta }_j)=\frac{(L{e}_j{)}^T\mathbf{q}}{||L{e}_j||\cdot ||\mathbf{q}||}$

其中，$e_j$是单位矩阵的第$j$列。

• 奇异值分解 ： 对词-文档矩阵$L$进行奇异值分解：

$L=U\Sigma V^T$
其中，矩阵$U$的列表示“超级文档”，这些超级文档通过结合实际文档捕捉到大多数术语使用的变化。选择一个秩$k$，使得$L$的秩$k$近似能够捕捉大部分变化。通常选择$k$满足${\sigma }_1^2+\dots +{\sigma }_k^2$与总变化之比大约为0.9。

• 截断矩阵： 计算截断矩阵$L_k=U_k{\Sigma }_k{V}_k^T$，其中，$U_k{\Sigma }_k{V}_k^T$表示前$k$个超级文档的近似。

使用截断矩阵计算新的余弦相似度：$\cos (\widehat{{\theta }_j})=\frac{(L_k{e}_j{)}^T\mathbf{q}}{||L_k{e}_j||\cdot ||\mathbf{q}||}$

• 投影与优化： 对查询向量$\mathbf{q}$进行投影：$\hat{\mathbf{q}}=U_k^T\mathbf{q}$投影结果表示为$U_k\hat{\mathbf{q}}$。

最终得到优化后的余弦相似度：$\cos (\widehat{{\theta }_j})\le \frac{s_j^T\hat{\mathbf{q}}}{||s_j||\cdot ||\hat{\mathbf{q}}||}=\cos ({\theta }_j^{\ast })$
其中，$s_j={\Sigma }_k{V}_k^{\text{T}}{e}_j$是独立于查询的预计算向量。
通过这些步骤，LSI 能够识别与查询相关的更多文档，且通常能够比之前的分析返回更多的文档，从而增强搜索结果。下面以一个具体的例子进行展示：
假设我们有一个$5\times 4$的矩$L$阵表示5个词在4个文档中的出现情况：
$L=\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 1& 1\\ 1& 0& 0& 1\\ 1& 0& 1& 1\end{array}\right]$
假设查询向量$\mathbf{q}$为$\left[\begin{array}{c}1;0;1;0;0\end{array}\right]$表示查询包含词1和词3；
对词-文档矩阵$L$进行分解进行奇异值分解：
$$\begin{gathered} L=US{V}^{\mathrm{T}}U=\left[\begin{array}{cccc}-0.43& 0.50& 0.50& -0.26\\ -0.43& 0.50& -0.50& -0.26\\ -0.43& -0.50& -0.50& -0.26\\ -0.43& -0.50& 0.50& -0.26\\ -0.53& 0& 0& 0.85\end{array}\right] \\ S=\left[\begin{array}{cccc}2.29& 0& 0& 0\\ 0& 1.41& 0& 0\\ 0& 0& 1.41& 0\\ 0& 0& 0& 0.8740\end{array}\right] \\ {V}^{\mathrm{T}}=\left[\begin{array}{cccc}-0.60& 0& 0.71& 0.37\\ -0.37& 0.71& 0& -0.60\\ -0.60& 0& -0.71& 0.37\\ -0.37& -0.71& 0& -0.60\end{array}\right] \end{gathered}$$
选取前$k=2$个奇异值及其对应的奇异值向量：
$U_k=\left[\begin{array}{cc}-0.43& 0.50\\ -0.43& 0.50\\ -0.43& -0.50\\ -0.43& -0.50\\ -0.53& 0\end{array}\right]$
$S_k=\left[\begin{array}{cc}2.29& 0\\ 0& 1.41\end{array}\right]$
$V_k^{\mathrm{T}}=\left[\begin{array}{cc}-0.60& 0\\ -0.37& 0.71\\ -0.60& 0\\ -0.37& -0.71\end{array}\right]$
将查询向量$q$投影到 $U_k$ 的列空间：
$\hat{\mathbf{q}}=U_k^T\mathbf{q}=\left[\begin{array}{c}0.85\\ 0\end{array}\right]$
计算文档的余弦相似度：
$s_1={\Sigma }_k{V}_k^{\text{T}}{e}_1=\left[\begin{array}{c}-1.3764\\ 0\end{array}\right]$
$\cos ({\theta }_1^{\ast })=\frac{s_1^T\hat{\mathbf{q}}}{||s_1||\cdot ||\hat{\mathbf{q}}||}=1.00$

同理，计算其它文档余弦相似度：
$\cos ({\theta }_2^{\ast })=0.65,\cos ({\theta }_3^{\ast })=1.00,\cos ({\theta }_4^{\ast })=0.65$
最终，最相关的文档顺序为：1、3、2、4

潜在语义索引MATLAB实现：

% 1. 创建词-文档矩阵 L
% 假设我们有一个 5x4 的矩阵表示 5 个词在 4 个文档中的出现情况
L = [
    1 1 0 0; % 词1
    0 1 1 0; % 词2
    0 0 1 1; % 词3
    1 0 0 1; % 词4
    1 0 1 0  % 词5
];

% 2. 对矩阵进行奇异值分解 (SVD)
[U, S, V] = svd(L, 'econ');

% 3. 选择前 k 个奇异值及其对应的奇异向量
k = 2; % 可以根据需要调整
Uk = U(:, 1:k);
Sk = S(1:k, 1:k);
Vk = V(:, 1:k);

% 4. 计算查询向量 q 的投影
% 假设查询向量 q 为 [1; 0; 1; 0; 0] 表示查询包含词1和词3
q = [1; 0; 1; 0; 0];
q_proj = Uk' * q; % 将 q 投影到 Uk 的列空间

% 5. 计算文档向量的投影和余弦相似度
numDocs = size(Vk, 1);
similarities = zeros(numDocs, 1);
for j = 1:numDocs
    % 计算文档 j 的投影
    doc_proj = Sk * Vk(j, :)';
    % 计算余弦相似度
    similarities(j) = dot(q_proj, doc_proj) / (norm(q_proj) * norm(doc_proj));
end

% 6. 输出相似度结果
disp('文档相似度：');
disp(similarities);

% 7. 按相似度排序，得到最相关的文档
[~, sorted_indices] = sort(similarities, 'descend');
disp('最相关的文档顺序：');
disp(sorted_indices);

主成分分析

主成分分析 (Principal Component Analysis, PCA) 是重要的降维工具。PCA 可以显著减少数据 的维数，同时保留数据中对方差贡献最大的成分。简单来说，PCA 的核心思想是通过线性变换将 高维数据映射到低维空间中，使得映射后的数据能够尽可能地保留原始数据的信息，同时去除噪 声和冗余信息，从而更好地描述数据的本质特征。以下展示了如何利用奇异值分解实现主成分分析：
首先中心化数据矩阵，对于已经去均值的矩阵${\mathbf{X}}_{\mathbf{n}\times \mathbf{D}}$进行奇异值分解，得到：
$\mathbf{X}=\mathbf{US}{\mathbf{V}}^{\mathrm{T}}(3.1)$
由于$\mathbf{U}$和$\mathbf{V}$均为正交矩阵，即满足：
$\begin{array}{rl}{\mathbf{UU}}^{\mathrm{T}}& ={\mathbf{U}}^{\mathrm{T}}\mathbf{U}=\mathbf{I}\\ {\mathbf{VV}}^{\mathrm{T}}& ={\mathbf{V}}^{\mathrm{T}}\mathbf{V}=\mathbf{I}\end{array}$
由于$\mathbf{X}$已经中心化，其协方差矩阵可以通过以下公式计算得到：
$\mathbf{\Sigma }=\frac{{\mathbf{X}}^{\mathrm{T}}\mathbf{X}}{\mathbf{n}-1}(3.2)$
将$(3.1)$代入$(3.2)$，可得：
$\mathbf{\Sigma }=\frac{{\left({\mathbf{USV}}^{\mathrm{T}}\right)}^{\mathrm{T}}{\mathbf{USV}}^{\mathrm{T}}}{n-1}=\frac{{\mathbf{VS}}^{\mathrm{T}}{\mathbf{SV}}^{\mathrm{T}}}{n-1}(3.3)$
对协方差矩阵进行特征值分解，结合公式$(1.2)$和协方差矩阵的求解过程，可知：
$\mathbf{\Sigma }=\mathbf{V\Lambda }{\mathbf{V}}^{\mathrm{T}}(3.4)$
联立$(3.3)$和$(3.4)$，可得：
$\frac{{\mathbf{VS}}^{\mathrm{T}}{\mathbf{SV}}^{\mathrm{T}}}{n-1}=\mathbf{V\Lambda }{\mathbf{V}}^{\mathrm{T}}(3.5)$
对于不含零矩阵的对角矩阵$\mathbf{S}$，整理$(3.5)$可得：
$\frac{{\mathbf{S}}^2}{\mathbf{n}-1}=\mathbf{\Lambda }$
即：
$\frac{1}{n-1}\left[\begin{array}{ccccc}{s}_1^2& & & & \\ & s_2^2& & & \\ & & \ddots & & \\ & & & s_D^2& \end{array}\right]=\left[\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_D\end{array}\right]$
其中，${\lambda }_1⩾{\lambda }_2⩾...⩾{\lambda }_D$，而奇异值和特征值存在如下关系：
$\frac{s_j^2}{n-1}={\lambda }_j$
$s_j$为第$j$个主成分的奇异值，${\lambda }_j$为协方差矩阵的第$j$个特征值。

背景删除

视频背景和前景的分离在计算机视觉和视频处理领域是一个重要的任务。背景矩阵通常是低秩的，因为背景在时间上变化很小；而前景矩阵包含了视频中所有的变化部分。因此，我们可以通过奇异值分解来分解原始矩阵，并从中提取出背景和前景，具体步骤如下：

• 创建矩阵$\mathbf{M}$：首先，从视频中定期采样图像帧，并将它们展平为一维数组，这些数组构成了矩阵M的列，如图3-3所示，对矩阵$\mathbf{M}$进行可视化，可以发现图上大部分是由水平的线组成，而这些水平的线条是由于背景重复造成的；

• 计算矩阵$\mathbf{M}$的奇异值分解： 使用奇异值分解对矩阵$\mathbf{M}$进行分解；
• 获取低秩近似： 选择矩阵的前k个奇异值来构建一个低秩近似矩阵，它表示背景；
• 计算前景矩阵 ： 通过从原始矩阵$\mathbf{M}$中减去背景矩阵，得到前景矩阵（图3-4所示）。

背景删除MATLAB代码实现：

% 指定图像文件夹路径
imageFolder = 'video'; % 替换为你的图像文件夹路径
imageFiles = dir(fullfile(imageFolder, '*.jpg')); % 读取所有.jpg文件
numImages = length(imageFiles);

% 初始化变量
sampleImage = imread(fullfile(imageFolder, imageFiles(1).name));
[height, width, numChannels] = size(sampleImage);

% 检查是否为灰度图像
if numChannels == 3
    sampleImage = rgb2gray(sampleImage);
end

% 初始化矩阵M
M = zeros(height * width, numImages);

% 填充矩阵M
for i = 1:numImages
    imageFile = fullfile(imageFolder, imageFiles(i).name);
    image = imread(imageFile);
    
    if size(image, 3) == 3 % 如果是彩色图像，转换为灰度图像
        image = rgb2gray(image);
    end
    
    M(:, i) = image(:); % 展平成列向量并存入矩阵M
end

% 使用imagesc可视化矩阵M
figure;
imagesc(M);
colormap('gray');
colorbar;
title('Matrix M Visualization');
xlabel('Image Index');
ylabel('Pixel Index');

% 进行SVD分解
[U, S, V] = svd(M, 'econ');

% 选择前k个奇异值来获得背景矩阵
k = 1; % 可以根据需要调整k值
U_k = U(:, 1:k);
S_k = S(1:k, 1:k);
V_k = V(:, 1:k);

% 构建低秩背景矩阵
M_background = U_k * S_k * V_k';

% 计算前景矩阵
M_foreground = M - M_background;

% 创建输出文件夹
outputFolder = 'output_images';
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

% 保存背景和前景图像
for i = 1:numImages
    backgroundFrame = reshape(M_background(:, i), [height, width]);
    foregroundFrame = reshape(M_foreground(:, i), [height, width]);

    % 将背景和前景保存为图像文件
    imwrite(uint8(backgroundFrame), fullfile(outputFolder, sprintf('background_%03d.png', i)));
    imwrite(uint8(foregroundFrame), fullfile(outputFolder, sprintf('foreground_%03d.png', i)));
end

disp('背景和前景图像已保存至文件夹: output_images');

参考文献

[1] Feeman, Timothy G. 2023. Applied Linear Algebra and Matrix Methods. Springer Undergraduate Texts in Mathematics and Technology. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-031-39562-8.

[2] Visualize-ML. 2024. Visualize-ML/Book4_Power-of-Matrix. Retrieved from [GitHub Repository https://github.com/Visualize-ML/Book4_Power-of-Matrix].