SLAM常见问题(五):Singular Value Decomposition(SVD)分解

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

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

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

计算 $A^TA$

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

计算 $AA^T$

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

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

降维

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

实战

数字例子

有矩阵A,对其进行SVD分解,已知:

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

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

1
2
3
4
5
6
7
8
9
10
11
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%)。由此,我们实现了图像压缩。

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

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
%% 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])

参考