DBMSolver: A Training-free Diffusion Bridge Sampler for High-Quality Image-to-Image Translation

A pretrained image-to-image diffusion model, sampled in 6 network calls instead of 119, with no retraining.

A diffusion bridge turns one image into another by running a noisy process that is pinned to the input at one end and to the output at the other. The equation that runs it backward is linear in the image and nonlinear only in the network's guess of the output image. DBMSolver solves the linear part exactly and approximates only that guess.

Hand a diffusion bridge model the edge sketch of a handbag and it paints the handbag. With the sampler it shipped with, the painting takes 119 passes through a 284-million-parameter network and scores an FID of 1.83. DBMSolver runs the same weights through six passes and scores 0.93.

Translation as a walk between two images

Image-to-image translation covers every task where the input is a picture and the output is another picture of the same scene: an edge map into a photograph (Edges2Handbags), a surface-normal map into an outdoor scene (DIODE), a face into its comic version, a segmentation mask into a face, a photo with a hole into the photo filled in. The paper studies all five. In each case the model sees the input and must produce one plausible output, and for most inputs there are many plausible outputs, so the model has to be generative rather than a fixed function.

An ordinary diffusion model is a poor fit for translation. It is built to start from pure Gaussian noise and walk toward data, so the input image has to be bolted on as extra conditioning while the walk still begins at noise. A diffusion bridge model (DBM), introduced in Denoising Diffusion Bridge Models (DDBM) by Zhou, Lou, Khanna and Ermon in 2023, changes the walk itself: the noisy process runs from the output image to the input image, and the network learns to walk it backward, from input to output.

The cost is counted the same way as for any diffusion model. Each step of the trajectory needs the network to look at the current image and make a prediction, and one such forward pass is one network function evaluation, NFE for short. NFE is the cost that matters because everything else in the sampler is arithmetic on tensors, and at 256×256 with a 534-million-parameter U-Net nearly all of the wall-clock time is spent inside the network. The sampler that DBMs shipped with, Hybrid Heun, takes 40 steps at up to 3 calls each, 119 NFE in total. A 2024 follow-up, DBIM, brought that to 20. DBMSolver brings it to 6 while improving the output, and it does so without touching the network: no distillation, no fine-tuning, no new training run. The Edges2Handbags and DIODE numbers in the paper come from the exact checkpoints the DDBM authors released, which is the sense in which the sampler is "training-free".

The argument runs in a few steps: what a bridge is and why it is pinned at both ends; the two equations that run it backward; the observation that those equations are linear in the image and nonlinear only in the network's prediction; the exact solution that follows from solving the linear part by hand, and the second-order approximation of what remains; the singular endpoint that forces the first step to be random; and the call-by-call algorithm with its budget of $2N-2$ network calls for $N$ steps.

A diffusion pinned at both ends

Start with an ordinary forward diffusion, the kind covered in the DDPM and score-based SDE explainers. A clean image $\mathbf{x}_0$ is scaled down and has Gaussian noise added, more of both as time $t$ runs from 0 to $T$:

Two summaries of the schedule appear everywhere below. The signal-to-noise ratio $\text{SNR}_t = \alpha_t^2/\sigma_t^2$ says how much of $\mathbf{x}_t$ is still image, and $\lambda_t = \log(\alpha_t/\sigma_t)$, half the log of the SNR, is the variable the solvers below integrate in: equal steps in $\lambda$ are wide in the early, high-noise part of the process and narrow in the late, nearly clean part. The bridge models in this paper use the "variance-preserving" schedule from the DDBM code, $\alpha_t = e^{-(0.05t + 0.5t^2)}$ and $\sigma_t^2 = 1-\alpha_t^2$ on $t\in[0,1]$. Its constants differ from the VP schedule of Song et al., which ends at nearly pure noise. Here, at $t=T=1$, the signal factor is $\alpha_T = 0.577$ and the SNR is 0.499: the end of this schedule is half image, half noise. A bridge does not need to reach noise, because its endpoint is going to be an image anyway. The half log-SNR runs from $\lambda_T = -0.35$ at the input end to $\lambda = 5.76$ at $t = 10^{-4}$, the smallest time the sampler visits.

To pin the process, take that forward diffusion and keep only the realizations that happen to arrive at a chosen point $\mathbf{x}_T$ at time $T$. Conditioning a diffusion on where it ends is a classical operation called Doob's h-transform, and its effect is a single extra drift term: at every time the process is pushed in the direction that makes arriving at $\mathbf{x}_T$ more likely, with strength $g(t)^2\,\nabla_{\mathbf{x}_t}\log p(\mathbf{x}_T\mid\mathbf{x}_t)$, where $g(t)$ is the diffusion coefficient (the rate at which noise is injected) and the gradient is of the log-probability of reaching the destination from the current point. Because the unconditioned process is Gaussian, that transition probability is Gaussian and its score has a closed form, the paper's equation (4):

Read it as a spring: it pulls $\mathbf{x}_t$ toward a rescaled copy of the destination, and the denominator shrinks to zero as $t\to T$, so the pull becomes infinitely stiff at the end. The stiffness guarantees arrival, and the same denominator makes the exact ODE solution singular at its starting point, the subject of the section on the first step. No network is involved in this term; it is a closed-form function of the schedule.

The process that results is a diffusion bridge: it leaves $\mathbf{x}_0$ at time 0 and arrives at $\mathbf{x}_T$ at time $T$ with certainty. Its state at an intermediate time is Gaussian, and the DDBM paper works out the mean and variance. With the shorthand $r_t = \text{SNR}_T/\text{SNR}_t$, a number between 0 and 1 that measures how far along the bridge the state is:

The mean is a blend of the two endpoints and the variance is the ordinary diffusion variance $\sigma_t^2$ shrunk by the factor $1-r_t$. At $t=0$ the blend is all $\mathbf{x}_0$ and the variance is zero because $\sigma_0 = 0$. At $t=T$ the ratio $r_T$ is 1, so the blend is all $\mathbf{x}_T$ and the variance is zero again. In between the band opens and closes. With the schedule above, at $t=0.5$ the weights are $a_t=0.260$ on the input and $b_t=0.711$ on the output, with standard deviation $c_t = 0.463$ on the $[-1,1]$ pixel scale. (The paper labels this Gaussian's score as the score of $p_t(\mathbf{x}_t\mid\mathbf{x}_T)$; the expression contains $\mathbf{x}_0$, so it is the score of the two-endpoint conditional written here. The distinction comes up again two paragraphs down.)

Drag the probe in the figure and read the three coefficients off the band. Then switch to the "variance-exploding" schedule the paper also supports, where $\alpha_t = 1$ and $\sigma_t = t$ up to $T=80$: the band changes shape but stays pinned at both ends.

Every sample path arrives at $\mathbf{x}_T$ itself, so when the process is run backward it can start from the input image itself rather than from noise. A process with zero variance at its start also has no room to differ from one run to the next, and a deterministic solver started there has nothing to fan out from; the section on the first step returns to this.

DDBM trains a network $\mathbf{D}_{\theta}(\mathbf{x}_t, t, \mathbf{x}_T)$ to predict the clean output $\mathbf{x}_0$ from a bridge sample: draw a training pair, draw a time, sample $\mathbf{x}_t$ directly from the Gaussian in (3) (no simulation of the process is ever needed), feed the network $\mathbf{x}_t$ with the input image concatenated along the channel axis, and regress onto $\mathbf{x}_0$ with a squared error. For DIODE at 256×256, $\mathbf{x}_t$ is a $3\times256\times256$ tensor of pixel values in $[-1,1]$, the U-Net receives six channels (the noisy image scaled by a preconditioning factor, then the raw input image), takes $250\ln t$ as its time embedding, and returns a $3\times256\times256$ prediction of the output. The paper calls this an $\mathbf{x}_0$-predicting DBM and writes it $\mathbf{D}_\theta(\mathbf{x}_s)$ for short. Predicting the clean image rather than the noise or the score is a choice inherited from DPM-Solver++: the clean image is bounded, so the weights the solver puts on it stay bounded even when a step covers a large range of noise levels, whereas weights on a noise prediction grow without limit toward the clean end.

The solver, though, is written in terms of a score, the gradient of the log-density of the bridge state given only the input, $\nabla\log p_t(\mathbf{x}_t\mid\mathbf{x}_T)$. The two are connected by the bridge version of Tweedie's formula. The score of the two-endpoint Gaussian (3) is $(a_t\mathbf{x}_T + b_t\mathbf{x}_0 - \mathbf{x}_t)/c_t^2$, which is affine in $\mathbf{x}_0$. The score given only $\mathbf{x}_T$ is the average of that expression over every $\mathbf{x}_0$ that could have produced this $\mathbf{x}_t$, and averaging an affine function of $\mathbf{x}_0$ replaces $\mathbf{x}_0$ by its posterior mean. The ideal network is that posterior mean, so:

The identity is exact when $\mathbf{D}_\theta$ equals the posterior mean, and off by the network's error otherwise. The two-endpoint Gaussian is the regression target, the one-endpoint score is what the trained network ends up approximating, and (6) converts between them. The $c_t^2$ in the denominator vanishes at both ends of the bridge, so any error in $\mathbf{D}_\theta$ is magnified by $1/c_t^2$ when converted to a score. The sampler below therefore carries the prediction $\mathbf{D}_\theta$ itself through every step and converts to a score once, in the final Euler step at $t=10^{-4}$.

Two ways to run the bridge backward

A forward diffusion can be reversed in two equivalent ways, and both carry over to the bridge. The first is a stochastic differential equation that injects fresh noise as it goes; the second is an ordinary differential equation with no noise at all, and the two produce the same distribution of images at every time even though any single run of one looks nothing like a single run of the other. For the bridge, with the drift $f(\mathbf{x}_t,t)$ and diffusion $g(t)$ of the underlying process, the paper's equations (2) and (5) are:

Both run from $t=T$ down to $t=0$, so the time increments are negative. Each has three drift pieces. The base drift $f$, for this schedule $f = -(0.05 + t)\,\mathbf{x}_t$, is known. The learned piece is the score of the bridge state, supplied through (6) by the network. The third piece is the h-transform drift toward the input from (4), known in closed form with no network. The SDE version, called the Bridge SDE, adds noise $g(t)\,\mathrm{d}\mathbf{w}_t$; the ODE version, the Bridge probability-flow ODE, has none.

Compare the two carefully: the score term carries a full $g^2$ in the SDE and half of it in the ODE, while the h-transform term carries $g^2$ in both. The half is not a convention. The noise in the SDE spreads the distribution, and to reproduce the same spreading with a noiseless flow you must convert the diffusion into a transport term, which costs half the score (the argument from Song et al. that gives every probability-flow ODE its $\tfrac12$). The h-transform term is part of the forward drift that defines the bridge, so the conversion leaves it untouched.

The fast samplers built for ordinary diffusion, the ones behind DDIM and DPM-Solver++, cannot be dropped onto a bridge model for the same reason. Their update rules are exact solutions of an ODE that has no h-transform term and a different linear coefficient. Fed a bridge network, they integrate the wrong equation: nothing steers the trajectory to $\mathbf{x}_T$, and the network is handed states it never saw during training. The paper describes the failure as those samplers assuming a Gaussian prior; the update formulas encode the wrong forward process, so handing them the input image as a starting point does not repair them. On CelebAMask-HQ at 6 NFE, DPM-Solver++ scores an FID of 169 against DBMSolver's 35.

Linear in the image, nonlinear only in the network

Substitute the closed forms (4) and (6) into the ODE (5) and collect terms in $\mathbf{x}_t$. Every piece except the network's prediction is a known scalar function of time multiplying $\mathbf{x}_t$ or $\mathbf{x}_T$, and the paper's equation (7) records the shape:

An equation of this form is called semi-linear: linear in the state, with a nonlinear remainder. The remainder here is affine in the network output, plus a known multiple of $\mathbf{x}_T$; all the nonlinearity lives inside $\mathbf{D}_\theta$, which depends on $\mathbf{x}_t$ in whatever way a U-Net does. For the VP schedule the linear coefficient works out to $L(t) = \tfrac{\mathrm{d}}{\mathrm{d}t}\log\sigma_t + \lambda_t'/\rho_t$ with $\rho_t = \text{SNR}_t/\text{SNR}_T - 1$, a quantity that vanishes at $t=T$, so $L$ blows up at the input end. DPM-Solver exploited the same structure for ordinary diffusion in 2022, and DBIM wrote it down for bridges in 2024; the next section covers what this paper adds.

Why does the split matter for a solver? Euler's method advances the state by the right-hand side times the step, treating the entire right-hand side as constant over the step. That is a poor approximation of the linear part, which is known in closed form and which, for this equation, is stiff: $L(t)$ contains the $1/c_t^2$ factors from (4) and (6), so it is enormous near both ends of the bridge, and a plain Euler step of any useful size overshoots and diverges. An exponential integrator instead treats $L(t)\mathbf{x}_t$ as the equation and $N$ as a forcing term, and the variation-of-constants formula gives the state at the end of a step from $s$ to $t$ exactly:

The first term propagates the current state with an exponential factor that is a closed-form number for this schedule, so no step-size error is ever attached to the linear part. The only thing left to approximate is how the network's prediction varies along the step inside the integral. Freeze it at its value at the start of the step and you have a first-order method that is exact whenever the prediction happens to be constant. Estimate its slope as well and you have a second-order method, with no error term that depends on the size of $L$.

The toy below is the same situation in one dimension: $\mathrm{d}x/\mathrm{d}t = -k\,x + N(t)$ with a stiff decay $k$ standing in for $L(t)$ and a smooth forcing $N(t)$ standing in for the network term. Watch what happens to Euler as you push $k$ past $2$ divided by the step, and what does not happen to the exponential integrators:

The paper uses "order" in the numerical-analysis sense. A method is k-th order when the error it commits in one step shrinks like the step size to the power $k+1$, so that after the $1/h$ steps needed to cross a fixed interval the accumulated error shrinks like $h^k$. Halve the step and a first-order method halves its final error; a second-order method quarters it. The paper's "error $O(\Delta t^2)$" for its first-order step is the per-step number; the error in the final image is one power lower. In the toy that is visible directly, and in the two-dimensional race later on the page it is measured: ratios of 2 per halving for the first-order step, ratios approaching 4 for the second.

The integrals above are taken in $\lambda$, the half log-SNR, rather than in $t$, a choice DPM-Solver introduced and DBIM and this paper keep. Both the propagator $e^{\int L}$ and the weight inside the integral become clean exponentials of $\lambda$, and the network's prediction is Taylor-expanded in $\lambda$ too. Since the network is fed $250\ln t$ as its time input, a linear change in $\lambda$ is closer to a linear change in what the network sees than a linear change in $t$ would be.

The exact solution and its second-order step

Carrying the variation-of-constants formula through for the bridge ODE, with the propagator worked out in closed form, gives the paper's Proposition 2. Write $\rho_t$ for $e^{2(\lambda_t-\lambda_T)}-1$, which is $\text{SNR}_t/\text{SNR}_T - 1$: zero at the input end, growing as the state moves toward the output. Then for a step from $s$ to $t<s<T$:

Each of the three terms has one job. The first carries the current state forward with a known weight. The second is a known weight on the input image. It exists because the score term enters the ODE at half strength, so it cancels only half of the h-transform's pull toward the input, whereas in the SDE the two cancel exactly; it has no counterpart in ordinary diffusion solvers. The third is an integral over the step of the network's prediction, weighted by $e^{2\lambda}/\sqrt{\rho_\lambda}$, and it is the only place approximation enters.

The approximation is a Taylor expansion of $\mathbf{D}_\theta$ in $\lambda$ around the start of the step, the paper's equation (10):

The derivatives $\mathbf{D}_\theta^{(n)}$ have to be estimated from network calls; the weight integrals depend only on the schedule and can be worked out once by hand. With $k=1$ the prediction is frozen, the single weight integral is $e^{2\lambda_T}\sqrt{\rho_t}\,(1-\sqrt{\rho_s/\rho_t})$, and the step becomes a fixed linear combination of the current state, the input, and one network prediction (the supplement's equation 21). That step is, coefficient for coefficient, the deterministic step of DBIM, which the paper acknowledges. The listing below gives the three weights in code and their values for the first step of a 6-NFE run, where a single step carries the state from $t=0.9999$ to $t=0.0412$, a span of 2.92 in $\lambda$:

# the three weights of the first-order step (supp. Eq. 21), s -> t < s < T
rs, rt = rho(lam(s)), rho(lam(t))
c_x = alpha(t)/alpha(s) * exp(2*(lam(s)-lam(t))) * sqrt(rt/rs)   # on x_s
c_T = alpha(t)/alpha(T) * exp(2*(lam(T)-lam(t))) * (1 - sqrt(rt/rs))  # on x_T
c_D = alpha(t) * exp(2*(lam(T)-lam(t))) * rt * (1 - sqrt(rs/rt))      # on D(x_s)
x_t = c_x * x_s + c_T * x_T + c_D * D_s

# 6-NFE run on DIODE, first ODE step  s = 0.9999 -> t = 0.0412:
#   c_x = 5.251   c_T = -5.245   c_D = 0.993
# so  x_t = 0.993 D_s + 5.25 (x_s - x_T) + 0.006 x_T   (to three decimals)

After the first ODE step the state is $0.993$ times the network's guess of the output, plus $5.25$ times the displacement of the current state from the input, plus a trace of the input. The guess dominates because $b_t = 0.994$ at $t = 0.041$: the bridge state there is almost entirely output. The section on the first step shows where the $5.25$ comes from.

With $k=2$ the expansion keeps the slope of the prediction along $\lambda$. The new weight integral has a closed form with an arctangent in it (supplement equation 24), and the full step, supplement equation 26, is the first-order step plus one correction:

The slope $\mathbf{D}_\theta^{(1)}$ has to be estimated from values, since the network returns only predictions. DBMSolver does it the way DPM-Solver's second-order singlestep method does: take the first-order step only as far as the midpoint of the interval in $\lambda$, call the network there, and divide the difference of the two predictions by the gap in $\lambda$. That midpoint call is the second network evaluation each step costs, and it is fresh information about the current step rather than a value remembered from the previous one. For the first ODE step of the 6-NFE run the midpoint sits at $t = 0.274$, the slope's weight comes out to $1.99$, and the finite difference divides by a gap of $1.46$. The supplement leaves the position of the midpoint unstated; the released code puts it at the $\lambda$-midpoint for DIODE and at the $t$-midpoint for Edges2Handbags, a per-dataset choice the paper does not mention. Both are valid slope estimates because the difference is always divided by the actual $\lambda$ gap, but they are different samplers.

The paper stops at $k=2$ because a third-order step would need the weight integral with $(\lambda-\lambda_s)^2$, and in this variable that integral has no elementary antiderivative: substituting $u=\sqrt{\rho_\lambda}$ turns it into an integral of $\log^2(1+u^2)$, which needs the dilogarithm. DBIM reaches second and third order a different way, by reusing the network's predictions from earlier steps in a multistep formula (one call per step, but extrapolating from stale information). The obstruction belongs to the expansion variable rather than to the bridge: DBIM expands in $\log(b_t/c_t)$ instead of $\lambda$, and in that variable every weight is elementary. DBIM's second-order slope is also a derivative with respect to that different variable, so the two second-order methods are different schemes rather than the same scheme with and without memory.

The Bridge SDE (2) has the same structure and yields to the same treatment. Its linear coefficient integrates to $(\alpha_t/\alpha_s)\,e^{2(\lambda_s-\lambda_t)}$, and after substituting (4) and (6) the terms in $\mathbf{x}_T$ cancel identically, so the SDE drift is linear in $\mathbf{x}_t$ plus a multiple of the network's prediction, the same shape as for an ordinary diffusion. Freezing the prediction and integrating the noise term exactly gives the paper's Proposition 1, its first-order SDE step:

With the prediction held fixed this is a draw from the Gaussian posterior of the forward process given the current state and a clean image equal to the prediction, and when the step starts at $s=T$ its three coefficients are the bridge coefficients $a_t, b_t, c_t$ of (3). The supplement also derives a second-order version of (8), and the coefficient it prints for the derivative term does not match the integral it comes from. Evaluating $2\alpha_t e^{-2\lambda_t}\int_{\lambda_s}^{\lambda_t} e^{2\lambda}(\lambda-\lambda_s)\,\mathrm{d}\lambda$ gives the first expression below; the supplement prints the second:

Expand both in the step $h=\lambda_t-\lambda_s$: the correct one begins at $h^2$, as a second-order correction must, while the printed one begins at $-h$ and would undo the method's order. The released code does not use the printed formula, and the paper's own default is the first-order SDE step, so nothing in the reported results depends on it. The next section covers why the SDE step exists at all.

Why the first step has to be random

Proposition 2 was stated for $s<T$, and the coefficients show why. At $s = T$ the quantity $\rho_s$ is zero, so $\sqrt{\rho_t/\rho_s}$ is infinite: the weight $c_x$ on the current state and the weight $c_T$ on the input both diverge, with opposite signs. The exact ODE solution cannot be started from the input image. In the geometry of Figure 1, a noiseless flow started at the pinched right end stays a single point, so it cannot produce the spread of outputs that the bridge has at every later time; the variation has to be injected before the ODE takes over.

The paper's answer is to take the first step, from $T$ to $T-\epsilon$ with $\epsilon = 10^{-4}$, with the SDE step (8), and every step after it with the ODE step. From $s=T$ the coefficients of (8) are $0.9998$ on the input, $0.00018$ on the network's guess, and $0.0145$ on the noise: the state moves by about one and a half percent of the pixel range. The very next step multiplies the displacement of the state from the input by $c_x = 5.25$, so a displacement with standard deviation $0.0145$ becomes one with standard deviation $0.076$, which is $\sigma_{0.041}$, the noise level of the bridge at the landing time, to three decimals. $c_x$ grows like one over the square root of the gap $T-s$, the noise standard deviation shrinks like the square root of the same gap, and their product is a schedule-determined constant.

So the seed of that one small draw decides which of the possible outputs the deterministic trajectory lands on. The supplement shows several faces generated from one segmentation map by changing only that seed. Check the balance in the figure: the two divergent weights, their finite sum, the shrinking noise, and the flat product, as the gap between the start and the pinned endpoint goes from one half down to one millionth:

Two consequences follow. Because only the first step is stochastic, DBMSolver has no knob for how random the sampler is; DBIM's $\eta$ parameter, which mixes noise into every step, has no counterpart here. And because the SDE step from $T$ is, with the prediction frozen, an exact draw from the bridge Gaussian around the network's guess, the same start is used by DBIM (its "booting" step) and by ODES3, the stochastic-start Heun sampler the paper compares against, whose own theorem shows that this draw beats a single Euler–Maruyama step of the SDE (Euler's method with a noise term added per step).

The algorithm, call by call

Algorithm 1 of the paper, with the details the released code settles, goes like this: pick $N$ step times running down from $t_N = T$ to $t_0 = 0$, with $t_{N-1} = T - 10^{-4}$ and the interior times on the Karras spacing (a power-law grid with exponent 7 that packs points near the clean end) down to $t_1 = 10^{-4}$; call the network once at the input to get $\mathbf{D}_\theta(\mathbf{x}_T)$ and take the SDE step (8) to $t_{N-1}$; take second-order ODE steps (S26) down to $t_1$, two calls each; and finish with one plain Euler step of the ODE (5) from $t_1$ to 0, converting the last prediction to a score with (6) and adding the h-transform term (4). The count of network calls is

Six calls means four steps on the grid 1, 0.9999, 0.0412, 0.0001, 0; twenty calls means eleven steps, with the first ODE step running from 0.9999 to 0.552.

Karras spacing puts many points near the clean end in $t$, yet in $\lambda$ the four-step grid is nearly uniform, with the first ODE step covering 2.92 of the 6.1 units of half log-SNR. The second-order correction matters most on that step: the network's guess changes most between the input end and the middle of the bridge, and freezing it across half the range is the error the midpoint call repairs.

Three implementation facts from the code that the paper leaves out. The final Euler step from $t=10^{-4}$ to 0 does not jump to the clean prediction: its coefficient on the move toward the bridge mean is $0.5005$, so it goes halfway, leaving about half of $\sigma_{10^{-4}} = 0.0032$ in residual noise, well below one 8-bit gray level. Every prediction is clamped to $[-1,1]$ before use, which is harmless because a posterior mean of pixels in that range lies in that range. And the network never receives $T$; it receives the scaled noisy image, the input image concatenated as three more channels, and $250\ln t$ as its time embedding. In pseudocode, the sampler is:

# DBMSolver, N steps, VP bridge (beta_min=0.1, beta_d=2, T=1).  NFE = 2N-2
# D(x, t) = network's x0 prediction, conditioned on x_T (concatenated channels)
lam  = lambda t: log(alpha(t) / sigma(t))          # half log-SNR
rho  = lambda l: exp(2 * (l - lam(T))) - 1          # rho(l, lam_T) of Eq. 9
grid = karras_rho7(n=N-1, t_min=1e-4, t_max=T-1e-4) + [0.0]

x = x_T                                             # the input image
D0 = D(x_T, T)                                      # call 1
x = sde_step(x_T, s=T, t=grid[0], D=D0, z=randn())  # Eq. 8: the only randomness
for s, t in zip(grid[:-1], grid[1:]):
    if t == 0:                                      # last step: plain Euler on Eq. 5
        x = x + (t - s) * bridge_ode_rhs(x, s, D(x, s))          # call
        break
    Ds  = D(x, s)                                                # call
    t_m = t_of_lam(0.5 * (lam(s) + lam(t)))         # midpoint in lambda (DIODE)
    x_m = k1_step(x, s, t_m, Ds)                    # supp. Eq. 21, first order
    D1  = (D(x_m, t_m) - Ds) / (lam(t_m) - lam(s))  # call; slope of D along lambda
    x   = k1_step(x, s, t, Ds) + c1(s, t) * D1      # supp. Eq. 26, second order
return clamp(x, -1, 1)

Same budget, three solvers

The claims so far can be checked on a problem small enough to solve to machine precision. Take a two-dimensional "image", a fixed input point $\mathbf{x}_T$, and let the possible outputs be a mixture of three Gaussians. For that problem the ideal predictor $\mathbb{E}[\mathbf{x}_0\mid\mathbf{x}_t,\mathbf{x}_T]$ has a closed form, so every sampler below runs on the same perfect network with the same VP schedule and the same Karras grid as the code, and the only thing that differs between them is the solver. Each seed has an exact answer, obtained by integrating the ODE with fifteen hundred tiny steps, marked by a ring; the segment from the ring to a sampler's endpoint is that seed's error. Hold the NFE budget fixed and switch methods:

Three things to take from the figure. Euler on the raw ODE diverges at every budget shown: its first step starts next to the pinned end, where the linear term is stiffest, and multiplies the small deviation from the input by hundreds instead of the exact factor, so the state leaves the canvas. The paper's Hybrid Heun and DPM-Solver++ numbers at 5 or 6 NFE (FID 442 and 169 on CelebAMask-HQ) are the image-scale version of a plain solver failing near that end. The first-order step is stable but pulls toward the average of the possible outputs when the steps are coarse, because a frozen prediction over a wide step is the network's guess from the noisy end of that step, which is a blend of all the outputs; in images this shows up as the blur and missing texture the paper reports for DBIM at 6 NFE. And the second-order step separates the modes with the same budget. Measured on this toy at equal NFE, its median endpoint error is 1.14 against 1.54 at 6 NFE, 0.49 against 1.02 at 10, 0.18 against 0.49 at 20, and 0.06 against 0.26 at 40, so the gap widens as the budget grows, as a difference of one order predicts.

The toy's predictor is exact, which no real network is. On real images the network's error in $\mathbf{D}_\theta$ is untouched by any solver, and the paper's own explanation for its weakest result, ImageNet inpainting, is that this error dominates there, since a better solver removes discretization error and nothing else.

What it buys on real images

The paper evaluates on five translation tasks with the standard metrics: FID, the Fréchet distance between Inception features of generated and reference images (lower is better, and it compares distributions, not paired images), plus Inception score, LPIPS (a learned perceptual distance between paired images) and mean squared error, and a classifier accuracy for inpainting. The comparison is against the sampler DDBM shipped with (Hybrid Heun, 119 NFE), DBIM at first, second and third order (20 NFE), a stochastic-start Heun sampler (28 NFE), and two distilled one- and two-step models that require training. Switch datasets in the figure and look for where the teal points sit:

On DIODE, the normals-to-outdoor-scene task at 256×256, six calls give FID 3.38 against Hybrid Heun's 4.43 at 119 calls and DBIM-3's 4.23 at 20. Twenty calls give 2.06, and the abstract's "53% drop" is this number against the authors' rerun of DBIM-2 at 4.40 (DBIM's own paper reports 4.34 for that setting, which would make it 52.5%). LPIPS and MSE, the paired metrics that measure faithfulness to the specific target, are tied between the samplers, so the gain is in how natural the set of outputs looks rather than in per-pixel closeness. On Edges2Handbags, six calls give 0.93 (Table 3; the text and Table 6 say 0.97 for the same run) and twenty give 0.53. These handbag FIDs are computed against the training set, following DDBM's protocol, which is why they are so much smaller than the FIDs reported for unconditional ImageNet generation; compare them with each other, not with those.

On Face2Comics at 20 NFE DBMSolver reaches 0.96 where DBIM at the same budget sits between 8.6 and 9.3 and Hybrid Heun at 119 gets 2.36. At 6 NFE it scores 3.04, better than any DBIM row but worse than Hybrid Heun's 2.36, so on Face2Comics DBMSolver matches Hybrid Heun's quality only at 20 calls, a saving of about 6 times rather than the headline 20. On CelebAMask-HQ, the segmentation-mask-to-face task, DBMSolver reaches 17.56 at 20 NFE against DBIM-3's 19.49, and the supplement adds 34.76 at 6 NFE against DBIM's 44.92. The GAN and other-diffusion rows in those tables are copied from the BBDM paper (Brownian Bridge Diffusion Models, an earlier bridge method that appears as a row in the same tables), which used its own reference split and pipeline, so the gap between 0.96 and BBDM's 23.20 is not a like-for-like measurement.

ImageNet inpainting, filling a 128×128 hole in the center of a 256×256 photo, is the counterexample, and the paper labels it as its limitation. At 20 NFE DBMSolver ties DBIM-2 and DBIM-3, FID 4.07 and classifier accuracy 72.0 for all three; at 6 NFE it scores 4.98, better than DBIM at 6 (5.36) but worse than DBIM at 20 (4.13). Table 5's caption claims superiority "across all NFEs", which its own rows do not support. The paper's explanation is the one from the toy above: on a harder, more varied dataset the network's prediction error dominates the discretization error, and a better solver can only fix the latter.

Against the stochastic-start Heun sampler ODES3 the margin is small on Edges2Handbags, 0.53 at 20 NFE against 0.54 at 28, with ODES3 marginally better on Inception score and LPIPS, so the advantage there is 20 calls against 28; on DIODE the gap is 2.06 against 2.29. Against distillation the paper quotes the two distilled models: consistency bridge distillation (CBD) at 2 NFE (1.30 on Edges2Handbags, 3.66 on DIODE) and Inverse Bridge Matching Distillation (IBMD) at 1 NFE (1.26, 4.07), both worse than DBMSolver's six-call numbers. CBD's own paper, Consistency Diffusion Bridge Models, also reports its stronger consistency bridge training (CBT) variant at 2 NFE, 0.80 and 2.93, which beats DBMSolver at 6; and IBMD at 4 NFE reaches 4.03 on ImageNet inpainting. If you can afford a training run, distillation still wins per call. DBMSolver costs nothing to adopt: the abstract's "up to 5×" fewer calls has no table pair behind it, but 3.3× fewer calls than DBIM at better FID on three of the five tasks, and 20× fewer than Hybrid Heun at better FID on four of them, are both in the tables.

For an implementer, the method reduces to a checklist. Take any DDBM checkpoint with its schedule $\alpha_t, \sigma_t$ and its $\mathbf{x}_0$-predicting denoiser. Build a Karras grid of $N-1$ times from $T-10^{-4}$ to $10^{-4}$ and append 0. Call the denoiser at the input image and take one SDE step (8) with a single Gaussian draw. For each interior interval, call the denoiser at the start, take the first-order step (supplement equation 21) to the $\lambda$-midpoint, call it again there, form the finite-difference slope, and apply (S26). Finish with one Euler step of (5) from $10^{-4}$ to 0 and clamp to $[-1,1]$. Four steps cost six calls and, on DIODE, score an FID of 3.38.