图片压缩:SVD(奇异值分解)之近似矩阵
写这篇博客的初衷的是想让读者浅尝一下SVD的魅力。 处于大数据时代,真的是,“SVD在手,天下我有”。最关键的是,SVD很容易理解,简单的线性代数知识就够了。
矩阵的秩
对于一个矩阵,比方如下:
A
=
[
6
6
−
2
4
−
8
−
15
−
15
5
−
10
20
12
12
−
4
8
−
16
15
15
−
5
10
−
20
6
6
−
2
4
−
8
]
5
×
5
\mathbf{A}= \begin{bmatrix} 6 & 6 & -2 & 4 & -8\\ -15 & -15 & 5 & -10 & 20\\ 12 & 12 & -4 & 8 & -16\\ 15 & 15 & -5 & 10 & -20\\ 6 & 6 & -2 & 4 & -8 \end{bmatrix}_{5\times5}
A=
6−15121566−1512156−25−4−5−24−108104−820−16−20−8
5×5
正常情况下,我们会以
25
25
25 个字节的形式存储下来。其实呢,这个矩阵的秩为
1
1
1,是我随机生成的两个长度为
5
5
5 向量相乘得来的:
a
1
=
[
2
,
−
5
,
4
,
5
,
2
]
T
,
a
2
=
[
3
,
3
,
−
1
,
2
,
−
4
]
,
A
=
a
1
×
a
2
\mathbf{a}_1 = [2, -5, 4, 5, 2]^{T}, \quad \mathbf{a}_2 = [3, 3, -1, 2, -4], \quad \mathbf{A} = \mathbf{a}_1\times\mathbf{a}_2
a1=[2,−5,4,5,2]T,a2=[3,3,−1,2,−4],A=a1×a2
所以存储起来,
10
10
10 个字节就够了。想象一个,如果是一个
n
×
m
n\times m
n×m 的矩阵呢,特别是处于大数据时代的我们,
n
n
n 和
m
m
m 可以随随便便就上万,甚至上亿,相似的方式存储这些矩阵,将会节省很多存储空间。
跟着这个思路,我们知道,对于一个大型的矩阵,如果它不是满秩,我们是可以用以上的方式来降低存储量。然而,现实生活中,这种矩阵, 比方说是每个列是某一个时刻的传感器数据,它虽然很可能不是满秩,但是也有可能接近满秩。 庆幸的是,这种大型的矩阵所包含的信息中,有一些是可以忽略的。如果我们可以求得这个矩阵的低阶近似矩阵,我们便可以达到最初的目的。
近似矩阵
为了降低存储量,我们的目标是找到一个大型矩阵的低阶近似矩阵。对于近似的理解,在向量的表示中,我们一般都是用欧几里德范数(Euclidean norm)或者欧几里德距离来衡量两个向量是否接近。在矩阵中,我们可以用弗罗贝尼乌斯范数(Frobenius norm)来衡量。所以,对于一个大型矩阵,比方说 X \mathbf{X} X,我们希望找到一个低阶的近似矩阵 X ~ \tilde{\mathbf{X}} X~。 ∥ X − X ~ ∥ F \left \| \mathbf{X}- \tilde{\mathbf{X}}\right \|_F X−X~ F 越小,两个矩阵越近似。
SVD
SVD 提供了一种系统的方法来根据主要特征确定高维数据的低阶近似。 这种技术是基于数据的,因为主要特征纯粹是从数据中发现的,不需要额外的信息。 SVD 在数值上是稳定的,并根据由数据内的主要相关性定义的新坐标系对数据进行分层。 最重要的一点是,SVD 对于任何矩阵都保证存在。
对于实数矩阵
X
\mathbf{X}
X (如果是复数的话,以下的转置都是共轭转置),我们把它表示为:
X
=
[
∣
∣
∣
x
1
x
2
⋯
x
m
∣
∣
∣
]
\mathbf{X}=\begin{bmatrix} | & | & & |\\ \mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_m\\ | & | & & | \end{bmatrix}
X=
∣x1∣∣x2∣⋯∣xm∣
其中,
x
i
\mathbf{x}_i
xi 是长度为
n
n
n 的列向量(通常我们遇到的数据中,
n
≫
m
n\gg m
n≫m,我们以下的表述中也基于这个假设)。存在且唯一存在一个SVD可以将矩阵
X
\mathbf{X}
X 分解为:
X
=
U
Σ
V
T
\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}
X=UΣVT
其中,
U
\mathbf{U}
U 和
V
\mathbf{V}
V 分别是
n
×
n
n\times n
n×n 和
m
×
m
m\times m
m×m 的正交矩阵 (
U
T
U
=
I
n
\mathbf{U}^{T}\mathbf{U} = I_n
UTU=In,
V
T
V
=
I
m
\mathbf{V}^{T}\mathbf{V} = I_m
VTV=Im);
Σ
\mathbf{\Sigma}
Σ 是
n
×
m
n\times m
n×m 的矩阵,且除了对角线上是非负数外,其余元素都是
0
0
0。
由于
n
⩾
m
n\geqslant m
n⩾m,
Σ
\mathbf{\Sigma}
Σ 对角线上最多有
m
m
m 个非零数值,我们也可以将
Σ
\mathbf{\Sigma}
Σ 表示为:
Σ
=
[
Σ
^
0
]
\mathbf{\Sigma}=\begin{bmatrix} \hat{\mathbf{\Sigma}}\\ \mathbf{0} \end{bmatrix}
Σ=[Σ^0]
进而,我们可以将
X
\mathbf{X}
X 表示为:
X
=
U
Σ
V
T
=
[
U
^
U
^
⊥
]
[
Σ
^
0
]
V
T
=
U
^
Σ
^
V
T
\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}=[\hat{\mathbf{U}} \quad \hat{\mathbf{U}}^{\bot}]\begin{bmatrix} \hat{\mathbf{\Sigma}}\\ \mathbf{0} \end{bmatrix}\mathbf{V}^{T}=\hat{\mathbf{U}}\hat{\mathbf{\Sigma}}\mathbf{V}^{T}
X=UΣVT=[U^U^⊥][Σ^0]VT=U^Σ^VT
也即简化的SVD(economy SVD)。值得注意的是,
Σ
\mathbf{\Sigma}
Σ 对角线上数值是从大到小排列的。我们将上式展开,得到:
X
=
∑
i
=
1
m
σ
i
u
i
v
i
T
=
σ
1
u
1
v
1
T
+
σ
2
u
2
v
2
T
+
⋯
+
σ
m
u
m
v
m
T
\mathbf{X}=\sum_{i=1}^{m}\sigma_i\mathbf{u}_i\mathbf{v}_i^T = \sigma_1\mathbf{u}_1\mathbf{v}_1^T+\sigma_2\mathbf{u}_2\mathbf{v}_2^T+\cdots +\sigma_m\mathbf{u}_m\mathbf{v}_m^T
X=i=1∑mσiuiviT=σ1u1v1T+σ2u2v2T+⋯+σmumvmT
近似定理(Eckart-Young theorem):
矩阵
X
\mathbf{X}
X 的最优
r
r
r 阶近似矩阵,
X
opt
r
=
arg
min
X
~
,
s.t. rank
(
X
~
)
=
r
∥
X
−
X
~
∥
F
\mathbf{X}^{r}_\text{opt}=\underset{\tilde{\mathbf{X}}, \text{ s.t. rank}(\tilde{\mathbf{X}})=r}{\arg\min}\left \| \mathbf{X} - \tilde{\mathbf{X}} \right \|_F
Xoptr=X~, s.t. rank(X~)=rargmin
X−X~
F, 是截断的
r
r
r 阶SVD,也即:
X
opt
r
=
∑
i
=
1
r
σ
i
u
i
v
i
T
=
σ
1
u
1
v
1
T
+
σ
2
u
2
v
2
T
+
⋯
+
σ
r
u
r
v
r
T
\mathbf{X}^{r}_\text{opt}=\sum_{i=1}^{r}\sigma_i\mathbf{u}_i\mathbf{v}_i^T = \sigma_1\mathbf{u}_1\mathbf{v}_1^T+\sigma_2\mathbf{u}_2\mathbf{v}_2^T+\cdots +\sigma_r\mathbf{u}_r\mathbf{v}_r^T
Xoptr=i=1∑rσiuiviT=σ1u1v1T+σ2u2v2T+⋯+σrurvrT
意味着,如果你期望近似矩阵为 10 10 10 阶的话,那么最近似的矩阵就是: ∑ i = 1 10 σ i u i v i T \sum_{i=1}^{10} \sigma_i\mathbf{u}_i\mathbf{v}_i^T ∑i=110σiuiviT。
图片压缩
把灰度(greyscale)图片看作是一个矩阵,我想大家都很容易理解这一点。下图中,图一是一个
1600
×
1200
1600\times1200
1600×1200 的灰度图片(丹麦的美人鱼)。
通过上述的近似矩阵方法,我们可以分别获取
20
20
20 阶、
100
100
100 阶 和
250
250
250 阶的低阶近似矩阵(基于MATLAB):
A = imread('./littlemermaid.jpeg');
X = double(rgb2gray(A)); %convert RGB to gray
nx = size(X,1); ny = size(X,2);
[U, S, V] = svd(X,'econ');
figure, subplot(1,4,1)
imagesc(X), axis off, colormap gray
title('Greyscale')
plotind = 2;
for r = [20 100 250] %truncation value
Xapprox = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; %approximate image
subplot(1,4,plotind), plotind = plotind+1;
imagesc(Xapprox), axis off
title(['r=',num2str(r,'%d'),', ', num2str(100*r*(nx+ny)/(nx*ny),'%2.2f'),'% storage'])
end
对比可以看出,低阶近似矩阵需要的存储量大大降低了,而且在实际中也经常有它的用处(比方说微信可以选择发送原图或者压缩的图片,这里就可以申明使用SVD压缩。当然,也可以用其它的方法压缩)。