0%

📝笔记: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}$都是正交矩阵,即$U^T U = I, V^T V = I$,即列向量是正交的单位向量,$U$称为left single vectors,$V$称为right single vectors;$\Sigma \in \mathbb{R}^{m \times n}$的对角阵(奇异值),奇异值$\geq 0$且按照顺序降序排列。

MIT Gilbert Strang 教授对 SVD 讲解得很清晰,如下:

此外,Barry讲解的SVD也很精彩。

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

形式: skinny SVD 与 full SVD

一张图可以解释了,上面介绍的SVD形式为skinny SVD 或者 thin SVD,即仅保留了奇异值矩阵 $\Sigma$ 的非零元素所在的部分,所以它的形状是个方阵。而full SVD保留了奇异值矩阵的非零部分,对应的 $U$ 和 $V$ 为方阵。

通常情况下,skinny SVD使用较多。

此外,还有向量表示形式,如上图最后一行的表示形式,此处不做赘述。

计算 $A^TA$

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

计算 $AA^T$

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

所以$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)结果是一致的(可能差个正负号)。

最小二乘 least square问题

MIT Gilbert Strang 教授介绍了4种方式来求解最小二乘问题。

问题:$ A x = b $,求解 $x$。

通常的解法如下:

其中$(A^T A)^{-1}A^T$称为$A$的伪逆$A^{\dagger}$,即

我们对A进行SVD分解,即 $ A = U \Sigma V^T $,则

其中$r = rank(A)$。 上面使用到了矩阵求逆交换律:$(EFG)^{-1} = G^{-1}F^{-1}E^{-1}$以及$U$和$V$的正交矩阵的性质。

所以以上最小二乘的解为:$\hat{x} = A^{\dagger} b = V\Sigma^{-1}U^T b$。

我们看一下最小二乘的error, 定义 $b$ 的估计值为 $\hat{b}$,它的具体形式如下:

所以$P_A = UU^T$,表示投影矩阵,它将$b$投影到了$A$的列空间的线性组合,或者说是值域即$span(A)$。所以:

$A\hat{x}$ 就是 $b$ 在 $span(A)$ 上的正交投影

再看一下error $e = b - \hat{b}$:

由于所奇异矩阵$U$的full形式为$U_{full} = [U | U_{\perp} ]$,则有:

所以: $P_{A\perp} = I - UU^T = U_{\perp}U_{\perp}^T$

所以最小二乘的误差为 $e = U_{\perp}U_{\perp}^Tb$

参考视频《Solving the Least-Squares Problem Using Geometry》

矩阵近似

祭上亲爱的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])

参考