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 model turns an input image into an output image by reversing a noisy process that is pinned to both. DBMSolver solves exactly the part of the reverse equation that does not involve the network, and approximates only the network's contribution, to second order.

Give a diffusion bridge model the edge sketch of a handbag and it paints the handbag. With the sampler it was released 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 6 passes and scores 0.93.

What DBMSolver speeds up

Image-to-image translation covers every task where the input is a picture and the output is another picture of the same scene. The paper studies five: an edge sketch into a handbag photo (Edges2Handbags), a surface-normal map into an outdoor scene (DIODE-Outdoor), a face into its comic version (Face2Comics), a segmentation mask into a face (CelebAMask-HQ), and a photo with a 128×128 hole in the middle into the filled-in photo (ImageNet inpainting). For most inputs many outputs are plausible, so the model has to sample one rather than compute a fixed function.

An ordinary diffusion model, the kind in the DDPM explainer, starts from pure Gaussian noise and walks toward an image. For translation the input image then has to be attached as extra conditioning while the walk still begins at noise. A diffusion bridge model (DBM), introduced as Denoising Diffusion Bridge Models (DDBM) by Zhou, Lou, Khanna and Ermon in 2023, changes the walk itself: a noisy process runs from the output image to the input image, and the network learns to run it backward, from input to output.

Sampling cost is counted in network function evaluations (NFE): one NFE is one forward pass of the network on the current image. Everything else a sampler does is a handful of elementwise operations on tensors, so at 256×256 with a 534-million-parameter U-Net nearly all of the wall-clock time is spent inside those passes. DDBM shipped with Hybrid Heun, which on each interval takes one Euler–Maruyama step of the SDE (one call) and then a Heun step of the ODE in the style of EDM (two calls): 40 steps, 119 NFE. A 2024 follow-up, Diffusion Bridge Implicit Models (DBIM), brought that to 20. DBMSolver brings it to 6, and at 6 calls it beats Hybrid Heun's FID on four of the five tasks.

"Training-free" means the network is untouched: no distillation and no fine-tuning. For Edges2Handbags and DIODE the paper samples the exact checkpoints the DDBM authors released (284M and 534M parameters). ImageNet uses DBIM's checkpoint, which was fine-tuned from a pretrained image generator; Face2Comics and CelebAMask-HQ have no public bridge checkpoint, so the authors trained bridges from scratch with the ADM U-Net of Diffusion Models Beat GANs and then sampled them with each sampler in turn.

The explanation needs six pieces, in this order: 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, and the second-order approximation of the one integral it leaves; the singular endpoint that forces the first step to be random; and the algorithm with its budget of $2N-2$ network calls for $N$ steps.

How a diffusion bridge works

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

Two summaries of the schedule appear in every formula below. The signal-to-noise ratio $\text{SNR}_t = \alpha_t^2/\sigma_t^2$ measures 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 solver integrates in. The bridges in this paper use the variance-preserving (VP) 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 ($\beta_{\min}=0.1$, $\beta_d=2$) are smaller than the VP schedule of Song et al., so at $t=T=1$ the signal factor is still $\alpha_T = 0.577$ and the SNR is 0.499: the far end of this schedule is about half image and half noise. A bridge has no need to reach pure noise, because its far end is going to be pinned to an image. 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, a span of 6.10.

To pin the process, condition that forward diffusion on arriving at a chosen point $\mathbf{x}_T$ at time $T$. Conditioning a diffusion on its endpoint is a classical operation called Doob's h-transform, and its whole effect is one extra drift term, $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 taken of the log-probability of reaching the destination from the current point, so the term pushes the state in the direction that makes arrival at $\mathbf{x}_T$ more likely. The unconditioned process is Gaussian, so that transition probability is Gaussian and its gradient has a closed form, the paper's equation (4):

It behaves like a spring that pulls $\mathbf{x}_t$ toward a rescaled copy of the destination. Its denominator goes to zero as $t\to T$, so the pull becomes infinitely stiff at the end and arrival is certain. No network appears in (4): it is a closed-form function of the schedule. The same vanishing denominator makes the exact solver singular at its starting point, which is the subject of the section on the first step.

The pinned process is a diffusion bridge: it leaves $\mathbf{x}_0$ at time 0 and arrives at $\mathbf{x}_T$ at time $T$ with probability one. (This is DDBM's construction. A Schrödinger bridge, the object of the DSBM explainer, solves a different, entropy-regularized transport problem, and the paper notes its solutions apply to DDBM's framework only.) The state at an intermediate time is Gaussian with a closed-form mean and variance. Write $r_t = \text{SNR}_T/\text{SNR}_t$, which runs from 0 at the output end to 1 at the input end:

The mean blends the two endpoints, and the variance is the ordinary diffusion variance $\sigma_t^2$ shrunk by $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$, $r_T = 1$, so the blend is all $\mathbf{x}_T$ and the variance is zero again. At $t=0.5$ on this schedule the weights are $a_t=0.260$ on the input and $b_t=0.710$ on the output, with standard deviation $c_t = 0.463$ on the $[-1,1]$ pixel scale: a pixel whose output value is 0.8 and input value is −0.2 sits, on average, at $0.260\cdot(-0.2) + 0.710\cdot 0.8 = 0.516$, give or take 0.463.

Figure 1 draws one pixel's bridge. Watch the amber band pinch to zero width at both ends, and move the probe to read the three coefficients at any time. The variance-exploding (VE) schedule, which the paper also supports, has $\alpha_t = 1$ and $\sigma_t = t$ up to $T=80$; the band changes shape but stays pinned at both ends.

Every path ends at $\mathbf{x}_T$, so the reverse process can start from the input image itself instead of from noise. A distribution with zero width at its start also leaves a deterministic solver nothing to spread 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. A training step draws an (output, input) pair and a time, samples $\mathbf{x}_t$ directly from the Gaussian (3) (the process is never simulated), and regresses the prediction onto $\mathbf{x}_0$ with a squared error. For DIODE, $\mathbf{x}_t$ is a $3\times256\times256$ tensor of pixel values in $[-1,1]$. The U-Net receives six channels, the noisy image times a preconditioning factor followed by the raw input image, takes $250\ln t$ as its time input, and returns a $3\times256\times256$ prediction of the output. The paper writes this $\mathbf{D}_\theta(\mathbf{x}_s)$ for short and calls the model $\mathbf{x}_0$-predicting. A solver could use a noise prediction instead, but the clean image is bounded, so the weights a solver puts on it stay bounded even when one step covers a large range of noise levels; weights on a noise prediction grow without limit toward the clean end, which is also why DPM-Solver++ predicts the clean image for ordinary diffusion.

The reverse equations are 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 bridge version of Tweedie's formula connects the two. 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. A network trained with squared error converges to that posterior mean, so the paper's equation (6) is:

The identity is exact when $\mathbf{D}_\theta$ equals the posterior mean, and off by the network's error otherwise. (The paper introduces (3)'s score as the score of $p_t(\mathbf{x}_t\mid\mathbf{x}_T)$; since the expression contains $\mathbf{x}_0$, it is the two-endpoint score, and (6) is the step that turns it into the one-endpoint score.) The $c_t^2$ in the denominator vanishes at both ends of the bridge, so an error of 0.01 in $\mathbf{D}_\theta$ becomes an error of $0.01\,b_t/c_t^2$ in the score, which grows without bound near either end. DBMSolver therefore carries the prediction $\mathbf{D}_\theta$ through every step and converts to a score only once, in a final Euler step at $t=10^{-4}$.

The reverse SDE and the probability-flow ODE

A forward diffusion can be reversed in two ways, and both carry over to the bridge. One is a stochastic differential equation that injects fresh noise as it runs. The other is an ordinary differential equation with no noise; it produces the same distribution of states at every time, although a single run of it looks nothing like a single run of the SDE. With drift $f(\mathbf{x}_t,t)$ and diffusion coefficient $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 terms. The base drift is known: for this schedule $f = -(0.05 + t)\,\mathbf{x}_t$ and $g^2 = 0.1 + 2t$. The score of the bridge state is the learned term, supplied by the network through (6). The h-transform drift toward the input is the closed form (4). The SDE (the paper's Bridge SDE) adds the noise term $g(t)\,\mathrm{d}\mathbf{w}_t$; the ODE (the Bridge probability-flow ODE, or PF ODE) has none.

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 noise in the SDE spreads the distribution; a noiseless flow reproduces the same spreading only by converting that diffusion into a transport term, and that conversion costs half the score; every probability-flow ODE has its $\tfrac12$ for this reason (Song et al.). The h-transform term belongs to the forward drift that defines the bridge, so the conversion leaves it at full strength.

The fast samplers built for ordinary diffusion, DDIM and DPM-Solver++, are exact or high-order solutions of an ODE that has no h-transform term and a different linear coefficient. Applied to a bridge network they integrate that other ODE: nothing steers the trajectory to $\mathbf{x}_T$, and the network receives states unlike the ones it was trained on. The paper attributes the failure to their Gaussian-prior assumption; starting them from the input image does not help, because their update formulas encode the wrong forward process. On CelebAMask-HQ at 6 NFE, DPM-Solver++ scores an FID of 169 against DBMSolver's 35, and its images keep the layout of the mask but lose colour and texture.

Why an exponential integrator fits the bridge ODE

Substitute the closed forms (4) and (6) into the ODE (5) and collect terms in $\mathbf{x}_t$. Every term except the network's prediction becomes 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 plus a remainder. Here the remainder is a known multiple of the network output plus a known multiple of $\mathbf{x}_T$, so all of the nonlinearity sits inside $\mathbf{D}_\theta$, which depends on $\mathbf{x}_t$ however a U-Net does. For the VP schedule the linear coefficient is $L(t) = \tfrac{\mathrm{d}}{\mathrm{d}t}\log\sigma_t + \lambda_t'/\rho_t$, where $\rho_t = \text{SNR}_t/\text{SNR}_T - 1$ is zero at $t=T$, so $L$ blows up at the input end. DPM-Solver used the same structure for ordinary diffusion in 2022, and DBIM wrote it down for bridges in 2024.

Euler's method advances the state by the right-hand side times the step, treating the whole right-hand side as constant across the step. For the linear part that throws away an exact answer, and the equation is stiff: $L(t)$ inherits the $1/c_t^2$ factors of (4) and (6), so it is enormous near both ends of the bridge, and an Euler step of any useful size there overshoots and diverges. An exponential integrator treats $L(t)\mathbf{x}_t$ as the equation and $N$ as a forcing term. The variation-of-constants formula then gives the state at the end of a step from $s$ to $t$ exactly:

For this schedule the exponential factor is a closed-form number, so no step-size error is attached to the linear part. What remains to approximate is how the network's prediction varies inside the second integral. Freeze it at its value at the start of the step and you get a first-order method that is exact whenever the prediction is constant. Estimate its slope as well and you get a second-order method. Neither has an error term that depends on the size of $L$.

Figure 2 is the same situation in one dimension: $\mathrm{d}x/\mathrm{d}t = -k\,x + N(t)$, with a decay rate $k$ in place of $L(t)$ and a smooth forcing $N(t)$ in place of the network term. Raise $k$ until the readout shows $k\cdot h$ above 2 and watch the gray Euler line start to zigzag and grow, while both teal lines stay on the amber curve.

In the default setting ($k=3$, 6 steps) the maximum errors are 17.2 for Euler, 0.393 for the first-order integrator and 0.070 for the second-order one. 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 after the $1/h$ steps needed to cross a fixed interval the accumulated error shrinks like $h^k$. Halving the step halves the final error of a first-order method and quarters that of a second-order one. The paper's "error $O(\Delta t^2)$" for its first-order step is the per-step figure; the error in the final image is one power lower.

The integrals are taken in $\lambda$ rather than in $t$, a choice DPM-Solver introduced and DBIM and this paper keep. In $\lambda$ the exponential factor and the weight inside the integral become exponentials of $\lambda$, and the network's prediction is Taylor-expanded in $\lambda$ too. Since the network's time input is $250\ln t$, a uniform step in $\lambda$ corresponds more closely to a uniform step in that input than a uniform step in $t$ does.

How the DBMSolver step is derived

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

The first term carries the current state forward with a known weight. The second is a known weight on the input image; ordinary diffusion solvers have no such term. It appears because the score term enters the ODE at half strength and cancels only half of the h-transform's pull toward the input (in the SDE the two cancel exactly). The third term is an integral of the network's prediction over the step, and it is the only place where anything is approximated.

The simplest approximation freezes the prediction at its value at the start of the step. The integral of the weight alone 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). Coefficient for coefficient, this is the deterministic step of DBIM, which the paper acknowledges. The listing gives the three weights in code and their values for the first ODE step of a 6-NFE run, which carries the state from $t=0.9999$ to $t=0.0412$, a span of 2.92 in $\lambda$ (from −0.35 to 2.57):

# 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 that step the state is 0.993 times the network's guess of the output, plus 5.25 times the current state's displacement 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 5.25 is taken up in the section on the first step.

To see what freezing the prediction gets wrong, divide the third term of (9) by its first-order coefficient $c_D$. What is left is a weighted average of the network's prediction over the step, with weight $e^{2\lambda}/\sqrt{\rho_\lambda}$ at each $\lambda$, and the first-order step replaces that average by the value at the start. The substitution $u = \sqrt{\rho_\lambda}$ makes the weight easy to read. Since $\mathrm{d}\rho = 2(\rho+1)\,\mathrm{d}\lambda$ and $e^{2\lambda} = e^{2\lambda_T}(\rho+1)$, the weight becomes a constant:

So the weighted average over the step is a plain average over $u$, and once $\lambda$ is a unit or so past $\lambda_T$, $u$ grows like $e^{\lambda-\lambda_T}$. On the first ODE step of the 6-NFE run, $u$ goes from 0.018 to 18.5; the upper half of that range is $\lambda$ above 1.88, the last 0.69 of the 2.92-unit step. The weight's centre sits 2.00 units past the start, and 52% of the weight lies in the last quarter of the step. The first-order step evaluates the network at the noisy end, where its guess is still a blend of every plausible output, and uses that blend for a step whose weight is concentrated near the output end.

The second-order step expands the prediction in $\lambda$ around the start of the step, the paper's equation (10), and keeps the slope:

The derivatives $\mathbf{D}_\theta^{(n)}$ must be estimated from network calls; the weight integrals depend only on the schedule and are worked out once by hand. In the variable $u$ the $n=1$ weight is $e^{2\lambda_T}\int(\lambda(u)-\lambda_s)\,\mathrm{d}u$, and the antiderivative of $\log(1+u^2)$ is $u\log(1+u^2) - 2u + 2\arctan u$; the arctangent in (S26) comes from that antiderivative. The full step, supplement equation 26, is the first-order step plus one correction:

The network returns predictions, not slopes, so $\mathbf{D}_\theta^{(1)}$ is estimated the way DPM-Solver's second-order singlestep method does it. Take the first-order step only as far as a midpoint of the interval, 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 measures the current step rather than reusing a value from the previous one. For the first ODE step of the 6-NFE run the $\lambda$-midpoint sits at $t = 0.274$, the bracket in (S26) times its prefactor comes to 1.99, and the finite difference divides by a gap of 1.46. The paper does not say where the midpoint goes; the released code puts it at the $\lambda$-midpoint for DIODE and at the $t$-midpoint for Edges2Handbags. Both give valid slope estimates, because the difference is always divided by the actual $\lambda$ gap, but they are different samplers.

Figure 3 runs this on a stand-in for one pixel of the network's guess, which sits at 0 (the average of two possible outputs) near the input end and rises to 0.75 around some later $\lambda$. The default step is the first ODE step of the 6-NFE run. Look at the amber strip at the bottom, which is the step's weight: most of it is at the right end. Then compare the readouts. At the default the exact weighted average is 0.521, the first-order estimate is 0.015 and the second-order estimate is 0.452. Switching to a $t$-midpoint puts the probe at $t=0.521$ ($\lambda=0.48$), before the guess has moved much, and the estimate drops to 0.236.

Push the step end to 5.76, a single step across the whole bridge, and the second-order estimate overshoots to 1.22 against an exact 0.74: a straight line fitted near the start cannot follow a guess that bends and flattens. Shrink the step to end at $\lambda=0.5$ and the guess barely moves within it: the first-order estimate is off by 0.032 and the second-order estimate by 0.007.

The paper stops at $k=2$ because the next weight, with $(\lambda-\lambda_s)^2$, becomes an integral of $\log^2(1+u^2)$ in $u$, which has no elementary antiderivative and needs the dilogarithm. DBIM reaches second and third order a different way, by reusing predictions from earlier steps in a multistep formula (one call per step, extrapolating from the previous interval). The obstruction comes from the expansion variable. DBIM expands in $\log(b_t/c_t)$, and from (3), $b_t/c_t = e^{\lambda_T}\sqrt{\rho_t} = e^{\lambda_T}u$, so DBIM's variable is $\lambda_T + \log u$. Polynomials in $\log u$ integrate against $\mathrm{d}u$ in elementary terms at every order. DBIM's slope is a derivative in that other variable, so the two second-order methods are different schemes, not one scheme with and without memory.

Why the first step has to be random

Proposition 2 holds only for $s<T$, and the coefficients in (9) show why. At $s = T$, $\rho_s = 0$, 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, and the exact ODE solution cannot start from the input image. In Figure 1's picture, a noiseless flow started at the pinched right end stays a single point and cannot produce the spread of outputs the bridge has at every later time.

DBMSolver takes its first step, from $T$ to $T-\epsilon$ with $\epsilon = 10^{-4}$, with the Bridge SDE (2) instead. That equation has the same semi-linear structure: 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 exactly. Freezing the prediction and integrating the noise exactly gives Proposition 1, the first-order SDE step:

With the prediction held fixed, (8) is a draw from the Gaussian of the forward process at time $t$ given the current state and a clean image equal to the prediction. Started at $s=T$, its three coefficients are the bridge coefficients $a_t, b_t, c_t$ of (3): $0.9998$ on the input, $0.00018$ on the network's guess, and $0.0145$ on the noise. The state moves by a standard deviation of about 0.7% of the $[-1,1]$ pixel range.

The next step is the first ODE step, and it multiplies the state's displacement from the input by $c_x = 5.25$. A displacement with standard deviation 0.0145 becomes one with standard deviation 0.0761, which is $c_{0.041}$, the bridge's own standard deviation at the landing time (there it agrees with $\sigma_{0.041}$ to three decimals). $c_x$ grows like one over the square root of the gap $T-s$, the SDE step's noise shrinks like the square root of the same gap, and as the gap goes to zero their product tends to $c_t$ exactly, which follows from the formulas for $c_x$ and $c_s$. Figure 4 plots both, with the landing time held at $t=0.5$, as the gap goes from one half down to one millionth.

At the default gap the readouts give $c_x = 31.9$, $c_T = -31.7$, a sum of 0.267, a noise standard deviation of 0.0145 and a product of 0.463, which is $c_{0.5}$ from Figure 1. Moving the slider changes $c_x$ by orders of magnitude and leaves the sum and the product in place. The seed of the one draw in (8) therefore decides which output the deterministic trajectory reaches, and the supplement shows several faces generated from one segmentation map by changing only that seed.

Because only the first step is random, DBMSolver has no knob for how stochastic sampling is; DBIM's $\eta$, which mixes noise into every step, has no counterpart. The same first draw appears in DBIM (its "booting" noise) and in ODES3, the stochastic-start sampler the paper compares against. ODES3's own Theorem 3 shows the draw is closer to the true distribution than one Euler–Maruyama step of the SDE (Euler's method with a noise term added).

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:

With $h=\lambda_t-\lambda_s$ and $\text{SNR}_s/\text{SNR}_t = e^{-2h}$, the correct expression expands to $\alpha_t(h^2 - \tfrac23h^3 + \dots)$, as a second-order correction must, while the printed one starts at $-h$ and would cancel the method's order. The released code does not use the printed formula, and the paper's default first step is the first-order SDE step, so no reported number depends on it. The paper's ablation shows the second-order SDE step would not help anyway: FID 3.38 either way at 4 steps, 2.06 against 2.07 at 11.

The algorithm, call by call

Algorithm 1 of the paper, with the details the released code settles, goes like this: pick $N$ step times from $t_N = T$ down to $t_0 = 0$, with $t_{N-1} = T - 10^{-4}$ and the interior times on the Karras grid (EDM's power-law spacing with exponent 7, dense near $t=0$) down to $t_1 = 10^{-4}$; call the network once at the input 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, and the first ODE step then runs from 0.9999 to 0.552. Figure 5 lays the grid out in $t$ and in $\lambda$ and marks every call; slide $N$ and watch the first ODE step's share of the $\lambda$ axis.

Karras spacing puts most points near $t=0$ in $t$, yet in $\lambda$ the four-step grid is close to uniform, and the first ODE step covers 2.92 of the 6.10 units of half log-SNR. The midpoint call matters most on that step, because the network's guess changes most between the input end and the middle of the bridge, which is the situation of Figure 3.

The code settles three details 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 and leaves about half of $\sigma_{10^{-4}} = 0.0032$ as residual noise, below one 8-bit gray level ($2/255 \approx 0.0078$). Every prediction is clamped to $[-1,1]$ before use. And the network never receives $T$; its inputs are the scaled noisy image, the input image as three more channels, and $250\ln t$. In pseudocode:

# 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]

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)

Three solvers at the same budget

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. The ideal predictor $\mathbb{E}[\mathbf{x}_0\mid\mathbf{x}_t,\mathbf{x}_T]$ then has a closed form, so every sampler in Figure 6 runs on the same perfect network, the VP schedule and the Karras grid of the code, and only the solver differs. Each seed's exact endpoint, from integrating the ODE with 1,500 small steps, is marked by a ring, and the segment from the ring to a sampler's endpoint is that seed's error.

Euler on the raw ODE diverges at every budget shown. Its first step starts next to the pinned end, where $L(t)$ is largest, and multiplies the small deviation from the input by hundreds instead of by the exact factor, so the state leaves the canvas. The paper's Hybrid Heun at 5 NFE (FID 442) and DPM-Solver++ at 6 NFE (FID 169) on CelebAMask-HQ are the image-scale version of a plain solver failing near that end. At 2 NFE every method is the same two calls and diverges.

The first-order step is stable but, at 6 NFE, herds most seeds toward the mean of the three outputs. Figure 3 shows the mechanism: a frozen prediction is the network's guess from the noisy end of a wide step, which is a blend of all the outputs. In images that appears as the blur and missing texture the paper reports for DBIM at low NFE. The second-order step separates the modes on the same budget. Its median endpoint error on this toy is 1.14 against the first-order step's 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. From 10 to 20 to 40 NFE the second-order error falls by 2.7× and then 3.0× per doubling, the first-order error by 2.1× and 1.9×.

The toy's predictor is exact, which no real network is. A solver removes discretization error and leaves the network's own error in $\mathbf{D}_\theta$ untouched, and the paper's explanation of its weakest result, ImageNet inpainting, is that the network error dominates there.

DBMSolver results: FID against network calls

The paper reports FID, the Fréchet distance between Inception-network features of generated and reference images (lower is better; it compares distributions, not paired images), along with Inception score, LPIPS (a learned perceptual distance between paired images), mean squared error, and classifier accuracy for inpainting. Baselines are Hybrid Heun (119 NFE), DBIM at first, second and third order (20 NFE), ODES3 (28 NFE), and two distilled models that need a training run. Figure 7 plots every row; switch datasets and look at where the teal points sit relative to the others.

On DIODE (normals to outdoor scene, 256×256), six calls give FID 3.38 against Hybrid Heun's 4.43 at 119 and DBIM-3's 4.23 at 20. Twenty calls give 2.06. The abstract's "53% drop" is 2.06 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, which measure closeness to the specific target image, are essentially tied across the bridge samplers (LPIPS 0.196 to 0.203, MSE 0.015 to 0.018), so the gain is in how realistic the set of outputs looks rather than in per-pixel closeness.

On Edges2Handbags (64×64), six calls give 0.93 (Table 3; the text and Table 6 give 0.97 for the same run) and twenty give 0.53. These FIDs are computed against the full training set, following the DDBM protocol, which is why they are much smaller than the FIDs usually reported for ImageNet generation; compare them only with each other.

On Face2Comics, DBMSolver reaches 0.96 at 20 NFE where DBIM at the same budget scores between 8.6 and 9.3 and Hybrid Heun at 119 scores 2.36. At 6 NFE it scores 3.04, better than every DBIM row but worse than Hybrid Heun, so on this task DBMSolver matches Hybrid Heun only at 20 calls, a saving of about 6× rather than 20×. On CelebAMask-HQ it reaches 17.56 at 20 NFE against DBIM-3's 19.49, and 34.76 at 6 NFE against DBIM-1's 44.92. The GAN and other-diffusion rows in Table 4 (Pix2Pix, CycleGAN, SPADE, OASIS, LDM, BBDM and others) match the Brownian Bridge Diffusion Models (BBDM) paper's Tables 1 and 2 to two decimals, so they come from BBDM's own runs, not from the pipeline that produced the bridge rows; 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 exception, and the paper names it as its limitation. At 20 NFE DBMSolver, DBIM-2 and DBIM-3 all score FID 4.07 with classifier accuracy 72.0. At 6 NFE DBMSolver scores 4.98, better than DBIM-1 at 6 (5.36) but worse than DBIM-1 at 20 (4.13). Table 5's caption claims superiority "across all NFEs", which its own rows contradict. The authors' hypothesis is the one from the toy: on a varied dataset the network's prediction error dominates the discretization error, and a better solver can only reduce the second.

Against ODES3 the margin on Edges2Handbags is 0.53 at 20 NFE against 0.54 at 28, with ODES3 marginally better on Inception score (3.65 against 3.64) and LPIPS (0.097 against 0.099); on DIODE it is 2.06 against 2.29. Against distillation, the paper quotes 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 and 4.07), both worse than DBMSolver at 6. The same consistency-bridge paper also reports a stronger variant, consistency bridge training (CBT), at 0.80 and 2.93 with 2 NFE, which beats DBMSolver at 6 (the Consistency Models explainer covers the idea these build on). On ImageNet inpainting the paper's Table 6 lists the consistency-bridge model at 2 NFE (5.65) and IBMD at 1 NFE (5.87), both behind DBMSolver's 4.98 at 6. IBMD at 4 NFE reaches 4.03 on ImageNet inpainting, so with a training run a distilled model wins per call.

The abstract's "up to 5×" fewer calls has no matching pair of rows in the tables. Two speedups are supported: 6 calls against DBIM's 20 (3.3×) at better FID on three of the five tasks (DIODE, Edges2Handbags, Face2Comics), and 6 against Hybrid Heun's 119 (about 20×) at better FID on four (all but Face2Comics). For an implementer the recipe is the code listing above applied to any DDBM checkpoint with its schedule and its $\mathbf{x}_0$-predicting denoiser; on DIODE, four steps cost six calls and score FID 3.38.