引
这篇文章,讨论的是这样的一个问题:
M
=
L
0
+
S
0
M = L_0 + S_0
M=L0+S0
有这样的一个矩阵 M ∈ R n 1 × n 2 M \in \mathbb{R}^{n_1 \times n_2} M∈Rn1×n2,它由一个低秩的矩阵 L 0 L_0 L0和稀疏矩阵 S 0 S_0 S0混合而成,那么,能否通过 M M M将 L 0 L_0 L0和 S 0 S_0 S0分别重构出来?而且,最关键的是,也是这篇文章的亮点所在:对于 L 0 L_0 L0和 S 0 S_0 S0仅有一丢丢的微弱的假设。
一些微弱的假设:
关于
L
0
L_0
L0:
低秩矩阵
L
0
L_0
L0不稀疏,这个假设很弱,但有意义。因为如果
M
=
e
1
e
1
∗
M=e_1e_1^*
M=e1e1∗(文章采用
∗
*
∗作为转置符号,这里沿用这种写法),即只有左上角第一个元素非零,那么这个矩阵如何分解为低秩和稀疏呢,因为压根没法区分。作者引入incoherence condition(不连贯条件?抱歉,不晓得如何翻译):
假设
L
0
L_0
L0的SVD分解为:
L
0
=
U
Σ
V
∗
=
∑
i
=
1
r
σ
i
u
i
v
i
∗
L_0 = U\Sigma V^*=\sum \limits_{i=1}^r \sigma_i u_i v_i^*
L0=UΣV∗=i=1∑rσiuivi∗
其中
r
=
r
a
n
k
(
L
0
)
r = \mathrm{rank}(L_0)
r=rank(L0),
σ
i
\sigma_i
σi为奇异值,并且,
U
=
[
u
1
,
…
,
u
r
]
U = [u_1, \ldots, u_r]
U=[u1,…,ur],
V
=
[
v
1
,
…
,
v
r
]
V = [v_1, \ldots, v_r]
V=[v1,…,vr]。
incoherence condition:
分析这些条件可以发现,
U
,
V
U,V
U,V的每一行的模都不能太大,
U
V
∗
UV^*
UV∗的每个元素同样不能太大,所以最后结果
L
0
L_0
L0的各个元素的大小会比较均匀,从而保证不稀疏?
关于 S 0 S_0 S0:
类似地,需要假设
S
0
S_0
S0不低秩,论文另辟蹊径,假设
S
0
S_0
S0的基为
m
=
C
a
r
d
(
S
0
)
m=\mathrm{Card}(S_0)
m=Card(S0),其支撑(support)定义为:
Ω
=
{
(
i
,
j
)
∣
[
S
0
]
i
,
j
≠
0
}
\Omega = \{(i,j)| [S_0]_{i,j} \ne 0\}
Ω={(i,j)∣[S0]i,j̸=0}
论文假设
S
0
S_0
S0的支撑
Ω
\Omega
Ω从一个均匀分布采样。即,
S
0
S_0
S0不为零的项为
m
m
m,从
n
1
×
n
2
n_1\times n_2
n1×n2个元素中挑选
m
m
m个非零的组合有
C
n
1
×
n
2
m
\mathrm{C}_{n_1 \times n_2}^m
Cn1×n2m,论文假设每一种组合被采样的概率都是相同的。事实上,关于
S
0
S_0
S0元素的符号,论文在后面说明,符号固定或者随机(服从伯努利分布),都能将
L
0
,
S
0
L_0,S_0
L0,S0恢复出来。
问题的解决
论文反复强调,多么令人惊讶,多么不可思议,事实上的确如此,上述问题只需求解一个形式非常简单的凸优化问题:
其中
∥
L
∥
∗
=
∑
i
=
1
r
σ
i
(
L
)
\|L\|_* = \sum \limits_{i=1}^r \sigma_i(L)
∥L∥∗=i=1∑rσi(L)为其核范数。
论文给出了以下结论:
即:令
n
(
1
)
=
max
(
n
1
,
n
2
)
,
n
(
2
)
=
min
(
n
1
,
n
2
)
n_{(1)}= \max(n_1, n_2), n_{(2)} = \min (n_1, n_2)
n(1)=max(n1,n2),n(2)=min(n1,n2),对于任意矩阵
M
=
L
0
+
S
0
∈
R
n
1
×
n
2
M = L_0 + S_0 \in \mathbb{R}^{n_1 \times n_2}
M=L0+S0∈Rn1×n2,只要
L
0
L_0
L0和
S
0
S_0
S0满足上面提到的假设,那么就存在常数
c
c
c使得,PCP的解(即(1.1)的解,且
λ
=
1
/
n
(
1
)
\lambda= 1/ \sqrt{n_{(1)}}
λ=1/n(1))
L
^
\hat{L}
L^和
S
^
\hat{S}
S^有
1
−
c
n
(
1
)
−
10
1-c n_{(1)}^{-10}
1−cn(1)−10的概率重构
L
0
L_0
L0,
S
0
S_0
S0,即
L
^
=
L
0
\hat{L}=L_0
L^=L0,
S
^
=
S
0
\hat{S}=S_0
S^=S0。并且有下列性质被满足
r
a
n
k
(
L
0
)
≤
ρ
r
n
(
2
)
μ
−
1
(
log
n
(
1
)
)
−
2
,
m
≤
ρ
s
n
1
n
2
\mathrm{rank}(L_0) \le \rho_r n_{(2)} \mu^{-1} (\log n_{(1)})^{-2},m \le \rho_s n_1 n_2
rank(L0)≤ρrn(2)μ−1(logn(1))−2,m≤ρsn1n2。
这个结果需要说明的几点是,常数 c c c是根据 S 0 S_0 S0的支撑 Ω \Omega Ω决定的(根据后面的理论,实际上是跟 ∥ P Ω P T ∥ \|\mathcal{P}_{\Omega} \mathcal{P}_{T}\| ∥PΩPT∥有关),另外一点是, λ = 1 / n ( 1 ) \lambda = 1 / \sqrt{n_{(1)}} λ=1/n(1)。是的,论文给出了 λ \lambda λ的选择(虽然这个选择不是唯一的),而且,这个值是通过理论分析得到的。
关于问题 (1.1)的求解有很多现成的算法,这篇论文大量篇幅用于证明上述的理论结果,只在最后提到利用ALM(Augmented Largrange multiplier)方法有很多优势,算法如下:
其中
S
τ
[
x
]
=
s
g
n
(
x
)
max
(
∣
x
∣
−
τ
,
0
)
\mathcal{S}_{\tau}[x] = \mathrm{sgn} (x) \max (|x| - \tau, 0)
Sτ[x]=sgn(x)max(∣x∣−τ,0),
D
τ
(
X
)
=
U
S
τ
(
Σ
)
V
∗
,
X
=
U
Σ
V
∗
\mathcal{D}_{\tau} (X)=U\mathcal{S}_{\tau}(\Sigma)V^*,X=U\Sigma V^*
Dτ(X)=USτ(Σ)V∗,X=UΣV∗。
ALM算法实际上是交替最小化下式:
其中
<
A
,
B
>
:
=
T
r
(
A
T
B
)
<A, B>:= \mathrm{Tr(A^TB)}
<A,B>:=Tr(ATB)。
固定
L
,
Y
L, Y
L,Y,实际上就是最小化:
min
λ
∥
S
∥
1
+
<
Y
,
−
S
>
+
μ
2
T
r
(
S
T
S
−
2
S
T
(
M
−
L
)
)
\min \quad \lambda \|S\|_1 + <Y, -S> + \frac{\mu}{2}\mathrm{Tr}(S^TS -2S^T(M-L))
minλ∥S∥1+<Y,−S>+2μTr(STS−2ST(M−L))
其在
S
S
S处的导数(虽然有些点并不可导)为:
λ
s
g
n
(
S
)
−
Y
+
μ
S
−
μ
(
M
−
L
)
\lambda \mathrm{sgn}(S) - Y + \mu S-\mu(M-L)
λsgn(S)−Y+μS−μ(M−L)
实际上,对于每个
S
i
,
j
S_{i,j}
Si,j,与其相关的为:
a
S
i
,
j
2
+
b
S
i
,
j
+
c
∣
S
i
,
j
∣
,
α
>
0
a S_{i, j}^2 + b S_{i,j}+c|S_{i,j}|, \alpha > 0
aSi,j2+bSi,j+c∣Si,j∣,α>0
关于最小化这个式子,我们可以得到:
S
i
,
j
=
{
−
b
+
c
2
a
−
b
+
c
2
a
≥
0
−
b
−
c
2
a
−
b
−
c
2
a
≤
0
0
e
l
s
e
S_{i,j}= \left \{ \begin{array}{ll} -\frac{b+c}{2a} & -\frac{b+c}{2a} \ge 0 \\ -\frac{b-c}{2a} & -\frac{b-c}{2a} \le 0 \\ 0 & else \end{array} \right.
Si,j=⎩⎨⎧−2ab+c−2ab−c0−2ab+c≥0−2ab−c≤0else
实际上就是
S
i
,
j
=
S
c
/
2
a
(
−
b
/
2
a
)
S_{i,j}= \mathcal{S}_{c/2a}(-b/2a)
Si,j=Sc/2a(−b/2a),对
S
S
S拆解,容易得到固定
L
,
Y
L,Y
L,Y的条件下,
S
S
S的最优解为:
S
λ
/
μ
(
M
−
L
+
μ
−
1
Y
)
\mathcal{S}_{\lambda / \mu}(M-L+\mu^{-1}Y)
Sλ/μ(M−L+μ−1Y)
当
S
,
Y
S, Y
S,Y固定的时候:
实际上就是最小化下式:
∥
L
∥
∗
+
<
Y
,
−
L
>
+
μ
2
T
r
(
L
T
L
−
L
T
(
M
−
S
)
)
\|L\|_* + <Y, -L> + \frac{\mu}{2} \mathrm{Tr(L^TL-L^T(M-S))}
∥L∥∗+<Y,−L>+2μTr(LTL−LT(M−S))
观测其次梯度:
U
V
∗
+
W
−
Y
+
μ
L
−
μ
(
M
−
S
)
UV^*+W-Y+\mu L - \mu(M-S)
UV∗+W−Y+μL−μ(M−S)
其中
L
=
U
Σ
V
∗
,
U
∗
W
=
0
,
W
V
=
0
,
∥
W
∥
≤
1
L=U\Sigma V^*,U^*W=0,WV=0,\|W\| \le 1
L=UΣV∗,U∗W=0,WV=0,∥W∥≤1,
∥
X
∥
\|X\|
∥X∥为
X
X
X的算子范数。
最小值点应当满足
0
0
0为其次梯度,即:
L
=
−
μ
−
1
U
V
∗
−
μ
−
1
W
+
μ
−
1
Y
+
M
−
S
L = -\mu^{-1}UV^*-\mu^{-1}W+\mu^{-1}Y + M-S
L=−μ−1UV∗−μ−1W+μ−1Y+M−S
有解。
令
A
=
M
−
S
+
μ
−
1
Y
A = M-S+\mu^{-1}Y
A=M−S+μ−1Y
则:
L
=
−
μ
−
1
U
V
∗
−
μ
−
1
W
+
A
⇒
Σ
=
−
μ
−
1
I
+
U
∗
A
V
L = -\mu^{-1}UV^*-\mu^{-1}W+A \\ \Rightarrow \Sigma=-\mu^{-1}I+U^*AV
L=−μ−1UV∗−μ−1W+A⇒Σ=−μ−1I+U∗AV
假设
A
=
U
A
Σ
A
V
A
∗
A=U_A\Sigma_A V_A^*
A=UAΣAVA∗,那么:
Σ
=
−
μ
−
1
I
+
U
∗
U
A
Σ
A
V
A
∗
V
\Sigma=-\mu^{-1}I+U^*U_A\Sigma_A V_A^*V
Σ=−μ−1I+U∗UAΣAVA∗V
如果我们令
U
=
U
A
,
V
=
V
A
U=U_A,V=V_A
U=UA,V=VA,则:
Σ
=
Σ
A
−
μ
−
1
I
\Sigma = \Sigma_A-\mu^{-1}I
Σ=ΣA−μ−1I
但是我们知道,
Σ
\Sigma
Σ的对角元素必须大于等于0,所以我可以这样构造解:
假设
Σ
A
\Sigma_A
ΣA的对角元素,即奇异值降序排列
σ
1
(
A
)
≥
σ
2
(
A
)
…
σ
k
(
A
)
≥
μ
−
1
>
σ
k
+
1
(
A
)
≥
…
≥
σ
r
(
A
)
\sigma_1(A) \ge \sigma_2(A) \ldots \sigma_k(A) \ge \mu^{-1} > \sigma_{k+1}(A) \ge \ldots \ge \sigma_r(A)
σ1(A)≥σ2(A)…σk(A)≥μ−1>σk+1(A)≥…≥σr(A),相对于的向量为
u
1
(
A
)
,
…
,
u
r
(
A
)
u_1(A),\ldots,u_r(A)
u1(A),…,ur(A)和
v
1
A
,
…
,
v
r
(
A
)
v_1{A},\ldots, v_r(A)
v1A,…,vr(A)。我们令
L
=
∑
i
=
1
k
(
σ
i
(
A
)
−
μ
−
1
)
u
i
(
A
)
v
i
T
(
A
)
W
=
∑
i
=
k
+
1
r
μ
σ
i
(
A
)
u
i
(
A
)
v
i
T
(
A
)
L=\sum \limits_{i=1}^k (\sigma_i(A)-\mu^{-1})u_i(A)v_i^T(A) \\ W = \sum \limits_{i=k+1}^{r} \mu \sigma_i(A)u_i(A)v_i^T(A)
L=i=1∑k(σi(A)−μ−1)ui(A)viT(A)W=i=k+1∑rμσi(A)ui(A)viT(A)
容易验证
U
∗
W
=
0
,
W
V
=
0
U^*W=0,WV=0
U∗W=0,WV=0,又
μ
−
1
>
σ
i
(
A
)
,
i
>
k
\mu^{-1} > \sigma_i(A), i > k
μ−1>σi(A),i>k,所以
μ
σ
i
(
A
)
<
1
\mu \sigma_i(A) < 1
μσi(A)<1,所以
∥
W
∥
≤
1
\|W\| \le 1
∥W∥≤1也满足。同样的,次梯度为0的等式也满足,所以
L
L
L就是我们要找到点(因为这个是凸函数,极值点就是最小值点)。
理论
这里只介绍一下论文的证明思路。
符号:
去随机
论文先证明,如果
S
0
S_0
S0的符号是随机的(伯努利分布),在这种情况能恢复,那么
S
0
S_0
S0的符号即便不随机也能恢复。
Dual Certificates(对偶保证?)
引理2.4
引理2.4告诉我们,只要我们能找到一对
(
W
,
F
)
(W, F)
(W,F),满足
U
V
∗
+
W
=
λ
(
s
g
n
(
S
0
)
+
F
)
UV^*+W=\lambda(sgn(S_0)+F)
UV∗+W=λ(sgn(S0)+F)。
(
L
0
,
S
0
)
(L_0,S_0)
(L0,S0)是优化问题(1.1)的唯一解。但是不够,因为这个条件太强了,
W
,
F
W, F
W,F的构造是困难的。
引理2.5
进一步改进,有了引理2.5。
引理2.5的意义在哪呢?它意味着,我们只要找到一个
W
W
W,满足:
那么,
(
L
0
,
S
0
)
(L_0,S_0)
(L0,S0)就是问题
(
1.1
)
(1.1)
(1.1)的最优解,而且是唯一的。
Golfing Scheme
接下来,作者介绍一种机制来构造
W
W
W,并且证明这种证明能够大概率保证
W
W
W能够满足上述的对偶条件。作者将
W
=
W
L
+
W
S
W = W^L + W^S
W=WL+WS分解成俩部分,俩部分用不同的方式构造:
接下的引理和推论,总结起来便是最后的证明:
通过引理2.8和2.9,再利用范数的三角不等式,容易验证
W
=
W
L
+
W
S
W = W^L+W^S
W=WL+WS满足对偶条件。
最后提一下,定理的概率体现在哪里,在对偶条件的引理中,有这么一个假设
∥
P
Ω
P
T
∥
≤
1
\|\mathcal{P}_{\Omega}\mathcal{P}_T\| \le 1
∥PΩPT∥≤1和
∥
P
Ω
P
T
∥
≤
1
/
2
\|\mathcal{P}_{\Omega}\mathcal{P}_T\| \le 1/2
∥PΩPT∥≤1/2,这些条件并不一定成立,根据
Ω
\Omega
Ω的采样,成立的大概率的。
数值实验
按照论文里,普通的人工生成矩阵,实验下来和论文的无异。也的确, S 0 S_0 S0的规模大小(每个元素的大小)对能否恢复几乎没有影响,有影响的是 S 0 S_0 S0的稀疏度。实验下来,稀疏度越高(0基越小),恢复的越好,而且差距蛮大的。
我们的实验,
L
=
X
Y
L = XY
L=XY,其中
X
,
Y
∈
R
500
×
25
X,Y \in \mathbb{R}^{500 \times 25}
X,Y∈R500×25依标准正态分布生成的,
S
0
S_0
S0是每个元素有
ρ
\rho
ρ的概率为
100
100
100,
ρ
\rho
ρ的概率为
−
1
-1
−1,
1
−
2
ρ
1-2\rho
1−2ρ的概率为0。
ρ
\rho
ρ我们选了
1
/
10
,
1
/
15
,
…
,
1
/
55
1/10, 1/15, \ldots, 1/55
1/10,1/15,…,1/55。
下面的图,红色表示
∥
L
^
−
L
0
∥
F
/
∥
L
0
∥
F
\|\hat{L}-L_0\|_F/\|L_0\|_F
∥L^−L0∥F/∥L0∥F,蓝色表示
∥
S
^
−
S
0
∥
F
/
∥
S
0
∥
F
\|\hat{S}-S_0\|_F/\|S_0\|_F
∥S^−S0∥F/∥S0∥F。因为每种情况我只做了一次实验,所以结果有起伏,但是趋势是在的。
比较有意思的是图片的拆解,我在马路上拍了31张照片,按照论文的解释,通过将照片拉成向量,再混合成一个矩阵 M M M,然后通过优化问题解出 L ^ \hat{L} L^和 S ^ \hat{S} S^,那么 L ^ \hat{L} L^的行向量(我做实验的时候是将一个照片作为一个行向量,论文里的是列向量)重新展成图片应该是类似背景的东西,而 S ^ \hat{S} S^中对应的行向量应该是原先图片中会动的东西。
照片我进行了压缩(
416
×
233
416 \times 233
416×233)和灰度处理,处理
M
M
M用了1412次迭代,时间大概15分钟?比我预想的快多了。
我挑选其中比较混乱的图(就是车比较多的):
第一组:
M:
L:
S:
第二组:
M:
L:
S:
第三组:
M:
L:
S:
比较遗憾的是,
S
^
\hat{S}
S^中都出现了面包车,我以为面包车会没有的,结果还是出现了。
代码
"""
Robust Principal Component Analysis?
"""
import numpy as np
class RobustPCA:
def __init__(self, M, delta=1e-7):
self.__M = np.array(M, dtype=float)
self.__n1, self.__n2 = M.shape
self.S = np.zeros((self.__n1, self.__n2), dtype=float)
self.Y = np.zeros((self.__n1, self.__n2), dtype=float)
self.L = np.zeros((self.__n1, self.__n2), dtype=float)
self.mu = self.__n1 * self.__n2 / (4 * self.norm_l1)
self.lam = 1 / np.sqrt(max(self.__n1, self.__n2))
self.delta = delta
@property
def norm_l1(self):
"""返回M的l1范数"""
return np.sum(np.abs(self.__M))
@property
def stoprule(self):
"""停止准则"""
A = (self.__M - self.L - self.S) ** 2
sum_A = np.sqrt(np.sum(A))
bound = np.sqrt(np.sum(self.__M ** 2))
if sum_A <= self.delta * bound:
return True
else:
return False
def squeeze(self, A, c):
"""压缩"""
if c <= 0:
return A
B1 = np.where(np.abs(A) < c)
B2 = np.where(A >= c)
B3 = np.where(A <= -c)
A[B1] = [0.] * len(B1[0])
A[B2] = A[B2] - c
A[B3] = A[B3] + c
return A
def newsvd(self, A):
def sqrt(l):
newl = []
for i in range(len(l)):
if l[i] > 0:
newl.append(np.sqrt(l[i]))
else:
break
return np.array(newl)
m, n = A.shape
if m < n:
C = A @ A.T
l, u = np.linalg.eig(C)
s = sqrt(l)
length = len(s)
u = u[:, :length]
D_inverse = np.diag(1 / s)
vh = D_inverse @ u.T @ A
return u, s, vh
else:
C = A.T @ A
l, v = np.linalg.eig(C)
s = sqrt(l)
length = len(s)
v = v[:, :length]
D_inverse = np.diag(1 / s)
u = A @ v @ D_inverse
return u, s, v.T
def update_L(self):
"""更新L"""
A = self.__M - self.S + self.Y / self.mu
u, s, vh = np.linalg.svd(A) #or self.newsvd(A)
s = self.squeeze(s, 1 / self.mu)
s = s[np.where(s > 0)]
length = len(s)
if length is 0:
self.L = np.zeros((self.__n1, self.__n2), dtype=float)
elif length is 1:
self.L = np.outer(u[:, 0] * s[0], vh[0])
else:
self.L = u[:, :length] * s @ vh[:length]
def update_S(self):
"""更新S"""
A = self.__M - self.L + self.Y / self.mu
self.S = self.squeeze(A, self.lam / self.mu)
def update_Y(self):
"""更新Y"""
A = self.__M - self.L - self.S
self.Y = self.Y + self.mu * A
def process(self):
count = 0
while not (self.stoprule):
count += 1
assert count < 10000, "something wrong..."
self.update_L()
self.update_S()
self.update_Y()
print(count)
注意到,我自己写了一个newsvd的方法,这是因为,在处理图片的时候,用numpy的np.linalg.svd会报memoryerror,所以就自己简单调整了一下。