數據降維之多維縮放MDS（Multiple Dimensional Scaling）

網上看到關於數據降維的文章不少，介紹MDS的卻極少，遂決定寫一寫。

考慮一個這樣的問題。我們有n個樣本，每個樣本維度為m。我們的目標是用不同的新的k維向量（k<<m）替代原來的n個m維向量，使得在新的低維空間中，所有樣本相互之間的距離等於（或最大程度接近）原空間中的距離（默認歐氏距離）。

舉個栗子：原來有3個4維樣本(1,0,0,3),(8,0,0,5),(2,0,0,4)，顯然我們可以用三個新的二維樣本(1,3),(8,5),(2,4)來保持維度變小並相互之間距離不變。

那么問題來了，如果不是這么明顯的數據該如何來處理？降維后的距離一定會相等嗎？

MDS算法給出了在給定k值條件下的最優解決方案。

首先我們計算所有原空間中樣本相互之間的距離平方矩陣Dist[][]，顯然這是一非負對稱實數矩陣。至此，其實我們要維護的就是Dist不變，與原樣本已經無關了。

接下來我們要根據Dist推算出目標降維后內積矩陣B，B[i][j]就是降維后第i,j個向量的內積。關於推導過程可以看相關書籍，這里給出一個優美的結論。

　　B[i][j]=-0.5(Dist[i][j] - avg(Disti[i]) - avg(Distj[j]) + avg_Dist)

有了B，只需要對B分解成B=U*U^T的形式就達到我們的目標了。

對B做特征分解（奇異分解也一樣），B=V*diag*V^T。

我們可以取最大的k個特征值及其對應的特征向量構成diag_k和V_k。

此時U=V_k*diag_k^0,5就是我們降維后的n個行向量組成的矩陣了。

如果還有疑惑，下面的代碼運行試試就明白了。

 

 

召喚算法君：

 

import numpy as np

# run this to get a test matrix
# A = np.random.randint(1,100,(5,20))
# np.save('mat.npy', A)
# exit()

A = np.load('mat.npy')

n,m = A.shape
Dist = np.zeros((n,n))
B = np.zeros((n,n))
for i in range(n):
    for j in range(n):
        Dist[i][j] = sum((ix-jx)**2 for ix,jx in zip(A[i], A[j]))

disti2 = np.array([0]*n)
distj2 = np.array([0]*n)

for x in range(n):
    disti2[x] = np.mean([Dist[x][j] for j in range(n)])
    distj2[x] = np.mean([Dist[i][x] for i in range(n)])

distij2 = np.mean([Dist[i][j] for i in range(n) for j in range(n)])

for i in range(n):
    for j in range(n):
        B[i][j] = -0.5*(Dist[i][j] - disti2[i] - distj2[j] + distij2)

w,v = np.linalg.eig(B)

v=v.transpose()

U = [{'eVal':w[i], 'eVec':v[i]} for i in range(n)]

U.sort(key = lambda obj:obj.get('eVal'), reverse = True)
k=4
w=np.array([0]*k)
v=np.zeros((k,n))

for i in range(k):
    w[i] = U[i].get('eVal')**0.5
    v[i] = U[i].get('eVec')

ans = np.dot(v.transpose(), np.diag(w))

ans_dist = np.zeros((n,n))
for i in range(n):
    ans_str=""
    for j in range(n):
        ans_dist[i][j] = sum((ix-jx)**2 for ix,jx in zip(ans[i], ans[j]))

print("Orign dis[][] is :")
print Dist
print("MDS dis[][] is :")
print(ans_dist)