低秩稀疏矩陣恢復|ADM(IALM)算法

---

一曲新詞酒一杯，去年天氣舊亭臺。夕陽西下幾時回？
迫不得已花落去，似曾相識燕歸來。小園香徑獨徘徊。
———《浣溪沙·一曲新詞酒一杯》——晏殊

更多精彩內容請關注微信公衆號 「優化與算法」

上一期介紹了低秩矩陣填充問題，這一期介紹一下低秩稀疏矩陣恢復問題。

1. 低秩矩陣恢復

將一個矩陣 \(\bf{D}~(\bf {D} = \bf {A_0} +\bf E_0)\) 分解爲一個低秩矩陣部分 \(\bf{A}\) 和一個獨立同分布的高斯矩陣 \(\bf{E}\) 的問題是經典的主成分分析（PCA）問題，能夠經過對矩陣 \(\bf{D}\) 進行奇異值分解獲得最優解。

然而，當矩陣 \(\bf{E_0}\) 爲稀疏的噪聲矩陣時，PCA再也不適用於解決這個問題。此時 ，將一個矩陣 \(\bf{D}~(\bf {D} = \bf {A_0} +\bf E)\) 分解爲一個低秩矩陣部分 \(\bf{A}\) 和一個稀疏矩陣部分 \(\bf{E}\) 的問題能夠建模爲下述優化問題：

\[\begin{array}{l} \mathop {\min }\limits_{{\bf{A}},{\bf{E}}} ~~~rank({\bf{A}}) + \lambda {\left\| {\bf{E}} \right\|_0} \\ s.t.~~~{\bf{D}} = {\bf{A}} + {\bf{E}} \\ \end{array}~~~~~(1)\]

其中 \({\bf{D}},{\bf{A}},{\bf{E}},{{\bf{A}}_0},{{\bf{E}}_0}{ \in \mathbb{R}^{m \times n}}\)，\(\bf D\) 是觀測矩陣。(1)式中 \(rank(\bf A)\) 和 \({\left\| {\bf{E}} \right\|_0}\) 都是非線性且非凸的，優化起來很是困難，這個問題也被稱爲主成分追蹤（Principal component pursuit, PCP）。幸運的是咱們提早知道一些先驗信息，即矩陣 \(\bf A\) 是低秩的且矩陣 \(\bf E\) 是稀疏的，從上一期介紹的關於矩陣填充理論中可知，矩陣的秩和 \(\ell_0\) 範數問題均可以進行凸鬆弛，從而爲求解上述問題提供了途徑。因爲矩陣的核範數是矩陣秩的凸包絡，矩陣的（1,1）範數是矩陣0範數的凸包絡，所以能夠將問題(1)鬆弛爲以下凸優化問題：

\[\begin{array}{l} \mathop {\min }\limits_{{\bf{A}},{\bf{E}}}~~~ {\left\| {\bf{A}} \right\|_*} + \lambda {\left\| {\bf{E}} \right\|_{1,1}} \\ s.t.~~~{\bf{D}} = {\bf{A}} + {\bf{E}} \\ \end{array}~~~~~~(2)\]

求解式(2)也稱爲魯棒主成分分析(RPCA)。

文獻[1]中指出，只要低秩矩陣 \(\bf{A_0}\) 的奇異值分佈合理且稀疏矩陣的非零元素均勻分佈，那麼凸優化問題PCP就可以以接近1的機率從未知的任意偏差中恢復出原始低秩矩陣 \(\bf A_0\) 來。

求解(2)式的算法能夠分爲以下幾大類：

1. 迭代閾值算法
   對於PCP問題時，迭代閾值算法(Iterative Thresholding, IT) 經過交替更新矩陣 \(\bf A\) 和 \(\bf E\) 來求解。IT算法的迭代形式簡單且收斂，但它的收斂速度比較慢，且難以選取合適的步長。所以，IT算法具備有限的應用範圍。
2. 加速近端梯度算法
   加速近端梯度算法(Accelerated Proximal Gradient, APG)的主要思想是利用了Nesterov加速，使算法可以快速收斂。
3. 對偶方法
   將原問題轉化爲其對偶問題(非線性、非光滑)，並使用最速上升法等能夠求解。對偶方法比APG算法具備更好的可擴展性，這是由於在每次迭代過程當中對偶方法不須要矩陣的徹底奇異值分解。
4. 增廣拉格朗日乘子法

這些方法都很是經典，這裏再也不細述，總的來講，只要將問題轉化爲凸問題，就有一大堆方法能夠用來求解。這裏僅介紹一種增廣拉格朗日乘子算法，即交替方向方法(Alternating direction methods, ADM)，也稱爲不精確拉格朗日乘子法(Inexact ALM, IALM)。
下面給出上述幾種算法的比較(數據來源於網絡)

2. 交替方向算法(ADM)

對於優化問題(2)，首先構造增廣拉格朗日函數：

\[L({\bf{A}},{\bf{E}},{\bf{Y}},u) = {\left\| {\bf{A}} \right\|_*} + \lambda {\left\| {\bf{E}} \right\|_{1,1}} + \left\langle {{\bf{Y}},{\bf{D}} - {\bf{A}} - {\bf{E}}} \right\rangle + \frac{u}{2}\left\| {{\bf{D}} - {\bf{A}} - {\bf{E}}} \right\|_F^2~~~(3) \]

當 \({\bf{Y}} = {{\bf{Y}}_k},u = {u_k}\) 時，使用交替方法求解塊優化問題：

\[\mathop {\min }\limits_{{\bf{A}},{\bf{E}}} L({\bf{A}},{\bf{E}},{{\bf{Y}}_k},{u_k})~~~(4) \]

使用精確拉格朗日乘子法(Exact ALM, EALM)交替迭代矩陣 \(\bf A\) 和 \(\bf E\)，直到知足終止條件爲止。若 \({\bf{E}} = {\bf{E}}_{k + 1}^j\)，則

\[\begin{array}{l} {\bf{A}}_{k + 1}^{j + 1} = \arg \mathop {\min }\limits_{\bf{A}} L({\bf{A}},{\bf{E}}_{k + 1}^j,{{\bf{Y}}_k},{u_k}) \\ = \arg \mathop {\min }\limits_{\bf{A}} {\left\| {\bf{A}} \right\|_*} + \frac{{{u_k}}}{2}\left\| {{\bf{A}} - ({\bf{D}} - {\bf{E}}_{k + 1}^j + \frac{{{{\bf{Y}}_k}}}{{{u_k}}})} \right\|_F^2 \\ = {D_{\frac{1}{{{u_k}}}}}({\bf D} - {\bf{E}}_{k + 1}^j + \frac{{{{\bf{Y}}_k}}}{{{u_k}}}) \\ \end{array}~~~(5)\]

再根據 \({\bf{A}}_{k + 1}^{j + 1}\) 更新矩陣 \(\bf E\)：

\[\begin{array}{l} {\bf{E}}_{k + 1}^{j + 1} = \arg \mathop {\min }\limits_{\bf{E}} L({\bf{A}}_{k + 1}^{j + 1},{\bf{E}},{{\bf{Y}}_k},{u_k}) \\ = \arg \mathop {\min }\limits_{\bf{E}} \lambda {\left\| {\bf{E}} \right\|_{1,1}} + \frac{{{u_k}}}{2}\left\| {{\bf{E}} - ({\bf{D}} - {\bf{A}}_{k + 1}^{j + 1} + \frac{{{{\bf{Y}}_k}}}{{{u_k}}})} \right\|_F^2 \\ = {S_{\frac{\lambda }{{{u_k}}}}}({\bf D} - {\bf{A}}_{k + 1}^{j + 1} + \frac{{{{\bf{Y}}_k}}}{{{u_k}}}) \\ \end{array}~~~(6)\]

記 \({\bf{A}}_{k + 1}^{\rm{*}}\) 和 \({\bf{E}}_{k + 1}^{\rm{*}}\) 分別爲 \({\bf{A}}_{k + 1}^{j + 1}\) 和 \({\bf{E}}_{k + 1}^{j + 1}\) 的精確值，則矩陣 \(\bf Y\) 的更新公式爲：

\[{{\bf{Y}}_{k{\rm{ + }}1}}{\rm{ = }}{{\bf{Y}}_k}{\rm{ + }}{u_k}({\bf{D}} - {\bf{A}}_{k + 1}^{\rm{*}} - {\bf{E}}_{k + 1}^{\rm{*}})~~~(7) \]

參數 \({u_k}\) 能夠更新以下：

\[{u_{k + 1}} = \left\{ {\begin{array}{*{20}{c}} {\rho {u_k}~~~\frac{{{u_k}{{\left\| {{\bf{E}}_{k + 1}^{\rm{*}}{\rm{ - }}{\bf{E}}_k^{\rm{*}}} \right\|}_F}}}{{{{\left\| {\bf{D}} \right\|}_F}}} < \varepsilon } \\ {{u_k}~~~~~~~~otherwise} \\ \end{array}} \right.~~~~(8)\]

其中 \(\rho>1\) 爲常數，\(\varepsilon>0\) 爲一個小的正數。

上述精確ALM方法在內循環中要屢次更新，進行屢次奇異值分解，爲此文獻[1]提出了非精確拉格朗日乘子法(Inecact ALM, IALM)，它不須要在每次外循環開始以前要求 \(\mathop {\min }\limits_{{\bf{A}},{\bf{E}}} L({\bf{A}},{\bf{E}},{{\bf{Y}}_k},{u_k})\) 的精確解，也就是去掉了ALM方法的內循環，其更新公式變成了以下形式：

\[\begin{array}{l} {{\bf{A}}_{k + 1}} = \arg \mathop {\min }\limits_{\bf{A}} L({\bf{A}},{{\bf{E}}_{k + 1}},{{\bf{Y}}_k},{u_k}) \\ ~~~~~~~~~= {D_{\frac{1}{{{u_k}}}}}({\bf{D}} - {{\bf{E}}_{k + 1}} + \frac{{{{\bf{Y}}_k}}}{{{u_k}}}) \\ \end{array}~~~(9)\]

\[\begin{array}{l} {{\bf{E}}_{k + 1}} = \arg \mathop {\min }\limits_{\bf{E}} L({{\bf{A}}_{k + 1}},{\bf{E}},{{\bf{Y}}_k},{u_k}) \\ {~~~~~~~~~=S_{\frac{\lambda }{{{u_k}}}}}({\bf{D}} - {{\bf{A}}_{k + 1}} + \frac{{{{\bf{Y}}_k}}}{{{u_k}}}) \\ \end{array}~~~~(10)\]

上面式子中的奇異值閾值算子 \({D_{\frac{1}{{{u_k}}}}}( \cdot )\) 和軟閾值算子 \({S_{\frac{\lambda }{{{u_k}}}}}( \cdot )\) 的定義參見上一期<低秩矩陣填充|奇異值閾值算法>。

4. 低秩矩陣恢復的應用

低秩矩陣恢復技術是一個很是有研究價值和實用價值的技術，它的應用也很是普遍，好比說：

1. 視頻背景建模。
   

2. 圖像恢復(去光照、陰影等)
   

3. 圖像類別標籤淨化

4. 文本主題分析

5. 音樂詞曲分離

6. 圖像矯正與去噪
   

7. 圖像對齊

5. 仿真

ADM算法matlab代碼以下：

function [L,S] = pcp_ad(M,u,lambda,itemax,tol)
% solve PCP problem by ADM algorithm
% v1.0 2020-1-1
% function：solve the following optimization problem
%                  min  ||X||*+lambda||E||_F
%                  s.t. M = A+E

% initialize S0 and Y0 and L0
[m,n] = size(M) ;
S = zeros(m,n) ;
Y = S ;
L = S ;

% the observed matrix can contain non number
unobserved = isnan(M);
M(unobserved) = 0;

if nargin < 2
    lambda = 1 / sqrt(max(m,n));
end
if nargin < 3
    u = 10*lambda;
end
if nargin < 4
    tol = 1e-6;
end
if nargin < 5
    itemax = 1000;
end

for ii = 1:itemax
    L = sig_thre(M-S+(1/u)*Y,(1/u)) ;
    S = soft_thre(M-L+(1/u)*Y, lambda/u) ;
    Z = M-L-S ;
    Y = Y+u*Z ;
    error = norm(M-L-S,'fro')/norm(M,'fro') ;
    if (ii == 1) || (mod(ii, 10) == 0) || (error < tol)
        fprintf(1, 'iter: %04d\terr: %f\trank(L): %d\tcard(S): %d\n', ...
            ii, error, rank(L), nnz(S));
    end
    if error<tol
        break;
    end
end

數值測試代碼：

% solve PCP problem by alternating direction method
clear
clc

m = 100 ;
n = 100 ;
r = 0.05*n ;
rate = 0.05 ;
% Generating a low rank matrix
LL = randn(m,r)/sqrt(m)*randn(r,n)/sqrt(n) ;
% Generating a large sparse noise matrix (Bernoulli matrix)
SS = randi([0,1],m,n) ;
SS(SS==0) = -1 ;

% sampling
ss = SS(:) ;
index = sort(randperm(m*n,ceil(rate*n*m))) ;
ss1 = zeros(m*n,1) ;
ss1(index) = ss(index) ;
SS = reshape(ss1,m,n) ;
M = LL+SS ;

lambda = 1/sqrt(max(m,n)) ;
u = 10*lambda ;

% [L,S] = pcp_ad(M,u,lambda,1000) ;
[L,S] = RobustPCA(M,lambda,u);
% [L,S] = pcp_ad(M,u,lambda,500,1e-8);
% [L,S] = adm_lrr(M);
MM = M-L-S ;

norm(M-MM,'fro')/norm(M,'fro')
norm(M-L-S,'fro')/norm(M,'fro')
norm(L-LL,'fro')/norm(LL,'fro')
norm(S-SS,'fro')/norm(SS,'fro')

function A = soft_thre(B,T) 
A = sign(B).*max(abs(B)-T,0) ;
end

function [A] = sig_thre(B,T)

[s,v,d] = svd(B,'econ') ;
%     v(v<T) = 0 ;
%     A = s*v*d' ;
A = s*soft_thre(v,T)*d' ;
end

運行上面程序，顯示結果norm(M-L-S,'fro')/norm(M,'fro')約爲9e-7，norm(L-LL,'fro')/norm(LL,'fro')約爲1e-5。

低秩圖像恢復仿真程序：

% low rank and sparse noise image recovery
clear
clc

A = imread('C:\xxx\xxx\xxx.bmp') ;

WW = double(A) ;
a1 = double(A(:,:,1)) ;
a2 = double(A(:,:,2)) ;
a3 = double(A(:,:,3)) ;
[M,N] = size(a1);
X = zeros(M,N,3);
    
for jj = 1:3
    lambda = 1*1 / sqrt(max(M,N)); 
    u =  1*lambda;
    [ X(:,:,jj),S(:,:,jj)] = RobustPCA(WW(:,:,jj),lambda,u,1e-8,200) ;
end

figure(1)
subplot(3,1,1)
imshow(A)
title("原圖",'fontsize',12)
subplot(3,1,2)
imshow(uint8(X))
title("低秩圖",'fontsize',12)
d = S ;
d(d<20) = 255 ;
subplot(3,1,3)
imshow(uint8(d))
title("噪聲圖",'fontsize',12)

低秩圖像恢復結果以下：

從上面圖像恢復結果來看，效果還不錯。

參考文獻

[1] Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11.
[2] 史加榮, 鄭秀雲, 魏宗田, & 楊威. (2013). 低秩矩陣恢復算法綜述. 計算機應用研究, 30(6), 1601-1605.
[3] Cui, X., Huang, J., Zhang, S., & Metaxas, D. N. (2012, October). Background subtraction using low rank and group sparsity constraints. In European Conference on Computer Vision (pp. 612-625). Springer, Berlin, Heidelberg.
[4] Wright, J., Ganesh, A., Rao, S., Peng, Y., & Ma, Y. (2009). Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization. In Advances in neural information processing systems (pp. 2080-2088).
[5] Peng, Y., Ganesh, A., Wright, J., Xu, W., & Ma, Y. (2012). RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE transactions on pattern analysis and machine intelligence, 34(11), 2233-2246.

更多精彩內容請關注訂閱號優化與算法

更多精彩內容請關注微信公衆號 「優化與算法」