1

我正在尝试使用 GPML 书中给出的 NumPy 实现 ARD 内核（来自公式 5.2 的 M3）。 在此处输入图像描述

我正在努力为 NxM 内核计算向量化这个方程。我尝试了以下非矢量化版本。有人可以帮助在 NumPy/PyTorch 中对此进行矢量化吗？

import numpy as np
N = 30 # Number of data points in X1
M = 40 # Number of data points in X2
D = 6  # Number of features (ARD dimensions)

X1 = np.random.rand(N, D)
X2 = np.random.rand(M, D)
Lambda = np.random.rand(D, 1)
L_inv = np.diag(np.random.rand(D))
sigma_f = np.random.rand()

K = np.empty((N, M))
for n in range(N):
  for m in range(M):
    M3 = Lambda@Lambda.T + L_inv**2
    d = (X1[n,:] - X2[m,:]).reshape(-1,1)
    K[n, m] = sigma_f**2 * np.exp(-0.5 * d.T@M3@d)
4

2 回答 2

1

我们可以使用广播规则和简洁的 NumPy 函数einsum来向量化数组操作。简而言之，广播允许我们通过向结果数组添加新维度来单行操作数组，同时einsum允许我们通过显式使用索引表示法（而不是矩阵）对多个数组执行操作。

幸运的是，计算内核不需要循环。请参阅下面的矢量化解决方案ARD_kernelfunction，它在我的机器中比原始的 loopy 版本快 30 倍。现在，einsum通常是尽可能快，但可能有更快的方法，但我没有检查其他任何东西（例如，通常@的运算符而不是einsum）。

此外，代码中缺少一个术语（Kronecker delta），我不知道它是否故意省略（如果您在实现它时遇到问题，请告诉我，我会编辑答案）。

import numpy as np
N = 300 # Number of data points in X1
M = 400 # Number of data points in X2
D = 6  # Number of features (ARD dimensions)

np.random.seed(1)  # Fix random seed for reproducibility
X1 = np.random.rand(N, D)
X2 = np.random.rand(M, D)
Lambda = np.random.rand(D, 1)
L_inv = np.diag(np.random.rand(D))
sigma_f = np.random.rand()

# Loopy function
def ARD_kernel_loops(X1, X2, Lambda, L_inv, sigma_f):
    K = np.empty((N, M))
    M3 = Lambda@Lambda.T + L_inv**2
    for n in range(N):
        for m in range(M):
            d = (X1[n,:] - X2[m,:]).reshape(-1,1)
            K[n, m] = np.exp(-0.5 * d.T@M3@d)
    return K * sigma_f**2
    
# Vectorized function
def ARD_kernel(X1, X2, Lambda, L_inv, sigma_f):
    M3 = Lambda.squeeze()*Lambda + L_inv**2  # Use broadcasting to avoid transpose
    d = X1[:,None] - X2[None,...]            #  Use broadcasting to avoid loops
    
    # order=F for memory layout (as your arrays are (N,M,D) instead of (D,N,M))
    return sigma_f**2 * np.exp(-0.5 * np.einsum("ijk,kl,ijl->ij", d, M3, d, order = 'F'))
于 2021-06-16T22:19:30.340 回答
0

也许还有一个额外的优化。给出的 M 个矩阵的例子都是正定的。这意味着可以应用 Cholesky 分解，我们可以找到上三角 U 使得

M = U'*U

重点是，如果我们将 U 应用于 xs，那么

y[p] = U*x[p] p=1..

然后

(x[p]-x[q])'*M*(x[p]-x[q]) = (y[p]-y[q])'*(y[p]-y[q])

因此，如果有 N 个向量 x 每个维度 d，我们将 LHS 上的 N 平方 O(d 平方) 操作转换为 RHS 上的 N 平方 O(d) 操作

这花费了额外的 Choleski 分解 (O(d 立方)) 和 NO(d 平方) U 对 xs 的应用。

于 2021-06-20T09:05:15.867 回答