Skip to content

BUG: Draw the feature subsets once per evaluation, and speed up the v4 metric - #6

Open
vboussot wants to merge 9 commits into
mainfrom
subset-once-per-evaluation
Open

vboussot wants to merge 9 commits into
mainfrom
subset-once-per-evaluation

Conversation

@vboussot

@vboussot vboussot commented Oct 18, 2026 •

Copy link
Copy Markdown
Collaborator

Fixes how the v4 metric and the ConvexAdam stages draw their random numbers, and makes the v4 metric faster: Static mode about 6x, the online (Jacobian) mode 4x to 50x on a GPU.

Fixes

  • The v4 metric drew its feature subset per point (Static) or per batch (Jacobian), from a generator per work unit. NCC sums its statistics channel by channel over the points, so those sums mixed channels, and the result depended on the number of threads and on the batch size. With 2 of 4 channels, the Jacobian derivative moved by 7 to 15% between batch sizes. The subset is now drawn once per evaluation.
  • The coarse and fine stages called torch::manual_seed(Seed), which reseeds the whole process, a Python caller's generator included. They now use a generator of their own.
  • The losses summed a batch's values in float32 (2e-8 to 6e-8 relative error over 20000 points). They now sum in double.

Speed, Static mode (lung CT test data, 20000 points, a 32-channel model)

main this PR
32 channels, 24 threads 120 ms 20 ms
32 channels, 1 thread 1281 ms 219 ms
8 of 32 channels, 24 threads 54 ms 6-9 ms
8 of 32 channels, 1 thread 514 ms 33 ms
  • Static now evaluates a work unit's points in batches instead of point by point.
  • The vector interpolator computes the B-spline weights once per point for all channels, instead of once per channel. It gives the same result to the bit as BSplineInterpolateImageFunction, which the new interpolator test checks.
  • The moving gradient image, which IMPACT never reads, is no longer built in Initialize().

Speed, online mode (same data and model, 2000 points unless noted, best of 3)

main this PR
L2, GPU 103 ms 13 ms
L2, GPU, 20000 points 430 ms 110 ms
NCC, GPU 782 ms 15 ms
NCC, CPU, 24 threads 901 ms 121 ms
L2, CPU, 24 threads, 8 of 32 channels (best of 10) 83 ms 85 ms
  • On a GPU, the models ran each work unit's points in batches of 64, one work unit after the other. The work units now only collect their points; all the threads then cut the patches, and the models run on batches as large as the device holds. BatchSize now defaults to 0: that measured bound on a GPU, 64 on the CPU.
  • The trial batch that measures the device budget is no larger than the largest batch an evaluation runs. Sized to the whole free memory, it took that memory from the other processes on a shared GPU.
  • NCC only knows its derivative with respect to one point's features once it has seen every point, so it took every channel's derivative instead, with one backward pass per channel. It now takes two passes: the first, without gradient, gives the value and that derivative at every point; the second needs one backward pass per batch, as L2 does.

Each fix comes with a test that fails without it. A new test checks that the GPU gives the CPU's value and derivative, for L2 and NCC; it is skipped without a CUDA device.

The v4 metric drew its feature subset (SubsetFeatures) from a generator
per work unit, reseeded at every evaluation: per point in Static mode,
per batch in Jacobian mode. A point's channels therefore depended on its
rank in its work unit, on the batch size, and on whether a batch was
replayed after running out of device memory. NCC accumulates its
statistics channel by channel over the points, so its sums mixed
channels. With 2 of 4 channels, the Jacobian derivative moved by 7 to
15% between batch sizes.

The subsets are now drawn once per evaluation, before the threads start,
from the seed the planes use, and every point compares the same channels.
Two evaluations of the same parameters still compare the same channels.

The batch-size test now keeps 2 of the 4 channels instead of all of
them, and a new test compares 1 and 4 work units in Static mode. Both
fail without this change.
Evaluate and EvaluateDerivative took both by value, so every call copied
the subset, two or three times per point and per layer in Static mode.
Both stages called torch::manual_seed(Seed), which reseeds every torch
generator of the process: a Python caller's random numbers restarted
from the stage's seed after each run. The fine stage also drew its
sampled points on the device, from the device's global generator, which
the batch-budget probe draws from too, so a probe could shift them.

Each stage now draws from a CPU generator of its own, seeded with Seed,
and moves the draws to the device. A seed gives the same channels and
points on the CPU and on a GPU. The new test checks that running either
stage leaves torch's global generator where it was; it fails without
this change.
The point-wise losses (L1, L2, Cosine, L1Cosine, Dice) summed a batch's
values in float32 before adding the total to their double accumulator,
so the value lost digits with the size of the batch: 2e-8 to 6e-8 of
relative error over 20000 points. The online mode of the v4 metric and
elastix's IMPACT metric both hand the losses whole batches, and finite
differences of such a value are mostly rounding. The batch is now summed
in double. The new test compares a batch of 20000 points with the same
points fed one by one, and fails without this change.
In Static mode the v4 threader let ITK walk the points one by one. Each
point went through ITK's ProcessVirtualPoint, which also computed an
image gradient IMPACT never reads, then through ProcessPoint, which
built a handful of one-row tensors to update the losses.

Static now gathers a work unit's points, as the online mode already
did, and evaluates them in batches: the feature vectors and the
transform Jacobians are written into contiguous tensors, and each loss
is updated once per batch. A batch holds at most about 4 MB of
[points, channels, parameters] Jacobian, so that a dense domain under a
B-spline transform stays bounded. The vector interpolator gains the
buffer variants this needs, and a value-and-derivative evaluation that
shares the B-spline weights between the two.

On the lung CT test data, with 20000 points and 8 of a model's 32
channels, an evaluation takes 11 ms instead of 40 ms on 24 threads, and
136 ms instead of 468 ms on one.
ImageToImageMetricv4 smooths the whole moving image into a gradient
image in Initialize(), for the metrics that read image gradients at
their points. IMPACT compares features: the online mode takes its patch
gradients from the moving interpolator, and Static mode differentiates
the feature maps. The filter is now off by default. On the lung CT,
Initialize() of an online metric takes 0.01 s instead of 0.08 s.
InterpolateVectorImageFunction held one BSplineInterpolateImageFunction
per channel. Every channel of every point converted the point to a
continuous index, computed the support, the weights and the mirrored
indices, and allocated two or three vnl matrices to hold them, although
all of these are the same for every channel.

It now computes each channel's coefficients as that class does, and
evaluates all the selected channels with one support and one set of
weights, in the same order of operations, so the result is the same to
the bit. SetActiveChannels() keeps a copy of the channels
an evaluation compares, the channels of a voxel side by side, so that a
point reads the fewest cache lines whether it compares all of a model's
channels or a few; the v4 threader sets it once per evaluation. The
per-channel images are gone, and a second SetInputImage() no longer
keeps the first image's channels.

The interpolator test now compares values and derivatives with
BSplineInterpolateImageFunction channel by channel: 70 channels, a
region that does not start at zero, anisotropic spacing, a rotated
direction, through the active channels and through another subset. On
32 channels, a value takes 1.2 us instead of 8.4 us, a value and its
derivative 2.5 us instead of 19.5 us.
On a GPU the online mode ran each work unit's points as batches of at
most BatchSize (64 by default), one work unit after the other under the
device lock. The model takes 20000 patches in about 0.08 s, but they
went through it as some 300 batches.

On a GPU the work units now only collect their points. Once they are
done, all the threads cut the patches of a batch, and the models run on
batches as large as the device holds. BatchSize now defaults to 0, which
means that measured bound on a GPU, and 64 on the CPU, where the work
units remain the parallelism and keep evaluating their own points. The
measurement tries no batch larger than the largest an evaluation can
run: a trial batch sized to the whole free memory took it from the other
processes on a shared GPU.

On the lung CT test data, with the 32-channel model, an evaluation of
20000 points takes 112 ms instead of 430 ms on a GPU; with 2000 points,
13 ms instead of 103 ms (L2) and 94 ms instead of 782 ms (NCC). The new
test compares the GPU's value and derivative with the CPU's, and is
skipped without a CUDA device.
NCC compares statistics taken over all the points, so its derivative
with respect to one point's features is only known once every point has
been seen. The online mode took every channel's d(feature)/d(moving
coordinate) instead, from one backward pass per channel: 32 passes for a
32-channel layer.

A global loss now takes two passes. The first, without gradient,
computes the features of every point, from which the loss takes its
value and its derivative with respect to each point's features. The
second seeds one backward pass per batch with that derivative, as for
L2. To see every point at once, the work units now only collect their
points on the CPU too; the threads then evaluate them in shares, 64
points at a time as before.

On the lung CT test data, with 2000 points and the 32-channel model, an
NCC evaluation takes 15 ms instead of 94 ms on a GPU, and 121 ms
instead of 1006 ms on the CPU. L2 is unchanged.
@vboussot vboussot changed the title BUG: Draw the feature subsets once per evaluation, and speed up Static mode BUG: Draw the feature subsets once per evaluation, and speed up the v4 metric Oct 18, 2026

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant