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

Posted on 2019-08-18 Edited on 2025-01-27 In SLAM Waline: Views:

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

$\begin{matrix} (1) & A = U Σ V^{T} \end{matrix}$

其中， $U \in R^{m \times m}$ ， $V \in R^{n \times n}$ 都是正交矩阵，即 $U^{T} U = I, V^{T} V = I$ ，即列向量是正交的单位向量， $U$ 称为left single vectors， $V$ 称为right single vectors； $Σ \in R^{m \times n}$ 的对角阵（奇异值），奇异值 $\geq 0$ 且按照顺序降序排列。

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

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

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

形式： skinny SVD 与 full SVD

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

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

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

计算 $A^{T} A$

$\begin{matrix} (2) & A^{T} A = (U Σ V^{T})^{T} U Σ V^{T} = V Σ^{T} U^{T} U Σ V^{T} = V Σ^{T} Σ V^{T} \end{matrix}$ 可见， $V$ 正是矩阵 $A^{T} A$ 的特征向量，而 $Σ^{T} Σ$ 为矩阵 $A^{T} A$ 的特征值。

计算 $A A^{T}$

$\begin{matrix} (3) & A A^{T} = U Σ V^{T} (U Σ V^{T})^{T} = U Σ V^{T} V Σ^{T} U^{T} = U Σ Σ^{T} U^{T} \end{matrix}$ 可见， $U$ 正是矩阵 $A A^{T}$ 的特征向量，而 $Σ Σ^{T}$ 为矩阵 $A A^{T}$ 的特征值。

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

降维

$\begin{matrix} (4) & Σ = {[\begin{array}{cccccccc} σ_{1} & 0 & \dots & 0 & 0 & 0 & \dots & 0 \\ 0 & σ_{2} & \dots & 0 & 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ & 0 & 0 & \dots & 0 \\ 0 & 0 & σ_{k} & 0 & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 & 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & \dots & 0 & 0 & 0 & \dots & 0 \end{array}]}_{m \times n} \Rightarrow {[\begin{array}{cccc} σ_{1} & 0 & \dots & 0 \\ 0 & σ_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & σ_{k} \end{array}]}_{k \times k} \end{matrix}$

其中 $σ_{1} \geq σ_{2} \geq . . . σ_{k} > 0$ ，将 $Σ$ 中主对角线为0的部分删去，同样的 $U, V$ 对应的部分删去，SVD分解就变成了下图的形式。

实战

数字例子

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

$\begin{matrix} (5) & A = [\begin{matrix} 1 & 4 & 3 & 5 & 6 \\ 2 & 3 & 4 & 5 & 0 \\ 7 & 4 & 0 & 9 & 1 \end{matrix}] \end{matrix}$ 计算 $A^{T} A$ 以及 $A A^{T}$ ： $\begin{matrix} (6) & A^{T} A = [\begin{matrix} 54 & 38 & 11 & 78 & 13 \\ 38 & 41 & 24 & 71 & 28 \\ 11 & 24 & 25 & 35 & 18 \\ 78 & 71 & 35 & 131 & 39 \\ 13 & 28 & 18 & 39 & 37 \end{matrix}] A A^{T} = [\begin{matrix} 87 & 51 & 74 \\ 51 & 54 & 71 \\ 74 & 71 & 147 \end{matrix}] \end{matrix}$

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

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

奇异值：

$Σ^{T} Σ = Diag (238.2878, 37.3715, 12.3407) \Rightarrow Σ = Diag (15.4366, 6.1132, 3.5129)$

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

最小二乘 least square问题

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

问题： $A x = b$ ，求解 $x$ 。

通常的解法如下：

$\begin{matrix} (7) & \begin{aligned} A x & = b \\ A^{T} A x & = A^{T} b \\ x & = (A^{T} A)^{- 1} A^{T} b \end{aligned} \end{matrix}$

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

$\begin{matrix} (8) & A^{†} = (A^{T} A)^{- 1} A^{T} \end{matrix}$

我们对A进行SVD分解，即 $A = U Σ V^{T}$ ，则

$\begin{matrix} (9) & \begin{aligned} A^{†} & = (A^{T} A)^{- 1} A^{T} \\ = (V Σ^{T} U^{T} U Σ V^{T})^{- 1} V Σ^{T} U^{T} \\ = (V Σ^{T} Σ V^{T})^{- 1} V Σ^{T} U^{T} \\ = (V Σ^{2} V^{T})^{- 1} V Σ^{T} U^{T} \\ = V^{- T} Σ^{- 2} V^{- 1} V Σ^{T} U^{T} \\ = V Σ^{- 1} U^{T} \\ = \sum_{i = 1}^{r} \frac{1}{σ_{i}} v_{i} u_{i}^{T} d \end{aligned} \end{matrix}$

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

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

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

$\begin{matrix} (10) & \begin{aligned} \hat{b} & = A \hat{x} b ≜ P_{A} b \\ = A A^{†} b \\ = A V Σ^{- 1} U^{T} b \\ = U Σ V^{T} V Σ^{- 1} U^{T} b \\ = U U^{T} b \end{aligned} \end{matrix}$

所以 $P_{A} = U U^{T}$ ，表示投影矩阵，它将 $b$ 投影到了 $A$ 的列空间的线性组合，或者说是值域即 $s p a n (A)$ 。所以:

$A \hat{x}$ 就是 $b$ 在 $s p a n (A)$ 上的正交投影

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

$\begin{matrix} (11) & \begin{aligned} e & = b - \hat{b} \\ = b - U U^{T} b \\ = (I - U U^{T}) b ≜ P_{A ⊥} b \end{aligned} \end{matrix}$

由于所奇异矩阵 $U$ 的full形式为 $U_{f u l l} = [U | U_{⊥}]$ ，则有：

$\begin{matrix} (12) & U_{f u l l}^{T} U_{f u l l} = U_{f u l l} U_{f u l l}^{T} = I = U U^{T} + U_{⊥} U_{⊥}^{T} \end{matrix}$

所以: $P_{A ⊥} = I - U U^{T} = U_{⊥} U_{⊥}^{T}$

进而最小二乘的误差为 $e = U_{⊥} U_{⊥}^{T} b$

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

矩阵近似

祭上亲爱的Battle Angel Alita。

原始图像尺寸 $1440$ ，我们可以对该图像做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
Singular Value Decomposition (SVD) video1，video2
李宏毅关于SVD的介绍，PPT,课程列表