基于核范數凸優化的微陣列缺失點重建

孟繁馳，李書琴，蔡 騁

(西北農林科技大學信息工程學院，陜西楊凌712100)

0 引言

基因微陣列技術在生物實驗中已經得到了廣泛的應用，從微陣列得到的基因表達值，通常以大規模矩陣的形式呈現［1］。由于微陣列實驗中各種因素的影響，如雜交失敗、圖像的噪聲、污染等，從微陣列得到的基因矩陣通常含有缺失點［2］。然而，基因下游的實驗，如常用的分類、聚簇方法等，通常需要完整的矩陣作為輸入。鑒于此，不同類型的基因微陣列缺失點的重建方法陸續產生。目前常用的基因微陣列缺失點的重建方法包括奇異值分解SVD［3］、K最近鄰KNN［3］、貝葉斯主成分分析BPCA［4］、最小二乘法LSimpute［5］、局部最小二乘法 LLSimpute［6］等。在以上各種方法中，LLSimpute被實驗證明可以在各種類型的基因數據集上取得較高的重建精確度［1，6-8］。

1 LLSimpute方法

LLSimpute是Hyunsoo Kim等人在2005年提出的方法［6］。該方法利用基因與基因之間的線性關系來重建目標基因中的缺失點。以一個m×n的基因微陣列矩陣為例，簡要介紹一下LLSimpute的流程:

(1)從其他m-1個基因行中找出與目標基因 (即當前需要重建的含缺失點的一行基因)g1距離最近的k個基因，組成列向量 (g1，gs1，…，gsk)T。

(2)找出目標基因中的缺失點g1×q和非缺失點g1×(n-q)。其中q表示該行中有q個缺失點。

(3)將缺失的點g1×q組成的向量表示為α，將非缺失的點g1×(n-q)組成的向量表示為wT。從gs1到gsk中找出與g1中的缺失點對應的列上的基因，組成矩陣 (或向量)Bk×q，從gs1到gsk中找出與g1中的非缺失點對應的列上的基因，組成矩陣Ak×(n-q)。假設目標基因上有兩個缺失點，則形成的矩陣如式 (1)所示

(4)在以上矩陣中，wT、B、A是已知的，α表示目標基因中待重建的缺失點組成的向量。由wT和A之間的線性關系，求系數向量x，可以看作一個最小二乘法問題 (2)

其中(AT)是AT的偽逆陣 (Pseudoinverse)。從而

實驗表明，在k選擇合適的時候，LLSimpute可以獲得非常好的重建效果。［6］中，為了找出最優的k值，作者采取了如下的方法:首先，用每行的均值將該行中的缺失點填充，得到一個完整的矩陣。然后在該完整的矩陣當中人為地去掉一些點，構成一個人造的含缺失點的矩陣。用不同數量的k個鄰居，使用LLSimpute去重建，得到重建出的缺失點的數值。由于這些缺失點是人為制造的，所以它們的真實值是已知的，這樣，不同k對應的重建錯誤率是可以計算的。最后，取對應重建錯誤率最低的k，在基因微陣列矩陣中作為選取的鄰居數。

2 矩陣填充

矩陣填充 (matrix completion，MC)是 Candès等人在2009年提出的方法 ［9］。在 ［9］中，作者證明了一個含有缺失點的矩陣，在低秩的情況下，如果其中所含的非缺失點足夠多，則完整的矩陣是可以在一定的概率下恢復出來的。因此，矩陣填充可以認為是一個 (4)所示的優化問題

其中M是觀測到的，含有缺失點的矩陣，Ω是非缺失點的元素對應的下標組成的集合，X是待重建的完整矩陣。但 (4)是一個NP難問題。然而矩陣X的秩與矩陣X的奇異值中非零元素的個數是一樣的，所以［9］中，X的秩被X的核范數 (nuclear norm)所代替，上式變為 (5)所示的優化問題

其中‖·‖*表示核范數，等于所有奇異值的和。式(4)和式 (5)的區別在于，式 (4)中X的秩對應的是其奇異值組成的向量的0范數，而式 (5)中X的核范數對應其奇異值組成的向量的1范數。式 (5)所對應的優化問題是凸的，具有唯一的解。

矩陣填充已經被成功用在Netflix視頻推薦系統［10］、視頻去噪上［11］。林宙辰等人提出的非精確增廣拉格朗日乘子(inexact augmented lagrange multiplier，IALM)［12］是解決式(5)所示的優化問題的一個有效方法，在重建精確度和運行速度上均占有很大優勢。

IALM采用拉格朗日乘子解決 (5)所示的核范數凸優化問題，這里，矩陣填充問題被描述為式 (6)所示的一個優化問題

式中:A——待重建的矩陣，D——觀測到的有缺失點的矩陣，πΩ——一個的線性變換。IALM在每次迭代中求奇異值分解，每次取較大的幾個奇異值來重建矩陣A和E，從而使A趨向于低秩。詳細的迭代過程，可參考［12］。本文中，將使用IALM進行矩陣填充，得到完整基因微陣列矩陣的方法稱為MC方法。該方法和其他方法的對比實驗結果將在第4章中給出。

3 矩陣填充對LLSimpute的改進

由LLSimpute中的第1步可以看出來，與目標基因對應的k最近鄰基因行必須是完整的，才可以構成式中的A矩陣和B矩陣。同樣，在自適應地尋找最優的k時，也需要完整的矩陣。然而實際的基因微陣列矩陣往往是多行都含有缺失點。一種解決辦法是在找k最近鄰的時候，忽略掉含有缺失點的基因，只在完整的基因中尋找鄰居，這種辦法在缺失率比較低的時候是可行的，但是當缺失率較高時，大多數基因都是不完整的，按照這個方法，會帶來鄰居數量過少的問題。實驗也表明，LLSimpute對k選取是非常敏感的，若只在完整的行中尋找鄰居，在缺失率大于等于5%的時候，LLSimpute就會因可選取的鄰居數量過少，而使重建錯誤率大大升高。

由此可以看出，LLSimpute需要一個事先填充好的“中間矩陣”才可以構造A矩陣和B矩陣，以及尋找最優的k。對于這個“中間矩陣”，文獻［6］中采用行均值的方法得到，即在含缺失點的矩陣當中，用每行的均值來代替缺失點，得到完整的矩陣。在本文的第2章中，已經提到了用矩陣填充得到完整矩陣的方法。矩陣填充得到的結果，可以作為最終重建后的微陣列矩陣，也可以作為一個“中間矩陣”，繼續使用LLSimpute來進行重建。這樣，即將LLSimpute中的行均值填充的過程，替換為矩陣填充的過程。在本文中，將這種改進的LLSimpute稱為MC_LLS方法。

用行均值填充可以構成一個完整的矩陣，但行均值只利用了每個目標基因本行的信息。MC_LLS方法中的MC過程是以矩陣的低秩為優化目標，利用的是矩陣的全局特性，而LLSimpute本身，是基于矩陣的局部特征的，所以MC_LLS方法，實際上同時利用了矩陣的全局特征和局部特征。MC_LLS方法的實驗結果，也將在第4章中給出。

4 實驗與結果

對目前常用的幾種方法 (KNN、BPCA、LLSimpute)，和本文提出的MC、MC_LLS方法進行對比試驗，并給出結果的分析。

4.1 數據集

試驗在4個公開的真實基因微陣列數據集上進行。前兩個數據都屬于時間序列。其中第一個數據集來自［13］中的酵母細胞的CDC15序列和CDC28序列，將其稱為CDC15_28。第二個來自 ［13］中的酵母細胞的alpha序列，稱為SP_APLHA。第三個數據集來自文獻［14］中的人類癌細胞株，是一個非時間序列，稱為NCI60。第四個數據集來自［15］中的淋巴瘤細胞，也是一個非時間序列，稱為lymphoma。其中CDC15_28在綜述［1］中也被用于對各種方法進行比較，SP_APLHA在 ［6］中被用于檢驗LLSimpute的重建精確度，NCI60與綜述［1］中的NTS來自相同的細胞株，lymphoma在LSimpute［5］中也被采用。各個數據集的總行列數、完整的行列數 (即除去含缺失點的行之后的大小)以及類型可見表1。

表1 實驗采用的數據集

這些數據集原本都是不完整的，為了評價重建后的錯誤率，需要完整的微陣列矩陣作為參考。為此，將矩陣中不完整的行去掉，只保留完整的行 (如表1中的第三列)。再在這些完整的行中按不同的比例隨機去除一些已知點，作為缺失點。同樣的方法也在 ［1，6］中被采用。這樣，可以同時得到完整的微陣列矩陣和它們對應的不同缺失率下的缺失矩陣。對不同缺失率的矩陣進行重建，將重建后的矩陣和完整的矩陣對比，可以得出重建誤差。

4.2 參數選擇

對于KNN方法，在CDC15_28和SP_ALPHA數據集上，鄰居數k值設為10;在NCI60和lymphoma數據集上，k設為5，這兩個k值處于KNN最優的參數范圍區間內［1］。BPCA的參數設為其默認值，即微陣列矩陣的列數減去1。LLSimpute中的參數k使用 ［6］中模擬缺失點篩選k的方法得到的值。MC_LLS中的k也采用模擬缺失點的方法得到。需要說明的是，對于LLSimpute和MC_LLS，在尋找鄰居時，若缺失率足夠低 (此時完整的行數相對比較多)，則只從完整的基因行中尋找，而不是從用行均值或MC填充好的“中間矩陣”中尋找，以有效利用矩陣本身的真實值。

4.3 重建評價標準

目前對于基因重建效果的評價，多采用的是標準化根均方差 (normalized root mean squared error)［1，4，6］，稱為NRMSE，定義如下

式中:yj——基因第j位上原始的值，j——該位置上重建后的值，N——要重建的基因點的總數，即缺失點的個數，σy——N個原始基因點的標準差 (standard deviation)。NRMSE值越低，表明重建錯誤率越低。

4.4 重建結果

對4.1中的4個數據集，按照從1%到25%的之間的8種不同的缺失率，構造缺失點。對于每種數據集，不同缺失率下的重建都進行10次重復實驗，取10次標準化根均方差的平均值作為最終的重建錯誤率。圖1—圖4分別給出了不同方法在CDC15_28、SP_ALPHA、NCI60和lymphoma數據集上的結果。

圖1 在數據集CDC15_28上的重建錯誤率

圖1中，KNN的曲線在較高的缺失率下被截斷，表示在這些錯誤率下，基因的每行都存在缺失點，從而無法從完整的行中找到最近鄰。在缺失率小于10%的情況下，MC的重建錯誤率明顯低于KNN、BPCA和LLSimpute。但是隨著缺失率的升高，MC的錯誤率也開始變大。然而MC_LLS表現很穩定，在任何缺失率下，重建錯誤率都是最低的。

圖2 在數據集SP_ALPHA上的重建錯誤率

由圖2可以看出，在數據集SP_ALPHA下，KNN的表現仍然是最差的，重建錯誤率始終高于0.9。LLSimpute在缺失率較低的情況下有較低的錯誤率，但當缺失率大于5%時，錯誤率明顯提高。MC在缺失率為8%的時候，錯誤率小于LLSimpute，在其他缺失率下均高于LLSimpute、BPCA和MC_LLS。但是總體上，仍然是MC_LLS取得了最低的重建錯誤率。

圖3 在數據集NCI60上的重建錯誤率

NCI60數據集是一個非時間序列，具有很強的局部相關性，所以圖3中基于數據全局特征的BPCA和MC在這里的重建錯誤率明顯比LLS、MC_LLS高出很多。而基于局部相關性的KNN，在這里的表現雖然仍然比較差，但是已經取得了比時間序列下更低的錯誤率。對于MC_LLS，在大多數的缺失率下，仍然具有最低的重建錯誤率。

圖4 在數據集lymphoma上的重建錯誤率

與數據集NCI60類似，圖4中的數據集lymphoma也具有較強的局部相關性。LLSimpute和BPCA的表現非常接近，KNN的結果仍然是最差的。雖然MC本身的NRMSE高出LLSimpute和BPCA，但MC_LLS仍然在除1%之外的其他所有缺失率下取得了最低的重建錯誤率。

由以上實驗結果可以得出結論，基于矩陣全局特征的MC方法在局部相關性較弱的數據集上具有優勢，但是不適用于局部相關性很強的數據集。而MC_LLS由于結合了矩陣的整體相關性和局部相關性，所以在任何類型的數據集上都表現出了最好的重建結果。

表2所示的是各種方法在CDC15_28數據集上，不同缺失率下的運算時間，單位是秒。每個時間都是在進行了10次重復實驗后求得的平均值。實驗在一臺裝有4核Intel Xeon E5504 2.0GHz處理器的計算機上進行，采用Matlab R2011a。

表2 各種方法在數據集CDC15_28上的運算時間 (秒)

雖然矩陣填充中的奇異值分解通常是一項比較耗時的計算，但由于IALM使用了部分奇異值分解［12］，所需時間大大減少，所以MC方法在運算時間上具有明顯的優勢。而相對于LLSimpute，MC_LLS在時間開銷增加不明顯的情況下，卻明顯降低了重建錯誤率。

5 結束語

使用非精確增廣拉格朗日乘子 (IALM)解決了矩陣填充 (MC)問題，實現了基因微陣列矩陣的核范數凸優化。含有缺失點的微陣列矩陣，在進行矩陣填充后，可以作為最終的重建結果，并且矩陣填充得到的結果作為中間矩陣，還可以替換LLSimpute中的行均值矩陣。在四個常用的基因微陣列數據集上的實驗結果表明矩陣填充方法 (MC)在某些時間序列微陣列矩陣的部分缺失率下可以取得比KNN、LLSimpute和BPCA更高的重建精確度，并且計算時間具有明顯優勢。矩陣填充和LLSimpute結合的方法 (MC_LLS)，在實驗的四個數據集中的幾乎所有缺失率下，取得了最高的重建精確度，并且相對LLSimpute時間增加不明顯。

［1］Bras L P，Menezes J C.Dealing with gene expression missing data［J］.Syst Biol(Stevenage)，2006，153(3):105-19.

［2］Liew Alan Wee-Chung，Law Ngai-Fong，Yan Hong.Missing value imputation for gene expression data:Computational techniques to recover missing data from available information［J］.Briefings in Bioinformatics，2011，12(5):498-513.

［3］Troyanskaya Olga，Cantor Michael，Sherlock Gavin，et al.Missing value estimation methods for DNA microarrays［J］.Bioinformatics，2001，17(6):520.

［4］Oba Shigeyuki，Sato Masa-aki，Takemasa Ichiro，et al.A bayesian missing value estimation method for gene expression profile data［J］.Bioinformatics，2003，19(16):2088.

［5］B Trond Hellen，Dysvik Bjarte，Jonassen Inge.LSimpute:Accurate estimation of missing values in microarray data with least squares methods［J］.Nucleic Acids Research，2004，32(3):e34.

［6］Kim Hyunsoo，Golub Gene H，Park Haesun.Missing value estimation for DNA microarray gene expression data:Local least squares imputation［J］.Bioinformatics，2005，21(2):187.

［7］Brock Guy N，Shaffer John R，Blakesley Richard E，et al.Which missing value imputation method to use in expression profiles:A comparative study and two selection schemes［J］.BMCbioinformatics，2008，9(1):12.

［8］Magalie C，Alain M，Ga lle L.Comparative analysis of missing value imputation methods to improve clustering and interpretation of microarray experiments［J］.BMCgenomics，2010，11(1):15.

［9］Candès Emmanuel J，Recht Benjamin.Exact matrix completion via convex optimization［J］.Foundations of Computational Mathematics，2009，9(6):717-772.

［10］Bennett James，Lanning Stan.The netflix prize［R］.California，USA:KDDCup，2007.

［11］JI Hui，LIU Chaoqiang，SHEN Zuowei，et al.Robust video denoising using low rank matrix completion［C］//IEEE Conference on Computer Vision and Pattern Recognition，2010.

［12］LIN Zhouchen，CHEN Minming，MA Yi.The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices［J］.Arxiv preprint arXiv:1009.5055，2010.

［13］Spellman Paul T，Sherlock Gavin，ZHANG Michael Q，et al.Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization［J］.Molecular biology of the cell，1998，9(12):3273.

［14］Scherf Uwe，Ross Douglas T，Waltham Mark，et al.A gene expression database for the molecular pharmacology of cancer［J］.Nature genetics，2000，24(3):236-244.

［15］Alizadeh A Alizadeh，Eisen Michael B，Davis R.Eric，et al.Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling［J］.Nature，2000，403(6769):503-511.