馬氏距離(Mahalanobis Distence) [python]

馬氏距離(Mahalanobis Distence)

是度量學習(metric learning)中一種常用的測度，所謂測度/距離函數/度量(metric)也就是定義一個空間中元素間距離的函數，所謂度量學習也叫做相似度學習。

---

什么是馬氏距離

似乎是一種更好度量相似度的方法。

馬氏距離是基於樣本分布的一種距離。

物理意義就是在規范化的主成分空間中的歐氏距離。

所謂規范化的主成分空間就是利用主成分分析對一些數據進行主成分分解。

再對所有主成分分解軸做歸一化，形成新的坐標軸。由這些坐標軸張成的空間就是規范化的主成分空間。

解釋第一句話: 馬氏距考慮到了樣本在不同類別的不同分布情況。第二句: 是基於歐氏距離的一種改進。后兩句針對高維線性分布的數據中各維度間非獨立同分布分布的問題，后文繼續探討。

距離公式:

• 單點長度:

  \(D_M(x) = \sqrt[]{(x-\mu)^{T} \Sigma^{-1}(x-\mu)}\)

• 兩點間長度:

  \(D_M(x, y) = \sqrt[]{(x-y)^{T} \Sigma_{x,y}^{-1}(x-y)}\)

  1. 其中\(\Sigma\)是多維隨機變量的協方差矩陣，若其是單位矩陣，則馬氏距離退化為歐式距，即各維度間獨立且同分布(因為只有同分布才會方差相同，對角線是相同的方差)
  2. (單點長度中的這個\(\mu\)可以看成是樣本點的中心，歐氏距離的一般定義是把原點看成中心點)

馬氏距離的作用(基於歐式距離)

歐氏距離對於不同的量綱一視同仁

舉個例子：在模式識別課上有個數據是學生的身高、鞋碼、體重來分類性別。三個數據的單位都不一樣，如果x1樣本和x2的鞋碼之間相差10(其他一樣)，x1和x3的體重相差10(其他一樣)，在歐式空間中就認為x1，x2的距離和x1，x3的距離是相同的，明顯x1和x2的性別差異就很大。

下面隨便舉一個例子，隨機生成一個維度間相關的分布

歸一化可以解決這個問題。歸一化之后維度的中心，就在坐標原點的地方了。

歸一化后的歐氏距離考慮方差的影響

• 舉個例子:一個樣本點分別到兩個類分布中心的距離相同(歐式)，可以通過觀察方差分布來分類。

• 判斷兩個點哪一個是該類的例子：歐氏距離對於方差視而不見，即使中心化之后兩個點到中心的相對距離的互相關系（誰近誰遠）是不會改變的。

  ​ 如上圖(模擬某一個類別的兩個屬性值的分布)，兩個點到其維度中心的距離都相同，但是明顯左側的紅點不屬於這一類。

  ​ 但是僅計算測試樣本(紅，黃兩點)與樣本中心的歐式距離來判斷，這兩個點會被分為同一類。

所以要本質上解決這個問題，就要針對主成分分析中的主成分來進行標准化。

馬氏距離——向量空間按照主成分旋轉后的歐氏距離

主成分分析：找到主成分方向(方差大的維度)，將整個樣本空間按照主成分方向旋轉，讓主成分方向作為新的軸，讓維度之間盡可能互相獨立。

馬氏距離：在主成分空間中，樣本到原點的歐式距離(因為PCA的時候已經中心化數據了，所以該類別樣本的中心就是原點)

計算步驟：

1. 找出主成分方向，中心化后的數據求協方差矩陣，經過特征值分解找出排序后的特征值矩陣U

2. U矩陣作用於數據矩陣，旋轉空間

3. 新數據空間標准化，減去均值除以標准差，標准化后，讓維度同分布，獨立

4. 計算新坐標的歐氏距離

   注意: 此時原點是類中心

可見現在黃色的點距離中心的距離比紅色的點更近了，紅點即為離群點。

可見在主成分空間能更好的看相似度。

數學推導

根據步驟:

假設數據沒有中心化后的\(X\)矩陣每一列是一個樣本，旋轉矩陣為\(U\)，新的坐標為:

\[\begin{aligned} F =& (F_1, F_2, ...,F_m) = U^T X,X\in \mathbb{R}^{n\times m} \\ \mu_F =& (\mu_1, \mu_2,..., \mu_m) \\ F - \mu_F =& U^T (X - \mu_X) \\ &\mu_X是np.mean(X, axis=1)\\ &這里其實挺亂的, 用到的就是新坐標和原始數據的關系U \end{aligned} \]

變換后每個維度都線性無關，各維度方差為協方差矩陣的特征值，滿足如下(最后要用到):

\(Cov_F = (F-\mu_F)(F-\mu_F)^T = U^T(X-\mu_X)(X-\mu_X)^TU\)
\(=U^TCov_XU\)

馬氏距離是旋轉變換加上縮放之后的歐氏距離(這里縮放其實就是對F在做一次標准化，使得維度同方差):

每一個維度的方差:\(\sigma_{Fi} = \sqrt{\lambda_{Fi}},Cov_F=diag(\lambda_{F1},..., \lambda_{Fn})\)

標准化之后(f是一個樣本，f_i是第i個維度):\(fi' = (\frac{f_i-\mu_{Fi}}{\sigma_{Fi}})\)

坐標到原點的歐氏距離(未開方):\(D_M=f_1'^2 + f_2'^2 + f_3'^2 + ... + f_n'^2\)

\(=(f_1-\mu_{F1}, f_2-\mu_{F2}, ...,f_n-\mu_{Fn})diag(\frac{1}{\lambda_1}, \frac{1}{\lambda_2}, ..., \frac{1}{\lambda_n})(f_1-\mu_{F1}, f_2-\mu_{F2}, ...,f_n-\mu_{Fn})^T\)

\(=(f-\mu_F)^T(U^TCov_XU)^{-1}(f-\mu_F),中間的對角陣是Cov_F^{-1}\)

\(= (x-\mu_X)^TUU^TCov_X^{-1}UU^T(x-\mu_X), 由前面x旋轉得到的f\)

\(= (x-\mu_X)^TCov_X^{-1}(x-\mu_X), U\)是單位正交的
最后在開平方就是馬氏距離

推導式中的U可以是一個非方陣，也就是用了PCA的方法，但是不影響結果，因為在中間步驟U會被消去。

馬氏距離的問題

• 協方差矩陣必須是滿秩，需要求逆，對稱矩陣不一定可逆哦，原因是其特征值會出現0，如果出現則可以考慮先PCA，這里PCA不會損失信息，原因上面提到了。
• 不能處理非線性流行(manifold)，流行問題還沒有了解。。
• 馬氏距離其實也是基於樣本的，自學習到一種空間變化，在主成分空間計算距離

附python作圖&實現代碼

class MahalanobisDistance():
    # 沒有考慮的很周全的代碼
    def __init__(self, X):
        """
        構造函數
        :param
            X:m * n維的np數據矩陣 每一行是一個sample 列是特征
        """
        self._PCA = None
        self._mean_x = np.mean(X, axis=0)
        mean_removed = X # - self._mean_x
        # cov = np.dot(mean_removed.T, mean_removed) / X.shape[0] # 計算協方差矩陣
        cov = np.cov(mean_removed, rowvar=False)
        # 判斷協方差矩陣是否可逆
        if np.linalg.det(cov) == 0.0:
            # 對數據做PCA 去掉特征值0的維度
            eig_val, eig_vec = np.linalg.eig(cov)
            e_val_index = np.argsort(-eig_val)  # 逆序排
            e_val_index = e_val_index[e_val_index > 1e-3]  # 需要特征值大於0的維度
            self._PCA = eig_vec[:, e_val_index]     # 降維矩陣 Z = XU
            PCA_X = np.dot(X, self._PCA)        # 降維
            self._mean_x = PCA_X.mean(axis=0) # 重新計算均值 去中心
            mean_removed = PCA_X # - self._mean_x
            # cov = np.dot(mean_removed.T, mean_removed) / PCA_X.shape[0] # 重新計算協方差矩陣
            cov = np.cov(mean_removed, rowvar=False)
            
        self._cov_i = np.linalg.inv(cov) # 協方差矩陣求逆
        
    def __call__(self, X, y=None):
        """
        計算x與y的馬氏距離 如果不傳入y則計算x到樣本中心的距離
        :param
            X:行向量/矩陣 樣本點特征維數必須和原始數據一樣
            y:行向量 樣本點特征維數必須和原始數據一樣
        :return
            distance 馬氏距離 如果出錯則返回-1
        """        
        # 不考慮出錯的情況 維度不符合
        if y is None:
            # 計算到樣本中心的距離
            y = self._mean_x
        X_data = X.copy()
        if self._PCA is not None:
            X_data = np.dot(X_data, self._PCA)    # 對數據降維
        X_data = X_data - y
        distance = np.dot(np.dot(X_data, self._cov_i), X_data.T)
        if len(X.shape) != 1:
            # X是一個矩陣
            distance = distance.diagonal()
            
        return np.sqrt(distance)

# m_distance = MahalanobisDistance(train_data)
# d = m_distance(x, y)

下面的作圖可看可不看，作為記錄

def plot_distri(x1, x2, rang, p1, p2):
    # 畫分布圖 p1 p2是測試點每一個維度的方差
    plt.figure()
    plt.scatter(x1, x2, s=10, marker='H', alpha=0.6)
    plt.scatter(p1[0], p2[0], s=20, color='r')
    plt.scatter(p1[1], p2[1], s=20, color='y')
    plt.xlim(rang)
    plt.ylim(rang)
    plt.axvline(x=0, color='k')
    plt.axhline(y=0, color='k')

def normalization(X):
    mean = X.mean(axis=1)
    std = X.std(axis=1)
    # print("mean", mean)
    return (X - np.tile(mean.reshape(-1, 1), X.shape[1])) / np.tile(std.reshape(-1, 1), X.shape[1]), std

x1 = np.random.normal(6, 3, 200)
x2 = x1 * 1.2 + 2.5 * np.random.randn(200)

rang =(-10, 20)
# 畫一個有維度相關性的圖
# 去中心化
x1_m = x1.mean()
x1_std = x1.std()
x2_m = x2.mean()
x2_std = x2.std()
x1s = (x1 - x1_m) / x1_std
x2s = (x2 - x2_m) / x2_std

px1 = np.array([-1, 1])     # 第二個點是 正例
px2 = np.array([1, 1])

pxx1 = px1 * x1_std + x1_m
pxx2 = px2 * x2_std + x2_m

plot_distri(x1, x2, rang, pxx1, pxx2)
plot_distri(x1s, x2s, (-3, 3), px1, px2)

X = np.vstack([x1s, x2s])

covX = np.cov(X, rowvar=True)
eigval, eigvec = np.linalg.eig(covX)
eig_arg = np.argsort(-eigval)
U = eigvec[:, eig_arg]
# 按照主成分方向旋轉
M = np.dot(U, X)
# 按照特征大小縮放
M, std = normalization(M)

M_p = np.dot(U, np.array([[-1, 1], [1, 1]]))
M_p[0, :] = M_p[0, :] / std[0]
M_p[1, :] = M_p[1, :] / std[1]

plot_distri(M[0, :], M[1, :], (-6, 6), M_p[:, 0], M_p[:, 1])
plt.show()