reasoning
This is a 3D transposed convolution. Let me analyze the operation and think about optimization strategies.
First, let me understand the reference:
- Input shape: (batch_size, in_channels, depth, height, width)
- Kernel shape: (out_channels/groups, in_channels/groups, kd, kh, kw)
- Output shape depends on stride, padding, output_padding
Given test parameters:
- batch_size = 8
- in_channels = 32
- out_channels = 32
- kernel_size = (3, 5, 7)
- depth = 12
- height = 24
- width = 48
- stride = (2, 2, 2)
- padding = (1, 2, 3)
- output_padding = (1, 1, 1)
- groups = 4
For ConvTranspose3d, the formula for output size is:
out_size = (in_size - 1) * stride - 2 * padding + kernel_size + output_padding
D: (12 - 1) * 2 - 2 * 1 + 3 + 1 = 22 - 2 + 3 + 1 = 24
H: (24 - 1) * 2 - 2 * 2 + 5 + 1 = 46 - 4 + 5 + 1 = 48
W: (48 - 1) * 2 - 2 * 3 + 7 + 1 = 94 - 6 + 7 + 1 = 96
So output shape: (8, 32, 24, 48, 96)
Total output elements: 8 * 32 * 24 * 48 * 96 = 28,311,552
Groups = 4, so each group has 8 in_channels and 8 out_channels.
For transposed convolution with stride > 1, a common optimization is to decompose it into a direct convolution. However, the most straightforward approach is:
1. Expand the input by stride (insert zeros between elements)
2. Apply a regular convolution (dilation/stride=1)
3. Crop/pad the result
But ConvTranspose3d is equivalent to:
- For each output element, we compute a linear combination of kernel weights and corresponding input elements
- With stride, each input element contributes to multiple output elements
Actually, another way to think about it: transposed convolution is equivalent to a normal convolution with dilated input. When stride=s, each input pixel is separated by s-1 zeros in the expanded input.
A more efficient way: instead of materializing the expanded input, we can directly map from output positions to input positions.
Let me think about the GPU kernel approach:
One approach is to compute the output directly:
- For each output pixel (b, oc, od, oh, ow)
- We need to find all input pixels that contribute to it
- In transposed conv with stride, each input pixel at (id, ih, iw) contributes to output pixels in the range [id*s, id*s + kd-1], etc.
- So for output pixel (od, oh, ow), input pixel (id, ih, iw) contributes iff:
id = floor((od - kd + 1 + 2*pd) / s) ... no, this is getting complicated.
Actually, the standard formulation:
- The expanded input has size D' = (D-1)*s_d + 1 + 2*pd
- Wait, for transposed conv:
- Expand input: size becomes (D-1)*stride + 1
- Pad this with `padding` on both sides
- Apply conv with kernel_size, stride=1
- The output size is (D-1)*stride - 2*padding + kernel_size + output_padding
Actually let me re-derive:
If we have input of size D, and want transposed conv output:
- The dilated (expanded) input has size D_dilated = (D - 1) * stride + 1
- After removing padding (equivalent to cropping): D_crop = D_dilated - 2*padding
- After conv with kernel size k, stride=1: D_out = D_crop + k - 1
- After adding output_padding: D_out = D_out + output_padding
So: D_out = (D - 1) * stride - 2*padding + k + output_padding
For our test:
D: (12-1)*2 - 2*1 + 3 + 1 = 22 - 2 + 4 = 24 ✓
Now, the operation can be viewed as:
1. Expand input with zeros (insert stride-1 zeros between each element)
2. Pad with `padding` amount of zeros on both sides
3. Perform a regular 3D convolution with stride=1
4. Pad the output with output_padding zeros (usually on the right/bottom)
So an efficient kernel could:
- Directly compute output without materializing the expanded input
- For each output position, find which input positions and kernel positions contribute
For output position (od, oh, ow), the contributing expanded input position would be:
- (id_dilated, ih_dilated, iw_dilated) where the kernel overlaps
- id_dilated = od - kd_start, where kd_start ranges from 0 to kd-1
- But we only care about positions where id_dilated corresponds to actual (non-zero) input
In the expanded input, non-zero positions are at: i * stride for i in [0, D-1]
So id_dilated = id * stride_d for some id in [0, D-1]
For a given output od and kernel offset kd_offset:
id_dilated = od - kd_offset
This is valid if id_dilated % stride_d == 0 and 0 <= id_dilated / stride_d < D
And we need: 0 <= id_dilated < (D-1)*stride + 1
Actually, let me be more careful. If we think of expanded input indices:
- Expanded input size: (D-1)*stride + 1
- Valid non-zero positions: 0, stride, 2*stride, ..., (D-1)*stride
- After padding with `padding` on each side, the expanded+pad input size is: (D-1)*stride + 1 + 2*padding
- The conv computes output at positions od where the kernel window [od, od+kd-1] is within the padded input
- Output position od in the conv result corresponds to expanded+pad input window starting at od
Wait, I need to be careful about indexing. Let me think again.
For a regular conv with input size I, kernel size k, padding p, stride 1:
Output size = I + 2p - k + 1
With I = (D-1)*stride + 1 (expanded input without padding)
And padding = p:
Padded input size = I + 2p
Output size = I + 2p - k + 1 = (D-1)*stride + 1 + 2p - k + 1
But we add output_padding to get the final size.
Actually for ConvTranspose3d, I think the implementation typically:
1. Creates expanded input with zeros inserted
2. Pads with `padding` zeros
3. Convolve with stride=1
4. The output is then trimmed/expanded by output_padding
For a simple direct kernel, I can:
- For each output element, iterate over all possible kernel positions that could contribute
- For each kernel position (kd_offset, kh_offset, kw_offset):
- Compute the corresponding expanded input position: (od - kd_offset, oh - kh_offset, ow - kw_offset)
- Check if this is a valid non-zero position in the expanded input (i.e., divisible by stride and in range)
- If so, map to actual input position and accumulate
Let me formalize this. With expanded+pad input:
- The padded positions range from -padding to (D-1)*stride + padding (conceptually)
- Actually more precisely, the padded input has indices from 0 to (D-1)*stride + 2*padding
- The original expanded input starts at index `padding` and ends at index (D-1)*stride + padding
For an output index od in [0, D_out-1]:
The conv window covers input positions [od, od+kd-1] in the padded input.
A kernel offset (k_d, k_h, k_w) means the kernel element at this offset multiplies the input at position (od + k_d, oh + k_h, ow + k_w) in the padded input.
Wait, no. For conv, output at position o is sum over k of input[o+k] * weight[k].
Hmm, but we need to be consistent. Actually in transposed conv, when viewed as conv on expanded input:
- output[o] = sum_k expanded_input[o + k] * weight[k]
Where expanded_input is the zero-padded, dilated input.
For expanded_input to be non-zero at position p:
- p must be in [padding, (D-1)*stride + padding] (the range of actual dilated input in the padded array)
- (p - padding) % stride == 0
And the actual input index is (p - padding) / stride.
So for output position (od, oh, ow), kernel offset (kd, kh, kw):
p_d = od + kd
p_h = oh + kh
p_w = ow + kw
This contributes if:
- padding_d <= p_d <= (D-1)*stride_d + padding_d
- (p_d - padding_d) % stride_d == 0
- And similarly for h and w
And the actual input index is ((p_d - padding_d) / stride_d, ...)
The weight used is weight[oc, ic, kd, kh, kw] for the appropriate group.
Wait, let me double-check the weight indexing. In ConvTranspose3d:
- weight shape is (in_channels, out_channels/groups, k_d, k_h, k_w)
Hmm, actually I should check. For nn.ConvTranspose3d, the weight shape is:
(in_channels, out_channels/groups, kd, kh, kw) when groups=1
For groups > 1: (in_channels, out_channels/groups, kd, kh, kw)
Actually, let me verify. PyTorch ConvTranspose3d weight is shape:
(in_channels, out_channels // groups, *kernel_size)
Wait no, actually I think for ConvTranspose3d:
- weight shape: (in_channels, out_channels // groups, *kernel_size)
Hmm that's different from Conv3d which is (out_channels, in_channels // groups, *kernel_size).
Actually, I should verify this in my head. For a transposed convolution, if we think of the operation as the transpose of a convolution matrix, the weight dimensions would be swapped.
Let me just check: for ConvTranspose3d, if input has C_in channels and output has C_out channels, then each group has C_in/groups input channels and C_out/groups output channels. The weight for each group has shape (C_in/groups, C_out/groups, *kernel_size).
So total weight shape: (groups, C_in/groups, C_out/groups, *kernel_size) or flattened as (C_in, C_out/groups, *kernel_size).
Actually, in PyTorch, ConvTranspose3d weight is stored as (in_channels, out_channels/groups, *kernel_size).
Wait, I'm second-guessing myself. Let me think about this more carefully.
For Conv3d with groups:
- weight shape: (out_channels, in_channels/groups, k_d, k_h, k_w)
- Each output channel oc corresponds to group g = oc // (out_channels/groups)
- The weight for oc is weight[oc, :, :, :, :] which has shape (in_channels/groups, k_d, k_h, k_w)
- This is convolved with input channels [g*(in_channels/groups) : (g+1)*(in_channels/groups)]
For ConvTranspose3d with groups, the operation is the "transpose" of the above. If the corresponding conv would map from in_channels to out_channels, then transposed conv maps from out_channels to in_channels.
Wait, actually the naming can be confusing. In ConvTranspose3d:
- The first argument is in_channels, which is the number of channels in the input
- The second argument is out_channels, which is the number of channels in the output
- So it maps from in_channels → out_channels, just like regular conv
But the weight shape for ConvTranspose3d is:
(in_channels, out_channels // groups, *kernel_size)
So for an output channel oc in group g:
- The relevant weights are weight[g*(in_channels/groups) : (g+1)*(in_channels/groups), oc - g*(out_channels/groups), :, :, :]
- This has shape (in_channels/groups, k_d, k_h, k_w)
- And it's applied to input channels in group g
Actually wait, that doesn't make sense dimensionally. Let me reconsider.
In ConvTranspose3d, if we have input x with shape (N, C_in, D, H, W):
- For a regular conv, output channel j is computed from all input channels in its group using weight[j, :, :, :, :]
- For transposed conv, output channel j should be computed similarly but from the transposed operation
I think a clearer way: ConvTranspose3d is equivalent to:
1. For each input channel i, we have a set of filters that produce all output channels in the same group
2. Specifically, input channel i (in group g) contributes to output channels in group g using the filter weight[i, :, :, :, :]
Yes, that matches the PyTorch weight shape (in_channels, out_channels//groups, *kernel_size).
So in the kernel, for each output channel oc in group g:
oc_group_local = oc - g * (out_channels/groups)
We sum over input channels ic in group g:
ic_group_local = ic - g * (in_channels/groups)
Accumulate input[g, ic_local, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
Let me now think about the kernel design. The output shape for our test is (8, 32, 24, 48, 96).
I can write a CUDA kernel where:
- Each thread computes one output element
- Or each thread block computes a tile of output elements
- We load input and weight into shared memory for reuse
Given that we have:
- batch_size = 8
- out_channels = 32
- D_out = 24, H_out = 48, W_out = 96
Total output elements: 8 * 32 * 24 * 48 * 96 = 28,311,552
If each thread computes one output element, we need ~28M threads. With 1024 threads per block, that's ~27K blocks. Actually that's fine.
But each output element requires:
- For groups=4, each group has 8 out_channels and 8 in_channels
- Kernel size = 3*5*7 = 105
- So each output element needs up to 8*105 = 840 multiply-adds
- Total FLOPs: 28,311,552 * 8 * 105 * 2 ≈ 47.6 GFLOPs
At 67 TFLOPS (FP32 peak), this would take ~0.7 ms if compute bound. Memory bandwidth is 3.35 TB/s.
Input size: 8 * 32 * 12 * 24 * 48 * 4 bytes ≈ 14.1 MB
Output size: 8 * 32 * 24 * 48 * 96 * 4 bytes ≈ 113.2 MB
Weight size: 32 * 8 * 3 * 5 * 7 * 4 bytes ≈ 107.5 KB
Total memory traffic if naive: ~127 MB, which at 3.35 TB/s is ~0.038 ms.
So this is definitely compute bound, not memory bound.
However, the naive PyTorch implementation might not be perfectly optimized for this specific case. A custom kernel could potentially:
1. Use shared memory for weight reuse
2. Use vectorized loads when possible
3. Optimize the specific case (groups, kernel size, etc.)
But wait - PyTorch's ConvTranspose3d might use cuDNN which is usually very well optimized. Let me think about how to beat it.
One important insight: since the stride is (2, 2, 2), each output pixel only receives contributions from a sparse set of input pixels. For each spatial dimension, out of kernel_size elements, only about ceil(kernel_size / stride) or similar number actually hit non-zero dilated input positions.
Let me compute more carefully. For output position od and kernel offset kd:
p_d = od + kd - padding_d ... wait, I need to be careful about the formula.
Let me re-derive with a concrete example. Suppose:
- Input size D = 12
- Stride = 2
- Padding = 1
- Kernel size = 3
- Output_padding = 1
Expanded input (0-indexed): positions 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22
Size = (12-1)*2 + 1 = 23
Pad with 1 on each side: positions -1 to 23 (with zeros at -1 and 23)
Size = 23 + 2 = 25
Conv with kernel size 3, stride 1:
Output size = 25 - 3 + 1 = 23
Add output_padding = 1: final output size = 24
Wait, that doesn't match. Let me recheck.
PyTorch formula: (D - 1) * stride - 2*padding + kernel_size + output_padding
= (12-1)*2 - 2*1 + 3 + 1 = 22 - 2 + 4 = 24 ✓
Hmm, so where does the 24 come from?
25 - 3 + 1 = 23, plus output_padding = 1 gives 24? But output_padding is added to the end.
Actually I think the output_padding means we pad the output on one side. So the conv produces 23 elements and we append 1 more? But that doesn't sound right for the formula to be exact.
Wait, actually for ConvTranspose3d with stride=2, the output should be roughly 2x the input size. 12*2 = 24, which matches.
Let me re-examine. If we think of the conv on expanded input:
Expanded input: 23 elements with non-zero at even positions (0, 2, 4, ..., 22)
Pad with 1: now we have 25 elements, with actual data from index 1 to 23
Wait no, padding is applied before expansion, or after? Let me check PyTorch docs.
According to PyTorch docs for ConvTranspose3d:
The output spatial dimensions are computed as:
out = (in - 1) * stride - 2 * padding + dilation * (kernel_size - 1) + output_padding + 1
Oh! I was using the wrong formula. Let me recalculate:
With dilation=1:
out = (12 - 1) * 2 - 2 * 1 + 1 * (3 - 1) + 1 + 1
= 11 * 2 - 2 + 2 + 1 + 1
= 22 - 2 + 2 + 2
= 24
So the general formula with dilation is:
out = (D - 1) * stride - 2 * padding + dilation * (kernel_size - 1) + output_padding + 1
For our specific case with dilation=1:
out = (D - 1) * stride - 2 * padding + kernel_size + output_padding
OK so my original calculation was right. But I need to understand the mechanism.
In ConvTranspose3d, the conceptual operation is:
1. The input is dilated with (stride - 1) zeros between elements
2. The dilated input is padded with `padding` zeros on both sides
3. A regular conv is applied with the given kernel_size, stride=1, dilation=1
4. The output is padded with `output_padding` zeros on one side
So for our example:
1. Dilate: 12 elements become 23 elements (positions 0, 2, 4, ..., 22)
2. Pad with 1: 25 elements, data from index 1 to 23
3. Conv with k=3, stride=1: output of size 25 - 3 + 1 = 23
4. Output padding 1: final size 24
Hmm, but where exactly is the output padding applied? I think it's added to the "right" side.
Actually, another way to think about it: the formula ensures that the output size is approximately input_size * stride when padding is small.
Now for the kernel computation. Let me think about how to map from output position to input.
For output position od in [0, 23]:
This came from the conv at position od in the 23-element conv output.
The conv at position od uses the padded input at positions [od, od+1, od+2].
Wait, if conv output at position o is sum over k of padded_input[o + k] * weight[k], then:
- Output 0: uses padded_input[0, 1, 2]
- Output 22: uses padded_input[22, 23, 24]
But our padded input has indices 0 to 24 (25 elements), with actual dilated input data at indices 1 to 23.
So for output 0:
- padded_input[0] = 0 (padding)
- padded_input[1] = input[0] (dilated position 0 -> input[0])
- padded_input[2] = 0 (dilated position 2 -> input[1], but wait position 2 should be dilated input[1]*stride=2)
Wait, I think I made an error. Dilated input positions:
input[0] → position 1 in padded array (after adding padding of 1)
input[1] → position 3
input[2] → position 5
...
input[i] → position 1 + i*2 = 1 + 2i
So padded positions with non-zero data: 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23
For output position od, conv uses padded_input[od + k] for k in [0, 2].
- od=0: padded[0, 1, 2] → only padded[1] is non-zero, which is input[0]
- od=1: padded[1, 2, 3] → padded[1]=input[0], padded[3]=input[1]
- od=2: padded[2, 3, 4] → only padded[3]=input[1]
- ...
So for a given output od, the non-zero contributions come from kernel offsets k where:
(od + k) is in the dilated input range [1, 23] AND (od + k - 1) % 2 == 0
More generally:
p_d = od + k_d
For p_d to correspond to actual input at index id:
- padding_d <= p_d <= (D-1)*stride_d + padding_d
- (p_d - padding_d) % stride_d == 0
- id = (p_d - padding_d) / stride_d
And then the weight used is weight[ic, oc_local, k_d, k_h, k_w] for the appropriate group.
But wait, this is for the conceptual "regular conv on padded dilated input" view. However, ConvTranspose3d weight might be stored/flipped differently.
Actually, in ConvTranspose3d, PyTorch stores the weight as (in_channels, out_channels//groups, *kernel_size). When used in a regular conv on dilated input, I believe the weight is applied directly without flipping (since transposed conv is already the transpose operation).
Hmm, but actually I need to be more careful. Let me think about what ConvTranspose3d really computes.
The simplest mental model: ConvTranspose3d is the gradient computation of Conv3d. If Conv3d computes y = conv(x, w), then ConvTranspose3d applied to dy computes dx = conv_transpose(dy, w).
In terms of actual computation for forward pass of ConvTranspose3d:
Given input x and weight w, output y is computed such that if we were to backprop through the corresponding Conv3d, we'd get x.
For a simple 1D case with no padding, stride=1:
Conv3d output at position i: y[i] = sum_k x[i+k] * w[k]
Then dy/dx[i] = sum_k dy[i-k] * w[k]
So transposed conv with input x (playing the role of dy) gives:
output[i] = sum_k x[i-k] * w[k] = sum_k x[i+k] * w[-k] if we flip
Actually, this gets complicated. Let me just test with PyTorch to verify the exact behavior, but I can't run arbitrary code easily. Let me think about this differently.
Actually, I can use the reference implementation. If I write a custom kernel, I need to match the PyTorch output exactly. One approach is:
1. Implement the kernel using what I believe is the correct formula
2. Test it for correctness
3. If it doesn't match, debug
But for efficiency, let me look up what ConvTranspose3d actually does.
PyTorch's ConvTranspose3d is mathematically equivalent to:
- For each output position
- The value is a sum over all input positions that "see" this output position in the corresponding forward conv
For a transposed conv with stride s, kernel k, padding p:
An input at position i "spreads" to output positions [i*s - p, i*s - p + k - 1] (for stride=1 conv view of transposed conv).
Wait, I think it's better to think in terms of the standard formula for transposed convolution. Given input x, the output at position o is:
y[o] = sum_{i, k} x[i] * w[k] * indicator(o == i*s - p + k)
Where the indicator says that output o receives contribution from input i with kernel element k if o = i*s - p + k.
Or equivalently, for each output o:
y[o] = sum_k x[(o + p - k) / s] * w[k]
Where the division is only valid when (o + p - k) is divisible by s and the result is a valid input index.
Hmm wait, this is if the transposed conv corresponds to a regular conv with stride s and padding p. But actually the relationship is more subtle.
For Conv3d with stride=s, padding=p:
y[o] = sum_k x[o*s - p + k] * w[k] (simplified 1D)
Then transposed conv with the same parameters should satisfy:
If y = conv(x, w, stride=s, padding=p), then for any z, dot(conv_transpose(z, w), x) = dot(z, conv(x, w)).
This means conv_transpose at position i receives contributions from z[o] where o*s - p + k = i, i.e., o = (i + p - k) / s.
Hmm, but this assumes the conv uses full output. In practice, conv produces output of a specific size.
Let me just use the expansion formulation which I know is correct:
1. Create dilated input by inserting (stride-1) zeros between input elements
2. Pad with `padding` zeros on both sides
3. Convolve with stride=1 using the transposed conv weight
4. Add output_padding zeros to the output
The key question is: does the conv in step 3 use the weight as-is, or flipped?
For a regular conv, stride=1, no padding, output at o is sum_k input[o+k] * w[k].
For dilated input padded with `padding`, the valid output range is such that o+k is in [padding, dilated_size + padding - 1].
Actually, I think the simplest approach is to write a kernel that directly matches PyTorch's output by computing the output as I described and then testing it.
Let me write a kernel with the following strategy:
- Each thread computes one output element
- For each output element (b, oc, od, oh, ow):
- Determine group g = oc / (out_channels / groups)
- oc_local = oc % (out_channels / groups)
- Initialize accumulator to 0
- For each kernel offset (kd, kh, kw):
- Compute dilated position pd = od + kd, ph = oh + kh, pw = ow + kw
- Check if padded position is valid: padding_d <= pd <= (D-1)*stride_d + padding_d
- Check if pd corresponds to actual input: (pd - padding_d) % stride_d == 0
- If valid, compute input index id = (pd - padding_d) / stride_d
- Check if id is in valid range [0, D-1]
- Similarly for h and w
- If all valid, accumulate over input channels in the group:
- For ic_local in [0, in_channels/groups - 1]:
- ic = g * (in_channels/groups) + ic_local
- accumulator += input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
- Add bias if present
- Store result
Wait, but this formula for pd might be wrong. Let me verify with a simple example.
Actually, now I'm confused because there are different formulations. Let me think about this more carefully using known results.
For nn.ConvTranspose2d with input size H, kernel k, stride s, padding p, output_padding op:
Output size = (H - 1) * s - 2p + k + op
One common implementation internally converts ConvTranspose to a regular conv by:
1. Creating a dilated input of size (H-1)*s + 1
2. Padding with p on each side
3. Running a regular conv with kernel k, stride 1
4. The weight used is the same weight tensor
The output of step 3 has size:
(H-1)*s + 1 + 2p - k + 1 = (H-1)*s - 2p + k + 2 - 1 = (H-1)*s - 2p + k + 1
Then with output_padding = op, the final size is:
(H-1)*s - 2p + k + 1 + op - 1 = (H-1)*s - 2p + k + op
Wait, that doesn't add up. If conv on padded dilated input gives size A, and we add output_padding, we get size A + op. But the formula says:
out = (H-1)*s - 2p + k + op
And conv on padded dilated input of size (H-1)*s + 1 + 2p with kernel k gives:
(H-1)*s + 1 + 2p - k + 1 = (H-1)*s + 2p - k + 2
So: (H-1)*s + 2p - k + 2 + op = (H-1)*s - 2p + k + op
This would require: 2p - k + 2 = -2p + k, i.e., 4p = 2k - 2, which is not generally true.
I made an error somewhere. Let me look up the actual formula.
OK I just realized: in ConvTranspose, the `padding` parameter does NOT add padding to the dilated input. Instead, it effectively removes padding from the output.
Actually, according to PyTorch documentation:
"The padding argument effectively adds dilation * (kernel_size - 1) - padding amount of zero padding to both sizes of the input."
Wait, what? Let me read this more carefully:
For ConvTranspose, the effect of `padding`:
"padding controls the amount of implicit zero padding on both sides for dilation * (kernel_size - 1) - padding number of points."
Hmm, this is confusing. Let me look at the actual formula from PyTorch docs:
For ConvTranspose1d:
L_out = (L_in - 1) * stride - 2 * padding + dilation * (kernel_size - 1) + output_padding + 1
OK so the formula I have is correct. But the implementation details matter.
Another way to think about it: ConvTranspose with parameters (stride, padding) is equivalent to convolving a dilated input with stride=1 and padding=0, but with a modified effective kernel.
Actually, I found a clearer description. In ConvTranspose:
- The `padding` parameter in ConvTranspose actually corresponds to cropping/removing padding from the dilated input before convolution, rather than adding padding.
- Specifically, dilated input is first computed, then `padding` amount is removed from both ends, then conv is applied.
Let me test this theory:
Dilated input size = (H-1)*s + 1 = (12-1)*2 + 1 = 23
Remove padding from both ends: 23 - 2*1 = 21
Conv with k=3, stride=1: 21 - 3 + 1 = 19
Add output_padding=1: 20
But expected is 24. So that's wrong.
Another theory: the dilated input is padded with `padding` zeros, but then the conv uses a different effective formula.
Wait, let's work backwards from the answer. Expected output size is 24.
If output = dilated_size + 2*padding - kernel_size + 1 + output_padding = 24
Then: 23 + 2*p - 3 + 1 + 1 = 24
22 + 2p = 24
p = 1
So: dilated size 23, pad by 1 on each side → 25, conv with k=3 → 23, add output_padding 1 → 24. ✓
But then the formula says:
dilated_size = (H-1)*s + 1 = 23
padded_size = dilated_size + 2*padding = 25
conv_output = padded_size - k + 1 = 25 - 3 + 1 = 23
final_output = conv_output + output_padding = 24
And the expected formula: (H-1)*s - 2p + k + output_padding = 22 - 2 + 3 + 1 = 24 ✓
So my original intuition was correct! The formula matches. I just had an algebra error earlier.
conv_output = (H-1)*s + 1 + 2p - k + 1 = (H-1)*s + 2p - k + 2
Expected: (H-1)*s - 2p + k + op
These are not the same without additional conditions. Wait...
(H-1)*s + 2p - k + 2 + op =? (H-1)*s - 2p + k + op
Only equal if 4p - 2k + 2 = 0, i.e., 2p - k + 1 = 0.
But I computed with specific numbers:
(H-1)*s + 2p - k + 2 + op = 22 + 2 - 3 + 2 + 1 = 24 ✓
(H-1)*s - 2p + k + op = 22 - 2 + 3 + 1 = 24 ✓
24 = 24! For these specific numbers both give 24. Let me check:
22 + 2 - 3 + 2 + 1 = 24 ✓
22 - 2 + 3 + 1 = 24 ✓
But according to my formula: conv_output + op = dilated_size + 2p - k + 2 + op = 23 + 2 - 3 + 2 + 1 = 25
Wait, that's 25, not 24! Let me recheck.
dilated_size = 23
padded_size = 25
conv_output = padded_size - k + 1 = 25 - 3 + 1 = 23
final = 23 + op = 23 + 1 = 24
My "formula": dilated_size + 2p - k + 2 + op = 23 + 2 - 3 + 2 + 1 = 25 ≠ 24
Where did I get "dilated_size + 2p - k + 2 + op"? I said conv_output = padded_size - k + 1. So conv_output = (dilated_size + 2p) - k + 1. Then final = conv_output + op = dilated_size + 2p - k + 1 + op = 23 + 2 - 3 + 1 + 1 = 24 ✓
I mistakenly added an extra +1 earlier. OK so everything is consistent!
Good. So the formulation is:
1. Dilate input with (stride-1) zeros between elements
2. Pad with `padding` zeros on both sides
3. Regular convolution with kernel, stride=1
4. Add `output_padding` zeros (to one side, typically the "end")
Now, the critical question is: what is the weight indexing in step 3?
PyTorch ConvTranspose3d stores weight as (in_channels, out_channels//groups, *kernel_size). In the dilated-input-conv formulation, this weight is used directly in a regular conv, where:
- Each output channel is computed from all input channels in its group
- But wait, the weight shape for regular conv is (out_channels, in_channels//groups, *kernel_size), which is different!
Actually, in the dilated-input-conv formulation of transposed conv:
- The "input" is the dilated tensor (with in_channels feature maps)
- The "output" is the conv result (with out_channels feature maps)
- For regular conv, we'd need weight (out_channels, in_channels//groups, *kernel_size)
But PyTorch stores ConvTranspose weight as (in_channels, out_channels//groups, *kernel_size).
So the relationship is: if you transpose the first two dimensions of the ConvTranspose weight, you get something like the Conv weight, but with in_channels and out_channels swapped.
Actually, I think the correct statement is:
- In PyTorch's ConvTranspose, the operation is equivalent to taking the dilated input, and for each spatial position, computing a linear transformation from in_channels/group to out_channels/group
- The weight is applied as: for output channel oc in group g, output[b, oc, ...] += dilated_input[b, ic, ...] * weight[ic, oc_local, ...]
- This is different from regular conv where: output[b, oc, ...] += input[b, ic, ...] * weight[oc, ic_local, ...]
So in the transposed conv, each input channel contributes to all output channels in its group, using the same spatial kernel. This is exactly the structure of the stored weight.
To write the kernel, I need to compute:
output[b, oc, od, oh, ow] = sum over (ic, kd, kh, kw) of
input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
Where (id, ih, iw) is derived from (od, oh, ow, kd, kh, kw) based on the dilation/stride/padding.
Actually, I'm still not 100% sure about the spatial mapping. Let me think about it with a concrete example.
Consider 1D case: input=[a, b, c] (size 3), stride=2, padding=1, kernel_size=3.
Dilated input: [a, 0, b, 0, c] (positions 0, 2, 4)
Pad with 1: [0, a, 0, b, 0, c, 0] (positions 0 to 6)
Weight: [w0, w1, w2]
Regular conv with stride=1:
output[0] = 0*w0 + a*w1 + 0*w2 = a*w1
output[1] = a*w0 + 0*w1 + b*w2 = a*w0 + b*w2
output[2] = 0*w0 + b*w1 + 0*w2 = b*w1
output[3] = b*w0 + 0*w1 + c*w2 = b*w0 + c*w2
output[4] = 0*w0 + c*w1 + 0*w2 = c*w1
output[5] = c*w0 + 0*w1 + 0*w2 = c*w0
With output_padding=0, final output has 6 elements.
Output size = (3-1)*2 - 2*1 + 3 + 0 = 4 - 2 + 3 = 5
Wait, but I got 6 elements above. What's wrong?
Let me recheck: padded dilated input has size (3-1)*2 + 1 + 2*1 = 5 + 2 = 7.
Conv with k=3 gives 7 - 3 + 1 = 5 elements.
Oh I see, I miscounted. Dilated input: [a, 0, b, 0, c] has 5 elements.
Pad with 1 on each side: 7 elements, indices 0 to 6, with data at indices 1, 3, 5.
output[0]: padded[0,1,2] = [0, a, 0] → a*w1
output[1]: padded[1,2,3] = [a, 0, b] → a*w0 + b*w2
output[2]: padded[2,3,4] = [0, b, 0] → b*w1
output[3]: padded[3,4,5] = [b, 0, c] → b*w0 + c*w2
output[4]: padded[4,5,6] = [0, c, 0] → c*w1
That's 5 outputs. With output_padding=1, we'd get 6 outputs: [a*w1, a*w0+b*w2, b*w1, b*w0+c*w2, c*w1, 0]
Let me verify with formula: (3-1)*2 - 2*1 + 3 + 1 = 4 - 2 + 3 + 1 = 6. ✓
Great, so now I understand the mapping. For output position o:
output[o] = sum_k padded_dilated_input[o + k] * w[k]
Where padded_dilated_input[i] is non-zero only for i in [padding, padding + (D-1)*stride] with step stride.
For our kernel, given output (od, oh, ow) and kernel offset (kd, kh, kw):
The padded dilated input position is (od + kd, oh + kh, ow + kw).
This position contains actual input if:
1. padding_d <= od + kd <= padding_d + (D-1)*stride_d
2. (od + kd - padding_d) % stride_d == 0
And the actual input index is:
id = (od + kd - padding_d) / stride_d
Similarly for h and w.
So the formula I had earlier was correct! I just had some initial confusion.
Now, let me verify: for the 1D example:
stride=2, padding=1, D=3, k=3
For output o=1, kernel k=0:
padded position = 1 + 0 = 1
1 in [1, 1+4] = [1, 5] ✓
(1 - 1) % 2 = 0 ✓
id = 0/2 = 0
weight[0] * input[0] = a * w0
For output o=1, kernel k=2:
padded position = 1 + 2 = 3
3 in [1, 5] ✓
(3 - 1) % 2 = 0 ✓
id = 2/2 = 1
weight[2] * input[1] = b * w2
So output[1] = a*w0 + b*w2, which matches!
Great, so my formula is correct. Now let me think about the output_padding.
With output_padding=1, the output is:
[o0, o1, o2, o3, o4, 0]
So output_padding just adds zeros at the end? Or does it expand the conv output?
Actually, it depends on how you think about it. The conv produces 5 elements. To get 6 elements with output_padding=1, you either:
1. Append 1 zero, or
2. The conv actually produces 6 elements somehow
I think in practice, output_padding means that the "effective" dilated input is larger by output_padding, or equivalently, the conv output indices go up to conv_output_size + output_padding - 1, with the extra positions being zero.
For our kernel, a simpler formulation: the output size is fixed. For output positions that exceed what the conv on padded dilated input would produce, we just output zero.
Wait, but in our 1D example, the conv produces outputs 0 to 4. With output_padding=1, the output is 0 to 5. Is output[5] necessarily zero?
Let's compute output[5]:
padded[5, 6, 7] but padded only goes up to index 6.
5+2=7 > 6, so this is out of bounds.
If we just say: for any kernel offset k where o+k is outside [0, padded_size-1], it contributes 0.
Then output[5] with k=0: padded[5]=c, but 5+0=5 <= 6, ok
k=1: padded[6]=0, 5+1=6 <= 6, ok
k=2: padded[7]=? out of bounds, 5+2=7 > 6, so 0
output[5] = c*w1 + 0*w2 = c*w1, not 0!
But wait, the formula says output size is 6, and my manual computation says output[5] = c*w1. But is this what PyTorch actually does?
Let me think about this differently. The standard way to think about output_padding is that it adds extra size to the output. In the context of ConvTranspose being the "inverse" of Conv:
If Conv with stride=2, padding=1, kernel=3 maps size 6 → size 3, then ConvTranspose with the same parameters should map size 3 → size 6.
For input of size 6, conv output = floor((6 + 2*1 - 3) / 2) + 1 = floor(5/2) + 1 = 3. ✓
So the transposed conv maps 3 → 6. And output_padding=1 is needed to get exactly size 6 instead of size 5.
But wait, for input size 7: conv output = floor((7 + 2 - 3)/2) + 1 = floor(6/2) + 1 = 4.
And transposed conv with output_padding=0 maps 4 → ?
(4-1)*2 - 2 + 3 = 6 - 2 + 3 = 7.
So output_padding is specifically needed to exactly invert the downsampling.
In any case, for our kernel:
- Output indices go from 0 to D_out - 1
- For each output position, we just compute the valid contributions
- If there are no valid contributions, output is 0 (or bias)
So output_padding doesn't need special handling in the kernel! The formula for D_out already accounts for it, and we just compute output at each position by summing valid contributions. Some positions at the boundary might happen to sum to 0 (especially if they're entirely in the "output_padding" region), but that's fine.
Wait, actually there might be some difference. Let me think about whether output_padding affects which input positions are accessed.
Actually, in the formulation where we do conv on padded dilated input, output_padding simply means the conv output is larger. How?
If output_padding = op, then the "effective" kernel position or dilated input size is adjusted. One way to implement this is:
- The conv computes at positions 0 to (D-1)*stride - 2p + k - 1 + op
- For positions beyond what the standard conv would give, the result depends on how exactly the padding is applied
But I think the easiest way is: output_padding simply allows the conv to be computed over more positions, and the extra positions include more of the padded dilated input on the right side.
Actually, another standard view: for ConvTranspose with stride s, the dilated input is placed starting at position 0, and then the conv with kernel k is computed. The output size without output_padding would be (D-1)*s + k - 2p if p is the padding removed. With output_padding, it's (D-1)*s + k - 2p + op.
In implementation, I think the most reliable approach is:
1. Compute output using the formula: for each output (od, oh, ow), iterate over kernel offsets
2. Check if od + kd is in the range where dilated input is present (considering padding and output_padding)
3. The actual condition depends on how PyTorch implements it
Let me just look at what Aten/PyTorch actually does. Actually, PyTorch's ConvTranspose uses the im2col + gemm approach or calls into cuDNN. For our custom kernel, the exact behavior matching is essential.
You know what, let me just write a test kernel with my understanding and see if it matches PyTorch. If not, I can debug.
Actually, I realize there's another subtlety. In the conv formulation of transposed conv, does the weight get flipped or not?
When we do regular conv on dilated input:
output[o] = sum_k dilated_input[o + k] * weight[k]
In transposed conv, is this the correct formula, or is it:
output[o] = sum_k dilated_input[o + k] * weight[k]
Actually for regular conv: y[n, c_out, i] = sum_{c_in, k} x[n, c_in, i + k] * w[c_out, c_in, k]
For transposed conv (the transpose operation of conv):
If y = conv(x, w), then dy/dx is computed by a transposed conv applied to dy.
For conv with stride=1:
dy[n, c_in, i] / dx would involve: for all c_out, k where i = j + k for some output j:
dy/dx[n, c_in, i] = sum_{c_out, k} dy[n, c_out, i - k] * w[c_out, c_in, k]
So transposed conv applied to z gives:
out[n, c_in, i] = sum_{c_out, k} z[n, c_out, i - k] * w[c_out, c_in, k]
Now with groups, if w has shape (out_channels, in_channels/groups, k), then for group g:
out[n, c_in, i] = sum_{c_out in group g, k} z[n, c_out, i - k] * w[c_out, c_in_local, k]
If we rename c_in → c_out (since it's the output of transposed conv) and z[n, c_out, :] is the input to transposed conv:
out[n, c_out_of_transpose, i] = sum_{c_in in group g, k} z[n, c_in, i - k] * w[c_in, c_out_local, k]
Wait, the weight indexing is confusing. Let me use different notation.
For Conv3d with groups:
- weight shape: (out_channels, in_channels_per_group, k_d, k_h, k_w)
- y[n, oc, i] = sum_{ic in group(oc), k} x[n, ic, i+k] * w[oc, ic_local, k]
For ConvTranspose3d:
- weight shape: (in_channels, out_channels_per_group, k_d, k_h, k_w)
- This is the "adjoint" operation. If y = conv(x, w), then backprop to x is: dx[i] = sum_{oc, k} dy[oc, i-k] * w[oc, ic_local, k]
But PyTorch ConvTranspose implements this as a forward operation. So if z is the input to ConvTranspose:
out[n, oc, i] = sum_{ic in group(oc), k} z[n, ic, i - k] * w[ic, oc_local, k]
This still assumes stride=1. For stride > 1, the input is dilated, so:
out[n, oc, i] = sum_{ic in group(oc), k} dilated_z[n, ic, i - k] * w[ic, oc_local, k]
Or equivalently:
out[n, oc, i] = sum_{ic in group(oc), k} z[n, ic, (i - k - p) / s] * w[ic, oc_local, k] for valid indices
Actually this is getting messy. Let me use the padded dilated formulation:
dilated_z is zero except at positions i where (i - p) % s == 0.
out[n, oc, o] = sum_{ic in group(oc), k} dilated_z[n, ic, o + k] * w[ic, oc_local, k]
Wait, I keep going back and forth. Let me re-derive from the standard definition.
Actually, I realize there's yet another formulation. Many sources say ConvTranspose is equivalent to:
1. Distribute each input value to a k_d × k_h × k_w block in the output
2. The blocks overlap and sum
Specifically, for each input position (id, ih, iw) and each kernel offset (kd, kh, kw):
output[id*s + kd - p, ih*s + kh - p, iw*s + kw - p] += input[id, ih, iw] * weight[..., kd, kh, kw]
This is the "scatter" or "deconvolution" view.
For the output position formula to match, we need:
- output indices in range [0, D_out-1]
- For a given output (od, oh, ow), it receives contributions from inputs where:
id*s + kd - p_d = od
ih*s + kh - p_h = oh
iw*s + kw - p_w = ow
So: id = (od + p_d - kd) / s, which requires (od + p_d - kd) % s == 0.
Let me verify with the 1D example:
stride=2, padding=1, D=3, k=3
For output o=1:
id = (1 + 1 - kd) / 2 = (2 - kd) / 2
- kd=0: id = 2/2 = 1, input[1] * w[0] = b * w0
- kd=1: id = 1/2, not integer
- kd=2: id = 0/2 = 0, input[0] * w[2] = a * w2
So output[1] = b*w0 + a*w2
But earlier with the conv-on-dilated-input formulation, I got output[1] = a*w0 + b*w2.
These are different! In one case it's a*w0 + b*w2, in the other it's a*w2 + b*w0.
So the weight is flipped in one formulation vs the other!
In the "scatter" formulation where each input creates a block in the output:
input a at id=0 creates block: for kd=0,1,2: output[0*2 + kd - 1] += a * w[kd]
= output[-1] += a*w0, output[0] += a*w1, output[1] += a*w2
input b at id=1 creates block: output[1*2 + kd - 1] += b * w[kd]
= output[1] += b*w0, output[2] += b*w1, output[3] += b*w2
input c at id=2 creates block: output[2*2 + kd - 1] += c * w[kd]
= output[3] += c*w0, output[4] += c*w1, output[5] += c*w2
So:
output[0] = a*w1
output[1] = a*w2 + b*w0
output[2] = b*w1
output[3] = b*w2 + c*w0
output[4] = c*w1
output[5] = c*w2 (but wait, does output go up to 5?)
Expected size is 5 with output_padding=0. So output[5] would be out of bounds. Actually with output_padding=1, size is 6, so output[5] = c*w2.
Let me verify formula: size = (3-1)*2 - 2*1 + 3 + 1 = 6. But wait, that included op=1. With op=0: size = 5.
So op=0: outputs 0 to 4
op=1: outputs 0 to 5
With op=0 and only outputs 0-4:
output[4] = c*w1
And c*w2 goes to output[5] which doesn't exist.
But what about the conv-on-dilated-input formulation?
Dilated input: [a, 0, b, 0, c] padded: [0, a, 0, b, 0, c, 0]
Regular conv:
output[0] = [0,a,0]·[w0,w1,w2] = a*w1
output[1] = [a,0,b]·[w0,w1,w2] = a*w0 + b*w2
output[2] = [0,b,0]·[w0,w1,w2] = b*w1
output[3] = [b,0,c]·[w0,w1,w2] = b*w0 + c*w2
output[4] = [0,c,0]·[w0,w1,w2] = c*w1
output[5] = [c,0,0]·[w0,w1,w2] = c*w0
OK so with op=1, output[5] = c*w0 in the conv formulation, but c*w2 in the scatter formulation.
These are different! So which one does PyTorch use?
I think the conv-on-dilated-input formulation must be wrong about the weight indexing. Maybe the weight is flipped before being used in the conv.
Actually, yes! Transposed convolution is the transpose of regular convolution. If regular conv uses weight w[k] at offset k, then transposed conv uses w[k] but with the spatial relationship reversed.
Actually, let's think about this more carefully. Regular conv:
y[i] = sum_k x[i+k] * w[k] (using valid region)
Transposed conv (the adjoint):
For the adjoint, if <y, conv(x, w)> = <transpose_conv(y, w), x>
LHS = sum_i y[i] * sum_k x[i+k] * w[k] = sum_i sum_k y[i] * x[i+k] * w[k]
Let j = i+k, then i = j-k:
= sum_j sum_k y[j-k] * x[j] * w[k]
= sum_j x[j] * sum_k y[j-k] * w[k]
So transpose_conv(y, w)[j] = sum_k y[j-k] * w[k]
If we compare with scatter formulation above where input z at position id contributes z[id] * w[kd] to output[id*s + kd - p], then for output position j = id*s + kd - p:
z[id] * w[kd] where id = (j + p - kd) / s
So transpose_conv(y, w)[j] = sum_{k, (j+p-k)/s valid} y[(j+p-k)/s] * w[k]
With dilated input: let dilated_y[j'] = y[(j'-p)/s] when (j'-p) % s == 0
Then transpose_conv(y, w)[j] = sum_k dilated_y[j + k] * w[k]
Wait! This is:
transpose_conv(y, w)[j] = sum_k dilated_y[j + k] * w[k]
But earlier I said regular conv is y[i] = sum_k x[i+k] * w[k].
And transpose_conv(y, w)[j] = sum_k dilated_y[j + k] * w[k]
These have the same form! So transposed conv = regular conv on dilated input, with the SAME weight orientation.
But this conflicts with the scatter formulation! Let me check the scatter formulation again.
In scatter formulation: for each id,kd, output[id*s + kd - p] += input[id] * w[kd]
In the conv formulation for output j:
output[j] = sum_k dilated_input[j+k] * w[k]
dilated_input[j+k] is non-zero only when (j+k-p) % s == 0.
Let j+k = id*s + p, so id = (j+k-p)/s.
Then output[j] = sum_{k: (j+k-p)%s==0} input[(j+k-p)/s] * w[k]
In scatter formulation:
output[j] = sum_{kd: j = id*s + kd - p for some id} input[id] * w[kd]
= sum_{kd, id: j = id*s + kd - p} input[id] * w[kd]
= sum_{kd: (j+p-kd)/s is integer} input[(j+p-kd)/s] * w[kd]
Wait, id = (j + p - kd) / s in scatter, but id = (j + k - p) / s in conv.
These are id = (j + p - kd) / s vs id = (j + k - p) / s.
For them to be the same: p - kd = k - p, so 2p = k + kd.
But k and kd are the same variable (kernel offset). So id = (j + p - k) / s in scatter, and id = (j + k - p) / s in conv.
These are DIFFERENT! One has (p - k), the other has (k - p).
Ah, I see the issue. In the conv formulation, I used dilated_input[j+k] which assumes the conv starts at position j and looks forward k positions. But standard conv is y[i] = sum_k x[i+k] * w[k] where k goes from 0 to K-1.
For transposed conv, the correct formula is:
transpose_conv(y, w)[j] = sum_k y[j-k] * w[k]
So with dilated y, this becomes sum_k dilated_y[j-k] * w[k].
Let me recompute:
dilated_y[j-k] is non-zero when ((j-k) - p) % s == 0, i.e., j-k = id*s + p, so k = j - id*s - p.
Then transpose_conv(y, w)[j] = sum_{id: k=j-id*s-p in [0,K-1]} y[id] * w[j - id*s - p]
This matches the scatter formulation! output[j] = sum_{kd=j-id*s-p} input[id] * w[kd]
So the correct formula is:
output[j] = sum_k dilated_input[j-k] * w[k]
NOT output[j] = sum_k dilated_input[j+k] * w[k].
So I had the sign wrong! The transposed conv convolves with FLIPPED weight, or equivalently, it looks backward in the dilated input.
With dilated input padded with p:
output[j] = sum_k padded_dilated_input[j - k + p] * w[k] ... wait, I need to be more careful.
Actually, let me re-derive with padding included.
For regular conv with padding p:
y[i] = sum_k padded_input[i+k] * w[k]
where padded_input[j] = input[j-p] for j in [p, p+D-1], and 0 otherwise.
= sum_k input[i+k-p] * w[k] for valid indices.
Transpose of this:
For adjoint: sum_i y[i] * input[i+k-p] = sum_j transpose(y)[j] * input[j]
So transpose(y)[j] = sum_{k, i: i+k-p = j} y[i] * w[k]
= sum_k y[j+p-k] * w[k]
With dilated y: y[i] is non-zero when i = id*s, so:
transpose(y)[j] = sum_{k, id: id*s + p - k = j} y[id*s] * w[k]
= sum_k dilated_y[j + k - p] * w[k]
Hmm wait, let me use the padded dilated input formulation more carefully.
dilated_input of size (D-1)*s + 1, padded to size (D-1)*s + 1 + 2p.
Call this P[j] for j in [0, (D-1)*s + 2p].
Regular conv on P:
output[i] = sum_{k=0}^{K-1} P[i+k] * w[k]
This gives output of size ((D-1)*s + 1 + 2p) - K + 1 = (D-1)*s + 2p - K + 2.
With op=0, expected size is (D-1)*s - 2p + K - 1 + 1 = (D-1)*s - 2p + K.
So: (D-1)*s + 2p - K + 2 =? (D-1)*s - 2p + K
Only if 4p - 2K + 2 = 0.
But with our numbers: D=3, s=2, p=1, K=3:
(D-1)*s + 2p - K + 2 = 4 + 2 - 3 + 2 = 5
Expected: (3-1)*2 - 2*1 + 3 = 4 - 2 + 3 = 5
They match! I was using the wrong formula for expected size. Let me recheck general formula.
PyTorch formula: L_out = (L_in - 1) * stride - 2 * padding + dilation * (kernel_size - 1) + output_padding + 1
With dilation=1:
= (D-1)*s - 2p + (K-1) + op + 1
= (D-1)*s - 2p + K + op
My conv formula: (D-1)*s + 2p - K + 2
Set equal: (D-1)*s - 2p + K + op = (D-1)*s + 2p - K + 2
-2p + K + op = 2p - K + 2
2K + op - 2 = 4p
K + op/2 - 1 = 2p
For D=3, s=2, p=1, K=3, op=0: 3 + 0 - 1 = 2 = 2*1 ✓
For these same values, both formulas give 5.
For our actual test values: D=12, s=2, p=1, K=3, op=1:
Expected: 11*2 - 2*1 + 3 + 1 = 22 - 2 + 4 = 24
My conv formula: 11*2 + 2*1 - 3 + 2 = 22 + 2 - 3 + 2 = 23 ✗
So they don't match in general! The conv formula gives 23 but expected is 24.
Wait, for the test values: (D-1)*s + 2p - K + 2 + op = 22 + 2 - 3 + 2 + 1 = 24. ✓
Actually I forgot to add output_padding in my conv formula. Let me redo:
conv_output = dilated_size + 2p - K + 1 = (D-1)*s + 1 + 2p - K + 1 = (D-1)*s + 2p - K + 2
final = conv_output + op = (D-1)*s + 2p - K + 2 + op
Expected: (D-1)*s - 2p + K + op
For these to be equal: 2p - K + 2 = -2p + K, i.e., 4p = 2K - 2, i.e., 2p = K - 1.
For our test: K=3, so 2p = 2, p=1. ✓
For D=3 example: K=3, p=1, 2*1 = 3-1 = 2 ✓
But what if p=0, K=3?
Expected: (D-1)*s + 3 + op
Conv: (D-1)*s - 3 + 2 + op = (D-1)*s - 1 + op
These are NOT equal! (D-1)*s + 3 ≠ (D-1)*s - 1
Hmm, but maybe p=0, K=3, s=2, D=3:
Expected: 4 + 3 + op = 7 + op
My conv formula: 4 - 1 + op = 3 + op
Remainder: this only works when 2p = K - 1, which is not general.
So the "conv on dilated input" formulation must be wrong about what padding means in ConvTranspose.
According to PyTorch docs for ConvTranspose:
"The padding argument effectively adds dilation * (kernel_size - 1) - padding amount of zero padding to both sizes of the input."
Oh! This is the key. `padding` in ConvTranspose is NOT adding padding to the dilated input. Instead, it's a "negative padding" or cropping.
Specifically, the effective padding added is: (K - 1) - p.
So for p=1, K=3: effective padding = 3 - 1 - 1 = 1. Which matches what I computed.
For p=0, K=3: effective padding = 3 - 1 - 0 = 2.
This means:
1. Dilate input with stride
2. Pad with (K-1-p) zeros on both sides
3. Apply regular conv with stride=1, padding=0
4. Add output_padding
Wait, but let me verify with p=0, K=3, s=2, D=3, op=0:
Dilated: size 5
Pad with K-1-p = 2: size 5 + 4 = 9
Conv with K=3, stride=1: size 9 - 3 + 1 = 7
Expected: (3-1)*2 - 0 + 3 + 0 = 7 ✓
And with p=1, K=3, s=2, D=3, op=0:
Dilated: size 5
Pad with K-1-p = 1: size 5 + 2 = 7
Conv: 7 - 3 + 1 = 5
Expected: 4 - 2 + 3 = 5 ✓
Great! So the correct formulation is:
Effective padding = kernel_size - 1 - padding
1. Dilate input with (stride-1) zeros
2. Pad with effective_padding zeros on both sides
3. Apply conv with stride=1, padding=0, kernel_size
4. Add output_padding
Now for the actual kernel computation. Given output position o:
In the conv formulation with padding=0 (after step 2):
output[o] = sum_k padded_dilated_input[o + k] * w[k]
Where padded_dilated_input has size dilated_size + 2*effective_padding.
And padded_dilated_input[j] is non-zero when:
effective_padding <= j < effective_padding + dilated_size
AND (j - effective_padding) % stride == 0
With actual input index = (j - effective_padding) / stride
So for output o and kernel offset k:
j = o + k
Condition: effective_padding <= o + k < effective_padding + dilated_size
And (o + k - effective_padding) % stride == 0
Let me verify with p=1, K=3, s=2, D=3, op=0:
effective_padding = 3 - 1 - 1 = 1
dilated_size = 5
For o=0:
k=0: j=0, not >= 1, skip
k=1: j=1, >=1 and < 1+5=6, (1-1)%2=0, id=0, input[0]*w[1]
k=2: j=2, in range, (2-1)%2=1, skip
So output[0] = a*w1. ✓
For o=1:
k=0: j=1, id=0, input[0]*w[0]
k=1: j=2, not divisible
k=2: j=3, id=1, input[1]*w[2]
output[1] = a*w0 + b*w2
For o=2:
k=0: j=2, skip
k=1: j=3, id=1, input[1]*w[1]
k=2: j=4, id=2, input[2]*w[2]
output[2] = b*w1 + c*w2
Wait, but earlier with the scatter formulation I had output[2] = b*w1. What's going on?
Hmm, for o=2, k=2, j=4, id=(4-1)/2=1.5, not integer!
Wait, (4-1) = 3, 3%2 = 1, not 0. So skip.
output[2] = b*w1. ✓
For o=3:
k=0: j=3, id=1, b*w0
k=1: j=4, skip
k=2: j=5, id=2, c*w2
output[3] = b*w0 + c*w2
For o=4:
k=0: j=4, skip
k=1: j=5, id=2, c*w1
k=2: j=6, skip (6 = 1+5, not < 6)
output[4] = c*w1. ✓
Great! This matches the scatter formulation! Let me recheck scatter:
output[0] = a*w1 ✓
output[1] = a*w0 + b*w2... wait, I previously computed a*w2 + b*w0.
In scatter, input a creates: output[-1]=a*w0, output[0]=a*w1, output[1]=a*w2
input b creates: output[1]=b*w0, output[2]=b*w1, output[3]=b*w2
input c creates: output[3]=c*w0, output[4]=c*w1, output[5]=c*w2 (but 5 is out of bounds for op=0)
So:
output[0] = a*w1
output[1] = a*w2 + b*w0
output[2] = b*w1
output[3] = b*w2 + c*w0
output[4] = c*w1
But in conv formulation, I got:
output[1] = a*w0 + b*w2
These are different! a*w2 + b*w0 vs a*w0 + b*w2.
Unless w is stored differently in the two formulations. In the scatter formulation, w[kd] at position kd means "multiply input by w[kd] and place at output offset kd". In conv formulation, w[k] is used at the same position.
Actually, there's a sign difference. In scatter: output[id*s + kd - p] += input[id] * w[kd]
For id=0, kd=0: output[0-1] = output[-1] += a*w[0]
For id=0, kd=1: output[1-1] = output[0] += a*w[1]
For id=0, kd=2: output[2-1] = output[1] += a*w[2]
In conv formulation: output[o] += dilated_input[j] * w[k] where j = o+k
For o=1, j=1, k=0: output[1] += dilated_input[1]*w[0] = a*w[0]
For o=1, j=3, k=2: output[1] += dilated_input[3]*w[2] = b*w[2]
So conv gives a*w[0] + b*w[2] for o=1.
Scatter gives a*w[2] + b*w[0] for o=1.
These are different! In scatter, w[0] is at the "start" of the block, while in conv, w[0] is at the "start" of the conv window.
For transposed conv, which one is correct?
I think the conv formulation might have the weight in a different order. If we think of conv_transpose as the transpose of conv, and conv uses w[0], w[1], w[2] at offsets 0, 1, 2, then conv_transpose should produce outputs that correspond to the adjoint.
Actually, let me think about this: if we have a 1D conv with kernel [w0, w1, w2] and input [a, b, c]:
output[0] = a*w0 + b*w1 + c*w2 (assuming full correspondence for simplicity)
Then if we backprop with dy = [d0, d1, d2]:
da = d0*w0 + d1*w1 + d2*w2
db = d0*w1 + d1*w2 + d2*w3... wait, that's with different kernel sizes.
For stride=1, no padding, kernel size 3:
y[0] = x[0]*w0 + x[1]*w1 + x[2]*w2
y[1] = x[1]*w0 + x[2]*w1 + x[3]*w2
Backprop:
dx[0] = dy[0]*w0
dx[1] = dy[0]*w1 + dy[1]*w0
dx[2] = dy[0]*w2 + dy[1]*w1
So with input z = dy:
transpose_conv(z, w)[0] = z[0]*w0
transpose_conv(z, w)[1] = z[0]*w1 + z[1]*w0
transpose_conv(z, w)[2] = z[0]*w2 + z[1]*w1
This is the "scatter" view where lower indices of w are used for lower offsets!
So for transposed conv of z=[a, b] with w=[w0, w1, w2]:
out[0] = a*w0
out[1] = a*w1 + b*w0
out[2] = a*w2 + b*w1
out[3] = b*w2
Hmm, but our transposed conv example had padding and stride. Let me apply the same logic with stride=2, padding=1, D=3.
In the forward conv (corresponding to our transposed conv), we'd have input of some size, conv with stride=2, padding=1, kernel=3, producing output of size 3.
For input of size 5:
y[0] = x[0]*w0 + x[1]*w1 + x[2]*w2 (positions 0, 1, 2 with padding=1 means padded input has x[0] at pos 1)
Wait, for input of size 5, padding=1, padded size=7.
y[0] = padded[0]*w0 + padded[1]*w1 + padded[2]*w2 = 0*w0 + x[0]*w1 + x[1]*w2
y[1] = padded[2]*w0 + padded[3]*w1 + padded[4]*w2 = x[1]*w0 + x[2]*w1 + x[3]*w2
y[2] = padded[4]*w0 + padded[5]*w1 + padded[6]*w2 = x[3]*w0 + x[4]*w1 + 0*w2
Backprop to x:
dx[0] = dy[0]*w1
dx[1] = dy[0]*w2 + dy[1]*w0
dx[2] = dy[1]*w1
dx[3] = dy[1]*w2 + dy[2]*w0
dx[4] = dy[2]*w1
So for z = dy of size 3:
transpose_conv(z, w)[0] = z[0]*w1
transpose_conv(z, w)[1] = z[0]*w2 + z[1]*w0
transpose_conv(z, w)[2] = z[1]*w1
transpose_conv(z, w)[3] = z[1]*w2 + z[2]*w0
transpose_conv(z, w)[4] = z[2]*w1
This exactly matches the scatter formulation! output[1] = a*w2 + b*w0.
So the correct formula is the scatter one:
output[j] = sum_{k, id: j = id*s + k - p} input[id] * w[k]
Or equivalently, in terms of conv on dilated input, the weight must be spatially flipped!
If regular conv is y[i] = sum_k dilated_input[i+k] * w[k], then transposed conv is:
output[j] = sum_k dilated_input[j-k] * w[k]
Or with flipped weight:
output[j] = sum_k dilated_input[j+k] * w[K-1-k]
Let me verify with K=3: w flipped is [w2, w1, w0]
output[1] = dilated[1]*w2 + dilated[2]*w1 + dilated[3]*w0 = 0*w2 + a*w1 + 0*w0... wait, that's not right.
Hmm, with j=1:
k=0: dilated[1+0]=dilated[1]=a, weight_flipped[0]=w2, product a*w2
k=1: dilated[1+1]=dilated[2]=0, weight_flipped[1]=w1, product 0
k=2: dilated[1+2]=dilated[3]=b, weight_flipped[2]=w0, product b*w0
output[1] = a*w2 + b*w0. ✓
Great! So if I use conv on dilated input with FLIPPED weight, it matches the scatter formulation. But since we read from weight directly, it's easier to just use the scatter formulation directly.
Actually, for the custom kernel, I'll use:
output[b, oc, od, oh, ow] = sum over contributing (ic, kd, kh, kw):
input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
Where:
id = (od + padding_d - kd) / stride_d (must be integer and in range)
Wait no, from scatter: od = id * stride_d + kd - padding_d
So: id = (od + padding_d - kd) / stride_d
OR equivalently: od + padding_d - kd must be divisible by stride_d, and id = (od + padding_d - kd) / stride_d must be in [0, D-1].
Similarly:
ih = (oh + padding_h - kh) / stride_h
iw = (ow + padding_w - kw) / stride_w
Let me verify with our 1D example (s=2, p=1, D=3, K=3):
od=0:
kd=0: (0+1-0)/2 = 0.5, skip
kd=1: (0+1-1)/2 = 0, id=0, input[0]*w[1] = a*w1
kd=2: (0+1-2)/2 = -0.5, skip
output[0] = a*w1 ✓
od=1:
kd=0: (1+1-0)/2 = 1, id=1, input[1]*w[0] = b*w0
kd=1: (1+1-1)/2 = 0.5, skip
kd=2: (1+1-2)/2 = 0, id=0, input[0]*w[2] = a*w2
output[1] = b*w0 + a*w2 ✓
Wait, this is a*w2 + b*w0. But in my conv backprop I got a*w2 + b*w0 too. ✓
Actually, I said scatter gave a*w2 + b*w0, which matches. In the conv-backprop derivation:
dx[1] = dy[0]*w2 + dy[1]*w0 = a*w2 + b*w0. ✓
Great!
Now, what about output_padding? The formula says out_size with op=1 is 6 instead of 5. How does op affect the computation?
In the backprop view, if the forward conv produced 3 elements with op=0, and we want 5 vs 6 or something... Actually, output_padding doesn't change the "core" computation. It just means the output tensor is larger, and the extra elements are computed if there are valid contributions, or set to zero if there are none.
Actually no, I need to think about this more carefully. Output_padding pads the output shape by adding extra "space" at the end. In the scatter view, it means some output positions that would not normally exist are now accessible.
For our 1D example with op=1, output size is 6.
od=5:
kd=0: (5+1-0)/2 = 3, id=3, out of range (D=3, so id in [0,2])
kd=1: (5+1-1)/2 = 2.5, skip
kd=2: (5+1-2)/2 = 2, id=2, input[2]*w[2] = c*w2
output[5] = c*w2
But in the conv-backprop view, the forward conv with input size 6 produces:
(6 + 2*1 - 3) / 2 + 1 = (7)/2 + 1 = 3 + 1 = 4 outputs
So if the transposed conv is the exact inverse, z should have size 4, and transposed conv gives size 6. But our test has z size 3 and gives size 5/6. So output_padding=1 just means the result is 1 larger than the "natural" inverse.
In any case, for the kernel:
- We compute output according to the scatter formula
- For output positions that have no valid contributions, the value is 0 (or bias)
Now I need to be careful about output_padding. In PyTorch, output_padding adds zeros (or more precisely, extends the output shape, with the extension filled by whatever the conv computes). Does it actually affect which (id, kd) pairs are valid?
Actually, yes! Because with output_padding=1, the output size is 6 instead of 5. For od=5, we computed c*w2 above. So it's not just "add zeros at the end". The conv computation is extended.
How does this happen in the "effective padding" formulation?
effective_padding = K - 1 - p = 3 - 1 - 1 = 1
dilated_size = 5
padded_dilated_size = 5 + 2*1 = 7
conv with K=3, op=1: output = 7 - 3 + 1 + 1 = 6
Wait, how does adding output_padding change the conv? I think it's effectively adding more padding to the dilated input, or allowing the conv to extend further.
If the padded dilated input has size 7, and we conv with K=3:
Valid conv outputs: 0 to 4 (5 outputs)
With op=1: output 0 to 5 (6 outputs)
For output 5, we'd need padded_input[5, 6, 7], but padded only goes up to index 6. So padded[7] = 0.
Actually, I think the typical implementation is that output_padding adds extra positions at the end of the output. For those extra positions, the conv is computed with the assumption that there's additional padding in the input.
So if op=1, it's like the padded dilated input has size 7 + 1 = 8 instead of 7. With padded positions 0 to 7, where original data is at 1, 3, 5.
For output 5: conv uses padded[5, 6, 7] = [0, c, 0] (wait, padded[5] should be 0 since data is at 1, 3, 5; padded[6] = 0 too)
Hmm, actually in the effective padding formulation with p=1, K=3:
effective_padding = K - 1 - p = 1
Dilated positions: 0, 2, 4 in the original dilated array.
After padding with 1: positions 0 to 6, with data at 1, 3, 5.
Wait, dilated of [a, b, c] with s=2: [a, 0, b, 0, c]
Pad by 1: [0, a, 0, b, 0, c, 0]
Indices: 0, 1, 2, 3, 4, 5, 6
Conv output 0: [0, a, 0] · [w0, w1, w2] = a*w1
Conv output 1: [a, 0, b] · [w0, w1, w2] = a*w0 + b*w2
Conv output 2: [0, b, 0] · [w0, w1, w2] = b*w1
Conv output 3: [b, 0, c] · [w0, w1, w2] = b*w0 + c*w2
Conv output 4: [0, c, 0] · [w0, w1, w2] = c*w1
Conv output 5: [c, 0, 0] · [w0, w1, w2] = c*w0
With op=1, output size is 6, so outputs 0 to 5. output[5] = c*w0.
But in the scatter formulation, I computed output[5] = c*w2.
These are DIFFERENT! c*w0 vs c*w2.
This is because of the weight flipping difference. In the effective-padding formulation, the conv uses weight [w0, w1, w2] directly. In the scatter formulation, the weight is used as-is. But we established earlier that scatter gives the correct result matching backprop.
Wait, let me re-examine. In the effective-padding formulation, I said conv output[5] uses [c, 0, 0]. But is this formulation actually correct for PyTorch?
Let me check: if PyTorch uses effective_padding = K - 1 - p, and then does a regular conv, then the weight would be used in its natural order. But we know from backprop analysis that transposed conv uses the weight "as-is" in the scatter sense, which gives a*w2 + b*w0 for output[1].
But in the effective-padding conv formulation, output[1] = [a, 0, b] · [w0, w1, w2] = a*w0 + b*w2.
These are different! So either:
1. The effective-padding formulation is not what PyTorch does, OR
2. The effective-padding formulation uses flipped weight
If PyTorch's ConvTranspose internally calls a conv routine, it might flip the weight first. Or it might use a completely different algorithm.
Let me check PyTorch source code mentally. PyTorch ConvTranspose uses the backward pass of Conv, which in cuDNN corresponds to the "backward data" algorithm. In cuDNN, the weight is NOT flipped; instead, the algorithm internally handles the transposed nature.
For our custom kernel, the safest approach is to compute exactly what PyTorch computes. Since I can't easily run PyTorch to inspect the exact output, I should use the scatter formulation which I derived from first principles as the mathematical transpose of convolution.
Actually, I CAN test this because the reference implementation is provided. I just need to write a kernel and see if it matches.
Let me write a simple but correct CUDA kernel using the scatter formulation.
Actually, on further reflection, I think there might be yet another subtlety. When we think about ConvTranspose3d in PyTorch, the weight tensor has shape (in_channels, out_channels/groups, k_d, k_h, k_w).
In the scatter formulation:
For each input channel ic, we scatter into output channels in the same group using weight[ic, :, kd, kh, kw].
In the conv-on-dilated-input formulation, if we view it as conv, the weight mapping would be different.
Let me just look at what PyTorch actually computes for a tiny example. But I can't because I don't have interactive PyTorch execution here. I'll just write a test kernel.
One important thing: PyTorch's ConvTranspose3d weight orientation. As discussed, weight is (in_channels, out_channels_per_group, *kernel_size). And for output channel oc in group g, the computation is:
output[b, oc, ...] = sum_{ic in group g} input[b, ic, ...] convolved with weight[ic, oc_local, ...]
Where "convolved with" means the transposed convolution operation.
So in the scatter formulation:
For each input position (id, ih, iw), each kernel offset (kd, kh, kw), and each output channel oc in the same group:
output[b, oc, id*s + kd - p, ih*s + kh - p, iw*s + kw - p] +=
input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
With the understanding that out-of-bound writes are ignored, and for (kd, kh, kw) where the output position is not an integer or out of bounds, nothing is written.
Actually wait, id*s + kd - p must be the output index. But s and p are per-dimension. So:
od = id * stride_d + kd - padding_d
oh = ih * stride_h + kh - padding_h
ow = iw * stride_w + kw - padding_w
And od must be in [0, D_out-1], etc.
For a reverse formulation (given od, find contributing id and kd):
od = id * s + kd - p
kd = od + p - id * s
Since kd ∈ [0, K-1], we have:
0 ≤ od + p - id*s ≤ K-1
id*s ≤ od + p ≤ id*s + K - 1
id ≤ (od + p) / s
id ≥ (od + p - K + 1) / s
So for a given od, id can range from ceil((od + p - K + 1) / s) to floor((od + p) / s).
And for each valid id, kd = od + p - id*s, which automatically satisfies 0 ≤ kd ≤ K-1.
This formulation is very clean! For each output position, we iterate over valid input positions id and compute the corresponding kd.
Similarly for ih, iw.
Then for each valid (id, ih, iw), we iterate over oc and ic in the group and accumulate.
This is very convenient for the kernel. Let me verify with the 1D example (s=2, p=1, K=3):
For od=0:
id ∈ [ceil((0+1-3+1)/2), floor((0+1)/2)] = [ceil(-0.5), floor(0.5)] = [0, 0]
id=0: kd = 0 + 1 - 0 = 1
output[0] += input[0] * w[1] = a*w1 ✓
For od=1:
id ∈ [ceil((1+1-3+1)/2), floor((1+1)/2)] = [ceil(-0.5), 1] = [0, 1]
id=0: kd = 1 + 1 - 0 = 2
output[1] += input[0] * w[2] = a*w2
id=1: kd = 1 + 1 - 2 = 0
output[1] += input[1] * w[0] = b*w0
output[1] = a*w2 + b*w0 ✓
For od=2:
id ∈ [ceil((2+1-3+1)/2), floor((2+1)/2)] = [ceil(0.5), 1] = [1, 1]
id=1: kd = 2 + 1 - 2 = 1
output[2] += input[1] * w[1] = b*w1 ✓
For od=3:
id ∈ [ceil((3+1-3+1)/2), floor((3+1)/2)] = [ceil(0.5), 2] = [1, 2]
id=1: kd = 3 + 1 - 2 = 2
output[3] += input[1] * w[2] = b*w2
id=2: kd = 3 + 1 - 4 = 0
output[3] += input[2] * w[0] = c*w0
output[3] = b*w2 + c*w0 ✓
For od=4:
id ∈ [ceil((4+1-3+1)/2), floor((4+1)/2)] = [ceil(1.5), 2] = [2, 2]
id=2: kd = 4 + 1 - 4 = 1
output[4] += input[2] * w[1] = c*w1 ✓
For op=1, od=5:
id ∈ [ceil((5+1-3+1)/2), floor((5+1)/2)] = [ceil(2), 3] = [2, 3]
id=2: kd = 5 + 1 - 4 = 2, output[5] += c*w2
id=3: out of range (D=3, id must be < 3)
So output[5] = c*w2. But in the effective-padding conv formulation, I got c*w0.
So which one does PyTorch produce?
Actually, I just realized: the "effective-padding conv formulation" I used earlier might have had the wrong weight indexing. If the weight [w0, w1, w2] in the conv formulation corresponds to [w2, w1, w0] in the scatter formulation, then they could match.
In PyTorch's ConvTranspose, does the weight get flipped spatially?
OK here is what actually happens in PyTorch internally. If you look at ConvTranspose2d in PyTorch, the operation is typically implemented as:
1. Flip the weight spatially (both H and W dimensions)
2. Use this flipped weight in a regular convolution on the dilated input
So PyTorch internally does: weight_flipped[k] = weight[K-1-k].
Let me verify: if PyTorch flips weight before conv, then in the conv formulation:
output[5] = [c, 0, 0] · [w2, w1, w0] = c*w2. ✓ matches scatter!
output[1] = [a, 0, b] · [w2, w1, w0] = a*w2 + b*w0. ✓ matches scatter!
Great! So both formulations are equivalent IF you flip the weight in the conv formulation. But PyTorch stores the weight in "unflipped" form, and the operation inherently uses it without flipping (from the user's perspective).
From the kernel writer's perspective, the scatter formulation is easier and doesn't require thinking about flipping. In scatter:
output[od] += input[id] * weight[kd] when od = id*s + kd - p
So for a custom kernel, I can just implement the scatter formulation directly.
Now let me think about output_padding again. With op=1:
For od=5, id=2, kd=2: output[5] += input[2] * w[2] = c*w2.
What if op were 2? Then output size = 7, and od=6:
id ∈ [ceil((6+1-3+1)/2), floor((6+1)/2)] = [ceil(2.5), 3] = [3, 3]
id=3: out of range, so output[6] = 0.
So output_padding=2 would add a zero at the end. PyTorch's output_padding parameter ensures that the transposed conv can exactly invert a corresponding conv, but only for certain combinations.
OK now I have a solid understanding. Let me design the kernel.
For our specific test:
- batch_size = 8
- in_channels = 32
- out_channels = 32
- kernel_size = (3, 5, 7)
- depth = 12
- height = 24
- width = 48
- stride = (2, 2, 2)
- padding = (1, 2, 3)
- output_padding = (1, 1, 1)
- groups = 4
Output sizes:
D_out = (12-1)*2 - 2*1 + 3 + 1 = 24
H_out = (24-1)*2 - 2*2 + 5 + 1 = 48
W_out = (48-1)*2 - 2*3 + 7 + 1 = 96
So output shape: (8, 32, 24, 48, 96)
Each group: 8 in_channels, 8 out_channels
Kernel: (3, 5, 7)
For each output position (od, oh, ow), valid inputs:
For od=0: id ∈ [ceil((0+1-3+1)/2), floor((0+1)/2)] = [0, 0], so id=0, kd=1
For od=1: id ∈ [0, 1], id=0(kd=2), id=1(kd=0)
For od=2: id ∈ [1, 1], id=1, kd=1
For od=3: id ∈ [1, 2], id=1(kd=2), id=2(kd=0)
etc.
Actually, for each od, there can be at most ceil(K/s) valid id's. Since K=3, s=2, at most 2 id's are valid for each od.
Similarly for h: K=5, s=2, at most 3 valid ih's per oh.
For w: K=7, s=2, at most 4 valid iw's per ow.
So for each output pixel, at most 2*3*4 = 24 input pixels contribute.
And with 8 input channels per group and 8 output channels per group, each output pixel requires at most 24*8 = 192 multiplies (before summing over the 8 output channels... wait, no).
Actually, for each output CHANNEL oc, we sum over all input channels ic in the group and all valid spatial positions:
output[b, oc, od, oh, ow] = sum_{ic in group} sum_{valid (id,ih,iw,kd,kh,kw)} input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
For each output element (b, oc, od, oh, ow), the computation involves:
- At most 24 spatial combinations
- For each spatial combo, 8 input channels
- Total: 24*8 = 192 multiply-adds
Wait, I think I miscounted. For each valid (id, kd) pair, there's exactly one kd given id and od. So the total number of (id, kd) pairs is at most 2 per od. Similarly (ih, kh) at most 3, (iw, kw) at most 4. So spatial combinations = 2*3*4 = 24.
For each spatial combination (id, ih, iw, and corresponding kd, kh, kw), we iterate over ic from 0 to 7:
sum += input[b, ic_g, id, ih, iw] * weight[ic_g, oc_local, kd, kh, kw]
So yes, 24*8 = 192 multiplies per output element.
Total output elements: 8*32*24*48*96 ≈ 28.3M
Total multiply-adds: 28.3M * 192 ≈ 5.4 GFLOPs
Wait, I calculated this wrong earlier. Let me recalculate:
The reference would compute: for each output element, sum over in_channels/groups * kernel_size elements.
With the naive formulation: in_channels/groups = 8, kernel_size = 3*5*7 = 105.
So 8*105 = 840 multiply-adds per output element.
But with the optimized formulation that skips the zeros in dilated input, we only do 192 multiply-adds per output element!
So the PyTorch reference using cuDNN might actually do 840 multiply-adds per output element if it uses an im2col + gemm approach that materializes the dilated input. Or cuDNN might be smart and skip zeros.
Actually, cuDNN's ConvTranspose backward_data algorithm is highly optimized and likely doesn't materialize all the zeros. But it's possible our custom kernel can still win with better data layout and shared memory usage for this specific configuration.
Wait, let me reconsider. cuDNN backward_data for ConvTranspose might actually use winograd or FFT for some cases, but probably not for 3D conv with small kernels. It likely uses some form of implicit GEMM that's quite efficient.
However, a custom kernel has the advantage of:
1. No overhead from cuDNN heuristics and algorithm selection
2. Can be perfectly tuned for the exact sizes
3. Can exploit the specific sparsity pattern explicitly
For example, since only 24 out of 105 kernel positions are "active" per output position, we can skip 77% of the multiplications!
Actually, is this right? Let me think again.
In PyTorch's ConvTranspose3d operation, the effective computation (mathematically) involves all kernel positions, but many of them multiply by zero because of the dilation. However, cuDNN doesn't know about the zeros in advance unless it's specifically optimized for transposed conv.
For transposed conv, cuDNN typically uses the backward_data algorithm, which might be implemented as:
1. The input is treated as a "dy" tensor
2. The output is "dx"
3. The algorithm computes dx = conv_transpose(dy, w)
In this case, cuDNN likely uses an implicit GEMM approach where each output element is computed by iterating over all kernel positions and checking some condition. Or it might use a direct approach similar to what we'd write.
Actually, for 3D convolutions, cuDNN's algorithms might not be as heavily optimized as 2D. And with the specific sizes here (small batch, groups=4), there might be room for improvement.
Let me design a kernel strategy:
Option 1: Output-stationary, one thread per output element
- For each output element, iterate over all valid (id, ih, iw, kd, kh, kw, ic)
- Pro: simple, good parallelism
- Con: each thread does 24*8 = 192 loads from input and weight, which might not be coalesced
Option 2: Tile-based with shared memory
- Load a tile of input into shared memory
- Load a tile of weight into shared memory
- Compute a tile of output
- Pro: better data reuse
- Con: more complex, shared memory might be limiting for 3D
For our test sizes:
- Input: (8, 32, 12, 24, 48) ≈ 14.1 MB
- Output: (8, 32, 24, 48, 96) ≈ 113.2 MB
- Weight: (32, 8, 3, 5, 7) ≈ 107.5 KB
Each output element is touched once. Each input element is read by multiple output elements (roughly (2*3*4)*8 = 192 times in the naive formulation, but actually each output element only reads 24 input positions, and each input position contributes to about 24*8 = 192 output channels... wait, no.
Actually, each input spatial position (id, ih, iw) contributes to multiple output spatial positions. With K=(3,5,7) and s=(2,2,2), each input spatial position contributes to at most 3*5*7 = 105 output positions, but because of stride, only certain combinations are valid. Actually in the scatter formulation, each input at (id, ih, iw) with kernel (kd, kh, kw) contributes to exactly one output position:
od = id*s_d + kd - p_d = id*2 + kd - 1
So each input element is used in exactly K_d*K_h*K_w = 105 output positions per output channel? No, wait.
For each input element at (id, ih, iw), and for each (kd, kh, kw), there's one output spatial position. So each input element contributes to 3*5*7 = 105 output spatial positions.
But since we have 8 output channels per group, and each input channel contributes to all 8 output channels in its group, each input element is read:
- 105 times for spatial positions
- For each spatial position, 8 times for the 8 output channels
Wait, that's not right either. Let me think more carefully.
For each group, we have 8 input channels and 8 output channels.
The weight is shape (8, 8, 3, 5, 7).
For input channel ic and output channel oc, the kernel is weight[ic, oc, :, :, :].
So for a fixed input element at (b, ic, id, ih, iw):
It contributes to output channels oc in the same group, at spatial positions:
(b, oc, id*2 + kd - 1, ih*2 + kh - 2, iw*2 + kw - 3) for all kd, kh, kw.
For each (kd, kh, kw), there's one output spatial position and 8 output channels.
So each input element is read once and then multiplied by 8*105 = 840 weights to produce 8*105 output values? No, the weight for each (ic, oc, kd, kh, kw) is different.
For a fixed input element, to compute all its contributions, we'd need to do 8 (output channels) * 105 (kernel positions) multiplications.
Actually, total FLOPs can be computed as:
For each input element (b, ic, id, ih, iw):
For each oc in same group:
For each (kd, kh, kw):
output[b, oc, od, oh, ow] += input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
Number of input elements: 8 * 32 * 12 * 24 * 48 = 3,538,944
For each: 8 * 105 = 840 multiplies
Total: 3,538,944 * 840 ≈ 2.97 GFLOPs
Wait, that seems low. Let me recalculate.
For each output element, we do 8*105 = 840 multiply-adds in the naive formulation.
Total output elements: 8*32*24*48*96 = 28,311,552
Total multiply-adds: 28,311,552 * 840 ≈ 23.8 GFLOPs
But in the optimized formulation where we skip zeros, each output element only does 8*24 = 192 multiply-adds.
Total: 28,311,552 * 192 ≈ 5.4 GFLOPs
That's a huge difference! But cuDNN might already skip the zeros implicitly. Let me think...
Actually, in cuDNN's backward_data (which is used for ConvTranspose), the algorithm knows that the input has stride > 1 in the corresponding forward conv. The backward_data algorithm doesn't need to materialize zeros; it can directly compute the gradient. So cuDNN likely already does something efficient.
But can we beat cuDNN for this specific case?
One potential issue with cuDNN for 3D convolutions with groups=4 and relatively small sizes is that the problem might not be large enough to fully utilize the GPU, and cuDNN's generic algorithms might have some overhead.
A custom kernel that explicitly iterates only over valid (id, ih, iw, kd, kh, kw) combinations could be faster because:
1. It avoids any unnecessary computation
2. It can be perfectly coalesced for our specific sizes
3. It avoids cuDNN's algorithm selection/setup overhead
Let me try to write a kernel. I'll use a tile-based approach where:
- Each block computes a tile of output elements
- Threads within a block collaborate to load input and weight data
- Compute the partial sums
Given the sizes, let me think about a good tiling strategy.
One simple approach: each thread computes one output element. This gives us 28M threads, which should be enough to saturate the GPU. With 1024 threads per block, we have ~27K thread blocks.
For an H100 with 132 SMs (I think), each SM can run multiple blocks. With good occupancy, we should be fine.
The key question is memory bandwidth and cache behavior. If each thread loads 192 input values and 192 weight values, and these loads are not cached well, we might be memory bound.
Actually, for 28M output elements, each doing ~192 loads of input and 192 loads of weight:
Total input loads: 28M * 192 * 4 bytes ≈ 21.4 GB
Total weight loads: 28M * 192 * 4 bytes ≈ 21.4 GB
Total output writes: 28M * 4 bytes ≈ 113 MB
If these all go to global memory, total traffic ≈ 43 GB. At 3.35 TB/s, that's ~13 ms. But the loads should be cached.
Actually, many threads will load the same input and weight values. For example, threads computing nearby output positions will need overlapping input regions. With the L2 cache on H100 (50 MB), we might get good hit rates.
But a better approach is to use shared memory explicitly. Let's think about this.
If we process output channels in tiles, we could load a subset of weights into shared memory and reuse them. Similarly for input spatial positions.
Actually, maybe a better approach is to think in terms of the scatter operation:
For each input element, we "scatter" its value to multiple output positions. This is more memory-write-friendly.
Consider a kernel where each thread processes one input element and updates all its contributing output positions:
For input (b, ic, id, ih, iw):
For oc in [0, out_channels/group-1]:
For kd, kh, kw:
od = id*2 + kd - 1
oh = ih*2 + kh - 2
ow = iw*2 + kw - 3
if valid:
atomicAdd(&output[b, oc_g, od, oh, ow], input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw])
But atomic adds to the same output position from different threads would be a bottleneck.
Another approach: process output channels in tiles, and for each tile, iterate over input channels.
Let me think about a more structured approach:
One block could compute a (B, OC, OD, OH, OW) tile of the output. For example:
- 4 batches
- 8 output channels (one full group's worth)
- 4 depth positions
- 4 height positions
- 4 width positions
Total: 4*8*4*4*4 = 2048 output elements. With 256 threads, each thread computes 8 output elements.
Within the block:
1. Load the relevant input region into shared memory
2. Load the relevant weights into shared memory
3. Compute
For 8 output channels in a group, we need all 8 input channels in that group.
For spatial tile (4, 4, 4) in output, what input region do we need?
Given od ∈ [OD_start, OD_start+3], and od = id*2 + kd - 1 with kd ∈ [0,2]:
id can range from floor((OD_start + 3 + 1)/2) down to ceil((OD_start + 0 + 1 - 2)/2)
Roughly, id covers about (4/2) + 1 = 3 input depth positions.
Similarly for h: about 3 input height positions.
For w: about 3 input width positions.
So input region: 8 channels * 3 * 3 * 3 ≈ 216 floats → ~864 bytes
Weight region: 8 input channels * 8 output channels * 3*5*7 = 64 * 105 = 6720 floats → ~27 KB
Total shared memory: ~28 KB, which fits in the 227 KB limit.
But the weight region is for ALL spatial kernel positions. If we want to load weights for only the "active" combinations, it gets more complex because different (od, oh, ow) positions within the tile need different subsets of kernel positions.
Actually, maybe it's simpler to load the whole weight tensor for the group. 27 KB is manageable.
Alternatively, we can have each thread compute one output channel within the spatial tile, and stream over input channels.
Let me design a specific kernel:
For each group g:
Input channels: [g*8 : (g+1)*8]
Output channels: [g*8 : (g+1)*8]
Weight: [g*8 : (g+1)*8, :, :, :, :] which is [8, 8, 3, 5, 7]
Grid dimensions:
- x: batch * groups (8*4 = 32)
- y: output depth * output height * output width / tile_size (24*48*96 / 256)
- z: or some other decomposition
Actually, for simplicity let me try a 1D output element per thread approach first, but with the optimized iteration over valid spatial positions. If that doesn't beat cuDNN, I'll try tiling.
Actually, wait. For output-stationary with one thread per output element:
Each thread needs to:
1. Determine its (b, oc, od, oh, ow)
2. Find g = oc / 8, oc_local = oc % 8
3. For id in valid range for od:
kd = od + p_d - id*s_d
For ih in valid range for oh:
kh = oh + p_h - ih*s_h
For iw in valid range for ow:
kw = ow + p_w - iw*s_w
For ic_local in [0, 7]:
ic = g*8 + ic_local
sum += input[b, ic, id, ih, iw] * weight[ic, oc_local, kd, kh, kw]
This is straightforward. The number of iterations is at most 2*3*4*8 = 192.
The memory access pattern:
- input[b, ic, id, ih, iw]: For a fixed (b, oc, od, oh, ow), as we vary id, ih, iw, ic, the input addresses might not be contiguous.
- weight[ic, oc_local, kd, kh, kw]: For fixed oc_local and varying ic, the weight addresses are strided by the output channel dimension. With 8 output channels, weight[ic, :, kd, kh, kw] are contiguous in ic? Actually weight is (32, 8, 3, 5, 7). So weight[ic, oc_local, kd, kh, kw] for varying ic has stride 8*3*5*7 = 840. Not great for coalescing.
Wait, but each thread is independent. Within a warp, different threads have different (b, oc, od, oh, ow). Do they access the same weight? No, they have different oc_local and different (kd, kh, kw), so weight accesses are likely uncoalesced.
Input accesses might have some coalescing if threads in a warp have similar spatial positions.
For a warp of 32 threads, if they're computing consecutive output elements in some dimension, their (od, oh, ow) might differ, leading to different (id, ih, iw), so input accesses are also uncoalesced.
This lack of coalescing could hurt performance.
To improve coalescing, maybe I should:
1. Reorder the loops so that threads in a warp access contiguous memory
2. Or use shared memory to stage data
One idea: have each warp compute 32 output elements for the SAME spatial position but DIFFERENT output channels. Since output channels in a warp are contiguous, and weight[ic, oc_local, ...] for contiguous oc_local might have some pattern...
Actually, weight shape is (32, 8, 3, 5, 7). For a fixed ic and varying oc_local from 0 to 7, weight[ic, oc_local, kd, kh, kw] are contiguous in memory! Because the second dimension (8) is the innermost dimension before the spatial dimensions.
Wait, weight layout in PyTorch is (in_channels, out_channels_per_group, kd, kh, kw).
So weight[ic, 0, 0, 0, 0] and weight[ic, 1, 0, 0, 0] are adjacent in memory. Great!
So if a warp computes the same spatial position for output channels 0-7 (which are in different groups... wait, out_channels=32, groups=4, so channels 0-7 are group 0), then:
- Each thread in the warp has different oc_local ∈ [0, 7]
- For fixed ic, kd, kh, kw, weight[ic, oc_local, kd, kh, kw] are contiguous
- So all 32 threads in the warp can load 32 contiguous floats with one or two memory transactions! (Assuming float4 loads or even float loads are coalesced)
But we only need 8 output channels per group. What about the other threads in the warp?
We can have the warp compute 4 different spatial positions, each with 8 output channels. Or compute 8 output channels for 4 batches. Or compute 32 consecutive output channels which spans 4 groups.
If we have 32 threads computing channels 0-31, that's 4 groups. Each group has different input channels and different weights. The weight accesses for contiguous oc_local within a group are contiguous, but across groups, there's a discontinuity.
Actually, weight[ic, oc_local, ...] where ic ranges over 32 (all groups) and oc_local ranges over 8 per group. The memory layout is:
[ic=0, oc_local=0-7, kd, kh, kw], [ic=1, oc_local=0-7, ...], etc.
For a warp computing oc=0,1,2,...,31:
- oc=0 (g=0, ocl=0): uses weight[ic∈[0,7], 0, kd, kh, kw]
- oc=1 (g=0, ocl=1): uses weight[ic∈[0,7], 1, kd, kh, kw]
- ...
- oc=7 (g=0, ocl=7): uses weight[ic∈[0,7], 7, kd, kh, kw]
- oc=8 (g=1, ocl=0): uses weight[ic∈[8,15], 0, kd, kh, kw]
Within the warp, threads are accessing different memory locations, but for threads 0-7 (same group), their weight accesses for a fixed ic, kd, kh, kw are contiguous.
So for better coalescing, maybe have each warp compute 8 consecutive output channels for the SAME spatial position, and use 4 warps per block to cover 32 channels. Or have 32 threads in a block, each handling one output channel for the same spatial position.
Wait, but 32 threads per block is very few. We want at least 128 or 256 threads for good occupancy.
Another idea: process output in tiles where each block handles a spatial tile and all output channels for a subset of groups.
Let me think about a more concrete kernel design:
Block dimensions: (8, 8, 4) = 256 threads
- x dimension: 8 output channels within a group
- y dimension: 8 spatial positions (e.g., 2x2x2 or some other decomposition)
- z dimension: 4 batches or something
Actually, given that we only have 8 output channels per group, if we have 256 threads, we need to distribute them across multiple groups, batches, and spatial positions.
Option: Block shape (8, 4, 8)
- 8 output channels (all channels in one group)
- 4 spatial positions in depth
- 8 spatial positions in height*width
Total: 8*4*8 = 256 threads.
Each block computes 8 output channels for a tile of 4*8 = 32 spatial positions.
Within the block, we can load the relevant input and weights into shared memory.
For a spatial tile of (4_d, 2_h, 4_w) = 32 positions:
- Input depth needed: roughly 4/2 + 1 = 3 positions
- Input height needed: 2/2 + 1 = 2 positions (with some slack)
- Input width needed: 4/2 + 1 = 3 positions
- Input channels: 8
So input shared memory: 8 * 3 * 2 * 3 = 144 floats ≈ 576 bytes
For weights: we need all weights for 8 input channels and 8 output channels: 8*8*105 = 6720 floats ≈ 27 KB
Total shared memory ≈ 27.6 KB, very reasonable.
But wait, different spatial positions in the output tile need DIFFERENT subsets of the kernel offsets. For example, for od=0, we need kd=1. For od=1, we need kd=0,2.
If we load all weights (27 KB) into shared memory, that's fine. The shared memory limit is 227 KB, so we have plenty of room.
For input, the needed region depends on the output spatial tile. Let me be more precise.
If the block computes output for od ∈ [OD0, OD0+3], oh ∈ [OH0, OH0+1], ow ∈ [OW0, OW0+3]:
For each od in [OD0, OD0+3]:
Valid id range: [ceil((od + p_d - K_d + 1)/s_d), floor((od + p_d)/s_d)]
id_min = min over od of lower bound
id_max = max over od of upper bound
The union of all valid id's forms the needed input depth range.
Similarly for h and w.
For d: each od needs at most 2 id's, and 4 consecutive od's need about 3-4 id's.
For h: each oh needs about 3 ih's (K_h=5, s_h=2), and 2 consecutive oh's need about 4 ih's.
For w: each ow needs about 4 iw's (K_w=7, s_w=2), and 4 consecutive ow's need about 5 iw's.
So input tile roughly: 8 channels * 4 * 4 * 5 = 640 floats ≈ 2.5 KB
Total shared memory: 27.6 + 2.5 ≈ 30 KB. Still very manageable.
But is this worth it? The computation per block:
256 threads * (4*4*5*8) = 256 * 640 multiply-adds ≈ 164K FLOPs
At 67 TFLOPS, this is 2.4 microseconds per block.
Total blocks: 32 (groups*batches) * (24/4 * 48/2 * 96/4) = 32 * (6 * 24 * 24) = 32 * 3456 = 110,592 blocks
Hmm, that's a lot of blocks but each does very little work. The overhead might be high.
Actually, let me recalculate the number of blocks. If each block handles 1 group, 1 batch, and (4, 2, 4) spatial positions:
Number of group-batch combinations: 4 * 8 = 32
Spatial blocks: (24/4) * (48/2) * (96/4) = 6 * 24 * 24 = 3456
Total blocks: 32 * 3456 = 110,592
With 132 SMs, each SM processes about 830 blocks. If each block takes ~5 microseconds (including overhead), total time is about 4.1 ms. But there might be better ways to organize.
Actually, I should also consider that the actual multiply-adds might be less because not all 4*4*5 input positions are used for each output position. The exact number depends on the overlap pattern.
Actually, wait. For a block computing 32 output spatial positions and 8 output channels, each output spatial position only uses a subset of the input tile. So the total number of valid (input_spatial, kernel) combinations across all 32 output spatial positions is approximately 32 * 24 = 768. Then times 8 input channels = 6144 multiply-adds per output channel. Times 8 output channels = 49,152 multiply-adds per block.
With 256 threads, that's 192 multiply-adds per thread. At peak, we want thousands of FLOPs per thread to amortize overhead.
So maybe larger tiles are better. What if the block computes:
- All 8 output channels
- 8x8x8 = 512 spatial positions
Total threads: 256 (each thread computes 2 output channels? or 2 spatial positions?)
Actually, 512 spatial positions is too many for a single block to manage in shared memory.
What about computing 32 output channels (4 groups' worth) for a (4, 4, 4) spatial tile?
32 * 64 = 2048 output elements.
With 256 threads, each thread computes 8 output elements.
2048 * 192 = 393K multiply-adds per block. That's better.
But then each block needs weights for 4 groups: 4 * 27 KB = 108 KB shared memory, which is still under 227 KB.
And input for 4 groups: 4 * 2.5 KB ≈ 10 KB.
Total: ~118 KB shared memory. Still fits!
Actually wait, for 4 groups, the input channels are different (groups don't share input channels), but the spatial region is the same. So we need 32 input channels * spatial region. If spatial region is 4*4*5 ≈ 80 positions, that's 32*80 = 2560 floats ≈ 10 KB.
Total: 108 + 10 = 118 KB. This is well within limits and should give good occupancy (at most 1 block per SM due to shared memory, but that's fine if compute is sufficient).
With 118 KB per block and 227 KB limit, we can fit 1 block per SM. With 132 SMs and enough compute per block, this could work well.
But maybe I'm overcomplicating this. Let me first try a simpler approach: a warp-specialized or simple tiled kernel, and measure performance.
Actually, let me reconsider the problem. The PyTorch reference runtime might be around 1-5 ms for this problem. A custom kernel needs to beat that.
Let me look at what limits performance:
1. If cuDNN is very well optimized, it might be hard to beat
2. But for 3D conv with groups=4 and relatively small sizes, cuDNN might choose an algorithm with suboptimal performance
One thing that could give us a big win: since stride=2 and kernel sizes > stride, there are "holes" (zeros) in the dilated input. A naive im2col approach would expand the input and waste computation on zeros. Even cuDNN might not fully exploit the sparsity pattern if it's using generic GEMM.
However, I just realized: cuDNN backward_data for transposed conv is specifically designed for this. It should handle the stride implicitly without materializing zeros. So cuDNN might already be quite efficient.
But let's try anyway. I'll write an explicit kernel that directly computes the output without any unnecessary operations.
One optimization I can try: use vectorized loads (float4) for both input and weights when possible.
Another optimization: since many output pixels share the same input pixels, I can process output pixels in a way that maximizes cache reuse. For example, compute all output channels for a given spatial position before moving to the next spatial position. This keeps the input pixel in registers/cache while applying different weights.
Wait, that's an interesting idea. If I process one output spatial position at a time but compute all 32 output channels for it:
- For each valid (id, ih, iw, kd, kh, kw):
- Load input[b, ic, id, ih, iw] for all ic in the group
- For each oc in the group:
- output[b, oc, od, oh, ow] += input_value[ic] * weight[ic, oc_local, kd, kh, kw]
This way, each input value is loaded once and used for all 8 output channels in the group. That reduces input bandwidth by 8x!
And for weights, we load weight[ic, oc_local, kd, kh, kw] for all 8 output channels, which are contiguous in memory. Perfect for vectorized loads.
Let me formalize this:
For each output spatial position (od, oh, ow) and batch b:
For each group g:
Initialize accumulators for 8 output channels: acc[8] = {0}
For each valid (id, ih, iw, kd, kh, kw):
For ic_local in [0, 7]:
input_val = input[b, g*8+ic_local, id, ih, iw]
weight_vec = load 8 floats: weight[g*8+ic_local, 0:8, kd, kh, kw]
For oc_local in [0, 7]:
acc[oc_local] += input_val * weight_vec[oc_local]
For oc_local in [0, 7]:
output[b, g*8+oc_local, od, oh, ow] = acc[oc_local] + bias
Number of multiply-adds per output spatial position and group:
For each valid spatial combination (at most 24), and each ic_local (8), and each oc_local (8):
24 * 8 * 8 = 1536 multiply-adds.
But the number of weight loads is 24 * 8 = 192 vectors of 8 floats each. That's 1536 float loads, or 384 float4 loads.
Input loads: 24 * 8 = 192 float loads. But each input float is used for 8 multiplies (one per oc_local), so this is efficient.
Total float loads: 192 (input) + 1536 (weight) = 1728 floats = 6.9 KB per output spatial position, per group, per batch.
Output writes: 8 floats = 32 bytes.
For all 28M output elements, if organized this way, we'd load each input element multiple times but each time use it efficiently.
Actually, this approach has a nice property: for each input channel, we load its value once per valid spatial position, and then multiply by 8 weights (for the 8 output channels) before moving to the next input channel or spatial position.
This is quite efficient! And if we process multiple output spatial positions in parallel using a block of threads, we can share the weight loads.
Let me design a block-level kernel:
Block dimensions: (8, 32) = 256 threads
- 8 threads per warp correspond to 8 output channels in a group
- 32 "warps" or thread groups correspond to 32 output spatial positions? No, 32*8=256 threads total.
Actually, let's use 256 threads organized as:
- 8 threads along oc_local dimension (output channels in group)
- 32 threads along some spatial dimension
Each block computes 32 output spatial positions, each with 8 output channels.
But the 32 spatial positions should be chosen to maximize weight sharing. If they're nearby spatial positions, they need different weight slices (because different (kd, kh, kw) are used for different od, oh, ow).
Actually, nearby output positions DO share some input positions. For example, od=0 and od=1 both might use id=0 or id=1, but with different kd values.
Hmm, the weight sharing for nearby output positions is tricky because each output position needs a different subset of kernel offsets.
Alternative: process the SAME output spatial position with multiple threads, where threads compute different groups.
For example, 256 threads computing:
- 32 output channels (4 groups)
- 8 spatial positions
Or 256 threads computing 32 output channels for 1 spatial position... but then 32*1 = 32 threads are enough per spatial position.
Let me think about this differently. The "all output channels for one spatial position" approach is very sequential in nature. To parallelize it:
Each block computes:
- 8 output channels for one group
- A tile of spatial positions
Threads within the block:
- Each thread computes one output channel for a subset of spatial positions
- Or each thread computes all output channels for one spatial position
Option A: 8 threads compute all 8 output channels for one spatial position (sequential over spatial positions), and we have 32 such teams in parallel (256/8 = 32).
So block computes 32 spatial positions, each with all 8 channels.
Within each team of 8 threads, they jointly compute one output spatial position by sharing the input values and each handling one output channel. But then each thread needs different weights.
Actually, within a team, thread t computes output channel t. For a given (id, ih, iw, kd, kh, kw, ic_local):
- All threads load the same input value
- Thread t loads weight[ic_local, t, kd, kh, kw]
Since weight[ic_local, 0:8, kd, kh, kw] are contiguous, the 8 threads in the team can load them cooperatively with one float42 (8 float) vector load, or two float4 loads.
This is excellent! Within each team of 8 threads:
- Thread 0 loads input value
- All 8 threads collaboratively load weight vector (4 float4 loads or similar)
- Each thread does one multiply-add
Then move to the next input channel/spatial combination.
For the block with 32 teams (256 threads total), each team handles one output spatial position. The block computes 32 output spatial positions in parallel.
But wait, the 32 teams are computing different spatial positions, so they need different input values. And because different spatial positions need different valid (kd, kh, kw) combinations, the control flow might diverge.
However, for our specific sizes with stride=2 and kernel sizes (3,5,7), all output positions have the SAME number of valid combinations: 2*3*4 = 24. The specific (kd, kh, kw) values differ, but the loop count is the same.
Wait, is the loop count exactly 24 for ALL output positions? Let's check boundaries.
For od=0: id can only be 0 (kd=1). That's 1 combination.
For od=1: id can be 0 (kd=2) or 1 (kd=0). That's 2 combinations.
Actually, the number of valid combinations varies near the boundaries!
For d: K_d=3, s_d=2, p_d=1.
Most positions have valid id count = ceil(K_d/s_d) = 2. But boundaries might have 1.
Actually, the number of valid id's for a given od is the number of integers in [ceil((od+p-K+1)/s), floor((od+p)/s)].
For od=0: [ceil(-0.5), 0] = [0,0] → 1 id
For od=1: [ceil(-0.5), 1] = [0,1] → 2 ids
For od=2: [ceil(0.5), 1] = [1,1] → 1 id
For od=3: [ceil(0.5), 2] = [1,2] → 2 ids
etc.
So the pattern alternates between 1 and 2 for d! Similarly for h and w, the count varies between 2 and 3, or 3 and 4.
The total number of spatial combinations varies! For example:
- od=0 (1) * oh=0 (let's say 2) * ow=0 (let's say 3) = 6 combinations
- od=1 (2) * oh=0 (2) * ow=0 (3) = 12 combinations
This means control flow can diverge between threads computing different output spatial positions.
To avoid divergence, maybe process output positions with the same "phase" together.
For d: od % 2 determines whether there are 1 or 2 valid ids.
For h: oh % 2 determines the count.
For w: ow % 2 determines the count.
So there are 2*2*2 = 8 "phase classes". If we process only output positions in the same phase class within a warp, the loop count is the same.
Actually, even within the same phase class, the specific (kd, kh, kw) values differ, so the weight indices differ. But the loop count is the same, which avoids divergent execution.
For our tiling, we can assign each block to one phase class, and within the block, process multiple spatial positions of that phase.
Let me count how many positions are in each phase class:
D_out=24, H_out=48, W_out=96
d phase 0 (even): 12 positions
d phase 1 (odd): 12 positions
h phase 0 (even): 24 positions
h phase 1 (odd): 24 positions
w phase 0 (even): 48 positions
w phase 1 (odd): 48 positions
Total positions per phase class: 12*24*48 = 13,824
With 32 teams per block, we can process 32 spatial positions per block. Number of blocks per phase class: 13,824 / 32 = 432.
Total blocks: 8 * 432 = 3456 spatial-blocks.
Times 4 groups * 8 batches = 32 group-batch combinations.
Total blocks: 3456 * 32 = 110,592 blocks.
This is the same total as before. With 132 SMs, about 838 blocks per SM.
Now, if each block computes 32 spatial positions * 8 channels = 256 output elements, and each output element needs on average maybe 100-200 multiply-adds, the total compute per block is 25K-50K FLOPs. At 67 TFLOPS, this is ~0.4-0.7 microseconds per block. With 838 blocks per SM, total time ~0.3-0.6 ms.
This seems very fast. But there are overheads:
1. Loading input and weights
2. Shared memory allocation
3. Thread synchronization
And memory bandwidth might be the bottleneck, not compute. Let me estimate memory traffic.
For each block (32 spatial positions, 8 output channels, 1 group, 1 batch):
Input needed: depends on the spatial tile. For 32 positions of the same phase, they're spread out. To share input loading, the positions should be nearby.
Actually, if we choose 32 positions that are spatially contiguous within the phase class, their union of input regions might be quite large.
For example, if we choose a 4x4x2 tile of positions all with the same phase:
- The od positions are 0, 2, 4, 6 (spacing of 2)
- The oh positions are 0, 2, 4, 6
- The ow positions are 0, 2 (or some spacing)
Actually, for positions to have the same phase, their coordinates must all be even or all be odd. If we take a dense block of same-phase positions, like od=0,2,4,6 and oh=0,2,4,6 and ow=0,2,4,6, that's only 4x4x4 = 64 positions, which is more than 32.
If we take od∈{0,2}, oh∈{0,2,4,6}, ow∈{0,2,4,6}: that's 2*4*4 = 32 positions.
The union of input depth needed:
For od=0: id=0
For od=2: id=1
So input depths: {0, 1}
For oh=0: ih=0 (maybe, need to compute exact)
For oh=2: ih=1
etc.
Input heights: maybe {0, 1}
Input widths: maybe {0, 1, 2}
So input tile: 8 channels * 2 * 2 * 3 = 96 floats ≈ 384 bytes.
Hmm, that's very small! Actually, I think I'm underestimating. Let me do a more careful calculation.
For a given od, the valid id's are:
If od is even (say od=0,2,4,...):
For od=0: id ∈ [0,0], id=0, kd=1
For od=2: id ∈ [1,1], id=1, kd=1
For od=4: id ∈ [2,2], id=2, kd=1
So for even od, id = od/2, kd = 1.
If od is odd (say od=1,3,5,...):
For od=1: id ∈ [0,1], id=0(kd=2), id=1(kd=0)
For od=3: id ∈ [1,2], id=1(kd=2), id=2(kd=0)
So for odd od, id = (od-1)/2 and (od+1)/2, with kd=2 and 0 respectively.
Interesting! For even od, there's exactly 1 valid id and kd=1.
For odd od, there are exactly 2 valid ids, with kd=2 or 0.
Similarly for h:
K_h=5, s_h=2, p_h=2.
For even oh:
oh=0: id ∈ [ceil((0+2-5+1)/2), floor((0+2)/2)] = [ceil(-1), 1] = [-1? wait let me recalculate]
Actually id must be >= 0.
oh=0: [ceil((0+2-5+1)/2), floor(2/2)] = [ceil(-1), 1] = [-1, 1]
But id >= 0, so [0, 1].
id=0: kh = 0+2-0 = 2
id=1: kh = 0+2-2 = 0
oh=2: [ceil((2+2-5+1)/2), floor((2+2)/2)] = [ceil(0), 2] = [0, 2]
id=0: kh=2+2-0=4
id=1: kh=2+2-2=2
id=2: kh=2+2-4=0
For odd oh:
oh=1: [ceil((1+2-5+1)/2), floor((1+2)/2)] = [ceil(-0.5), 1] = [0, 1]
id=0: kh=1+2-0=3
id=1: kh=1+2-2=1
oh=3: [ceil((3+2-5+1)/2), floor((3+2)/2)] = [ceil(0.5), 2] = [1, 2]
id=1: kh=3+2-2=3
id=2: kh=3+2-4=1
So for any oh, there are either 2 or 3 valid ih's.
For w: K_w=7, s_w=2, p_w=3.
ow=0: [ceil((0+3-7+1)/2), floor((0+3)/2)] = [ceil(-1.5), 1] = [-1, 1], clipped to [0, 1]
id=0: kw=3
id=1: kw=1
ow=1: [ceil((1+3-7+1)/2), floor((1+3)/2)] = [ceil(-1.5), 2] = [0, 2]
id=0: kw=4
id=1: kw=2
id=2: kw=0
So there are 2 or 3 valid iw's for each ow.
For our block computing 32 positions with the same parity (all even or all odd):
If all d positions are even, each has exactly 1 valid id.
If all d positions are odd, each has exactly 2 valid ids.
This means within a block, all threads have the same loop count if they're in the same phase class.
But the number of id's differs between even and odd phases for d. So we should definitely process even and odd d positions in separate blocks.
So 8 phase classes becomes... actually, for d we must separate even and odd, but for h and w we can mix? No, better to separate all phases for maximum uniformity.
Actually, for the most uniform blocks, we'd have:
- Phase where d is even: each od has 1 id
- Phase where d is odd: each od has 2 ids
- For h: each oh has either ceil(5/2)=3 or floor(5/2)=2 ids depending on the boundary
- For w: each ow has either ceil(7/2)=4 or floor(7/2)=3 ids
The exact count for h and w depends on oh and ow values, not just parity! For example:
oh=0: 2 valid ih's (0, 1)
oh=1: 2 valid ih's (0, 1)
oh=2: 3 valid ih's (0, 1, 2)
oh=3: 2 valid ih's (1, 2)
Wait, let me recheck:
oh=0: [0,1] → 2
oh=1: [0,1] → 2
oh=2: [0,2] → 3
oh=3: [1,2] → 2
oh=4: [1,3] → 3? Let me check:
[ceil((4+2-5+1)/2), floor((4+2)/2)] = [ceil(1), 3] = [1, 3] → 3
So the count alternates between 2 and 3, but with a pattern of [2, 2, 3, 2, 3, 2, 3, ...]?
Actually: oh=0→2, oh=1→2, oh=2→3, oh=3→2, oh=4→3, oh=5→2, oh=6→3...
For w:
ow=0: [0,1] → 2
ow=1: [0,2] → 3
ow=2: [1,2] → 2
ow=3: [1,3] → 3
ow=4: [2,3] → 2
ow=5: [2,4] → 3
...
For w: [2, 3, 2, 3, 2, 3, ...]
So h and w counts depend on the exact coordinate modulo some pattern. For perfect uniformity, we'd need very fine-grained phase classes.
Actually, for h:
Lower bound: ceil((oh + 2 - 5 + 1) / 2) = ceil((oh - 2) / 2)
Upper bound: floor((oh + 2) / 2)
For oh = 0: ceil(-1) = -1 → 0 (clipped), floor(1) = 1 → [0,1], count=2
oh = 1: ceil(-0.5) = 0, floor(1.5) = 1 → [0,1], count=2
oh = 2: ceil(0) = 0, floor(2) = 2 → [0,2], count=3
oh = 3: ceil(0.5) = 1, floor(2.5) = 2 → [1,2], count=2
oh = 4: ceil(1) = 1, floor(3) = 3 → [1,3], count=3
So pattern: 2, 2, 3, 2, 3, 2, 3, ...
Actually: 0→2, 1→2, 2→3, 3→2, 4→3, 5→2, 6→3, 7→2, 8→3...
For even oh >= 2: count is 3 if oh/2 >= 1, i.e., oh >= 2.
Wait, oh=2: 3, oh=4: 3, oh=6: 3. All even oh >= 2 have count 3?
oh=2: [0,2]→3 ✓
oh=4: [1,3]→3 ✓
oh=6: [2,4]→3 ✓
For odd oh: oh=1→2, oh=3→2, oh=5→2. All odd oh have count 2.
For w:
ow=0: [0,1]→2
ow=1: [0,2]→3
ow=2: [1,2]→2 (wait, let me recheck)
[ceil((2+3-7+1)/2), floor((2+3)/2)] = [ceil(-0.5), 2] = [0, 2] → 3
ow=3: [ceil((3+3-7+1)/2), floor((3+3)/2)] = [ceil(0), 3] = [0, 3] → 4? Let me recheck.
(3+3-7+1)/2 = 0/2 = 0. ceil(0) = 0.
floor((3+3)/2) = floor(3) = 3.
[0,3] → 4 values!
Hmm wait, I need to check if id=0,1,2,3 are all valid for D=24 inputs (actually D=12 for input depth, but W=48 for input width).
ow=3, id=3: kw = 3 + 3 - 3*2 = 6 - 6 = 0. ✓
id must be < W_in = 48. id=3 is valid.
ow=4: [ceil((4+3-7+1)/2), floor((4+3)/2)] = [ceil(0.5), 3] = [1, 3] → 3
ow=5: [ceil((5+3-7+1)/2), floor((5+3)/2)] = [ceil(1), 4] = [1, 4] → 4
So for w:
even ow: 2, 3, 3, 3, ...
ow=0: 2
ow=2: 3
ow=4: 3
ow=6: [ceil((6+3-7+1)/2), floor((6+3)/2)] = [ceil(1.5), 4] = [2, 4] → 3
odd ow: 3, 4, 4, 4, ...
ow=1: 3
ow=3: 4
ow=5: 4
Wait, ow=3 gives count 4? Let me double check:
[0, 3] for iw, that's iw=0,1,2,3.
kw = ow + p - iw*s = 3 + 3 - iw*2 = 6 - 2*iw
iw=0: kw=6 (but K_w=7, so kw=6 is valid: 0-6)
iw=1: kw=4
iw=2: kw=2
iw=3: kw=0
All valid! So 4 combinations.
For ow=5:
[1, 4]: iw=1,2,3,4
kw = 5+3-2*iw = 8-2*iw
iw=1: kw=6
iw=2: kw=4
iw=3: kw=2
iw=4: kw=0
All valid! Count=4.
For ow=7:
[ceil((7+3-7+1)/2), floor((7+3)/2)] = [ceil(2), 5] = [2, 5] → 4
iw=2,3,4,5: 4 values.
So for w, the counts are:
ow=0: 2
ow=1: 3
ow=2: 3
ow=3: 4
ow=4: 3
ow=5: 4
ow=6: 3
ow=7: 4
...
Pattern: ow=0→2, ow>=1→3 or 4.
Specifically, for even ow>=2: count=3
For odd ow>=1: count=4? Let's check ow=1: 3, not 4.
Hmm, let me recheck ow=1:
[ceil((1+3-7+1)/2), floor((1+3)/2)] = [ceil(-1), 2] = [-1, 2], clipped to [0, 2] → 3
ow=3: [0, 3] → 4
ow=5: [1, 4] → 4
So:
ow=0: [0,1] → 2
ow=1: [0,2] → 3
ow=2: [0,2] → 3
ow=3: [0,3] → 4
ow=4: [1,3] → 3
ow=5: [1,4] → 4
ow=6: [2,4] → 3
ow=7: [2,5] → 4
Pattern for even ow:
ow=0: 2
ow=2: 3
ow=4: 3
ow=6: 3
...
ow=2k for k>=1: 3
Pattern for odd ow:
ow=2k+1:
ow=1: 3
ow=3: 4
ow=5: 4
ow=7: 4
...
ow=2k+1 for k>=1: 4
So for w, only ow=0 is special with count 2. All others are 3 or 4.
Actually, let me verify ow=2:
[ceil((2+3-7+1)/2), floor((2+3)/2)] = [ceil(-0.5), 2] = [0, 2] → 3
iw=0: kw=5
iw=1: kw=3
iw=2: kw=1
All valid (0-6). ✓
And ow=96-1=95:
[ceil((95+3-7+1)/2), floor((95+3)/2)] = [ceil(46), 49] = [46, 49] → 4
iw=46: kw=95+3-92=6
iw=47: kw=4
iw=48: kw=2
iw=49: kw=0
But wait, W_in=48, so iw must be < 48. iw=48 and 49 are out of range!
Ah! I forgot to clip to the actual input size. So for ow=95:
Upper bound floor((95+3)/2) = 49, but W_in=48, so upper bound is 47.
[46, 47] → 2 values.
iw=46: kw=6
iw=47: kw=4? 95+3-94=4. But kw=4 is valid.
Actually wait, iw=47 gives kw=95+3-47*2=98-94=4 ✓
So for large ow near the boundary, the count drops back down.
This means the loop count varies significantly across the output domain, making warp-level uniformity challenging.
Given all this complexity, maybe the simplest approach is actually best: one thread per output element, with the exact iteration count for that element. The hardware can handle moderate divergence, and for most of the output domain, the counts are near the maximum.
For the interior regions:
- d: 2 for odd, 1 for even
- h: 3 for even >= 2, 2 for odd
- w: 3 for even >= 2, 4 for odd >= 3
Interior combinations:
For even d (1), even h (3), even w (3): 9 combinations
For odd d (2), even h (3), even w (3): 18 combinations
For even d (1), odd h (2), even w (3): 6 combinations
For even d (1), even h (3), odd w (4): 12 combinations
...
Average: roughly between 9 and 24, maybe around 15 for most positions.
So each thread does about 15 * 8 * 8 = 960 multiply-adds on average. Not too bad.
Actually, there's an even better approach. Instead of iterating over valid (id, ih, iw) and deriving kd, kh, kw, we can iterate over id, ih, iw and derive kd, kh, kw. Then we check if kd, kh, kw are in [0, K-1].
For a given od:
id_min = max(0, ceil((od + p - K + 1) / s))
id_max = min(D-1, floor((od + p) / s))
For each id in [id_min, id_max]:
kd = od + p - id * s
if 0 <= kd < K:
valid!
This gives the same result and might be simpler to compute in the kernel.
Similarly for h and w.
For our test case, since the range is small, we can just iterate id from 0 to D-1 and check the condition. But that would be 12 iterations for d instead of ~2. However, if we unroll or use the exact range, it's fine.
Wait, actually computing id_min and id_max once per thread, then iterating over that small range, is quite efficient.
OK, let me write the kernel. I'll start with a relatively simple approach:
- Each thread computes one output channel for one spatial position
- Use 256 threads per block
- Block computes a tile of output
Grid organization:
- x: batches * out_channels (8*32 = 256)
- y: spatial positions (24*48*96 = 110,592)
- Total threads: 256 * 110,592 = 28,311,552 (one per output element)
Wait, that's exactly one thread per output element. But with 256 threads per block, we'd have 110,592 blocks in the y dimension. That's a lot.
Actually, we can combine spatial dimensions:
- x: out_channels (32)
- y: batch * D_out * H_out * W_out / some_tile_size
If we use 256 threads per block, and each thread handles one output element, then:
- blocks along x: 32 / 32 = 1 warp for channels? No, 32 channels per block means 32 threads for channels.
- Remaining 224 threads for spatial positions.
- We can have 7 warps, each handling 32 spatial positions? No, 224 threads / 32 = 7, but that's not clean.
Better: 256 threads per block, each computing 1 output element.
Block dimensions:
- x: 32 threads for 32 output channels... but groups=4, so each warp would span multiple groups.
- y: 8 threads for something
- z: not needed, or could use for batches
Actually, let me try a different layout. Each block computes:
- 16 output channels (2 groups' worth)
- 16 spatial positions
Total: 16*16 = 256 threads.
Block dimensions: (16, 16) = 256 threads.
- x: 16 output channels (2 groups of 8)
- y: 16 spatial positions
For x=16 channels: channels 0-7 are group 0, channels 8-15 are group 1.
Within the block, each thread computes one output element. The spatial positions can be organized as a 4x4 tile, or 2x2x4, etc.
Let me think about memory access within a warp. A warp of 32 threads consists of 2 "rows" in the block if we use (16, 16) layout, meaning:
Threads (0-15, 0) and (0-15, 1) form the first warp.
These threads all have x ∈ [0,15] (channels 0-15) and y ∈ {0, 1} (2 spatial positions).
So within a warp, 16 threads compute channels 0-7 and 8-15 for spatial position 0, and 16 threads compute the same channels for spatial position 1.
For spatial position 0 and 1, the valid (id, ih, iw) combinations might differ, so the loop count differs. But within each half-warp (16 threads), they all compute the same spatial position with different channels.
For the 16 threads computing the same spatial position and different channels (0-15 spanning 2 groups):
- They need different weights (different oc_local and different group)
- For a fixed ic and spatial combination, weight[ic, :, kd, kh, kw] for channels 0-7 are contiguous, but channels 8-15 are in a different group.
So loading weights for the first 8 channels is coalesced, but channels 8-15 need a different weight region.
This is getting complex. Let me simplify.
How about: each block computes 8 output channels (1 group) for 32 spatial positions?
Block: (8, 32) = 256 threads.
A warp (32 threads) would be:
- 8 threads for channels 0-7 at spatial position y=0 to 3 (if we tile y as 32 positions)
Or more precisely, with block shape (8, 32):
Warp 0: threads (0-7, 0-3) = 8*4 = 32 threads.
All threads in warp 0 compute spatial positions y=0,1,2,3 with channels 0-7.
But wait, with block shape (8, 32), the thread index is threadIdx.x + threadIdx.y * 8.
Warp 0: threads 0-31, which are (x=0-7, y=0-3).
So within warp 0, 4 threads compute spatial position y=0 with channels 0,4,8,12? No, all 8 threads compute all channels 0-7 for their respective y.
Actually, within warp 0:
- Thread (0,0), (1,0), ..., (7,0): spatial pos y=0, channels 0-7
- Thread (0,1), (1,1), ..., (7,1): spatial pos y=1, channels 0-7
- Thread (0,2), ..., (7,2): spatial pos y=2, channels 0-7
- Thread (0,3), ..., (7,3): spatial pos y=3, channels 0-7
So there are 4 groups of 8 threads within the warp, each group computing a different spatial position but the same set of channels (0-7).
For each spatial position, the 8 threads computing channels 0-7 can load weights coalesced for a fixed ic, kd, kh, kw.
And the 4 spatial positions within the warp might have different loop counts (due to boundary effects), but the counts are small and the hardware handles warp divergence.
The key benefit: for a fixed spatial combination (id, ih, iw, kd, kh, kw) and fixed ic:
- 4 threads (one per spatial position) load the same input value? No, different spatial positions need different input values.
- But 8 threads for channels load contiguous weights.
So the weight loading is perfectly coalesced within each sub-warp of 8 threads. The input loading is not coalesced within the same sub-warp because they need different spatial positions, but across the 4 sub-warps, they might have some locality.
Actually, for spatial positions y=0,1,2,3, if these correspond to different od values with the same oh and ow, then the input id values differ, so input loads are to different addresses.
But if y=0,1,2,3 correspond to consecutive ow with the same od and oh, then the input iw values might be overlapping (because consecutive ow values share some iw values).
This suggests that organizing spatial positions along the w dimension might be best for locality, since w has the most positions and the most overlap.
Actually, it might be even better to use a 3D spatial tile for y dimension.
Let me try: block shape (8, 4, 8) = 256 threads.
- x: 8 channels
- y: 4 depth positions
- z: 8 (height*width) positions, e.g., 2x4 or something
A warp would be 32 threads. With (8, 4, 8), the warp could cover:
- All 8 channels (x)
- 1 depth position (y)
- 4 height*width positions (z)
Or with the hardware warp layout, it depends on how thread IDs map.
Actually, in CUDA, threads in a warp have consecutive threadIdx.x values if the block dimensions are organized appropriately. But the exact mapping of 3D thread indices to 1D warp lanes is:
tid = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y
For (8, 4, 8):
tid = x + y*8 + z*32, where x∈[0,7], y∈[0,3], z∈[0,7].
Warp 0: tid 0-31 → x∈[0,7], y∈[0,3], z=0.
So warp 0 covers channels 0-7, depth 0-3, and spatial position 0.
This means all 32 threads in warp 0 compute the same spatial position (z=0) but different channels and depth. Wait no, z=0 is one height*width position, but y=0,1,2,3 are 4 different depth positions.
So warp 0 computes 4 depth positions (y=0,1,2,3) at the same (h, w) position z=0, for 8 channels each.
For a fixed channel and fixed (h,w):
The 4 threads for different depths load input from different id's. These might be nearby or far apart in memory depending on depth stride.
Input memory layout in PyTorch is (N, C, D, H, W). So input[b, ic, id, ih, iw] and input[b, ic, id+1, ih, iw] are separated by H*W floats.
For H=24, W=48, H*W=1152 floats ≈ 4.5 KB. So consecutive depth positions are 4.5 KB apart. In L1 cache, this might be ok if the cache line is large enough, but they're not contiguous.
Hmm. What if we organize the block as (8, 1, 32)?
Warp 0: tid 0-31 → x∈[0,7], y=0, z∈[0,3] (since 8*1*4=32).
Warp 0 covers 8 channels and 4 spatial positions (all at y=0).
Or (32, 1, 8):
Warp 0: tid 0-31 → x∈[0,31], y=0, z=0.
All 32 threads compute the same spatial position with 32 different channels. But we only have 32 channels total, and groups=4. So x=0-7 is group 0, x=8-15 is group 1, x=16-23 is group 2, x=24-31 is group 3.
For a fixed spatial position and 32 channels (4 groups), all threads load the same input values but different weights. The weights for channels 0-7 (group 0) are contiguous, but then there's a gap to group 1, etc.
Actually, with groups=4, the weight tensor is (32, 8, 3, 5, 7). The layout is:
[group 0: ic=0-7, oc=0-7], [group 1: ic=8-15, oc=0-7], etc.
For a fixed oc_local=0 and varying ic:
weight[0, 0, ...], weight[1, 0, ...], ..., weight[7, 0, ...] are contiguous.
Then weight[8, 0, ...] comes next.
So for a warp computing channel 0 (oc=0) and channel 8 (oc=8), thread 0 loads weight[ic=0, oc_local=0, ...] and thread 8 loads weight[ic=8, oc_local=0, ...]. These are separated by 8*3*5*7 = 840 floats.
If ic varies and oc_local varies too, the pattern is more complex.
Given all this analysis, I think the most important optimization is to:
1. Minimize the number of multiply-adds by only iterating over valid combinations
2. Ensure coalesced weight loading when possible
3. Use shared memory if beneficial
For a first attempt, let me write a simple kernel where each thread handles one output element, and we iterate over the valid combinations. This avoids complex shared memory management and is likely to give decent speedup over a generic cuDNN call for this problem size.
Actually, I just realized something. For ConvTranspose3d, cuDNN might choose an algorithm that's optimized for larger batch sizes or different configurations. With batch=8, channels=32, and 3D spatial (12*24*48), the problem might be too small for cuDNN to fully utilize the GPU, and the algorithm selection overhead might be significant.
A custom kernel that directly implements the operation with minimal overhead and explicit loop bounds might indeed be faster.
Let me write the kernel. I'll organize threads as one per output element, with the grid covering all output elements.
Grid: (batch_size, out_channels, D_out * H_out * W_out)
Block: 256 threads, mapping to a subset of output elements.
More precisely:
Total output elements per batch and channel: 24*48*96 = 110,592
With 256 threads per block, that's 432 blocks.
Total threads needed: 8 * 32 * 110,592 = 28,311,552.
CUDA grid dimensions can be up to 2^31-1 in x, 65535 in y and z. So I could use:
- x: 110,592 spatial positions (too large for some GPUs, but should be ok on H100)
- y: 32 channels
- z: 8 batches
Or combine into fewer dimensions:
- x: D_out * H_out * W_out = 110,592
- y: out_channels = 32
- z: batch_size = 8
With block size 256 in x dimension:
Blocks in x: 110,592 / 256 = 432
So grid = (432, 32, 8), block = (256, 1, 1).
Total: 432 * 32 * 8 = 110,592 blocks.
Actually... with block=(256,1,1) and grid=(432,32,8):
Total threads = 256 * 432 * 32 * 8 = 28,311,552. Correct.
Within each block:
- threadIdx.x determines the spatial position within the block's tile.
- blockIdx.x determines which spatial tile (0 to 431).
- blockIdx.y determines output channel (0 to 31).
- blockIdx.z determines batch (0 to 7).
Each thread computes one output element: output[b, oc, od, oh, ow].
Wait, but blockIdx.y determines the output channel, so ALL threads in the block compute the SAME output channel. Is this good?
If all threads compute the same output channel, they all load the same weights for a given ic and spatial combination. But they load different input values for different spatial positions.
For weights: if all threads need weight[ic, oc_local, kd, kh, kw] for the SAME oc_local, then a single load can be broadcast to all threads in the block! Actually, the compiler might cache the weight in a register and reuse it, but since oc_local is the same for all threads in the block, there's potential for constant memory caching or L1 broadcast.
Wait, but different threads in the block compute different spatial positions, so for the SAME oc_local but DIFFERENT (kd, kh, kw). So they need different weights.
Hmm. Let me think. If two threads compute spatial positions that share the same (kd, kh, kw) for a given (id, ih, iw), they need the same weight. But in general, different spatial positions need different (kd, kh, kw).
Actually, if all threads in the block have the SAME oc but different spatial positions, then:
- For a fixed ic, the weight is weight[ic, oc_local, kd, kh, kw] where kd, kh, kw depend on the thread's spatial position.
- These are different for different threads, so weight accesses are not uniform.
But if we organize threads to have DIFFERENT oc but same spatial position:
- For a fixed ic and fixed (kd, kh, kw), weight[ic, oc_local, kd, kh, kw] for different oc_local are contiguous.
- All threads load the same input value.
This seems much better for coalescing!
So let me reorganize:
- Threads in the SAME warp/block compute the SAME spatial position but DIFFERENT output channels.
- We process spatial positions sequentially across warps/blocks.
Block dimensions: (32, 8) = 256 threads.
- x: 32 output channels (all 32 channels for the same spatial position)
- y: 8 spatial positions (but warps would span this)
Wait, with block=(32,8):
tid = x + y*32.
Warp 0: tid 0-31 → x∈[0,31], y=0.
Warp 1: tid 32-63 → x∈[0,31], y=1.
...
Warp 7: tid 224-255 → x∈[0,31], y=7.
Each warp computes 32 output channels for ONE spatial position.
Within the warp, 32 threads load weights for 32 different output channels. Are these contiguous?
Weight layout: (32, 8, 3, 5, 7).
For fixed ic, fixed kd,kh,kw, and varying oc:
weight[ic, 0, kd, kh, kw], weight[ic, 1, kd, kh, kw], ..., weight[ic, 7, kd, kh, kw] are contiguous.
But we have 32 output channels, which span 4 groups. For group 0, oc_local=0-7 with ic∈[0,7]. For group 1, oc_local=0-7 with ic∈[8,15].
So for a fixed ic and varying oc:
- If ic is in group 0 (0-7), only oc in group 0 (0-7) use this ic. Channels 8-31 are in different groups and don't use ic∈[0,7].
- So if a warp computes channels 0-31, and we're iterating over ic, only some channels participate for each ic.
This means we should process one group at a time per block, or at least per warp.
Block dimensions: (8, 32) = 256 threads.
- x: 8 output channels in one group
- y: 32 spatial positions
With block=(8,32):
tid = x + y*8.
Warp 0: tid 0-31 → x∈[0,7], y∈[0,3].
So 4 spatial positions, each with 8 channels.
Within warp 0, for a fixed ic and fixed spatial combination:
- 4 threads (one per spatial position, same channel offset) load the same weight value? No, 4 threads per channel means 4 spatial positions, and threads with the same x but different y load different weights if kd,kh,kw differ.
- But 8 threads with different x and same y load contiguous weights for the same spatial position.
Actually, within warp 0, for fixed y=0 (spatial pos 0), threads x=0-7 load weight[ic, 0:7, kd, kh, kw] which are contiguous! Perfect coalescing!
And for the other 3 spatial positions (y=1,2,3), similar coalescing occurs within each group of 8 threads.
But the 4 spatial positions might have different numbers of valid combinations, causing divergence.
Given that the valid combinations for nearby spatial positions are highly correlated and the counts don't vary too much (mostly between 9 and 24), the divergence might be acceptable.
Actually, for maximum efficiency within a warp, all 32 threads should have the same control flow. So ideally, all 32 threads in a warp compute spatial positions with the same number of valid combinations.
One way: within a warp, use 8 threads for 8 channels and 4 spatial positions that are guaranteed to have the same loop count.
We can choose spatial positions with the same (od % 2, oh % 2, ow % 2) parity if the counts depend on parity. But as I showed earlier, for h and w, the count depends on the exact position, not just parity.
However, the counts are the same for most positions. For w, only positions near the boundary have count=2 instead of 3 or 4. For d, even positions have count=1 and odd have count=2. For h, even>=2 have count=3, odd have count=2, and some boundaries have 1.
To simplify, let me just not worry too much about divergence and write the kernel. H100 has good divergence handling.
Actually, another idea: process the SAME spatial position with a whole block, where 256 threads compute 256 output channels. But we only have 32 output channels total! So that's not enough parallelism.
Unless we process the same spatial position for multiple batches simultaneously. 8 batches * 32 channels = 256. Perfect!
Block dimensions: (32, 8) = 256 threads.
- x: 32 output channels (0-31)
- y: 8 batches (0-7)
All threads in the block compute the SAME spatial position (od, oh, ow) but different (batch, channel) combinations.
Now, for a fixed spatial combination (id, ih, iw, kd, kh, kw) and fixed ic:
- All threads load the SAME input value (but from different batches). If batches are stored consecutively in memory, input[b, ic, id, ih, iw] for b=0-7 are separated by C*D*H*W. So consecutive batch indices are 32*12*24*48 = 442,368 floats apart ≈ 1.7 MB. That's not contiguous at all.
Hmm. With PyTorch layout (N, C, D, H, W), elements from the same spatial position in different batches are very far apart. That's bad for coalescing if threads load different batches.
So processing multiple batches in parallel is bad for input coalescing.
What if a block processes multiple spatial positions for ONE batch?
Block=(32, 8), but x=32 channels, y=8 spatial positions.
Warp 0: x=0-31, y=0. Computes 32 channels for spatial pos 0.
Warp 1: x=0-31, y=1. Computes 32 channels for spatial pos 1.
...
For each warp, 32 threads compute 32 channels for the same spatial position.
Within the warp, for fixed ic and fixed spatial combination:
- Threads load weight[ic, oc_local, kd, kh, kw] for oc_local=0-7 in each group.
- But since channels 0-7, 8-15, 16-23, 24-31 are in different groups, and each group uses different ic ranges, this is complex.
Actually, for group g, channels are g*8 to g*8+7, and input channels are also g*8 to g*8+7.
So when processing channel 0 (group 0), we iterate over ic=0-7.
When processing channel 8 (group 1), we iterate over ic=8-15.
If a warp processes channels 0-31 simultaneously, it would need to iterate over ic=0-31, but only some channels participate for each ic. Specifically:
- For ic=0, channels 0-7 participate (group 0), channels 8-31 do nothing.
- For ic=8, channels 8-15 participate (group 1), channels 0-7 and 16-31 do nothing.
This would cause significant warp divergence! Not good.
So within a warp, we should compute channels from the SAME group.
Block=(8, 32):
- x: 8 channels in one group
- y: 32 spatial positions
Warp 0: x=0-7, y=0-3. All threads in warp 0 compute group 0 channels for 4 spatial positions.
For a fixed spatial combination, all 32 threads in warp 0 use the same ic range (0-7) and different spatial positions/channels.
When iterating over ic from 0 to 7:
- For each ic, all 32 threads are active (computing different spatial positions and channels within group 0).
- They load input[b, ic, id, ih, iw] for different spatial positions. Depending on how spatial positions are assigned, the id, ih, iw values might have some locality.
If we assign y=0-31 to consecutive spatial positions in row-major order (ow varies fastest):
y=0: (od=0, oh=0, ow=0)
y=1: (od=0, oh=0, ow=1)
...
y=31: (od=0, oh=0, ow=31)
Then nearby y values have the same od and oh, but different ow. For these, the iw values are nearby (e.g., ow=0→iw∈[0,1], ow=1→iw∈[0,1,2], so some overlap).
But the input loads for different y values are to different iw's, so not perfectly coalesced. However, nearby iw values might fall in the same cache line.
Actually, for consecutive ow values, the iw values overlap significantly. For example:
ow=0: iw∈[0,1] (2 values)
ow=1: iw∈[0,1,2] (3 values)
ow=2: iw∈[0,1,2] (3 values)
ow=3: iw∈[0,1,2,3] (4 values)
For a fixed combination (say the first valid iw for each ow):
ow=0, iw=0
ow=1, iw=0
ow=2, iw=0
ow=3, iw=0
So all 4 threads load input at the same iw! That's great.
But for a different combination within the iteration:
ow=0 has 2 iw's, ow=1 has 3 iw's, etc.
When we loop over all valid (id, ih, iw) combinations for each spatial position, the exact combinations differ, so perfect alignment is hard.
OK, I think I'm overthinking this. Let me just write a reasonable kernel and test it. I'll start with a simple output-per-thread kernel with block=(256) and exploit the fact that we skip zero multiplications.
Let me draft the CUDA kernel:
```cuda
__global__ void conv_transpose3d_kernel(
const float* __restrict__ input,
const float* __restrict__ weight,
const float* __restrict__ bias,
float* __restrict__ output,
int batch_size, int in_channels, int out_channels,
int depth, int height, int width,
int out_depth, int out_height, int out_width,
int kernel_d, int kernel_h, int kernel_w,
int stride_d, int stride_h, int stride_w,
int padding_d, int padding_h, int padding_w,
int output_padding_d, int output_padding_h, int output_padding_w,
int groups
) {
int out_ch_per_group = out_channels / groups;
int in_ch_per_group = in_channels / groups;
int total_out = batch_size * out_channels * out_depth * out_height * out_width;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= total_out) return;
// Decode idx to (b, oc, od, oh, ow)
int tmp = idx;
int ow = tmp % out_width; tmp /= out_width;
int oh = tmp % out_height; tmp /= out_height;
int od = tmp % out_depth; tmp /= out_depth;
int oc = tmp % out_channels; tmp /= out_channels;
int b = tmp;
int g = oc / out_ch_per_group;
int oc_local = oc % out_ch_per_group;
float acc = 0.0f;
// For each valid input depth
int id_min = max(0, (od + padding_d - kernel_d + 1 + stride_d - 1) / stride_d);
int id_max = min(depth - 1, (od + padding_d) / stride_d);
for (int id = id_min; id <= id_max; ++id) {
int kd = od + padding_d - id * stride_d;
if (kd < 0 || kd >= kernel_d) continue;
for (int ih = max(0, (oh + padding_h - kernel_h + 1 + stride_h - 1) / stride_h);
ih <= min(height - 1, (oh + padding_h) / stride_h); ++ih) {
int kh = oh + padding_h - ih * stride_h;
if (kh < 0 || kh >= kernel_h) continue;
for (int iw = max(0, (ow + padding_w - kernel_w + 1 + stride_w - 1) / stride_w);
iw <= min(width - 1, (ow + padding_w) / stride_w); ++iw) {
int kw = ow + padding_w - iw * stride_w;
if (kw < 0 || kw >= kernel_w) continue;
for (int ic_local = 0; ic_local < in_ch_per_group; ++ic_local) {
int ic = g * in_ch_per_group + ic_local;
float in_val = input[((b * in_channels + ic) * depth + id) * height + ih) * width + iw];
float w_val = weight[(((ic * out_ch_per_group + oc_local) * kernel_d + kd) * kernel_h + kh) * kernel_w + kw];
acc += in_val * w_val;
}
}
}
}
if (bias != nullptr) {
acc += bias[oc];
}
output[((b * out_channels + oc) * out_depth + od) * out_height + oh) * out_width + ow] = acc;
}
```
Wait, I need to be careful about integer division rounding. The formula `ceil(a/b)` for positive integers is `(a + b - 1) / b`.
In my code:
id_min = max(0, ceil((od + padding_d - kernel_d + 1) / stride_d))
id_max = min(depth - 1, floor((od + padding_d) / stride_d))
For positive values, floor division is just integer division `/` in C++ for positive numbers.
Ceil division is `(numerator + denominator - 1) / denominator` for positive numerator.
But `od + padding_d - kernel_d + 1` can be negative! If it's negative, ceil(negative / positive) might be negative, and we clip to 0 anyway.
So:
id_min = max(0, (od + padding_d - kernel_d + 1 + stride_d - 1) / stride_d)
= max(0, (od + padding_d - kernel_d + stride_d) / stride_d)
Let me verify with the example: od=0, p=1, K=3, s=2:
id_min = max(0, (0 + 1 - 3 + 2) / 2) = max(0, 0 / 2) = 0 ✓
od=0 requires id=0, kd=1.
id_max = min(D-1, (0 + 1) / 2) = min(2, 0) = 0
So id=0 is the only value.
kd = 0 + 1 - 0*2 = 1. ✓
od=1:
id_min = max(0, (1 + 1 - 3 + 2) / 2) = max(0, 1 / 2) = 0
id_max = min(2, (1+1)/2) = min(2, 1) = 1
id=0: kd = 1+1-0 = 2 ✓
id=1: kd = 1+1-2 = 0 ✓
od=2:
id_min = max(0, (2+1-3+2)/2) = max(0, 2/2) = 1
id_max = min(2, (2+1)/2) = min(2, 1) = 1
id=1: kd = 2+1-2 = 1 ✓
Good. But for od=5 with op=1:
id_min = max(0, (5+1-3+2)/2) = max(0, 5/2) = 2
id_max = min(2, (5+1)/2) = min(2, 3) = 2
id=2: kd = 5+1-4 = 2 ✓
output[5] = c*w2.
Now, what if id_max computes to a value >= depth and id_min also >= depth? Then the loop doesn't execute, and output is 0. That corresponds to output positions where there are no valid contributions.
Wait, but I need to be careful about the output_padding region. For the test, D_out=24 with D_in=12. od=23:
id_min = max(0, (23+1-3+2)/2) = max(0, 23/2) = 11
id_max = min(11, (23+1)/2) = min(11, 12) = 11
id=11: kd = 23+1-22 = 2
So there IS a valid input position for od=23. Good.
For the "extra" output_padding positions, are they all zero?
od=23 is actually the last output position (since D_out=24, indices 0-23). And it has a valid input. So output_padding doesn't create all-zero positions in this case.
OK, the kernel formula looks correct.
Now, about performance. The kernel as written has the loops over id, ih, iw in the inner loops, and ic in the innermost loop. This means:
For each valid spatial combination (id, ih, iw, kd, kh, kw):
For each ic_local:
Load input[b, ic, id, ih, iw]
Load weight[ic, oc_local, kd, kh, kw]
Multiply and accumulate
With 256 threads per block, if all threads are computing the same oc but different spatial positions, the weight loads for the same (ic, oc_local, kd, kh, kw) differ between threads. That's not ideal.
If all threads compute different oc but the same spatial position, the input loads are the same (broadcast), and weight loads for the same (ic, kd, kh, kw) but different oc_local are contiguous. That's better.
But as discussed, having different oc in the same warp is problematic because they might be in different groups.
What if we restructure the kernel to process one group at a time, and within a group, have contiguous threads handle contiguous spatial positions and contiguous output channels?
Let me try: each block handles 1 batch, 1 group (8 in-channels, 8 out-channels), and some number of spatial positions.
Block = (8, 32) → 256 threads.
- x: 8 output channels (all channels in the group)
- y: 32 spatial positions
All threads in the block are in the same group. Grid:
- x: not used (or 1)
- y: not used (or 1)
- z: spatial blocks
Actually, let me change the grid layout:
Grid = (spatial_blocks, groups * batch_size)
Block = (8, 32)
Each block computes 8 output channels for 32 spatial positions, for 1 group and 1 batch.
For spatial decomposition: total spatial positions = 24*48*96 = 110,592.
With 32 positions per block: 110,592 / 32 = 3456 spatial blocks.
Total blocks: 3456 * (4 * 8) = 3456 * 32 = 110,592 blocks.
Each thread computes 1 output element.
Within a warp (32 threads): 8 channels for 4 spatial positions.
For warp 0: spatial positions 0, 1, 2, 3 (in some ordering).
For each valid spatial combination, 8 threads load contiguous weights for the 8 channels.
And 4 groups of 8 threads (one per spatial position) load input from different spatial positions.
The input loads might not be coalesced, but they can be cached.
The weight loads are coalesced within each sub-group of 8 threads.
This might be good enough. Let me also think about how to improve input cache hits.
If the 4 spatial positions in a warp are consecutive in some dimension where input overlaps, we get better cache reuse.
For example, if consecutive spatial positions differ in ow by 1, they share many iw values. So processing spatial positions in row-major order (ow fastest) is good.
Let me assign the 32 spatial positions in a block as a 1D tile in row-major order. For example, the block handles a region like:
(od_start, oh_start, ow_start) to (od_start, oh_start, ow_start+31) (wrapping as needed).
Or a 2x2x8 tile, etc. A 1x1x32 tile (varying ow) would maximize w-dimension overlap.
Actually, let me think about what happens when we process 32 consecutive ow values with the same od and oh.
For ow=0 to 31:
The valid iw ranges overlap significantly. iw typically ranges from 0 to about 16 for these ow values.
But different threads in the warp are processing different ow values, so their iw loops have different bounds and iterations. This causes divergence.
If I use a tile of 2x4x4 for example:
od∈{0,1}, oh∈{0,1,2,3}, ow∈{0,1,2,3}
For d even (0): 1 id. For d odd (1): 2 ids. So there's some divergence.
For h: oh=0→2, oh=1→2, oh=2→3, oh=3→2.
For w: ow=0→2, ow=1→3, ow=2→3, ow=3→4.
The loop counts vary, but not wildly. And nearby positions share some input values.
I think for the first version, I'll just use a simple 1D spatial decomposition where each block gets a contiguous chunk of spatial positions. The exact details might not matter too much if the kernel is memory/compute efficient.
Actually, there's another approach that might be even simpler and more efficient: unroll the spatial combinations explicitly.
Since the kernel sizes are small (3, 5, 7) and the number of valid combinations is small (at most 24), we can precompute which (kd, kh, kw) are valid for each (od, oh, ow) phase.
But this would require a lot of precomputation. For a simpler first attempt, I'll stick with the loop-over-id approach.
One more optimization: I should handle the loop ordering to maximize reuse. The innermost loop should be over ic_local since the input and weights for consecutive ic_local are memory-adjacent or at least have predictable strides.
Wait, input[b, ic, id, ih, iw] for consecutive ic are separated by depth*height*width. That's not contiguous. But with groups=4 and in_ch_per_group=8, we only iterate over 8 channels, which is a small range.
Weight[ic, oc_local, kd, kh, kw] for consecutive ic are separated by out_ch_per_group * kd * kh * kw = 8 * 3 * 5 * 7 = 840 floats.
Neither is particularly contiguous. But if we load weight for a fixed ic and all oc_local (8 values), that's a contiguous vector of 8 floats. Good for float4 or float2 loads.
Similarly, if we load input for a fixed (id, ih, iw) and all ic_local (8 values), that's 8 floats separated by D*H*W. Not contiguous at all.
So it makes sense to have the oc_local loop inside the ic_local loop, so that for each ic_local, we load the input once and then use it for all oc_local? No, that would mean loading 8 weights per input value, but the weights are contiguous.
Actually, let me think about this differently. For a fixed spatial combination:
For ic_local = 0:
Load input_0
Load weight_0[0:7] (8 floats, contiguous)
For oc_local = 0..7: acc[oc_local] += input_0 * weight_0[oc_local]
For ic_local = 1:
Load input_1
Load weight_1[0:7]
For oc_local = 0..7: acc[oc_local] += input_1 * weight_1[oc_local]
This processes one ic_local at a time, loading one input value and one weight vector of 8. The weight vector is contiguous and can be loaded with two float4 instructions.
If a thread computes all 8 output channels, it needs 8 accumulators and does 8 multiply-adds per ic_local. That's 8*8 = 64 operations per spatial combination.
If a thread computes only 1 output channel, it loads 1 weight per ic_local, but the weights for different threads (different oc_local) are contiguous. A warp of 8 threads computing 8 channels could load all 8 weights with two float4 loads and broadcast the input value.
The second approach (one thread per channel, cooperative loading) might have better register usage but requires warp-level cooperation.
For simplicity, let me have each thread compute all 8 channels for its spatial position. This requires 8 accumulators (can fit in registers) and loads an 8-float weight vector per ic_local per spatial combination.
With at most 24 spatial combinations and 8 ic_local values, each thread loads 24*8 = 192 weight vectors of 8 floats. That's 1536 floats per thread. In 32 registers (512 floats), we can't store all weights, so they must be loaded from global memory.
But if the weights fit in L1 cache or are read repeatedly by many threads, the effective bandwidth is lower.
Actually, all threads in the block that are in the same group load the same weight values (since they all compute channels in the same group, just different spatial positions). So with 256 threads per block and 32 spatial positions, each weight is loaded 32 times (once per spatial position).
With L1 cache of 128 KB per SM on H100, the weights for one group are 27 KB. This easily fits in L1 cache. So after the first few reads, weights are served from L1 cache at very high bandwidth.
Similarly, input values for different spatial positions overlap, so L1 cache helps there too.
So the memory bandwidth should not be a bottleneck if L1 caching works well.
OK, I think the kernel design is reasonable. Let me write a concrete CUDA kernel and test it.
One more thought: I can also try using shared memory for weights. If I load all weights for the group (27 KB) into shared memory once per block, all threads can access them at shared memory speed. This avoids L1 cache misses and contention.
For 256 threads, loading 27 KB of weights into shared memory is quite fast. And 27 KB fits in the 227 KB limit.
With shared memory for weights:
1. At the start of the block, all threads cooperatively load weight[g*8:(g+1)*8, 0:8, 0:3, 0:5, 0:7] into shared memory.
2. Each thread then accesses weights from shared memory.
This should be faster than loading from global memory, especially since the weight accesses are random (different kd, kh, kw for different threads).
Similarly, I could also load input into shared memory. But since different spatial positions need different input regions, and the total input needed for 32 spatial positions might be large, this is trickier.
Let me estimate: for 32 spatial positions in a tile, if they're all with the same od and oh but different ow (e.g., ow=0 to 31), the iw range is roughly 0 to 17 (since ow=31→iw up to floor((31+3)/2)=17). And ih is roughly 0 to 13 or so (depending on oh).
So input for one batch, one group (8 channels), id range (say 0-1 or similar), ih range (~14 positions), iw range (~18 positions):
8 * 2 * 14 * 18 ≈ 4032 floats ≈ 16 KB.
This also fits in shared memory!
So we could load both weights (27 KB) and relevant input (16 KB) into shared memory, totaling 43 KB per block. Still well within limits.
But the exact input region depends on the spatial tile. And loading input into shared memory requires knowing which (id, ih, iw) are needed and having all threads participate.
Actually, this is getting complex. For the first iteration, let me just write a kernel with shared memory for weights, and see if it compiles and passes correctness. Then I can optimize further.
Let me draft the final kernel code. I'll write it carefully.
Actually, before I commit to a complex kernel, let me just try a simple version first to establish a baseline.
Simple version: each thread computes one output element, with all loops in the kernel and global memory loads. No shared memory.
Then I can measure and see if it's correct. If the speedup is already good, great. If not, I'll add shared memory.
I'll write the simple version first.