基于對(duì)數(shù)行列式實(shí)現(xiàn)的矩陣補(bǔ)全算法

秦國峰 彭沖 魏計(jì)鵬

摘要：傳統(tǒng)的基于低秩假設(shè)的矩陣補(bǔ)全模型常常對(duì)目標(biāo)矩陣采用核范數(shù)的約束，由于核范數(shù)對(duì)秩函數(shù)的近似不夠精確，基于核范數(shù)的低秩模型可能無法產(chǎn)生最優(yōu)的效果。為此，采用對(duì)數(shù)行列式代替核范數(shù)，提出基于最小化矩陣對(duì)數(shù)行列式的矩陣補(bǔ)全模型。研究結(jié)果表明，基于最小化對(duì)數(shù)行列式實(shí)現(xiàn)的矩陣補(bǔ)全算法能夠有效地恢復(fù)矩陣的低秩信息，能夠有效地補(bǔ)全圖像的缺失信息。

關(guān)鍵詞：矩陣補(bǔ)全;低秩結(jié)構(gòu);對(duì)數(shù)行列式;機(jī)器學(xué)習(xí)

中圖分類號(hào)：TP181

文獻(xiàn)標(biāo)志碼：A

收稿日期：2020-11-06

通信作者：彭沖，男，博士，副教授，主要研究方向?yàn)闄C(jī)器學(xué)習(xí)。E-mail： Pchong1991@163.com

在現(xiàn)實(shí)世界中，數(shù)據(jù)量變得越來越大，越來越多的數(shù)據(jù)開始采用矩陣的形式存儲(chǔ)。由于存儲(chǔ)設(shè)備的損壞，網(wǎng)絡(luò)中信息傳輸不穩(wěn)定導(dǎo)致傳輸過程中數(shù)據(jù)丟包等原因，出現(xiàn)了大量數(shù)據(jù)缺失的問題[1]，因此從非常有限的信息中估計(jì)缺失值這項(xiàng)技術(shù)變得尤為重要。矩陣補(bǔ)全技術(shù)被廣泛應(yīng)用在很多領(lǐng)域。例如，圖像恢復(fù)[2]，視頻去噪[3]和推薦系統(tǒng)[4-5]等。在近二十年中，矩陣補(bǔ)全算法得到了長足的發(fā)展，其中基于低秩的模型具有顯著的性能[6-9]。基于低秩的模型通常對(duì)目標(biāo)矩陣做出一個(gè)合理的假設(shè)，即目標(biāo)矩陣是低秩的或近似低秩的，基于此假設(shè)，研究人員提出了矩陣補(bǔ)全問題的基礎(chǔ)建模[6]。然而，矩陣補(bǔ)全基礎(chǔ)模型秩最小化問題很難求解，為了解決這個(gè)問題，現(xiàn)有算法通常使用核范數(shù)來代替秩函數(shù)。理論研究表明核范數(shù)是最貼近秩函數(shù)的凸下界[10];采用核范數(shù)來近似代替秩函數(shù)的現(xiàn)有方法有很多，例如奇異值閾值（SVT）[11]，核規(guī)范正則最小二乘法（NNLS）[12]和魯棒性主成分分析（Robust PCA）[13]。研究發(fā)現(xiàn)，當(dāng)矩陣存在較大奇異值時(shí)，核范數(shù)對(duì)矩陣秩的近似往往不夠精確，因此在實(shí)際應(yīng)用中可能會(huì)獲得次優(yōu)性能[14]。本文針對(duì)核范數(shù)對(duì)于秩函數(shù)的擬合不夠精確這一問題，提出了改進(jìn)后的矩陣補(bǔ)全算法。相對(duì)于將矩陣所有奇異值直接相加的核范數(shù)，本文采用基于矩陣對(duì)數(shù)行列式的方式近似矩陣的秩，從而有效地減小較大的奇異值對(duì)矩陣低秩性質(zhì)的影響。相比于核范數(shù)而言，對(duì)數(shù)行列式對(duì)于秩函數(shù)的近似更為精確。因此，本文通過采用對(duì)數(shù)行列式代替秩函數(shù)，得到了新的矩陣補(bǔ)全模型。

1 相關(guān)工作介紹

近些年，矩陣補(bǔ)全技術(shù)被廣泛應(yīng)用于各種領(lǐng)域，包括機(jī)器學(xué)習(xí)，計(jì)算機(jī)視覺，信號(hào)處理等。最近的研究進(jìn)展表明，基于核范數(shù)的矩陣補(bǔ)全算法對(duì)于低秩部分的近似不夠精確，針對(duì)這個(gè)問題本研究考慮采用對(duì)數(shù)行列式來解決。

1.1 對(duì)數(shù)行列式

對(duì)于任意矩陣X∈Rm×n，矩陣的對(duì)數(shù)行列式定義為

F（X）=logdet（I+XTX）（1）

顯然，logdet（I+XTX）=∑ni=1log（1+σ2i，X），其中，σ1，X>σ2，X>…>σn，X≥0，σi，X是矩陣X的第i個(gè)奇異值。

1.2 奇異值閾值算法（SVT）

SVT[11]算法是采用核范數(shù)近似代替秩函數(shù)從而實(shí)現(xiàn)矩陣補(bǔ)全的一種算法，模型為

minXX*+αX2Fs.t. PΩ（X）=PΩ（M）（2）

其中，X∈Rm×n是恢復(fù)出的矩陣;M∈Rm×n是給定的低秩的不完整的數(shù)據(jù)矩陣;α為平衡參數(shù);Ω是一組對(duì)應(yīng)于能被觀察到的項(xiàng)目的位置信息。PΩ（X）定義為：當(dāng)（i，j）∈Ω時(shí)，PΩ（X）ij=Xij; 當(dāng)（i，j）∈Ωc時(shí)，PΩ（X）ij=0。

通過使用對(duì)數(shù)行列式來對(duì)函數(shù)的低秩部分進(jìn)行擬合，這種效果要好于核范數(shù)的擬合效果，本文采用這種方式成功的解決了核范數(shù)作為秩函數(shù)近似不夠精確的問題[14]。

2 本文方法

2.1 目標(biāo)函數(shù)

考慮到核范數(shù)作為矩陣低秩部分的表現(xiàn)不夠理想，本文采用對(duì)數(shù)行列式作為對(duì)秩函數(shù)的非凸近似

min Xlogdet（I+XTX）s.t. PΩ（X）=PΩ（M）（3）

其中，I為單位陣;M∈Rm×n; X∈Rm×n是修復(fù)后的新矩陣。

2.2 模型的優(yōu)化

由于難以直接優(yōu)化目標(biāo)函數(shù)，本文通過引入輔助變量來簡化模型，從而得到新模型

min logdet（I+XTX）s.t. W=X，PΩ（X）=PΩ（M）（4）

優(yōu)化階段，本文采用了ALM（增廣拉格朗日法）[15]的優(yōu)化算法，交替迭代優(yōu)化每一個(gè)變量。給出模型的增廣拉格朗日方程式

L=F（X）+ρ2W-X+θρ2F（5）

其中，ρ 是平衡參數(shù)，θ 是拉格朗日乘子。

2.2.1 最小化X 關(guān)于X的子問題

Xt+1=argminX F（X）+ρt2Wt-X+θtρt2F（6）

其中，t為迭代次數(shù)。對(duì)于使用對(duì)數(shù)行列式的子問題，通過文獻(xiàn)[14]的證明，可知式（6）等價(jià)于求解

σ*i=argminσ*i∑ig（σ*i）+ρt2（σ*i-σi，H）2s.t. σ*i≥0，（i=1，…，n）（7）

其中，g（x）=log（1+x2）; σ*i是所求更新后的最優(yōu)解X的第i個(gè)奇異值，H=Wt+θt/Pt，σi，H≥0，（i=1，…，n）是矩陣H的第i個(gè)奇異值。對(duì)關(guān)于σi的函數(shù)f（σi）=log（1+σ2i）+ρt2（σi-σi，H）2分別求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，得

f′（σi ） = 2σi 1 + σ2i+ ρt （σi -σi，H ）（8）

f''（σi ） = ρt σ4i+ （2ρt -2）σ2i+ （2 + ρt ）（1 + σ2i）2（9）

令f′（σi）=0，得

ρtσ3i-ρtσi，Hσ2i+（ρt+2）σi-ρtσi，H=0（10）

其中，σi≥0，（i=1，…，n）; 矩陣H的SVD是U diag（{σi，H}ni=1）VT; U，V均為單位正交矩陣，且UUT=I，VVT=I，U為左奇異矩陣，V為右奇異矩陣;diag（{σi，H}ni=1）是主對(duì)角線元素為σi，H，（i=1，…，n），其余元素為0的主對(duì)角矩陣;σi，H≥0，（i=1，…，n）。

本文通過兩種情況進(jìn)行分析。

情況1：如果σi，H=0，由于σi≥0，可知f′（σi）≥0，故f（σi）是非減函數(shù)，因此，σi=0時(shí)，函數(shù)f（σi）最小。

情況2：如果σi，H>0，此時(shí)有f′（0）<0，f''（0）>0，因此至少存在一個(gè)根∈（0，σi，H）; 當(dāng)ρt>0.25時(shí)，可以看出f''（σi）>0，因此f′（σi）是有唯一全局最優(yōu)解的凸函數(shù)，此時(shí)σ*i可以通過f′（σi）=0求解得到。當(dāng)0<ρt<0.25時(shí)，可以通過以下方式進(jìn)行求解：刪掉式（10）求解中負(fù)根的集合，獲得解集Ο+，根據(jù)一階必要條件定理，可知σ*i=argminσi∈{0}∪Ο+f（σi）。

設(shè)置初始化ρt=0.000 1，且ρt隨著迭代次數(shù)不斷增大，在迭代到17次時(shí)，滿足ρt>0.25，通過這種方式，來保證當(dāng)σi，H>0時(shí)，有σ*i∈（0，σi，H）; 當(dāng)σi，H=0時(shí)，σ*i=0。

綜上，得到Xt+1的更新方式為

Xt+1=U diag（σ*1，…，σ*n）VT（11）

2.2.2 最小化W 關(guān)于W的子問題為

Wt+1=argminWρt2W-Xt+1+θtρt2F（12）

針對(duì)上述最小化問題，通過對(duì)式（10）中W求導(dǎo)，得到Wt+1的更新方式為

Wt+1=Xt+1-θtρt（13）

為了保證原圖像已知的像素值不變，加入約束條件PΩ（W）=PΩ（M），得到Wt+1的更新方式

Wt+1=PΩc（Wt+1）+PΩ（M）（14）

其中，Ωc是一組對(duì)應(yīng)于矩陣中缺失項(xiàng)的位置信息。

2.2.3 其他變量的更新 關(guān)于變量ρ與θ的更新方式為

ρt+1=ρt×k（15）

θt+1=θt+ρt+1（Wt+1-Xt+1）（16）

其中，k>1，保證每次迭代ρ都是增大的。

上述優(yōu)化算法可以歸為算法1。

算法1

輸入：需要處理的數(shù)據(jù)矩陣，以及缺失位置的下標(biāo);初始化ρ和θ。

（1） 給定最大的迭代次數(shù)matrix，迭代進(jìn)行優(yōu)化;

（2） 設(shè)置t=1，計(jì)算 max{Wt-Xt2F，Xt+1-Xt2F，Wt+1-Wt2F}≤0.001，判斷是否收斂;

（3） 根據(jù)式（11）計(jì)算Xt+1;

（4） 根據(jù)式（14）計(jì)算Wt+1;

（5） 根據(jù)式（15）計(jì)算ρt+1，根據(jù)式（16）計(jì)算θt+1;

（6） 跳轉(zhuǎn)至第二步直到收斂;

（7） Return Xt+1。

3 實(shí)驗(yàn)及結(jié)果分析

實(shí)驗(yàn)采用和SVT[11]算法相同的參數(shù)初始值，令ρ=0.000 1。當(dāng)?shù)M(jìn)行到17次時(shí)，ρ>0.25，此時(shí)滿足實(shí)驗(yàn)方法所需條件，可以認(rèn)為本文算法在前16次的迭代中，提供了滿足實(shí)驗(yàn)條件的初始化。

3.1 對(duì)比方法及實(shí)驗(yàn)數(shù)據(jù)

為保證實(shí)驗(yàn)的有效性，選取實(shí)驗(yàn)數(shù)據(jù)時(shí)，借鑒文獻(xiàn)[16]，選用對(duì)比度較為明顯的圖1～圖4，能夠更直觀的看到對(duì)比實(shí)驗(yàn)效果。實(shí)驗(yàn)設(shè)計(jì)時(shí)，分別進(jìn)行了25%，50%，75%的隨機(jī)缺失值處理，并記錄缺失值的位置。在對(duì)比實(shí)驗(yàn)部分，選用經(jīng)典的SVT（奇異值閾值算法）。

3.2 恢復(fù)圖像對(duì)比

圖1～圖4，展示了矩陣補(bǔ)全算法的對(duì)比實(shí)驗(yàn)效果圖（原數(shù)據(jù)圖片節(jié)選自文獻(xiàn)[16]，下載鏈接為https：//github.com/xueshengke/TNNR）。實(shí)驗(yàn)選用隨機(jī)添加了75%缺失值的噪聲圖像的矩陣補(bǔ)全效果圖作為展示。從圖1～圖4的對(duì)比圖展示部分，不難看出，新算法恢復(fù)后的圖像噪聲信息更少，色彩飽和度和局部圖像質(zhì)量更加接近原圖像。

3.3 PSNR值比較

實(shí)驗(yàn)中選用PSNR值來衡量恢復(fù)后的圖像質(zhì)量，PSNR值越大，表示恢復(fù)出的圖像質(zhì)量越好。對(duì)三種不同噪聲水平的圖片，分別用SVT算法和本文的算法分別進(jìn)行圖像恢復(fù)，并記錄修復(fù)后圖像的PSNR值，記錄在表1～表4中。

從表1～表4展示的實(shí)驗(yàn)結(jié)果很容易發(fā)現(xiàn)，在相同條件下，修復(fù)后的圖像在PSNR值的對(duì)比上，使用本文方法獲得的PSNR值較大，矩陣補(bǔ)全效果明顯好于經(jīng)典的SVT算法。實(shí)驗(yàn)驗(yàn)證了本文的猜想，相比于SVT算法采用核范數(shù)近似秩函數(shù)的做法，采用對(duì)數(shù)行列式的方式，能夠更好的作為秩函數(shù)的近似，充分利用了數(shù)據(jù)矩陣的低秩性質(zhì)，提高了矩陣恢復(fù)的效果。

3.4 收斂性分析

為了測試本文方法的收斂性，實(shí)驗(yàn)設(shè)計(jì)中本文針對(duì)圖1分別計(jì)算了每次迭代中的PSNR值，W-X2F，Xt+1-Xt2F，以及Wt+1-Wt2F的值，收斂圖分別如圖5～圖8所示?？芍诘?0次左右，模型已經(jīng)收斂。

4 結(jié)論

本文提出了一種新的矩陣補(bǔ)全算法的模型。在目標(biāo)函數(shù)低秩性部分改進(jìn)了核范數(shù)替代秩函數(shù)不夠精確的問題，采用對(duì)數(shù)行列式代替核范數(shù)作為秩函數(shù)的近似。在實(shí)驗(yàn)優(yōu)化部分，本文采用ALM的優(yōu)化方式對(duì)模型進(jìn)行優(yōu)化，解決了非凸問題優(yōu)化困難的技術(shù)難點(diǎn)。在所設(shè)計(jì)的圖像補(bǔ)全實(shí)驗(yàn)中，本文算法表現(xiàn)良好，證明了其有效性。

參考文獻(xiàn)

[1]王興宏.大數(shù)據(jù)應(yīng)用及新時(shí)期所面臨的挑戰(zhàn)研究[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2020，33（3）：22-27.

[2]NIKOS K. Image completion using global optimization[C]//Computer Vision and Pattern Recognition. New York， 2006： 442-452.

[3]JI H， LIU C， SHEN Z， et al. Robust video denoising using low rank matrix completion[C]// 23rd IEEE Conference on Computer Vision and Pattern Recognition， CVPR. San Francisco， 2010： 13-18.

[4]STECK H. Training and testing of recommender systems on data missing not at random[C]//Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Washington， 2010，713-722.

[5]KANG Z， PENG C， CHENG Q. Top-N recommender system via matrix completion[C]// 30th Association for the Advancement of Artificial Intelligence （AAAI） Conference on Artificial Intelligence. Phoenix， 2016，179-184.

[6]CANDES E J， RECHT B. Exact matrix completion via convex optimization[J]. Computational Math， 2009， 9： 717-772.

[7]CANDES E J， TAO T. The power of convex relaxation： Near-optimal matrix completion[J]. IEEE Transactions on Information Theory， 2010， 56（5）：2053-2080.

[8]LIU J， MUSIALSKI P， WONKA P， et al. Tensor completion for estimating missing values in visual data[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2013， 35（1）：208-220.

[9]ERIKSSON A， HENGEL A V D. Efficient computation of robust low-rank matrix approximations in the presence of missing data using the L1 Norm[C]//2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. San Francisco， 2010： 771-778.

[10] RECHT B， FAZEL M， PARRILO P A. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[J]. Siam Review， 2010， 52（3）：471-501.

[11] CAI J， CAND E J， SHEN Z. A singular value thresholding algorithm for matrix completion[J]. Siam Journal on Optimization， 2008， 20（4）：1956-1982.

[12] TOH K C， YUN S. An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[J]. Pacific Journal of Optimization， 2010， 6（3）：615-640.

[13] PENG C， CHEN Y， KANG Z， et al. Robust principal component analysis： A factorization-based approach with linear complexity[J]. Information Sciences， 2019， 513： 581-599.

[14] PENG C， KANG Z， LI H， et al. Subspace clustering using log-determinant rank approximation[C]//Proceedings of the 21th ACM SIGKDD international conference on Knowledge Discovery and Data Ming.Sydney，2015：925-934.

[15] YANG J， YUAN X. Linearized augmented lagrangian and alternating direction methods for nuclear norm minimization[J]. Mathematics of Computation， 2012， 82（281）：301-329.

[16] HU Y， ZHANG D B， YE J P， et al. Fast and accurate matrix completion via truncated nuclear norm regularization[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence， 2013， 35（9）：2117-2130.

Matrix Completion Algorithm Based On Log-determinant

QIN Guo-feng， PENG Chong， WEI Ji-peng

（College of Computer Science and Technology， Qingdao University， Qingdao 266071， China）

Abstract：

The nuclear norm constraint is mainly used for the target matrix in traditional matrix completion models based on low-rank assumptions. However， the nuclear norm is not accurate enough to approximate the rank function. Thus it may lead to degraded performance in low-rank recovery problem. To resolve the drawback of the nuclear norm in low-rank recovery， the log-determinant is used to replace nuclear norm and the matrix completion model based on the log-determinant minimization is proposed. Experimental results show that proposed method can effectively recover the low-rank structure of the target matrix and complete the missing entries of corrupted images in image completion application.

Keywords：

matrix completion; low-rank structure; log-determinant; machine learning