---

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

更多精彩內容請關注微信公眾號 “優化與演算法”

低秩矩陣恢復是稀疏向量恢復的拓展，二者具有很多可以類比的性質。首先，稀疏是相對於向量而言，稀疏性體現在待恢復向量中非零元素的數量遠小於向量長度；而低秩是相對於矩陣而言，低秩體現在矩陣的秩遠小於矩陣的實際尺寸。其次，稀疏向量恢復問題可以轉化為基於 \(\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.​

更多精彩內容請關注微信公眾號 “優化與演算法”

相關推薦