Numpy多元正态分布向量化

Question

我有两个二维numpy数组A和B。我想使用scipy.stats.multivariate_normal来计算每行A的联合logpdf，使用每行B作为协方差矩阵。是否有一种不需要显式循环遍历行的方法？将scipy.stats.multivariate_normal直接应用于A和B会计算出A中每行的logpdf（这正是我想要的），但使用整个2D数组A作为协方差矩阵，而这不是我想要的（我需要每行B创建一个不同的协方差矩阵）。我正在寻找一种使用numpy向量化并避免显式循环遍历两个数组的解决方案。

- W.Ali

---

你现在的代码是什么样子？ - FHTMitchell

scipy.stats.multivariate_normal( cov=B).logpdf(A) - W.Ali

我认为你必须循环，因为我所能想到的唯一替代方案就是传入一个三维数组（协方差矩阵堆栈）。然而，这将导致错误* ValueError：数组“cov”必须至多为二维*。 - MB-F

如果numpy函数（或scipy）不能处理更高维度的输入，那么避免显式循环的numpy向量化就会变得困难。 - hpaulj

是的，我也这么想，只是想确认一下我没有漏掉什么显而易见的东西。 - W.Ali

Answer 1

我也试图实现类似的功能。以下是我的代码，它接收三个NxD矩阵。X的每一行是一个数据点，means的每一行都是均值向量，covariances的每一行都是对角协方差矩阵的对角线向量。结果是一个长度为N的对数概率向量。

def vectorized_gaussian_logpdf(X, means, covariances):
    """
    Compute log N(x_i; mu_i, sigma_i) for each x_i, mu_i, sigma_i
    Args:
        X : shape (n, d)
            Data points
        means : shape (n, d)
            Mean vectors
        covariances : shape (n, d)
            Diagonal covariance matrices
    Returns:
        logpdfs : shape (n,)
            Log probabilities
    """
    _, d = X.shape
    constant = d * np.log(2 * np.pi)
    log_determinants = np.log(np.prod(covariances, axis=1))
    deviations = X - means
    inverses = 1 / covariances
    return -0.5 * (constant + log_determinants +
        np.sum(deviations * inverses * deviations, axis=1))

请注意，此代码仅适用于对角协方差矩阵。在这种特殊情况下，以下数学定义被简化：行列式变为元素乘积，逆变为逐元素倒数，矩阵乘法变为逐元素乘法。

多元正态分布概率密度函数的图片如下：

进行正确性和运行时间的快速测试：

def test_vectorized_gaussian_logpdf():
    n = 128**2
    d = 64

    means = np.random.uniform(-1, 1, (n, d))
    covariances = np.random.uniform(0, 2, (n, d))
    X = np.random.uniform(-1, 1, (n, d))

    refs = []

    ref_start = time.time()
    for x, mean, covariance in zip(X, means, covariances):
        refs.append(scipy.stats.multivariate_normal.logpdf(x, mean, covariance))
    ref_time = time.time() - ref_start

    fast_start = time.time()
    results = vectorized_gaussian_logpdf(X, means, covariances)
    fast_time = time.time() - fast_start

    print("Reference time:", ref_time)
    print("Vectorized time:", fast_time)
    print("Speedup:", ref_time / fast_time)

    assert np.allclose(results, refs)

我获得了大约250倍的加速。(是的，我的应用需要计算16384个不同的高斯函数。)

- Daniel W