馮雅莉 孫為軍

(廣東工業(yè)大學(xué) 廣東 廣州 510006)

0 引 言

近年來，矩陣的使用越來越廣泛和越來越受重視，而矩陣分析方法的重要性也顯得越來越突出，例如矩陣分解、非負矩陣分解、低秩矩陣近似、矩陣填充等。而本文主要是關(guān)注其中的矩陣填充問題，在實際問題中，矩陣填充問題主要表現(xiàn)在圖像和視頻處理、推薦系統(tǒng)、文本分析和機器學(xué)習(xí)[1-2]等方面。為了方便對數(shù)據(jù)進行處理，一般我們會把數(shù)據(jù)轉(zhuǎn)化成矩陣的表示形式，但這些數(shù)據(jù)轉(zhuǎn)為矩陣后，經(jīng)常會面臨缺失，損毀和噪聲污染等問題，所以如何恢復(fù)數(shù)據(jù)和對數(shù)據(jù)進行去噪是我們需要解決和研究的問題，而矩陣填充技術(shù)就是解決這類問題的其中一種有效的方法。矩陣填充技術(shù)在圖像修復(fù)[3]、核磁共振圖像分割[4]、EEG信號處理、推薦系統(tǒng)[5]等領(lǐng)域中被廣泛應(yīng)用。

在近二十年里，矩陣填充技術(shù)越來越受關(guān)注。它的主要目的是通過矩陣中已知的元素來恢復(fù)未知的元素，通常利用矩陣中行與行或者列與列之間的線性關(guān)系對未知元素進行估計和恢復(fù)，所以需要被填充的矩陣是低秩的，這是矩陣填充的重要前提之一。近年來，在矩陣填充領(lǐng)域中涌現(xiàn)了許多有效的算法，例如近鄰梯度下降算法[6](proximal gradient descent,PGD)、分塊坐標(biāo)下降算法[7]、矩陣空間Bregman迭代算法[8]、交替方向乘子法[9](alternating direction method of multipliers,ADMM),還有一些隨機優(yōu)化方法，例如隨機近鄰梯度下降算法(stochastic proximal gradient descent,SPGD)。

在解決實際的矩陣填充的問題中，對矩陣進行降維能夠大大減少矩陣填充在計算中的成本，常用的方法有主成分分析(principal component analysis,PCA)和奇異值分解(singular value decomposition,SVD)。大多矩陣填充方法都是利用傳統(tǒng)的SVD方法對原始矩陣進行降維的，但是傳統(tǒng)的SVD的操作計算成本很高，所以不適用于大規(guī)模大尺寸的矩陣中。本文針對這個問題，提出了一種快速的隨機投影方法，同時利用人工仿真數(shù)據(jù)和實際圖片像素數(shù)據(jù)進行了多個實驗，驗證了本算法的可行性。

1 基礎(chǔ)知識

對于秩為r的矩陣X∈Rm×n，它的奇異值分解為：X=UΣVT，其中U∈Rm×r和V∈Rn×r滿足UTU=VTV=I。(I為r階單位矩陣)，Σ=diag(σ1,σ2,…,σr)且其對角線上的元素滿足σ1≥σ2≥…≥σr>0，σi表示X的第i個奇異值。

X的Frobenious范數(shù)滿足：

式中：〈X,X〉表示X與X的內(nèi)積，trace(·)表示矩陣的跡算子。

2 快速矩陣填充方法

2.1 矩陣填充模型

秩是矩陣的一個重要的全局信息[9]，它能體現(xiàn)出矩陣的信息冗余度，而矩陣能夠被充分地填充的前提條件是矩陣的信息是冗余的，其冗余度越低，填充難度就越高。所以矩陣填充的問題一般是最小化矩陣的秩的問題，其模型如下：

(1)

s.t.XΩ=MΩ

式中：X,M∈Rm×n，M可觀察的元素的索引存放在Ω集合中，而索引不在Ω集合中的元素均為未知的。由于函數(shù)rank(X)是非凸的，在實際優(yōu)化中求解非常困難。文獻[10]提出對rank(X)的最小化問題可以轉(zhuǎn)化成對X的核范數(shù)最小化問題，模型如下：

(2)

s.t.XΩ=MΩ

在過去使用核范數(shù)優(yōu)化模型的矩陣填充方法很多，例如奇異值閾值算法(Singular Value Thresholding,SVT)[11], 加速近端梯度算法(Accelerated Proximal Gradient,APG)[12]，增廣拉格朗日乘子法(Augmented Lagrange Multipliers,ALM)。

本文引入一個新的變量Z，并把式(2)改寫成：

(3)

rank(Z)≤rank(X)

2.2 隨機快速矩陣填充方法

解式(3)，一般采用SVD的方法，但是采用傳統(tǒng)的SVD的方法，計算的時間和計算復(fù)雜度都非常大和高，這不利于對大尺寸的矩陣進行填充和處理，所以本文提出了一種快速的隨機投影方法。

假設(shè)一個矩陣X∈Rm×n，Ω和ΩC分別表示X矩陣中的已知和未知元素的索引集合。為了解式(3)，我們可以采用交替最小二乘法對Z進行更新。

本文方法分為以下幾個部分：

第一部分是對數(shù)據(jù)進行預(yù)處理，初始化帶未知元素的矩陣X，把矩陣X的未知元素用獨立同分布的均值為0方差為1的高斯隨機分布賦值，如下所示：

第二部分，如果對X進行傳統(tǒng)的SVD操作，則計算復(fù)雜度為Ο(min(m2n,mn2)),隨著矩陣的尺寸增大，計算的復(fù)雜度呈指數(shù)增長，所以需要對矩陣X進行降維，方法如下：

C=X×H

(4)

式中：H∈Rn×r，且H是由獨立同分布的高斯隨機分布組成，至于r是如何選擇的可以參考文獻[13]。

第三部分，就是求出X的近似左奇異值矩陣：

U=C×V×Σ-1

(5)

式中：V是利用對CTC求特征向量獲得，Σ是通過對CTC求特征值開平方根獲得。

第四部分，對X進行矩陣填充，或者是低秩近似。因為U在正常情況下是C的一個正交矩陣，所以UUT=I，而C是X的一個近似矩陣，所以我們以如下方式更新X：

X≈U×UT×X

(6)

第五部分是更新停止的條件的設(shè)定，本文采用均方差(mean squared error,MSE)來作為停止條件之一，定義如下：

(7)

當(dāng)MSE小于或者等于我們設(shè)定的閾值時，或者迭代次數(shù)達到了我們設(shè)定的最大迭代次數(shù)時，迭代停止。

本文的隨機投影算法的詳細步驟如算法1所示。

算法1基于隨機投影的矩陣填充算法

輸入：帶有未知元素的X矩陣(X∈Rm×n)，估計的秩r，閾值t，最大迭代次數(shù)M

輸出：填充后的矩陣X

初始化：X中的未知元素都用獨立同分布的均值為0方差為1的高斯隨機分布賦值

循環(huán)： 迭代次數(shù)n=1:M

步驟1C=X×H，其中H∈Rn×r

i.i.d.N(0,1)

步驟2 [V,d]=eig(CTC)，其中V是C的右奇異矩陣，d是CTC的特征值組

步驟3U=C×V×Σ-1，其中

Σ=diag(diag(d).^0.5)

步驟4X=U×UT×X

當(dāng)MSE≤t或者n=M時，迭代結(jié)束。

2.3 復(fù)雜度分析

算法1中，步驟1的計算復(fù)雜度為Ο(mnr),步驟2計算復(fù)雜度為Ο(nr+r2)，步驟3 的計算復(fù)雜度為Ο(2nr2)，步驟4的計算復(fù)雜度為Ο(m2r+m2n)，所以本文算法迭代一次的計算復(fù)雜度大約為Ο((n+r)m2+mnr+2nr2+nr+r2)。

3 實驗分析

將本文算法與SVT[11]、APG[12]和ALM[14]三個經(jīng)典的算法進行比較，由于精確增廣拉格朗日乘子法(Exact Augmented Lagrange Multipliers，EALM)的計算成本較高，所以本文僅對ALM算法中的非精確增廣拉格朗日乘子法(Inexact Augmented Lagrange Multipliers，IALM)作比較。實驗中，我們假設(shè)需要填充的矩陣M∈Rm×n，其秩為r=(r1,r2)，實驗的收斂條件為[15]：

(8)

或者是迭代次數(shù)達到設(shè)置的最大迭代次數(shù)。式中Vi表示我們需要迭代更新的第i個變量，ε表示收斂或者迭代停止的閾值，一般設(shè)定為10-5，本文中的最大迭代次數(shù)設(shè)置為100。

3.1 仿真數(shù)據(jù)實驗

本節(jié)包含兩個部分：第一部分是需要被恢復(fù)的矩陣的秩不同的情況下，四個算法的恢復(fù)情況；第二部分是在不同采樣率下四個算法的恢復(fù)情況。對矩陣恢復(fù)有重要影響的兩個參數(shù)為恢復(fù)矩陣的秩和采樣率。我們用均方差(MSE)和運行時間來作為對矩陣恢復(fù)效果的衡量標(biāo)準(zhǔn)。

第一部分，假設(shè)矩陣M∈R100×100，在其秩從2遞增至20間設(shè)計了10組實驗，采樣率p=0.5。實驗結(jié)果如圖1所示。由圖1(a)可知，當(dāng)矩陣的秩在2～20之間時，F(xiàn)RPMC方法的對矩陣的恢復(fù)的效果僅次于ALM算法，但圖1(b)中，F(xiàn)RPMC的恢復(fù)的速度比ALM快了一倍。在仿真實驗中雖然APG和IALM的運行速度相對比較快，但是它們矩陣恢復(fù)的相對誤差較高，即恢復(fù)的精度相對較低。

第二部分，假設(shè)矩陣M∈R100×100，矩陣的秩為5，在采樣率分別為0.4、0.6、0.8的條件下，用FRPMC、APG、IALM分別對矩陣進行恢復(fù)，結(jié)果如圖2所示。由圖2(a)可知，F(xiàn)RPMC方法的精度相對另外的兩種算法要高，而且隨著采樣率的增加，恢復(fù)的精度越高，恢復(fù)效果越好，即FRPMC在高采樣率的情況下，恢復(fù)效果比較顯著。

(a)

(b)圖1 不同算法在矩陣的秩不同的情況下矩陣恢復(fù)的結(jié)果

(a)

3.2 真實數(shù)據(jù)實驗

本節(jié)主要利用對實際的缺失的黑白圖片的恢復(fù)來說明FRPMC能夠?qū)崿F(xiàn)圖像恢復(fù)的功能。因為黑白圖片的像素相當(dāng)于一個二維的矩陣。實驗將FRPMC算法與SVT、APG、EALM、IALM進行比較，并用峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)[17]作為衡量圖像恢復(fù)的效果的完整性的標(biāo)準(zhǔn)，定義如下：

(9)

圖3中的6幅圖是分別在不同隨機采樣下生成的，采樣率從左至右分別是0.3、0.4、0.5、0.6、0.7、0.8。圖4的6幅圖則是圖3中相對應(yīng)的用FRPMC算法恢復(fù)后的圖像。圖4(a)是采樣率為0.3的情況下恢復(fù)后的圖像，從圖中可以清晰地看到房子的輪廓窗、房子在水中的倒影，說明FPRMC方法能夠有效恢復(fù)帶有丟失數(shù)據(jù)的圖像。

(a) (b) (c) (d) (e) (f)圖4 經(jīng)過FRPMC處理后與圖3相對應(yīng)的恢復(fù)后的圖像

接下來，隨機抽取來自Berkeley Segmentation數(shù)據(jù)集中的50幅圖來進行實驗。首先把這50幅彩圖進行灰度化，使其變成黑白圖片，再進行實驗。

對這50幅處理后的黑白圖片進行隨機采樣，采樣率分別是0.3、0.4、0.5、0.6、0.7、0.8，分別用FRPMC、SVT、APG和IALM四種方法進行實驗。實驗假設(shè)最大迭代次數(shù)為100。實驗結(jié)果如表1和表2所示。

表1 在不同采樣率下四個算法的PSNR值

表2 在不同采樣率下四個算法的運行時間

由表1可知，在圖像的尺寸為512×512的情況下，即FPRMC算法對黑白圖像恢復(fù)的PSNR是最高的，即效果最好的。由表2可知，F(xiàn)RPMC算法在實際圖像恢復(fù)中的速度也是最快的。

圖5為在采樣率為0.7的實驗中隨機選取的5幅圖的實驗情況。表3是圖5中對應(yīng)的5幅圖的恢復(fù)結(jié)果，分別是PSNR值和各個算法對應(yīng)的運行時間。

(a) 原始圖片 (b) 加噪后圖片(P=0.7) (c) FRPMC算法 (d) SVT算法 (e) APG算法 (f) IALM算法圖5 隨機選取5幅圖的實驗情況

表3 圖5中5幅圖相對應(yīng)的PSNR值和對應(yīng)算法的運行時間

4 結(jié) 語

傳統(tǒng)的矩陣填充方法很難避免標(biāo)準(zhǔn)的SVD的計算操作，而隨著矩陣維度的增加，SVD操作的計算復(fù)雜度呈指數(shù)增加，不適用于大規(guī)模、高維和高階的矩陣的處理。如何設(shè)計可擴展的快速的算法是矩陣填充算法研究的核心。為了解決這個問題，本文提出了一種快速的隨機投影方法，相對比于一般傳統(tǒng)的算法。FRPMC算法主要是利用隨機投影的方法對經(jīng)過隨機初試化的矩陣(帶有缺失值)進行降維，然后再對其進行特征值和特征向量的計算，利用SVD的計算模型來設(shè)計了一個矩陣的低秩近似的模型，對矩陣進行恢復(fù)，大大減少了計算的成本。

在今后的矩陣填充算法研究中，還將關(guān)注以下幾個方面：(1) 在確保精度的情況下，提高矩陣填充的速度；(2) 設(shè)計適合不同應(yīng)用場景的矩陣填充的模型；(3) 將矩陣填充擴展到高維數(shù)組中，例如張量填充問題上；(4) 研究可以利用已知元素來確定或者估計矩陣的秩的方法。