王 蕾，李因果，夏利宇

(1.中國(guó)人民大學(xué) a.應(yīng)用統(tǒng)計(jì)科學(xué)研究中心，b.統(tǒng)計(jì)學(xué)院，北京 100872；2.江蘇師范大學(xué) 商學(xué)院，江蘇 徐州 221116)

一、引 言

主成分分析是經(jīng)典的數(shù)據(jù)處理和降維的工具，在工程、生物以及社會(huì)科學(xué)等領(lǐng)域都有著十分廣泛的應(yīng)用。主成分分析旨在尋找原始變量的一系列線性組合，使得由主成分張成的子空間能夠盡可能地刻畫(huà)原始數(shù)據(jù)的變動(dòng)情況。然而，在求解主成分的過(guò)程中并未考慮到響應(yīng)變量Y的信息，使其有效性受到限制，正如Cox所言：“主成分是由自變量X的邊際分布計(jì)算得來(lái)，我們沒(méi)有理由認(rèn)為主成分一定包含了跟Y有關(guān)的信息。”作為一種無(wú)監(jiān)督的降維方法，通過(guò)直接分解自變量X的協(xié)方差矩陣得到的主成分，并不一定是所感興趣的投影方向。在回歸和分類(lèi)問(wèn)題中，如果能將原始變量投影到與響應(yīng)變量有關(guān)的方向上去，則更可能得到較合理的結(jié)果，這也促使了一些有監(jiān)督的主成分分析方法的出現(xiàn)。

Oman在給定Y的條件下對(duì)X(即X|Y)建立了逆回歸模型，并認(rèn)為在給定X的條件下對(duì)Y(即Y|X)的回歸模型中系數(shù)向量應(yīng)當(dāng)被包含在或者靠近X前幾個(gè)主成分張成的空間[1]；Bair等在求解主成分之前，先用自變量與因變量的皮爾遜相關(guān)系數(shù)進(jìn)行了變量篩選[2]；Cook等提出了關(guān)于主成分的逆回歸模型及主成分?jǐn)M合回歸[3-4]；Barshan等用希爾伯特-施密特獨(dú)立性準(zhǔn)則度量主成分與因變量的相依性，旨在尋找一個(gè)使X與Y的相依程度最大化的主成分子空間[5]；Fuentes等將稀疏偏最小二乘的思想引入時(shí)間序列數(shù)據(jù)的預(yù)測(cè)中，即對(duì)同時(shí)包含X和Y的矩陣進(jìn)行特征根分解[6]；Giovannelli等在宏觀經(jīng)濟(jì)時(shí)間序列數(shù)據(jù)的預(yù)測(cè)問(wèn)題中，考慮了用標(biāo)準(zhǔn)化的主成分與響應(yīng)變量Y的相關(guān)性來(lái)選擇保留的主成分[7]。

上述文獻(xiàn)在估計(jì)主成分時(shí)結(jié)合了因變量的信息，但得到的每一個(gè)主成分都是所有原始變量的線性組合，而在自變量個(gè)數(shù)較多時(shí)則難以解釋每一個(gè)主成分的意義，這也促使了一些稀疏主成分分析方法的出現(xiàn)。Zou等將求解主成分的問(wèn)題轉(zhuǎn)化為線性回歸的最優(yōu)化問(wèn)題，并將最小絕對(duì)收縮和選擇算子以及彈性網(wǎng)的罰函數(shù)引入目標(biāo)函數(shù)中，使某些較小的載荷直接變?yōu)?[8]；Johnstone等證明了一定條件下稀疏主成分估計(jì)具有一致性[9]；Ma提出用迭代閾值法估計(jì)主成分子空間[10]；Kawano等在Zou的研究基礎(chǔ)上，在目標(biāo)函數(shù)中加入了主成分回歸的擬合誤差項(xiàng)[11]；Yi等提出聯(lián)合稀疏主成分分析，放松了各主成分的正交限制，且在目標(biāo)函數(shù)中直接對(duì)特征向量矩陣的稀疏性進(jìn)行懲罰[12]。國(guó)內(nèi)也有將稀疏主成分應(yīng)用于實(shí)際數(shù)據(jù)分析的文獻(xiàn)，如喻勝華[13]。

本文主要研究基于正交迭代的有監(jiān)督稀疏主成分分析方法(SSPCA)。當(dāng)前幾個(gè)特征值相等或很接近時(shí)，分別求解單個(gè)主成分可能出現(xiàn)不可識(shí)別的問(wèn)題，因此本文使用正交迭代方法求解由前幾個(gè)特征向量張成的特征根空間。正交迭代是計(jì)算特征根空間的常用方法，Ma基于此算法提出了計(jì)算稀疏主成分的ITSPCA及DTSPCA，前者是在正交迭代過(guò)程中增加了閾值步(Thresholding Step)，而后者選取了方差較大的一部分自變量計(jì)算主成分，再進(jìn)行正交分解，可用DTSPCA的結(jié)果作為ITSPCA的初值。ITSPCA和DTSPCA都是無(wú)監(jiān)督的方法，若用主成分作為新的自變量進(jìn)行預(yù)測(cè)，則無(wú)法識(shí)別出方差很大卻與Y無(wú)關(guān)的自變量。有監(jiān)督的主成分方法主要有：BSPCA(Bair’s SPCA)、ESPCA(Elnaz Barshan’s SPCA)、FSPCA(Fuentes’s SPCA)、GSPCA(Giovannelli’s SPCA)等，但FSPCA和GSPCA主要針對(duì)時(shí)間序列數(shù)據(jù)，因此在第三部分?jǐn)?shù)值模擬中不予比較。

本文使用正交迭代方法求解主成分子空間，并利用距離相關(guān)系數(shù)度量各自變量與因變量之間的相關(guān)性，將與Y的相關(guān)性較弱且方差較小的自變量對(duì)應(yīng)的載荷變?yōu)?，以得到稀疏的特征向量空間。相比ITSPCA和DTSPCA，本文給出的SSPCA更傾向于留下與Y有較強(qiáng)相關(guān)性的自變量，其預(yù)測(cè)效果也得到極大提高；對(duì)于BSPCA和ESPCA，SSPCA更易于識(shí)別X與Y的非線性關(guān)系，且能得到更稀疏的結(jié)果，使主成分更易解釋。

二、研究方法

(一)主成分逆回歸模型

考慮以下逆回歸模型：

X=μ+Γvy+σε

(1)

其中X為p維列自變量；μ為p維常數(shù)列向量；Γ為p×s實(shí)數(shù)矩陣，s

Oman將主成分逆回歸模型用于描述校準(zhǔn)問(wèn)題，如Y為某些量的精確值，通常難以獲得或者需要耗費(fèi)大量財(cái)力物力，而X為易得的近似測(cè)量結(jié)果，一般用X來(lái)推斷Y即建立關(guān)于Y|X的向前回歸模型，但近似測(cè)量結(jié)果往往也與精確值有關(guān)，即X的取值受到Y(jié)的影響，故在此類(lèi)問(wèn)題中同時(shí)建立關(guān)于X|Y的逆回歸模型也十分合理。

Cook認(rèn)為此逆回歸模型可以應(yīng)用于測(cè)量誤差、圖像識(shí)別、微陣列數(shù)據(jù)等領(lǐng)域，且在引理中闡述了逆回歸模型的系數(shù)向量Γ是向前回歸模型中的一個(gè)充分降維的方向，即Y|X與給定ΓX的條件下Y(即Y|ΓX)的分布相同[3]。此外，Γ列空間的極大似然估計(jì)恰為X的協(xié)差陣的前s個(gè)特征向量張成的空間。不難看出，在X與Y滿足此逆回歸模型的條件下，若使用ΓX代替原始自變量X來(lái)預(yù)測(cè)Y，在降維的同時(shí)還可以保持預(yù)測(cè)效果不變。但在實(shí)際應(yīng)用中，收集到的變量中只有一小部分是與Y相關(guān)的，即高維問(wèn)題中的稀疏性假設(shè)，因此若能在求解Γ時(shí)將那些沒(méi)有預(yù)測(cè)能力的自變量對(duì)應(yīng)的載荷變?yōu)?，則可以利用較少的變量得到較好的預(yù)測(cè)效果。

(二)正交迭代

設(shè)有p×s的列滿秩矩陣Z，p>s，存在p×s的列正交矩陣Q和s×s的上三角矩陣R，使得Z=QR為Z的輕正交分解。

給定p×p正定矩陣A、列相互正交的p×s矩陣Q0，則正交迭代按如下步驟產(chǎn)生正交矩陣序列{Qk}?p×s：

1.Zk=AQk-1

2.QkRk=Zk(輕正交分解)

且Qk的列空間會(huì)收斂到A的前s個(gè)特征向量張成的空間[14]367-368，即主成分子空間，其中Qk、Rk、Zk分別代表第k次迭代得到的矩陣Q、R、Z。

(三)基于正交迭代的有監(jiān)督的稀疏主成分分析

算法1輸入：自變量樣本協(xié)方差陣S；主成分子空間維度s；初始正交陣Q0及用Q0建模擬合Y的擬合誤差r0；閾值函數(shù)η，臨界值γ1j，γ2k，γ3k，j=1，2，…，s；稀疏容忍度Sparsity。

(1)從k=1開(kāi)始做。

(8)r0=rk

算法2輸入：自變量樣本協(xié)方差陣S；臨界值αn。

(1)選取自變量子集XB={Xi∶i∈B}，B為與Y相關(guān)性大于臨界值的自變量的下標(biāo)集合，B={i∶ωi≥αn}。

三、數(shù)值模擬

考慮用模型(1)產(chǎn)生數(shù)據(jù)，令μ=0、σ=1、s=1，2，3；Γ為服從[-10，10]均勻分布的p×s維矩陣列空間的一組正交基；vy是由{exp(Y)、I(Y<0)、Y}構(gòu)成的向量，即s=1時(shí)vy=exp(Y)，s=2時(shí)vy={exp(Y)，I(Y<0)}，依此類(lèi)推。Y和ε皆來(lái)自標(biāo)準(zhǔn)正態(tài)分布N(0，1)。設(shè)p=1 024，觀測(cè)個(gè)數(shù)n=200，其中0.7n=140個(gè)觀測(cè)為訓(xùn)練集，其余0.3n=60個(gè)觀測(cè)為測(cè)試集。假定1 024個(gè)自變量中分別只有前m=5，10，20，30，40個(gè)自變量與Y相關(guān)，由上述模型(1)產(chǎn)生，其余變量為服從t(3)分布的假變量。在Bair的數(shù)值模擬中，僅有1%的自變量與Y相關(guān)，故本文m的設(shè)置是合理的。用五種方法得到的主成分代替原始自變量，與Y建立四次多項(xiàng)式回歸模型。在固定m=20時(shí)，令n=120、160、200、240，隨觀測(cè)個(gè)數(shù)增加而研究觀察模擬結(jié)果的變化情況。

重復(fù)模擬100次，并將本文提出的方法SSPCA與ITSPCA、DTSPCA、BSPCA、ESPCA從擬合誤差、預(yù)測(cè)誤差和特征向量稀疏程度三方面進(jìn)行比較，其中擬合誤差用訓(xùn)練集的殘差平方和來(lái)衡量，預(yù)測(cè)誤差用測(cè)試集的殘差平方和來(lái)衡量，特征向量稀疏程度則是計(jì)算平均每個(gè)特征向量不為0的元素個(gè)數(shù)。BSPCA的臨界值選取使用R包superpc中的相關(guān)函數(shù)，ITSPCA和DTSPCA的有關(guān)參數(shù)均按原文推薦的設(shè)定，算法2中αn取0.03，模擬結(jié)果見(jiàn)表1～4所示。

表1～4中數(shù)據(jù)均為100次模擬結(jié)果的中位數(shù)，可看出在不同的設(shè)定下SSPCA和BSPCA的擬合誤差和預(yù)測(cè)誤差都優(yōu)于其他方法，且在SSPCA保留的自變量個(gè)數(shù)遠(yuǎn)小于BSPCA的情況下，其預(yù)測(cè)誤差仍然小于BSPCA。由表1表2易見(jiàn)，隨著s的增大，五種方法的擬合誤差基本都在減小，然而只有SSPCA和BSPCA的預(yù)測(cè)誤差也呈下降趨勢(shì)，其余三種方法的預(yù)測(cè)誤差大多有所增加。

因SSPCA與ITSPCA、DTSPCA皆為稀疏主成分方法，得到的特征子空間都比較稀疏，每個(gè)特征向量平均非零元素個(gè)數(shù)均在10～30之間。ESPCA的結(jié)果最不稀疏，主成分是全部原始變量的線性組合，保留了最多的信息，同時(shí)也保留了所有的假變量，因此其擬合和預(yù)測(cè)效果并不盡如人意。BSPCA在求解主成分之前利用皮爾遜相關(guān)系數(shù)進(jìn)行了變量選擇，使得特征向量具有一定的稀疏性，但相比較稀疏主成分方法，保留的自變量個(gè)數(shù)仍然較多。由表1～2知，盡管測(cè)試集的樣本量小于訓(xùn)練集，BSPCA的預(yù)測(cè)誤差大多大于擬合誤差，說(shuō)明BSPCA在本文設(shè)定下，可能存在過(guò)擬合的問(wèn)題。

表1 n為200時(shí)不同m的擬合誤差與預(yù)測(cè)誤差表

表2 m為20時(shí)不同n的擬合誤差與預(yù)測(cè)誤差表

表3 n為200時(shí)不同m的特征向量稀疏程度表

表4 m為20時(shí)不同n的特征向量稀疏程度表

四、實(shí)例分析

下面將SSPCA應(yīng)用于體脂數(shù)據(jù)body fat的分析，body fat數(shù)據(jù)集來(lái)自http://lib.stat.cmu.edu/datasets/bodyfat，旨在利用年齡、體重、身高等基本數(shù)據(jù)估算人體去脂體重，數(shù)據(jù)共有252個(gè)觀測(cè)、15個(gè)變量，其中前143個(gè)觀測(cè)為訓(xùn)練集，其余為測(cè)試集。人體密度和體脂率有近似一一對(duì)應(yīng)的關(guān)系，因此將人體密度這一變量去掉，同時(shí)還向數(shù)據(jù)中添加了1 011個(gè)服從t(3)分布的假變量，以體脂率為因變量，其余為自變量建立回歸模型，體脂率的數(shù)據(jù)用百分?jǐn)?shù)表示，變動(dòng)范圍在0～100之間，訓(xùn)練集體脂率的箱線圖見(jiàn)圖1，從圖1中可以看到體脂率主要在5～40之間波動(dòng)。將數(shù)值模擬中的五種方法應(yīng)用于body fat數(shù)據(jù)集中，分別選取前1～6個(gè)主成分，建立二次多項(xiàng)式回歸模型(結(jié)果見(jiàn)表5表6所示)，易見(jiàn)在用前3～6個(gè)主成分建模時(shí)SSPCA都能達(dá)到較優(yōu)的擬合誤差與預(yù)測(cè)誤差，在用前5個(gè)主成分建模時(shí)SSPCA的預(yù)測(cè)誤差達(dá)到最小。再結(jié)合五種方法各特征向量平均保留的變量個(gè)數(shù)來(lái)看，SSPCA方法優(yōu)于其他四種方法。

圖1 訓(xùn)練集體脂率的箱線圖

sSSPCAITSPCADTSPCABSPCAESPCA擬合誤差預(yù)測(cè)誤差擬合誤差預(yù)測(cè)誤差擬合誤差預(yù)測(cè)誤差擬合誤差預(yù)測(cè)誤差擬合誤差預(yù)測(cè)誤差15 808.84 003.65 826.84 017.45 831.04 023.35 807.74 031.68 673.17 622.524 635.83 511.64 613.33 507.24 618.43 516.94 587.53 468.28 620.87 133.732 676.22 440.74 134.13 332.04 109.33 324.34 073.93 303.88 610.27 132.342 572.22 368.74 052.13 260.34 093.73 296.34 037.53 289.78 531.17 242.752 485.42 266.62 758.22 498.32 831.72 504.22 812.92 742.68 441.77 165.062 426.12 286.12 692.92 563.32 770.92 566.32 719.62 844.28 089.27 764.1

表6 不同方法分析體脂數(shù)據(jù)的特征向量稀疏程度表

五、結(jié) 論

本文基于正交迭代和距離相關(guān)系數(shù)，提出一種有監(jiān)督的稀疏主成分分析方法SSPCA。該方法考慮了自變量與因變量之間的相關(guān)性，并在迭代求解的過(guò)程中將一些與Y相關(guān)性很弱的自變量對(duì)應(yīng)的系數(shù)變?yōu)?，使得所求的特征向量只保留預(yù)測(cè)能力較強(qiáng)的信息。相比ITSPCA、DTSPCA、BSPCA和ESPCA方法，SSPCA方法在數(shù)值模擬和實(shí)例分析中都表現(xiàn)出了顯著的優(yōu)勢(shì)。

參考文獻(xiàn)：

[1] Oman S D.Random Calibration with Many Measurements：An Application of Stein Estimation[J].Technometrics，1991(2).

[2] Bair E，Hastie T，Paul D，Tibshirani R.Prediction by Supervised Principal Components[J].Journal of the American Statistical Association，2006，101(473).

[3] Cook R D.Fisher Lecture：Dimension Reduction in Regression[J].Statistical Science，2007，22(1).

[4] Cook R D，F(xiàn)orzani L.Principal Fitted Components for Dimension Reduction in Regression[J].Statistical Science，2008(4).

[5] Barshan E，Ghodsi A，Azimifar Z，Jahromi M Z.Supervised Principal Component Analysis：Visualization，Classification and Regression on Subspaces and Submanifolds[J].Pattern Recognition，2011，44(7).

[6] Fuentes J，Poncela P，Rodríguez，J.Sparse Partial Least Squares in Time Series for Macroeconomic Forecasting[J].Journal of Applied Econometrics，2015，30(4).

[7] Giovannelli A，Proietti T.On the Selection of Common Factors for Macroeconomic Forecasting[J].MPRA Paper，2015.

[8] Zou H，Hastie T，Tibshirani R.Sparse Principal Component Analysis[J].Journal of Computational and Graphical Statistics，2006，15(2).

[9] Johnstone I M，Lu A Y.On Consistency and Sparsity for Principal Components Analysis in High Dimensions[J].Journal of the American Statistical Association，2009，104(486).

[10] Ma Z.Sparse Principal Component Analysis and Iterative Thresholding[J].The Annals of Statistics，2013，41(2).

[11] Kawano S，F(xiàn)ujisawa H，Takada T.Sparse Principal Component Regression with Adaptive Loading[J].Computational Statistics and Data Analysis，2015，89(C).

[12] Yi S，Lai Z，He Z，Cheung Y，Liu Y.Joint Sparse Principal Component Analysis[J].Pattern Recognition，2017，61.

[13] 喻勝華.基于穩(wěn)健稀疏主成分的經(jīng)濟(jì)增長(zhǎng)影響因素分析[J].統(tǒng)計(jì)與信息論壇，2017，32(3).

[14] Golub G H，Van Loan C F.Matrix Computation[M].Baltimore：Johns Hopkins University Press，2013.

[15] Li R，Zhong W，Zhu L.Feature Screening Via Distance Correlation Learning[J].Journal of the American Statistical Association，2012，107(499).