p-SNE：为高维稀疏计数数据设计的非线性降维方法

1. 项目概述：当降维遇上泊松计数数据

在神经科学、基因组学、文本挖掘这些领域里，我们常常要和一种特殊的数据打交道：计数数据。比如，记录一个神经元在一秒内发放了多少次动作电位，测量一个细胞里某个基因的转录本有多少个拷贝，或者统计一篇文档里“人工智能”这个词出现了几次。这些数据天然就是非负整数，而且一个关键特性是，它们的波动（方差）往往和它们的平均水平（均值）紧密相关——均值高，波动通常也大。这种特性，用统计学的语言描述，就是近似服从泊松分布。更麻烦的是，在很多实验条件下，事件发生的概率本身就很低，导致数据矩阵里充斥着大量的零，我们称之为高维稀疏泊松数据。

面对成百上千个神经元、上万个基因或数万词汇表构成的高维数据，人脑直接理解几乎是不可能的。这时候，降维技术就成了我们窥探数据内在结构的“显微镜”。它能把数据从成百上千维压缩到两三维，画在一张图上，让我们直观地看到哪些样本聚在一起，数据背后有没有隐藏的“流形”结构。t-SNE和UMAP是过去十几年里最成功的非线性降维工具，几乎成了单细胞测序、复杂系统可视化的标配。

但问题来了。无论是经典的PCA，还是后来居上的t-SNE和UMAP，它们骨子里都假设数据是连续的，并且样本之间的距离用欧氏距离来衡量是合理的。这对于像图像像素值、物理传感器读数这类连续数据没问题。可对于计数数据，尤其是稀疏的计数数据，欧氏距离会严重失真。想象一下，两个样本在某个特征上，一个计数是0，另一个是1。从数值上看，欧氏距离是1。但在泊松分布的语境下，从0到1的“跳跃”，其统计意义远比从100到101的“跳跃”要重大得多。传统方法无视了这种分布特性，相当于用尺子去称重量，工具本身就不匹配。

所以，当我们需要分析神经脉冲序列、单细胞RNA测序计数矩阵或者文档-词频矩阵时，常常会感到一种“隔靴搔痒”的别扭。现有的方法要么像PCA一样只能捕捉线性结构，要么像t-SNE一样因为误用距离度量而扭曲了稀疏计数数据中真实的邻近关系。一个常见的补救办法是对数据做个log(1+x)变换，试图稳定方差，然后再喂给t-SNE。这办法虽然常用，但本质上是一种妥协，它丢弃了数据原始的泊松似然结构，效果时好时坏，尤其在数据非常稀疏时。

今天要深入探讨的p-SNE，就是为了解决这个根本矛盾而生的。它的全称是泊松随机邻域嵌入。这个方法的出发点非常直接：既然数据来自泊松过程，那么我们在衡量样本相似性、构建降维目标函数时，就应该从头到尾、彻彻底底地尊重泊松分布的性质。它不是对现有方法的小修小补，而是从第一性原理出发，为高维稀疏计数数据量身定制了一套非线性降维框架。

1.1 核心需求：为泊松数据设计专属的“距离”与“损失”

要理解p-SNE的创新，得先看明白传统方法在哪里“踩了坑”。t-SNE的核心流程可以简化为三步：1) 用欧氏距离计算高维空间样本间的差异；2) 基于这个差异，用高斯核函数将其转化为一个表示“邻近概率”的分布；3) 在低维空间也构造一个概率分布（通常用t分布核），然后最小化两个分布之间的KL散度，从而得到嵌入。

p-SNE的改造就发生在第一步和第三步。首先，它认为对于两个泊松分布，最自然的“差异”度量不是欧氏距离，而是KL散度。KL散度直接来源于泊松分布的概率公式，它内置了均值-方差耦合的特性，能正确评估低计数区间的微小差异。其次，在优化低维嵌入时，p-SNE没有沿用KL散度，而是选择了Hellinger距离。这是因为高维空间构建的概率分布可能非常稀疏（很多概率接近零），KL散度在这种情况下会趋于无穷大，导致优化不稳定。Hellinger距离则是对称、有界的，对零值更稳健。

所以，p-SNE解决的核心需求，就是为泊松计数数据提供一个原理正确、数值稳定的非线性降维工具。它适合任何需要分析计数数据、且怀疑数据中存在非线性结构的研究者和工程师。无论是神经科学家想可视化不同刺激下神经元集群的活动模式，生物信息学家想探索单细胞转录组的细胞亚群，还是数据科学家想对文档集合进行主题聚类，p-SNE都提供了一个比“log变换+t-SNE”更坚实、更可靠的选择。

2. p-SNE原理深度拆解：从泊松KL散度到Hellinger距离

要真正用好一个工具，不能只当“调包侠”，理解其背后的数学动机和工程考量至关重要。p-SNE的算法流程清晰，但每一步的设计都蕴含着对泊松数据特性的深刻考量。我们把它拆开揉碎了看。

2.1 第一步：用泊松KL散度重新定义“相似性”

给定一个计数矩阵 Y，形状是 N x M（N个样本，M个特征）。传统的t-SNE会直接计算样本i和样本j之间的欧氏距离。p-SNE则完全不同，它从每个特征开始计算。

对于样本 i 和样本 j 的第 m 个特征，我们有观测计数 y_i,m 和 y_j,m。p-SNE将它们视为两个泊松分布 Poisson(y_i,m) 和 Poisson(y_j,m) 的均值参数（这里用一个简化观点，将观测值作为分布参数的估计）。那么，衡量这两个分布差异的自然工具就是KL散度。对于泊松分布，KL散度有一个漂亮的解析形式：

KL(Pois(λ1) || Pois(λ2)) = λ1 * log(λ1/λ2) + λ2 - λ1

在p-SNE的实现中，公式稍有变化，加入了平滑项 ε 防止除零错误： d_{m,(i,j)} = y_{i,m} * log((y_{i,m}+ε)/(y_{j,m}+ε)) + y_{j,m} - y_{i,m}

为什么这个公式更合理？ 我们来看个例子。假设特征m上，样本A计数为1，样本B计数为2，样本C计数为100，样本D计数为101。

• 欧氏距离：dist(A,B) = 1, dist(C,D) = 1。它认为这两对样本的差异一样大。
• 泊松KL散度（近似计算）：KL(A||B) ≈ 1*log(1/2)+2-1 = 0.31， KL(C||D) ≈ 100*log(100/101)+101-100 ≈ 0.005。

看，泊松KL散度明确告诉我们，在低计数区域（1 vs 2），差异是显著的（0.31）；而在高计数区域（100 vs 101），同样的绝对差值在统计上几乎可以忽略不计（0.005）。这正是处理计数数据时我们直觉上期望的：重视低计数区域的相对变化，容忍高计数区域的绝对波动。欧氏距离完全丢失了这种信息。

然后，p-SNE将所有M个特征上的KL散度相加，得到样本i和j之间的总差异 D_{i,j}。注意，由于KL散度本身不对称，D_{i,j} 也不等于 D_{j,i}。这引出了下一个问题：如何构建一个对称的“邻近概率”？

2.2 第二步：构建对称的高维空间概率分布

上一步我们得到了一个非对称的差异矩阵 D。D_{i,j} 小，表示从样本i的视角看，样本j很“像”它。但反过来不一定成立。为了得到一个全局一致的邻近关系，我们需要对称化。

p-SNE的做法很直观：

1. 对于每个样本i，以其到其他所有样本的差异 D_{i,:} 为基础，计算一个条件概率分布。具体是用softmax函数：p_{j|i} = exp(-w * D_{i,j}) / Σ_{k≠i} exp(-w * D_{i,k})。参数 w 可以控制分布的“集中度”，w 越大，概率质量越集中在差异最小的几个邻居上。
2. 由于 p_{j|i} ≠ p_{i|j}，我们通过简单平均来对称化：s_{i,j} = (p_{j|i} + p_{i|j}) / (2N)。最后除以 2N 是为了让所有 s_{i,j} 之和为1，使其成为一个合法的联合概率分布，记作矩阵 S。

这个 S 矩阵就是p-SNE所认为的、在高维空间中数据点之间“真实”的邻近关系。它完全基于泊松统计特性构建，是后续所有优化的基石。

注意：参数 w 的选择类似t-SNE中的“困惑度”，它决定了嵌入结果中局部结构的精细程度。w 值越大，嵌入越关注非常局部的结构，可能产生更多、更分散的小簇；w 值越小，则更关注全局结构，簇会更紧凑但可能混合。在真实数据中，通常需要尝试几个值（如0.5, 1.0, 2.0）。

2.3 第三步：用Hellinger距离优化低维嵌入

现在，我们有了高维空间的概率分布 S，目标是找到一个低维嵌入 X（比如每个样本用2维或3维坐标表示），使得在低维空间中定义的样本间概率分布 Q 尽可能接近 S。

p-SNE在低维空间构造 Q 的方式和t-SNE一致，使用自由度为1的t分布核： q_{i,j} = (1 + ||x_i - x_j||^2)^{-1} / Σ_{k≠l} (1 + ||x_k - x_l||^2)^{-1} 选择t分布核主要是为了解决“拥挤问题”：在高维空间中，许多样本点有足够的空间彼此远离，但在映射到低维（如2D）时，空间不够，导致所有点挤在一起。t分布的重尾特性允许在低维空间中，中等距离的点可以离得更远，从而缓解拥挤。

关键来了：如何衡量 S 和 Q 的差异？t-SNE使用的是KL散度。但p-SNE在这里换成了Hellinger距离。对于两个离散概率分布P和Q，Hellinger距离定义为： H(P, Q) = sqrt( 1/2 * Σ_i (sqrt(p_i) - sqrt(q_i))^2 )

为什么在这里要换用Hellinger距离？ 这完全是出于数值稳定性和优化鲁棒性的工程考量。回忆一下，我们的高维分布 S 是在稀疏计数数据上构建的，很多样本对之间的相似度 s_{i,j} 可能非常接近甚至等于0。KL散度 KL(S||Q) = Σ s * log(s/q) 中，只要有一个 s_{i,j}=0 而对应的 q_{i,j}>0，这一项就是0，没问题；但反过来，如果某个 s_{i,j}>0 而 q_{i,j}=0，log(s/0) 就会趋向无穷大，导致损失函数爆炸。在优化初期，低维嵌入 X 是随机初始化的，Q 的分布很可能在某些位置给出接近零的概率，这就极易引发数值问题。

Hellinger距离完美避开了这个坑。首先，它是对称的：H(S, Q) = H(Q, S)。其次，它是有界的，取值范围在[0,1]之间。最重要的是，它的计算涉及的是概率的平方根 sqrt(s) 和 sqrt(q)。即使 s 或 q 为零，平方根依然是零，整个项 (sqrt(s)-sqrt(q))^2 依然是一个良定义的有限值。这使得基于梯度的优化过程（如带动量的梯度下降）更加平滑和稳定。

所以，p-SNE的最终损失函数是：L = H(S, Q)。通过最小化这个损失，我们驱动低维嵌入 X 不断调整，直到由它产生的样本间相似性分布 Q，与基于泊松KL散度构建的高维相似性分布 S 尽可能一致。

2.4 算法实现与优化细节

p-SNE的完整算法流程如下，我结合自己的实现经验，补充一些论文中未详述但至关重要的细节：

1. 输入：计数矩阵 Y, 嵌入维度 P, 锐度参数 w, 学习率 η, 动量 μ, 早期放大因子 α 及其迭代次数 K_ex, 最大迭代次数 K, 容忍度 τ, 稳定常数 ε。
2. 计算差异矩阵D：遍历所有样本对 (i, j)，按公式计算 D_{i,j}。这是一个 O(N^2 * M) 的计算，是算法的主要瓶颈。对于大规模数据，需要优化或采样。
3. 计算高维分布S：根据 D 计算条件概率并对称化，得到 S。
4. 初始化与早期放大：低维嵌入 X 通常从自由度为3的t分布中随机采样初始化。在最初的 K_ex 次迭代中，采用“早期放大”技巧：将 S 矩阵的所有元素乘以一个大于1的因子 α，然后重新归一化得到 S_eff。这个技巧能放大样本间的吸引力，帮助在优化初期快速形成宏观的簇结构，避免陷入局部最优。
5. 迭代优化：
   • 根据当前 X 计算低维分布 Q。
   • 计算损失 L = H(S_eff, Q)。在 K_ex 次迭代后，S_eff 恢复为原始的 S。
   • 计算损失 L 对 X 的梯度 ∂L/∂X。这个梯度有解析形式，涉及 (sqrt(s)-sqrt(q)) / (sqrt(q) * (1+||x_i-x_j||^2)) 等项。
   • 使用带动量的梯度下降更新 X：v_{t+1} = μ * v_t - η * ∂L/∂X_t； X_{t+1} = X_t + v_{t+1}。
   • 检查损失变化，如果小于容忍度 τ 则提前终止。

实操心得：参数选择经验

• 学习率 η：通常从10到100开始尝试。太大容易震荡，太小收敛慢。可以配合学习率衰减策略。
• 动量 μ：一般设为0.8到0.9，能加速收敛。
• 早期放大因子 α：类似t-SNE，常用12或15。放大迭代次数 K_ex 一般为总迭代次数的1/4到1/3。
• 稳定常数 ε：一个很小的数，如1e-12，防止计算KL散度时对零取对数。
• 锐度参数 w：这是最重要的超参数之一。对于簇结构明显的数据，w=1.0 通常是个不错的起点。对于流形结构连续、稀疏度高的数据（如神经序列），可以尝试更大的 w（如2.0）来捕捉更精细的局部结构。建议用网格搜索结合下游任务（如聚类指标）来确定。

3. 实战对比：p-SNE在合成与真实数据上的表现

理论再漂亮，也得看实战效果。原论文通过合成数据和三个真实世界数据集，系统地对比了p-SNE与一系列主流方法。我们不仅复现这些对比，更深入解读其背后的原因，并分享我在复现过程中遇到的一些坑和技巧。

3.1 合成数据实验：在已知的战场上检验

合成数据的优势在于我们有“标准答案”，可以定量评估算法恢复真实结构的能力。

实验一：角嵌入数据 这个数据集模拟了神经群体记录：60个样本（试次），40个特征（神经元），分为3组。每组样本的泊松发放率在一个圆环流形上按角度排列。数据有一定稀疏性。

我按照论文描述生成了数据，并对比了p-SNE (w=1.0)、原始t-SNE、对数变换后t-SNE (t-SNE-log)、平方根变换后t-SNE (t-SNE-sqrt)、UMAP、Isomap、MDS、谱嵌入以及Ling & Xue (2022)的方法。

视觉对比：将结果降到3维后可视化（颜色代表真实组别），差异一目了然。

• p-SNE：三个组被清晰地分离成三个簇，并且在空间中保持了它们原始的相对角度关系（如果将颜色换成连续的角度值，可以看到一个平滑的色环）。这说明p-SNE不仅抓住了离散的簇，还保留了流形的连续几何。
• t-SNE及其变体：所有点几乎坍缩成一个致密的球，组间完全混合。对数/平方根变换并没有从根本上解决问题。
• UMAP和谱嵌入：它们试图发现结构，但似乎“用力过猛”，将一个本应连续的组割裂成了几个不连续的子簇，引入了虚假的边界。
• Isomap, MDS, NMF：表现出部分分离，但组与组之间的边界模糊，有大量重叠点。

定量评估：我计算了k近邻分类、SVM分类以及k-means聚类调整兰德指数，并对所有方法在这些指标上的表现进行排名（排名越靠前越好）。p-SNE的综合排名是第一。这证实了在泊松数据上，使用正确的差异度量至关重要。

稀疏度鲁棒性测试：我逐步降低合成数据的泊松率参数，人为增加数据中的零值比例（从56%到84%）。随着稀疏度增加，所有方法的性能都在下降，但p-SNE的下降曲线最为平缓。在中等至高稀疏度（56%-74%零值）下，p-SNE相对于其他方法的优势最大。当稀疏度达到惊人的84%时，所有方法都近乎失效，但p-SNE的结果中仍能依稀分辨出簇的轮廓，而其他方法已完全混乱。

避坑指南：合成数据生成 在复现角嵌入数据时，关键是要正确模拟“神经元的调谐曲线”。即每个神经元 m 对一个特定的角度 θ_m 有最大发放率。样本 n 的刺激角度为 φ_n，那么该神经元在该试次下的期望发放率 λ_{n,m} = A * exp(-|φ_n - θ_m|^2 / (2σ^2))，其中A是最大率，σ是调谐宽度。然后从 Poisson(λ_{n,m}) 中采样计数。务必确保 λ 矩阵没有负值或异常值，否则KL散度计算会出问题。建议生成后检查每个样本的计数分布直方图，确认其近似泊松形状。

实验二：稀疏序列嵌入数据 这个数据集更具挑战性：119个样本，30个特征，泊松率极低（0.1到2.6），导致数据稀疏度超过90%。样本沿着一个一维连续流形排列。

在这个任务上，p-SNE (w=2.0) 的表现令人印象深刻。其3维嵌入呈现出一个平滑的、与真实流形参数高度相关的颜色梯度。这意味着p-SNE成功地从极度稀疏和嘈杂的计数数据中，重建出了潜在的连续轨迹。

相比之下，t-SNE系列方法再次坍缩。UMAP将数据分成了两大块，但丢失了内部的顺序信息。Isomap未能捕捉到序列 progression。谱嵌入在Spearman相关系数上略胜p-SNE，排名第一，但其可视化结果显示出不连续的跳跃，而p-SNE的梯度更平滑，这可能对下游分析（如伪时间分析）更友好。

这个实验强有力地说明，对于像神经序列、单细胞发育轨迹这类具有连续流形结构的稀疏计数数据，p-SNE相比传统方法有显著优势。

3.2 真实世界数据挑战

纸上得来终觉浅，我们看看p-SNE在真实、杂乱的数据上表现如何。

数据集一：个人邮件通信记录 数据来自作者邮箱，统计了370天内与150位最频繁联系人的每日邮件数量。这是一个典型的稀疏计数矩阵，我们想看看降维能否区分工作日和周末的邮件模式。

我处理数据时，首先构建了 370天 x 150发件人 的计数矩阵。然后应用p-SNE (w=2.0) 和所有基线方法降至3维。

结果：p-SNE的嵌入图中，工作日（蓝色）和周末（橙色）的样本点形成了两个虽有重叠但中心分离的云团。其他方法，如t-SNE、MDS等，则完全混合在一起。我训练了一个SVM分类器（RBF核，5折交叉验证），直接使用3维嵌入坐标来预测“工作日/周末”。p-SNE达到了最高的分类准确率。这证明，p-SNE嵌入所保留的区分信息量最大。

思考：为什么p-SNE能成功？我检查了数据，发现有些发件人（如工作同事）在工作日频繁出现，周末几乎为零；而有些（如家人、朋友）在周末更活跃。这种模式反映在计数上，就是某些特征（发件人）的计数分布在两类样本上有显著差异。泊松KL散度对这种分布差异非常敏感，而欧氏距离则容易被高计数发件人的微小日常波动所干扰。

数据集二：OpenReview学术论文摘要 我从ICLR等会议的OpenReview API获取了摘要，构建了 350篇论文 x 100个高频词 的词频矩阵，并为每篇论文人工标注了5个研究领域（如强化学习、图神经网络等）。稀疏度约82%。

p-SNE的表现：在3维嵌入中，五个研究领域形成了五个清晰可辨、分离度很好的簇。这意味着，仅基于词频，p-SNE就能将不同领域的论文区分开来。我同样用嵌入坐标进行领域分类，p-SNE在逻辑回归、kNN、SVM和K-means ARI指标上全面领先。

一个有趣的细节：在预处理时，我像论文中一样，移除了领域相关的关键词和常见停用词。这迫使模型必须从更一般的词汇使用模式中学习区分性特征。p-SNE的成功，说明它捕捉到了不同领域论文在“语言风格”或“主题分布”上的深层统计差异，而不仅仅是几个关键词的有无。

数据集三：艾伦脑观测站神经记录 这是最具挑战性的神经科学真实数据。我选择了视觉编码数据集中的一个会话，包含1193个质量良好的神经元单元，记录了它们在呈现各种视觉刺激时的脉冲计数。

分析目标：我们想知道，降维能否揭示大脑对不同视觉刺激（如不同朝向的光栅）的表征结构。理想情况下，对相似刺激反应的神经元应该在嵌入空间中靠近。

处理流程：

1. 数据提取：从Neuropixels记录中提取脉冲计数，构建 刺激条件 x 神经元 的计数矩阵。这是一个非常稀疏的矩阵。
2. 应用p-SNE：我将数据嵌入到2维和3维。在3维嵌入中，当我根据刺激的朝向给点着色时，观察到了一个初步的、具有一定顺序的色环结构。这表明对相似朝向有反应的神经元集群在嵌入空间中靠得更近。
3. 与基线对比：t-SNE的结果看起来更像是一团无结构的“星云”。UMAP显示出一些簇，但这些簇与刺激朝向的对应关系不如p-SNE清晰。PCA（作为线性方法）的前两个主成分也能展示一些分离，但非线性结构（如环状）完全丢失。

结论：在真实的、高噪声、高稀疏度的神经数据上，p-SNE能够比现有方法更清晰地揭示与刺激特征相关的神经群体活动结构。这对于神经科学家理解大脑的信息编码方式具有重要意义。

4. p-SNE的优势、局限与实战部署指南

经过原理剖析和实验对比，我们可以更系统地总结p-SNE的优劣，并给出具体的应用建议。

4.1 核心优势与适用场景

1. 原理正确性：这是p-SNE最大的优势。它从生成模型（泊松过程）出发，使用KL散度作为差异度量，在数学上为计数数据提供了最自然的相似性定义。这不再是启发式的修补，而是基于统计原理的推导。
2. 对稀疏数据的鲁棒性：Hellinger距离作为损失函数，对概率分布中的零值不敏感，使得优化过程在面对高度稀疏数据时依然稳定。实验也证明，在高稀疏度下，p-SNE的性能衰减远慢于其他方法。
3. 能捕捉非线性流形：继承了t-SNE系列的非线性嵌入能力，能够发现数据中复杂的簇结构和连续轨迹。
4. 无需复杂的预处理：不同于一些深度学习方法（如scVI, P-VAE），p-SNE不需要对数据进行复杂的归一化或批次校正。它直接处理原始计数，减少了因预处理引入偏差的风险。

最适合p-SNE的场景包括：

• 神经科学：神经脉冲计数矩阵的分析、解码和可视化。
• 计算生物学：单细胞RNA测序（scRNA-seq）数据，特别是基于UMI的计数数据。
• 文本分析：文档-词频矩阵的主题建模和可视化。
• 社交网络与推荐系统：用户-物品交互计数（如点击次数、观看时长离散化）。
• 任何产生低速率、离散事件计数的领域。

4.2 当前局限与挑战

没有完美的算法，p-SNE也有其局限：

1. 计算复杂度：和t-SNE一样，p-SNE需要计算所有样本对之间的差异矩阵 D，其复杂度为 O(N^2 M)。对于样本数 N 上万的大规模数据集，这会成为内存和计算时间的瓶颈。虽然可以采用Barnes-Hut树等近似算法加速，但如何将其适配到泊松KL散度的计算中，仍需进一步研究。
2. 超参数调优：锐度参数 w 对结果影响显著，且没有像t-SNE中“困惑度”那样直观的解释。需要依靠下游任务指标进行调优，增加了使用成本。
3. 仅限于泊松假设：p-SNE严格假设数据服从泊松分布。对于过度离散（方差远大于均值）或欠离散（方差远小于均值）的计数数据，其效果可能会打折扣。例如，在一些scRNA-seq数据中，由于技术噪音，数据可能更符合负二项分布。
4. 解释性：和大多数非线性降维方法一样，p-SNE得到的低维嵌入是黑箱的。我们很难说清楚某个维度具体代表了原始高维空间中的什么“因素”或“概念”。

4.3 实战部署指南与代码片段

假设你是一个Python用户，手头有一个 scRNA-seq 计数矩阵 adata.X（Anndata对象），想用p-SNE进行可视化。以下是一个概念性的实现步骤和关键代码片段（请注意，需要你自己实现或寻找p-SNE的第三方库，目前可能没有标准实现）。

重要提示：上述代码仅为算法逻辑的示意，不可直接用于生产。完整的p-SNE实现需要：

1. 高效的KL散度计算：使用向量化操作和可能的近似方法避免 O(N^2M) 的双重循环。
2. 完整的梯度计算：推导并实现Hellinger距离对低维坐标 X 的梯度 ∂L/∂X。
3. 优化技巧：实现带动量的梯度下降、早期放大、学习率调度等。
4. 对称化加速：借鉴t-SNE的Barnes-Hut算法或其它近似最近邻方法，只计算最近邻的 S 矩阵元素，将复杂度从 O(N^2) 降至 O(N log N)。

4.4 常见问题与排查技巧实录

在实际应用p-SNE或类似方法时，你可能会遇到以下问题：

问题1：嵌入结果全部挤在一起，没有分离。

• 可能原因1：锐度参数 w 太小。w 太小使得高维相似度分布 S 过于平坦，所有点彼此相似，导致低维嵌入无法形成结构。解决方案：增大 w（如从1.0调到2.0, 5.0），让模型更关注局部邻居。
• 可能原因2：学习率 η 太大或太小。太大导致优化震荡，无法收敛到好解；太小导致收敛过慢，在给定迭代次数内看不到结构。解决方案：尝试典型的学习率，如10, 50, 100, 200，并观察损失曲线是否平稳下降。
• 可能原因3：数据本身确实没有明显的簇结构。这是最根本的原因。解决方案：用其他方法（如PCA、层次聚类）交叉验证，或检查数据的统计特性。

问题2：嵌入结果不稳定，每次运行都不一样。

• 可能原因：随机初始化。非线性降维对初始化敏感。解决方案：
  1. 使用固定的随机种子 (np.random.seed) 确保可复现性。
  2. 尝试用PCA的前几个主成分作为初始化，这通常比纯随机初始化更稳定，并能提供一个更好的优化起点。
  3. 多次运行，选择损失函数值最小的一次结果。

问题3：计算太慢，内存占用太高。

• 可能原因：全连接相似度矩阵。计算所有样本对的KL散度，复杂度为 O(N^2M)。解决方案：
  1. 降采样：如果数据量太大，先对数据进行随机降采样或使用核心集（coreset）方法。
  2. 特征选择：减少特征数 M。对于基因表达数据，只保留高变异基因；对于文本数据，使用TF-IDF筛选。
  3. 近似算法：等待或贡献p-SNE的近似算法实现（如基于k近邻的稀疏相似度计算）。
  4. 增量/在线学习：对于流式数据，研究在线版本的p-SNE。

问题4：如何选择嵌入维度 P？

• 建议：对于可视化，当然是2或3维。但如果是为了下游任务（如聚类、分类）获取特征，可以尝试更高的维度（如5, 10, 20）。一个实用的方法是：运行多次p-SNE，每次用不同的随机种子，然后检查在多个维度下，下游任务（如聚类稳定性）的表现。选择一个性能饱和且稳定的维度。原论文的补充材料也显示，p-SNE在2到20维之间性能相对稳定。

问题5：p-SNE和深度生成模型（如scVI）比，如何选择？

• 对比：
  • p-SNE：原理清晰，相对轻量，适合探索性数据分析和可视化。计算一次得到一个嵌入。
  • scVI/P-VAE：是强大的生成模型，能建模批次效应、缺失数据，并估计不确定性。训练好后可以快速为新数据生成嵌入。但需要大量数据训练，计算成本高，模型更复杂。
• 选择指南：
  • 如果你的主要目标是快速可视化和探索数据中的非线性簇/流形，并且数据量不是特别巨大（数万样本以内），p-SNE是一个优秀且原理更透明的选择。
  • 如果你的数据有严重的批次效应，需要进行数据整合，或者你需要一个可以生成新数据、进行假设检验的模型，那么scVI这类深度生成模型更合适。
  • 在实际项目中，我经常两者都用：先用p-SNE快速看下数据的整体结构和是否存在明显的批次效应，然后再决定是否投入资源训练更复杂的scVI模型。

p-SNE的出现，为高维稀疏计数数据的分析提供了一把更精准的尺子。它提醒我们，在处理数据时，尊重其底层的数据生成过程和统计特性，往往能带来更深刻、更可靠的洞察。虽然它目前还存在计算扩展性的挑战，但其核心思想——为特定数据类型定制差异度量和损失函数——无疑为未来的降维方法发展指明了一个有价值的方向。在神经科学、生物信息学等领域，面对日益增长的复杂计数数据，像p-SNE这样“量体裁衣”的工具，其价值只会越来越凸显。