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

Posted on Edited on In Waline: Views:

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}\)都是正交矩阵，即\(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\)

\[ 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] \]

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

1
2
3
4
5
6
7
8
9
10
11
12
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\)。

通常的解法如下：

\[ \begin{aligned} A x &= b \\ A^T A x &= A^T b \\ x &= (A^T A)^{-1}A^T b \end{aligned} \]

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

\[ A^{\dagger} = (A^T A)^{-1}A^T \]

我们对A进行SVD分解，即 \(A = U \Sigma V^T\)，则

\[ \begin{aligned} A^{\dagger} &= (A^T A)^{-1}A^T \\ &= (V \Sigma^T U^T U \Sigma V^T)^{-1}V \Sigma^T U^T \\ &= (V \Sigma^T \Sigma V^T)^{-1}V \Sigma^T U^T \\ &= (V \Sigma^2 V^T)^{-1}V \Sigma^T U^T \\ &= V^{-T}\Sigma^{-2}V^{-1}V \Sigma^T U^T \\ &= V\Sigma^{-1}U^T \\ &= \sum_{i=1}^{r}\frac{1}{\sigma_i}v_{i}u_{i}^{T}d \end{aligned} \]

其中\(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}\)，它的具体形式如下：

\[ \begin{aligned} \hat{b} &= A\hat{x} b \triangleq P_A b \\ &= AA^{\dagger} b \\ &= AV\Sigma^{-1}U^T b \\ &= U \Sigma V^TV\Sigma^{-1}U^T b \\ &= UU^T b \end{aligned} \]

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

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

再看一下error \(e = b - \hat{b}\)：

\[ \begin{aligned} e &= b - \hat{b} \\ &= b - UU^T b \\ &= (I - UU^T) b \triangleq P_{A\perp}b \end{aligned} \]

由于所奇异矩阵\(U\)的full形式为\(U_{full} = [U | U_{\perp} ]\)，则有：

\[U_{full}^TU_{full} = U_{full}U_{full}^T = I = UU^T + U_{\perp}U_{\perp}^T\]

所以: \(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 $，我们可以对该图像做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
103
%% 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])

参考