学习日记_241027_高维数据可视化，t-分布邻域嵌入（t-SNE）

学习日记，高维数据可视化，t-sne

行然梦实

1219人浏览 · 2024-10-28 21:43:18

行然梦实 · 2024-10-28 21:43:18 发布

前言

提醒：
文章内容为方便作者自己后日复习与查阅而进行的书写与发布，其中引用内容都会使用链接表明出处（如有侵权问题，请及时联系）。
其中内容多为一次书写，缺少检查与订正，如有问题或其他拓展及意见建议，欢迎评论区讨论交流。

相关链接：
t-SNE高维数据可视化（python）
t-SNE使用过程中的一些坑
 代码抄录源

代码抄录

import numpy as np 
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt
import time

def cal_pairwise_dist(x):
    sum_x=np.sum(np.square(x),1)
    dist=np.add(np.add(-2*np.dot(x,x.T),sum_x).T,sum_x)
    return dist
def cal_perplexity(dist,idx=0,beta=1.0):
    prob=np.exp(-dist*beta)
    prob[idx]=0
    sum_prob=np.sum(prob)
    if sum_prob<1e-12:
        prob=np.maximum(prob,1e-12)
        perp=-12
    else:
        perp=np.log(sum_prob) + beta*np.sum(dist*prob) / sum_prob
        prob=prob/sum_prob
    return perp,prob
def search_prob(x,tol=1e-5,perplexity=30.0):
    print("Computing pairwise distances...")
    (n,d)=x.shape
    dist=cal_pairwise_dist(x)
    dist[dist<0]=0
    pair_prob=np.zeros((n,n))
    beta=np.ones((n,1))
    base_perp=np.log(perplexity)
    for i in range(n):
        if i%500==0:
            print("Computing pair_prob for point %s of %s..." % (i,n))
        betamin=-np.inf
        betamax=np.inf
        perp,this_prob=cal_perplexity(dist[i],i,beta[i])
        perp_diff=perp-base_perp
        tries=0
        while np.abs(perp_diff)>tol and tries<50:
            if perp_diff>0:
                betamin=beta[i].copy()
                if betamax==np.inf or betamax==-np.inf:
                    beta[i]=beta[i]*2
                else:
                    beta[i]=(beta[i]+betamax)/2
            else:
                betamax=beta[i].copy()
                if betamin==np.inf or betamin==-np.inf:
                    beta[i]=beta[i]/2
                else:
                    beta[i]=(beta[i]+betamin)/2
            perp,this_prob=cal_perplexity(dist[i],i,beta[i])
            perp_diff=perp-base_perp
            tries=tries+1
        pair_prob[i,]=this_prob
    print("Mean value of sigma: %f" % np.mean(np.sqrt(1/(beta))))
    return pair_prob
def tsne(x,no_dims=2,perplexity=30.0,max_iter=1000):
    if isinstance(no_dims,float):
        print("Error: array x should have type float")
        return -1
    (n,d)=x.shape
    initial_momentum=0.5
    final_momentum=0.8
    eta=500
    min_gain=0.01
    y=np.random.randn(n,no_dims)
    dy=np.zeros((n,no_dims))
    iy=np.zeros((n,no_dims))
    gains=np.ones((n,no_dims))
    P=search_prob(x,1e-5,perplexity)
    P=P+np.transpose(P)
    P=P/np.sum(P)
    print("T_SNE DURING:%s" % time.process_time())
    P=P*4
    P=np.maximum(P,1e-12)
    for iter in range(max_iter):
        sum_y=np.sum(np.square(y),1)
        num=1/(1+np.add(np.add(-2*np.dot(y,y.T),sum_y).T,sum_y))
        num[range(n),range(n)]=0
        Q=num/np.sum(num)
        Q=np.maximum(Q,1e-12)
        PQ=P-Q
        for i in range(n):
            dy[i,:]=np.sum(np.tile(PQ[:,i]*num[:,i],(no_dims,1)).T*(y[i,:]-y),0)
        if iter<20:
            momentum=initial_momentum
        else:
            momentum=final_momentum
        gains=(gains+0.2)*((dy>0)!=(iy>0))+(gains*0.8)*((dy>0)==(iy>0))
        gains[gains<min_gain]=min_gain
        iy=momentum*iy-eta*(gains*dy)
        y=y+iy
        y=y-np.tile(np.mean(y,0),(n,1))
        if (iter+1)%100==0:
            C=np.sum(P*np.log(P/Q))
            print("Iteration",(iter+1),": error is",C)
            if (iter+1)!=100:
                ratio=C/oldC   
                print("ratio ",ratio)
                if ratio>=0.95:
                    break
            oldC=C
        if iter==100:
            P=P/4
    print("finished training ")
    return y

if __name__=="__main__":
    digits=load_digits()
    X=digits.data
    Y=digits.target
    data_2d=tsne(X,2)
    plt.scatter(data_2d[:,0],data_2d[:,1],c=Y)
    plt.show()

代码分析

数据集介绍

load_digits() 是 scikit-learn 库中的一个函数，用于加载手写数字数据集。这个数据集常用于图像分类和机器学习算法的测试和演示。

详细介绍

load_digits(): 加载手写数字数据集。
digits.data: 一个二维数组，形状为 (1797, 64)，每一行表示一个 8x8 的图像展开成的 64 维特征向量。
digits.images: 一个三维数组，形状为 (1797, 8, 8)，表示原始的 8x8 图像。
digits.target: 一个一维数组，长度为 1797，表示图像对应的数字标签（0 到 9）。
digits.target_names: 数组，表示可能的目标名称（数字 0 到 9）。

代码展示

from sklearn.datasets import load_digits
import matplotlib.pyplot as plt

# 加载数据集
digits = load_digits()

# 数据集的基本信息
print(f"数据集形状: {digits.data.shape}")
print(f"目标类别数量: {len(digits.target_names)}")

# 显示前5个数据点的图像和标签
fig, axes = plt.subplots(1, 5, figsize=(10, 3))
for ax, image, label in zip(axes, digits.images, digits.target):
    ax.set_axis_off()
    ax.imshow(image, cmap=plt.cm.gray_r, interpolation='nearest')
    ax.set_title(f'Target: {label}')
plt.show()

# 打印第一张图片的矩阵
print(digits.images[0])

运行结果为：

注意 $d i g i t s . d a t a . s ha p e$ d的形状为（1797，64），即数据集有1797条数据，每条数据为64个数（8*8矩阵摊平）的行向量。

t-SNE算法简介

t-SNE（t-分布邻域嵌入）是一种非线性降维技术，尤其适合用于高维数据的可视化。它通过将数据点嵌入到二维或三维空间，保留高维空间中相似点的距离关系。

数学表达

相似度定义
- 在高维空间中，对于两个数据点 $x_i$ 和 $x_j$ ，定义其相似度为条件概率 $p_{j|i}$ ，表示在给定 $x_i$ 的情况下选择 $x_j$ 的概率。通常使用高斯分布来计算：
  $p_{j|i} = \frac{\exp(-\|x_i - x_j\|^2 / 2\sigma_i^2)}{\sum_{k \neq i} \exp(-\|x_i - x_k\|^2 / 2\sigma_i^2)}$
- 最终相似度 $p_{ij}$ 是对称的，定义为：
  $p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}$
低维空间相似度
- 在低维空间中，使用 t-分布计算相似度 $q_{ij}$ ：
  $q_{ij} = \frac{(1 + \|y_i - y_j\|^2)^{-1}}{\sum_{k \neq l} (1 + \|y_k - y_l\|^2)^{-1}}$
Kullback-Leibler 散度
- t-SNE 通过最小化高维和低维分布之间的 Kullback-Leibler 散度来找到最优嵌入：
  $\text{KL}(P || Q) = \sum_{i \neq j} p_{ij} \log \frac{p_{ij}}{q_{ij}}$
梯度下降
- 使用梯度下降法最小化 KL 散度，更新低维空间坐标 $y_i$ 。

优势和应用

保留局部结构: t-SNE 能够很好地保留高维空间中局部结构的关系。
数据可视化: 被广泛用于图像、文本和基因数据的降维与可视化。

注意事项

计算复杂度高: t-SNE 对于大数据集可能较慢。
参数选择: 需要调节参数如perplexity来获得最佳结果。

cal_pairwise_dist(x)

def cal_pairwise_dist(x):
    sum_x=np.sum(np.square(x),1)
    dist=np.add(np.add(-2*np.dot(x,x.T),sum_x).T,sum_x)
    return  dist

代码展示的结果为：
（在 np.sum(np.square(x), 1) 中，1 表示沿着矩阵的行方向（即第2个维度）进行求和。）

sum_x=np.sum(np.square(x),1)：

dist=np.add(np.add(-2*np.dot(x,x.T),sum_x).T,sum_x)：

最终，得到的矩阵 $\mathbf{D}$ 的每个元素 $\mathbf{D}_{ij}$ 代表数据点 $x_i$ 和 $x_j$ 之间的欧几里得距离的平方。

cal_perplexity(dist,idx=0,beta=1.0)

def cal_perplexity(dist,idx=0,beta=1.0):
    prob=np.exp(-dist*beta)
    prob[idx]=0
    sum_prob=np.sum(prob)
    if sum_prob<1e-12:
        prob=np.maximum(prob,1e-12)
        perp=-12
    else:
        perp=np.log(sum_prob) + beta*np.sum(dist*prob) / sum_prob
        prob=prob/sum_prob
    return perp,prob

prob=np.exp(-dist*beta)
prob[idx]=0
sum_prob=np.sum(prob)
~~概率（probability）困惑度（perplexity）~~
计算未归一化的相似性概率:
对于给定的点 $i$ ，计算与其他点 $j$ 的未归一化相似性概率：
$\text{prob}_j = \exp(-\beta \cdot \text{dist}(x_i, x_j))$
这里， $\text{dist}(x_i, x_j)$ 是点 $i$ 和点 $j$ 之间的距离， $\beta$ 是控制高斯宽度的参数。
排除自身概率:
设置 $\text{prob}_i = 0$ 以避免自我相似性计算：
$\text{prob}_i = 0$
计算概率和:
计算所有未归一化概率的和：
$sum_prob = ∑ j prob j \text{sum\_prob} = \sum_{j} \text{prob}_j$
这个过程为每个点计算与其他点的相对相似性，除了自身。

    if sum_prob<1e-12:
        prob=np.maximum(prob,1e-12)
        perp=-12
    else:
        perp=np.log(sum_prob) + beta*np.sum(dist*prob) / sum_prob
        prob=prob/sum_prob

这段代码处理相似性概率的数学表达如下：
情况 1: $sum_prob < 1 × 1 0 − 12 \text{sum\_prob} < 1 \times 10^{-12}$

调整概率:
$\text{prob}_j = \max(\text{prob}_j, 1 \times 10^{-12})$

设定困惑度:
$\text{perp} = -12$

情况 2: $sum_prob ≥ 1 × 1 0 − 12 \text{sum\_prob} \geq 1 \times 10^{-12}$

计算困惑度:
$sum_prob ) + β ∑ j dist ( x i , x j ) ⋅ prob j sum_prob \text{perp} = \log(\text{sum\_prob}) + \frac{\beta \sum_{j} \text{dist}(x_i, x_j) \cdot \text{prob}_j}{\text{sum\_prob}}$

归一化概率:
$sum_prob \text{prob}_j = \frac{\text{prob}_j}{\text{sum\_prob}}$

解释

调整概率: 当 $sum_prob \text{sum\_prob}$ 非常小时，概率调整到一个最小值以避免数值不稳定性。

困惑度定义: 困惑度（perplexity）提供了一种衡量信息熵的方式。它反映了概率分布的复杂性。

归一化: 确保所有概率之和为1，从而得到有效的概率分布。

search_prob(x,tol=1e-5,perplexity=30.0)

def search_prob(x,tol=1e-5,perplexity=30.0):
    print("Computing pairwise distances...")
    (n,d)=x.shape
    dist=cal_pairwise_dist(x)
    dist[dist<0]=0
    pair_prob=np.zeros((n,n))
    beta=np.ones((n,1))
    base_perp=np.log(perplexity)
    for i in range(n):
        if i%500==0:
            print("Computing pair_prob for point %s of %s..." % (i,n))
        betamin=-np.inf
        betamax=np.inf
        perp,this_prob=cal_perplexity(dist[i],i,beta[i])
        perp_diff=perp-base_perp
        tries=0
        while np.abs(perp_diff)>tol and tries<50:
            if perp_diff>0:
                betamin=beta[i].copy()
                if betamax==np.inf or betamax==-np.inf:
                    beta[i]=beta[i]*2
                else:
                    beta[i]=(beta[i]+betamax)/2
            else:
                betamax=beta[i].copy()
                if betamin==np.inf or betamin==-np.inf:
                    beta[i]=beta[i]/2
                else:
                    beta[i]=(beta[i]+betamin)/2
            perp,this_prob=cal_perplexity(dist[i],i,beta[i])
            perp_diff=perp-base_perp
            tries=tries+1
        pair_prob[i,]=this_prob
    print("Mean value of sigma: %f" % np.mean(np.sqrt(1/(beta))))
    return pair_prob

dist=cal_pairwise_dist(x)
计算点对之间的距离平方的矩阵 dist

pair_prob=np.zeros((n,n))
概率矩阵pair_prob为（n,n）的全零矩阵，最终得到 $pair_prob ( i , j ) = P ( j ∣ i ) \text{pair\_prob}(i, j) = P(j|i)$
$P (j ∣ i)$ 表示在给定数据点 $i$ 的情况下选择数据点 $j$ 的条件概率。

beta=np.ones((n,1))
初始化 $\beta$ ：
$\beta = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}$

base_perp=np.log(perplexity)
确定目标困惑度

for i in range(n):
……
其余代码含义：使用二分法进行困惑度计算以确定pair_prob
对于每个数据点 $i$ ，调整参数 $\beta_i$ 使得该点的条件概率分布的困惑度接近目标困惑度 perplexity。数学上，目标是找到：
$Perplexity_target ∣ \beta_i = \arg\min_{\beta_i} | \text{Perplexity}(P(i|j; \beta_i)) - \text{Perplexity\_target} |$
其中，条件概率计算为
$\frac{\exp(-\beta_i \cdot \text{dist}(x_i, x_j))}{\sum_{k \neq i} \exp(-\beta_i \cdot \text{dist}(x_i, x_k))}$

具体步骤为：

计算当前困惑度差异：
$perp_diff = Perplexity ( P ( i ∣ j ; β i ) ) − Perplexity_target \text{perp\_diff} = \text{Perplexity}(P(i|j; \beta_i)) - \text{Perplexity\_target}$

调整 $\beta_i$ 使用二分法：

如果 $perp_diff > 0 \text{perp\_diff} > 0$ ：

则增加 $\beta_i$ ：

如果 $\beta_{\max}$ 未定义，则 $\beta_i = 2 \beta_i$

否则， $\beta_i = \frac{\beta_i + \beta_{\max}}{2}$

如果 $perp_diff < 0 \text{perp\_diff} < 0$ ：

则减少 $\beta_i$ ：

如果 $\beta_{\min}$ 未定义，则 $\beta_i = \frac{\beta_i}{2}$

否则， $\beta_i = \frac{\beta_i + \beta_{\min}}{2}$

更新上限或下限：

如果 $perp_diff > 0 \text{perp\_diff} > 0$ ，更新 $\beta_{\min} = \beta_i$

如果 $perp_diff < 0 \text{perp\_diff} < 0$ ，更新 $\beta_{\max} = \beta_i$
这段过程重复，直到 $perp_diff ∣ ≤ tol |\text{perp\_diff}| \leq \text{tol}$ 或达到最大尝试次数（50次）。
最终更新 $pair_prob [ i , : ] = this_prob \text{pair\_prob}[i,:] = \text{this\_prob}$ ，作为点 $i$ 的条件概率分布。

tsne(x,no_dims=2,perplexity=30.0,max_iter=1000)

def tsne(x,no_dims=2,perplexity=30.0,max_iter=1000):
    if isinstance(no_dims,float):
        print("Error: array x should have type float")
        return -1
    (n,d)=x.shape
    initial_momentum=0.5
    final_momentum=0.8
    eta=500
    min_gain=0.01
    y=np.random.randn(n,no_dims)
    dy=np.zeros((n,no_dims))
    iy=np.zeros((n,no_dims))
    gains=np.ones((n,no_dims))
    P=search_prob(x,1e-5,perplexity)
    P=P+np.transpose(P)
    P=P/np.sum(P)
    print("T_SNE DURING:%s" % time.process_time())
    P=P*4
    P=np.maximum(P,1e-12)
    for iter in range(max_iter):
        sum_y=np.sum(np.square(y),1)
        num=1/(1+np.add(np.add(-2*np.dot(y,y.T),sum_y).T,sum_y))
        num[range(n),range(n)]=0
        Q=num/np.sum(num)
        Q=np.maximum(Q,1e-12)
        PQ=P-Q
        for i in range(n):
            dy[i,:]=np.sum(np.tile(PQ[:,i]*num[:,i],(no_dims,1)).T*(y[i,:]-y),0)
        if iter<20:
            momentum=initial_momentum
        else:
            momentum=final_momentum
        gains=(gains+0.2)*((dy>0)!=(iy>0))+(gains*0.8)*((dy>0)==(iy>0))
        gains[gains<min_gain]=min_gain
        iy=momentum*iy-eta*(gains*dy)
        y=y+iy
        y=y-np.tile(np.mean(y,0),(n,1))
        if (iter+1)%100==0:
            C=np.sum(P*np.log(P/Q))
            print("Iteration",(iter+1),": error is",C)
            if (iter+1)!=100:
                ratio=C/oldC   
                print("ratio ",ratio)
                if ratio>=0.95:
                    break
            oldC=C
        if iter==100:
            P=P/4
    print("finished training ")
    return y

tsne(x,no_dims=2,perplexity=30.0,max_iter=1000)

$x$ : 原始数据矩阵，形状为 $ n \times d $（$ n $ 个样本，每个样本有 $ d $ 个特征）。

$no_dims \text{no\_dims}$ : 目标低维空间的维度，默认为 2。

$\text{perplexity}$ : 用于计算条件概率的困惑度。

$max_iter \text{max\_iter}$ : 最大迭代次数，默认为 1000。

initial_momentum=0.5
final_momentum=0.8
eta=500
min_gain=0.01
这些是 t-SNE 算法中的一些超参数，控制着梯度下降的行为：

initial_momentum = 0.5:
初始动量，用于在前几次迭代中加速收敛。动量是梯度下降中的一项技术，帮助跳出局部最小值。

final_momentum = 0.8:
最终动量。在算法运行一段时间后，增加动量的值以帮助更稳定地收敛。

eta = 500:
学习率，控制每次迭代时更新的步长大小。较大的学习率可能导致不稳定，较小的学习率可能导致收敛缓慢。

min_gain = 0.01:
增益的最小值，用于防止更新步长过小。增益用于调整学习率，以适应梯度变化。

这些参数需要根据具体数据集进行调整，以获得最佳结果。

y=np.random.randn(n,no_dims)
dy=np.zeros((n,no_dims))
iy=np.zeros((n,no_dims))
gains=np.ones((n,no_dims))
P=search_prob(x,1e-5,perplexity)
P=P+np.transpose(P)
P=P/np.sum(P)

$y$ : 随机初始化的低维表示，形状为 $no_dims n \times \text{no\_dims}$ 。

$\text{dy}$ , $\text{iy}$ : 梯度和动量，初始化为零矩阵。

$\text{gains}$ : 增益矩阵，初始化为全 1。

$P$ : 高维空间中点对的联合概率矩阵，经过对称化和归一化处理。

P=P*4；P=np.maximum(P,1e-12)
这段代码的作用是对概率矩阵 ( P ) 进行调整：
$\times 4$ :

将矩阵 $P$ 中的每个元素乘以 4。这个步骤在 t-SNE 中是为了在前几次迭代中提高 $P$ 的权重，使低维表示更快地形成。
$\text{np.maximum}(P, 1e-12)$ :

将 $P$ 中的每个元素与 $1 e - 12$ 进行比较，确保所有元素至少为 $1 e - 12$ 。这是为了避免对数计算中的零值问题（如计算 KL 散度时 $\log(0)$ 是未定义的）。这样可以提高数值稳定性。

for iter in range(max_iter)：

计算相似性矩阵：
$\text{num} = \frac{1}{1 + \|y_i - y_j\|^2}$
$\frac{\text{num}}{\sum \text{num}}$
$Q$ 是低维空间中的相似性矩阵。
关于P=np.maximum(P,1e-12)
np.maximum()函数的作用是逐元素比较两个数组，并返回两者中较大的值。表达式Q=np.maximum(Q,1e-12)的意思是将数组Q中的每个元素与1e-12进行比较，并将Q中每个小于1e-12的元素替换为1e-12。这样做通常是为了避免在后续的计算中出现数值上的不稳定。

计算梯度：
$\text{dy}[i,:] = \sum_j \left((P_{ij} - Q_{ij}) \cdot \text{num}_{ij} \cdot (y_i - y_j)\right)$
关于*np.tile(PQ[:,i]num[:,i],(no_dims,1))
将数组 PQ 的第 i 列与数组 num 的第 i 列对应的元素相乘，得到一个一维数组。然后将这个一维数组在垂直方向上重复 no_dims 次形成一个新的二维数组，而水平方向上保持不变。

更新动量和增益：
$\text{gains} = (\text{gains} + 0.2) \cdot ((\text{dy} > 0) \neq (\text{iy} > 0)) + (\text{gains} \cdot 0.8) \cdot ((\text{dy} > 0) == (\text{iy} > 0))$
$min_gain ] = min_gain \text{gains}[\text{gains} < \text{min\_gain}] = \text{min\_gain}$
$\text{iy} = \text{momentum} \cdot \text{iy} - \eta \cdot (\text{gains} \cdot \text{dy})$

更新低维表示：
$\text{iy}$
并使 $y$ 的均值为零。

计算和输出误差（KL散度）：
$\sum_{i,j} P_{ij} \log\left(\frac{P_{ij}}{Q_{ij}}\right)$

输出

返回：低维表示 $y$ 。
通过这些步骤，函数在高维数据中寻找低维表示，尽可能保留数据的局部结构。