算法

如前所述，t-SNE采用高維數據集，並將其簡化為保留了許多原始信息的低維圖。

假設我們有一個由3個不同的類組成的數據集。

我們希望將2D圖簡化為1D圖，同時保持群集之間的清晰邊界。

回想一下，僅將數據投影到軸上是降低維數的一種較差的方法，因為我們會丟失大量信息。

相反，我們可以使用降維技術(提示：t-SNE)來實現所需的功能。 t-SNE算法的第一步涉及測量從一個點到另一個點的距離。而不是直接處理距離，我們將它們映射到概率分布。

在分布中，相對於當前點距離最小的點的可能性很高，而遠離當前點的點的可能性很低。

再看一下2D圖，請注意藍色簇比綠色簇更分散。如果我們不解決這種尺度差異，則綠點的可能性將大於藍點的可能性。為了說明這個事實，我們用可能性之和做歸一化。

因此，盡管兩點之間的絕對距離不同，但它們被認為是相似的。

讓我們嘗試將這些概念與基礎理論聯係起來。在數學上，我們將正態分布的方程寫為如下形式：

如果我們將所有項都做指數計算，並使用另一個點代替均值，同時用求和規歸一化解決前麵討論的尺度問題，得到如下公式（參考論文)。

接下來，我們考慮降維到低維空間的情形。首先，我們創建一個n_samples x  n_components的矩陣(在這種情況下為9×1)並用隨機值(即位置)填充。

如果我們對上述情況采取類似的方法(測量點之間的距離並將其映射到概率分布)，則可以得到以下等式。

請注意，就像以前一樣，我們采用正態分布方程式，將所有內容放在前麵，使用其他點代替均值，然後通過除以所有其他點的似然之和來解決尺度問題(這裏忽略了標準差)。

如果我們能使降維後特征空間中的點的概率分布近似於原始特征空間中的點的概率分布，則可以得到定義良好的聚類。

為此，我們使用了稱為Kullback-Leiber散度的東西。 KL散度是一個概率分布與另一個概率分布之間差異的度量。

KL散度的值越小，兩個分布之間的距離越近。 KL散度為0表示所討論的兩個分布是相同的。

回想一下在線性回歸的情況下，我們如何通過使用梯度下降來最小化損失函數(即均方誤差)來確定最佳擬合曲線。在t-SNE中，我們同樣使用梯度下降法將所有數據點上Kullback-Leiber散度的總和最小化。

針對每個點取成本函數的偏導數，以便提供每次更新的方向。

t-SNE的Python代碼

在接下來的部分中，我將嘗試將算法和相關的數學方程式實現為Python代碼。為了完成該過程，我從scikit-learn 源代碼的TSNE類 中借鑒了一些東西。

首先，我們將導入以下庫並設置一些繪圖屬性，這些屬性將在我們繪製數據時發揮作用。

import numpy as np
from sklearn.datasets import load_digits
from scipy.spatial.distance import pdist
from sklearn.manifold.t_sne import _joint_probabilities
from scipy import linalg
from sklearn.metrics import pairwise_distances
from scipy.spatial.distance import squareform
from sklearn.manifold import TSNE
from matplotlib import pyplot as plt
import seaborn as sns
sns.set(rc={'figure.figsize':(11.7,8.27)})
palette = sns.color_palette("bright", 10)

在此示例中，我們將使用手寫數字數據。scikit-learn庫提供了一種將它們導入我們程序的方法。

X, y = load_digits(return_X_y=True)

人能夠理解或者可見的最多維度是3，所以這裏為t-SNE選擇2或3作為組件數(也就是聚類數量)。另一方麵，perplexity與算法中使用的最近鄰居的數量有關。不同的perplexity可能會導致最終結果發生巨大變化。在我們的情況下，我們將其設置為scitkit-learn默認值t-SNE(30)的實現。注：根據numpy文檔，機器精度(epsilon)是最小的可表示正數，因此1.0 + eps != 1.0。換句話說，由於缺少必要的位，因此計算機無法操縱機器精度(epsilon)以下的任何數字。正如我們將看到的，使用np.maximum檢查矩陣中的值是否小於機器的epsilon並在出現時替換它們。

MACHINE_EPSILON = np.finfo(np.double).eps
n_components = 2
perplexity = 30

接下來，我們定義fit功能。我們調用fit做數據轉換。

def fit(X):
    n_samples = X.shape[0]
    # Compute euclidean distance
    distances = pairwise_distances(X, metric='euclidean', squared=True)

    # Compute joint probabilities p_ij from distances.
    P = _joint_probabilities(distances=distances, desired_perplexity=perplexity, verbose=False)

    # The embedding is initialized with iid samples from Gaussians with standard deviation 1e-4.
    X_embedded = 1e-4 * np.random.mtrand._rand.randn(n_samples, n_components).astype(np.float32)

    # degrees_of_freedom = n_components - 1 comes from
    # "Learning a Parametric Embedding by Preserving Local Structure"
    # Laurens van der Maaten, 2009.
    degrees_of_freedom = max(n_components - 1, 1)

    return _tsne(P, degrees_of_freedom, n_samples, X_embedded=X_embedded)

fit函數實現了很多功能，下麵分解介紹。

1.將樣本數存儲在變量中，以便後續使用。

2.計算每個數據點之間的歐式距離。這對應於||xi — xj||^2。

3.將上一步中計算出的歐幾裏德距離作為參數傳遞給_join_probabilities函數，然後計算並返回一個矩陣p_ji(使用相同的方程式)。

4.使用標準差為1e-4的高斯分布隨機選擇值來創建縮減特征空間。

5.定義degrees_of_freedom。源代碼中有一條注釋， 在這篇論文中有解釋他們推理過程。基本上，從經驗上可以看出，當將degrees_of_freedom設置為組件數減去一個時，我們會得到更好的結果(粗體)。

6.最後，我們調用tsne函數，該函數的實現如下。

def _tsne(P, degrees_of_freedom, n_samples, X_embedded):
    params = X_embedded.ravel()
    obj_func = _kl_divergence

    params = _gradient_descent(obj_func, params, [P, degrees_of_freedom, n_samples, n_components])

    X_embedded = params.reshape(n_samples, n_components)
    return X_embedded

此函數實際上沒有太多內容。首先，我們使用np.ravel將向量展平為一維數組。

>>> x = np.array([[1, 2, 3], [4, 5, 6]])
>>> np.ravel(x)
array([1, 2, 3, 4, 5, 6])

然後我們使用梯度下降來最小化KL散度。完成後，我們將嵌入更改回2D數組並返回它。

接下來，看看計算細節。以下代碼塊負責基於kl發散和梯度計算誤差。

def _kl_divergence(params, P, degrees_of_freedom, n_samples, n_components):
    X_embedded = params.reshape(n_samples, n_components)

dist = pdist(X_embedded, “sqeuclidean”)
dist /= degrees_of_freedom
dist += 1.
dist **= (degrees_of_freedom + 1.0) / -2.0
Q = np.maximum(dist / (2.0 * np.sum(dist)), MACHINE_EPSILON)

# Kullback-Leibler divergence of P and Q
kl_divergence = 2.0 * np.dot(P, np.log(np.maximum(P, MACHINE_EPSILON) / Q))

# Gradient: dC/dY
grad = np.ndarray((n_samples, n_components), dtype=params.dtype)
PQd = squareform((P – Q) * dist)
for i in range(n_samples):
grad[i] = np.dot(np.ravel(PQd[i], order=’K’), X_embedded[i] – X_embedded)
grad = grad.ravel()
c = 2.0 * (degrees_of_freedom + 1.0) / degrees_of_freedom
grad *= c

return kl_divergence, grad

同樣，讓我們​​一步一步地講解代碼。

1.第一部分計算低維映射中各點的概率分布。

實際上使用的是上麵方程式的變體，其中包括自由度。

其中α表示Student-t分布的自由度數

2.計算KL散度(np.dot)。

3.計算梯度(偏導數)。dist實際上是yi — yj：

同樣，他們使用自由度以上方程的變化形式。

其中α表示Student-t分布的自由度數

梯度下降函數通過最小化KL散度來更新嵌入中的值。當梯度範數低於閾值或達到最大迭代次數而沒有任何進展時，我們提前停止。

def _gradient_descent(obj_func, p0, args, it=0, n_iter=1000,
                      n_iter_check=1, n_iter_without_progress=300,
                      momentum=0.8, learning_rate=200.0, min_gain=0.01,
                      min_grad_norm=1e-7):
    p = p0.copy().ravel()
    update = np.zeros_like(p)
    gains = np.ones_like(p)
    error = np.finfo(np.float).max
    best_error = np.finfo(np.float).max
    best_iter = i = it

    for i in range(it, n_iter):
        error, grad = obj_func(p, *args)
        grad_norm = linalg.norm(grad)
        inc = update * grad < 0.0
        dec = np.invert(inc)
        gains[inc] += 0.2
        gains[dec] *= 0.8
        np.clip(gains, min_gain, np.inf, out=gains)
        grad *= gains
        update = momentum * update - learning_rate * grad
        p += updateprint("[t-SNE] Iteration %d: error = %.7f,"
                      " gradient norm = %.7f"
                      % (i + 1, error, grad_norm))

        if error < best_error:
            best_error = error
            best_iter = i
        elif i - best_iter > n_iter_without_progress:
            break
        if grad_norm <= min_grad_norm:
            break

      return p

到這裏，已經做好了在數據上執行fit了。

X_embedded = fit(X)

如我們所見，該模型在根據像素位置分離不同數字方麵表現很好。

sns.scatterplot(X_embedded[:,0], X_embedded[:,1], hue=y, legend='full', palette=palette)

t-SNE的Sklearn代碼

下麵是scikit-learn的t-SNE實現。

tsne = TSNE()
X_embedded = tsne.fit_transform(X)

如我們所見，該模型成功地獲取了64維數據集，並將其投影到2維空間中，從而使相似的樣本聚在一起。

sns.scatterplot(X_embedded[:,0], X_embedded[:,1], hue=y, legend='full', palette=palette)