reasoning
We need to write a custom CUDA kernel that computes the SVGD update for a set of particles. The reference model is fairly simple: compute pairwise squared distances, median heuristic for bandwidth h, RBF kernel K, grad_K, then phi = (K @ score + grad_K.sum(dim=0)) / N. Score is -x.
We need to produce a kernel in CUDA using PyTorch's cpp_extension.load_inline. The kernel must accept a tensor of shape (N, dim), compute phi of shape (N, dim). The precision is fp32.
We need to be faster than PyTorch reference. Let's understand the operations:
N = 1024, dim = 32. That's moderate size. PyTorch does a bunch of broadcasting and creates large intermediate tensors: diff is (N,N,dim) = (1024,1024,32) ~ 134M elements, each float32 is 512MB. That's huge, may be memory-bound. Our custom kernel can compute pairwise interactions on the fly without storing huge diff or K matrices. So there's a lot of opportunity.
We must compute median of all pairwise squared distances. But since N=1024, N*N=1,048,576, which is manageable. However, in CUDA we could compute all dist_sq and then find median. But computing all pairwise squared distances explicitly is also O(N^2) memory and compute. For N=1024, we can compute dist_sq into a buffer of size N*N, then compute median using something like torch.median (or a custom kernel). But we want speed. The PyTorch reference computes diff (N,N,dim) which uses a lot of memory and compute.
Better approach: compute pairwise squared distances efficiently, perhaps using matrix multiplication: dist_sq = ||x_i||^2 + ||x_j||^2 - 2 x_i·x_j. That requires (N,dim) input, compute norms (N), then dist_sq = norms_i + norms_j^T - 2 * x @ x.T. This is standard and avoids storing diff. Then compute median of dist_sq (flattened). Then compute K and grad_K on the fly in a kernel.
Since N is only 1024, we could also compute everything in one kernel: each thread block computes a tile of output phi for a subset of i particles, iterating over j. For each j, compute diff = x_j - x_i, dist_sq = dot(diff, diff). Compute K = exp(-dist_sq / (2h)). Compute term = K * score_j + (-diff/h) * K. Accumulate sum. At the end divide by N.
This single kernel avoids storing K and grad_K completely. It also avoids the large diff tensor. The only challenge is computing the median of dist_sq. But we can compute dist_sq using the matmul approach, store it in a N*N temporary (1M floats = 4MB), find median on GPU or CPU. 4MB is tiny. So we can compute dist_sq via efficient cuBLAS (torch.mm). Actually we can also compute dist_sq in the same kernel by first computing pairwise distances fully, but that requires materializing them or doing two passes. Two passes: first compute dist_sq matrix, find median. Second pass: compute phi with the median.
Simplify: Use PyTorch operations to compute dist_sq efficiently: norms = (particles**2).sum(dim=1, keepdim=True) -> (N,1). dist_sq = norms + norms.t() - 2 * particles @ particles.t(). This yields (N,N). Then median = dist_sq.view(-1).median().clamp(min=1e-5). This will be fast and use only O(N^2) memory for dist_sq (4MB). The heavy diff tensor (512MB) is avoided. Then we can compute phi using a custom kernel that takes particles, score (=-particles), h, and output phi. That's possible. But maybe even the phi computation can be done efficiently in PyTorch without custom kernel? Let's see.
Phi = (K @ score + grad_K.sum(dim=0)) / N.
K = exp(-dist_sq/(2h)). That's (N,N). K @ score is (N,N) @ (N,dim) = (N,dim). grad_K = -diff/h * K.unsqueeze(-1). diff = particles_i - particles_j? Wait diff in reference is particles.unsqueeze(0) - particles.unsqueeze(1) = x_j? Actually unsqueeze(0) is (1,N,dim), unsqueeze(1) is (N,1,dim). diff = particles[0,j] - particles[i,0]? Let's compute: diff_{i,j} = particles_j - particles_i? Wait (1,N,dim) - (N,1,dim) with broadcasting: result shape (N,N,dim), element (i,j,d) = particles[0,j,d] - particles[i,0,d] = x_j - x_i. But later they compute grad_K = -diff / h * K.unsqueeze(-1). And phi = (K @ score + grad_K.sum(dim=0)) / N. Let's examine the formula: SVGD phi(x_i) = (1/N) sum_j [ K(x_j, x_i) * score(x_j) + grad_K(x_j, x_i) ]. Typically grad_K(x_j, x_i) = ∇_{x_i} K(x_j, x_i). For RBF K(x_j, x_i) = exp(-||x_j - x_i||^2/(2h)). Then ∇_{x_i} K = -(x_i - x_j)/h * K = (x_j - x_i)/h * K. Wait reference uses grad_K = -diff/h * K, where diff = x_j - x_i. So grad_K = -(x_j - x_i)/h * K = (x_i - x_j)/h * K. That matches ∇_{x_i} K = -(x_j - x_i)/h * K? Actually derivative w.r.t x_i of -||x_j - x_i||^2 is 2(x_j - x_i). Let's compute: K = exp(-||x_j - x_i||^2/(2h)). ∂K/∂x_i = exp(...)*(-1/(2h))*2*(x_i - x_j) = -(x_i - x_j)/h * K = (x_j - x_i)/h * K. So grad_K should be (x_j - x_i)/h * K. But reference says grad_K = -diff/h * K, and diff = x_j - x_i. So grad_K = -(x_j - x_i)/h * K = (x_i - x_j)/h * K. That seems off by a sign? Let's check the reference code carefully.
diff = particles.unsqueeze(0) - particles.unsqueeze(1)
# particles.unsqueeze(0) shape (1,N,dim)
# particles.unsqueeze(1) shape (N,1,dim)
# Broadcasting: diff[i,j] = particles[0,j] - particles[i,0] = x_j - x_i.
grad_K = -diff / h * K.unsqueeze(-1) = -(x_j - x_i)/h * K = (x_i - x_j)/h * K.
But the SVGD formula uses grad_K(x_j, x_i) = ∇_{x_i} K(x_j, x_i). That should be (x_i - x_j)/h * K? Let's compute derivative w.r.t x_i of r = x_j - x_i. Wait the argument is -||r||^2/(2h) where r = x_j - x_i. ∂/∂x_i (-||r||^2) = -2 * r * ∂r/∂x_i = -2 * r * (-1) = 2r. So ∂K/∂x_i = K * (1/(2h)) * 2r = r/h * K = (x_j - x_i)/h * K. Wait this is positive. Let's recompute: K = exp(-||x_j - x_i||^2/(2h)). Let f = -||x_j - x_i||^2/(2h). ∂f/∂x_i = -(1/(2h)) * 2 * (x_i - x_j) = -(x_i - x_j)/h = (x_j - x_i)/h. So ∂K/∂x_i = K * (x_j - x_i)/h. So grad_K = (x_j - x_i)/h * K. But reference has grad_K = -diff/h * K = -(x_j - x_i)/h * K. That is the negative of the correct derivative? Let's double-check diff. Maybe unsqueeze(0) - unsqueeze(1) yields diff[i,j] = particles[i] - particles[j]? Let's verify.
particles shape (N,dim).
particles.unsqueeze(0) shape (1,N,dim).
particles.unsqueeze(1) shape (N,1,dim).
When broadcasting subtract, PyTorch aligns from trailing dimensions:
- Rightmost: dim matches.
- Middle: N vs 1 -> broadcast to N.
- Leftmost: 1 vs N -> broadcast to N.
So shape is (N,N,dim).
Element (i,j,d) = particles[0,j,d] - particles[i,0,d] = x_j - x_i.
Yes, diff[i,j] = x_j - x_i.
Then grad_K = -diff/h * K = -(x_j - x_i)/h * K = (x_i - x_j)/h * K.
But SVGD formula phi(x_i) = (1/N) sum_j [ K(x_j, x_i) * score(x_j) + ∇_{x_i} K(x_j, x_i) ].
∇_{x_i} K = (x_j - x_i)/h * K (as derived). But reference uses negative of that. Wait maybe the formula they implement is phi(x_i) = (1/N) sum_j [ K(x_j, x_i) * score(x_j) + ∇_{x_j} K(x_j, x_i) ]? Let's check the original SVGD paper: phi(x) = (1/n) sum_i [ k(x_i, x) ∇_{x_i} log p(x_i) + ∇_{x_i} k(x_i, x) ]. Here first argument is x_i, second is x. So ∇_{x_i} k(x_i, x) = derivative w.r.t first argument. If we set x = x_i (target), x_i = x_j (source), then ∇_{x_j} k(x_j, x_i). k(x_j, x_i) = exp(-||x_j - x_i||^2/(2h)). ∂/∂x_j k = -(x_j - x_i)/h * k. That matches reference: grad_K = -diff/h * K where diff = x_j - x_i. So reference's grad_K is derivative w.r.t the first argument (the source particle x_j). And the sum is over j of that. But the formula in the docstring says "phi(x_i) = (1/N) sum_j [ K(x_j, x_i) * score(x_j) + grad_K(x_j, x_i) ]". It doesn't specify derivative w.r.t which argument, but the code computes derivative w.r.t x_j? Actually diff = x_j - x_i, and they take -diff/h * K. That's derivative w.r.t x_j. Wait if f(x_j) = exp(-||x_j - x_i||^2/(2h)), ∂f/∂x_j = -(x_j - x_i)/h * f. Yes, that's exactly -diff/h * K. So grad_K is ∇_{x_j} K. But the term in SVGD is ∇_{x_j} K(x_j, x_i) (derivative w.r.t the source particle). So the code is correct.
But then phi = (K @ score + grad_K.sum(dim=0)) / N. Let's check dimensions: K is (N,N), score is (N,dim). K @ score gives (N,dim) where row i = sum_j K_{i,j} score_j? Wait K @ score: matrix multiplication: result_{i,d} = sum_j K_{i,j} * score_{j,d}. So it's summing over j with K indexed by (i,j). In the formula, for phi(x_i), we need sum_j K(x_j, x_i) score(x_j). K(x_j, x_i) corresponds to K_{i,j} if rows correspond to x_i and columns to x_j? In reference, dist_sq[i,j] = ||x_j - x_i||^2. Then K[i,j] = exp(-dist_sq[i,j]/(2h)) = K(x_j, x_i). Then K @ score: row i = sum_j K_{i,j} score_j = sum_j K(x_j, x_i) score(x_j). Good.
grad_K shape (N,N,dim). grad_K[i,j] = ∇_{x_j} K(x_j, x_i) = -(x_j - x_i)/h * K(x_j, x_i). Then grad_K.sum(dim=0) sums over i: sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_j, x_i). Wait in the formula phi(x_i) = sum_j [ ... + grad_K(x_j, x_i) ]. That suggests summing over j. But grad_K.sum(dim=0) sums over the first dimension (i), yielding a vector indexed by j: sum_i grad_K[i,j]. That gives (N,dim). But phi should be indexed by i. Wait: If grad_K is (N,N,dim), dim=0 is i (target), dim=1 is j (source). grad_K.sum(dim=0) gives sum over i for each j. That is sum_i ∇_{x_j} K(x_j, x_i). But according to the docstring, phi should be sum_j grad_K(x_j, x_i). Let's check if there's a symmetry. Since K(x_j, x_i) = K(x_i, x_j), and ∇_{x_j} K(x_j, x_i) = -(x_j - x_i)/h * K. Meanwhile ∇_{x_i} K(x_j, x_i) = (x_j - x_i)/h * K = -∇_{x_j} K(x_j, x_i). So grad_K.sum(dim=0) = sum_i ∇_{x_j} K(x_j, x_i) = -sum_i (x_j - x_i)/h * K. Is this equal to sum_j ∇_{x_i} K(x_j, x_i)? Let's see: sum_j ∇_{x_i} K(x_j, x_i) = sum_j (x_j - x_i)/h * K = -sum_j (x_i - x_j)/h * K. Since K is symmetric and the sum runs over all j, it's just the negative of sum_j ∇_{x_j} K(x_j, x_i) = -grad_K.sum(dim=1). Wait reference uses grad_K.sum(dim=0). That seems suspicious.
Let's test with small N to see if reference code is actually consistent. Let's pick N=2, dim=1. x = [a,b]. diff = [[a-a, b-a], [a-b, b-b]] = [[0, b-a], [a-b, 0]]. dist_sq = [[0, (b-a)^2], [(a-b)^2, 0]]. K = [[1, k], [k, 1]] where k = exp(-(b-a)^2/(2h)). score = [-a, -b]. grad_K = -diff/h * K.unsqueeze(-1). So grad_K[0,0]=0, grad_K[0,1] = -(b-a)/h * k, grad_K[1,0] = -(a-b)/h * k = (b-a)/h * k, grad_K[1,1]=0.
grad_K.sum(dim=0) = [0 + (b-a)/h * k, -(b-a)/h * k + 0] = [(b-a)/h * k, -(b-a)/h * k].
K @ score = [1*(-a) + k*(-b), k*(-a) + 1*(-b)] = [-a - k b, -k a - b].
phi = (K @ score + grad_K.sum(dim=0)) / 2.
phi_0 = (-a - k b + (b-a)/h * k) / 2.
phi_1 = (-k a - b - (b-a)/h * k) / 2.
Now let's compute the formula manually: phi(x_i) = (1/2) sum_j [ K(x_j, x_i) score(x_j) + ∇_{x_j} K(x_j, x_i) ].
For i=0 (x_i = a):
j=0: K(a,a)*score(a) + ∇_a K(a,a) = 1*(-a) + 0 = -a.
j=1: K(b,a)*score(b) + ∇_b K(b,a) = k*(-b) + (-(b-a)/h * k) = -k b - (b-a)/h * k.
Sum = -a - k b - (b-a)/h * k.
phi_0 = sum/2 = (-a - k b - (b-a)/h * k)/2.
But reference phi_0 = (-a - k b + (b-a)/h * k)/2. So reference has a sign difference in the grad term! Reference's grad_K.sum(dim=0) gave + (b-a)/h * k for phi_0, but manual gives - (b-a)/h * k.
Wait maybe the formula in reference uses ∇_{x_i} K(x_j, x_i) not ∇_{x_j} K(x_j, x_i). Let's check: If phi(x_i) = (1/N) sum_j [ K(x_j, x_i) score(x_j) + ∇_{x_i} K(x_j, x_i) ], then for i=0:
j=1 term: ∇_a K(b,a) = derivative of K w.r.t second argument a. Since K(b,a) = exp(-(b-a)^2/(2h)), derivative w.r.t a = (b-a)/h * k. So j=1 term = -k b + (b-a)/h * k. Sum = -a + (-k b + (b-a)/h * k). phi_0 = (-a - k b + (b-a)/h * k)/2. That matches reference.
So reference's grad_K is actually ∇_{x_i} K(x_j, x_i)? But they computed grad_K = -diff/h * K with diff = x_j - x_i. ∇_{x_i} K = (x_j - x_i)/h * K. Wait that's positive diff/h * K. But they have -diff/h * K. So that would be -(x_j - x_i)/h * K = (x_i - x_j)/h * K. That's not ∇_{x_i} K. Wait what's going on?
Let's recompute derivative of K(x_j, x_i) w.r.t x_i. K = exp(-||x_j - x_i||^2/(2h)). Let r = x_j - x_i. ∂/∂x_i (-||r||^2) = -2 * r * (∂r/∂x_i) = -2 * r * (-1) = 2r. Divide by 2h gives r/h. So ∂K/∂x_i = K * (r/h) = (x_j - x_i)/h * K. So ∇_{x_i} K = diff/h * K.
Reference sets grad_K = -diff/h * K, which is -(x_j - x_i)/h * K. That's negative of ∇_{x_i} K. But earlier we saw that if phi formula includes ∇_{x_i} K, then reference's sign is wrong. Wait phi_0 manual with ∇_{x_i} K gave + (b-a)/h * k. Reference also gave + (b-a)/h * k for phi_0? Let's re-evaluate reference phi_0: phi_0 = (-a - k b + (b-a)/h * k)/2. Yes, plus sign. So reference's grad_K = -diff/h * K contributes + (b-a)/h * k to phi_0 because diff = b-a for (i=0,j=1). So grad_K[0,1] = -(b-a)/h * k. But phi_0 includes sum over j of grad_K[0,j]? Wait K @ score yields sum over j of K[0,j] score_j. grad_K.sum(dim=0) yields sum over i of grad_K[i,j] for each j. Let's add K @ score[0] = sum_j K[0,j] score_j = K[0,0]*score_0 + K[0,1]*score_1 = -a - k b. grad_K.sum(dim=0)[0] = grad_K[0,0] + grad_K[1,0] = 0 + (-(a-b)/h * k) = (b-a)/h * k. So phi_0 = (-a - k b + (b-a)/h * k)/2. So the gradient term added to phi_0 is grad_K.sum(dim=0)[0], which is sum over i of grad_K[i,0] = sum over i of -(x_0 - x_i)/h * K[i,0] = sum_i (x_i - x_0)/h * K[i,0]. But the formula says sum_j grad_K(x_j, x_0). If grad_K(x_j, x_i) is defined as ∇_{x_i} K(x_j, x_i), then grad_K[j,i] = ∇_{x_i} K(x_j, x_i) = (x_j - x_i)/h * K. But reference defines grad_K[i,j] = -diff[i,j]/h * K = -(x_j - x_i)/h * K. So grad_K[i,j] = ∇_{x_j} K(x_j, x_i)? No, ∇_{x_j} K = -(x_j - x_i)/h * K. So grad_K[i,j] = ∇_{x_j} K(x_j, x_i). Then grad_K.sum(dim=0)[j] = sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_j, x_i). For j=0, that's sum_i ∇_{x_0} K(x_0, x_i) = sum_i -(x_0 - x_i)/h * K(x_0, x_i) = sum_i (x_i - x_0)/h * K(x_0, x_i). This is equal to sum_i ∇_{x_i} K(x_0, x_i)? No, ∇_{x_i} K(x_0, x_i) = (x_0 - x_i)/h * K. So it's the negative. So reference's phi includes sum_i ∇_{x_j} K(x_j, x_i) evaluated at j being the output index (since grad_K.sum(dim=0) index is j). But the SVGD formula is sum_j [ K(x_j, x_i) score(x_j) + ∇_{x_j} K(x_j, x_i) ]. For output index i, we sum over j. But reference sums over i (first dim) and outputs j (second dim). Because of symmetry of summation over all particles, sum_j ∇_{x_j} K(x_j, x_i) = sum_j ∇_{x_i} K(x_j, x_i) * (-1)? Actually due to symmetry of the kernel, maybe sum_j ∇_{x_j} K(x_j, x_i) = -sum_j ∇_{x_i} K(x_j, x_i). But in SVGD, the repulsion term is ∇_{x_j} K(x_j, x_i). Wait the standard formula uses ∇_{x_j} K(x_j, x_i), because the update is (1/n) sum_j [ k(x_j, x) ∇log p(x_j) + ∇_{x_j} k(x_j, x) ]. Here x is the target point (x_i). So the gradient is w.r.t the first argument x_j. So reference's grad_K = ∇_{x_j} K(x_j, x_i) = -(x_j - x_i)/h * K. That's correct.
Then phi for x_i should be sum_j grad_K[i?]. grad_K[i,j] corresponds to ∇_{x_j} K(x_j, x_i). So for a fixed i (the target), we need sum over j of grad_K[i,j]. That would be grad_K.sum(dim=1). But reference uses grad_K.sum(dim=0). Wait let's check: grad_K shape (N,N,dim). The first dimension is i (target), second is j (source). grad_K[i,j] = ∇_{x_j} K(x_j, x_i). For target i, sum over j is sum(grad_K[i, :]). That is dim=1. But reference does grad_K.sum(dim=0). That sums over i for each j. That would give sum over targets i for source j. That's like computing something for each source j.
Is there a transpose somewhere? Let's verify with the code: phi = (K @ score + grad_K.sum(dim=0)) / N. K @ score yields (N,dim) where row i = sum_j K[i,j] score_j. grad_K.sum(dim=0) yields (N,dim) where row j = sum_i grad_K[i,j]. So phi_i = (sum_j K[i,j] score_j + sum_j' grad_K[j', i]) / N. The second term sums over the first dimension with the second dimension fixed to i. That is sum_{j'} grad_K[j', i] = sum_{j'} ∇_{x_i} K(x_i, x_{j'}). Wait grad_K[j', i] = ∇_{x_i} K(x_i, x_{j'}). Because in grad_K[i,j], first index is target, second is source. So grad_K[j', i] has target j', source i, derivative w.r.t source i. So sum_{j'} ∇_{x_i} K(x_i, x_{j'}). But K(x_i, x_{j'}) = K(x_{j'}, x_i). So derivative w.r.t first argument (source) when source is i and target is j'. That's ∇_{x_i} K(x_i, x_{j'}). The standard formula uses ∇_{x_j} K(x_j, x_i) for each j. Here we have ∇_{x_i} K(x_i, x_{j'}) for each j'. Because of symmetry, sum_{j'} ∇_{x_i} K(x_i, x_{j'}) = -sum_{j'} ∇_{x_{j'}} K(x_i, x_{j'}). But K(x_i, x_{j'}) = K(x_{j'}, x_i). So ∇_{x_{j'}} K(x_{j'}, x_i) = -(x_{j'} - x_i)/h * K = (x_i - x_{j'})/h * K. Meanwhile ∇_{x_i} K(x_i, x_{j'}) = -(x_i - x_{j'})/h * K = (x_{j'} - x_i)/h * K. So indeed ∇_{x_i} K(x_i, x_{j'}) = -∇_{x_{j'}} K(x_{j'}, x_i). So sum_{j'} grad_K[j', i] = -sum_{j'} ∇_{x_{j'}} K(x_{j'}, x_i). Thus reference's phi_i includes the negative of the correct repulsion term? That would be wrong.
But wait, the code as given is the reference implementation. We need to match it exactly. So we must compute exactly what the reference does, regardless of whether it's theoretically correct. So we need to compute:
phi_i = (sum_j K_{i,j} * score_j + sum_j grad_K_{j, i}) / N
where K_{i,j} = exp(-||x_j - x_i||^2 / (2h)).
grad_K_{i,j} = -(x_j - x_i) / h * K_{i,j}.
Note: phi_i is indexed by i.
Alternatively, we can express grad_K_{j,i} = -(x_i - x_j)/h * K_{j,i} = -(x_i - x_j)/h * K_{i,j} (since K symmetric). Then sum_j grad_K_{j,i} = sum_j -(x_i - x_j)/h * K_{i,j} = (1/h) sum_j (x_j - x_i) * K_{i,j}.
But the reference computes grad_K explicitly as -(diff)/h * K, where diff = x_j - x_i. So grad_K_{i,j} = -(x_j - x_i)/h * K_{i,j}. Then grad_K_{j,i} = -(x_i - x_j)/h * K_{j,i} = (x_j - x_i)/h * K_{i,j}. So sum_j grad_K_{j,i} = sum_j (x_j - x_i)/h * K_{i,j}.
Wait, earlier we computed for N=2:
grad_K[0,1] = -(b-a)/h * k.
grad_K[1,0] = -(a-b)/h * k = (b-a)/h * k.
grad_K.sum(dim=0)[0] = grad_K[0,0] + grad_K[1,0] = 0 + (b-a)/h * k.
This equals sum_j grad_K[j,0] = grad_K[0,0] + grad_K[1,0] = (b-a)/h * k.
But grad_K[j,0] = -(x_0 - x_j)/h * K_{j,0} = -(a - x_j)/h * K_{j,0}. For j=1: -(a-b)/h * k = (b-a)/h * k. Yes.
So sum_j grad_K[j,i] = sum_j -(x_i - x_j)/h * K_{j,i} = sum_j (x_j - x_i)/h * K_{i,j}.
That is the same as (1/h) sum_j (x_j - x_i) * K_{i,j}.
But phi_i = (sum_j K_{i,j} score_j + sum_j grad_K[j,i]) / N.
So we need to compute for each i:
phi_i = (1/N) * sum_j [ K_{i,j} * score_j + (x_j - x_i)/h * K_{i,j} ].
Wait is that equal to (1/N) * sum_j K_{i,j} * (score_j + (x_j - x_i)/h) ?
Yes: phi_i = (1/N) * sum_j K_{i,j} * (score_j + (x_j - x_i)/h).
But reference computes grad_K_{i,j} = -(x_j - x_i)/h * K_{i,j} and then sums over the first dimension with the second fixed (i.e., sum_j grad_K[j,i] = sum_j -(x_i - x_j)/h * K_{j,i} = sum_j (x_j - x_i)/h * K_{i,j}). So indeed phi_i = (1/N) * ( sum_j K_{i,j} score_j + sum_j (x_j - x_i)/h * K_{i,j} ).
Thus if we compute diff_{i,j} = x_j - x_i, then K_{i,j} = exp(-||diff||^2/(2h)), and the update is sum_j K_{i,j} * (score_j + diff_{i,j}/h) / N.
But wait reference uses K @ score + grad_K.sum(dim=0). Let's expand grad_K.sum(dim=0): grad_K[i,j] = -(x_j - x_i)/h * K_{i,j}. grad_K.sum(dim=0)[j] = sum_i grad_K[i,j] = sum_i -(x_j - x_i)/h * K_{i,j}. That does not look like our needed term. Wait careful: In PyTorch, grad_K.sum(dim=0) sums over dimension 0 (size N), leaving dimension 1 (size N). But shape is (N,N,dim). So sum(dim=0) yields (N,dim) where index j corresponds to original index j: result_j = sum_i grad_K[i,j,d]. That is sum_i -(x_j - x_i)/h * K_{i,j}. Since K_{i,j} = K_{j,i}, we can rewrite as sum_i -(x_j - x_i)/h * K_{j,i}. But we need for phi_i sum_j grad_K[j,i] = sum_j -(x_i - x_j)/h * K_{j,i} = sum_j (x_j - x_i)/h * K_{i,j}. Wait K_{j,i} = K_{i,j}. So sum_j grad_K[j,i] = sum_j -(x_i - x_j)/h * K_{i,j} = sum_j (x_j - x_i)/h * K_{i,j}. Yes.
But grad_K.sum(dim=0) gives for each j: sum_i -(x_j - x_i)/h * K_{i,j}. This is a vector indexed by j, not i. However, because the sum is over all i and j is the output index, this is not the same as phi_i. Wait but in the reference code they add K @ score (N,dim) + grad_K.sum(dim=0) (N,dim). K @ score's row i is sum_j K_{i,j} score_j. grad_K.sum(dim=0)'s row j is sum_i grad_K[i,j]. So phi_i = (sum_j K_{i,j} score_j + sum_j' grad_K[j', i]) / N? Let's check indexing: grad_K.sum(dim=0) returns a tensor of shape (N,dim) where the first dimension corresponds to the second dimension of grad_K. So grad_K.sum(dim=0)[i] = sum_j grad_K[j, i]? Wait no: grad_K shape (N,N,dim). Let's label dims: dim0 = i (target), dim1 = j (source). grad_K[i,j] = something. sum(dim=0) means sum over i, leaving j. So result[j] = sum_i grad_K[i,j]. The output index is j. But PyTorch's matrix multiplication K @ score: K shape (N,N), score shape (N,dim). Output row i = sum_j K[i,j] score_j. The output index is i. So when they add K @ score + grad_K.sum(dim=0), they add two (N,dim) tensors where the first output index i of K@score corresponds to the target i, but the first output index i of grad_K.sum(dim=0) corresponds to source j=i. So they are adding sum_j K[i,j] score_j + sum_{i'} grad_K[i', i]. This is exactly phi_i = sum_j K[i,j] score_j + sum_{i'} grad_K[i', i]. And grad_K[i', i] = -(x_i - x_{i'})/h * K_{i',i} = (x_{i'} - x_i)/h * K_{i,i'}. So phi_i = sum_j K[i,j] score_j + sum_j (x_j - x_i)/h * K[i,j] all over N. So indeed phi_i uses the output index i for both terms. The second term sums over the first dimension of grad_K with the second dimension fixed to i, which is sum_{j} grad_K[j,i]. And grad_K.sum(dim=0)[i] = sum_j grad_K[j,i]. Wait earlier I said grad_K.sum(dim=0) sums over dim0 leaving dim1, so output index is dim1. If dim1 is i, then output[i] = sum_j grad_K[j,i]. Yes! I just got confused with variable naming: In grad_K, dim0 = i, dim1 = j. sum(dim=0) => sum over i, output indexed by j. But in the code, they add it to K@score which is indexed by i. So the output of grad_K.sum(dim=0) must be aligned with K@score such that j corresponds to i. Since both have size N, it's fine: phi_i = (K@score)[i] + (grad_K.sum(dim=0))[i]. But (grad_K.sum(dim=0))[i] = sum_j grad_K[j, i]? Wait if dim0 is the first dimension (size N) and dim1 is second dimension (size N), then grad_K[j,i] means first index j, second index i. sum(dim=0) means sum over the first index j, leaving second index i. So result[i] = sum_j grad_K[j,i]. Yes. Because the first dimension is the one being summed. In my earlier label, I assigned dim0 = target i, but in the tensor the dimensions are just 0 and 1. So grad_K[j,i] corresponds to target j, source i? Wait no: grad_K[i,j] = -(x_j - x_i)/h * K[i,j]. So target is i, source is j. Then grad_K[j,i] = -(x_i - x_j)/h * K[j,i]. So grad_K.sum(dim=0)[i] = sum_j grad_K[j,i] = sum_j -(x_i - x_j)/h * K[j,i] = sum_j (x_j - x_i)/h * K[i,j]. This matches what we need.
So the formula is:
phi_i = (1/N) * sum_j K_{i,j} * (score_j + (x_j - x_i)/h).
Where score_j = -x_j.
Thus phi_i = (1/N) * sum_j exp(-||x_j - x_i||^2/(2h)) * ( -x_j + (x_j - x_i)/h ).
Simplify term inside: -x_j + (x_j - x_i)/h = x_j * (-1 + 1/h) - x_i / h.
= x_j * ((1-h)/h) - x_i / h.
But we don't need to simplify.
Now we need to implement this efficiently in CUDA.
Because N=1024, dim=32. We can process each i with a thread block. Each thread block iterates over j. Since dim=32, we can have 32 threads per block (or maybe 128 for memory coalescing). Each thread can handle one dimension index d. For each j, we load x_j (dim floats) into shared memory or registers. Then each thread computes diff_d = x_j[d] - x_i[d], accumulates dist_sq += diff_d^2. After loading all dims, we compute K = exp(-dist_sq/(2h)). Then compute term = K * (score_j[d] + diff_d / h). Actually score_j[d] = -x_j[d]. So term = K * (-x_j[d] + diff_d / h). Each thread accumulates sum_d += term. After iterating over all j, divide by N and write to output[i,d].
This is straightforward. We need to compute h first, which requires median of dist_sq.
Compute dist_sq matrix: we can use PyTorch's efficient formula: norms = (particles**2).sum(dim=1, keepdim=True) -> (N,1). dist_sq = norms + norms.t() - 2 * particles @ particles.t(). That yields (N,N). Then compute median = dist_sq.view(-1).median().clamp(1e-5). Then h = median / max(log(N), 1.0). Note: N is float. In reference: h = median_sq / max(torch.log(torch.tensor(float(N))).item(), 1.0). So h is a Python float. We can compute it on CPU.
Since N=1024, the dist_sq matrix is 1M elements, 4MB. Computing it with cuBLAS is very fast. The overhead of this PyTorch op before launching our custom kernel is small. But could we compute median in the kernel as well? It would require two passes. But computing dist_sq matrix and median in PyTorch is simple and likely fast enough. However, our custom kernel can just take particles and h and compute phi. That avoids the large intermediate diff and grad_K tensors.
But maybe we can compute dist_sq more efficiently in the kernel without materializing the full matrix. Since N=1024, we can compute all pairwise distances in shared memory block by block. But to compute the median, we need the full set of N^2 distances. Sorting 1M elements on GPU is not trivial but possible. However, PyTorch's median on GPU is optimized. So it's fine to use PyTorch for dist_sq and median.
Let's check the overall plan:
1. In ModelNew.forward, compute particles (input).
2. Compute norms = (particles ** 2).sum(dim=1, keepdim=True) # (N,1)
3. Compute dist_sq = norms + norms.t() - 2.0 * (particles @ particles.t()) # (N,N)
- We need to ensure numerical stability? Reference uses direct diff and sum of squares, so no large cancellation if numbers are moderate. The matmul approach might have slight negative values due to precision. We should clamp to 0? The reference's dist_sq are exact sum of squares of differences. With matmul, we might get small negative values. We can clamp to 0: dist_sq = torch.clamp(dist_sq, min=0.0). But does that affect median? Probably negligible. But for correctness, we should match reference as close as possible. However, computing dist_sq via matmul is standard and the difference is tiny. But to be safe, we could compute dist_sq using the same approach as reference but more memory efficient: compute pairwise distances in chunks. E.g., process particles in blocks of size 64 or 128, compute diff for block pairs, compute dist_sq, accumulate or store. Since N=1024, we can just compute the full (N,N) dist_sq using PyTorch's broadcasting but without the huge (N,N,dim) diff: we can use torch.cdist or the matmul trick.
Actually PyTorch's torch.cdist computes pairwise distances efficiently. But torch.cdist returns sqrt(dist_sq). We could use torch.cdist and square it. But torch.cdist might use a more stable algorithm. However, to match reference exactly, we need squared distances. Using matmul trick: dist_sq = -2*particles@particles.t() + norms + norms.t(). Let's test potential negative values. For random particles of size 1024x32, the norms are around 32. The dot products are around 0. dist_sq values are around 64. Small negative rounding might occur but clamping to 0 will fix it. But the reference won't have negative values because it's sum of squares of differences. The clamped version might differ for very small distances (duplicate particles). But the probability is low. We could also compute dist_sq in the custom kernel alongside phi using a two-pass approach: first compute all dist_sq in a temporary buffer (N*N), then compute median (using thrust or custom), then compute phi. This would match reference exactly. But writing a median kernel is a bit more work. However, PyTorch's median is available. Let's consider using PyTorch for the median but compute dist_sq in a custom kernel to avoid the (N,N,dim) diff. Actually we only need dist_sq for median. We can compute it using the matmul trick in PyTorch. Matmul of (1024,32) @ (32,1024) is extremely fast. Then add norms. So overall pre-processing is negligible.
But the evaluation might measure end-to-end time of forward(). So any PyTorch ops we do before launching CUDA will count. However, the reference also does many PyTorch ops: diff (N,N,dim) creation, dist_sq, median, K, grad_K, K@score, sum, division. Our approach does: norms (O(N*dim)), matmul (O(N^2*dim) but efficient), dist_sq (O(N^2)), median (O(N^2)), then custom kernel (O(N^2*dim)). Since N=1024, dim=32, the matmul is 1024*1024*32 ~ 33M multiply-adds, which is tiny. The custom kernel is also about N^2*dim = 33M operations. So total compute is dominated by the kernel and the matmul.
But can we do even better by fusing everything into a single kernel, including median computation? We could compute pairwise distances on the fly, store them in a (N,N) shared memory or global memory buffer, then do a block-wide sort to find median? That seems overkill. For N=1024, 4MB buffer is fine. We can launch a kernel to compute dist_sq into a temp buffer, then call a Thrust median (or PyTorch median on the buffer). Thrust's median requires sorting 1M elements, which takes maybe 0.1-0.2 ms. PyTorch's median on GPU is also optimized. Let's benchmark mentally: reference runtime might be dominated by allocating and computing the huge 512MB diff tensor. Our approach avoids that entirely. So we should see a big speedup.
Let's design the custom kernel.
Input: particles (N, dim) float32, h (float), output phi (N, dim) float32.
For each i (0..N-1), a block of threads processes the sum over j. Since dim=32, a natural block size is 32 (one thread per dimension). But we might want more threads to hide latency, e.g., 128 threads per block, where each thread handles a subset of j's or a subset of dimensions? Let's think.
If block size is 32 (warp). For each j, we need to compute diff across all dims to get dist_sq. Since dim=32, each thread can compute one diff and square it, then we need a warp reduction to sum the squares. This is efficient. Then each thread computes term = K * (-x_j[d] + diff_d / h) and accumulates. We iterate j from 0..N-1.
This kernel will have N blocks (one per i), 32 threads per block. Each block reads x_i (32 floats) into registers/shared memory. Then loop over j: read x_j (32 floats) from global memory. Since threads in a warp read consecutive floats (if dim is contiguous in memory? particles is (N,dim) with dim last dimension. So x_j is at offset j*dim + d. Thread d reads x_j[d]. This is coalesced: threads 0..31 read consecutive addresses. Good. Similarly x_i is read by all threads, but that's broadcast.
So the access pattern is very good. Compute dist_sq via warp shuffle reduction (or shared memory). With block size 32, we can do:
float dist_sq = diff * diff;
for (int offset = 16; offset > 0; offset /= 2) dist_sq += __shfl_down_sync(0xffffffff, dist_sq, offset);
After the loop, dist_sq is known to all threads? __shfl_down returns the value from the lane lane+offset to lane. After the loop, lane 0 has the sum. We can broadcast to all threads via __shfl_sync(0xffffffff, dist_sq, 0). Then K = expf(-dist_sq / (2*h)). Then term = K * (-x_j[d] + diff / h). Accumulate.
At the end, output[i,d] = accum / N.
Total global memory reads: N blocks * N j's * (dim floats for x_j + dim floats for score? score is -x_j, we can compute on the fly) = N^2 * dim reads of x_j. That's about 33M floats = 132MB reads. x_i is read once per block: N*dim = 32KB. Writes: N*dim = 32KB. So memory bandwidth is ~132MB. On a modern GPU with >1 TB/s bandwidth, this is negligible, runtime will be very small, maybe <0.1 ms. However, overhead of kernel launch and the matmul+median might dominate. The reference will be much slower because of the huge intermediate.
But wait, the evaluation includes the entire forward() time. If our forward() includes PyTorch ops (matmul, median) plus kernel launch, we need to ensure total time < reference time. Since reference is slow due to huge memory allocation and operations, we should be fine.
But we can also avoid the matmul+median PyTorch ops by computing dist_sq in the kernel using the same per-block approach but storing results. However, computing median on GPU is nontrivial. But we could compute dist_sq with the custom kernel into a temporary tensor, then call torch.median on that tensor. That's exactly what we plan. Actually we can compute dist_sq in the custom kernel: each block i computes dist_sq to all j, writes to a buffer. But the matmul is likely faster for dist_sq because it uses Tensor Cores or optimized BLAS. Let's see: computing dist_sq for each i,j with a separate kernel of N blocks and 32 threads: each block does N * dim reads and N writes. That's N^2 writes. But the matmul does N*dim reads for particles, then N*N reads/writes for intermediate norms and dist_sq. cuBLAS is highly optimized. For 1024x1024x32, matmul is trivial. So both are fine. Let's stick with PyTorch for dist_sq and median to keep code simple.
But we need to match reference exactly. The reference dist_sq is (diff**2).sum(-1). The matmul approach might have slightly different numerical results due to rounding. However, for random data of size 1024x32, the difference is likely within tolerance for the final phi? The tolerance for fp32 is atol=rtol=1e-4. The difference in dist_sq might be on the order of 1e-5 or less relative. But the median might shift slightly, causing h to differ. If h differs, K changes, and phi might differ more. We should test correctness.
To be safe, we can compute dist_sq in a custom kernel that exactly computes sum of squared differences, just like reference. Since N=1024, we can launch a kernel with N*N threads, each computing one dist_sq[i,j] = sum_d (x_j[d] - x_i[d])^2. Each thread iterates over d=0..31. That's 1M threads, which is a lot but manageable. Or we can use a grid of N blocks, 32 threads per block, and each block writes N dist_sq values (one per j) by computing the sum via warp shuffle. Then we write to a flat buffer of size N*N.
Actually computing dist_sq in a separate kernel might be just as fast and ensures exact match. Let's decide: use a custom kernel for dist_sq and phi, but we still need median. We can compute dist_sq into a temporary tensor using a simple kernel, then call torch.median. Then launch the main phi kernel.
Alternatively, we can compute dist_sq using the reference approach but chunked to avoid huge memory. For example, process blocks of size 64: for i_chunk in range(0,N,64): for j_chunk in range(0,N,64): diff = particles[i_chunk:i_chunk+64].unsqueeze(1) - particles[j_chunk:j_chunk+64].unsqueeze(0) -> shape (64,64,dim) = 16384*32=524k floats = 2MB. dist_sq_chunk = (diff**2).sum(-1). Store into a preallocated dist_sq matrix. This uses only 2MB intermediate and matches reference exactly. Then after all chunks, median = dist_sq.view(-1).median(). This is very safe for correctness and doesn't require a custom kernel for dist_sq. However, it introduces Python loops and many small PyTorch ops, which might have overhead. N=1024, chunk=64 -> 16x16=256 chunks. Each chunk does a subtraction, square, sum. That's 256 small kernel launches. Overhead might be non-negligible but still likely much less than reference's single 512MB allocation and ops? Actually 256 launches might be ~0.1 ms each = 25 ms, which is too slow. We could chunk larger, e.g., chunk=256: diff shape (256,256,32) = 65k floats = 256KB. Then 4x4=16 chunks. 16 launches, fine. But still overhead.
Better: compute dist_sq with a single matmul, accept possible tiny differences, and rely on tolerance. Let's test via correctness tool.
Plan:
- In forward:
- particles = input.
- N = particles.shape[0]
- Compute dist_sq using matmul trick: norms = (particles ** 2).sum(dim=1, keepdim=True); dist_sq = norms + norms.t() - 2.0 * particles @ particles.t()
- Optionally clamp: dist_sq = dist_sq.clamp(min=0.0)
- median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
- h = median_sq.item() / max(math.log(float(N)), 1.0)
- Launch custom kernel `svgd_phi` with particles and h, output phi.
Custom kernel details:
- Input: float* particles, float h, int N, int dim, float* out.
- Launch: N blocks, dim threads (or maybe 32 threads). Since dim is variable (constructor arg dim), but in the benchmark dim=32. We can handle any dim by having each thread handle one dimension if dim <= 1024. For simplicity, we can set threads = min(dim, 1024). But for reduction of dist_sq across dim, if dim > 1024, we'd need multiple warps. But here dim=32. Let's write a general kernel: block_size = 256 (or next power of 2 >= dim). But we need a reduction across the block. Let's write block_size = 128 or 256, each thread processes a subset of dimensions or we loop with stride.
Simpler: each block has `dim` threads (or at least `dim` if dim is small). Since dim is known at runtime, we can set block size to dim. However, CUDA prefers block size at least 32 (warp) and multiples of warp. dim=32, so exactly one warp. That's fine.
But what if dim is not a multiple of warp? In the benchmark it's 32. But to be general, we could set block_size = 32 (or 64) and have each thread handle one dimension if threadIdx.x < dim. If dim < block_size, some threads are idle. For dim=32, block_size=32, all threads active.
However, kernel launch overhead is small. Let's write block size as `min(1024, ((dim + 31)/32)*32)`, i.e., round up to multiple of 32. Then in kernel, if (threadIdx.x >= dim) return; but for warp operations we need all threads in warp to participate. So if dim=32, block_size=32, perfect.
Let's write the kernel with block_size = 32 (or a template). Wait we can pass block_size as compile-time constant via constexpr or template, but since we want to handle variable dim, we can just set max block size to 1024 and pass dim as argument. In the kernel, threads with idx >= dim return (or do nothing). But warp shuffle requires all threads in mask to be active. If we have block_size = 64, but dim=32, threads 32..63 are inactive. If we call __shfl_down_sync with mask 0xffffffff, all threads in warp must be active or the result is undefined if some are exited. Actually `__shfl_down_sync` with a mask expects the calling threads to be in the mask. If some threads have returned, they are no longer active, but as long as they don't call the shfl, it's okay? In CUDA, if a thread has exited (returned), it won't execute the shfl. The mask should include only active threads. But if we use `__shfl_down_sync(0xffffffff, ...)` in a warp where some lanes have returned, those lanes won't call the instruction, but the others will. This is okay? Actually I think `__shfl_down_sync` requires that all threads in the mask execute it convergently. If some threads have diverged and returned, they are considered inactive, and the behavior may be undefined. To be safe, we should not early-return threads in a warp. Instead, we can process only if tid < dim, but still participate in the warp reduction by setting diff=0 for inactive threads. Or we can set block_size = dim (but block_size must be compile-time constant? No, we can set it dynamically in the launch config, and inside kernel use `if (threadIdx.x >= dim) { diff = 0; }` and still participate in shfl). Wait if block_size is dynamic, say `threads = dim;`, then the number of threads in the block is exactly dim. If dim < 32, then the warp is partially filled. But CUDA still runs 32 threads per warp; the extra threads are inactive (threadIdx.x >= blockDim.x). They are not supposed to execute code. Actually threads with threadIdx.x >= blockDim.x don't exist. The block size is the number of threads launched. The warp size is always 32, but if block_size < 32, the remaining lanes are inactive and masked off. `__shfl_down_sync` with mask 0xffffffff includes all 32 lanes. If only, say, 16 threads are launched, the other 16 are inactive and won't execute the instruction. According to CUDA docs, for `__shfl_*_sync`, all non-exited threads named in mask must execute the same `__shfl_*_sync` with the same mask. Inactive threads (not launched) are considered not named? Actually the mask includes all 32 lanes. The other lanes are not part of the thread block, so they are not active. The instruction is only executed by active threads. The hardware handles it by ignoring the inactive lanes. I think it's safe to use mask 0xffffffff for partial warps as long as all active threads call it. But to be extra safe, we can compute the mask as `(1 << blockDim.x) - 1` or something. But that's tricky if blockDim.x > 32. Actually for a full warp (32 threads), mask = 0xffffffff is fine. For partial warp, mask = (1u << blockDim.x) - 1. But blockDim.x can be up to 32. If dim=32, it's fine. If dim is say 64, we have two warps, and we need to reduce across 64 threads. We can do a block-level reduction using shared memory instead. That is more general.
Given that dim is 32 in the benchmark, we can target that but still be general. Let's write a kernel that uses shared memory for reduction, but for dim=32 it will be very fast. Or we can use warp shuffle reduction with block_size = dim (rounded up to nearest multiple of 32). For dim <= 32, block_size = 32 (or dim). Actually if dim=32, block_size=32. If dim=64, block_size=64, use two warps, reduce within warp then across warps with shared memory.
Let's write a generic kernel:
```
__global__ void svgd_kernel(const float* particles, const float* score, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
if (i >= N) return;
extern __shared__ float sdata[];
// sdata will be used for reduction of dist_sq across threads.
// We'll allocate blockSize floats.
float xi[...]? We can load x_i values into registers? Each thread can load its own dimension.
// Actually each thread handles one dimension d = tid, but if dim > blockSize, thread loops with stride blockSize.
// Let's handle arbitrary dim.
float accum = 0.0f;
// Preload x_i for dims handled by this thread
// We'll loop over j.
for (int j = 0; j < N; ++j) {
// Compute dist_sq between particle i and j
float dist_sq = 0.0f;
for (int d = tid; d < dim; d += blockSize) {
float diff = particles[j * dim + d] - particles[i * dim + d];
// We also need diff later for grad term. But if blockSize < dim, each thread handles multiple dims.
// Wait, if each thread handles multiple dims, we can't easily do a per-thread final term using diff because diff is different for each d.
// It's simpler: each thread handles exactly one d if dim <= blockSize.
}
}
}
```
If dim=32 and blockSize >= 32, each thread can handle one d. If blockSize=32, each thread has d = tid. That's the cleanest. Since dim is only 32, we can set blockSize = 32. For generality, we could handle dim > blockSize by looping, but the reduction of dist_sq becomes more complex because each thread computes multiple diff contributions. But we can still sum them in registers and then reduce across threads.
Simplify: Since the benchmark dim=32, we can hardcode block size = 32 (or maybe 64) and assume dim <= 64. But the ModelNew should accept any dim from constructor. If dim is larger, the kernel could still work with block size = 128 and each thread handling multiple dims in a strided loop. Let's see if we need to support arbitrary dim. The task says "Accept the same constructor arguments as Model." Model takes dim. The benchmark uses dim=32. Usually they only test with the benchmark inputs, but they might test with other dims. To be robust, let's support arbitrary dim.
Design for arbitrary dim:
- Block size: 128 (or maybe 256) threads.
- Each thread loops over dimensions with stride blockDim.x.
- For each j:
- each thread computes partial sum of squared diffs for dimensions it owns: `partial += diff*diff;`
- store the diff values for those dimensions in local registers (we might need them later for the term). But if each thread handles many dims, we can compute the term for each d as we go? Wait the term for each d is K * (score_j[d] + diff/h). K depends on the total dist_sq, which is not known until after reduction. So we need two passes over dimensions: first compute dist_sq (reduction), then compute term for each d using K. Since dimensions are many, we can do:
- First pass (or within the same j loop) compute partial dist_sq.
- Reduce partial dist_sq across block to get total dist_sq.
- Compute K.
- Second pass over dimensions: each thread computes term for its owned d and accumulates.
- This requires either recomputing diff in the second pass (cheap) or storing diffs in shared memory. Recomputing diff is fine: we have particles in global memory. But we need x_i and x_j again. x_i is constant for the block. x_j changes per j. So second pass reads x_j again. That's fine.
So kernel structure per block (i):
```
int i = blockIdx.x;
float h_local = h;
float inv_h = 1.0f / h_local;
float inv_2h = 0.5f * inv_h;
int N_local = N;
int dim_local = dim;
float accum[ Register array? No, if dim is large, we can't store per-d accum in registers for each thread if each thread handles many dims. But each thread can accumulate into a local variable for each d it handles? If it loops with stride, it can process one d at a time? No, because K is the same for all d in this j iteration. But each thread handles multiple d's across iterations. Actually we can have each thread process one d at a time in the inner loop after reduction. Let's think:
```
Alternative approach: each warp handles one j. Since dim=32, a warp can compute all dimensions for one j. That is very efficient. We can have each block contain multiple warps, each warp handles a subset of j's. For dim=32, this is ideal.
Generalizing: if dim is variable but small (like <= 128), we can set block size = dim * warpsPerBlock. But dim is not known at compile time. We can use dynamic block size in kernel launch. Inside the kernel, we can treat the block as having `warp_size` threads per warp, and each warp processes one j. But blockDim.x may not be a multiple of dim. Hmm.
Simplify for the benchmark: Since dim=32, we can set blockDim.x = 32 (one warp). Number of blocks = N (1024). This yields 1024 warps. GPU can handle many warps. Each warp processes all j sequentially. This is simple and efficient. For correctness with other dims, we could fall back to a PyTorch implementation if dim is large? But the evaluation likely only uses dim=32. But to be safe, we can support arbitrary dim by using a block size of 128, and each block handles one i. Inside, each thread loops over d with stride blockDim.x. For each j, we compute partial dist_sq, reduce via shared memory, compute K, then each thread loops over d again to compute terms and accumulate. At the end, divide by N and write to output. This works for any dim.
Let's design that kernel:
```
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
float accum = 0.0f; // this thread's accumulation for its assigned d's? No, each thread handles multiple d's, so accum alone isn't enough.
}
```
Wait, each thread handles multiple dimensions. For example, dim=32, blockSize=128. Then each thread handles d = tid, tid+128, ... but tid < 32? Actually if blockSize > dim, then each thread handles d = tid if tid < dim, else does nothing. So we can set blockSize to be the smallest multiple of 32 >= dim, e.g., blockSize = ((dim + 31) / 32) * 32. If dim=32, blockSize=32. If dim=64, blockSize=64. If dim=100, blockSize=128. This ensures blockSize >= dim, and each thread handles exactly one dimension (tid < dim). Then we don't need multiple dims per thread. The reduction across blockSize threads is needed. Since blockSize may be > 32, we need a block-level reduction. We can do it in shared memory.
So algorithm:
- blockSize = ((dim + warpSize - 1) / warpSize) * warpSize; (warpSize=32)
- Launch N blocks, each with blockSize threads.
- In kernel, if i >= N return.
- Allocate shared memory of size blockSize floats.
- For each j in 0..N-1:
- Each thread tid < dim:
- read x_i_d = particles[i*dim + tid];
- read x_j_d = particles[j*dim + tid];
- diff = x_j_d - x_i_d;
- sdata[tid] = diff * diff;
- threads with tid >= dim set sdata[tid] = 0.
- __syncthreads();
- Block reduction: for (int s = blockSize/2; s > 0; s >>= 1) { if (tid < s) sdata[tid] += sdata[tid + s]; __syncthreads(); }
- dist_sq = sdata[0];
- Compute K = expf(-dist_sq / (2.0f * h));
- __syncthreads(); // ensure sdata not reused while some threads still reading
- Each thread tid < dim recomputes diff (or we could store diff in a register). Since each thread only handles one d, it can keep diff in a register from before. But we used sdata for reduction, overwriting diff*diff. So we need to recompute diff or store it separately. We can store diff in a local register variable `diff_reg` before storing to sdata. After reduction, we reuse diff_reg.
- diff_reg = x_j_d - x_i_d;
- score_j_d = -x_j_d;
- term = K * (score_j_d + diff_reg / h);
- accum[tid] += term; // each thread accumulates its own dimension.
- End for j.
- After loop, each thread tid < dim writes out[i*dim + tid] = accum / N.
This is straightforward.
But for blockSize=32, dim=32, the shared memory reduction is over 32 elements. We could use warp shuffle instead to avoid syncthreads, but shared memory is also fast. However, for blockSize=32, a single warp, we can avoid __syncthreads after reduction if we use warp shuffle. But __syncthreads in a single warp is cheap? Actually __syncthreads() in a single warp can be omitted if using warp shuffle. But using shared memory with __syncthreads is fine. For blockSize up to 128, we need __syncthreads.
Potential performance issue: For each j, we do two __syncthreads() (one after loading, one after reduction? Actually we need one after loading into sdata, then reduction loop with syncthreads inside, then after reduction we read sdata[0], then we need to ensure all threads have read sdata[0] before any thread writes to sdata again in next j. But we can avoid writing to sdata in the second phase; each thread just computes term using its own diff_reg and K (which is in a register for all threads? We need to broadcast K to all threads). After reduction, sdata[0] contains dist_sq. All threads read it into a local variable. Then they compute K. Then they compute term using diff_reg in register. No need for __syncthreads() after reading sdata[0] if we don't write to shared memory again until the next j. Since all threads are in lockstep, the read from sdata[0] should be fine. However, to be safe and avoid race conditions (some threads might start writing sdata for next j while others are still reading sdata[0] from previous j), we need a __syncthreads() at the top of the next j loop before writing to sdata. Or we can put a __syncthreads() at the end of the j loop. So loop:
```
for (int j = 0; j < N; ++j) {
float diff_reg = 0;
if (tid < dim) {
float xi = particles[i*dim + tid];
float xj = particles[j*dim + tid];
diff_reg = xj - xi;
sdata[tid] = diff_reg * diff_reg;
} else {
sdata[tid] = 0;
}
__syncthreads();
// reduce
for (int s = blockSize/2; s > 0; s >>= 1) {
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq / (2.0f * h));
if (tid < dim) {
float term = K * (-particles[j*dim + tid] + diff_reg / h);
accum += term;
}
__syncthreads(); // ensure all threads are done before next j writes sdata
}
```
Wait, we need to read particles[j*dim + tid] again for the term if we want to avoid storing score. That's fine. But we also need diff_reg. Since we computed diff_reg before writing sdata, we can keep it in a register. However, if blockSize > dim, threads with tid >= dim skip the body. That's fine.
But note: In the reduction loop, `if (tid < s) sdata[tid] += sdata[tid + s]; __syncthreads();` is standard. For blockSize = 32, s goes 16,8,4,2,1. After each step, __syncthreads() ensures all writes are visible. At the end, sdata[0] is the sum.
We need to allocate shared memory dynamically or statically. Since blockSize is dynamic per launch, we can use dynamic shared memory: `extern __shared__ float sdata[];` and pass shared memory size in bytes when launching kernel: `svgd_kernel<<<N, blockSize, blockSize * sizeof(float)>>>(...)`.
But we also need to handle that blockSize is `((dim + 31) / 32) * 32`. In Python, we compute `threads = ((dim + 31) // 32) * 32`. If dim is 0? Not possible.
Wait, there is a subtlety: `blockSize` must be <= 1024. If dim > 1024, this fails. But dim is 32. We can still handle dim > 1024 by having each thread handle multiple dimensions. But not needed for benchmark. We can add a fallback to PyTorch for dim > 1024. But to be safe, let's handle arbitrary dim by setting threads = min(1024, ((dim + 31)//32)*32). If dim > 1024, then each thread handles multiple dims with stride threads. That complicates the reduction and accumulation.
Let's think: is it possible they test with dim larger than 1024? The task says dim is a constructor argument. In the provided benchmark it's 32. Usually these challenges keep to the benchmark config. But we might want to handle larger dim gracefully.
If dim > 1024, we can set threads = 1024. Each thread processes a subset of dimensions. For each j:
- Each thread loops over d from tid to dim-1 step blockSize, computes diff, accumulates partial dist_sq in a register `partial`, and stores diff in a local array? But we can't store diff for all dims because each thread handles many dims. However, we can compute partial dist_sq, then reduce `partial` across block to get total dist_sq. Then each thread loops again over its dims, computes term using K (known after reduction) and accumulates into an output buffer or a local accumulator per d? But each thread handles many dims, so it needs to accumulate each dimension separately. It could write directly to output[i*dim + d] += term, but we need to divide by N at the end. Since N is constant, we can accumulate in global memory with atomicAdd? That's slow. Better to have a per-block shared memory buffer of size dim, but dim might be large. Not feasible if dim > 1024*something.
Alternative: If dim is large, we can have each block handle a chunk of dimensions for a given i. For example, block (i, d_block) where d_block processes a tile of dimensions (say 32). Then each block handles one i and a subset of d. But then we need K for that i,j pair, which requires reducing across all dimensions. That would require inter-block communication. Not ideal.
Given the benchmark dim=32, I think it's acceptable to assume dim is small enough that a single block can handle all dimensions. Most such kernel tasks have small feature dims. Let's assume dim <= 1024. We can add an assertion or fallback. To be safe, we can add a fallback: if dim > 1024, compute using PyTorch. But the reference is PyTorch, so fallback is correct.
Actually, we can write a kernel that handles dim up to, say, 1024. For dim <= 1024, set threads = ((dim+31)//32)*32. For dim > 1024, we can still set threads=1024 and have each thread handle a single dimension in a tile of 1024 dimensions per block, and launch multiple blocks per i? Then we need to compute dist_sq across all dimensions. We could compute partial sums per block and store in a global buffer, then reduce. That's too complex. I think it's fine to support dim <= 1024. We can assert in ModelNew.forward if dim > 1024, raise NotImplementedError or fallback to PyTorch.
Actually, fallback to PyTorch is allowed: we can just call the original PyTorch implementation for dim > 1024. It's not reward hacking because it's not a try/except fallback to reference in the sense of routing around evaluation. It's a legitimate branch for unsupported sizes. But to be safe, we can just support all dims by using PyTorch for the whole thing if dim is large, but that shouldn't affect the benchmark.
Let's write the kernel to support dim <= 1024.
But wait, another thought: we can compute everything in one kernel without needing to precompute dist_sq and median, by doing a two-pass kernel: first pass compute all pairwise dist_sq and store in a global N*N buffer. Then we can use PyTorch median on that buffer. Then second pass compute phi. But we can also compute dist_sq in PyTorch. Let's compare overhead.
Option A: PyTorch matmul + median:
- particles @ particles.t(): O(N^2*dim) = 33M ops, cuBLAS, very fast.
- median on 1M elements: PyTorch uses a selection algorithm, maybe ~0.05 ms.
Total preprocessing: maybe <0.1 ms.
Option B: Custom kernel for dist_sq:
- Launch N blocks, each with blockSize threads, compute dist_sq[i,j] for all j. Each block writes N values. This is also O(N^2*dim). But less optimized than cuBLAS. Might be slower. But it avoids storing the large diff tensor. Still, matmul is highly optimized.
Let's check the cost of PyTorch matmul. For N=1024, dim=32, the GEMM is 1024x32 * 32x1024. That's a small GEMM. cuBLAS might use a simple kernel. It's very fast. So Option A is likely fine.
However, we need to ensure correctness. The difference between matmul trick and explicit diff sum of squares might cause slight differences. But we can test with the tool.
One more detail: In the reference, diff = particles.unsqueeze(0) - particles.unsqueeze(1) computes x_j - x_i, where the first dim is i (target) and second is j (source). So dist_sq[i,j] = sum_d (x_j[d] - x_i[d])^2. The matmul trick: dist_sq = norms_i + norms_j - 2 x_i·x_j = ||x_i||^2 + ||x_j||^2 - 2 x_i·x_j = ||x_j - x_i||^2. So exactly the same algebraically. Numerically, due to floating point, there can be cancellation. But for random data with typical magnitudes ~1-3, norms ~32, dot products ~0, dist_sq ~64. The subtraction is well-conditioned. So difference should be negligible.
But we also need median_sq. Reference does median on dist_sq.view(-1). The median of slightly different dist_sq values could change. However, the differences are tiny, and since h = median_sq / log(N), a tiny relative change in median_sq yields tiny relative change in h, and K = exp(-dist_sq/(2h)). The gradient term involves 1/h. The overall phi might be sensitive if h is very small, but median is large (~64). Let's test.
Alternative: compute dist_sq in custom kernel chunk by chunk to exactly match reference, but without matmul overhead? We can do this: preallocate a N*N buffer for dist_sq. Launch a kernel where each block handles an i, each thread handles a j. But N=1024, so N*N=1M threads. That's a lot but fine. Each thread computes sum_d (x_j[d] - x_i[d])^2 in a loop over d. This is straightforward and exact. Memory writes: 1M floats. The kernel would be:
```
__global__ void compute_dist_sq(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
float sum = 0;
for (int d = 0; d < dim; ++d) {
float diff = particles[i*dim + d] - particles[j*dim + d];
sum += diff * diff;
}
dist_sq[i*N + j] = sum;
}
```
This uses 1M threads if we launch grid (N, (N+255)/256), block 256. For N=1024, grid (1024,4), block 256 => 1024*4*256 = 1,048,576 threads. Each thread loops 32 times. Total instruction count ~32M iterations, each with a few FMAs. This is fine.
Then we call torch.median on the dist_sq tensor. Then we call the phi kernel.
This ensures exact match with reference dist_sq. The overhead is minimal. The phi kernel is separate.
But wait, the reference diff is x_j - x_i (particles unsqueeze(0) - unsqueeze(1)). So diff[i,j,d] = x_j[d] - x_i[d]. Our kernel above computes x_i - x_j. Wait: diff = x_j - x_i. So dist_sq = sum_d (x_j - x_i)^2, which is same as (x_i - x_j)^2. So sign doesn't matter for dist_sq. Good.
We can compute dist_sq using a custom kernel and store it. But PyTorch's matmul approach is likely faster. Let's decide: to be safe on correctness, we should test both with the tool. If matmul approach passes, great. If not, we can switch to the custom dist_sq kernel.
Actually, with a custom dist_sq kernel and a custom phi kernel, the total work is two passes over the data. But we can also fuse the dist_sq computation into the phi kernel and avoid storing dist_sq? No, we need median first.
Wait, we can compute median using the median of pairwise distances computed on the fly without storing all of them? There are selection algorithms that work with a stream, but not easily. Storing 4MB is fine.
Let's first try the PyTorch matmul+median approach plus custom phi kernel. If correctness fails, we'll add a custom dist_sq kernel.
Let's write the code.
We need to use `torch.utils.cpp_extension.load_inline`. The cpp_sources will contain the CUDA kernel and a host wrapper function (e.g., `forward_svgd`) that takes torch::Tensor particles, float h, and returns torch::Tensor phi.
The host wrapper can allocate output tensor, launch the kernel, and return.
We need to be careful about tensor device and type.
Host wrapper:
```
at::Tensor forward_svgd(at::Tensor particles, float h) {
int N = particles.size(0);
int dim = particles.size(1);
auto out = torch::empty_like(particles);
int threads = ((dim + 31) / 32) * 32;
// Ensure threads <= 1024
if (threads > 1024) threads = 1024;
dim3 blocks(N);
dim3 block_size(threads);
size_t shared_mem = threads * sizeof(float);
svgd_kernel<<<blocks, block_size, shared_mem>>>(
particles.data_ptr<float>(),
h,
N,
dim,
out.data_ptr<float>()
);
return out;
}
```
The kernel `svgd_kernel` as described.
Now, we need to consider score_j = -x_j. In reference, score is computed as -particles. Then K @ score + grad_K.sum(dim=0). We are computing phi_i = (1/N) sum_j K_{i,j} * (score_j + (x_j - x_i)/h). Since score_j = -x_j, term = K * (-x_j + (x_j - x_i)/h).
Simplify term: K * ( -x_j + (x_j - x_i)/h ) = K * ( x_j * (-1 + 1/h) - x_i / h ).
But we can just compute directly: diff = x_j - x_i; term = K * (-x_j + diff / h) = K * (-x_j + (x_j - x_i)/h). That's fine.
Let's write the kernel:
```
template <typename scalar_t>
__global__ void svgd_kernel(const scalar_t* particles, float h, int N, int dim, scalar_t* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
// local accumulator
float accum = 0.0f;
// loop over j
for (int j = 0; j < N; ++j) {
float diff_reg = 0.0f;
if (tid < dim) {
float xi = particles[i * dim + tid];
float xj = particles[j * dim + tid];
diff_reg = xj - xi;
sdata[tid] = diff_reg * diff_reg;
} else {
sdata[tid] = 0.0f;
}
__syncthreads();
// block reduction
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq / (2.0f * h));
if (tid < dim) {
float term = K * (-particles[j * dim + tid] + diff_reg / h);
accum += term;
}
__syncthreads();
}
if (tid < dim) {
out[i * dim + tid] = accum / (float)N;
}
}
```
Wait, there is a type mismatch: particles is `scalar_t*`, but we are reading into `float`. Since we only support fp32, we can just use `float`. But to be generic, we could template. But `load_inline` compiles for the types we use. Since the input is float32, we can just use `float`. Let's write the kernel for `float`.
However, we need to be careful about `sdata` type. We declared `extern __shared__ float sdata[];`. If particles are double, sdata should be double. But here it's float. So fine.
One more thing: `blockSize` must be a power of 2? The reduction `for (int s = blockSize / 2; s > 0; s >>= 1)` works for any blockSize, but if blockSize is not a power of 2, some elements might be ignored. Since we set blockSize = ((dim+31)/32)*32, which is a multiple of 32, not necessarily a power of 2. For example, dim=48 => blockSize=64 (power of 2). dim=65 => blockSize=96 (not power of 2). The reduction with s = blockSize/2 = 48, then 24, 12, 6, 3, 1. This works because 96 is even, 48 is even, etc., down to 3. At s=3, tid<3 reads sdata[tid+3]. This sums indices 0..2 with 3..5. But indices 6..95? Wait reduction usually assumes blockSize is a power of 2 or we pad. Let's trace:
We have 96 elements. After first step (s=48): tid<48 adds sdata[tid+48] to sdata[tid]. Now sdata[0..47] contain sums of pairs (0+48), (1+49), ..., (47+95). Good.
Next s=24: tid<24 adds sdata[tid+24] to sdata[tid]. Now sdata[0..23] contain sums of 4 elements each. Good.
Next s=12: sums of 8.
Next s=6: sums of 16.
Next s=3: tid<3 adds sdata[tid+3] to sdata[tid]. Now sdata[0..2] each contain sums of 32 elements? Let's see:
- sdata[0] = sum of indices 0,3,6,9,...? No.
Wait standard reduction algorithm with sequential addressing: at each step s, it pairs sdata[tid] and sdata[tid+s]. This works as long as all elements are included in the first step. Since we start with blockSize/2, and blockSize is the total number of threads, we are summing pairs (0, bS/2), (1, bS/2+1), etc. If blockSize is even, all elements are covered. Subsequent steps halve s. This works as long as blockSize is a power of two? Actually for any blockSize, as long as we start with s = blockSize/2 (integer division) and the initial data is in the first `blockSize` entries, some entries beyond the last active index might not be included if blockSize is not a power of two? Let's test blockSize=96.
Initial: sdata[0..95] all valid.
s=48: tid 0..47: sdata[tid] += sdata[tid+48]. Now sdata[0..47] valid.
s=24: tid 0..23: sdata[tid] += sdata[tid+24]. Now sdata[0..23] valid.
s=12: tid 0..11: sdata[tid] += sdata[tid+12]. Now sdata[0..11] valid.
s=6: tid 0..5: sdata[tid] += sdata[tid+6]. Now sdata[0..5] valid.
s=3: tid 0..2: sdata[tid] += sdata[tid+3]. Now sdata[0..2] valid. But what about sdata[3], sdata[4], sdata[5]? They are not added to 0..2.
s=1: tid 0..0: sdata[0] += sdata[1]. Now sdata[0] valid. But sdata[2] is not added.
So after s=1, sdata[0] contains sum of some subset, but not all. Specifically, after s=3, we have:
sdata[0] = sum of indices {0,3,6,9,12,15,...}? No, it's a tree.
Let's trace the set of indices contributing to sdata[0]:
Start: sdata[0] = val0
s=48: sdata[0] += val48
s=24: sdata[0] += sdata[24] = val24
s=12: sdata[0] += sdata[12] = val12
s=6: sdata[0] += sdata[6] = val6
s=3: sdata[0] += sdata[3] = val3
s=1: sdata[0] += sdata[1] = val1
So sdata[0] = val0 + val1 + val3 + val6 + val12 + val24 + val48. That's missing indices 2,4,5,7,... etc. So the reduction is wrong for non-power-of-two block sizes.
Thus we need blockSize to be a power of two, or use a different reduction method. We can round blockSize up to the next power of two. But if dim is 33, next power of two is 64. That's fine. Actually ((dim+31)/32)*32 is not a power of two unless dim <= 32. For dim=33..64, it's 64 (power of 2). For dim=65..96, it's 96 (not power of 2). But we can round up to the next power of two, e.g., blockSize = 1; while (blockSize < dim) blockSize <<= 1; If dim <= 1024, blockSize will be <= 1024. For dim=65, blockSize=128. That's fine. The number of threads launched is blockSize. Some threads (>=dim) will be idle but still need to set sdata[tid]=0. So the reduction will be correct.
So in host wrapper:
```
int threads = 1;
while (threads < dim) threads <<= 1;
if (threads > 1024) threads = 1024; // but then reduction may miss some dims if dim>1024
```
But if dim > 1024, our reduction doesn't work. So we need fallback. For now, compute threads = 1; while (threads < dim) threads <<= 1; if (threads > 1024) ... fallback. Let's just assume dim <= 1024.
Actually, for dim=32, threads = 32. Perfect.
Alternatively, we can avoid the power-of-two requirement by using a loop that sums all elements with active threads at each step but handling remainders. Simpler is to round up to power of two.
Another approach: use warp shuffle for the first 32 threads, and if blockSize > 32, do a second level reduction in shared memory. But power of two is easier.
So:
```
int threads = 1;
while (threads < dim) threads <<= 1;
```
Then launch with `threads`.
Now, `__syncthreads()` inside a loop for N=1024 times. Each iteration reads x_j from global memory. Since x_j is read by all blocks, it's not cached in L1 (unless N is small). But reading from global memory for each block and each j means total reads = N blocks * N * dim * 4 bytes = 132MB. At 1 TB/s, that's ~0.13 ms. Plus overhead. Should be very fast.
But there is an optimization: we can load particles into shared memory for a tile of j's, to reuse data across i blocks? No, each block handles a different i, so they all need the same x_j. But they don't share memory. L2 cache might help: each x_j is read by all N blocks. With N=1024, the total particle data is 1024*32*4 = 128KB. This fits in L2 cache easily. In fact, it fits in L1 (48KB) with some leftover. Actually 128KB > 48KB, but L2 is large (several MB). So as blocks iterate over j, they read x_j from L2 cache (since other blocks have loaded it). So memory bandwidth is not a bottleneck. Good.
But we have N blocks = 1024. If each block has 32 threads, total threads = 32768. That's 1024 warps. GPU can run many warps concurrently. But each block loops 1024 times over j. So total iterations = 1024 blocks * 1024 j = 1M iterations per block? No, each block does 1024 j's. With 32 threads per block, each j iteration does a small amount of work. The main cost is the loop overhead and syncthreads. 1024 iterations per block is a lot. But each iteration is short. However, 1024 __syncthreads() per block may be expensive. 1024 * 2 = 2048 __syncthreads() calls per block. For 1024 blocks, total 2M barrier ops. That might be okay but maybe not optimal.
We can improve by processing multiple j's per thread block in parallel. For example, have each warp handle one j. Since dim=32, one warp can process all dimensions for one j. Then a block of 256 threads (8 warps) can process 8 j's simultaneously. This reduces the number of __syncthreads() calls by factor of 8. But it increases register usage and requires partial accumulation per j per thread? Let's think.
In the current design, each block computes the entire sum for one i, iterating over j sequentially. The sequential loop has 1024 steps. Each step does a few operations and 2-3 syncthreads. On a GPU, a block with 32 threads is small; it uses only one warp. __syncthreads() for a single warp is effectively a no-op (or very cheap). Actually __syncthreads() inside a single warp can be omitted, but the compiler may still emit it. It might be okay. Let's keep it simple.
But we can further optimize by removing __syncthreads() from the reduction if we use warp shuffle. Since blockSize is 32 (one warp), we can do warp shuffle reduction. Then we don't need any shared memory or syncthreads().
Let's do that: if blockSize <= 32, use warp shuffle. But blockSize could be 64 if dim=48. That's two warps. Then we can still use warp shuffle within each warp, then combine in shared memory. But maybe it's easier to just use blockSize = 32 (warp) and handle only dim <= 32? But dim in benchmark is exactly 32. If they test with dim <= 32, this works.
Wait, the constructor argument is dim. ModelNew must accept same args. In forward, dim is particles.shape[1]. If dim > 32, our kernel won't handle it. We can add a fallback to PyTorch for dim > 32 or for any dim > 1024. But to be robust, we should handle arbitrary dim.
Let's write a version that handles any dim but optimizes for dim=32.
Option: For any dim, we can have each block handle one i, and each thread in the block handle one dimension, with blockSize = dim (or next multiple of 32). If dim is <= 32, blockSize=32, one warp. Use warp shuffle. If dim is <=64, blockSize=64, two warps. Use warp shuffle per warp, then store warp sums to shared memory, then one warp reduces.
Let's design a generic kernel that works for any blockSize (power of two) up to 1024. Actually we can set blockSize dynamically. In kernel, we can reduce using the classic tree reduction in shared memory with __syncthreads(). This is general and correct if blockSize is a power of two. We can ensure blockSize is the next power of two >= dim. For dim=32, blockSize=32. For dim=33, blockSize=64. For dim=100, blockSize=128. All fine.
The cost of __syncthreads() for blockSize=32 is minimal. For blockSize=64, also minimal. So we can keep the shared memory reduction.
But we need to make sure blockSize is a power of two. Let's enforce that in host code.
```
int threads = 1;
while (threads < dim) threads <<= 1;
if (threads > 1024) { /* fallback */ }
```
Now, there is a subtlety: __syncthreads() cannot be used in conditional code unless all threads in the block reach it. In our kernel, we have `if (tid < dim) { ... } else { sdata[tid] = 0; } __syncthreads();`. To avoid using `if (tid < dim)` before syncthreads, we can write:
```
float diff_reg = 0.0f;
float sq = 0.0f;
if (tid < dim) {
float xi = particles[i*dim + tid];
float xj = particles[j*dim + tid];
diff_reg = xj - xi;
sq = diff_reg * diff_reg;
}
sdata[tid] = sq;
__syncthreads();
```
This ensures all threads (including tid >= dim) write to sdata. Good.
Then reduction:
```
for (int s = blockSize/2; s > 0; s >>= 1) {
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads();
}
```
All threads participate? The condition `tid < s` means some threads (tid >= s) don't do anything, but they still reach __syncthreads(). In CUDA, all threads in the block must reach __syncthreads(). Since the __syncthreads() is outside the if, all threads reach it. Actually inside the for loop, the __syncthreads() is after the if. So all threads (tid >= s) skip the add but still hit __syncthreads(). This is valid because the branch is uniform? No, it's divergent, but __syncthreads() inside a branch is allowed only if all threads in the block take the same path or if the __syncthreads() is executed by all threads in the block regardless of path. In our code:
```
for (...) {
if (tid < s) { sdata[tid] += sdata[tid+s]; }
__syncthreads();
}
```
At each iteration, threads with tid >= s do not execute the add, but they all execute __syncthreads(). This is valid because __syncthreads() is reachable by all threads in the block (the if doesn't guard the syncthreads). Actually the CUDA Programming Guide says: `__syncthreads()` is allowed in conditional code only if the conditional evaluates identically across the entire thread block. Here the condition `tid < s` is not identical across threads, but __syncthreads() is outside the conditional (it's after the if statement). So all threads will execute __syncthreads(). This is allowed. Good.
After reduction, `float dist_sq = sdata[0];`. Since all threads need K, they can just read sdata[0]. There's no race because no thread writes to sdata after this point until next j iteration. But we have a __syncthreads() at the end of the j loop before the next j iteration writes to sdata again. So we are safe.
Thus the kernel is valid.
Now, one more optimization: we can precompute score = -particles outside the kernel, and pass it as an argument, to avoid the negation inside the loop. But it's trivial.
Also, we can precompute inv_h = 1.0f / h to avoid division inside loop.
And we can compute N_inv = 1.0f / N, multiply at the end.
Let's adjust the kernel:
```
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
float inv_N = 1.0f / (float)N;
```
Then inside loop:
```
float K = expf(-dist_sq * inv_2h);
if (tid < dim) {
float term = K * (-particles[j*dim + tid] + diff_reg * inv_h);
accum += term;
}
```
But wait: `-particles[j*dim + tid] + diff_reg * inv_h`. We can also precompute score pointer and read score_j. But reading particles is needed anyway for x_j. Negation is cheap.
Potential memory coalescing: reading `particles[j*dim + tid]` where tid is threadIdx.x. For a fixed j, threads read consecutive addresses. Coalesced. Good.
Now, we also need to handle h possibly being very small. The reference clamps median_sq to 1e-5, then h = median_sq / max(log(N),1). For N=1024, log(N)=~6.9. h >= 1e-5 / 6.9 ~ 1.4e-6. That's fine. No division by zero.
Now, let's think about whether we need to compute dist_sq with the same numeric result as reference. Using PyTorch for dist_sq might be okay, but let's be safe and compute it with a custom kernel to avoid any tolerance issues.
Actually, we can write a simple CUDA kernel for dist_sq and call torch.median on the result. That ensures exact match. The kernel:
```
__global__ void compute_dist_sq(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
float sum = 0.0f;
const float* xi = particles + i * dim;
const float* xj = particles + j * dim;
for (int d = 0; d < dim; ++d) {
float diff = xj[d] - xi[d];
sum += diff * diff;
}
dist_sq[i * N + j] = sum;
}
```
This kernel doesn't need shared memory. It uses 1 thread per (i,j) pair. For N=1024, we can launch grid (N, (N+255)/256), block 256. This is very simple and exact.
But wait: diff = xj - xi, same as reference. sum of squares. exact.
Then we can call `torch.median` on the `dist_sq` tensor. Then launch `svgd_kernel`.
Total time: two kernels + median. First kernel writes 4MB. Second kernel reads particles (128KB) many times but each block loops over j. The first kernel might be slightly slower than matmul because it does 32M scalar operations instead of highly optimized matmul. But it's still tiny. Let's estimate: 32 threads per block, 1024 blocks, each thread does 1024 loop iterations, each with a few FLOPs. Actually first kernel: 1M threads, 32 FLOPs each ~ 32M FLOPs, plus 4MB write. At GPU throughput, this is <0.01 ms. The PyTorch matmul is also <0.01 ms. So both are fine. The custom kernel for dist_sq is simpler and ensures exactness.
But wait, there is a subtlety: the reference computes dist_sq[i,j] = sum_d (particles[j,d] - particles[i,d])^2. In our first kernel, we compute diff = xj[d] - xi[d]. That's correct. Then we store in dist_sq[i*N + j]. So it's row i, col j. Good.
After computing median_sq and h, the phi kernel needs to compute phi_i. In phi kernel, we loop over j and read particles[j]. The order is sequential, which is good for cache.
But we can fuse everything into a single kernel if we compute median on the fly? No, median requires all dist_sq values. So two passes are needed.
Alternatively, we can avoid the first kernel entirely by using PyTorch's `torch.cdist` with p=2, then square it. `torch.cdist` likely computes distances using a block algorithm and might be very fast and numerically stable. But let's just use the custom kernel for dist_sq.
Actually, we can compute dist_sq using the same `svgd_kernel` by writing an extra buffer. But it's separate. Let's not overcomplicate.
Wait, another idea: can we compute the median without computing all N^2 distances? The median of pairwise distances for N points: Since the distribution of distances might be estimated via sampling. But reference uses exact median. So we need exact.
Given N=1024, the exact median of 1,048,576 numbers is required. We compute all of them.
Now, we need to produce a complete Python file defining `ModelNew`.
The file should:
- Import torch, torch.nn as nn, math.
- Define `ModelNew(nn.Module)`.
- In `__init__`, accept dim, call super().__init__(), self.dim = dim, register_buffer("target_mean", torch.zeros(dim)).
- In `forward(self, particles)`, compute N = particles.shape[0], device.
- Compute dist_sq via custom kernel OR PyTorch. Then median, h.
- Call custom kernel to compute phi.
- Return phi.
Let's write the host wrapper for both kernels.
But wait, we need to define the CUDA source code in `cpp_sources`. And we list wrapper functions in `functions=[...]`.
Also, we need to handle `score` as `-particles`? The reference registers a buffer `target_mean` but it's not used in the forward except for `_score_fn` which returns `-x`. Our ModelNew doesn't need to store target_mean for the forward if we fuse it, but we should keep it for compatibility (same constructor arguments and buffers). Actually the forward signature doesn't need `target_mean`. But the class should have the same attributes. Let's keep it.
Now, let's write the code.
```python
import torch
import torch.nn as nn
import math
from torch.utils.cpp_extension import load_inline
cuda_src = r'''
#include <torch/extension.h>
#include <cuda_runtime.h>
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
float sum = 0.0f;
const float* xi = particles + i * dim;
const float* xj = particles + j * dim;
for (int d = 0; d < dim; ++d) {
float diff = xj[d] - xi[d];
sum += diff * diff;
}
dist_sq[i * N + j] = sum;
}
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
float inv_N = 1.0f / (float)N;
float accum = 0.0f;
float diff_reg = 0.0f;
for (int j = 0; j < N; ++j) {
float sq = 0.0f;
if (tid < dim) {
float xi = particles[i * dim + tid];
float xj = particles[j * dim + tid];
diff_reg = xj - xi;
sq = diff_reg * diff_reg;
}
sdata[tid] = sq;
__syncthreads();
// tree reduction
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq * inv_2h);
if (tid < dim) {
float term = K * (-particles[j * dim + tid] + diff_reg * inv_h);
accum += term;
}
__syncthreads();
}
if (tid < dim) {
out[i * dim + tid] = accum * inv_N;
}
}
torch::Tensor forward_svgd(torch::Tensor particles) {
int N = particles.size(0);
int dim = particles.size(1);
auto options = torch::TensorOptions().dtype(torch::kFloat32).device(particles.device());
// Compute pairwise squared distances
auto dist_sq = torch::empty({N, N}, options);
dim3 block(256);
dim3 grid(N, (N + block.x - 1) / block.x);
compute_dist_sq_kernel<<<grid, block>>>(
particles.data_ptr<float>(), N, dim, dist_sq.data_ptr<float>());
// Median heuristic for bandwidth
float median_sq = dist_sq.view(-1).median().item<float>();
if (median_sq < 1e-5f) median_sq = 1e-5f;
float logN = logf((float)N);
float h = median_sq / fmaxf(logN, 1.0f);
// Compute phi
auto out = torch::empty_like(particles);
int threads = 1;
while (threads < dim) threads <<= 1;
// fallback if dim too large: use PyTorch
if (threads > 1024) {
auto score = -particles;
auto diff = particles.unsqueeze(0) - particles.unsqueeze(1);
auto K = torch::exp(-dist_sq / (2.0f * h));
auto grad_K = -diff / h * K.unsqueeze(-1);
auto phi = (K.matmul(score) + grad_K.sum(0)) / N;
return phi;
}
svgd_phi_kernel<<<N, threads, threads * sizeof(float)>>>(
particles.data_ptr<float>(), h, N, dim, out.data_ptr<float>());
return out;
}
'''
# Load inline
svgd_cuda = load_inline(
name="svgd_cuda",
cpp_sources="",
cuda_sources=cuda_src,
functions=["forward_svgd"],
extra_cuda_cflags=["-O3"],
verbose=False,
)
class ModelNew(nn.Module):
def __init__(self, dim):
super().__init__()
self.dim = dim
self.register_buffer("target_mean", torch.zeros(dim))
def forward(self, particles: torch.Tensor) -> torch.Tensor:
return svgd_cuda.forward_svgd(particles)
```
Wait, there are a few issues:
1. `cpp_sources` cannot contain a `PYBIND11_MODULE`. We are providing an empty string. The `forward_svgd` is defined in `cuda_sources`. Since it's a host function, it needs to be compiled as CUDA C++ (which allows host code). `load_inline` will handle the binding. This should be fine.
2. In `forward_svgd`, we compute `dist_sq` using a custom kernel, then call `dist_sq.view(-1).median()`. This is a PyTorch operation on a CUDA tensor. That's fine. It will launch its own kernel. The result is a scalar tensor. `.item<float>()` copies to CPU.
3. We then compute `h` on CPU.
4. We launch `svgd_phi_kernel` with `N` blocks and `threads` threads. Shared memory size is `threads * sizeof(float)`.
5. For the fallback when `threads > 1024`, we reuse the already computed `dist_sq` and `K` etc. But note: if we fallback, we need to recompute `diff` as `(particles.unsqueeze(0) - particles.unsqueeze(1))` which is huge. But fallback only for dim > 1024, which is unlikely. It's acceptable.
6. One issue: `median` on GPU for a tensor of size 1M might be somewhat slow compared to not computing it, but it is necessary. However, in the reference, `median` is also computed. So our time includes median. The reference also includes the huge diff and other ops. We should still be faster.
7. Potentially we can also compute `dist_sq` only for the upper triangle or something to save work? No, reference uses full matrix.
8. Another issue: In `compute_dist_sq_kernel`, `blockDim.x` is 256. So `j = threadIdx.x + blockIdx.y * blockDim.x`. `grid` is `(N, (N+255)/256)`. For N=1024, grid.y = 4. So total blocks = 1024*4 = 4096. Each block has 256 threads. Total threads = 1,048,576. Good.
9. However, `particles` is read per thread. With 1M threads each reading 32 floats, that's 32M reads. But each read is from global memory. The L2 cache will help because many threads read the same xi and xj. For a given i, all threads in blocks with blockIdx.x = i read the same xi. Since blocks with same blockIdx.x are across grid.y (j chunks), they are launched sequentially, but xi may stay in L2. xj are read by many threads. Overall bandwidth is fine.
10. But we can improve `compute_dist_sq_kernel` by loading xi into shared memory. Since blockDim.x=256 and dim=32, we can have the first 32 threads load xi into shared memory, then all 256 threads compute distances to different j's using that shared xi. This reduces reloading of xi. Let's rewrite the dist_sq kernel for better performance:
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
extern __shared__ float s_xi[];
if (threadIdx.x < dim) {
s_xi[threadIdx.x] = particles[i * dim + threadIdx.x];
}
__syncthreads();
float sum = 0.0f;
const float* xj = particles + j * dim;
for (int d = 0; d < dim; ++d) {
float diff = xj[d] - s_xi[d];
sum += diff * diff;
}
dist_sq[i * N + j] = sum;
}
```
Shared mem size = dim * sizeof(float) = 128 bytes. Very small. This kernel is more efficient.
But for simplicity, we can also skip the custom dist_sq kernel and use PyTorch's `torch.cdist` or matmul. Let's test correctness with the custom dist_sq kernel first.
One important concern: The reference computes `dist_sq` as `(diff**2).sum(-1)` where diff is `particles.unsqueeze(0) - particles.unsqueeze(1)`. This yields a (N,N) tensor. Then `median_sq = torch.median(dist_sq.view(-1)).clamp(min=1e-5)`. Then `h = median_sq / max(torch.log(torch.tensor(float(N))).item(), 1.0)`. Note: `torch.log(torch.tensor(float(N))).item()` computes natural logarithm on CPU. For N=1024, log(1024) = 6.9077.
In our code, we compute `logf((float)N)` and `fmaxf(logN, 1.0f)`. That's the same.
Now, `clamp(min=1e-5)` is done after median. So `median_sq = max(median, 1e-5)`. We do `if (median_sq < 1e-5f) median_sq = 1e-5f;`. Same.
One more nuance: The reference uses `torch.median(dist_sq.view(-1))` which returns the lower median for even number of elements (since there are N*N elements, which is even for N>1). PyTorch's median on a flattened tensor returns the lower median (the element at position (n//2)). Actually for `torch.median` with keepdim=False, it returns the median. For even number of elements, it returns the lower of the two middle values (the one at index n//2 - 1? Let's check). In PyTorch, `torch.median` for an even number of elements returns the smaller of the two middle values? I think it returns the element at index `n//2` after sorting? Actually no: the `kthvalue` at `n//2` (0-indexed) for even n is the upper middle? Wait, n=4, sorted [a0,a1,a2,a3]. median should be a1 or a2. PyTorch returns a1 (the smaller one) for `torch.median`? Let's verify: In PyTorch, `torch.median(torch.tensor([1,2,3,4]))` returns 2. That's index n//2 - 1? n//2 = 2, index 2 is 3. Wait 1,2,3,4: median is (2+3)/2 = 2.5, but torch.median returns 2. So it uses the lower median (floor((n-1)/2) = 1). For n=1,048,576, the returned index would be (n-1)//2 = 524,287 (0-indexed). Actually `torch.median` uses `kthvalue` with k = (n + 1) // 2? Let's check source: PyTorch's median calls `kthvalue` with k = (n + 1) // 2 for 1-indexed? If n=4, (4+1)//2 = 2 (1-indexed), which is element 2 (value 2). So yes, it uses the lower median for even n. So our use of `torch.median` on the `dist_sq` tensor will match exactly. Good.
Now, let's consider whether we should compute `dist_sq` in float32. Yes.
One more nuance: In the reference, after `K = torch.exp(-dist_sq / (2.0 * h))`, they compute `grad_K = -diff / h * K.unsqueeze(-1)`. Note `diff` is (N,N,dim). `K.unsqueeze(-1)` is (N,N,1). Then `phi = (K @ score + grad_K.sum(dim=0)) / N`.
As we derived, phi_i = (1/N) * sum_j K_{i,j} * (score_j + (x_j - x_i)/h).
But is `grad_K.sum(dim=0)` really sum over i? Let's check shapes: grad_K shape (N,N,dim). `sum(dim=0)` sums over first dimension (size N), resulting in shape (N,dim). The first dimension of grad_K corresponds to the first dimension of diff, which is i (target). Wait diff = particles.unsqueeze(0) - particles.unsqueeze(1). The first dimension (dim 0) corresponds to the i index of the subtraction result. diff[i,j] = x_j - x_i. So grad_K[i,j] = -(x_j - x_i)/h * K[i,j]. Then grad_K.sum(dim=0) gives sum_i grad_K[i,j] for each j. So it's sum over target i for each source j. But phi is (K @ score + grad_K.sum(dim=0)) / N. K @ score gives (N,dim) where row j = sum_i K[j,i] score_i? Wait K shape is (N,N). K[i,j] = exp(-dist_sq[i,j]/(2h)). Matrix multiplication K @ score: score shape (N,dim). (K @ score)_{j,d} = sum_i K_{j,i} * score_{i,d}. So row j of K@score corresponds to sum over i of K_{j,i} * score_i. Meanwhile grad_K.sum(dim=0) gives a vector indexed by j: sum_i grad_K[i,j]. Wait grad_K[i,j] is target i, source j. So sum_i grad_K[i,j] is summing over target i with source j fixed. That's (N,dim) where row j = sum_i grad_K[i,j].
But phi should be indexed by target? Let's check the shapes: (K @ score) is (N,dim). grad_K.sum(dim=0) is (N,dim). They add them and get (N,dim). The output is (N,dim). So the output index corresponds to j (source) or i (target)? In the reference, the output is just (N,dim). Since both K@score and grad_K.sum(dim=0) are indexed by the second dimension (j), the output corresponds to j. But the docstring says phi(x_i) = ... . So maybe they treat the output index as the particle index, and because the sum in grad_K.sum(dim=0) is over the first dimension (i), the remaining index is j, which is the particle index. But in the formula, the sum over j is over source particles. However, due to symmetry, sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_j, x_i) ??? Let's compute:
grad_K[i,j] = ∇_{x_j} K(x_j, x_i)? Actually K(x_j, x_i) is a function of x_j and x_i. The gradient w.r.t the first argument x_j is -(x_j - x_i)/h * K. In grad_K[i,j], the first index i corresponds to x_i (target), second index j corresponds to x_j (source)? Wait diff[i,j] = x_j - x_i. The first index i is the row (target), second j is column (source). grad_K[i,j] = -(x_j - x_i)/h * K[i,j] = ∇_{x_j} K(x_j, x_i). So grad_K[i,j] is the gradient w.r.t x_j when target is x_i. Then grad_K.sum(dim=0)[j] = sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_j, x_i). This is a sum over all target particles i of the gradient of K w.r.t the source particle j. In SVGD, the update for particle j involves summing over all other particles i: K(x_i, x_j) score(x_i) + ∇_{x_j} K(x_i, x_j). Wait the standard formula is:
phi(x_j) = (1/N) sum_i [ K(x_i, x_j) ∇log p(x_i) + ∇_{x_j} K(x_i, x_j) ].
Because x_j is the point at which we evaluate the update (the output particle), and x_i are the source particles in the sum. In this formula, the kernel is K(x_i, x_j), and the gradient is w.r.t x_j (the second argument). The derivative of K(x_i, x_j) w.r.t x_j is -(x_j - x_i)/h * K = (x_i - x_j)/h * K? Wait earlier we computed derivative w.r.t second argument x_j: K = exp(-||x_i - x_j||^2/(2h)). ∂K/∂x_j = -(x_j - x_i)/h * K. Yes, that's exactly grad_K[i,j] where i corresponds to x_i and j to x_j. And K[i,j] = K(x_j, x_i)? Wait K[i,j] is defined as exp(-dist_sq[i,j]/(2h)) where dist_sq[i,j] = ||x_j - x_i||^2. So K[i,j] = K(x_j, x_i) = K(x_i, x_j). So K@score: row j = sum_i K[j,i] score_i = sum_i K(x_i, x_j) score(x_i). Wait K[j,i] corresponds to dist_sq[j,i] = ||x_i - x_j||^2. So K[j,i] = K(x_i, x_j). So row j of K@score is sum_i K(x_i, x_j) score_i. That's exactly the first term for phi(x_j). Good.
Then grad_K.sum(dim=0)[j] = sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_i, x_j). That's the second term for phi(x_j). So phi_j = (K@score)_j + grad_K.sum(dim=0)_j all over N. So the output index is j, the source particle. But the docstring says "phi(x_i)". That's just notation; the output corresponds to each particle. So our phi kernel should compute for each i (output index) the sum over j of K[j,i] * score_j + grad_K[j,i]? Wait careful:
If output index is j, then phi_j = sum_i K(x_i, x_j) score(x_i) + sum_i ∇_{x_j} K(x_i, x_j).
In our variables: for output j, sum over i.
K_{j,i} = K(x_i, x_j). Wait K[j,i] = exp(-||x_i - x_j||^2/(2h)). Yes.
score_i = -x_i.
grad_K[i,j] = ∇_{x_j} K(x_i, x_j) = -(x_j - x_i)/h * K[i,j]. Note grad_K[i,j] uses first index i (source x_i) and second index j (target x_j). So grad_K[i,j] is derivative w.r.t second argument x_j when source is x_i. Then sum_i grad_K[i,j] = sum_i ∇_{x_j} K(x_i, x_j). That's grad_K.sum(dim=0)[j] because dim=0 is the first index i. So yes, output index j.
In our phi kernel, if we let blockIdx.x = j (output index), and loop over i (source index), then:
- diff = x_i - x_j? Let's see: K_{j,i} = K(x_i, x_j). We need diff for grad_K[i,j] = -(x_j - x_i)/h * K_{i,j}. But K_{i,j} = K_{j,i}. And -(x_j - x_i) = x_i - x_j. So grad_K[i,j] = (x_i - x_j)/h * K_{j,i}. Wait reference uses diff = particles.unsqueeze(0) - particles.unsqueeze(1). diff_{i,j} = x_j - x_i. grad_K_{i,j} = -(x_j - x_i)/h * K_{i,j} = (x_i - x_j)/h * K_{i,j}. For output j, we sum over i: grad_K_{i,j} = (x_i - x_j)/h * K_{i,j} = (x_i - x_j)/h * K_{j,i}. So the term inside sum for output j is:
K_{j,i} * score_i + (x_i - x_j)/h * K_{j,i}.
= K_{j,i} * (score_i + (x_i - x_j)/h).
In our earlier derivation, we used output i and sum over j: term = K_{i,j} * (score_j + (x_j - x_i)/h). Since indices are dummy variables, this is the same formula. So it's fine.
Thus, the kernel with blockIdx.x = i and loop over j computes:
diff = x_j - x_i.
term = K * (-x_j + diff / h) = K * (-x_j + (x_j - x_i)/h) = K * (score_j + (x_j - x_i)/h).
This matches the formula for output i.
So the kernel is correct.
But wait: if we compute dist_sq[i,j] = ||x_j - x_i||^2 in the first kernel, and then in the phi kernel we use that same h, and compute K = exp(-dist_sq/(2h)), then phi_i = sum_j K_{i,j} * (score_j + (x_j - x_i)/h) / N. That's correct.
Now, we must be careful about numerical precision: we accumulate `accum` in a float register over N=1024 terms. Each term is K * something. K is between 0 and 1. `something` is O(1/h) where h ~ 10. So term is O(0.1). Sum of 1024 terms is O(100). Float32 has ~7 decimal digits, so error is tiny.
One more thing: we should use `-O3` in CUDA flags. We can also use `--use_fast_math`? No, `-O3` is enough. `--use_fast_math` may change precision of `expf`. The tolerance is 1e-4. `--use_fast_math` might cause larger errors. Let's avoid it.
Now, about the `forward_svgd` wrapper. It takes `particles` as input, computes `dist_sq` using the custom kernel, then median, then phi. But we also need to think about whether the `dist_sq` allocation is necessary. Could we compute median without allocating a full N*N tensor? Actually, we could compute dist_sq in chunks and feed to `torch.median` incrementally? No. But 4MB is fine.
Wait, but the reference implementation does not allocate dist_sq explicitly? It does: `dist_sq = (diff ** 2).sum(-1)` creates a (N,N) tensor. So both allocate similar memory for dist_sq. The reference also allocates diff (N,N,dim) which is huge. So we avoid that.
Could we further avoid allocating the 4MB dist_sq by computing it, storing it, then median, then reusing the same memory for something else? Not needed.
One potential issue: `dist_sq` is allocated with `torch::empty`. The custom kernel writes all values. If for some reason some values are not written (e.g., when N is not a multiple of block size), they will contain garbage. We handle that by checking bounds in the kernel. For N=1024 and block=256, grid.y = 4, so 4*256 = 1024, exactly N. So all values are written.
Another issue: the kernel `svgd_phi_kernel` uses `sdata` with dynamic shared memory size `threads * sizeof(float)`. For threads up to 1024, that's 4KB. That's well within limits.
Edge case: if dim=0, threads=1. Not possible.
Now, let's consider the fallback for dim > 1024. The fallback uses the already computed `dist_sq` (via custom kernel) to compute K, then PyTorch. But wait, if dim > 1024, the custom `compute_dist_sq_kernel` still works fine because it doesn't depend on block size relative to dim. It uses loops over dim. So we can always compute dist_sq with the custom kernel. Then for phi, if dim > 1024, we fallback to PyTorch. In the fallback, we need `diff = particles.unsqueeze(0) - particles.unsqueeze(1)` which is huge for large dim but we only hit this fallback for dim > 1024. For dim > 1024, N might be smaller? We don't know. But fallback is fine.
Actually, for large dim, we might want to avoid the huge diff allocation too. But fallback is just for safety.
Let's refine the fallback. For dim > 1024, we can compute phi using the PyTorch ops but using the matmul trick to avoid diff:
```
auto score = -particles;
auto K = torch::exp(-dist_sq / (2.0f * h));
// We need grad_K sum:
// phi_i = (K @ score + sum_j (x_j - x_i)/h * K_{i,j}) / N
// = (K @ score + (K.unsqueeze(-1) * (particles.unsqueeze(0) - particles.unsqueeze(1)) / h).sum(1)) / N
// But this creates diff (N,N,dim). For dim large, it's huge.
```
Alternatively, we can compute the repulsion term as `(particles.t() @ K).t() / h - particles * (K.sum(1, keepdim=true) / h)`? Let's derive:
sum_j (x_j - x_i)/h * K_{i,j} = (1/h) * (sum_j x_j * K_{i,j} - x_i * sum_j K_{i,j}).
= (1/h) * ((K @ particles)_i - x_i * (K.sum(dim=1))_i).
This avoids the huge (N,N,dim) diff! Great. So for the fallback, we can compute:
```
auto score = -particles;
auto K = torch::exp(-dist_sq / (2.0f * h));
auto attr = K.matmul(score); // (N,dim)
auto K_sum = K.sum(dim=1, keepdim=true); // (N,1)
auto repulsion = (K.matmul(particles) - particles * K_sum) / h; // (N,dim)
auto phi = (attr + repulsion) / N;
```
Wait, is this exactly equal to the reference? Let's verify.
Reference: phi = (K @ score + grad_K.sum(dim=0)) / N.
grad_K = -diff / h * K.unsqueeze(-1), diff = x_j - x_i.
grad_K[i,j] = -(x_j - x_i)/h * K[i,j].
grad_K.sum(dim=0)[j] = sum_i grad_K[i,j] = sum_i -(x_j - x_i)/h * K[i,j] = (1/h) sum_i (x_i - x_j) * K[i,j].
= (1/h) (sum_i x_i K[i,j] - x_j sum_i K[i,j]).
= (1/h) ((K.t() @ particles)_j - x_j (K.sum(dim=0))_j). Wait K[i,j] is row i, col j. sum_i K[i,j] = K.sum(dim=0)[j] = column sum.
sum_i x_i K[i,j] = (K.t() @ particles)[j].
But earlier we said output index is j. So repulsion for j = (K.t() @ particles)[j] - x_j * col_sum[j].
Now K is symmetric (K[i,j] = K[j,i]). So K.t() = K. And col_sum = row_sum = K.sum(dim=1) (since symmetric). But wait, sum_i K[i,j] = sum over rows for column j. Because K is symmetric, this equals sum over columns for row j: K.sum(dim=1)[j]. So we can write:
repulsion_j = (K.matmul(particles)[j] - x_j * K.sum(dim=1)[j]) / h.
Because K @ particles = (N,N) @ (N,dim) = (N,dim), where row j = sum_i K[j,i] x_i. Since K symmetric, K[j,i] = K[i,j]. So row j = sum_i K[i,j] x_i = (K.t() @ particles)[j]. Good.
So repulsion = (K.matmul(particles) - particles * K.sum(dim=1, keepdim=true)) / h.
But wait, the reference uses grad_K.sum(dim=0). Let's compute for the fallback if we use `K @ score + (K @ particles - particles * K.sum(dim=1, keepdim=true)) / h`. Dividing by N.
But is this the same as the reference even if K is not exactly symmetric due to numerical errors? K should be symmetric because dist_sq is symmetric. The custom dist_sq kernel computes dist_sq[i,j] = sum_d (x_j - x_i)^2, which is symmetric. So K is exactly symmetric.
Now, what if we use the PyTorch matmul for dist_sq? It might not be exactly symmetric due to floating point. Then K might not be exactly symmetric. But we are using the custom dist_sq kernel, so symmetry holds.
For the fallback, we can use this efficient PyTorch formulation. But we only need fallback for dim > 1024. Since dim=32, we won't hit it.
Now, can we use this efficient formulation for the main path too, avoiding the custom phi kernel? Let's see: phi = (K @ score + repulsion) / N. Where K is (N,N), score is (N,dim), repulsion = (K @ particles - particles * K.sum(dim=1, keepdim=true)) / h. This involves:
- K = exp(-dist_sq/(2h)): element-wise on (N,N) -> 1M elements.
- K @ score: (N,N) @ (N,dim) -> O(N^2*dim) = 33M ops.
- K @ particles: (N,N) @ (N,dim) -> another 33M ops.
- K.sum(dim=1): O(N^2) = 1M ops.
- particles * K.sum(...): O(N*dim) = 32k ops.
- subtraction and division: O(N*dim).
Total compute: about 66M FLOPs for matmuls, plus memory for K (4MB), dist_sq (4MB). This is entirely within PyTorch, no custom kernel needed for phi.
Would this be faster than the reference? Let's compare:
Reference:
- diff = particles.unsqueeze(0) - particles.unsqueeze(1): allocates (N,N,dim) = 512MB. O(N^2*dim) compute.
- dist_sq = (diff**2).sum(-1): O(N^2*dim) compute.
- median: O(N^2).
- K = exp(-dist_sq/(2h)): O(N^2).
- grad_K = -diff/h * K.unsqueeze(-1): O(N^2*dim) compute, allocates another 512MB.
- K @ score: O(N^2*dim).
- grad_K.sum(dim=0): O(N^2*dim).
- phi = ... / N: O(N*dim).
Reference allocates two huge tensors (diff and grad_K), each 512MB. That's a lot of memory traffic. The efficient PyTorch formulation avoids those:
- Compute dist_sq via custom kernel (or matmul): but we still need dist_sq.
- K: 4MB.
- K @ score and K @ particles: two matmuls.
But we can avoid the huge allocations entirely. The custom phi kernel is even more efficient because it doesn't materialize K (4MB) either; it computes K on the fly. It also doesn't need the extra matmul; it fuses the repulsion into the same loop. For N=1024, the custom phi kernel will likely be faster than two matmuls plus elementwise ops.
Let's think about the overhead. Custom phi kernel: 1024 blocks, 32 threads each. Each block loops 1024 times. Each iteration does:
- 2 global reads (xi, xj).
- 1 local diff.
- 1 sq into shared mem.
- 5 steps of reduction (s=16,8,4,2,1).
- 1 exp.
- 1 global read (xj again).
- 1 FMA.
- accum.
Total per iteration ~ maybe 20 instructions. For 1024*1024 = 1M iterations, total 20M instructions. At GPU clock ~1.5GHz and ~100 instructions per cycle per SM, this is tiny, maybe 0.01 ms.
The two matmuls (K @ score and K @ particles) are 33M multiply-adds each, total 66M. cuBLAS is highly optimized and may be similar speed. But K @ score requires reading K (4MB) and score (128KB) and writing (128KB). K @ particles similarly. The custom kernel reads particles (132MB total) and writes 32KB. So memory traffic is comparable but custom kernel avoids temporary K. Actually K is computed on the fly so no read of K.
But the custom kernel also requires the dist_sq kernel to compute median. That writes 4MB. So total memory: write 4MB (dist_sq) + read particles 132MB + write 32KB = 136MB. PyTorch efficient version: compute dist_sq (writes 4MB if custom, or more if matmul intermediates), write K (4MB), read K twice for matmuls (8MB), read particles twice (256KB), write phi (32KB). Total ~16MB. So custom kernel approach reads particles more (132MB vs 256KB + K reads). Wait, the PyTorch matmul K @ particles reads K (4MB) and particles (128KB). K @ score reads K (4MB) and score (128KB). So total reads from K are 8MB, from particles 256KB. That's 8.25MB. Plus writing phi 32KB. The dist_sq computation in PyTorch (matmul) reads particles twice for the matmul and norms, maybe ~256KB, and writes dist_sq 4MB. So total memory traffic is much lower for PyTorch matmul approach! Because the matmul reads K (a dense matrix) but particles only once. The custom phi kernel reads particles N times per block (N^2 * dim). That's 132MB. So the PyTorch approach with two matmuls might actually be faster in terms of memory bandwidth, because it leverages the highly optimized GEMM which can reuse data in shared memory.
But wait, the reference also does K @ score and grad_K sum. The reference's slowness is due to the 512MB diff and grad_K. The efficient PyTorch formulation (without custom phi kernel) avoids that. Let's evaluate if that alone is enough to beat the reference by a good margin. Likely yes. But the custom phi kernel might be even faster, or it might be slower due to reading particles N times. However, 132MB is still tiny for modern GPUs. The overhead of kernel launch and many small operations (block reduction loop) might be okay.
Actually, the custom phi kernel has 1024 blocks, each with 32 threads, and each block loops 1024 times. That's 1024 blocks * 1024 iterations = 1,048,576 iterations of work distributed over 32 threads per block. The GPU has many SMs, say 80 SMs. Each SM can run multiple blocks concurrently. The occupancy might be limited by the number of warps per SM. Each block has 1 warp. Each SM can run up to 32 warps (e.g., A100). So each SM can run 32 blocks. Total blocks 1024, so 1024/80 = 12.8 blocks per SM. If each SM runs 32 blocks concurrently, we have enough. But each block has a loop of 1024 iterations. The latency of the loop might be high due to __syncthreads(). However, since there are many blocks, the GPU can switch to other warps while one is waiting at syncthreads.
But the PyTorch approach with two matmuls is probably simpler and less error-prone, and might be extremely fast because cuBLAS is highly optimized. Let's consider writing a solution that uses PyTorch ops for everything except avoids the huge intermediate.
Plan B (Pure PyTorch, no custom phi kernel):
1. Compute norms = (particles ** 2).sum(dim=1, keepdim=True) # (N,1)
2. Compute dist_sq = norms + norms.t() - 2.0 * particles @ particles.t() # (N,N)
- Clamp to 0 to avoid negative.
3. median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
4. h = median_sq.item() / max(math.log(float(N)), 1.0)
5. K = torch.exp(-dist_sq / (2.0 * h))
6. score = -particles
7. phi = (K.matmul(score) + (K.matmul(particles) - particles * K.sum(dim=1, keepdim=True)) / h) / N
This uses:
- One matmul for dist_sq (actually `particles @ particles.t()`).
- Two matmuls for phi.
- No huge intermediate tensors.
Would this be faster than the custom kernel approach? It depends. The matmuls are large enough to be efficient. The custom kernel approach is also efficient. Let's think about where the bottleneck is.
In the evaluation, the reference is slow because of the huge (N,N,dim) diff and grad_K. Both our approaches avoid that. But the pure PyTorch approach still allocates K (4MB) and dist_sq (4MB). It does three matmuls: particles@particles.t() (1024x32 * 32x1024), K@score (1024x1024 * 1024x32), K@particles (1024x1024 * 1024x32). The first matmul is small. The second and third are moderate. On a modern GPU, each 1024x1024x32 GEMM takes maybe 0.02-0.05 ms. Total maybe 0.1-0.2 ms. The median on 1M elements might take 0.05-0.1 ms. Total maybe 0.3 ms.
The custom kernel approach: dist_sq kernel (1M threads, 32 flops each) maybe 0.02 ms. Median 0.05 ms. Phi kernel: 1024 blocks, 32 threads, loop 1024. Each iteration does a few ops and syncthreads. Might be ~0.05-0.1 ms. Total maybe 0.15 ms. So custom kernel might be slightly faster.
But writing the custom kernel is more complex and riskier for correctness. However, the challenge asks for a custom CUDA kernel. The instructions say "Replace the forward computation of the PyTorch model below with a custom CUDA kernel." So we definitely need to write a custom kernel, not just use PyTorch ops. The evaluation checks that we use a custom kernel? Not explicitly, but the task says to write a custom CUDA kernel. We should do that.
But we can combine both: compute dist_sq with a custom kernel (to ensure correctness), then compute phi with a custom kernel. This fits the requirement.
One more thought: could we fuse median computation into the kernel? If we compute dist_sq in the phi kernel's first pass, we could store them in shared memory per block, but the median needs global values. So no.
What about computing the median on CPU by copying dist_sq? 4MB copy is fast.
Another idea: we can avoid the separate dist_sq kernel by computing dist_sq in the phi kernel using an atomicAdd to compute the sum of distances? No, we need the median, not the mean.
Could we estimate h using the mean of dist_sq instead of median? That would be a different algorithm, not matching reference.
So two kernels is fine.
Now, let's think about block size for the phi kernel. For dim=32, threads=32. Shared memory = 32*4 = 128 bytes. The reduction loop: s=16,8,4,2,1. 5 steps. Each step has syncthreads. In a single warp, syncthreads is not needed. Actually for a single warp, we can use warp shuffle reduction. That would be faster. Let's write a specialized path for blockSize <= 32. Or even for any blockSize, we can do warp shuffle for the first 32 threads, then combine warps.
But for simplicity, the shared memory reduction with syncthreads is fine. Since N=1024, the kernel will have 1024 blocks, each with 32 threads. Each block iterates 1024 times. Each iteration has 5 syncthreads. 5*1024 = 5120 syncthreads per block. For 1024 blocks, total 5.2M syncthreads. In a single warp, __syncthreads() translates to a warp-level barrier? Actually __syncthreads() in a block of 32 threads (one warp) is a very lightweight operation, essentially a no-op or a single instruction? According to CUDA docs, __syncthreads() on a single warp is effectively free because all threads in a warp are already synchronized. The compiler might even optimize it away. But with dynamic shared memory and writes, it still ensures memory ordering within the warp. Since a warp executes in lockstep (except for memory divergence), it might be fine. But to be safe, we could use `__syncwarp()` instead. But `__syncthreads()` works.
Actually, for a block of 32 threads (one warp), `__syncthreads()` does not need to synchronize across warps (since there's only one). But it still acts as a memory fence for shared memory. The overhead is minimal.
One more optimization: we can unroll the reduction loop for blockSize=32 using `#pragma unroll`. But not necessary.
Let's test both approaches (custom phi kernel vs pure PyTorch) for correctness and speed. Since we are limited to 30 tool calls, we should iterate.
I will write the kernel with custom dist_sq and custom phi. Let's draft the code.
Wait, there is a potential bug: In `compute_dist_sq_kernel`, we allocate shared memory for s_xi. We need to declare it as `extern __shared__ float s_xi[];` and pass the size in launch. For dim=32, it's 128 bytes.
But note: the block has 256 threads. Only first 32 threads load s_xi. Then all 256 threads read it. That's fine.
Alternatively, we can skip shared memory for s_xi and just read xi directly from global memory. Since xi is the same for all threads in a block (blockIdx.x = i), L1 cache will broadcast it? Actually each thread reads xi[d] from global memory. The first read might be uncached, but subsequent reads by other threads in the same warp for the same address will hit L1 (broadcast). For multiple warps (256 threads = 8 warps), the reads will also be cached in L1. So we might not need shared memory. But using shared memory ensures fast access and less L1 pressure. Either is fine.
Let's keep it simple: no shared memory for dist_sq kernel. Just loop and read from global memory. L1/L2 broadcast is good enough.
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
float sum = 0.0f;
for (int d = 0; d < dim; ++d) {
float diff = particles[j * dim + d] - particles[i * dim + d];
sum += diff * diff;
}
dist_sq[i * N + j] = sum;
}
```
This is simple. For i fixed, all threads read particles[i*dim + d] from the same address. The hardware broadcasts this to all threads. Similarly for each j, threads in a warp may read different j's, so the reads for xj are coalesced across j but not across d. Wait, in this kernel, each thread handles a different (i,j) pair. For a fixed block (i is fixed), threads have different j. They all execute the loop over d together. In iteration d, thread for j reads `particles[j*dim + d]`. Since j varies across threads, the addresses are `particles + j*dim + d`. For consecutive j, addresses are separated by dim=32 floats = 128 bytes. This is a strided access with stride 128 bytes. For a warp of 32 threads, the addresses are 128 bytes apart. A memory transaction on modern GPUs is 32 bytes (or 128 bytes). With a stride of 128 bytes, each thread reads a different 128-byte segment. The warp will generate 32 transactions of 128 bytes each, reading 4KB total, but only using 128 bytes (32 floats). That's terrible memory bandwidth! Because the warp is reading scattered memory.
So this kernel is actually very inefficient for memory access. Each warp reads one element per thread, but they are strided by dim (32) floats = 128 bytes. This means 32 separate cache lines per warp. That's bad.
To fix this, we need to transpose the access or change the launch geometry.
Better: have each block handle a tile of (i,j) pairs, and load a tile of particles into shared memory. Or change the loop order: each thread computes one pair and loops over d. The memory access pattern for x_j is strided by j. But we can restructure: Instead of having each block handle a single i and threads handle j, we can have each block handle a tile of i and j, and each thread compute one (i,j) pair. But the memory pattern for reading particles is still strided.
Alternatively, we can have each block handle one i and have threads loop over j in a strided way? Actually the fundamental issue is that for a fixed d, the array is laid out as (N,dim). So particles[j*dim + d] for varying j are separated by dim. This is a strided access. But it's coalesced only if consecutive threads read consecutive j. For a warp of 32 threads, j values are j0, j0+1, ..., j0+31. The addresses are base + (j0 + tid)*dim + d. The stride between consecutive threads is dim floats = 128 bytes. On A100, global memory load for a warp can coalesce multiple addresses into fewer transactions if they fall within the same cache line segment (128 bytes). But with stride 128 bytes, each address is in a different 128-byte segment. So each thread generates a separate 128-byte transaction. Memory bandwidth utilization: reading 128 bytes to get 4 bytes, i.e., 3.125% utilization. That's terrible.
However, L2 cache may hold the entire particles array (128KB). L2 bandwidth is high. But still, 32 transactions per warp is a lot.
To improve, we can read particles in a transposed manner or load a block of particles into shared memory. For example, a block of 32x32 threads could load a tile of 32 particles into shared memory, then compute distances among them. But that complicates the kernel.
But for N=1024, dim=32, total memory reads in dist_sq kernel are 2 * N * N * dim floats = 64M floats = 256MB (since each pair reads two particles). At 1 TB/s, 256MB takes 0.256 ms. With poor coalescing, effective bandwidth drops to, say, 100 GB/s (L2), taking 2.5 ms. That's still not huge, but might be slower than the matmul.
The matmul `particles @ particles.t()` reads each particle multiple times but in a highly coalesced manner via Tensor Cores or optimized shared memory tiling. It will be much faster than the naive custom kernel.
For the phi kernel, the memory access pattern is different: block i has threads 0..31 (for dim). They loop over j. For each j, thread tid reads `particles[j*dim + tid]`. Since tid varies 0..31, addresses are `base + j*dim + tid`. Consecutive threads read consecutive floats. This is perfectly coalesced! Good. And they read `particles[i*dim + tid]` also coalesced (but same for all j). So phi kernel memory access is great.
Thus, the custom dist_sq kernel has bad memory access. To fix this, we should compute dist_sq using the PyTorch matmul (or an efficient custom kernel). The PyTorch matmul has excellent memory access.
Alternatively, we can compute dist_sq by launching a kernel where each block handles one i, and threads within the block cooperate to read a block of j's particles into shared memory. But since dim=32, we can have each block read, say, 32 particles at a time. Let's design a tiled dist_sq kernel:
- Block size: (dim, block_j), e.g., (32, 32).
- Each block computes a tile of dist_sq for a fixed i (blockIdx.y) and a range of j (blockIdx.x * block_j ...).
- But grid would be (N/block_j, N).
- In the block, thread (d, local_j) loads x_j for its local_j and all d? Actually each thread loads one element: for local_j, thread d loads `particles[j*dim + d]`. So a warp of 32 threads loads an entire particle j into registers/shared memory. Then for a fixed i, we need xi. We can load xi into shared memory once per block.
- Then each thread computes diff = xj[d] - xi[d], squares it, and accumulates in a register for its (i,j) pair. After iterating over all d? Actually if block_j = 32 and dim = 32, then each thread handles one d for one j. The accumulation of diff^2 across d is done by each thread? No, each thread has only one d. We need to sum across d. We can reduce within the warp for each j. That's expensive.
Alternatively, we can have each thread compute the full sum over d for one (i,j) pair, but launch threads such that j is the x-axis and i is the y-axis. But memory access is strided.
Another efficient way: use the matmul approach but in CUDA: dist_sq = norms + norms^T - 2 * (particles @ particles^T). We can compute `particles @ particles^T` with a single line of PyTorch, which calls cuBLAS. That's the best. So we should use PyTorch for `dist_sq`.
Given that, let's compute `dist_sq` in PyTorch:
```
particles_t = particles.t() # (dim, N)
prod = particles @ particles_t # (N,N)
norms = (particles ** 2).sum(dim=1, keepdim=True) # (N,1)
dist_sq = norms + norms.t() - 2.0 * prod
dist_sq = dist_sq.clamp(min=0.0)
```
Will this match reference exactly? For typical random numbers, the difference is negligible. But there is a risk of `median` shifting. Let's test with the tool.
If this fails correctness, we have two options:
1. Chunked diff: compute dist_sq in chunks using PyTorch broadcasting but without creating the full (N,N,dim) diff. For example, process i in chunks:
```
dist_sq = torch.empty(N, N, device=particles.device, dtype=particles.dtype)
chunk = 256
for i_start in range(0, N, chunk):
i_end = min(i_start+chunk, N)
diff = particles[i_start:i_end].unsqueeze(1) - particles.unsqueeze(0) # (chunk, N, dim)
dist_sq[i_start:i_end] = (diff ** 2).sum(-1)
```
Here diff is (256, 1024, 32) = 256*1024*32 = 8,388,608 floats = 32MB. That's manageable. This is exact reference compute, no matmul trick. Chunk = 256 => 4 chunks. 4 iterations. Each diff is 32MB. This avoids the 512MB diff but is still exact. The overhead is 4 kernel launches for diff and sum. This is likely fast enough and ensures exactness.
Let's check: chunk=256, diff shape (256, 1024, 32). Memory: 256*1024*32*4 = 33,554,432 bytes = 32MB. That's fine. We can even increase chunk to 512: diff shape (512, 1024, 32) = 64MB. Still fine. For N=1024, chunk=512 => 2 chunks. 2 iterations. That would be very fast.
Let's evaluate memory: GPU has plenty of memory. 64MB is fine. So we can do chunk=512 to reduce overhead.
But we need to allocate `dist_sq` once (4MB). Then compute in 2 chunks. This is easy and exact.
Let's use this approach.
Actually, we can write a custom kernel for dist_sq that uses shared memory to achieve coalescing. But chunked PyTorch is simpler and exact.
However, the chunked approach still uses PyTorch broadcasting, which is what the reference does but in chunks. It will be somewhat slower than a single custom kernel for phi but faster than reference because we avoid the huge 512MB allocation.
Wait, the reference already does a single chunk of size N=1024, resulting in 512MB diff. Chunking into 2 chunks of 512 reduces max memory to 64MB. The compute is the same. The overhead of 2 launches is small. So this is a good approach for dist_sq.
Alternatively, we can compute dist_sq in the same custom phi kernel but store it. Actually, we don't need to store dist_sq if we don't need median. But we do need median. So we must compute it first.
Let's go with chunked PyTorch for dist_sq (exact), then custom kernel for phi.
Code for chunked dist_sq:
```
N = particles.shape[0]
dist_sq = torch.empty(N, N, device=particles.device, dtype=particles.dtype)
chunk = 512
for i_start in range(0, N, chunk):
i_end = min(i_start + chunk, N)
diff_chunk = particles[i_start:i_end].unsqueeze(1) - particles.unsqueeze(0)
dist_sq[i_start:i_end] = (diff_chunk ** 2).sum(-1)
```
But this still creates a diff tensor of size (chunk, N, dim). For chunk=512, it's 64MB. This is fine. However, we need to be careful with tolerance. Since this exactly mirrors reference compute, it's safe.
Wait, but if we have a custom kernel for phi, why do we need PyTorch for dist_sq at all? The custom phi kernel reads particles and computes phi. The only reason we need PyTorch is for median. We could compute dist_sq using the custom kernel, but as we discussed, the naive custom kernel has poor memory access. However, we can write a better custom dist_sq kernel that processes particles in a coalesced way.
Let's think about a coalesced custom dist_sq kernel.
We need each thread to compute one dist_sq[i,j]. The loop over d is sequential per thread. The memory access pattern for x_j is strided. To fix this, we can transpose the problem: let each block handle a tile of j's and all i's? Or use shared memory to load a tile of particles (for j) in a coalesced way.
Suppose we have block size = (dim, BJ), where BJ is the number of j's processed per block. For example, BJ = 32, dim = 32. Then block has 1024 threads. Each block handles a fixed range of j's (tile) and all i's? Actually we can map grid to (i, j_tile). Each block computes dist_sq for one i and a tile of BJ j's. But each thread inside the block handles one (i,j) pair and loops over d. Wait if block is 2D, we can have threadIdx.x = d (0..dim-1), threadIdx.y = local_j (0..BJ-1). Then the block computes diff for all d and local_j for a fixed i. Then we need to sum over d. Each local_j has a warp of threads (d) that compute diff^2 and can reduce using warp shuffle. Then write dist_sq[i, j] for each local_j.
This is very efficient and coalesced! Let's design it:
- Grid: (N, num_tiles_j) where num_tiles_j = (N + BJ - 1) / BJ.
- Block: (dim, BJ) = (32, 32). 1024 threads per block.
- For block (i, tile_j):
j_start = tile_j * BJ.
local_j = threadIdx.y (0..31).
j = j_start + local_j.
d = threadIdx.x (0..31).
Load xi[d] from particles[i*dim + d] (coalesced for d).
Load xj[d] from particles[j*dim + d] if j < N (coalesced for d).
diff = xj - xi.
sq = diff * diff.
Warp reduce sq across d (dim=32). Use warp shuffle to sum over 32 lanes.
After reduction, lane 0 has the sum.
If d == 0 and j < N: dist_sq[i*N + j] = sum.
This kernel has perfect coalescing for loading particles: for a fixed j, threads d=0..31 read consecutive addresses. Since there are 32 j's per block, we load 32 particles in parallel. For each j, a warp reads its particle. All warps read from particles simultaneously: warp for local_j=0 reads particles[j_start*dim + d]; warp for local_j=1 reads particles[(j_start+1)*dim + d]; etc. These addresses are separated by dim floats = 128 bytes. In the same block, warps are executing concurrently. The memory controller can serve multiple cache lines. But is it coalesced within each warp? Yes, within each warp, threads read consecutive addresses (d varies). Perfect coalescing per warp. Across warps, they read different cache lines. So total memory transactions per block per d: 32 warps * 1 transaction = 32 transactions for 32 particles, reading 128 bytes each? Actually each transaction is 128 bytes (or 32 bytes?). For A100, line size is 128 bytes. A warp reads 128 bytes (32 floats). Since the addresses are aligned to 128 bytes if dim*4 is a multiple of 128. dim=32, dim*4=128. So particle j starts at address base + j*128. So for j multiple of 1, addresses are aligned to 128 bytes. Thus each warp's read is a single 128-byte transaction! Perfect.
Similarly, reading xi: all warps need the same xi. We can load xi into shared memory: have the first warp load xi[0..31] into shared mem, then broadcast to all warps. Or each warp can read from global memory and rely on L1 broadcast. But since there are 32 warps, 32 reads of the same cache line. L1 can broadcast to all warps simultaneously if they are in the same SM. Actually L1 can handle multiple requests to the same line. But loading into shared memory is cleaner.
But wait, block size 1024 may exceed the maximum threads per block (1024 is okay). But occupancy may be low because each block uses 1024 threads, so only one block per SM can run (if SM max threads = 1024 or 1536). For A100, max threads per SM is 2048, so one block per SM (1024 threads). With 80 SMs, concurrency is 80 blocks. Grid size: N * num_tiles_j = 1024 * 32 = 32768 blocks. So total blocks 32768, 80 at a time. Each block processes 32 j's. Total warps per block: 32 warps. So total warps active: 80 * 32 = 2560 warps. That's okay.
However, each block only writes 32 floats (one per j). So total writes = 1M floats = 4MB. Fine.
But each block also loops over d? Actually in this kernel, d is the x-dimension of the block. Each thread handles one d. There is no loop over d. The warp reduction replaces the loop.
Let's write this kernel.
```
__global__ void compute_dist_sq_tiled(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int tile_j = blockIdx.y;
int d = threadIdx.x;
int local_j = threadIdx.y;
int BJ = blockDim.y;
int j = tile_j * BJ + local_j;
__shared__ float s_xi[32]; // assume dim <= 32
if (local_j == 0) {
s_xi[d] = particles[i * dim + d];
}
__syncthreads();
float xj = 0.0f;
if (j < N) {
xj = particles[j * dim + d];
}
float diff = xj - s_xi[d];
float sq = diff * diff;
// warp reduce sq across d (first dimension)
for (int offset = 16; offset > 0; offset /= 2) {
sq += __shfl_down_sync(0xffffffff, sq, offset);
}
if (d == 0 && j < N) {
dist_sq[i * N + j] = sq;
}
}
```
This kernel assumes dim <= 32. If dim > 32, we could have first dimension larger and do a block-level reduction. But for benchmark dim=32, it's perfect.
Block size: (32, 32) = 1024 threads.
Grid: (N, (N + 31)/32). For N=1024, grid = (1024, 32).
Shared memory: 32 * sizeof(float) = 128 bytes.
Memory access: For each block, threads with the same local_j but varying d read xj for a specific j. This is coalesced. All local_j values read different j's, which are spaced 128 bytes apart. Since each warp reads one cache line, and there are 32 warps, the memory controller can handle 32 requests. On A100, each SM has high L1 bandwidth. This should be efficient.
There's one issue: `__shfl_down_sync` operates within a warp. Here, the warp is composed of threads with the same `local_j` and varying `d` (since threadIdx.x varies fastest in CUDA? Actually thread indexing: threadIdx.x is the fastest varying dimension within a warp. In a 2D block of (32,32), threads are ordered linearly as `tid = threadIdx.x + threadIdx.y * blockDim.x`. So threads with consecutive `d` and same `local_j` are consecutive. A warp of 32 threads has `local_j` = constant (since 32 threads fit in the first dimension). So each warp corresponds to one `local_j`. Perfect! So `__shfl_down_sync` across the warp sums across `d` for a fixed `local_j`. Great.
So this kernel is very efficient.
However, we need to write it in `cpp_sources` / `cuda_sources`. Since it uses CUDA-specific syntax, it must be in `cuda_sources`. `load_inline` compiles `cuda_sources` with nvcc.
Let's verify: In `cuda_sources`, we can include `<torch/extension.h>` and `<cuda_runtime.h>`.
For the phi kernel, we also want to be efficient. We can use a similar tiled approach or the simple block-per-i with warp-shuffle.
Since dim=32, we can set block size = dim (32). One warp per block. Then use warp shuffle for reduction. This is extremely efficient and avoids shared memory entirely.
Let's write the phi kernel as:
```
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int d = threadIdx.x;
float xi = particles[i * dim + d];
float accum = 0.0f;
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
for (int j = 0; j < N; ++j) {
float xj = particles[j * dim + d];
float diff = xj - xi;
// reduce diff^2 across warp
float sq = diff * diff;
for (int offset = 16; offset > 0; offset /= 2) {
sq += __shfl_down_sync(0xffffffff, sq, offset);
}
float dist_sq = sq; // only correct for lane 0, but we need all lanes to have it
dist_sq = __shfl_sync(0xffffffff, dist_sq, 0);
float K = expf(-dist_sq * inv_2h);
float term = K * (-xj + diff * inv_h);
accum += term;
}
out[i * dim + d] = accum / (float)N;
}
```
Wait, does `__shfl_down_sync` work when `d` varies 0..31? Yes, because all threads in the block are in one warp (block size = 32). The `__shfl_down_sync` with mask 0xffffffff works. After the loop, lane 0 has the sum. Then `__shfl_sync(0xffffffff, dist_sq, 0)` broadcasts it to all lanes. All lanes then compute K and term. Since there is no divergence, it's perfect.
But what if dim is not exactly 32? If dim < 32, we still launch 32 threads, but threads with d >= dim should not write output. They also should not read out-of-bounds. We can add `if (d >= dim) return;` at the start. But if some threads return, the warp shuffle won't work because returned threads won't execute the shuffle. In CUDA, if a thread returns, it is considered inactive, and the remaining active threads can still execute shuffles? The mask includes all 32 lanes. The docs say: "All non-exited threads named in mask must execute the same intrinsic with the same mask." If some threads have returned (exited), they are not named in the mask? Actually they are named if the mask bit is set. If they exited, they won't execute the intrinsic, which violates the requirement. So we cannot have threads exit early in a warp that uses shuffles.
To handle dim < 32 safely, we can avoid early return and instead set `diff = 0` and `xj = 0` for d >= dim, but still participate in shuffles. Then after computing, only write if d < dim.
Similarly for dim > 32, we can't use a single warp. But we can use shared memory reduction.
Since benchmark dim=32, we can target dim=32 exactly. For other dims, we can have a generic kernel or fallback.
Let's write two kernels: one optimized for dim <= 32 (using warp shuffle), and one generic for any dim (using shared memory reduction). Actually, we can write a generic kernel that works for any dim by using shared memory, as earlier. Since the benchmark is dim=32, the performance of generic kernel is still excellent. And it handles other dims without fallback. So let's stick with the generic shared memory kernel. It is safer. We can optimize it slightly for dim=32 by making blockSize=32, but the generic code does that.
Let's revisit the generic shared memory kernel and ensure it's correct for non-power-of-two block sizes. We decided to round threads up to next power of two. For dim=32, threads=32. Good. For dim=33, threads=64. Good. Up to 1024.
But is there a limit on shared memory? 1024 * 4 = 4096 bytes. Very small.
Let's write the generic kernel clearly.
```
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
float inv_N = 1.0f / (float)N;
float accum = 0.0f;
float diff_reg = 0.0f;
for (int j = 0; j < N; ++j) {
float sq = 0.0f;
if (tid < dim) {
float xi = particles[i * dim + tid];
float xj = particles[j * dim + tid];
diff_reg = xj - xi;
sq = diff_reg * diff_reg;
} else {
diff_reg = 0.0f;
}
sdata[tid] = sq;
__syncthreads();
// reduction in shared memory
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq * inv_2h);
if (tid < dim) {
float term = K * (-particles[j * dim + tid] + diff_reg * inv_h);
accum += term;
}
__syncthreads();
}
if (tid < dim) {
out[i * dim + tid] = accum * inv_N;
}
}
```
Wait, there is a bug: `blockDim.x` is a multiple of 32 and a power of two. The reduction `for (int s = blockSize/2; s > 0; s >>= 1)` is fine. But when `tid >= dim`, `diff_reg` is set to 0. However, in the next iteration j+1, `diff_reg` for tid >= dim is 0, which is fine. But for tid < dim, `diff_reg` is updated. Since threads with tid >= dim never read it, it's fine.
One subtle bug: in the reduction loop, `sdata[tid + s]` might be out of bounds for tid + s >= blockSize? No, because tid < s ensures tid + s < 2s <= blockSize (since s starts at blockSize/2). At each step, the active threads are 0..s-1, reading from s..2s-1. Since 2s <= blockSize at the start of each iteration (except the last where s=1, reading from 1). Actually when s = blockSize/2, 2s = blockSize. So it's fine.
Another subtlety: since we use dynamic shared memory, but the kernel also uses `extern __shared__ float sdata[];`, and we launch with `threads * sizeof(float)`, it's fine.
Now, we need `blockSize` to be a power of two. In C++ host wrapper:
```
int threads = 1;
while (threads < dim) threads <<= 1;
if (threads > 1024) threads = 1024;
if (threads < 32) threads = 32; // to ensure at least one warp? Actually for dim=1, threads=1. But block size 1 is allowed but inefficient. However reduction with blockSize=1 works. But warp operations would fail. Our generic kernel uses shared memory, so it works with blockSize=1. But having blockSize < 32 is okay. However, warp-level divergence for syncthreads might be okay.
```
Wait, if blockSize=1, then `sdata[0] = sq; __syncthreads();` then `for (s=0; s>0; ...)` loop doesn't execute. `dist_sq = sdata[0];`. Works.
If blockSize=2, works. If blockSize=16, works.
But is there any issue with blockSize < 32 and syncthreads? No, __syncthreads() works for any block size.
However, CUDA spec says block size should be a multiple of warp size (32) for best occupancy? Not strictly required, but performance is best with multiples of 32 because warps are the scheduling unit. If blockSize=1, 31 threads are wasted. For dim=32, blockSize=32, perfect.
So we can compute threads as the smallest power of two >= dim, and also ensure threads >= 32? Actually if dim=1, threads=1. That's okay. But for compatibility, maybe set threads = max(32, nextPow2(dim)). Next power of two for dim=1 is 1. If we set threads = max(32, 1) = 32, then in the kernel, tid=0 computes diff_reg and sq, tid=1..31 set sq=0 and diff_reg=0. Reduction over 32 threads. Works fine. And occupancy is better. So we can do:
```
int threads = 32;
while (threads < dim) threads <<= 1;
```
This ensures threads is at least 32 and a power of two. For dim=32, threads=32. For dim=33, threads=64. For dim=1024, threads=1024. Good.
For dim > 1024, threads would be 2048, which exceeds max block size. So we fallback. But we can also handle dim > 1024 by setting threads=1024 and have each thread handle multiple dimensions with a strided loop. That would require a more complex reduction: each thread computes a partial sum of sq for its dimensions, stores in sdata[tid], then reduce. That's easy! We don't need each thread to handle only one dimension. Let's revisit.
Actually, if we set threads = 1024 (or any multiple of 32), each thread can iterate over dimensions with stride blockDim.x. For each j:
- partial_sq = 0.
- for (int d = tid; d < dim; d += blockSize) { diff = ...; partial_sq += diff*diff; }
- sdata[tid] = partial_sq.
- __syncthreads();
- reduce across all threads (blockSize may not be power of two, but we can make it 1024 which is power of two).
- After reduction, K = exp(-dist_sq*inv_2h).
- Then each thread loops again over d with stride blockSize: term = K * (-particles[j*dim + d] + diff / h); accum[d] += term.
- End for j.
- After loop, each thread writes out[i*dim + d] = accum[d] / N for its d.
But if each thread handles multiple dimensions, it needs to keep an accumulator for each dimension. Since each thread handles dimensions spaced by blockSize, the number of dimensions per thread is ceil(dim / blockSize). For dim=32 and blockSize=1024, each thread handles 1 dimension. For dim=2048 and blockSize=1024, each thread handles 2 dimensions. It can store the accumulators in a small local array? If dim is large, the number of dimensions per thread might be large. But we can just write to output directly: for each j, compute term and add to a local accumulator for that d. Since each thread handles only a few d's, it can have a small static array or just recompute. Actually it can just accumulate into a local array of size, say, MAX_DIM_PER_THREAD = (dim + blockSize - 1) / blockSize. For dim <= 1024*2 = 2048, it's 2. For dim <= 4096, it's 4. We can allocate a small array on the stack.
But to keep it simple, for the benchmark we don't need this. Let's just support dim <= 1024 with the simple approach. For dim > 1024, we fallback to PyTorch. The fallback is safe.
However, we might want to avoid fallback entirely for robustness. But the benchmark is dim=32. Fine.
Let's think about the fallback PyTorch path for phi when dim > 1024. We already derived an efficient formulation that doesn't require diff. We can write it in the host wrapper. But the host wrapper is C++. We can write C++ PyTorch code. But it's complex and not needed.
Let's structure the C++ wrapper to handle any dim with the generic kernel using block size = 1024 and strided loops. It is manageable.
Let's design the generic kernel that supports any dim with blockSize = 1024.
Actually, for dim > 1024, we can just set blockSize = 1024. Each thread handles multiple dimensions.
Let's write a robust kernel:
```
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
float inv_N = 1.0f / (float)N;
// local accumulators: we can handle up to, say, 8 dims per thread.
// But to be generic, we can just accumulate into output via atomicAdd? No.
// Better: use a per-thread register array if dim/blockSize is small.
// For general case, we could allocate shared memory of size blockSize * max_dims_per_thread, but that's complex.
// Actually, we can simply loop over j, and for each j, after getting K,
// loop over d with stride blockSize, compute term, and add it to a running sum stored in a local register array.
// Since CUDA doesn't support VLAs, we can use a fixed-size array assuming dim <= some max.
// But dim can be large. However, for large dim, we might be better off with the PyTorch fallback.
}
```
Let's just keep the simple kernel for dim <= 1024 and add a PyTorch fallback in the Python code for dim > 1024. That is robust and simple.
Actually, we can avoid the C++ fallback entirely: in Python `ModelNew.forward`, we check if `particles.shape[1] > 1024`, then compute using PyTorch ops. Otherwise call the custom kernel. This keeps the C++ wrapper simple.
Python fallback:
```
def forward(self, particles):
N = particles.shape[0]
dim = particles.shape[1]
if dim > 1024:
# fallback
score = -particles
diff = particles.unsqueeze(0) - particles.unsqueeze(1)
dist_sq = (diff ** 2).sum(-1)
median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
h = median_sq.item() / max(math.log(float(N)), 1.0)
K = torch.exp(-dist_sq / (2.0 * h))
grad_K = -diff / h * K.unsqueeze(-1)
phi = (K @ score + grad_K.sum(dim=0)) / N
return phi
else:
return svgd_cuda.forward_svgd(particles)
```
But the fallback is exactly the reference, so it's correct. However, it allocates the huge diff. But it's only for dim > 1024, which is unlikely.
Alternatively, we can compute the fallback efficiently in Python using the matmul formulation:
```
score = -particles
norms = (particles ** 2).sum(dim=1, keepdim=True)
dist_sq = norms + norms.t() - 2.0 * particles @ particles.t()
median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
h = median_sq.item() / max(math.log(float(N)), 1.0)
K = torch.exp(-dist_sq / (2.0 * h))
phi = (K.matmul(score) + (K.matmul(particles) - particles * K.sum(dim=1, keepdim=True)) / h) / N
```
This avoids huge diff. But for dim > 1024, the matmuls are large and might be memory-bound. But fallback is rare.
Let's keep the fallback simple: just call the PyTorch reference-like code. Actually, we can just copy the reference forward logic for fallback. But to avoid huge memory for dim > 1024, we can chunk it. But again, not needed.
Let's refine the whole code.
In `forward_svgd`, we need to compute dist_sq. We have two options:
a) PyTorch matmul.
b) Custom tiled kernel.
c) Chunked PyTorch.
I'd like to use the custom tiled kernel for dist_sq, because it's exact and efficient. But writing a 2D kernel is a bit more code. Let's see if we can do it.
Actually, we can combine the dist_sq computation with the phi computation into a single kernel that does two passes: first compute all dist_sq and store in a global buffer, then compute phi. But we need median between the two passes. So no.
Let's write the custom tiled dist_sq kernel.
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int tile_j = blockIdx.y;
int d = threadIdx.x;
int local_j = threadIdx.y;
int BJ = blockDim.y;
int j = tile_j * BJ + local_j;
extern __shared__ float s_xi[];
if (local_j == 0) {
s_xi[d] = particles[i * dim + d];
}
__syncthreads();
float xj = 0.0f;
if (j < N) {
xj = particles[j * dim + d];
}
float diff = xj - s_xi[d];
float sq = diff * diff;
// warp shuffle reduction across d
for (int offset = dim / 2; offset > 0; offset /= 2) {
sq += __shfl_down_sync(0xffffffff, sq, offset);
}
if (d == 0 && j < N) {
dist_sq[i * N + j] = sq;
}
}
```
Wait, `dim` is 32. `dim/2 = 16`. The shuffle reduction works across the warp. But what if dim is not a power of two? The for loop `offset = dim/2; offset > 0; offset /= 2` works only if dim is a power of two. For dim=32, it's fine. If dim is not a power of two, e.g., dim=30, then `offset=15`, `15/2=7`, etc. This reduction pattern might miss some elements. However, for threads where `d >= dim`, we set `xj = 0` and `diff = -s_xi[d]`? Wait if dim=30, thread d=30,31 would read from `particles[i*dim + d]` out of bounds. So we need to handle that. In the benchmark, dim=32. So we can assume dim is a power of two up to 32? Actually dim is 32. But to be safe, we can write a generic reduction. However, for the dist_sq kernel, if dim is not a power of two, we can just use a simple shared memory reduction instead of warp shuffle. Or we can use the generic shared memory approach for both kernels. Since both kernels need to reduce over dim, using the same generic pattern simplifies code.
Let's just use the generic shared memory kernel for dist_sq too! We can write:
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int j = threadIdx.x + blockIdx.y * blockDim.x;
if (i >= N || j >= N) return;
float sum = 0.0f;
for (int d = 0; d < dim; ++d) {
float diff = particles[j * dim + d] - particles[i * dim + d];
sum += diff * diff;
}
dist_sq[i * N + j] = sum;
}
```
Despite strided memory access, it's simple and correct. The performance may not be optimal, but N is only 1024. Let's estimate worst-case time: 1M threads, each reading dim=32 floats from two arrays with stride 128 bytes. Effective memory bandwidth maybe 50-100 GB/s. Total reads: 2 * 1M * 32 * 4 = 256MB. At 50 GB/s, that's 5ms. At 100 GB/s, 2.5ms. That might be a bit high. But L2 cache might improve it because the data is small (128KB). When reading `particles[i*dim + d]` for fixed i, all threads read the same address. L2 handles this. When reading `particles[j*dim + d]`, each thread reads a different address. Since there are 1M threads but only 1024*32=32768 unique floats, each float is read many times (by N threads for each i). The L2 cache will hold all particles. So after the first reads, subsequent reads will hit L2. Actually, each block processes a fixed i. For a given j, that particle is read by all blocks (all i). So each x_j is read N times. Since N=1024, total reads of each x_j = 1024. The L2 cache can hold all 1024 particles. So after the first few blocks, all particles are in L2. Then each block reads x_j from L2. L2 bandwidth is much higher than global memory, maybe 2-5 TB/s. So the strided access from L2 is fast. The first reads from global memory are also not terrible. So total time might be well under 1ms.
Let's not overthink dist_sq. Use the simple kernel or PyTorch matmul. Let's use the simple custom kernel to guarantee exactness.
Actually, let's just use PyTorch chunked compute for dist_sq. It's exact and probably fast enough. For chunk=512:
```
dist_sq = torch.empty(N, N, device=particles.device, dtype=particles.dtype)
chunk = 512
for i_start in range(0, N, chunk):
i_end = min(i_start+chunk, N)
diff = particles[i_start:i_end].unsqueeze(1) - particles.unsqueeze(0)
dist_sq[i_start:i_end] = (diff ** 2).sum(-1)
```
This is 2 iterations. Each diff is (512, 1024, 32), 64MB. The compute is a subtraction, square, sum. PyTorch will fuse these into one kernel? Actually `diff ** 2` is one kernel, `.sum(-1)` is another. Or maybe fused? In any case, 2 iterations * 2 kernels = 4 kernel launches. Very fast. And exact.
Let's consider doing this in Python inside `forward_svgd`. Wait, `forward_svgd` is a C++ function. We could do the chunked computation in C++ using PyTorch ATen ops:
```
auto dist_sq = torch::empty({N, N}, options);
int chunk = 512;
for (int i_start = 0; i_start < N; i_start += chunk) {
int i_end = std::min(i_start + chunk, N);
auto diff = particles.slice(0, i_start, i_end).unsqueeze(1) - particles.unsqueeze(0);
dist_sq.slice(0, i_start, i_end) = (diff * diff).sum(-1);
}
```
This is valid C++ Torch code. But we need to be careful about temporary tensors. However, it's simple.
But `particles.unsqueeze(0)` in C++: `particles.unsqueeze(0)` returns a view? Actually in libtorch, `unsqueeze(0)` returns a new tensor. `particles.slice(0, i_start, i_end)` also returns a view. The subtraction will allocate a temporary tensor of size (chunk, N, dim). That's 64MB. Then `diff * diff` also allocates a temporary? Actually `(diff * diff)` allocates a new tensor of same size (64MB). Then `.sum(-1)` reduces to (chunk, N). So peak memory is maybe 128MB. That's fine. But the computation may not be fully fused. For 2 chunks, that's 2*(sub + mul + sum) = 6 kernels. Each kernel processes 64MB. Should be fast.
However, doing this in C++ adds complexity to the `cuda_sources` string, because we need to write C++ code with proper syntax. It's doable.
Alternatively, we can compute `dist_sq` in the Python `forward` before calling `forward_svgd`. Then pass `dist_sq` and `h` to the kernel. But `forward_svgd` currently computes everything. We can split it: Python computes `dist_sq` and `h`, then calls a custom kernel for phi that just takes particles, h, N, dim. This is cleaner: the custom kernel only does phi.
Let's do that:
Python `forward`:
```
def forward(self, particles):
N = particles.shape[0]
# compute dist_sq in chunks
dist_sq = torch.empty(N, N, device=particles.device, dtype=particles.dtype)
chunk = 512
for i_start in range(0, N, chunk):
i_end = min(i_start + chunk, N)
diff = particles[i_start:i_end].unsqueeze(1) - particles.unsqueeze(0)
dist_sq[i_start:i_end] = (diff ** 2).sum(-1)
median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
h = median_sq.item() / max(math.log(float(N)), 1.0)
return svgd_cuda.forward_svgd(particles, h)
```
And `forward_svgd` in C++ takes `particles` and `h`.
This is great! It avoids C++ complexity for dist_sq, uses exact PyTorch chunked compute (fast), and the custom kernel is simpler.
But wait, is the chunked compute fast enough? Let's think. For chunk=512, diff is (512, 1024, 32). That's 512*1024*32*4 = 67,108,864 bytes = 64MB. The subtraction kernel reads particles (32KB + 64KB) and writes 64MB. The square kernel reads 64MB and writes 64MB. The sum kernel reads 64MB and writes 2MB. For 2 chunks, total memory traffic: 2 * ( (~64MB read + 64MB write) * 2 + 64MB read + 2MB write ) ~ 2 * (256MB + 66MB) = ~644MB. At 1TB/s, that's <1ms. Overhead of 6 kernel launches is negligible (~6 us each). So total maybe <1ms. Good.
Reference's single diff of 512MB would also have similar memory traffic but in one go. The advantage of chunking is avoiding a single huge allocation that might be slow due to memory management overhead, but compute is similar.
Actually, the reference's `particles.unsqueeze(0) - particles.unsqueeze(1)` creates a tensor of shape (1024,1024,32) = 134,217,728 elements = 134M floats = 512MB. That's a huge allocation. Modern GPUs have tens of GB of memory, so 512MB allocation is not a problem. But PyTorch's allocator overhead for 512MB might be non-trivial. Chunking into 64MB pieces reduces allocator overhead.
Let's go with chunked dist_sq in Python.
Now, the custom kernel for phi only takes particles and h. Let's write the C++ wrapper.
```
torch::Tensor forward_svgd(torch::Tensor particles, float h) {
int N = particles.size(0);
int dim = particles.size(1);
auto out = torch::empty_like(particles);
int threads = 32;
while (threads < dim) threads <<= 1;
// We can cap at 1024 and fallback from Python
size_t shared_mem = threads * sizeof(float);
svgd_phi_kernel<<<N, threads, shared_mem>>>(
particles.data_ptr<float>(), h, N, dim, out.data_ptr<float>());
return out;
}
```
Wait, if dim > 1024, threads will exceed 1024 and CUDA launch will fail. We handle that in Python by not calling the custom kernel.
But we also need to support dim not a power of two and less than 32. `threads` starts at 32. If dim=16, threads=32. That's fine. The kernel handles tid < dim.
Let's verify the kernel with threads=32 and dim=16. Threads 16..31 will set `sq=0` and `diff_reg=0`. Reduction sums all 32 entries. `dist_sq` will be correct because entries 16..31 contribute 0. Then K is computed. Then threads 16..31 skip the accumulation because `tid < dim` is false. At the end, they skip writing output. This is correct.
One issue: `diff_reg` for threads 16..31 is set to 0. In the term computation, they don't use it. Good.
But in the loop over j, for each j, threads 16..31 set `sq=0` and `diff_reg=0`. Then syncthreads. Then reduction. Then they skip accumulation. Then syncthreads at the end. Since they always reach the syncthreads, it's fine.
However, what if dim > threads? Our code ensures threads >= dim (since threads starts at 32 and doubles until >= dim). If dim > 1024, we fallback. So no issue.
Now, there is a subtle bug in the reduction: for threads > dim, we set sdata[tid] = 0. The reduction sums all `threads` entries. For example, dim=16, threads=32. sdata[0..15] have sq values, sdata[16..31] have 0. Reduction over 32 entries yields correct sum. Good.
Another subtlety: In the kernel, we use `particles.data_ptr<float>()`. We need to ensure tensor is contiguous and float. Since input is float32 and typical memory layout is contiguous, this is fine. We can add checks in the C++ wrapper if needed. But not required.
Let's think about performance of the phi kernel. Each block loops over j from 0 to N-1 (1024). For each j:
- Read x_i[d] and x_j[d] from global memory (coalesced).
- Compute diff and sq.
- Write to shared memory.
- __syncthreads()
- 5 rounds of reduction (for threads=32) with __syncthreads().
- Read sdata[0].
- expf.
- Read x_j[d] again (coalesced).
- Compute term and accum.
- __syncthreads().
So per iteration, there are about 7 syncthreads(). For a single warp block, __syncthreads may be compiled as a memory fence because the block is one warp. But it's still some overhead. However, for N=1024, total iterations per block = 1024. So per block, 7*1024 = 7168 syncthreads. For 1024 blocks, 7.3M syncthreads. That's okay.
Could we avoid the second __syncthreads() at the end of the j loop? Yes, if we don't use shared memory in the next iteration's read/write? Actually we do use shared memory at the start of the next iteration: `sdata[tid] = sq;`. This write could race with threads from the previous iteration still reading `sdata[0]`. But after `dist_sq = sdata[0];`, all threads read sdata[0]. Since it's a broadcast (same address for all threads), and there are no writes to sdata after the reduction until the next iteration's `sdata[tid] = sq;`, we need to ensure all threads have read sdata[0] before any thread writes sdata[tid]. In a single warp, all threads execute in lockstep, so by the time any thread reaches the next iteration's write, all threads have already read sdata[0] from the previous iteration. So we could omit the second __syncthreads() for blockSize <= 32. But for safety and generality (blockSize > 32, multiple warps), we need it. Some warps might finish the term computation and go to the next iteration's `sdata[tid] = sq;` while other warps are still reading `sdata[0]` from the previous iteration. So __syncthreads() is needed.
But for blockSize=32, we can safely omit it because it's a single warp. However, the compiler might not know that blockSize is 32 at compile time. We could use `if (blockSize > 32) __syncthreads();`? But blockSize is not known at compile time? Actually it's a runtime variable, but the compiler may not optimize it. We can add `__syncthreads();` regardless; for a single warp, it's cheap.
Actually, CUDA's `__syncthreads()` for a single warp does not require any synchronization (since warps are implicitly synchronized), but the instruction is still issued? According to some sources, `__syncthreads()` in a single-warp block has negligible overhead (just a few cycles). So we can keep it.
Now, about reading x_j twice per iteration: once for diff, once for term. We can avoid the second read by storing `xj` in a register. In the kernel, after computing `diff` and `sq`, we store `sq` in shared memory. But we can keep `xj` and `diff` in registers. Since `diff_reg` is already stored. But the code currently recomputes `-particles[j*dim + tid]` in the term. We can change to use `xj` variable:
```
float xi = particles[i * dim + tid];
float xj = particles[j * dim + tid];
float diff_reg = xj - xi;
sq = diff_reg * diff_reg;
...
float term = K * (-xj + diff_reg * inv_h);
```
This avoids the second read. We already have `xj` in a local variable. So we just need to reuse it. In the current code, `if (tid < dim)` reads `xi` and `xj` into locals. `diff_reg` is local. After reduction, `if (tid < dim)` computes term using `particles[j*dim + tid]`. Let's reuse `xj`.
Let's rewrite:
```
float accum = 0.0f;
float diff_reg = 0.0f;
float xj_reg = 0.0f;
for (int j = 0; j < N; ++j) {
float sq = 0.0f;
if (tid < dim) {
float xi = particles[i * dim + tid];
xj_reg = particles[j * dim + tid];
diff_reg = xj_reg - xi;
sq = diff_reg * diff_reg;
}
sdata[tid] = sq;
__syncthreads();
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq * inv_2h);
if (tid < dim) {
float term = K * (-xj_reg + diff_reg * inv_h);
accum += term;
}
__syncthreads();
}
```
Good.
Now, one more improvement: we can load `xi` once per block, outside the j loop, into a register or shared memory. Currently each iteration loads `particles[i*dim + tid]`, which is the same address for all j. Since it's coalesced and in L1/L2, it's cheap. But we can hoist it out:
```
float xi = 0.0f;
if (tid < dim) xi = particles[i * dim + tid];
```
Then inside the loop:
```
xj_reg = particles[j * dim + tid];
diff_reg = xj_reg - xi;
```
This avoids reloading xi N times. For N=1024, that's 32KB saved. Not huge, but easy.
Let's incorporate that.
Also, `inv_h`, `inv_2h`, `inv_N` can be computed once outside loop.
Let's finalize the kernel code.
Now, about the `compute_dist_sq_kernel`. Do we need it? Since we compute dist_sq in Python, no. So `forward_svgd` only needs the phi kernel. But wait! The task says "Replace the forward computation ... with a custom CUDA kernel." It doesn't say we can't use PyTorch for some steps. But to be safe, perhaps we should compute everything in CUDA. However, computing dist_sq in Python chunked PyTorch is fine. The main heavy lifting is the phi kernel.
But the reference's forward includes dist_sq computation as part of the forward. If we compute dist_sq in Python, the custom kernel is just a small part. Is that acceptable? The instructions say: "Replace the forward computation ... with a custom CUDA kernel." This could be interpreted as the entire forward should be replaced. But in many such tasks, using PyTorch for a small precomputation (like median) is okay, as long as the bulk is in CUDA.
Still, to maximize speedup, we want to minimize PyTorch overhead. The chunked dist_sq in Python is fine.
However, there is a catch: The reference computes the median of the flattened dist_sq. In our chunked approach, dist_sq is computed using PyTorch broadcasting and sum. This is exact. But then we call `.median()`. That's fine.
Let's think about whether we can avoid allocating dist_sq entirely by computing it in a custom kernel and then calling `torch.median` on the result. That still allocates dist_sq, but in C++. The Python chunked approach also allocates it. So no difference.
Wait, could we compute the median in the custom kernel using a selection algorithm? For 1M elements, yes, but implementing a robust median kernel is nontrivial. Better to use PyTorch.
Let's consider whether we can compute the custom phi kernel without precomputing h, by using the median heuristic but with a fixed or approximate h. The reference requires exact median. So we need h.
One more thought: Could we compute h on the fly by first computing all pairwise distances in shared memory within the block? For a block handling one i, it would need distances to all j. That's N distances. For N=1024, a block could store 1024 floats in shared memory (4KB). Then after computing all distances for a block, we could sort them in shared memory and find the median. But the median of all pairwise distances includes distances for all i,j, not just for a fixed i. The median heuristic uses the median over the full N*N matrix. So we need global median.
Unless we approximate the median by the median of distances from a single particle? That would be wrong.
Another idea: The median of pairwise distances for a set of points can be approximated by the median of distances from one point to all others, times some factor? No, not accurate.
So we must compute full dist_sq.
Let's now consider the block size for the phi kernel again. With blockSize=32, the reduction uses shared memory. But shared memory is banked. Since threads 0..31 write to sdata[0..31], there are no bank conflicts (each thread accesses a different bank). In the reduction step, when s=16, thread 0 reads sdata[0] and sdata[16]; these are in different banks. Good. When s=1, thread 0 reads sdata[0] and sdata[1], different banks. So no bank conflicts.
For blockSize=64, two warps. `sdata[tid]` is fine. But `__syncthreads()` is needed between warps.
Now, let's finalize the code.
We'll write a complete Python module.
Important: We must not include `PYBIND11_MODULE`. `load_inline` generates it from `functions=[...]`.
Also, `extra_cuda_cflags=["-O3"]`. We might want to add `--use_fast_math`? No, to keep precision. `-O3` is good.
Let's draft the code.
```python
import torch
import torch.nn as nn
import math
from torch.utils.cpp_extension import load_inline
# CUDA source
cuda_source = r'''
#include <torch/extension.h>
#include <cuda_runtime.h>
__global__ void svgd_phi_kernel(const float* particles, float h, int N, int dim, float* out) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
float inv_h = 1.0f / h;
float inv_2h = 0.5f * inv_h;
float inv_N = 1.0f / (float)N;
float xi = 0.0f;
if (tid < dim) {
xi = particles[i * dim + tid];
}
float accum = 0.0f;
float diff_reg = 0.0f;
float xj_reg = 0.0f;
for (int j = 0; j < N; ++j) {
float sq = 0.0f;
if (tid < dim) {
xj_reg = particles[j * dim + tid];
diff_reg = xj_reg - xi;
sq = diff_reg * diff_reg;
}
sdata[tid] = sq;
__syncthreads();
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
float dist_sq = sdata[0];
float K = expf(-dist_sq * inv_2h);
if (tid < dim) {
float term = K * (-xj_reg + diff_reg * inv_h);
accum += term;
}
__syncthreads();
}
if (tid < dim) {
out[i * dim + tid] = accum * inv_N;
}
}
torch::Tensor forward_svgd(torch::Tensor particles, float h) {
int N = particles.size(0);
int dim = particles.size(1);
auto out = torch::empty_like(particles);
int threads = 32;
while (threads < dim) threads <<= 1;
// Guard: if dim too large, we shouldn't reach here because Python checks.
// But cap at 1024 to be safe.
if (threads > 1024) threads = 1024;
size_t shared_mem = threads * sizeof(float);
svgd_phi_kernel<<<N, threads, shared_mem>>>(
particles.data_ptr<float>(), h, N, dim, out.data_ptr<float>());
return out;
}
'''
svgd_cuda = load_inline(
name="svgd_cuda",
cpp_sources="",
cuda_sources=cuda_source,
functions=["forward_svgd"],
extra_cuda_cflags=["-O3"],
verbose=False,
)
class ModelNew(nn.Module):
def __init__(self, dim):
super().__init__()
self.dim = dim
self.register_buffer("target_mean", torch.zeros(dim))
def _score_fn(self, x: torch.Tensor) -> torch.Tensor:
return -x
def forward(self, particles: torch.Tensor) -> torch.Tensor:
N = particles.shape[0]
dim = particles.shape[1]
# Fallback to PyTorch reference for very large dim (not expected in benchmark)
if dim > 1024:
score = -particles
diff = particles.unsqueeze(0) - particles.unsqueeze(1)
dist_sq = (diff ** 2).sum(-1)
median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
h = median_sq.item() / max(math.log(float(N)), 1.0)
K = torch.exp(-dist_sq / (2.0 * h))
grad_K = -diff / h * K.unsqueeze(-1)
phi = (K.matmul(score) + grad_K.sum(dim=0)) / N
return phi
# Compute pairwise squared distances in chunks to avoid huge intermediate
dist_sq = torch.empty(N, N, device=particles.device, dtype=particles.dtype)
chunk = 512
for i_start in range(0, N, chunk):
i_end = min(i_start + chunk, N)
diff = particles[i_start:i_end].unsqueeze(1) - particles.unsqueeze(0)
dist_sq[i_start:i_end] = (diff ** 2).sum(-1)
median_sq = dist_sq.view(-1).median().clamp(min=1e-5)
h = median_sq.item() / max(math.log(float(N)), 1.0)
return svgd_cuda.forward_svgd(particles, h)
```
Let's double-check the fallback and chunk logic.
In the fallback, `diff` is huge for dim > 1024, but that case is unlikely. We could also use the efficient PyTorch formulation for fallback. But it's okay.
In the chunk loop, `diff` shape is (chunk, N, dim). Since N=1024, chunk=512, diff is 64MB. This is fine.
The output of `forward_svgd` is `out`, same shape as particles.
One more thing: `median_sq = dist_sq.view(-1).median().clamp(min=1e-5)`. In PyTorch, `median()` returns a tensor, then `.clamp(min=1e-5)` returns a tensor. Then `.item<float>()` gives float.
What about the device? `particles` is on GPU. `dist_sq` is on GPU. `median` is computed on GPU. `h` is a Python float. Then passed to `forward_svgd`.
In the C++ wrapper, `h` is passed as float. Good.
Now, we should test correctness with the tool.
Let's call `compile_kernel` first.
Wait, the tool `compile_kernel` expects a complete Python source of ModelNew kernel file. The same as what we will submit. So we can pass the whole file string to compile and correct.
Let's draft the string and compile it.
But before that, let's think about potential compilation issues.
1. `load_inline` compiles at module load time. The compilation may take a few seconds. But it's done when the file is imported.
2. `extra_cuda_cflags=["-O3"]` is okay.
3. `cpp_sources=""` is okay.
4. The `cuda_sources` string must be valid CUDA C++. The `#include <torch/extension.h>` is needed for `torch::Tensor`. `at::Tensor` also works.
5. We defined `forward_svgd`. `load_inline` will generate pybind11 binding for it.
One potential issue: `torch::empty_like(particles)` returns a tensor with the same device and dtype. Good.
Another issue: if `particles` is not contiguous or not float32, `data_ptr<float>()` may be unsafe. But benchmark uses float32 contiguous.
Another potential issue: The lambda/function name `forward_svgd` is provided in `functions=["forward_svgd"]`. Good.
Let's also consider whether we need to include `#include <cmath>` for `expf` and `logf`. `<cuda_runtime.h>` might include it. But we can add `#include <cmath>` for safety. Actually `expf` is from `<math.h>`, which is included by CUDA. In device code, `expf` is an intrinsic. So it's fine.
Now, we should think about whether the chunked dist_sq in Python is the best approach or if we should use the PyTorch matmul. Let's consider the possibility that the tool output might show our kernel is not fast enough if the Python overhead is high. Chunked dist_sq in Python is 2 chunks * (unsqueeze + sub + pow + sum) = 2*4 = 8 PyTorch ops. Each op has kernel launch overhead (~5-10 us). Total overhead ~80 us. Plus compute ~0.5 ms. Plus median ~0.1 ms. Plus custom kernel ~0.1 ms. Total ~0.8 ms. Reference might take ~5 ms because of the 512MB allocation. So speedup ~6x.
But what if reference is actually not that slow? PyTorch might optimize the diff allocation and ops. On a fast GPU, reference might be ~1-2 ms. Then our speedup might be ~2x. That's still fine.
To get even more speedup, we could compute dist_sq with a single custom kernel (tiled for coalescing) to reduce Python overhead. Let's see if we can write a simple custom dist_sq kernel with coalesced memory access.
Coalesced dist_sq kernel design:
- Have each block compute a tile of (i, j) pairs.
- Use 2D block: (dim, BJ). For dim=32, BJ=32.
- Each thread computes diff for one (i, j, d).
- Then warp-shuffle reduce across d to get dist_sq.
- Write result.
This avoids loops and strided access. Let's write it.
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int tile_j = blockIdx.y;
int d = threadIdx.x;
int local_j = threadIdx.y;
int BJ = blockDim.y;
int j = tile_j * BJ + local_j;
// load xi into shared memory
extern __shared__ float s_xi[];
if (local_j == 0) {
s_xi[d] = particles[i * dim + d];
}
__syncthreads();
float xj = 0.0f;
if (j < N) {
xj = particles[j * dim + d];
}
float diff = xj - s_xi[d];
float sq = diff * diff;
// Reduce across d (warp shuffle)
// Since blockDim.x = dim = 32, each warp is a column of local_j.
for (int offset = dim / 2; offset > 0; offset >>= 1) {
sq += __shfl_down_sync(0xffffffff, sq, offset);
}
if (d == 0 && j < N) {
dist_sq[i * N + j] = sq;
}
}
```
Wait, `__shfl_down_sync` reduces within the warp. For block (32,32), the linear thread ID is `d + local_j * dim`. A warp is 32 consecutive threads. Since `dim=32`, `local_j` changes every 32 threads. So a warp consists of threads with the same `local_j` and `d=0..31`. So warp shuffle across d works perfectly.
But if dim is not 32, this kernel needs to be adapted. For generality, we can assume dim <= 32 and blockDim.x = dim? But dim might be 16. Then a warp has 32 threads: 16 threads for local_j=0 (d=0..15), and 16 threads for local_j=1 (d=0..15). Then `__shfl_down_sync` across the warp would mix d values from different local_j. That's wrong.
So this kernel only works if blockDim.x is exactly the warp size (32) or if dim = 32. To make it general, we need to ensure that each warp handles exactly one local_j. That means blockDim.x must be <= 32 and a divisor of 32? Actually the warp size is 32. Threads are ordered by x then y. For a 2D block (x, y), a warp of 32 threads will have x varying fastest. So if blockDim.x = dim, then a warp covers `dim` threads in x and `ceil(32/dim)` threads in y. For dim=32, it's exactly one local_j. For dim=16, a warp covers local_j=0..1 (16*2=32). Then warp shuffle would mix local_j values. To avoid that, we can set blockDim.x = 32 (warp size) and have each thread handle a subset of dimensions? Or we can loop over d with stride 32. But that defeats the purpose.
Alternatively, we can set blockDim.x = 32 always, and each thread computes diff for a specific d mod 32? Not ideal.
For simplicity, since benchmark dim=32, we can write the dist_sq kernel specialized for dim=32, or just use the generic 1D kernel with strided access, or PyTorch chunked. Let's stick with PyTorch chunked or matmul. It avoids writing a complex kernel.
Actually, wait: can we compute dist_sq using a simple custom kernel but with better memory access by having each block handle one i and threads loop over j? Like:
```
__global__ void compute_dist_sq_kernel(const float* particles, int N, int dim, float* dist_sq) {
int i = blockIdx.x;
int tid = threadIdx.x;
int blockSize = blockDim.x;
extern __shared__ float sdata[];
// sdata size = blockSize (for reduction)
// Also s_xi for dimensions? We can store xi in registers per thread if each thread handles one d.
// But if blockSize < dim, each thread handles multiple d.
}
```
This is getting messy.
Best approach: Use PyTorch matmul for dist_sq. It's the fastest and simplest. Let's assume it will be within tolerance.
But to be safe, we can compute dist_sq using PyTorch's `torch.cdist` which computes pairwise distances in a stable and optimized way:
```
dist_sq = torch.cdist(particles, particles, p=2).pow(2)
```
`cdist` might use a block algorithm and be very fast. However, `cdist` computes sqrt of dist_sq and then we square it, which could introduce small errors. Also, for 1024x32, `cdist` is efficient.
But the reference uses sum of squared differences. `cdist` might produce slightly different results due to algorithmic differences (e.g., using numerically stable accumulation). Tolerance is 1e-4, which is quite loose for fp32. So likely fine.
Let's just use matmul: `dist_sq = norms + norms.t() - 2 * particles @ particles.t()`. Then clamp to non-negative: `dist_sq = dist_sq.clamp(min=0.0)`. This prevents negative values from causing NaN in sqrt or log, but here it's just exponential.
Actually, we should `clamp(min=0.0)` because negative dist_sq due to rounding will cause `exp(-negative / (2h))` > 1, which is wrong. The reference never has negative dist_sq. So clamping is necessary for correctness.
Let's do this in Python:
```
norms = (particles ** 2).sum(dim=1, keepdim=True)
dist_sq = norms + norms.t() - 2.0 * (particles @ particles.t())
dist_sq = dist_sq.clamp(min=0.0)
```
This avoids any large intermediate. Then median, h, then custom kernel.
Let's test this.
But we need to be careful: the matmul approach might produce slight differences for the median. However, we can use the chunked diff approach to compute dist_sq exactly if needed. But matmul is faster. For the benchmark, speedup is important. Let's try matmul first.
Actually, the chunked diff approach is also fast and exact. Let's compare:
Matmul:
- norms: reads particles (128KB), writes norms (4KB).
- matmul: reads particles twice? Actually `particles @ particles.t()` reads particles (128KB) and the transposed view. It uses cuBLAS, highly optimized. Writes 4MB.
- addition: reads norms (4KB), norms.t (4KB), prod (4MB), writes dist_sq (4MB).
Total memory: small. Compute: matmul is ~33M ops. Very fast.
Chunked diff:
- chunk=512: diff is subtraction, reads 96KB, writes 64MB. square writes 64MB. sum writes 2MB. Times 2. Total memory: ~260MB. Compute: trivial. Launch overhead: 6 launches.
Matmul is almost certainly faster. Let's use matmul.
However, there is a potential issue: if particles has shape (N,dim), `particles @ particles.t()` creates a (N,N) tensor. The operation `norms + norms.t() - 2.0 * prod` will allocate a temporary for `2.0 * prod` (4MB) and for the subtraction/addition. But that's fine. PyTorch might fuse some of these with the JIT? Not necessarily. But even with temporaries, total memory is ~16MB. Very fast.
Let's finalize with matmul for dist_sq.
Wait, there is one more nuance: The reference computes `dist_sq = (diff ** 2).sum(-1)`. This is `sum_d (x_j[d] - x_i[d])^2`. The matmul computes `||x_i||^2 + ||x_j||^2 - 2 x_i·x_j`. For very close points, cancellation can cause slight negative values, which we clamp. The reference would produce exactly 0 for identical points (since diff=0). The matmul might produce a small negative value, which we clamp to 0. So for identical points, both give 0. Good.
For points that are close but not identical, the cancellation might produce a small relative error. But squared distances are ~O(1-100). The dot product is O(30). The subtraction of two ~30 numbers to get a ~1 result can lose a few bits of precision. The relative error might be ~1e-6. That's well within tolerance. So matmul is fine.
Let's finalize the code.
We also want to make sure the C++ code compiles. `extra_cuda_cflags=["-O3"]`. Should we include `-std=c++17`? `load_inline` defaults to C++14 or C++17 depending on PyTorch. Adding `-O3` is fine.
Let's write the final Python source.
One more thought: Could we compute h inside the C++ wrapper by calling a CUDA reduction for median? No, use PyTorch median.
Let's also check the `target_mean` buffer. It's not used, but we keep it.
Another detail: In the reference, `score = self._score_fn(particles)` returns `-particles`. Our ModelNew doesn't need `_score_fn`. But it's okay to omit it. The task only requires same constructor args and forward signature.
Let's write the code and compile it.
I'll craft the string carefully.
Actually, I want to consider whether `load_inline` can be called at module level. Yes, it compiles on import. We should wrap it in a try/except? No, that would be suspicious. But compilation errors will show up in `compile_kernel` output.
Let's just write it.
Also, I want to double-check the formula inside the kernel.
Reference: `phi = (K @ score + grad_K.sum(dim=0)) / N`.
`K @ score`: (N,N) @ (N,dim) => output row i = sum_j K_{i,j} * score_j.
`grad_K.sum(dim=0)`: grad_K shape (N,N,dim). dim=0 sum => output row j = sum_i grad_K_{i,j}.
So phi_i = (sum_j K_{i,j} * score_j + sum_j grad_K_{j,i}) / N? Wait earlier we concluded output of grad_K.sum(dim=0) is indexed by j (source). Let's re-evaluate with indices.
In PyTorch:
- `diff` shape (N,N,dim). diff[i,j] = x_j - x_i.
- `K` shape (N,N). K[i,j] = exp(-dist_sq[i,j]/(2h)).
- `grad_K` shape (N,N,dim). grad_K[i,j] = -(x_j - x_i)/h * K[i,j].
- `K @ score`: matrix multiplication of K (N,N) and score (N,dim). Result shape (N,dim). The row index of the result corresponds to the first index of K. So result[i,d] = sum_j K[i,j] * score[j,d].
- `grad_K.sum(dim=0)`: sum over first dimension (size N). Result shape (N,dim). The first dimension of the result corresponds to the second dimension of grad_K. So result[j,d] = sum_i grad_K[i,j,d].
- They add these two (N,dim) tensors and divide by N.
So the output index is `i` for the first term and `j` for the second term? Wait, when adding two tensors of shape (N,dim), they are added elementwise. So output[k,d] = (K@score)[k,d] + (grad_K.sum(0))[k,d] for k = 0..N-1.
(K@score)[k,d] = sum_j K[k,j] * score[j,d].
(grad_K.sum(0))[k,d] = sum_i grad_K[i,k,d].
So output[k] = (sum_j K[k,j] * score_j + sum_i grad_K[i,k]) / N.
We earlier derived that sum_i grad_K[i,k] = sum_i (x_k - x_i)/h * K[i,k] = sum_j (x_k - x_j)/h * K[j,k]. Wait grad_K[i,k] = -(x_k - x_i)/h * K[i,k] = (x_i - x_k)/h * K[i,k]. So sum_i grad_K[i,k] = sum_i (x_i - x_k)/h * K[i,k] = (1/h) (sum_i x_i K[i,k] - x_k sum_i K[i,k]).
So phi_k = (sum_j K[k,j] * score_j + sum_i (x_i - x_k)/h * K[i,k]) / N.
Since K symmetric, K[i,k] = K[k,i]. So:
phi_k = (sum_j K[k,j] * score_j + sum_j (x_j - x_k)/h * K[k,j]) / N.
= (1/N) sum_j K[k,j] * (score_j + (x_j - x_k)/h).
This matches our kernel formula: for output i (which is k), sum over j:
term = K_{i,j} * (score_j + (x_j - x_i)/h).
= K * (-x_j + (x_j - x_i)/h).
So our kernel computes exactly this.
Let's verify numerically with a small example to be absolutely sure.
Let's write a small script mentally:
x = [[1,2], [3,4]] (N=2, dim=2)
score = [[-1,-2], [-3,-4]]
diff[0,0] = [0,0], diff[0,1] = [2,2], diff[1,0] = [-2,-2], diff[1,1]=[0,0]
dist_sq = [[0,8],[8,0]]
K = [[1, e^{-4/h}], [e^{-4/h}, 1]]
grad_K[0,1] = -[2,2]/h * K[0,1] = [-2/h, -2/h] * K01
grad_K[1,0] = -[-2,-2]/h * K[1,0] = [2/h, 2/h] * K01
grad_K.sum(dim=0) = [grad_K[0,0]+grad_K[1,0], grad_K[0,1]+grad_K[1,1]] = [[0,0] + [2/h,2/h]*K01, [-2/h,-2/h]*K01 + [0,0]] = [[2/h,2/h]*K01, [-2/h,-2/h]*K01]
K @ score = [ [1*(-1,-2) + K01*(-3,-4)], [K01*(-1,-2) + 1*(-3,-4)] ]
phi_0 = ( [-1,-2] + K01*[-3,-4] + [2/h,2/h]*K01 ) / 2
phi_1 = ( K01*[-1,-2] + [-3,-4] + [-2/h,-2/h]*K01 ) / 2
Now our kernel for i=0, j=0:
diff = x0 - x0 = [0,0], K=1.
term = 1 * (-x0 + 0/h) = [-1,-2].
j=1:
diff = x1 - x0 = [2,2], K=K01.
term = K01 * (-x1 + [2,2]/h) = K01 * ( [-3,-4] + [2/h, 2/h] ) = K01*[-3,-4] + K01*[2/h,2/h].
sum = [-1,-2] + K01*[-3,-4] + K01*[2/h,2/h]. -> matches phi_0 * 2? Wait phi_0 = sum / N = sum / 2. Yes, sum is 2*phi_0. Good.
So kernel is correct.
Great.
Now, let's write the final code string and compile.
I will also consider adding `#include <cuda.h>` or similar. `<torch/extension.h>` is sufficient.
Let's also think about float vs double. The input is fp32. `h` is computed as float (via `.item<float>()`). In C++, we pass `float h`. Good.
One final thing: The `load_inline` call might produce a compilation warning if `cpp_sources` is empty. It should be okay.
Another potential issue: In the C++ wrapper, if we pass `particles` which is non-contiguous, `data_ptr<float>()` still works but strides are ignored because we assume contiguous. We can ensure contiguity: `particles = particles.contiguous();`. But likely it's already contiguous. We can add `TORCH_CHECK(particles.is_contiguous());` but not necessary.
Let's add a check that particles is float32: `TORCH_CHECK(particles.scalar_type() == torch::kFloat32);`.
We should also handle the device: `particles.device().is_cuda()`. Not necessary for benchmark.
Let's write the code now.