破解密度代码-为什么-MAF-流在-KDE-停滞不前

破解密度代码：为什么 MAF 流在 KDE 停滞不前

原文链接

第一部分：简介

在高维密度估计中出现的一个主要问题是，随着我们的维度增加，我们的数据变得更加稀疏。因此，对于依赖于局部邻域估计的模型，随着维度的增加，我们需要指数级更多的数据来继续获得有意义的成果。这被称为维度诅咒。

在我关于密度估计的前一篇文章中，我展示了核密度估计器（KDE）如何有效地用于一维数据。然而，在更高维度的性能会显著下降。为了说明这一点，我进行了一个模拟，以确定在估计多变量高斯分布的密度时，KDE 需要多少样本才能达到 0.2 的平均相对误差。带宽是使用 Scott 规则选择的。结果如下：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
np.random.seed(42)

# Gaussian sample generator
def generate_gaussian_samples(n_samples, dim, mean=0, std=1):
    return np.random.normal(mean, std, size=(n_samples, dim))

def compute_bandwidth(samples):
    # Scott method
    n, d = samples.shape
    return np.power(n, -1./(d + 4))

# KDE error computation
def compute_kde_error(samples, dim, n_test=1000):
    bandwidth = compute_bandwidth(samples)
    kde = KernelDensity(bandwidth=bandwidth).fit(samples)
    test_points = np.random.normal(0, 1, size=(n_test, dim))
    kde_density = np.exp(kde.score_samples(test_points))
    true_density = np.exp(-np.sum(test_points**2, axis=1) / 2) / ((2 * np.pi)**(dim / 2))
    error = np.mean(np.abs(kde_density - true_density) / true_density)
    return error, bandwidth

# Determine required samples for a target error
def find_required_samples(dim, target_error=0.2, max_samples=500000, start_samples=10, n_experiments=5):
    samples = start_samples
    while samples <= max_samples:
        errors = [compute_kde_error(generate_gaussian_samples(samples, dim), dim)[0] for _ in range(n_experiments)]
        avg_error = np.mean(errors)
        if avg_error <= target_error:
            return samples, avg_error
        samples = int(samples * 1.5)
    return max_samples, avg_error

# Main
def analyze_kde(dims, target_error):
    results = []
    for dim in dims:
        samples, error = find_required_samples(dim, target_error)
        results.append((dim, samples))
        print(f"Dim {dim}: {samples} samples")
    return results

# Visualization
def plot_results(dims, results,target_error=.2):
    samples = [x[1] for x in results]
    plt.figure(figsize=(8, 6))
    plt.plot(dims, samples, 'o-', color='blue')
    plt.yscale('log')
    plt.xlabel('Dimension')
    plt.ylabel('Required Number of Samples (log scale)')
    plt.title(f'Samples Needed for a Mean Relative Error of {target_error}')
    plt.grid(True)

    for i, sample in enumerate(samples):
        plt.text(dims[i], sample * 1.15, f'{sample}', fontsize=10, ha='right', color='black')  
    plt.show()

# Run the analysis
dims = range(1, 7)
target_error = 0.2
results = analyze_kde(dims, target_error)
plot_results(dims, results) 

对了：在我的模拟中，为了匹配一维中仅 22 个数据点的准确性，你需要在六个维度中需要超过 360,000 个数据点！更令人惊讶的是，在 David W. Scott 的《多元密度估计》一书中，他表明，根据度量标准，在八个维度中需要超过一百万个数据点才能达到一维中仅 50 个数据点的相同准确性。

希望这足以说服你，核密度估计器在估计高维密度时并不理想。但替代方案是什么呢？

---

第二部分：归一化流简介

一个有希望的替代方案是归一化流，而我将重点关注的特定模型是掩码自回归流（MAF）。

本节部分借鉴了 George Papamakarios 和 Balaji Lakshminarayanan 的工作，这些工作在 Kevin P. Murphy 所著的《概率机器学习：高级主题》（Probabilistic Machine Learning: Advanced Topics）的第二十三章中有所介绍（详见该书以获取更多细节）。

正态化流背后的核心思想是，可以通过从简单的基分布（例如高斯分布）中抽取随机变量开始，然后通过一系列可微、可逆的变换（微分同胚）来模拟分布 p(x)。每个变换逐步重塑分布，逐渐将基分布映射到目标分布。下面是这个过程的一个可视化说明。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
np.random.seed(42)

#Sample from a standard normal distribution
n_points = 1000
initial_dist = np.random.normal(loc=[0, 0], scale=1.0, size=(n_points, 2))

#Generate target distribution
theta = np.linspace(0, np.pi, n_points//2)
r = 2
x1 = r * np.cos(theta)
y1 = r * np.sin(theta)
x2 = (r-0.5) * np.cos(theta)
y2 = (r-0.5) * np.sin(theta) - 1
target_dist = np.vstack([
    np.column_stack([x1, y1 + 0.5]),
    np.column_stack([x2, y2 + 0.5])
])
target_dist += np.random.normal(0, 0.1, target_dist.shape)

def f1(x, t):
    """Split transformation"""
    shift = 2 * t * np.sign(x[:, 1])[:, np.newaxis] * np.array([1, 0])
    return x + shift

def f2(x, t):
    """Curve transformation"""
    theta = t * np.pi / 2
    r = np.sqrt(x[:, 0]**2 + x[:, 1]**2)
    phi = np.arctan2(x[:, 1], x[:, 0]) + theta * (1 - r/4)
    return np.column_stack([r * np.cos(phi), r * np.sin(phi)])

def f3(x, t):
    """Fine-tune to target"""
    return (1 - t) * x + t * target_dist

# Create figure
fig, ax = plt.subplots(figsize=(10, 10))
scatter = ax.scatter([], [], alpha=0.6, s=10)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)

def sigmoid(x):
    """Smooth transition function"""
    return 1 / (1 + np.exp(-(x - 0.5) * 10))

def get_title(t):
    if t < 0.33:
        return f'Applying Split Transformation (f₁)'
    elif t < 0.66:
        return f'Applying Curve Transformation (f₂)'
    else:
        return f'Fine-tuning to Target Distribution (f₃)'

def init():
    scatter.set_offsets(initial_dist)
    ax.set_title('Initial Gaussian Distribution', pad=20, fontsize=18)
    return [scatter]

def update(frame):
    #Normalize frame to [0, 1]
    t = frame / 100

    #Apply transformations sequentially
    points = initial_dist

    #f1: Split the distribution
    t1 = sigmoid(t * 3) if t < 0.33 else 1
    points = f1(points, t1)

    #f2: Create curves
    t2 = sigmoid((t - 0.33) * 3) if 0.33 <= t < 0.66 else (0 if t < 0.33 else 1)
    points = f2(points, t2)

    #f3: Fine-tune to target
    t3 = sigmoid((t - 0.66) * 3) if t >= 0.66 else 0
    points = f3(points, t3)

    #Update scatter plot
    scatter.set_offsets(points)
    colors = points[:, 0] + points[:, 1]
    scatter.set_array(colors)

    #Update title
    ax.set_title(get_title(t), pad=20, fontsize=18)

    return [scatter]

#Create animation
anim = FuncAnimation(fig, update, frames=100, init_func=init,
                    interval=50, blit=True)
plt.tight_layout()
plt.show()

#Save animation as a gif
anim.save('normalizing_flow_single.gif', writer='pillow')

更正式地说，假设以下内容：

然后，我们的目标分布由以下变量变换公式定义：

其中 J_{f^{-1}}(x) 是在 x 处评估的 f^{-1} 的雅可比矩阵。

由于我们需要计算行列式，也存在计算上的考虑；我们的变换函数的理想情况是具有行列式易于计算的雅可比矩阵。设计一个既能模拟复杂分布又能产生可处理行列式的同胚函数是一个具有挑战性的任务。在实践中，这是通过通过更简单函数的流来构建目标分布来解决的。因此，f 定义如下：

然后，由于同胚函数的复合也是同胚的，f 将是可逆且可微的。

对于 f 有几个典型的候选者。以下是一些流行的选择。

仿射流

仿射流由以下函数给出：

我们需要将 A 限制为可逆的方阵，以便 f 是可逆的。仿射流本身在建模数据方面并不很好，但与其他函数混合使用时很有用。

元素流

元素流逐元素地变换向量 u。设 h 是一个标量双射，我们可以创建一个如下定义的向量值双射 f：

雅可比矩阵的行列式如下给出：

与仿射流类似，元素流本身在建模数据方面并不非常有效，因为它们没有捕捉到维度之间的交互。然而，它们经常与其他变换结合使用来构建更复杂的流。

配对流

配对流，由 Dinh 等人于 2015 年提出，与之前讨论的流不同，因为它们允许使用非线性函数更好地捕捉数据的结构。在这里使用图片表示，但为了避免混淆，我需要使用内联 LaTeX。

在这里，f-hat 的参数是通过将 u 的子集 b 通过 Θ 计算得到的，其中 Θ 是一个称为条件器的通用函数。这种设置与仿射流形成对比，仿射流只线性混合维度，与元素流形成对比，元素流保持每个维度独立。配对流通过条件器允许维度的非线性混合。如果您对已提出的耦合层类型感兴趣，请参阅 Kobyzev, Ivan & Prince, Simon & Brubaker, Marcus. (2020)。

由于 x-b 对 u-b 的偏导数为 0，行列式的计算相当简单。因此，雅可比矩阵是以下上三角矩阵：

雅可比矩阵的行列式如下给出：

以下展示了这些函数如何影响分布的视觉效果。

掩码自回归流

假设 u 是一个包含 d 个元素的向量。一个自回归双射函数，输出一个包含 d 个元素的向量 x，定义为以下：

在这里，h 是一个由 Θ 参数化的标量双射参数，其中 Θ 是一个任意的非线性函数，通常是神经网络。由于自回归结构，每个元素 x_i 仅依赖于 u 的前 i 个索引的元素。因此，雅可比矩阵将是三角矩阵，其行列式将是主对角线元素的乘积，如下所示：

如果我们使用多个自回归双射函数作为我们的 f，我们需要训练 d 个不同的神经网络，这可能会非常计算量大。因此，为了解决这个问题，实践中更有效的方法是将条件器组合成一个单一模型 Θ，该模型接受共享输入 x 并输出参数集（Θ1、Θ2、…、Θd）。然而，为了保持自回归结构，我们必须确保每个 Θi 仅依赖于 x1、x2、…、xi−1。

掩码自回归流（MAF）使用多层感知器作为非线性函数，然后应用掩码来消除任何违反自回归结构的计算路径。通过这样做，MAF 确保每个输出 Θi 仅条件依赖于之前的输入 x1、x2、…、xi−1，从而允许高效训练。

---

第三部分：对决

为了确定 KDE 或 MAF 在更高维度的分布建模中哪个更好，我设计了一个实验，这个实验与我对 KDE 的初步分析类似。我逐渐在更大的数据集上训练这两个模型，直到每个模型达到 0.5 的 KL 散度。

对于那些不熟悉这个度量的人来说，KL 散度量化了一个概率分布与参考分布之间的差异。具体来说，它衡量使用一个分布来近似另一个分布时的预期“惊喜”的过剩。KL 散度为 0.0 表示分布之间完美对齐，而更高的值表示更大的差异。为了提供直观的视觉理解，下面的图示展示了比较两个三维分布时 0.5 KL 散度的样子：

0.5 KL 散度可视化

实验设计涵盖了三个不同的分布家族，每个家族都被选择来测试模型能力的不同方面。首先，我检查了条件高斯分布，这是最简单的情况，具有单峰、对称的概率质量。其次，我测试了条件高斯混合分布，引入了多峰性来挑战模型捕捉数据中多个不同模式的能力。最后，我包括了条件偏斜正态分布，以评估在非对称分布上的性能。

对于核密度估计模型，选择合适的带宽参数对于较大的维度来说是一项挑战。我最终采用了留一法交叉验证（LOOCV），这是一种技术，其中每个数据点被排除在外，而剩余的点被用来估计最佳带宽。这个过程虽然计算成本高，需要为 n 个数据点进行 n 个单独的模型拟合，但在高维度上获得可靠结果却是必要的。在我之前使用替代带宽选择方法的实验版本中，所有方法都表现出较差的性能，需要大量更多的训练数据才能达到相同的 KL 散度阈值。

Masked Autoregressive Flow 模型需要不同的优化策略。像大多数基于神经网络的模型一样，MAF 依赖于多个超参数。我开发了一种缩放策略，其中这些超参数与输入的维度成比例调整。重要的是要注意，这种缩放是基于合理的启发式方法，而不是彻底的优化。超参数搜索被保持在最小，以建立基线性能，更复杂的调整可能会为 MAF 模型带来大的性能提升。

包括数据生成、模型实现、训练过程和评估指标在内的完整代码库，可在这个仓库中找到，以实现可重复性和进一步的实验。以下是结果：

实验结果展示了 KDE 和 MAF 在相对性能上的显著差异！如图表所示，在第五维左右发生转变。低于这个阈值，KDE 表现出更好的性能，然而，超过五个维度，MAF 开始以越来越大的幅度大幅超越 KDE。

在维度 7 时，这种差异的幅度变得非常明显，我们的结果表明在数据效率方面存在深刻的差异。在所有三种分布类型测试中，KDE 始终需要超过 100,000 个数据点才能达到令人满意的性能。相比之下，MAF 在所有分布中仅需最多 2,000 个数据点就达到了相同的性能阈值。这代表着从 50 倍到 100 倍的改进因子！

除了样本效率之外，计算性能的差异同样引人注目，因为在这些更高维度上，KDE 的训练时间大约是 MAF 的 12 倍。

优秀的数据处理效率和更快的训练时间使得 MAF 在处理高维任务时成为当之无愧的赢家。KDE 对于低维问题仍然是一个非常有价值的工具，但如果你正在处理涉及超过五个维度的应用，我强烈建议尝试 MAF 方法。

---

第四部分：为什么 MAF 能击败 KDE？

要理解为什么 KDE 在高维空间中表现不佳，我们首先必须检查 KDE 实际上是如何工作的。正如我在之前的文章中讨论的，核密度估计使用局部邻域估计，对于任何我们想要评估密度的点，KDE 会查看附近的数据点，并使用它们的邻近性来估计局部概率密度。每个核函数在每个数据点周围创建一个邻域，任何位置的密度估计是所有包含该位置的核函数贡献的总和。

这种局部方法在低维空间中表现良好。然而，随着维度的增加，数据变得稀疏，导致估计器需要指数级更多的数据点来填充相同密度的空间。

相反，MAF 不使用基于局部邻域的估计。MAF 不是通过查看附近的点来估计密度，而是学习将先前变量映射到条件分布参数的函数。神经网络的权重在整个输入空间中共享，这使得它可以从训练数据中泛化，而无需填充局部邻域。这种架构差异使得 MAF 在维度上的扩展能力远优于 KDE。

这种局部和全局方法之间的区别解释了我实验中观察到的显著性能差距。虽然 KDE 必须用数据点填充指数级扩展的空间以保持准确的局部邻域，但 MAF 可以利用神经网络的组合结构，使用远少于样本的数量来学习全局模式。

结论

核密度估计器在低维数据中的非参数数据分析方面表现优秀；它直观、快速，并且需要很少的调整。然而，对于高维数据或当计算时间成为关注点时，我建议尝试正常化流。虽然该模型还没有像 KDE 那样经过实战检验，但它们是一个值得尝试的可靠替代品，可能最终会成为你新的首选密度估计器。

除非另有说明，所有图像均由作者提供。主要实验的代码位于此 仓库*。 *