图像的低秩近似与Eckart–Young–Mirsky定理

引子

Hello 大家好呀，今天我们从理论的角度来分析一下为什么奇异值分解有助于降噪。

我们先回忆一下使用SVD分解进行降噪的过程：首先我们有一张被噪音污染的图片$A+E$

from  torchvision.io import read_image
import torch
import torchvision.transforms as transforms
import matplotlib.pyplot as plt

# Read in Images
img_path="cat.jpg"
cat=read_image(img_path)
cat=cat.float()
cat/=255

#Add noise
noise = (0.1**0.5)*torch.randn(cat.shape)
plt.imshow((cat+noise).permute(1,2,0))

我们对这张图片使用SVD分解 \(A+E=USV^T\)

选取一个阈值$r$,使得对角阵$S$中只保留比$r$大的那些值,记这个矩阵为$S_{clean}$: \(S_{clean}=\left[\begin{array}{cc} S^{\# r} & 0 \\ 0 & 0 \end{array}\right]\)

然后我们再把原先的U,V乘上去,即可获得去噪之后的图片

#Perform SVD
def fil(S,thres):
    s1=torch.tensor([ss if ss.item()>thres else 0. for ss in S[0,]])
    s2=torch.tensor([ss if ss.item()>thres else 0. for ss in S[1,]])
    s3=torch.tensor([ss if ss.item()>thres else 0. for ss in S[2,]])
    S_c=torch.stack([s1,s2,s3])
    print('Rankd for channel1: {}, channel2: {}, channel3: {}'.format(sum(S_c[0]>0).item(),sum(S_c[1]>0).item(),sum(S_c[2]>0).item()))
    return S_c

def get(S):
    a=U[0,]@torch.diag(S[0,])@V[0,]
    b=U[1,]@torch.diag(S[1,])@V[1,]
    c=U[2,]@torch.diag(S[2,])@V[2,]
    recover=torch.stack([a,b,c],dim=2)
    return recover
    
#Recover Clean Img
def whole(thres):
    U,S,V=torch.linalg.svd(cat+noise)
    S_c=fil(S,thres)
    recover=get(S_c)
    return recover
thres=1e1 
recover=whole(thres)
plt.imshow(recover)

奇异值的大小可以解读为矩阵$u_i v_i^T$的重要性，值越大说明这一部分在图像中越重要，值比较小的部分就对应无意义的噪音。

通过保留前k个较大的主要成分，我们可以在保留原始图片的同时达到一定程度的去噪效果。 \(A_k:=\sum_{i=1}^k \sigma_i u_i v_i^{\top}\)

那么这样做为什么会有效呢？我们其实需要证明通过奇异值的截断操作，我们能够使下式最小化，从而从数学上证明上述的操作有理论支撑。

这就引出了我们今天的主题：低秩逼近与Eckart–Young–Mirsky定理。

低秩逼近

我们先来讲一下低秩逼近是一个什么样的问题：

问题就是我们希望用一个秩不超过$r$的矩阵$\widehat{D}$去逼近原始矩阵$D$,使得他们俩之间的F范数最小。

考虑到图像本来就具有低秩的性质，在上述操作中我们使用截尾SVD来重建被噪音破坏的图像，也可以视为一种低秩逼近问题。

Eckart–Young–Mirsky定理就是证明了这个问题刚好可以被奇异值分解所解决。那么我们就一起来看看吧！

一句话来说，就是上述低秩逼近的解就是其SVD分解的秩$r$截尾解。

接下来我们就来看一下如何证明这个结论，从2范数和F范数的角度我们都可以证明上述截尾解就是最优解。

从2范数角度证明

首先我们先利用前面分好块的形式做一些计算：

\(=\sum_{i=1}^{\color{red}r} \sigma_i u_i v_i^{\top}+\sum_{i={\color{red}r+1}}^n \sigma_i u_i v_i^{\top}\) 我们记第一部分为$D_r$，则我们需要证明对任意秩小于等于r的$m\times n$矩阵$B_r$,都有

则此时秩r截断矩阵$D_r$为原优化问题的最优解。

我们先来计算等式左边的具体表达

所以我们只需证明右边不超过原矩阵$D$的第r+1个奇异值即可。

我们先对任意秩小于等于r的矩阵$B_r$做一个关于维数$r$的分解，$B_r=XY^T$，其中$X\in\mathbb{R}^{m\times r},Y^T\in\mathbb{R}^{r\times n}$。

如果$B_r$的秩刚好为$r$,直接用它的满秩分解即可，若秩小于$r$，则在它对应的满秩分解后补$\boldsymbol{0}$凑够$r$列即可。

对于$Y^T w=0$这个齐次方程组来说，由于这里至多只有r个线性无关的向量，根据施密特正交化的原理，很容易在正交矩阵$V\in\mathbb{R}^{n\times n}$的前r+1个线性无关且相互正交的正交基${v_1,\cdots,v_{r+1}}$中找到一个非零线性组合$w$，

使得$Y^Tw=0$。

并且我们可以对$w$的系数做一个规范化，使之落在单位球上: \(\|w\|_2=\gamma_1^2+\cdots+\gamma_{r+1}^2=1\)

由Cauchy-Schwarz不等式

由于 \(B_rw=X(Y^Tw)=X\mathbf{0}=\mathbf{0}\) 所以右边只剩下 \(\|Dw-\mathbf{0}\|^2_2\)

接下来我们来展开算一下这一块是多少，我们先把$V$和$V^T$用列向量的方式表示出来：

\(V_{n\times n}^T=\left[\begin{array}{ll} v_1^T\\ \vdots\\v_{r+1}^T\\ \hline \vdots\\ v_{n}^T \end{array}\right]\) 记$\boldsymbol{\gamma}=\left[\begin{array}{ll} \gamma_1,\cdots,\gamma_{r+1} \end{array}\right]^T\in\mathbb{R}^{(r+1)\times 1}$ \(w_{n\times 1}=V_{n\times (r+1)}^{\#}\boldsymbol{\gamma}_{(r+1)\times 1}=\left[\begin{array}{ll} v_1,\cdots,v_{r+1} \end{array}\right]\left[\begin{array}{ll} \gamma_1\\\vdots\\\gamma_{r+1} \end{array}\right]\)

由前面的SVD分解，$D=U\Sigma V^T$

通过使用结合律，我们先来计算$V^T V^{#}$,利用正交矩阵列向量互相正交且模长为1的特性，我们有

再把这一部分和$\boldsymbol{\gamma}$乘起来

$\left[\begin{array}{ll} \boldsymbol{I}{r+1}\ \boldsymbol{0}{(n-r-1)\times (r+1)} \end{array}\right]{n\times (r+1)}\left[\begin{array}{ll} \gamma_1\\vdots\\gamma{r+1} \end{array}\right]{(r+1)\times 1}=\left[\begin{array}{lll} \gamma_1
\vdots
\gamma{r+1}
\boldsymbol{0} \end{array}\right]_{n \times 1}$

现在$V^TV^{#}\boldsymbol{\gamma}$的部分都计算完了，我们就可以把计算的结果乘上前面的大$\Sigma$

现在原来的$Dw=U\Sigma V^T V^{#}\boldsymbol{\gamma}$就变成了

记$U=\left[\begin{array}{lll} u_1 \cdots u_m \end{array}\right]$

则$u_i$互相垂直且模长为1

所以有 \(\|Dw\|_2^2=\sigma_1^2\gamma_1^2 +\cdots+\sigma_{r+1}^2\gamma_{r+1}^2\)

由于$\gamma_1^2 +\cdots+\gamma_{r+1}^2=1$，且$\sigma_1\geq\sigma_2\geq\cdots\geq\sigma_{r+1}$,

相当于上式是关于奇异值的平方的一个凸组合，系数非负，且各成分都大于等于$\sigma_{r+1}^2$，

所以 \(\sigma_1^2\gamma_1^2 +\cdots+\sigma_{r+1}^2\gamma_{r+1}^2\geq \sigma_{r+1}^2\)

由此，我们就通过细致的计算证明了对任意的秩不超过$r$的矩阵$B_r$,

从而证明了SVD分解$\Sigma$矩阵秩r截尾的组合$U\Sigma^{#}_r V^T$就是原低秩优化问题的最优解(矩阵的2范数意义下)。

接下来我们再从F范数的意义下进行证明，表明SVD的这种截断解在F范数的意义下依然是最优解。

从F范数角度证明

证明的思路还是一样的，先计算

因为后半部分矩阵的奇异值就是$\sigma_{r+1}\cdots\sigma_{n}$，F范数平方的等价定义是奇异值的平方的和，所以上式就直接等于$\sum_{i=r+1}^n \sigma_{i}^2$。我们的目标就是证明对任意秩不超过r的矩阵$B_r$,有

我们先可以算出左边的展开式为

若想证明原式成立，需要证明

借助Weyl定理的奇异值版本，我们有如下定理可用

借用这个结论，我们可以令$A=D-B_r,B=B_r$，则有

取$i=1,j=r+1$，则$i+j-1=r+1$

由于$B_r$的秩不超过r，所以$\sigma_{r+1}(B_r)=0$

上式变为 \(\sigma_{r+1}(D)\leq \sigma_1(D-B_r)\)

同样地再取$i=2,\cdots,(n-r-2);j=r+1$ $\Rightarrow i+j+1=r+2,\cdots,n$

\begin{equation} \begin{aligned} \Rightarrow \sigma_{r+2}(D)&\leq \sigma_2(D-B_r)\\ &\cdots \\ \sigma_{n}(D)\leq& \sigma_{n-r-2}(D-B_r) \end{aligned} \end{equation}

把各项平方再相加，我们就证明了

所以SVD分解后的秩r截断解在F范数的意义下依然是原优化问题的最优解。

安可环节

上面用到的这个神奇的Weyl定理的奇异值版本怎么证呢？我们可以一起来看一看。我们要证明

先对AB各自做SVD分解，记为$A=U\Sigma V^T,B=X\Gamma Y^T$，因为我们研究的对象关于A是第i个特征值，关于B是第j个特征值，所以我们分别取出V和Y在i和j之后的列向量，记由这一部分正交基张成的子空间分别为 \(S'=\operatorname{Span}\{w_i,\cdots,w_n\}\) \(S''=\operatorname{Span}\{y_j,\cdots,y_n\}\) $dim(S’)=n-i+1,dim(S^{‘’})=n-j+1$

由维数定理 \(dim(S'\cap S^{''})=dim(S')+dim( S^{''})-dim(S'+ S^{''})\) \(=n-i+1+n-j+1-dim(S'+ S^{''})\) 由于$dim(S’+ S^{‘’})\leq n$,所以上式 \begin{equation} \begin{aligned} &\geq n-i+1+n-j+1-n\
&=n-(i+j-1)+1 \end{aligned} \end{equation} 由于我们指标$i+j-1$不能超出$n$，所以上式大于等于1，记 \begin{equation} \begin{aligned} \nu&=dim(S’\cap S^{‘’})\
&=n-(i+j-1)+1\geq 1 \end{aligned} \end{equation} 移项可得 \(n-\nu+1\leq i+j-1\) 由于指标小的奇异值大，所以 \begin{equation} \begin{aligned} \sigma_{i+j-1}(A+B)\leq \sigma_{n-\nu+1}(A+B) \end{aligned} \end{equation}

我们借助奇异值的一个等价定义

上式右边可以接着写成 \begin{equation} \begin{aligned} = \min _{\mathcal{S} \subset \mathbb{R}^n \atop dim(\mathcal{S})=\nu}\max _{x \in \mathcal{S}\atop\|x\|_2=1}\| (A+B) x\|_2 \end{aligned} \end{equation}

前面已经证明了$S’$与$S^{‘’}$相交的子空间$dim(S’\cap S^{‘’})=\nu$，所以它属于第一个$\min$考虑的集合。

由于$S’\cap S^{‘’}$算第一个$\min$下条件的特例，在它上面取max一定大于等于全局最小值，所以上式可以放大为

\begin{equation} \begin{aligned} \leq \max _{x \in S’\cap S^{‘’}\atop\|x\|_2=1}\| (A+B) x\|_2 \end{aligned} \end{equation}

再由范数的三角不等式放缩为 \begin{equation} \begin{aligned} \leq \max _{x \in S’\cap S^{‘’}\atop\|x\|_2=1}\| A x\|_2+\max _{x \in S’\cap S^{‘’}\atop\|x\|_2=1}\| B x\|_2 \end{aligned} \end{equation}

由于$(S’\cap S^{‘’})\subset S’,(S’\cap S^{‘’})\subset S^{‘’}$，在一个更大的子空间里面取最大值一定比在其所包含的子空间里面取最大值大。

把第一项的max的空间范围放大到$S’$，第二项的空间范围放大到$S^{‘’}$ \begin{equation} \begin{aligned} \leq \max _{x \in S’\atop\|x\|_2=1}\| A x\|_2+\max _{x \in S^{‘’}\atop\|x\|_2=1}\| B x\|_2 \end{aligned} \end{equation}

我们看第一项$x \in S’=\operatorname{Span}{w_i,\cdots,w_n}$,由于${w_i,\cdots,w_n}$是对应奇异值$\sigma_{i}$之后的奇异值向量，他们互相正交且奇异值依次减小。

我们容易得知在这个空间里能取到的最大值就是$\sigma_i(A)$，同理有$\sigma_j(B)$，所以上式 \begin{equation} \begin{aligned} = \sigma_i(A)+\sigma_j(B) \end{aligned} \end{equation}

串联起最头和最尾就得到了

定理得证。