📝笔记：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}$都是正交矩阵，即$U^{T}U=I,V^{T}V=I$，即列向量是正交的单位向量，$U$称为left single vectors，$V$称为right single vectors；$\Sigma \in {\mathbb{R}}^{m\times n}$的对角阵（奇异值），奇异值$\ge 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^{T}A$

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

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

计算 $A{A}^T$

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

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

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

降维

$${\Sigma }_k={\left[\begin{array}{cccc}{\sigma }_1& 0& \dots & 0\\ 0& {\sigma }_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& & {\sigma }_k\end{array}\right]}_{k\times k}$$

实战

数字例子

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

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

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

计算$A^{T}A$以及$A{A}^T$：

$$A^{T}A=\left[\begin{array}{ccccc}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{array}\right]$$ $$A{A}^T=\left[\begin{array}{ccc}87& 51& 74\\ 51& 54& 71\\ 74& 71& 147\end{array}\right]$$

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

奇异值：

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

最小二乘 least square问题

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

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

通常的解法如下：

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)$

$$\begin{array}{rl}Ax& =b\\ A^{T}Ax& =A^{T}b\\ x& =(A^{T}A{)}^{-1}{A}^{T}b\end{array}$$

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

$$A^{†}=(A^{T}A{)}^{-1}{A}^T$$

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

$$\begin{array}{rl}{A}^{†}& =(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\end{array}$$

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

所以以上最小二乘的解为：

$$\hat{x}=A^{†}b=V{\Sigma }^{-1}{U}^{T}b$$

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

$$\begin{array}{rl}\hat{b}& =A\hat{x}\triangleq P_{A}b\\ & =A{A}^{†}b\\ & =AV{\Sigma }^{-1}{U}^{T}b\\ & =U\Sigma V^{T}V{\Sigma }^{-1}{U}^{T}b\\ & =U{U}^{T}b\end{array}$$

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

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

再看一下error：

$$\begin{array}{rl}e& =b-\hat{b}\\ & =b-U{U}^{T}b\\ & =(I-U{U}^T)b\triangleq P_{A\perp }b\end{array}$$

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

$$U_{\text{full}}^T{U}_{\text{full}}=U_{\text{full}}{U}_{\text{full}}^T=I=U{U}^T+U_{\perp }{U}_{\perp }^T$$

进而最小二乘的误差为：

$$e=P_{A\perp }b=U_{\perp }{U}_{\perp }^{T}b$$

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

矩阵近似

祭上亲爱的Battle Angel Alita。

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

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

所以：

$$P_{A\perp }=I-U{U}^T=U_{\perp }{U}_{\perp }^T$$

代码

%% 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

%% 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

%% 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;
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,课程列表