本發(fā)明涉及大數(shù)據(jù)分析技術領域,特別涉及一種處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法。
背景技術:
主成分分析,即PCA(Principal Component Analysis),是一種常用的數(shù)據(jù)分析方法。PCA通過矩陣計算提取出原始數(shù)據(jù)在線性空間中的一組主要基向量(即主要特征分量),然后將原始數(shù)據(jù)在這組基上投影,實現(xiàn)高維數(shù)據(jù)的降維。對經(jīng)過降維后的數(shù)據(jù),可進一步做聚類、分類等運算,實現(xiàn)特征提取、自動分類、識別等人工智能應用。當前,主成分分析作為一種重要的無監(jiān)督學習方法,已廣泛用于數(shù)據(jù)挖掘、機器學習有關的各種應用問題中。
在實際問題中,數(shù)據(jù)往往可以表示為一個矩陣。不失一般性,將每個數(shù)據(jù)看成矩陣A的一行,那么矩陣的列數(shù)就是每個數(shù)據(jù)的維度。主成分分析計算的目標是原始數(shù)據(jù)的若干個主要特征分量,可通過矩陣的特征值分解或奇異值分解得到。基于矩陣特征值分解的方法是先計算矩陣ATA,然后對ATA進行特征值分解,得到最大的若干特征值對應的特征向量就是要求的“主成分”?;诰仃嚻娈愔捣纸獾姆椒ㄖ苯訉仃嘇做奇異值分解:A=UΣVT,其中U和V均為正交陣,Σ為對角元從大到小排列的對角陣,得到的V矩陣的前若干列就是要求的“主成分”。若數(shù)據(jù)維度不太高,即A的列數(shù)遠小于行數(shù),基于特征值分解的方法計算效率比較高,因為其處理的ATA矩陣是一個階數(shù)較小的矩陣。
另一方面,隨著移動設備、互聯(lián)網(wǎng)、傳感器網(wǎng)絡、基因工程的迅速發(fā)展,產(chǎn)生數(shù)據(jù)的來源變得多樣化,同時數(shù)據(jù)量也呈現(xiàn)出指數(shù)級的增長趨勢。也就是說,當前正處在所謂的“大數(shù)據(jù)”時代。如何在可承受的時間空間限制下存儲、分析和管理日益增長的數(shù)據(jù)集成為傳統(tǒng)的數(shù)據(jù)處理手段面臨的一個難題。研究表明,目前85%的數(shù)據(jù)都可以直接或通過轉換后表示為數(shù)值型數(shù)據(jù),即常見的整型、浮點型數(shù)據(jù),而數(shù)據(jù)庫中存儲數(shù)值型數(shù)據(jù)構造的“表”結構通常被當作矩陣進行處理。因此,如何針對這些大數(shù)據(jù)在產(chǎn)生、存儲、應用等方面的特點,研究出有效的“大矩陣”數(shù)據(jù)分析方法變得異常重要。具體來說,由于數(shù)據(jù)規(guī)模太大,它們可能是分布式存儲的(即在網(wǎng)絡上不同的計算機節(jié)點上)、或存儲在計算機硬盤上且無法完整地載入到內存中(由于內存容量限制)。在其他一些應用場景中,這些數(shù)據(jù)也可能是按“數(shù)據(jù)流”的方式逐漸生產(chǎn)、獲取到的,不適合采用傳統(tǒng)的先存儲下來再計算的方式對它們進行處理。考慮到傳統(tǒng)的計算主成分分析的方法需對整個矩陣進行特征值分解或奇異值分解,其內在算法需要反復的讀取、遍歷數(shù)據(jù)矩陣的元素(若要計算前k個主成分,至少要完整地讀取矩陣元素k遍),顯然它們不適合對上述場景中讀取開銷巨大的大數(shù)據(jù)進行分析。
考慮到上述背景,基于隨機化的矩陣計算方法,包括特征值分解、奇異值分解的算法,在近年來備受人們關注。在文獻:N.Halko,P.-G.Martinsson and J.A.Tropp,Finding structure with randomness:Probabilistic algorithms for constructing approximate matrix decompositions,SIAM Review,53(2011),no.2,pp.217-288(以下簡寫作SIAM2011)中,提出了一種對矩陣數(shù)據(jù)遍歷次數(shù)較少的隨機奇異值分解算法。該方法通過將原始矩陣A乘以一個僅含k列的隨機矩陣,得到原始矩陣列空間的k維特征子空間,然后求出該子空間的正交基向量矩陣Q,以及A的近似分解:A≈QB,其中B為一個只有k行的矩陣。最后對B這個較小的矩陣進行傳統(tǒng)奇異值分解計算,可近似得到原始矩陣A的前k個奇異值和相應的左、右奇異向量。在文獻SIAM2011中,還對上述近似算法的準確度進行了理論分析,結果表明它在很大的概率上能使誤差落在很小的限度內,同時也提出了幾種提高結果準確度的技巧。
應當指出,文獻SIAM2011所提方法雖然相比傳統(tǒng)的奇異值分解算法大大減少了對矩陣元素的遍歷次數(shù),但它至少需要遍歷矩陣元素兩遍,從計算效率上仍有提升空間,并且無法適應數(shù)據(jù)流式大數(shù)據(jù)的處理要求。
技術實現(xiàn)要素:
本發(fā)明旨在至少解決上述技術問題之一。
為此,本發(fā)明的目的在于提出一種處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法,該方法適合于多種大數(shù)據(jù)分析場景,具有較高的計算效率和實用性。
為了實現(xiàn)上述目的,本發(fā)明的實施例公開了一種處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法,包括以下步驟:S1:在內存中生成一個n行、l列的隨機數(shù)矩陣Ω;S2:選取原始數(shù)據(jù)矩陣A,并根據(jù)所述原始數(shù)據(jù)矩陣A計算矩陣G和H,并將矩陣G和H存儲于內存中,其中,G=AΩ,H=ATG,所述原始數(shù)據(jù)矩陣A為m×n矩陣;S3:初始化變量j=1,并初始化m×l矩陣Q和l×n矩陣B均為零矩陣;S4:設定G[j,j+b]和Ω[j,j+b]分別為矩陣G和矩陣Ω的第j到j+b列,且當j>1時,計算G[j,j+b]-QBΩ[j,j+b],并將計算結果覆蓋G[j,j+b],其中,b為不超過l-j的非負整數(shù);S5:對矩陣G[j,j+b]做簡化QR分解,得到m×(b+1)的列正交矩陣Q[j,j+b]和上三角方陣R,其中,Q[j,j+b]為存儲在矩陣Q的第j到j+b列;S6:如果j>1,則計算矩陣Q[j,j+b]-Q(QTQ[j,j+b])的簡化QR分解,將得到的m×(b+1)列正交矩陣覆蓋Q[j,j+b],以得到上三角陣為并計算矩陣乘法并將計算結果覆蓋R;S7:設H[j,j+b]表示矩陣H的第j到j+b列,如果j=1,計算否則計算得到結果為(b+1)×n的矩陣Btemp,并將Btemp存儲在矩陣B的第j到j+b行;S8:將變量j+b+1的值賦值給變量j;S9:如果j≤l,則返回執(zhí)行所述S4,否則執(zhí)行所述S10;S10:對矩陣B做奇異值分解:B=UΣVT,其中,矩陣V的前k列為所述前k個主成分向量,Σ的前k個對角元為所述對應的奇異值。
另外,根據(jù)本發(fā)明上述實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法還可以具有如下附加的技術特征:
在一些示例中,在所述S1中,所述參數(shù)l為至少比k大5的整數(shù)。
在一些示例中,所述S1,進一步包括:S11:根據(jù)隨機數(shù)生成器軟件生成一個n×l隨機數(shù)矩陣Ω;S12:初始化變量i=0,變量P為小于10的非負整數(shù);S13:如果i=P,則結束執(zhí)行,否則轉到所述S14繼續(xù)執(zhí)行;S14:計算矩陣乘法AΩ,并對計算結果進行簡化QR分解,將得到的m×l列正交陣賦值給矩陣G;S15:計算矩陣乘法ATG,并對計算結果進行簡化QR分解,將得到的n×l列正交陣賦值給矩陣Ω;S16:將i的值加1,并轉到所述S13繼續(xù)執(zhí)行。
在一些示例中,在所述S2中,根據(jù)所述原始數(shù)據(jù)矩陣A的不同產(chǎn)生方式或來源,通過遍歷一遍所述原始數(shù)據(jù)矩陣A中的元素來計算出矩陣G=AΩ和H=ATG。
在一些示例中,所述S2,進一步包括:S21:在內存中開辟二維數(shù)組空間存儲n×l的矩陣H,并將所述矩陣H的數(shù)據(jù)初始化為0;S22:獲取原始數(shù)據(jù)矩陣A的預設行的數(shù)據(jù)并存于內存中,并設定所述預設行形成s×n的矩陣Ai,計算矩陣乘運算Gi=AiΩ,其中,所述Gi為矩陣G對應的行;S23:計算并將計算結果賦值給矩陣H;S24:判斷是否獲取原始數(shù)據(jù)矩陣A的所有行,如果是,則停止執(zhí)行,否則返回執(zhí)行所述S22。
根據(jù)本發(fā)明實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法,基于目前的隨機奇異值分解算法,但通過改進將算法主體部分中對數(shù)據(jù)矩陣的遍歷次數(shù)由兩次減少為一次,并且保持原算法的準確度不變,從而可以適合于多種大數(shù)據(jù)分析場景,具有較高的計算效率和實用性。
本發(fā)明的附加方面和優(yōu)點將在下面的描述中部分給出,部分將從下面的描述中變得明顯,或通過本發(fā)明的實踐了解到。
附圖說明
本發(fā)明的上述和/或附加的方面和優(yōu)點從結合下面附圖對實施例的描述中將變得明顯和容易理解,其中:
圖1是根據(jù)本發(fā)明實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法的流程圖。
具體實施方式
下面詳細描述本發(fā)明的實施例,所述實施例的示例在附圖中示出。
以下結合附圖描述根據(jù)本發(fā)明實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法。
圖1是根據(jù)本發(fā)明一個實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法的流程圖。如圖1所示,該方法包括以下步驟:
步驟S1:在內存中生成一個n行、l列的隨機數(shù)矩陣Ω。其中,參數(shù)l為至少比k大5的整數(shù),k為一個隨機矩陣的列數(shù)。
在本發(fā)明的一個實施例中,隨機數(shù)矩陣Ω中的元素為均勻分布隨機數(shù)或標準正態(tài)分布隨機數(shù),或者通過更復雜的方式得到,此處不再一一列列舉贅述。
在本發(fā)明的一個實施例中,步驟S1進一步包括:
S11:根據(jù)隨機數(shù)生成器軟件生成一個n×l隨機數(shù)矩陣Ω。
S12:初始化變量i=0,變量P為小于10的非負整數(shù)。
S13:如果i=P,則結束執(zhí)行,否則轉到步驟S14繼續(xù)執(zhí)行。
S14:計算矩陣乘法AΩ,并對計算結果進行簡化QR分解,將得到的m×l列正交陣賦值給矩陣G。
S15:計算矩陣乘法ATG,并對計算結果進行簡化QR分解,將得到的n×l列正交陣賦值給矩陣Ω。
S16:將i的值加1,并轉到S13繼續(xù)執(zhí)行,從而提高結果的準確度。
步驟S2:選取原始數(shù)據(jù)矩陣A,并根據(jù)原始數(shù)據(jù)矩陣A計算矩陣G和H,并將矩陣G和H存儲于內存中,其中,G=AΩ,H=ATG,原始數(shù)據(jù)矩陣A為m×n矩陣,其每行代表一個數(shù)據(jù),總共m行。
在本發(fā)明的一個實施例中,在步驟S2中,根據(jù)原始數(shù)據(jù)矩陣A的不同產(chǎn)生方式或來源,通過遍歷一遍原始數(shù)據(jù)矩陣A中的元素來計算出矩陣G=AΩ和H=ATG。
基于此,步驟S2進一步包括:
S21:在內存中開辟二維數(shù)組空間存儲n×l的矩陣H,并將矩陣H的數(shù)據(jù)初始化為0。
S22:獲取原始數(shù)據(jù)矩陣A的預設行的數(shù)據(jù)并存于內存中,并設定預設行形成s×n的矩陣Ai,計算矩陣乘運算Gi=AiΩ,其中,Gi為矩陣G對應的行。
S23:計算并將計算結果賦值給(覆蓋)矩陣H。
S24:判斷是否獲取原始數(shù)據(jù)矩陣A的所有行,如果是,則停止執(zhí)行,否則返回執(zhí)行步驟S22。
步驟S3:初始化變量j=1,并初始化m×l矩陣Q和l×n矩陣B均為零矩陣。
步驟S4:設定G[j,j+b]和Ω[j,j+b]分別為矩陣G和矩陣Ω的第j到j+b列,且當j>1時,計算G[j,j+b]-QBΩ[j,j+b],并將計算結果覆蓋G[j,j+b],其中,設定b為不超過l-j的非負整數(shù)。
步驟S5:對矩陣G[j,j+b]做簡化QR分解,得到m×(b+1)的列正交矩陣Q[j,j+b]和上三角方陣R,其中,Q[j,j+b]存儲在矩陣Q的第j到j+b列。
步驟S6:如果j>1,則計算矩陣Q[j,j+b]-Q(QTQ[j,j+b])的簡化QR分解,將得到的m×(b+1)列正交矩陣覆蓋Q[j,j+b],以得到上三角陣為并計算矩陣乘法并將計算結果覆蓋R。
步驟S7:設H[j,j+b]表示矩陣H的第j到j+b列,如果j=1,計算否則計算得到結果為(b+1)×n的矩陣Btemp,并將Btemp存儲在矩陣B的第j到j+b行。
步驟S8:將變量j+b+1的值賦值給變量j。
步驟S9:如果j≤l,則返回執(zhí)行S4,否則執(zhí)行S10。
步驟S10:對矩陣B做奇異值分解,以得到前k個主成分向量和對應的奇異值。具體地,對矩陣B做奇異值分解的公式為:
B=UΣVT,
其中,矩陣V的前k列為前k個主成分向量,Σ的前k個對角元為對應的奇異值。
需要說的是,在本發(fā)明的上述實施例中,假設原始數(shù)據(jù)矩陣A為m×n矩陣,其每行代表一個數(shù)據(jù),總共m行,基于此,本發(fā)明的目的是計算出該數(shù)據(jù)的前k個主成分向量和對應的奇異值。另一方面,如果A的每列代表一個數(shù)據(jù),可對其轉置AT進行上述本發(fā)明上述實施例所描述的計算過程。
為了便于更好地理解本發(fā)明實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法,以下結合具體的實施例對本發(fā)明進行進一步詳細描述。
在本實施例中,本發(fā)明實施例的方法例如可用任何一種編程語言實現(xiàn),在具有CPU和內存的計算設備上執(zhí)行。本實施例中使用的隨機數(shù)生成器、對矩陣執(zhí)行加/減法、乘法、轉置、矩陣求逆(或線性方程組求解),以及QR分解和奇異值分解,均為現(xiàn)有技術,通過調用相應編程語言的數(shù)值計算函數(shù)庫可以實現(xiàn)。
在本實施例中,考慮一個存儲在硬盤上的大規(guī)模數(shù)據(jù)矩陣A,它是某種時間序列數(shù)據(jù),每個數(shù)據(jù)特征數(shù)比較多,比如100,000個,而數(shù)據(jù)總共有500,000個,每個數(shù)值采用雙精度浮點數(shù)存儲,需要8個字節(jié)。所以,整個數(shù)據(jù)需要約400GB的存儲量。假設需要提取出數(shù)據(jù)的1000個主要特征,也就是要算前1000個主成分向量??蓤?zhí)行下述步驟,只遍歷數(shù)據(jù)一遍即完成計算(參數(shù)值k=1000,m=500,000,n=100,000),具體步驟如下:
步驟1:根據(jù)參數(shù)k確定l的值,令l=k+5。
步驟2:在內存中生成一個n行、l列的標準正態(tài)分布隨機數(shù)矩陣Ω。
步驟3:通過遍歷硬盤上的矩陣A一遍來計算出矩陣G=AΩ和H=ATG。具體步驟如下:
步驟3.1:在內存中開辟二維數(shù)組空間存儲n×l的矩陣H和m×l矩陣G,將其數(shù)據(jù)初始化為0,打開存儲矩陣A的硬盤文件,將讀指針置于文件頭。
步驟3.2:從文件指針位置開始讀取矩陣A的1000行數(shù)據(jù)存于內存中,設它們形成矩陣Ai,計算Gi=AiΩ,結果Gi存于矩陣G對應的行上。
步驟3.3:計算將結果賦值給(覆蓋)H。
步驟3.4:如果沒有取完A對應的文件,返回步驟3.2執(zhí)行,否則執(zhí)行步驟4。
步驟4:初始化變量j=1,初始化m×l矩陣Q和l×n矩陣B均為零矩陣。
步驟5:設定b的值為19。
步驟6:記G[j,j+b]和Ω[j,j+b]分別為矩陣G和Ω的第j到j+b列;如果j>1,計算G[j,j+b]-QBΩ[j,j+b],將結果覆蓋G[j,j+b]。
步驟7:對矩陣G[j,j+b]做簡化QR分解,得到m×(b+1)的列正交矩陣Q[j,j+b]和上三角方陣R,Q[j,j+b]即存儲在矩陣Q的第j到j+b列。
步驟8:如果j>1,計算矩陣Q[j,j+b]-Q(QTQ[j,j+b])的簡化QR分解,將得到的m×(b+1)列正交矩陣覆蓋Q[j,j+b],而得到上三角陣為計算矩陣乘法結果覆蓋R。
步驟9:設H[j,j+b]表示矩陣H的第j到j+b列;如果j=1,計算否則計算得到結果為(b+1)×n的矩陣Btemp,將它存儲在矩陣B的第j到j+b行。
步驟10:將變量j+b+1的值賦值給變量j。
步驟11:如果j≤l,返回步驟5執(zhí)行,否則執(zhí)行步驟12。
步驟12:對矩陣B做奇異值分解,即B=UΣVT,則矩陣V的前k列就是所要的前k個主成分向量,Σ的前k個對角元就是相應的奇異值。
由于從硬盤讀數(shù)據(jù)的時間遠大于在內存中進行計算的時間,上述僅讀取了硬盤數(shù)據(jù)一遍的算法其執(zhí)行速度是文獻SIAM2011中相應算法的兩倍,大大節(jié)省了整個大數(shù)據(jù)分析的時間。
上述過程得的前1000個主成分向量組成的矩陣V,注意它是一個100,000×1000的列正交陣。然后計算AV,得到一個500,000×1000的降維數(shù)據(jù)矩陣,它仍表示500,000個數(shù)據(jù),但數(shù)據(jù)維度大大縮減。利用降維后的數(shù)據(jù)可進一步進行聚類、分類,等等,做數(shù)據(jù)挖掘與分析。
綜上,根據(jù)本發(fā)明實施例的處理大規(guī)模矩陣數(shù)據(jù)的主成分分析方法,基于目前的隨機奇異值分解算法,但通過改進將算法主體部分中對數(shù)據(jù)矩陣的遍歷次數(shù)由兩次減少為一次,并且保持原算法的準確度不變,從而可以適合于多種大數(shù)據(jù)分析場景,具有較高的計算效率和實用性。
在本說明書的描述中,參考術語“一個實施例”、“一些實施例”、“示例”、“具體示例”、或“一些示例”等的描述意指結合該實施例或示例描述的具體特征、結構、材料或者特點包含于本發(fā)明的至少一個實施例或示例中。在本說明書中,對上述術語的示意性表述不一定指的是相同的實施例或示例。而且,描述的具體特征、結構、材料或者特點可以在任何的一個或多個實施例或示例中以合適的方式結合。
盡管已經(jīng)示出和描述了本發(fā)明的實施例,本領域的普通技術人員可以理解:在不脫離本發(fā)明的原理和宗旨的情況下可以對這些實施例進行多種變化、修改、替換和變型,本發(fā)明的范圍由權利要求及其等同限定。