📝笔记：SLAM常见问题(五)：Singular Value Decomposition（SVD）分解

SVD分解就是一种矩阵拆解术，它能够把任意矩阵$A \in \mathbb{R}^{m \times n}$拆解成3个矩阵的乘积形式，即：

$$A = U \Sigma V^T$$

其中，$U \in \mathbb{R}^{m \times m}$，$V \in \mathbb{R}^{n \times n}$都是正交矩阵，即列向量是正交的单位向量，$\Sigma \in \mathbb{R}^{m \times n}$的对角阵（奇异值）。搬运了来自MIT OpenCourseWare的在线课程并放在了B站，讲解得很清晰。

刚才说了矩阵$U, \Sigma, V$的形式，视频中还提到了这三个矩阵的物理意义，即SVD分解可以理解为：任意矩阵都可以分解为(rotation)*(Stretch)*(rotation)的形式。接下来说明一下这三个矩阵是如何来的。

计算 $A^TA$

$$A^TA = (U \Sigma V^T)^TU \Sigma V^T = V{\Sigma}^TU^TU \Sigma V^T = V{\Sigma}^T \Sigma V^T$$

可见，$V$正是矩阵$A^TA$的特征向量，而${\Sigma}^T \Sigma $为矩阵$A^TA$的特征值。

计算 $AA^T$

$$AA^T = U \Sigma V^T(U \Sigma V^T)^T = U \Sigma V^TV{\Sigma}^TU^T = U\Sigma {\Sigma}^T U^T$$

可见，$U$正是矩阵$AA^T$的特征向量，而$\Sigma {\Sigma}^T $为矩阵$AA^T$的特征值。

所以$U, \Sigma, V$都可以通过上述方式来计算。

降维

$$\Sigma = \left[ \begin{array} {cccc|cccc} {\sigma_{1}} & {0} & {\dots} & {0} & {0} & {0} & {\dots} & {0} \\ {0} & {\sigma_{2}} & {\dots} & {0} & {0} & {0} & {\dots} & {0} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} & {0} & {0} & {\dots} & {0} \\ {0} & {0} & {} & {\sigma_{k}} & {0} & {0} & {\dots} & {0} \\ \hline {0} & {0} & {\dots} & {0} & {0} & {0} & {\dots} & {0} \\ {\vdots} & {\vdots} & {} & {\vdots} & {\vdots} & {\vdots} & {} & {\vdots} \\ {0} & {0} & {\dots} & {0} & {0} & {0} & {\dots} & {0} \end{array} \right]_{m \times n} \Rightarrow \left[ \begin{array} {cccc} {\sigma_{1}} & {0} & {\dots} & {0} \\ {0} & {\sigma_{2}} & {\dots} & {0} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {0} & {0} & {} & {\sigma_{k}} \end{array} \right]_{k \times k}$$

其中${\sigma}_1 \geq {\sigma}_2 \geq … {\sigma}_k > 0 $，将$\Sigma$中主对角线为0的部分删去，同样的$U,V$对应的部分删去，SVD分解就变成了下图的形式。

实战

数字例子

有矩阵A，对其进行SVD分解，已知：

$$A = \left[ \begin{matrix} ​ 1 & 4 & 3 & 5 & 6 \cr ​ 2 & 3 & {4} & 5 & 0 \cr ​ 7 & 4 & 0 & 9 & 1 \cr \end{matrix} \right]$$

计算$A^TA$以及$AA^T$：

$$A^TA = \left[ \begin{matrix} ​ {54} & {38} & {11} & {78} & {13} \cr ​ {38} & {41} & {24} & {71} & {28} \cr ​ {11} & {24} & {25} & {35} & {18} \cr ​ {78} & {71} & {35} & {131} & {39} \cr ​ {13} & {28} & {18} & {39} & {37} \end{matrix} \right] \\ AA^T = \left[ \begin{matrix} ​ {87} & {51} & {74} \cr ​ {51} & {54} & {71} \cr ​ {74} & {71} & {147} \cr \end{matrix} \right]$$

对以上两式做特征值分解得到：

V =
    0.4269    0.5222    0.1760   -0.5292   -0.4839
    0.4087   -0.1757   -0.0655   -0.5258    0.7221
    0.2100   -0.4474   -0.7536   -0.1512   -0.4062
    0.7389    0.1520   -0.0603    0.6481    0.0853
    0.2464   -0.6879    0.6271   -0.0258   -0.2687

U =
    0.5095    0.7999    0.3171
    0.4285    0.0838   -0.8997
    0.7462   -0.5942    0.3001

奇异值$\Sigma ^T \Sigma = \text{Diag}(238.2878, 37.3715, 12.3407) \Rightarrow \Sigma = \text{Diag}(15.4366, 6.1132, 3.5129)$

这与直接调用svd(A)结果是一致的（可能差个正负号）。

图像处理

祭上亲爱的Battle Angel Alita。

原始图像尺寸$1440\times 2560 $，我们可以对该图像做SVD分解，然后仅保留奇异值的前10，50，100重构图像，比较重构图像与原始图像的质量差异。可见仅仅保留其前10个奇异值时，图像质量遭到了极大破坏（此时仅保留原始图像信息的58.864%），随着奇异值数量的增多，图像质量也会逐渐提升，可以看到当奇异值个数为100时，基本上已经看不出与原图的差异（此时仅保留原始图像信息的87.37%）。由此，我们实现了图像压缩。

下图是保留的奇异值数量与图像质量的关系图，保留的奇异值越多，图像质量越高，图像压缩效果越不明显；反之，奇异值越少，图像质量越差，图像压缩效果越明显。这只是一种非常简单的图像压缩算法，仅作原理验证使用，在实际中用到的概率不是很大。

代码

%% simple test of using SVD decomposistion
clear all;
close all;
clc;

% A = [1 4 3 5 6;
%          2 3 4 5 0;
%          7 4 0 9 1]
% A'*A;
% A*A';
% 
% [V,Dv] = eig(A'*A);
% 
% lambda = wrev(diag(Dv));
% V = fliplr(V)
% 
% [U,Du] = eig(A*A');
% 
% lambda = wrev(diag(Du));
% U = fliplr(U)

%% LOADING IMAGE
img = imread('alita_origin.png');



ENERGE = 0;
for i = 1:3
    [U(:,:,i) D(:,:,i) V(:,:,i)] = svd(double(img(:,:,i)))  ;
    ENERGE = ENERGE +sum(diag(D(:,:,i)));
end

%% 10
DIM = 10;
ENERGE10 = 0;
for i = 1:3
    img_recons10(:,:,i) = U(:,1:DIM,i)*D(1:DIM,1:DIM,i)*V(:,1:DIM,i)';
     ENERGE10 = ENERGE10 +sum(diag(D(1:DIM,1:DIM,i)));
end
% figure;
% imshow(mat2gray(img_recons10))
% imwrite(mat2gray(img_recons10),'alita_10.png');

%% 50
DIM = 50;
ENERGE50 = 0;
for i = 1:3
    img_recons50(:,:,i) = U(:,1:DIM,i)*D(1:DIM,1:DIM,i)*V(:,1:DIM,i)';
     ENERGE50 = ENERGE50 +sum(diag(D(1:DIM,1:DIM,i)));
end
% figure;
% imshow(mat2gray(img_recons50))
% imwrite(mat2gray(img_recons50),'alita_50.png');

%% 100
DIM = 100;
ENERGE100 = 0;
for i = 1:3
    img_recons100(:,:,i) = U(:,1:DIM,i)*D(1:DIM,1:DIM,i)*V(:,1:DIM,i)';
    ENERGE100 = ENERGE100 +sum(diag(D(1:DIM,1:DIM,i)));
end
% figure;
% imshow(mat2gray(img_recons100))
% imwrite(mat2gray(img_recons100),'alita_100.png');

figure;
set(gcf,'pos',[ 986 414 1274 826])

FONTSIZE = 15;
h(1) = subplot(221);imshow(mat2gray(img)); 
xlabel('origin Alita');set(gca,'fontsize',FONTSIZE)

h(2) = subplot(222);imshow(mat2gray(img_recons10));
xlabel(['Using 10 singular values: ' num2str(ENERGE10/ENERGE)]);set(gca,'fontsize',FONTSIZE)

h(3) = subplot(223);imshow(mat2gray(img_recons50));
xlabel(['Using 50 singular values: ' num2str(ENERGE50/ENERGE)]);set(gca,'fontsize',FONTSIZE)

h(4) = subplot(224);imshow(mat2gray(img_recons100));
xlabel(['Using 100 singular values: ' num2str(ENERGE100/ENERGE)]);set(gca,'fontsize',FONTSIZE)
set(gcf,'color',[1 1 1])



%% SHOW ENERGY
ENERGY_tmp = zeros(size(img,1),1);
for DIM_ = 1:size(img,1)
   for i = 1:3
     ENERGY_tmp(DIM_,1) = ENERGY_tmp(DIM_,1) +sum(diag(D(1:DIM_,1:DIM_,i)));
   end
end

figure;
FONTSIZE = 30;
ratio = ENERGY_tmp/ENERGE;
X =   1:size(img,1);
plot(X,ratio,'linewidth',5,'color','r');
set(gcf,'color',[1 1 1])
xlabel('Number of Singular values');
ylabel('Image Quality');
set(gca,'fontsize',FONTSIZE)
set(gcf,'pos',[ 986 414 1274 826])

参考

• A Singularly Valuable Decomposition The SVD of a Matrix
• 李宏毅关于SVD的介绍，PPT,课程列表