低秩矩陣填充|奇異值閾值算法

---

斜風細雨做小寒，淡煙疏柳媚晴灘。入淮清洛漸漫漫。
雪沫乳花浮午盞，蓼茸蒿筍試春盤。人間有味是清歡。
---- 蘇軾

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

低秩矩陣恢復是稀疏向量恢復的拓展，兩者具備不少能夠類比的性質。首先，稀疏是相對於向量而言，稀疏性體如今待恢復向量中非零元素的數量遠小於向量長度；而低秩是相對於矩陣而言，低秩體如今矩陣的秩遠小於矩陣的實際尺寸。其次，稀疏向量恢復問題能夠轉化爲基於 \(\ell _1\) 範數的凸優化問題，由於 \(\ell _1\) 範數是 \(\ell _0\) 的最佳凸包絡；而矩陣的核範數在必定條件下也是矩陣秩的最佳凸近似，所以，也能夠利用這一性質將低秩矩陣恢復問題鬆弛爲一個凸問題來求解。本文重點介紹低秩矩陣恢復的一些經常使用算法，並給出了奇異值閾值（STV）算法用於低秩矩陣恢復的仿真分析。

1. 低秩矩陣恢復問題背景

隨着5G技術和物聯網（IoT）技術的發展，做爲5G三大應用場景之一的大規模機器類通訊（eMTC）將使「萬物互聯」成爲現實。大規模機器類通訊的部署勢必致使海量數據的產生，而實時、精確地處理這些大規模數據對傳統的信息處理技術帶來了挑戰。因爲數據維度的上升，在信號處理、圖像處理、計算機視覺、機器學習和數據挖掘等領域，所需處理的數據量呈幾何級數增加，面對海量的數據，不少傳統的信號處理技術已不堪重負，沒法適應於將來爆炸式增加的數據環境，所以急需開發更加先進的信號處理技術。

隨着數據維度的增長，高維數據之間每每存在較多的相關性和冗餘度，並且數據自己信息量的增加比數據維度的增加一般要慢得多。例如，視頻信號的可壓縮空間比單幅圖像的可壓縮空間大得多。一幅天然圖片通過小波變換以後，只有少許的係數在數值大小上是相對顯著的。若是將一幅圖片當成一個像素矩陣，對其進行奇異值分解後，其前10%奇異值包含的信息量每每佔了整幅圖像的90%。這些實例都說明在高維數據中，每每存在不一樣程度的相關性，利用這些相關性能夠大幅減少數據的存儲空間和處理複雜度。

此外，現實生活中的大規模數據經常會有部分數據缺失、數據偏差、損壞等問題，這將進一步加大數據處理和分析難度。這種如今在實際生活中很常見，例如在人臉識別中,訓練集中的或是待識別的人臉圖像每每含有陰影、高光、遮擋、變形等;在運動恢復結構(Structure from motion,SFM)問題中,進行特徵提取和特徵匹配時每每存在大的匹配偏差. 這些因素的存在使得不少傳統的分析和處理方式失效, 須要新的理論和實用的算法爲相關的應用提供理論支撐和有力的求解工具.正確並高效地從不完整的、帶有損毀的數據中 恢復和利用它們, 對現代大規模數據的分析與處理相當重要。

假定原始數據矩陣是低秩的，可是矩陣中含有不少未知的元素。從一個不完整的矩陣中恢復出一個完整的低秩陣，即是低秩矩陣填充問題。例如，著名的Netflix問題即是一個典型的低秩矩陣填充問題。Netflix是美國的一家影片租賃公司，其推薦系（recommendation system）要從用戶僅有的對少數的電影打分中爲用戶推薦影片，若是這種推薦越符合用戶的喜愛，也就越能提升該公司租賃電影的業務，爲此，該公司設立了百萬美圓的獎金用於懸賞可以最好地提升該公司推薦系統質量的解決方法。這個問題能夠用矩陣填充來進行建假設矩陣的每一行分別表明同一用戶對不一樣電影的打分, 每一列表明不一樣用戶對同一電影的打分，用戶數量巨大, 電影數目巨大，所以這個矩陣的維度十分大。因爲用戶所打分的電影有限，這個矩陣中只有很小一部分的元素值已知，並且可能含有噪聲或偏差，那麼Netflix問題就是如何從這個不完整的矩陣中推測其中未知元素的值。矩陣填充的越準確，爲用戶推薦的電影也就越符合用戶的喜愛。

2. 低秩矩陣補全

矩陣補全是推薦系統、計算機視覺、圖像處理等領域常常遇到的問題，具備很高的研究價值。矩陣補全問題旨在經過真實未知矩陣 \(\bf M\) 的部分觀測 \({\bf{M}}_{i,j},~~(i,j)\in\Omega\) 來估計 \(\bf M\) 中未觀測到的其餘元素。若是不加其餘約束，這樣的問題是徹底不可解的。可是，若是數據矩陣具備一些特殊的性質，這些特殊的性質將使得矩陣補全問題成爲可能。低秩性就是這樣的性質。若是事先知道矩陣 \(\bf M\) 是低秩的，那麼矩陣補全問題能夠形式爲爲以下優化問題：

\[\begin{array}{l} \mathop {\min }\limits_{\bf{X}}~~~ rank({\bf{X}}) \\ s.t.~~~{{\bf{X}}_{i,j}} = {{\bf{M}}_{i,j}},~~(i,j) \in \Omega \\ \end{array}~~~~~~(1)\]

其中 \(\Omega\) 爲觀測樣本下標的集合，\(\bf X\) 爲優化變量，\(\bf M\) 爲真實的未知矩陣。定義投影算子 \(P_{\Omega}\):
\[{\left[ {{P_\Omega }({\bf{X}})} \right]_{i,j}} = \left\{ {\begin{array}{*{20}{c}} {{{\bf{X}}_{i,j}}~~~(i,j) \in \Omega } \\ {0~~~~~{\rm{otherwise}}} \\ \end{array}} \right.~~~~(2)\]

從而(1)式能夠簡潔地表述爲：
\[\begin{array}{l} \mathop {\min }\limits_{\bf{X}} ~~~rank({\bf{X}}) \\ s.t.~~~{P_\Omega }({\bf{X}}) = {P_\Omega }({\bf{M}}) \\ \end{array}~~~~~~~(3)\]

顯然，上述問題是一個NP-hard問題，且其複雜度隨矩陣維度的增城呈平方倍指數關係增加。爲了求解問題(1)，有學者對其進行凸鬆弛，轉化爲一個凸優化問題：
\[\begin{array}{l} \mathop {\min }\limits_{\bf{X}}~~~ {\left\| {\bf{X}} \right\|_*} \\ s.t.~~{{\bf{X}}_{i,j}} = {{\bf{A}}_{i,j}},~~(i,j) \in \Omega \\ \end{array}~~~~~(4)\]

文獻[2]中指出：問題（4）的解在知足強非相干性條件下，很大機率等價於問題(1)的解。在某種意義上，這是NP-難秩最小化問題的緊凸鬆弛。由於核範數球是譜範數爲1的rank one矩陣集合的凸包。具體定理及證實參見文獻[1]。
對於(4)中的優化問題，當採樣點數知足 \(m \ge Cn^{6/5}r\log n\) 時，\(\bf{M}\) 能以很大機率求得精確解，其中 \(r\) 爲矩陣 \(\bf{M}\) 的秩，\(C\) 爲一個正常數。

求解(4)式可使用一些凸鬆弛方法，好比半定規劃法，然而半定規劃法一般使用內點法來求解，其求解上述問題的算法複雜度爲 \(O(p{(m + n)^3} + {p^2}{(m + n)^2} + {p^3})\)。能夠經過matlab軟件包CVX等來求解。

3. 奇異值閾值算法（SVT）

2009年，Cai Jian-feng 和 Candes（e上面有個四聲符號）等人提出了求解核範數最小化問題的奇異值閾值算法（Singular Value Thresholding, SVT）。該算法能夠類比爲求解向量 \(\ell _0\) 範數最小化的軟閾值迭代算法。 SVT算法先將最優化問題(4)正則化，即有：
\[\begin{array}{l} \mathop {\min }\limits_{\bf{X}} ~~~~\tau {\left\| {\bf{X}} \right\|_*} + \frac{1}{2}\left\| {\bf{X}} \right\|_F^2 \\ s.t.~~~{_\Omega }({\bf{X}}) = {P_\Omega }({\bf{M}}) \\ \end{array}~~~~~(5)\]

其中，\(\tau > 0\)。當 $\tau \to + \infty $ 時，上述最優化問題的最優解收斂到(4)式的最優解。構造最優化問題(5)的拉格朗日函數：
\[L({\bf{X}},{\bf{Y}}) = {\left\| {\bf{X}} \right\|_*} + \frac{1}{2}\left\| {\bf{X}} \right\|_F^2 + \left\langle {{\bf{Y}},{P_\Omega }({\bf{M}} - {\bf{X}})} \right\rangle ~~~~(6)\]

其中，拉格朗日乘子 \({\bf{Y}}{ \in {\mathbb{R}}^{m \times n}}\)。若是 \(({{\bf{X}}^*},{{\bf{Y}}^*})\) 爲優化問題(4)的原-對偶問題的最優解，則有：
\[\mathop {\sup }\limits_{\bf{Y}} \mathop {\inf }\limits_{\bf{X}} L({\bf{X}},{\bf{Y}}) = \mathop {\inf }\limits_{\bf{X}} \mathop {\sup }\limits_{\bf{Y}} L({\bf{X}},{\bf{Y}}) = L({{\bf{X}}^*},{{\bf{Y}}^*})~~~(7)\]

SVT算法使用交替迭代方法求解優化問題(4)，其迭代格式能夠簡單表述以下：
\[\left\{ {\begin{array}{*{20}{c}} {{{\bf{X}}^k} = {D_\tau }({{\bf{Y}}^{k - 1}} )} \\ {{{\bf{Y}}^k} = {{\bf{Y}}^{k - 1}} + {\delta _k}{P_\Omega }({\bf{M}} - {{\bf{X}}^k})} \\ \end{array}} \right. ~~~~(8)\]

其中 \({D_\tau }({\bf{W}})\) 爲奇異值閾值軟閾值操做相似，不過這裏的對象是矩陣。能夠具體描述爲：
\[{\bf{W}}' = {D_\tau }({\bf{W}}) = \left\{ {\begin{array}{*{20}{c}} {(1)~~~~~~~~~~~~~~[{\bf{U}},{\bf{S}},{\bf{V}}] = {\mathop{\rm svd}\nolimits} (W)} \\ {(2)~~{\bf{S}} = {\mathop{\rm sgn}} ({\bf{S}}).*\max ({\mathop{\rm abs}\nolimits} ({\bf{S}}) - \tau ,0)} \\ {(3)~~~~~~~~~~~~~~~~~{\bf{W}}' = {\bf{U}}*{\bf{S}}*{{\bf{V}}^T}} \\ \end{array}} \right.~~~~(9)\]

其中(8)式中的第一步是這樣來的，初始化 \({{\bf{Y}}^0} = 0\)，當 \({{\bf{Y}}^{k-1}}\) 固定時：
\[\begin{array}{l} {{\bf{X}}^k} = \mathop {\arg \min }\limits_{\bf{X}} \tau {\left\| {\bf{X}} \right\|_*} + \frac{1}{2}\left\| {\bf{X}} \right\|_F^2 - \left\langle {{{\bf{Y}}^{k - 1}},{P_\Omega }({\bf{X}})} \right\rangle \\ ~~~~~= \mathop {\arg \min }\limits_{\bf{X}} \tau {\left\| {\bf{X}} \right\|_*} + \frac{1}{2}\left\| {\bf{X}} \right\|_F^2 - \left\langle {{\bf{X}},{P_\Omega }({{\bf{Y}}^{k - 1}})} \right\rangle \\ ~~~~~= {D_\tau }({{\bf{Y}}^{k - 1}}) \\ \end{array}~~~~(10)\]

(8)式中的第(2)式，當 \({{\bf{X}}^k}\) 給定時，用梯度降低來更新 \(\bf{Y}\)。

4. 仿真

SVT算法的matlab代碼以下：

function [ X,iterations ] = SVT(M,P,T,delta,itermax,tol)
%Single value thresholding algorithm，SVT
% function：solve the following optimization problem
%                  min  ||X||*
%               s.t. P(X-M) = 0
% X: recovered matrix
% M: observed matrix
% T: single value threshold
% delta: step size
% output：X,iterations

% initialization
Y = zeros(size(M));
iterations = 0;

if nargin < 3
    T =  sqrt(n1*n2);
end
if nargin < 4
    delta = 1;
end
if nargin < 5
    itermax = 200 ;
end
if nargin < 6
    tol = 1e-7;
end

for ii = 1:itermax
    % singular value threshold operation
    [U, S, V] = svd(Y, 'econ') ;
    S = sign(S) .* max(abs(S) - T, 0) ;
    X = U*S*V' ;
    % update auxiliary matrix Y
    Y = Y + delta* P.* (M-X);
    Y = P.*Y ;
    % computer error
    error= norm( P.* (M-X),'fro' )/norm( P.* M,'fro' );
    if error<tol
        break;
    end
    % update iterations
    iterations = ii ;
end
end

數值驗證代碼爲以下：

% Numerical verification for SVT algorithm
clear
clc

r = 2 ; % rank(M) = 2 ;
n1 = 200 ; % row length of M
n2 = 300 ; % column length of M
sample_rate = 0.5 ; 
% real
% M = rand(n1,r)*rand(r,n2) ;
% complex
M = (rand(n1,r)+1i*rand(n1,r))*(rand(r,n2)+1i*rand(r,n2))/2 ;

% construct index matrix
P = zeros(n1*n2,1) ;  
MM = reshape(M,n1*n2,1) ;
pos = sort(randperm(n1*n2,n1*n2*sample_rate))' ;
P(pos) = MM(pos) ;
index1 = find(P) ;
P(index1) = 1 ;
P = reshape(P,n1,n2) ;

% set threshold & step size
T = sqrt(n1*n2); 
delta = 2 ;

[ X,iterations ] = SVT(M,P,T,delta) ;

% norm(P.*(X-M),'fro')
norm(P.*(X-M),'fro')/norm(P.*M,'fro')

% norm(X-M,'fro')
% norm(X-M,'fro')/norm(M,'fro')

上述測試迭代200次中止時，norm( P.* (M-X),'fro' )/norm( P.* M,'fro' )偏差大約爲1e-7。

% low rank image completion by SVT algorithm
clear
clc

A = imread('C:\Users\pc\Desktop\CS_Rectest\MC\ce_svt\channel.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);
% sampling
pos = sort(randperm(M*N,M*N*0.5))' ;
    
for jj = 1:3
    P = zeros(M*N,1) ;
    MM = reshape(double(WW(:,:,jj)),M*N,1) ;
    P(pos) = MM(pos) ;
    index1 = find(P) ;
    P(index1) = 1 ;
%     P(index1) = rand(size(index1,1),1) ;
    P = reshape(P,M,N) ; 
    T = 50000 ;
    delta_t = 1 ; 
    [ X(:,:,jj),iterations ] = SVT(WW(:,:,jj),P,T,delta_t) ;
%     [ X(:,:,jj)] = pcp_ad(WW(:,:,jj)) ;
end

DD = P.*WW(:,:,1) ;
DD1 = P.*WW(:,:,2) ;
DD2 = P.*WW(:,:,3) ;
ff(:,:,1) = DD;
ff(:,:,2) = DD1;
ff(:,:,3) = DD2;

figure(1)
subplot(1,3,1)
imshow(A)
title("原圖")

subplot(1,3,2)
imshow(uint8(ff))
title("採樣後的圖")
subplot(1,3,3)
imshow(uint8(X))
title("恢復的圖")

仿真效果以下：

能夠看出，在採樣率爲0.5時，圖像恢復效果很不錯。

因爲時間關係，下次再介紹低秩與稀疏矩陣恢復算法，敬請期待！

5. 參考文獻

[1] Cai, J. F., Candès, E. J., & Shen, Z. (2010). A singular value thresholding algorithm for matrix completion. SIAM Journal on optimization, 20(4), 1956-1982.
[2] Davenport, M. A., & Romberg, J. (2016). An overview of low-rank matrix recovery from incomplete observations. IEEE Journal of Selected Topics in Signal Processing, 10(4), 608-622.​

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