快速多維標度算法研究*

屈太國，蔡自興

1.衡陽師范學院 計算機科學與技術學院，湖南 衡陽 421002

2.智能信息處理與應用湖南省重點實驗室，湖南 衡陽 421002

3.中南大學 信息科學與工程學院，長沙 410083

1　引言

隨著大數據時代的到來，數據規模越來越大[1]，數據高維化趨勢越來越明顯，這帶來了一系列問題，如“維數災難”[2]。數據降維是應對數據高維化的一種有效方法。

經典多維標度法（classicalmultidimensionalscaling，CMDS）基于樣本之間的距離，求它們在低維歐氏空間中的坐標，使它們在歐氏空間的距離盡量逼近原來的距離[3]。CMDS是數據降維和數據可視化的常用方法[4-8]，被廣泛應用于Matlab、SAS、SPSS、Statistica、S-Plus等計算機語言中。

CMDS的優點是具有解析解，缺點是速度慢。它的處理時間為Θ(N3)，其中N表示樣本個數。隨著數據（樣本）規模越來越大，提高CMDS的速度成為一個迫切需要解決的問題。人們提出了多種CMDS的快速算法。文獻[9-11]基于彈簧-質點模型（springmass model），通過最小化代價函數（cost function），求樣本在低維歐氏空間的坐標。Chalmers算法[9]的時間為Θ(N2)；Morrison算法[10]將時間減少到Θ(NlgN)，但它只能給出樣本在二維歐氏空間的坐標；Williams等人[11]對Chalmers算法進行了一些改進，使之能適應高維大規模樣本集。這3種算法都屬于迭代算法，容易陷入局部最小。另一類快速算法[12-15]屬于Nystrom算法[16]。其中FastMap[12]將樣本投影到一組相互正交的樞軸（pivot）上，這些投影構成了樣本在低維歐氏空間中的坐標。HyperMap[13]對FastMap進行了推廣，提供了“多視角”分析數據的靈活性。MetricMap[14]將目標空間由FastMap中的歐氏空間推廣到偽歐氏（pseudo-Euclidean）空間。LMDS（landmark multidimensional scaling）[15]首先指定部分樣本為標志點（landmark point），利用CMDS求解標志點之間的距離矩陣，得到標志點在低維空間的坐標，其余點的坐標由一個線性公式給出。上述所有算法在速度上都比CMDS快，但給出的都是近似解。因此，如何快速求得與CMDS完全一致的解仍有待研究。

作者先后提出了一種改進的FastMap算法（improved FastMap，iFastMap）[17]和基于分治策略的多維標度算法（divide-and-conquer based MDS，dcMDS）[18]。這兩種算法和本文提出的iLMDS（improved LMDS）算法都能快速得到與CMDS一致的解。

2　預備知識

CMDS從樣本集的距離矩陣出發，求樣本在歐氏空間中的坐標，使它們在歐氏空間的距離盡量逼近原來的距離。

本文將距離限定為歐氏距離，即存在正整數m和x1,x2,…,xN∈Rm，滿足：

其中，xi=(xi1,xi2,…,xim)T，i=1,2,…,N表示樣本的坐標，N表示樣本集的大小；dij表示樣本間的歐氏距離。

本文采用X=(xij)=(x1,x2,…,xN)、D=(dij)分別表示樣本集的坐標矩陣和距離矩陣。在CMDS問題中，已知的是距離矩陣，而坐標矩陣往往是未知的，但這個未知量在快速多維標度算法的討論中起到非常關鍵的作用。

基于距離矩陣，可以構造下面兩個矩陣：

其中：

通過對矩陣B進行譜分解可以得到CMDS的解析解（如表1），本文稱之為樣本的CMDS坐標。類似的，本文中其他算法的結果，都稱為相應算法的坐標。為區別起見，稱X為樣本集的原始坐標矩陣。

CMDS需要對矩陣B進行譜分解，其時間復雜度為Θ(N3)。隨著樣本規模的擴大，CMDS的時間急劇增加，從而限制了其在大規模樣本集上的應用。

2.1　CMDS和PCA的等價關系

主成分分析（principal component analysis，PCA）是另一種常用的多元數據分析方法。它從樣本集的坐標矩陣出發，旨在找到樣本方差最大的方向，然后把樣本投影到這些方向上。這些投影構成PCA坐標。

CMDS與PCA算法歸納如下，見表1。

Table 1　CMDS and PCA表1　CMDS和PCA

由表1可知，PCA坐標是通過對協方差矩陣Σ=(XH)(XH)T進行譜分解得到的。

令Λm=diag(λ1,λ2,…,λm)，UX=(u1,u2,…,um)，則協方差矩陣可以寫成如下譜分解的形式：

X的PCA坐標矩陣可以寫成如下矩陣形式：

其中，u1,u2,…,um是方差最大的方向，稱為樣本集的主軸向量；UX=(u1,u2,…,um)稱為X的主軸矩陣。

根據文獻[3]，對歐氏距離，樣本集的CMDS坐標矩陣與PCA坐標矩陣相等，即：

這表明：CMDS坐標從本質上可視為樣本在各主軸向量上的投影，也可視為對原始坐標的一種變換。這種變換屬于下面要介紹的保距變換。

2.2　保距變換

定義1（保距變換）給定常向量r∈Rm和m階正交矩陣O，如下變換稱為保距變換[20]：

用二元組表示。y稱為x在下的像。不難驗證，任意兩點間的歐氏距離在變換前后保持不變。

對樣本集X=(x1,x2,…,xN)，令Y=(y1,y2,…,yN)，其中yi=O(xi-r)，可得X在下的像矩陣Y=

關于保距變換，有如下定理。

定理1坐標矩陣與其像矩陣具有相同的PCA坐標矩陣。

證明假定X為m×N矩陣，r∈Rm，O為m階正交矩陣，Y為X在下的像矩陣，即Y=O(X-由1TNH=0可知：

由此可得Y的協方差矩陣：

在上式證明中用到了式（4）。令UY=OUX，有：

上式就是ΣY的譜分解。根據式（5），可得Y的PCA坐標矩陣：

定理1表明，如能得到X的像矩陣，則只需對它進行PCA就可以得到樣本集的CMDS坐標。

與CMDS基于整個樣本集的距離信息不同，本文介紹的3類算法利用樣本集上滿足特定條件的某個子集的距離信息，得到X的像矩陣，從而求得與CMDS完全一致的解，并且提高了速度。下面介紹的內在維數規定了樣本子集所需滿足的條件。

2.3　樣本集的內在維數

定義2（內在維數[21]）樣本集內在維數指的是包含樣本集上所有點的最小歐氏空間的維數。

根據文獻[21]，rank(XN-1-xN1TN-1)就是樣本集內在維數，其中XN-1=(x1,x2,…,xN-1)，1N-1表示分量均為1的N-1維列向量。

不難驗證：

根據式（3），有：

本文采用m表示樣本集的內在維數，即：

為了得到一個內在維數等于m的子集，可以從樣本集中隨機抽取部分點。本文采用最簡單的策略，即抽取前n個樣本。顯然，內在維數等于m的子集至少包含m+1個樣本，即n≥m+1。如果子集的內在維數小于m，逐步增加n的值，直至子集內在維數等于m。這就是本文要采用的子集選擇算法。

算法1子集選擇算法

輸入：距離矩陣D=(dij),內在維數m。

輸出：n。

1.令n=m+1；

2.利用式（6）計算前n個點的內在維數m1；

3.如果m1=m，返回n，否則，n=n+1，轉2。

3　快速多維標度算法

首先回顧兩種快速多維標度算法iFastMap和dcMDS，然后介紹一種新的快速算法iLMDS。iFast-Map、iLMDS分別為FastMap和LMDS的改進算法。

3.1　iFastMap算法

FastMap算法包含多次投影，每次投影包含3步：首先，選取兩個相距較遠的樣本構成一個樞軸；然后，利用余弦定理，將各樣本投影到樞軸上；最后，修改樣本之間的距離。

第k次投影后，各點之間的距離為:

FastMap算法的時間主要用于距離的計算，每次投影后，要重新計算所有樣本之間的距離，因此s次投影需要的時間為Θ(sN2)。

FastMap算法可以概括為：在整個樣本集上找到m個樞軸，投影、修改所有樣本之間的距離。

其中，UF=(e1,e2,…,em)，e1,e2,…,em為各樞軸上的單位向量，彼此正交；r取決于各樞軸的起點坐標。

式（9）表明，尋找樞軸的本質是得到m個彼此正交的單位向量。根據文獻[17]，尋找樞軸的范圍可以從整個樣本集縮小到一個內在維數等于m的子集上。

其次，根據式（7），在計算坐標時，只用到了各點與樞軸點之間的距離。換句話說，每次投影后，修改所有樣本間的距離是不必要的。因此，如果能事先確定各樞軸，可以大大減少距離的計算。

最后，根據式（9），YFastMap是X的像矩陣，因此對YFastMap進行PCA，其結果YiFastMap等于CMDS坐標。

基于上述討論，文獻[17]提出了iFastMap算法。首先確定一個內在維數等于m的樣本子集（稱之為樞軸子集），然后在這個子集上確定m個樞軸，將各樣本投影到這些樞軸上，最后對投影結果進行PCA。

算法2iFastMap算法

輸入：距離矩陣D=(dij),內在維數m。

輸出：iFastMap坐標。

1.利用子集選擇算法確定n,取前n個樣本構成樞軸子集；

2.在這個樣本子集上運行FastMap算法

2.1初始化

令k=0；，i,j∈{1,2,…,n}

2.2重復以下投影過程m次；

2.2.1k=k+1；

3.對樞軸子集以外的樣本執行以下操作

3.1k=0；

3.2重復以下投影過程m次；

3.2.1k=k+1；

3.2.3利用式（8）計算i與{1,2,…,n}中樣本點之間的距離；

4.對YFastMap進行主成分分析，返回其結果YiFastMap。

由于每次投影只需計算各點與樞軸點之間的距離，當m?N時，iFastMap算法的運算時間為Θ(m2N)，它具有線性時間復雜度。

3.2　dcMDS算法

在iFastMap算法中，選定的子集與樣本上其他所有點的距離必須已知。但在實際應用中，經常出現只知道局部距離信息的情形。下面討論如何基于局部信息求出整個樣本集的CMDS坐標。

每個局部對應一個樣本子集。根據前面的討論，如果對每個子集直接用CMDS求解，會得到各子集原始坐標矩陣在不同保距變換下的像矩陣。因此，必須將這些像矩陣整合成同一保距變換下的像矩陣。

首先考慮兩個樣本子集的整合。假定子集A在下的像矩陣為YA，子集B在下的像矩陣為ZB。A、B的像矩陣能夠整合的前提是二者之間存在一個內在維數等于m的交集，用C表示。根據文獻[18]，C在、下的像矩陣YC、ZC之間存在保距變換，滿足：

其中，NC表示C中樣本個數。因為C的內在維數等于m，所以可以利用線性回歸求出。

文獻[18]進一步指出，如果將作用于YA，有：

這表明，利用，把A在保距變換下的像矩陣YA轉換成在下的像矩陣，使得A與B具有相同的保距變換，即A向B“對齊”。

對多個子集的情形：首先，從每個子集中隨機抽取部分點，構成一個內在維數等于m的“基準”子集；然后，將各個“基準”子集合并成一個“基準”集，并求出“基準”集的CMDS坐標；最后，將各子集向“基準”集“對齊”。

在像矩陣對齊的基礎上，結合分而治之策略[22-23]，作者提出了dcMDS算法[18]。dcMDS算法包含兩個過程：首先將樣本集自上而下逐級劃分，即將樣本集分成若干子集，如果子集的規模仍很大，則將子集進一步劃分，直到每個子集足夠小。第二個過程是自下而上逐級求子集的解。最底層的子集采用CMDS直接求解；底層子集的解逐級整合，直到得到整個樣本集的像矩陣。根據定理1，對該像矩陣進行PCA，其結果YdcMDS與CMDS坐標完全一致。

根據文獻[18]，當m?N時，dcMDS的運算時間為Θ(NlgN)。

3.3　iLMDS算法

在LMDS中，首先選取部分點作為標志點，然后計算標志點的CMDS坐標，最后給出其余點的歐氏坐標。LMDS算法歸納如下。

算法3LMDS算法

輸入：距離矩陣D=(dij),內在維數m。

輸出：LMDS坐標。

1.不失一般性，取前n個點作為標志點；

2.利用式（2）求標志點的中心化內積矩陣Bn；

令m1=rank(Bn)，則Bn有m1個正特征值對應的單位正交特征向量為

3.計算LMDS坐標：

由于將特征值的計算僅限于標志點，而在計算各點坐標時采用的是線性公式，當m?N時，LMDS具有線性時間復雜度Θ(N)。

根據文獻[15]，LMDS坐標是各樣本在標志點主軸向量上的投影，即：

從上式可以看出，如果m1=m，YLMDS是X的像矩陣。根據定理1，只需對YLMDS進行PCA，其結果YiLMDS等于CMDS坐標。

在LMDS算法中，對標志點的選取，只強調了n至少取為m+1，這不能確保m1=m。其實，只需采用子集選擇算法，就可以確保標志點集的內在維數等于m，從而得到CMDS的一致解。

基于上述分析，本文提出iLMDS算法，它由以下步驟組成：

（1）利用子集選擇算法確定n值，采用前n個樣本構成標志點集；

（2）調用LMDS算法求樣本集的CMDS坐標矩陣YLMDS；

（3）對YLMDS進行PCA，其結果YiLMDS等于樣本集的CMDS坐標，即：

4　實驗

實驗采用Acer TM4330筆記本電腦，CPU為雙核1.6 GHz，內存為2 GB；采用Matlab 2010a編程。針對USPS[24]和UCI[25]上的多組數據進行了實驗。原始數據對應樣本的各種屬性，換句話說，它們都屬于原始坐標，因此實驗之前，先把它們轉換成樣本間的歐氏距離。每個實驗都與CMDS算法進行對比，而CMDS算法對計算機的CPU、內存要求都很高，因此對于數據規模超過4 000的數據，都只選擇其前4 000個樣本。第1個實驗針對一組基準數據集比較了快速多維標度算法解與CMDS解的差異；第2個實驗比較了它們在該基準數據集上的運算時間；第3個實驗研究了它們的運算時間隨樣本個數的變化規律。

4.1　快速多維標度算法與CMDS算法解的一致性

為了驗證快速算法與CMDS算法解的一致性，本文采用以下兩個指標衡量快速算法與CMDS算法結果的差異。

（1）“標準化”誤差

（2）stress值

下式定義的stress[26]值常用來衡量多維標度算法的性能。

其中，?kl表示依據某快速算法坐標得到的歐氏距離；dkl表示依據CMDS坐標得到的歐氏距離。stress值衡量了這兩種距離的差，并用CMDS距離的和進行歸一化。

從表2、表3中可以看出，各種快速算法的“標準化”誤差和stress值都幾乎為0，表明快速算法能得到與CMDS算法一致的解。

Table 2　Normalized error of algorithms表2　各種算法的“標準化”誤差

Table 3　stressof algorithms表3　各種算法的stress值

4.2　各種算法的時間

Table 4　Running time of algorithms表4　各種算法的運算時間

各種算法的運算時間如表4所示。可以看出，3種快速算法在速度上都有了明顯改善。對大多數樣本集，它們的時間都在1 s以內。對usps這樣的數據集，由于樣本間的相關性較強，計算標志點集（樞軸點集）的時間增大，影響了速度，但即便如此，它們仍比CMDS算法快。

4.3　運算時間與樣本個數的關系

為了定量分析各種算法運算時間與樣本個數的關系，本文利用Matlab生成40組20維的隨機數，各組隨機數的個數為N=100,200,…,4 000。各算法的運算時間如圖1所示。從圖中可以看出，3種快速算法的速度較CMDS有了很大的提高。

Fig.1　Time of algorithms vs.number of samples圖1　各種算法運算時間與樣本個數的關系

進一步的分析表明：CMDS的運算時間在2.5×10-8×N3和 9.2×10-9×N3之間；iLMDS在N/140 000和N/240 000之間；iFastMap在N20 000和N/35 000之間；dcMDS在2.5×10-5×NlgN和1.5×10-5×NlgN之間。實驗結果驗證了前面的分析，即CMDS、iLMDS、iFastMap、dcMDS的時間復雜度分別為Θ(N3)、Θ(N)、Θ(N)、Θ(NlgN)。

5　討論

本文回顧了iFastMap和dcMDS算法，并提出了iLMDS算法。與CMDS算法不同的是，它們都是將樣本投影到一個子集上，從而提高了速度。

下面首先對比這3種算法的特點，然后討論各種算法中子集的選擇。

（1）3種算法的特點分析

在CMDS實際問題中，樣本集的距離信息可能以多種方式呈現，如圖2所示。圖（a）各樣本間的距離已知；圖（b）存在一個特殊的子集，該子集與樣本集上各點的距離已知；圖（c）只知道相鄰樣本之間的距離。

Fig.2　Three kinds of given distances圖2　3種距離情形

顯然，CMDS算法只能求解第一種情形；iLMDS和iFastMap算法可用于前兩種情形；dcMDS算法可用于任何一種情形。由此可見，快速算法拓展了多維標度算法的應用領域。

另外，CMDS是一種“批量”算法，即如果樣本集新增了樣本，必須重新計算，才能得到各點的坐標；而iFastMap和iLMDS都屬于“增量”式算法，可以在求出標志點集（樞軸子集）后，逐一計算新樣本的坐標。dcMDS算法也可以通過整合計算新樣本的坐標。

（2）子集的選擇

快速算法都涉及子集的選擇。確定iFastMap中的樞軸子集、iLMDS中的標志點集以及dcMDS中各“基準”子集，其判斷準則都是一樣的，即它們的內在維數等于m。

如前所述，在dcMDS算法中，對樣本集采取自上而下逐級劃分，直至每個子集足夠小。在實際應用中，對規模小于4m的子集停止劃分，以確保從該子集中能抽取出內在維數等于m的“基準”子集。

對樞軸子集、標志點集以及“基準”子集的抽取，都采用子集選擇算法。表5列出了這些子集的n值，其中dcMDS給出的是所有“基準”子集的最大n值。

從表5中可以看出，對多數樣本集，隨機選取m+1個樣本所得到的樣本子集，其內在維數等于m。

（3）未來的工作

CMDS是一種常用的數據降維和可視化方法，應用領域非常廣泛。本文側重闡述3種快速算法的基本理論，從理論上闡明它們與CMDS解的一致性。為了驗證它們的高效性以及與CMDS解的一致性，所有的實驗都與CMDS算法進行對比，這限制了實驗中樣本集的規模。將這些算法應用于大數據，是下一階段的工作。

Table 5　Size of chosen subsets表5　各算法選取的子集大小

6　結束語

本文提出了一種新的快速多維標度算法，并介紹了原來提出的另兩種快速多維標度算法。這3種方法都能得到與CMDS一致的解，而在速度上較CMDS有很大提高。下一階段的工作是考慮將這些算法用于解決數據降維和可視化問題。

[1]Mayer-Sch?nberger V,Cukier K.Big data:a revolution that will transform how we live,work,and think[M].New York:Houghton Mifflin Harcourt,2013:6-15.

[2]Verleysen M,Fran?ois D.The curse of dimensionality in data mining and time series prediction[C]//LNCS 3512:Proceedings of the 8th International Work-Conference on Artificial Neural Networks,Barcelona,Jun 8-10,2005.Berlin,Heidelberg:Springer,2005:758-770.

[3]Cox T F,Cox Michael A A.Multidimensional scaling[M].Boca Raton:Chapman&Hall/CRC,2001:38-44.

[4]Wu Chunqing,Ren Peige,Wang Xiaofeng.Survey on semanticbased organization and search technologies for network big data[J].Chinese Journal of Computers,2015,38(1):1-17.

[5]Ingram S,Munzner T.Dimensionality reduction for documents with nearest neighbor queries[J].Neurocomputing,2015,150(1):557-569.

[6]Patel A P,Tirosh I,Trombetta J J,et al.Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma[J].Science,2014,344(6190):1396-1401.

[7]Wang Shaowei,Zhuo Zhizheng,Yang Hongyu,et al.An approach to facial expression recognition integrating radial basis function kernel and multidimensional scaling analysis[J].Soft Computing,2014,18(7):1363-1371.

[8]Vogelstein J T,Park Y,Ohyama T,et al.Discovery of brainwide neural-behavioral maps via multiscale unsupervised structure learning[J].Science,2014,344(6182):386-392.

[9]Chalmers M.Alinear iteration time layout algorithm for visualising high-dimensional data[C]//Proceedings of the 7th IEEE Visualization Conference,San Francisco,Oct 27-Nov 1,1996.Piscataway:IEEE,1996:127-131.

[10]Morrison A,Ross G,Chalmers M.Fast multidimensional scaling through sampling,springs and interpolation[J].Information Visualization,2003,2(1):68-77.

[11]Williams M,Munzner T.Steerable,progressive multidimensional scaling[C]//Proceedings of the 10th IEEE Symposium on Information Visualization,Austin,Oct 10-12,2004.Washington:IEEE Computer Society,2004:57-64.

[12]Faloutsos C,Lin K.FastMap:a fast algorithm for indexing,data-mining and visualization of traditional and multimedia datasets[C]//Proceedings of the 1995 ACM SIGMOD International Conference on Management of Data,San Jose,May 22-25,1995.New York:ACM,1995:163-174.

[13]An Jiyuan,Yu J X,Ratanamahatana C A,et al.A dimensionality reduction algorithm and its application for interactive visualization[J].Journal of Visual Languages&Computing,2007,18(1):48-70.

[14]Wang J T L,Wang Xiong,Shasha D,et al.MetricMap:an embedding technique for processing distance-based queries in metric spaces[J].IEEE Transactions on Systems,Man,and Cybernetics:Part B Cybernetics,2005,35(5):973-987.

[15]de Silva V,Tenenbaum J B.Sparse multidimensional scaling using landmark points[R].Palo Alto:Stanford University,2004.

[16]Platt J C.FastMap,MetricMap,and Landmark MDS are all Nystr?m algorithms[C]//Proceedings of the 10th International Workshop on Artificial Intelligence and Statistics,Bridgetown,Jan 6-8,2005:261-268.

[17]Qu Taiguo,Cai Zixing.An improved FastMap algorithm[J].Journal of Nanjing University:Natural Sciences,2016,52(4):682-692.

[18]Qu Taiguo,Cai Zixing.A divide-and-conquer based multidimensional scaling algorithm[J].Pattern Recognition andArtificial Intelligence,2014,27(11):961-969.

[19]Zhang Runchu.Multivariate statistical analysis[M].Beijing:Science Press,2006.

[20]You Chengye.Analytic geometry[M].Beijing:Peking University Press,2004.

[21]Young G,Householder A S.Discussion of a set of points in terms of their mutual distances[J].Psychometrika,1938,3(1):19-22.

[22]Tian Wenbiao,Rui Guosheng,Kang Jian.Divide and conquer method for sparsity estimation within compressed sensing framework[J].Electronics Letters,2014,50(9):677-678.

[23]Li Shengguo,Gu Ming,Cheng Lizhi,et al.An accelerated divide-and-conquer algorithm for the bidiagonal SVD problem[J].SIAM Journal on Matrix Analysis and Applications,2014,35(3):1038-1057.

[24]Hull J J.A database for handwritten text recognition research[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1994,16(5):550-554.

[25]Lichman M.UCI machine learning repository[EB/OL].(2013).[2016-05-03].http://archive.ics.uci.edu/ml.

[26]Kruskal J B.Nonmetric multidimensional scaling:a numerical method[J].Psychometrika,1964,29(2):115-129.

附中文參考文獻：

[4]吳純青,任沛閣,王小峰.基于語義的網絡大數據組織與搜索[J].計算機學報,2015,38(1):1-17.

[17]屈太國,蔡自興.一種改進的FastMap算法[J].南京大學學報:自然科學,2016,52(4):682-692.

[18]屈太國,蔡自興.基于分而治之的多維標度算法[J].模式識別與人工智能,2014,27(11):961-969.

[19]張潤楚.多元統計分析[M].北京:科學出版社,2006.

[20]尤承業.解析幾何[M].北京:北京大學出版社,2004.