詳解視覺化利器 t-SNE 演算法：數無形時少直覺

本文介紹了 T 分佈隨機近鄰嵌入演算法，即一種十分強大的高維資料降維方法。我們將先簡介該演算法的基本概念與直觀性理解，再從詳細分析與實現該降維方法，最後我們會介紹使用該演算法執行視覺化的結果。

T 分佈隨機近鄰嵌入（T-Distribution Stochastic Neighbour Embedding）是一種用於降維的機器學習方法，它能幫我們識別相關聯的模式。t-SNE 主要的優勢就是保持區域性結構的能力。這意味著高維資料空間中距離相近的點投影到低維中仍然相近。t-SNE 同樣能生成漂亮的視覺化。

當構建一個預測模型時，第一步一般都需要理解資料。雖然搜尋原始資料並計算一些基本的統計學數字特徵有助於理解它，但沒有什麼是可以和圖表視覺化展示更為直觀的。然而將高維資料擬合到一張簡單的圖表（降維）通常是非常困難的，這就正是 t-SNE 發揮作用的地方。

在本文中，我們將探討 t-SNE 的原理，以及 t-SNE 將如何有助於我們視覺化資料。

t-SNE 演算法概念

這篇文章主要是介紹如何使用 t-SNE 進行視覺化。雖然我們可以跳過這一章節而生成出漂亮的視覺化，但我們還是需要討論 t-SNE 演算法的基本原理。

t-SNE 演算法對每個資料點近鄰的分佈進行建模，其中近鄰是指相互靠近資料點的集合。在原始高維空間中，我們將高維空間建模為高斯分佈，而在二維輸出空間中，我們可以將其建模為 t 分佈。該過程的目標是找到將高維空間對映到二維空間的變換，並且最小化所有點在這兩個分佈之間的差距。與高斯分佈相比 t 分佈有較長的尾部，這有助於資料點在二維空間中更均勻地分佈。

控制擬合的主要引數為困惑度（Perplexity）。困惑度大致等價於在匹配每個點的原始和擬合分佈時考慮的最近鄰數，較低的困惑度意味著我們在匹配原分佈並擬合每一個資料點到目標分佈時只考慮最近的幾個最近鄰，而較高的困惑度意味著擁有較大的「全域性觀」。

因為分佈是基於距離的，所以所有的資料必須是數值型。我們應該將類別變數通過二值編碼或相似的方法轉化為數值型變數，並且歸一化資料也是也十分有效，因為歸一化資料後就不會出現變數的取值範圍相差過大。

T 分佈隨機近鄰嵌入演算法（t-SNE）

Jake Hoare 的部落格並沒有詳細解釋 t-SNE 的具體原理和推導過程，因此下面我們將基於 Geoffrey Hinton 在 2008 年提出的論文和 liam schoneveld 的推導與實現詳細介紹 t-SNE 演算法。如果讀者對這一章節不感興趣，也可以直接閱讀下一章節 Jake Hoare 在實踐中使用 t-SNE 進行資料視覺化。

• liam schoneveld 推導與實現地址：https://nlml.github.io/in-raw-numpy/in-raw-numpy-t-sne/
  
• 論文地址：http://www.jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf
  

因為 t-SNE 是基於隨機近鄰嵌入而實現的，所以首先我們需要理解隨機近鄰嵌入演算法。

隨機近鄰嵌入（SNE）

假設我們有資料集 X，它共有 N 個資料點。每一個資料點 x_i 的維度為 D，我們希望降低為 d 維。在一般用於視覺化的條件下，d 的取值為 2，即在平面上表示出所有資料。

SNE 通過將資料點間的歐幾里德距離轉化為條件概率而表徵相似性（下文用 p_j|i 表示）：

如果以資料點在 x_i 為中心的高斯分佈所佔的概率密度為標準選擇近鄰，那麼 p_j|i 就代表 x_i 將選擇 x_j 作為它的近鄰。對於相近的資料點，條件概率 p_j|i 是相對較高的，然而對於分離的資料點，p_j|i 幾乎是無窮小量（若高斯分佈的方差σ_i 選擇合理）。

其中σ_i 是以資料點 x_i 為均值的高斯分佈標準差，決定σ_i 值的方法將在本章後一部分討論。因為我們只對成對相似性的建模感興趣，所以可以令 p_i|i 的值為零。

現在引入矩陣 Y，Y 是 N*2 階矩陣，即輸入矩陣 X 的 2 維表徵。基於矩陣 Y，我們可以構建一個分佈 q，其形式與 p 類似。

對於高維資料點 x_i 和 x_j 在低維空間中的對映點 y_i 和 y_j，計算一個相似的條件概率 q_j|i 是可以實現的。我們將計算條件概率 q_i|j 中用到的高斯分佈的方差設定為 1/2。因此我們可以對對映的低維資料點 y_j 和 y_i 之間的相似度進行建模：

我們的總體目標是選擇 Y 中的一個資料點，然後其令條件概率分佈 q 近似於 p。這一步可以通過最小化兩個分佈之間的 KL 散度（損失函式）而實現，這一過程可以定義為：

因為我們希望能最小化該損失函式，所以我們可以使用梯度下降進行迭代更新，我們可能對如何迭代感興趣，但我們會在後文討論與實現。

使用 NumPy 構建歐幾里德距離矩陣

計算 p_i|j 和 q_i|j 的公式都存在負的歐幾里德距離平方，即-||x_i - x_j||^2，下面可以使用程式碼實現這一部分：

def neg_squared_euc_dists(X):
    """Compute matrix containing negative squared euclidean
    distance for all pairs of points in input matrix X

    # Arguments:
        X: matrix of size NxD
    # Returns:
        NxN matrix D, with entry D_ij = negative squared
        euclidean distance between rows X_i and X_j
    """
    # Math? See https://stackoverflow.com/questions/37009647
    sum_X = np.sum(np.square(X), 1)
    D = np.add(np.add(-2 * np.dot(X, X.T), sum_X).T, sum_X)
    return -D

為了更高的計算效率，該函式使用矩陣運算的方式定義，該函式將返回一個 N 階方陣，其中第 i 行第 j 列個元素為輸入點 x_i 和 x_j 之間的負歐幾里德距離平方。

使用過神經網路的讀者可能熟悉 exp(⋅)/∑exp(⋅) 這樣的表達形式，它是一種 softmax 函式，所以我們定義一個 softmax 函式：

def softmax(X, diag_zero=True):
    """Take softmax of each row of matrix X."""

    # Subtract max for numerical stability
    e_x = np.exp(X - np.max(X, axis=1).reshape([-1, 1]))

    # We usually want diagonal probailities to be 0.
    if diag_zero:
        np.fill_diagonal(e_x, 0.)

    # Add a tiny constant for stability of log we take later
    e_x = e_x + 1e-8  # numerical stability

    return e_x / e_x.sum(axis=1).reshape([-1, 1])

注意我們需要考慮 p_i|i=0 這個條件，所以我們可以替換指數負距離矩陣的對角元素為 0，即使用 np.fill_diagonal(e_x, 0.) 方法將 e_x 的對角線填充為 0。

將這兩個函式放在一起後，我們能構建一個函式給出矩陣 P，且元素 P(i,j) 為上式定義的 p_i|j：

def calc_prob_matrix(distances, sigmas=None):
    """Convert a distances matrix to a matrix of probabilities."""
    if sigmas is not None:
        two_sig_sq = 2. * np.square(sigmas.reshape((-1, 1)))
        return softmax(distances / two_sig_sq)
    else:
        return softmax(distances)

困惑度

在上面的程式碼段中，Sigmas 引數必須是長度為 N 的向量，且包含了每一個σ_i 的值，那麼我們如何取得這些σ_i 呢？這就是困惑度（perplexity）在 SNE 中的作用。條件概率矩陣 P 任意行的困惑度可以定義為：

其中 H(P_i) 為 P_i 的夏農熵，即表示式如下：

在 SNE 和 t-SNE 中，困惑度是我們設定的引數（通常為 5 到 50 間）。我們可以為矩陣 P 的每行設定一個σ_i，而該行的困惑度就等於我們設定的這個引數。直觀來說，如果概率分佈的熵較大，那麼其分佈的形狀就相對平坦，該分佈中每個元素的概率就更相近一些。

困惑度隨著熵增而變大，因此如果我們希望有更高的困惑度，那麼所有的 p_j|i（對於給定的 i）就會彼此更相近一些。換而言之，如果我們希望概率分佈 P_i 更加平坦，那麼我們就可以增大σ_i。我們配置的σ_i 越大，概率分佈中所有元素的出現概率就越接近於 1/N。實際上增大σ_i 會增加每個點的近鄰數，這就是為什麼我們常將困惑度引數大致等同於所需要的近鄰數量。

搜尋σ_i

為了確保矩陣 P 每一行的困惑度 Perp(P_i) 就等於我們所希望的值，我們可以簡單地執行一個二元搜尋以確定σ_i 能得到我們所希望的困惑度。這一搜尋十分簡單，因為困惑度 Perp(P_i) 是隨σ_i 增加而遞增的函式，下面是基本的二元搜尋函式：

def binary_search(eval_fn, target, tol=1e-10, max_iter=10000, 
                  lower=1e-20, upper=1000.):
    """Perform a binary search over input values to eval_fn.
    
    # Arguments
        eval_fn: Function that we are optimising over.
        target: Target value we want the function to output.
        tol: Float, once our guess is this close to target, stop.
        max_iter: Integer, maximum num. iterations to search for.
        lower: Float, lower bound of search range.
        upper: Float, upper bound of search range.
    # Returns:
        Float, best input value to function found during search.
    """
    for i in range(max_iter):
        guess = (lower + upper) / 2.
        val = eval_fn(guess)
        if val > target:
            upper = guess
        else:
            lower = guess
        if np.abs(val - target) <= tol:
            break
    return guess

為了找到期望的σ_i，我們需要將 eval_fn 傳遞到 binary_search 函式，並且將σ_i 作為它的引數而返回 P_i 的困惑度。

以下的 find_optimal_sigmas 函式確實是這樣做的以搜尋所有的σ_i，該函式需要採用負歐幾里德距離矩陣和目標困惑度作為輸入。距離矩陣的每一行對所有可能的σ_i 都會執行一個二元搜尋以找到能產生目標困惑度的最優σ。該函式最後將返回包含所有最優σ_i 的 NumPy 向量。

def calc_perplexity(prob_matrix):
    """Calculate the perplexity of each row 
    of a matrix of probabilities."""
    entropy = -np.sum(prob_matrix * np.log2(prob_matrix), 1)
    perplexity = 2 ** entropy
    return perplexity


def perplexity(distances, sigmas):
    """Wrapper function for quick calculation of 
    perplexity over a distance matrix."""
    return calc_perplexity(calc_prob_matrix(distances, sigmas))


def find_optimal_sigmas(distances, target_perplexity):
    """For each row of distances matrix, find sigma that results
    in target perplexity for that role."""
    sigmas = [] 
    # For each row of the matrix (each point in our dataset)
    for i in range(distances.shape[0]):
        # Make fn that returns perplexity of this row given sigma
        eval_fn = lambda sigma: \
            perplexity(distances[i:i+1, :], np.array(sigma))
        # Binary search over sigmas to achieve target perplexity
        correct_sigma = binary_search(eval_fn, target_perplexity)
        # Append the resulting sigma to our output array
        sigmas.append(correct_sigma)
    return np.array(sigmas)

對稱 SNE

現在估計 SNE 的所有條件都已經宣告瞭，我們能通過降低成本 C 對 Y 的梯度而收斂到一個良好的二維表徵 Y。因為 SNE 的梯度實現起來比較難，所以我們可以使用對稱 SNE，對稱 SNE 是 t-SNE 論文中一種替代方法。

在對稱 SNE 中，我們最小化 p_ij 和 q_ij 的聯合概率分佈與 p_i|j 和 q_i|j 的條件概率之間的 KL 散度，我們定義的聯合概率分佈 q_ij 為：

該表示式就如同我們前面定義的 softmax 函式，只不過分母中的求和是對整個矩陣進行的，而不是當前的行。為了避免涉及到 x 點的異常值，我們不是令 p_ij 服從相似的分佈，而是簡單地令 p_ij=(p_i|j+p_j|i)/2N。

我們可以簡單地編寫這些聯合概率分佈 q 和 p：

def q_joint(Y):
    """Given low-dimensional representations Y, compute
    matrix of joint probabilities with entries q_ij."""
    # Get the distances from every point to every other
    distances = neg_squared_euc_dists(Y)
    # Take the elementwise exponent
    exp_distances = np.exp(distances)
    # Fill diagonal with zeroes so q_ii = 0
    np.fill_diagonal(exp_distances, 0.)
    # Divide by the sum of the entire exponentiated matrix
    return exp_distances / np.sum(exp_distances), None


def p_conditional_to_joint(P):
    """Given conditional probabilities matrix P, return
    approximation of joint distribution probabilities."""
    return (P + P.T) / (2. * P.shape[0])

同樣可以定義 p_joint 函式輸入資料矩陣 X 並返回聯合概率 P 的矩陣，此外我們還能一同估計要求的σ_i 和條件概率矩陣：

def p_joint(X, target_perplexity):
    """Given a data matrix X, gives joint probabilities matrix.

    # Arguments
        X: Input data matrix.
    # Returns:
        P: Matrix with entries p_ij = joint probabilities.
    """
    # Get the negative euclidian distances matrix for our data
    distances = neg_squared_euc_dists(X)
    # Find optimal sigma for each row of this distances matrix
    sigmas = find_optimal_sigmas(distances, target_perplexity)
    # Calculate the probabilities based on these optimal sigmas
    p_conditional = calc_prob_matrix(distances, sigmas)
    # Go from conditional to joint probabilities matrix
    P = p_conditional_to_joint(p_conditional)
    return P

所以現在已經定義了聯合概率分佈 p 與 q，若我們計算了這兩個聯合分佈，那麼我們能使用以下梯度更新低維表徵 Y 的第 i 行：

在 Python 中，我們能使用以下函式估計梯度，即給定聯合概率矩陣 P、Q 和當前低維表徵 Y 估計梯度：

def symmetric_sne_grad(P, Q, Y, _):
    """Estimate the gradient of the cost with respect to Y"""
    pq_diff = P - Q  # NxN matrix
    pq_expanded = np.expand_dims(pq_diff, 2)  #NxNx1
    y_diffs = np.expand_dims(Y, 1) - np.expand_dims(Y, 0)  #NxNx2
    grad = 4. * (pq_expanded * y_diffs).sum(1)  #Nx2
    return grad

為了向量化變數，np.expand_dims 方法將十分有用，該函式最後返回的 grad 為 N*2 階矩陣，其中第 i 行為 dC/dy_i。一旦我們計算完梯度，那麼我們就能利用它執行梯度下降，即通過梯度下降迭代式更新 y_i。

估計對稱 SNE

前面已經定義了所有的估計對稱 SNE 所需要的函式，下面的訓練函式將使用梯度下降演算法迭代地計算與更新權重。

def estimate_sne(X, y, P, rng, num_iters, q_fn, grad_fn, learning_rate,
                 momentum, plot):
    """Estimates a SNE model.

    # Arguments
        X: Input data matrix.
        y: Class labels for that matrix.
        P: Matrix of joint probabilities.
        rng: np.random.RandomState().
        num_iters: Iterations to train for.
        q_fn: Function that takes Y and gives Q prob matrix.
        plot: How many times to plot during training.
    # Returns:
        Y: Matrix, low-dimensional representation of X.
    """

    # Initialise our 2D representation
    Y = rng.normal(0., 0.0001, [X.shape[0], 2])

    # Initialise past values (used for momentum)
    if momentum:
        Y_m2 = Y.copy()
        Y_m1 = Y.copy()

    # Start gradient descent loop
    for i in range(num_iters):

        # Get Q and distances (distances only used for t-SNE)
        Q, distances = q_fn(Y)
        # Estimate gradients with respect to Y
        grads = grad_fn(P, Q, Y, distances)

        # Update Y
        Y = Y - learning_rate * grads
        if momentum:  # Add momentum
            Y += momentum * (Y_m1 - Y_m2)
            # Update previous Y's for momentum
            Y_m2 = Y_m1.copy()
            Y_m1 = Y.copy()

        # Plot sometimes
        if plot and i % (num_iters / plot) == 0:
            categorical_scatter_2d(Y, y, alpha=1.0, ms=6,
                                   show=True, figsize=(9, 6))

    return Y

為了簡化表達，我們將使用 MNIST 資料集中標籤為 0、1 和 8 的 200 個資料點，該過程定義在 main() 函式中：

# Set global parameters
NUM_POINTS = 200            # Number of samples from MNIST
CLASSES_TO_USE = [0, 1, 8]  # MNIST classes to use
PERPLEXITY = 20
SEED = 1                    # Random seed
MOMENTUM = 0.9
LEARNING_RATE = 10.
NUM_ITERS = 500             # Num iterations to train for
TSNE = False                # If False, Symmetric SNE
NUM_PLOTS = 5               # Num. times to plot in training


def main():
    # numpy RandomState for reproducibility
    rng = np.random.RandomState(SEED)

    # Load the first NUM_POINTS 0's, 1's and 8's from MNIST
    X, y = load_mnist('datasets/',
                      digits_to_keep=CLASSES_TO_USE,
                      N=NUM_POINTS)

    # Obtain matrix of joint probabilities p_ij
    P = p_joint(X, PERPLEXITY)

    # Fit SNE or t-SNE
    Y = estimate_sne(X, y, P, rng,
             num_iters=NUM_ITERS,
             q_fn=q_tsne if TSNE else q_joint,
             grad_fn=tsne_grad if TSNE else symmetric_sne_grad,
             learning_rate=LEARNING_RATE,
             momentum=MOMENTUM,
             plot=NUM_PLOTS)

構建 t-SNE

前面我們已經分析了很多關於隨機近鄰嵌入的方法與概念，並推匯出了對稱 SNE，不過幸運的是對稱 SNE 擴充套件到 t-SNE 是非常簡單的。真正的區別僅僅是我們定義聯合概率分佈矩陣 Q 的方式，在 t-SNE 中，我們 q_ij 的定義方法可以變化為： 

上式通過假設 q_ij 服從自由度為 1 的學生 T 分佈（Student t-distribution）而推匯出來。Van der Maaten 和 Hinton 注意到該分佈有非常好的一個屬性，即計數器（numerator）對於較大距離在低維空間中具有反平方變化規律。本質上，這意味著該演算法對於低維對映的一般尺度具有不變性。因此，最優化對於相距較遠的點和相距較近的點都有相同的執行方式。

這就解決了所謂的「擁擠問題」，即當我們試圖將一個高維資料集表徵為 2 或 3 個維度時，很難將鄰近的資料點與中等距離的資料點區分開來，因為這些資料點都聚集在一塊區域。

我們能使用以下函式計算新的 q_ij：

def q_tsne(Y):
    """t-SNE: Given low-dimensional representations Y, compute
    matrix of joint probabilities with entries q_ij."""
    distances = neg_squared_euc_dists(Y)
    inv_distances = np.power(1. - distances, -1)
    np.fill_diagonal(inv_distances, 0.)
    return inv_distances / np.sum(inv_distances), inv_distances

注意我們使用 1. - distances 代替 1. + distances，該距離函式將返回一個負的距離。現在剩下的就是重新估計損失函式對 Y 的梯度，t-SNE 論文中推導該梯度的表示式為：

同樣，我們很容易按照計算對稱 SNE 梯度的方式構建 t-SNE 的梯度計算方式：

def tsne_grad(P, Q, Y, inv_distances):
    """Estimate the gradient of t-SNE cost with respect to Y."""
    pq_diff = P - Q
    pq_expanded = np.expand_dims(pq_diff, 2)
    y_diffs = np.expand_dims(Y, 1) - np.expand_dims(Y, 0)

    # Expand our inv_distances matrix so can multiply by y_diffs
    distances_expanded = np.expand_dims(inv_distances, 2)

    # Multiply this by inverse distances matrix
    y_diffs_wt = y_diffs * distances_expanded

    # Multiply then sum over j's
    grad = 4. * (pq_expanded * y_diffs_wt).sum(1)
    return grad 

以上我們就完成了 t-SNE 的具體理解與實現，那麼該演算法在具體資料集中的視覺化效果是怎樣的呢？Jake Hoare 給出了實現視覺化的效果與對比。

t-SNE 視覺化

下面，我們將要展示 t-SNE 視覺化高維資料的結果，第一個資料集是基於物理特徵分類的 10 種不同葉片。這種情況下，t-SNE 需要使用 14 個數值變數作為輸入，其中就包括葉片的生長率和長寬比等。下圖展示了 2 維視覺化輸出，植物的種類（標籤）使用不同的顏色表達。

物種 Acer palmatum 的資料點在右上角形成了一個橙色叢集，這表明它的葉子和其他物種有很大的不同。該示例中類別通常會有很好的分組，相同物種的葉子（同一顏色的資料點）趨向於彼此靠緊聚集在一起。左下角有兩種顏色的資料點靠近在一起，說明這兩個物種的葉子形狀十分相似。

最近鄰準確度表明給定兩個隨機點，它們是相同物種的概率是多少。如果這些資料點完美地根據不同物種而分類，那麼準確度就會非常接近 100%，高準確度意味著資料能被幹淨地分為不同的叢集。

調整困惑度

下面，我們對可樂品牌做了類似的分析。為了演示困惑度（perplexity）的影響，我們首先需要將困惑度設定為較低的值 2，每個資料點的對映只考慮最近鄰。如下，我們將看到許多離散的小叢集，並且每一個叢集只有少量的資料點。

下面我們將 t-SNE 的困惑度設定為 100，我們可以看到資料點變得更加擴散，並且同一類之間的聯絡變弱。

在該案例中，可樂本身就要比樹葉更難分割，即使一類資料點某個品牌要更集中一些，但仍然沒有明確的邊界。

在實踐中，困惑度並沒以一個絕對的標準，不過一般選擇 5 到 50 之間會有比較好的結果。並且在這個範圍內，t-SNE 一般是比較魯棒的。

預測的解釋

度量資料點之間的角度或距離並不能推斷出任何資料上的具體或定量的資訊，所以 t-SNE 的預測更多的是用於視覺化資料。

在模型搭建前期直觀地挖掘資料模式有助於指導資料科學下一步程式。如果 t-SNE 能很好地分割資料點，那麼機器學習同樣也能找到一種將未知新資料投影到相應類別的好方法。給定一種預測演算法，我們就能實現很高的準確度。

上例中每一個類別都是孤立分類的，因此簡單的機器學習模型就能將該類別與其他類別分離開。但如果類別重疊，我們可能就要構建更精細的模型做出預測。如下所示，如果我們按照某個品牌的偏好從 1 到 5 進行分類，那麼類別可能更加離散、更難以預測，最近鄰精度也會相對較低。

對比 PCA

很自然我們就希望將 t-SNE 和其他降維演算法做比較。降維演算法中比較流行的是主成分分析法（PCA），PCA 會尋找能儘可能保留資料方差的新維度。有較大的方差的資料保留的資訊就越多，因為彼此不同的資料可以提供不同的資訊，所以我們最好是保留彼此分離的資料點，因為它們提供了較大的方差。

下面是採用 PCA 演算法對上文的樹葉類別進行降維的結果，我們看到雖然左側的橙色是分離的，但其它類別都有點混合在一起。這是因為 PCA 並不關心區域性的最近鄰關係，此外 PCA 是一種線性方法，所以它表徵非線性關係可能效果並不是很好。不過 PCA 演算法在壓縮資料為更小的特徵向量而投入到預測演算法中有很好地表現，這一點它要比 t-SNE 效果更好。

結語

t-SNE 是一種視覺化高維資料的優秀演算法，它經常要比其它降維演算法生成更具特點的視覺化結果。在資料分析中，獲得資料的先驗知識總是很重要的，正如華羅庚先生說過：數無形時少直覺，形少數時難入微，我們只有先理解了資料的大概分佈，然後再能選擇具體的演算法對這些資料進一步分析。數形結合百般好，隔離分家萬事休，也許高維資料的視覺化與機器學習演算法的結合才是資料分析的正確開啟方式。