AdjointDEIS: Efficient Gradients for Diffusion Models

Authors: zander w blasingame
chen  liu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:49:19

Abstract: The optimization of the latents and parameters of diffusion models with respect to some differentiable metric defined on the output of the model is a challenging and complex problem. The sampling for diffusion models is done by solving either the *probability flow* ODE or diffusion SDE wherein a neural network approximates the score function allowing a numerical ODE/SDE solver to be used. However, naive backpropagation techniques are memory intensive, requiring the storage of all intermediate states, and face additional complexity in handling the injected noise from the diffusion term of the diffusion SDE. We propose a novel family of bespoke ODE solvers to the continuous adjoint equations for diffusion models, which we call *AdjointDEIS*. We exploit the unique construction of diffusion SDEs to further simplify the formulation of the continuous adjoint equations using *exponential integrators*. Moreover, we provide convergence order guarantees for our bespoke solvers. Significantly, we show that continuous adjoint equations for diffusion SDEs actually simplify to a simple ODE. Lastly, we demonstrate the effectiveness of AdjointDEIS for guided generation with an adversarial attack in the form of the face morphing problem. Our code will be released on our project page [https://zblasingame.github.io/AdjointDEIS/](https://zblasingame.github.io/AdjointDEIS/)

Main paper: 
 introduction: 
 Diffusion models are a large family of state-of-the-art generative models which learn to map samples drawn from white Gaussian noise into the data distribution [1, 2] . These diffusion models have achieved state-of-the-art performance on prominent tasks such as image generation [3] [4] [5] , audio generation [6, 7] , or video generation [8] . Often, the state-of-the-art models are quite large and training them is prohibitively expensive [9] . As such, it is fairly common to adapt a pre-trained model to a specific task for post-training. This way, the generative model can learn new concepts, identities, or tasks without needing to train the whole model [10] [11] [12] . Additional work has also proposed algorithms for guiding the generative process of the diffusion models [13, 14] .One method to guide or direct the generative process is to solve an optimization problem w.r.t. some guidance function L defined on the image space R d . This guidance function works on the output of the diffusion model and assesses how "good" the output is. However, the diffusion model works by iteratively removing noise until a clean sample is reached. As such, we need to be able to efficiently backpropagate gradients through the entire generative process. As Song et al. [15] showed, the diffusion SDE can be simplified to an associated ODE, and as such, many efficient ODE/SDE solvers have been developed for diffusion models [16] [17] [18] . However, naïvely applying backpropagation to the diffusion model is inflexible and memory intensive; moreover, such an approach is not trivial to apply to the diffusion models that used an SDE solver instead of an ODE solver.AdjointDEISx t N ϵ θ z ϵ t N x t N -1 ϵ θ z ϵ t1x t0. . .x t0 L ϵ θ x t1x tN ϵ θ a x ( t0 ) a x ( t1 ) a x ( tM ) . . . . . .Figure 1 : A high-level overview of the AdjointDEIS solver to the continuous adjoint equations for diffusion models. The sampling schedule consists of {t n } N n=0 timesteps for the diffusion model and { tn } M n=0 timesteps for AdjointDEIS. The gradients a x (T ) can be used to optimize x T to find some optimal x * T .Inspired by the work of [19, 20] we study the application of the continuous adjoint equations to diffusion models, with a focus on training-free guided generation with diffusion models. We introduce several theoretical contributions and technical insights to both improve the ability to perform certain guided generation tasks and to gain insight into guided generation with diffusion models.First, we introduce AdjointDEIS a bespoke family of ODE solvers which can efficiently solve the continuous adjoint equations for both diffusion ODEs and SDEs. To the best of our knowledge, this is the first general backpropagation technique designed for diffusion SDEs. Moreover, we show that the continuous adjoint equations for diffusion SDEs simplify to a mere ODE.Next, we show how to calculate the continuous adjoint equation for conditional information which evolves with time (rather than being constant). To the best of our knowledge, we are the first to consider conditional information which evolves with time for neural ODEs.Overall, multiple theoretical contributions and technical insights are provided to bring a new family of techniques for the guided generation of diffusion models, which we evaluate experimentally on the task of face morphing. 
diffusion models: 
 In this section, we provide a brief overview of diffusion models. Diffusion models learn a generative process by first perturbing the data distribution into an isotropic Gaussian by progressively adding Gaussian noise to the data distribution. Then, a neural network is trained to perform denoising steps, allowing for sampling of the data distribution via sampling of a Gaussian distribution [2, 9] . Assume we have an n-dimensional random variable x ∈ R n with some distribution p data (x). Then diffusion models begin by diffusing p data (x) according to the diffusion SDE [2] , an Itô SDE given asdx t = f (t)x t dt + g(t) dw t (2.1)where t ∈ [0, T ] denotes time with fixed constant T > 0, f (•) and g(•) denote the drift and diffusion coefficients, and w t denotes the standard Wiener process. The trajectories of x t follow the distributions p t (x t ) with p 0 (x 0 ) ≡ p data (x) and p T (x T ) ≈ N (0, I). Under some regularity conditions Song et al. [15] showed that Equation (2.1) has a reverse process as time runs backwards from T to 0 with initial marginal distribution p T (x T ) governed bydx t = [f (t)x t -g 2 (t)∇ x log p t (x t )] dt + g(t) d wt (2.2)where wt is the standard Wiener process as time runs backwards. Solving Equation (2.2) is what allows diffusion models to draw samples from p data (x) by sampling p T (x T ). The unknown term in Equation (2.2) is the score function ∇ x log p t (x t ), which in practice is modeled by a neural network that estimates the scaled score function, ϵ θ (x t , t) ≈ -σ t ∇ x log p t (x t ), or some closely related quantity like x 0 -prediction [1, 2, 21] .The practical choice of a step size when discretizing SDEs is limited by the randomness of the Wiener process as a large step size, i.e., a small number of steps, can cause non-convergence, particularly in high-dimensional spaces [16] . Sampling an equivalent Ordinary Differential Equation (ODE) over an SDE would enable faster sampling. Song et al. [15] showed there exists an Ordinary Differential Equation (ODE) whose marginal distribution at time t is identical to that of Equation (2.2) given asdx t dt = f (t)x t - 1 2 g 2 (t)∇ x log p t (x t ).(2.3)The ODE in Equation (2.3) is known as the probability flow ODE [15] . As the noise prediction network, ϵ θ (x t , t), is trained to model the scaled score function, Equation (2.3) can be parameterized asEQUATIONw.r.t. the noise prediction network. For brevity, we refer to this as a diffusion ODE.While there exist several popular choices for the drift and diffusion coefficients, we opt to use the de facto choice which is known as the Variance Preserving (VP) type diffusion SDE [1, 15, 22] . The coefficients for VP-type SDEs are given asEQUATIONwhich corresponds to sampling x t from the distribution q(x t | x 0 ) = N (α t x 0 , σ 2 t I). 
adjoint diffusion odes: 
 Problem statement. Given the diffusion ODE in Equation (2.4), we wish to solve the following optimization problem:EQUATIONI.e., we wish to find the optimal x T , z, and θ which satisfy our guidance function L.Unlike GANs which can update the latent representation through GAN inversion [23, 24] , as seen in Equation (3.1) diffusion models require more care as they model an ODE or SDE and require numerical solvers. Therefore, to update the latent representation, model parameters, and conditional information we must backpropagate the gradient of loss defined on the output, ∂L(x 0 )/∂x 0 through the whole ODE or SDE.A key insight of this work is the connection between the adjoint ODE used in neural ODEs by Chen et al. [19] and specialized ODE/SDE solvers by [16] [17] [18] for diffusion models. It has been well observed that diffusion models are a type of neural ODE [15, 25] . Since a diffusion model can be thought of as a neural ODE, then we can solve the continuous adjoint equations [20] to find useful gradients for guided generation. We can then exploit the unique structure of diffusion models to develop efficient bespoke ODE solvers for the continuous adjoint equations.Let f θ describe a parameterized neural field of the probability flow ODE, i.e., the R.H.S of Equation (2.4), defined asEQUATIONThen f θ (x t , z, t) describes a neural ODE which admits an adjoint state, a x := ∂L/∂x t (and likewise for a z (t) and a θ (t)), which solve the continuous adjoint equations [20, Theorem 5.2] in the form of the following Initial Value Problem (IVP):EQUATIONWe refer to this system of equations in Equation (3.3) as the adjoint diffusion ODE1 as it describes the continuous adjoint equations for the empirical probability flow ODE.N.B., in the literature of diffusion models the sampling process is often done in reverse-time, i.e., the initial noise is x T and the final sample is x 0 . Due to this convention solving the adjoint diffusion ODE backwards actually means integrating forwards in time. Thus while diffusion models learn to compute x t from x s with s > t, the adjoint diffusion ODE seeks to compute a x (s) from a x (t).We show that rather than treating f θ as a black box, the specific structure of the probability flow ODE is carried over to the adjoint probability flow ODE, allowing the adjoint probability flow ODE to be simplified into a special exact formulation.By evaluating the gradient of f θ w.r.t. x t for each term in Equation (3.2) we can rewrite the adjoint diffusion ODE for a x (t) in Equation (3.3) asEQUATIONDue to the gradient of the drift term in Equation (3.4), further manipulations are required to put the empirical adjoint probability flow ODE into a sufficiently "nice" form. We follow the approach used by [16, 18] to simplify the empirical probability flow ODE with the use of exponential integrators and a change of variables. By applying the integrating factor exp t 0 f (τ ) dτ to Equation (3.4), the following ODE is expressed asEQUATIONThen, the exact solution at time s given time t < s is found to beEQUATIONLike with solvers for diffusion models which leverage exponential integrators, we are able to transform the adjoint diffusion ODE into a non-stiff form by separating the linear and non-linear component. Moreover, we can compute the linear in closed form, thereby eliminating the discretization error in the linear term. However, we still need to approximate the non-linear term which consists of a difficult integral about our complex neural network. This is where the insight of Lu et al. [16] to integrate in the log-SNR domain becomes invaluable. Let λ t := log(α t /σ t ) be one half of the log-SNR. Then, with using this new variable and computing the drift and diffusion coefficients in closed form, we can rewrite Equation (3.6) as EQUATIONEQUATIONEQUATION)a θ (s) = a θ (t) + λs λt α λ e -λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂θ dλ. (3.10)The full derivations of Proposition 3.1 can be found in Appendix A.1.There is a nice symmetry between Equations (3.8) to (3.10), while the adjoint of the solution trajectories evolves with a weighting of α t /α s in the linear term and the integral term is weighted by α 2 t /α 2 s reflecting the double partial ∂x t in the adjoint and Jacobian terms. Conversely, the adjoint state for the conditional information and model parameters evolves with no weighting on the linear term and the integral is only weighted by α t /α s . This follows from the vector fields being independent of a z and a θ . These equations while reflecting the special nature of this formulation of diffusion models also have an appealing parallel with the exact solution for diffusion ODEs Lu et al. [16, Proposition 3.1] .The numerical solver for the adjoint empirical probability flow ODE, now in light of Equation (3.8) , only needs to focus on approximating the exponentially weighted integral of ϵ θ from λ t to λ s , a well-studied problem in the literature on exponential integrators [26, 27] . To approximate this integral, we evaluate the Taylor expansion of the vector Jacobian product to further simplify the ODE. For k ≥ 1, the (k -1)-th Taylor expansion at λ t isEQUATIONPlugging this expansion into Equation (3.8) and letting h = λ s -λ t yieldEQUATIONWith this expansion, the number of terms which need to be estimated is further reduced as the exponentially weighted integral λs λt(λ-λt) n n!e -λ dλ can be solved analytically by applying n times integration-by-parts [16, 28] . Therefore, the only errors in solving this ODE occur in the approximation of the n-th order total derivatives of the vector-Jacobian product and the higher-order error terms O(h k+1 ). By dropping the O(h k+1 ) error term and approximating the first (k -1)-th derivatives of the vector-Jacobian product, we can derive k-th order solvers for adjoint diffusion ODEs. We decide to name such solvers as Adjoint Diffusion Exponential Integrator Sampler (AdjointDEIS) reflecting our use of the exponential integrator to simplify the ODEs and pay homage to DEIS from [18] that explored the use of exponential integrators for diffusion ODEs. Consider the case of k = 1, by dropping the error term O(h 2 ) we construct the AdjointDEIS-1 solver with the following algorithm.AdjointDEIS-1. Given an initial augmented adjoint state [a x (t), a z (t), a θ (t)] at time t ∈ (0, T ), the solution [a x (s), a z (s), a θ (s)] at time s ∈ (t, T ] is approximated byEQUATIONhigher-order expansions of Equation (3.11) require estimations of the n-th order derivatives of the vector Jacobian product which can be approximated via multi-step methods, such as Adams-Bashforth methods [29] . This has the added benefit of reduced computational overhead, as the multi-step method just reuses previous values to approximate the higher-order derivatives. Moreover, multi-steps are empirically more efficient than single-step methods [29] . Combining the Taylor expansions in Equation (3.11) with techniques for designing multi-step solvers, we propose a novel multi-step second-order solver for the adjoint empirical probability flow ODE which we call AdjointDEIS-2M. This algorithm combines the previous values of the vector Jacobian product at time t and time r to predict a s without any additional intermediate values.AdjointDEIS-2M. We assume having a previous solution a x (r) and model output ϵ θ (x r , z, r) at time r < t < s, let ρ denote ρ = λt-λr h . Then the solution a s at time s to Equation (3.4) is estimated to beEQUATIONEQUATIONFor brevity, we omit the details of the AdjointDEIS-2M solver for a z (t) and a θ (t); rather, we provide the full derivation and details in Appendix A. Likewise, the full algorithm can be found in Appendix F.1. The advantage of a higher-order solver is that it is generally more efficient, requiring fewer steps due to its higher convergence order. We show that AdjointDEIS-k is a k-th order solver, as stated in the following theorem. The proof is in Appendix B. Theorem 3.1 (AdjointDEIS-k as a k-th order solver). Assume the function ϵ θ (x t , z, t) and its associated vector-Jacobian products follow the regularity conditions detailed in Appendix B, then for k = 1, 2, AdjointDEIS-k is a k-th order solver for adjoint diffusion ODEs, i.e., for the sequence {ã x (t i )} M i=1 computed by AdjointDEIS-k, the global truncation error at time T satisfies ãx (t M )a x (T ) = O(h 2 max ), where h max = max 1≤j≤M (λ ti -λ ti-1 ). Likewise, AdjointDEIS-k is a k-th order solver for the estimated gradients w.r.t. z and θ.As previous work has shown that higher-order solvers may be unsuitable for large guidance scales [16] [17] [18] we do explicitly construct or analyze any solvers for k > 2 and leave such explorations for future study.Thus far, we have held the conditional information constant across time, i.e., the same text conditioning like a prompt "fire dragon" can be fed to the noise prediction network. In guided generation tasks it is not uncommon to take advantage of the iterative nature of diffusion models by scheduling the conditional information to exhibit different values at different timesteps. E.g., blending concepts by alternating between the prompt for "fire dragon" and the prompt for "ice dragon" at each timestep to create a picture of a dragon with both the qualities of ice and fire-a technique which has been popularized by Stable Diffusion [3] .We show that converting a constant z into a scheduled z t does not actually change the continuous adjoint equation for z t , meaning we can still apply the AdjointDEIS-k solvers derived above to find a z (t) by simply replacing z with z t . We state this observation more formally in the following theorem. The proof is in Appendix C. Theorem 3.2 (Continuous adjoint equations for time-dependent conditional information). Suppose there exists a function z : R → R z which is continuously differentiable in t and is the conditional input into a parameterized vector field f θ : R d × R z × R → R d be continuous in t, uniformly Lipschitz in x, and continuously differentiable in x. Let x : R → R d be the unique solution for the ODEdx t dt = f θ (x t , z t , t)with initial condition x 0 . Then there exists a unique solution a z : R → R z to the following IVP:a z (0) = 0, da z dt (t) = -a x (t) ⊤ ∂f θ (x t , z t , t) ∂z t .While motivated by the case of scheduled conditional information in guided generation with diffusion models, this result applies to neural ODEs more generally, which could open future research directions. 
adjoint diffusion sdes: 
 As recent work [30, 31] i.e., any gap in the mismatched prior distributions p t (x t ) and pt (x t ) for the true distribution p t and edited distribution pt will remain between p 0 (x 0 ) and p0 (x 0 ), whereas for diffusion SDEs the gap can be reduced between pt (x t ) and p t (x t ) as t tends towards 0. Motivated by this reasoning, we present a framework for solving the adjoint diffusion SDE using exponential integrators.The diffusion SDE with noise prediction model is given byEQUATIONwhere 'dt' is an infinitesimal negative timestep. Note how the drift term of the SDE looks remarkably similar to the probability flow ODE sans a missing factor of 1/2 in front of the noise prediction model. This is due to differing manipulations of the forward Kolomogorov equations-which describe the evolution of p t (x t )-used by Anderson [32] to derive the reverse-time SDE and later by Song et al. [15] to derive the probability-flow ODE. This connection is very important as it enables one to simplify the AdjointDEIS solvers for the adjoint diffusion SDE.We show that for the special case of Stratonovich SDEs2 with a diffusion coefficient g(t) which does not depend on the process state x t , then the adjoint process has a unique strong solution that evolves with what is essentially an ODE. Intuitively, this tracks as the stochastic term g(t) • dw t has nothing to do with x t . We state this observation somewhat informally in the following theorem. The proof can be found in Appendix D.Theorem 4.1. Let f : R d × R → R d be in C ∞,1 b and g : R → R d×w be in C 1 b . Let L : R d → R be a scalar-valued differentiable function. Let w t : [0, T ] → R w be a w-dimensional Wiener process. Let x : [0, T ] → R d solve the Stratonovich SDE dx t = f (x t , t) dt + g(t) • dw t ,with initial condition x 0 . Then the adjoint process a x (t) := ∂L(x T )/∂x t is a strong solution to the backwards-in-time ODEda x (t) = -a x (t) ⊤ ∂f ∂x t (x t , t) dt. (4.2)This is a boon for us as diffusion models model use only a mere scalar diffusion coefficient, g(t). Therefore, the continuous adjoint equations for the diffusion SDE just simplify to an ODE. Not only that, but as mentioned before, the drift term of the diffusion SDE and probability flow ODE differ only by a factor of 2 in the term with the noise prediction network. As only the drift term of the diffusion SDE is used when constructing the continuous adjoint equations, it follows that the only difference between the continuous adjoint equations for the probability flow ODE and diffusion SDE is a factor of 2. Therefore, the exact solutions are given by: Proposition 4.1 (Exact solution of adjoint diffusion SDEs). Given initial values [a x (t), a z (t), a θ (t)] at time t ∈ (0, T ), the solution [a x (s), a z (s), a θ (s)] at time s ∈ (t, T ] of adjoint diffusion SDEs isEQUATIONEQUATION)EQUATIONRemark 4.1. While the adjoint diffusion SDEs evolve with an ODE the same cannot be said for the underlying state, x t . Rather this evolves with a backwards SDE (more details in Appendix D) which requires the same realization of the Wiener process used to sample the image as the one used in the backwards SDE.Lu et al. [17] propose the following first-order solver for diffusion SDEsEQUATIONwhere ϵ s ∼ N (0, I). To solve the SDE backwards in time we follow the approach initially proposed by Wu and la Torre [33] and used by later works [31] . Given a particular realization of the Wiener process that admits x t ∼ N (α t x 0 | σ 2 t I), then for two samples x t and x s the noise ϵ s can be calculated by rearranging Equation (4.6) to findϵ s = x t -αt αs x s + 2σ t (e h -1)ϵ θ (x s , z, s) σ t √ e 2h -1 (4.7)With this the sequence {ϵ ti } N i=1 of added noises can be calculated which will exactly reconstruct the original input from the initial realization of the Wiener process. This technique is referred to as Cycle-SDE after the CycleDiffusion paper [33] . 
experiments: 
 To illustrate the efficacy of our technique, we examine the application of guided generation for the face morphing attack. The face morphing attack is a new emerging attack on Face Recognition (FR) systems. This attack works by creating a singular morphed face image x (ab) 0 that shares biometric information with the two contributing faces x (a) 0 and x (b) 0 [34] [35] [36] . A successfully created morphed face image can trigger a false accept with either of the two contributing identities in the targeted Face Recognition (FR) system, see Figure 2 for an illustration. Recent work in this space has explored the use of diffusion models to generate these powerful attacks [34, 37, 38] . All prior work on diffusionbased face morphing used a pre-trained diffusion autoencoder [39] trained on the FFHQ [40] dataset at a 256×256 resolution. We illustrate the use of the AdjointDEIS solvers by modifying the Diffusion Morph (DiM) architecture proposed by Blasingame and Liu [34] to use the AdjointDEIS solvers to find the optimal initial noise x (ab) T and conditional z ab . The AdjointDEIS solvers are used to calculate the gradients with respect to the identity loss [38] defined as ), and F : X → V is an FR system which embeds images into a vector space V which is equipped with a measure of distance, d. We used the ArcFace [41] FR system for the identity loss.EQUATIONWe compare against three preexisting DiM methods, the original DiM algorithm [34] , Fast-DiM [37] , and Morph-PIPE [38] as well as a GAN-inversion based face morphing attack, MIPGAN-I and MIPGAN-II [42] based on the StyleGAN [40] and StyleGAN2 [43] architectures respectively. Fast-DiM improves DiM by using high-order ODE solvers to decrease the number of sampling steps required to create a morph. Morph-PIPE performs a very simple version of guided generation by generating a large batch of morphed images derived from a discrete set of interpolations between x T , and z a and z b . For reference purposes, we compare against a reference GAN-based method [42] which uses GAN-inversion w.r.t.to the identity loss to find the optimal morphed face, and we include prior state-of-the-art Webmorph, a commercial off-the-shelf system [44] .We run our experiments on the SYN-MAD 2022 [44] morphed pairs which are constructed from the Face Research Lab London dataset [45] , more details in Appendix G.4. The morphed images are evaluated against three FR systems, the ArcFace [41] , ElasticFace [46] , and AdaFace [47] models, further details are found in Appendix G.5. To measure the efficacy of a morphing attack, the Mated Morph Presentation Match Rate (MMPMR) metric [48] is used. The MMPMR metric as proposed by Scherhag et al. [48] is defined asM (δ) = 1 M M m=1 min n∈{1,...,Nm} S n m > δ (5.2)where δ is the verification threshold, S n m is the similarity score of the n-th subject of morph m, N m is the total number of contributing subjects to morph m, and M is the total number of morphed images.In our experiments, we used a learning rate of 0.01, N = 20 sampling steps, M = 20 steps for AdjointDEIS, and 50 optimization steps for gradient descent. In Table 1 we present the effectiveness of the morphing attacks against the three FR systems. Guided generation with AdjointDEIS massively increases the performance of DiM, supplanting the old state-of-the-art for face morphing. Interestingly, the SDE variant did not fare as well as the ODE variant. This is likely due to the difficulty in discretizing SDEs with large step sizes [15] [16] [17] . We present further results in Appendix E which explore the impact of the choice of learning rate and number of discretization steps for AdjointDEIS. 
conclusion: 
 We present a unified view on guided generation by updating latent, conditional, and model information of diffusion models with a guidance function using the continuous adjoint equations. We propose AdjointDEIS, a family of solvers for the continuous adjoint equations of diffusion models. We exploit the unique construction of diffusion models to create efficient numerical solvers by using exponential integrators. We prove the convergence order of solvers and show that the continuous adjoint equations for diffusion SDEs evolve with an ODE. Furthermore, we show how to handle conditional information that is scheduled in time, further expanding the generalizability of the proposed technique. Our results in face morphing show that the gradients produced by AdjointDEIS can be used for guided generation tasks.Limitations. There are several limitations. Empirically, we only explored a small subset of the true potential AdjointDEIS by evaluating on a single scenario, i.e., face morphing. Likewise, we only explored a few different hyperparameter options. In particular, we did not explore much the impact of the number of optimization steps and the number of sampling steps for diffusion SDEs on the visual quality of the generated face morphs.Broader Impact. Guided generation techniques can be misused for a variety of harmful purposes. In particular, our approach provides a powerful tool for adversarial attacks. However, better knowledge of such techniques should hopefully help direct research in hardening systems against such kinds of attacks. 
 
 Appendix: In this section we provide the full derivations for the family of AdjointDEIS solvers. First recall the full definition of the continuous adjoint equations for the empirical probability flow ODE:a x (0) = ∂L ∂x 0 , da x dt (t) = -a x (t) ⊤ ∂f θ (x t , z, t) ∂x t , a z (0) = 0, da z dt (t) = -a x (t) ⊤ ∂f θ (x t , z, t) ∂z , a θ (0) = 0, da θ dt (t) = -a x (t) ⊤ ∂f θ (x t , z, t) ∂θ . (A.1)We can simplify the equations by explicitly solving gradients of the neural vector field f θ for the drift term to obtainda x dt (t) =-f (t)a x (t) - g 2 (t) 2σ t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x t , da z dt (t) = 0 - g 2 (t) 2σ t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂z , da θ dt (t) = 0 - g 2 (t) 2σ t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂θ . (A.2)Remark A.1. The last two equations in Equations (A.2) the vector fields are independent of (a z , a θ ), reducing these equations to mere integrals; however, it is often useful to compute the whole system a aug = (a x , a z , a θ ) as an augmented ODE. Remark A.2. Likewise, the last two equations in Equations (A.2) are functionally identical with a simple swap of z for θ or vice versa.As such, for the sake of brevity, the derivations for the AdjointDEIS solvers for (a z , a θ ) will only explicitly include the derivations for a z .Focusing first on the continuous adjoint equation for a x we apply the integrating factor expt 0 f (τ ) dτ to Equation (A.2) to find d dt e t 0 f (τ ) dτ a x (t) = -e t 0 f (τ ) dτ g 2 (t) 2σ t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x t . (A.3)Then, the exact solution at time s given time t < s is found to bee s 0 f (τ ) dτ a x (s) = e t 0 f (τ ) dτ a x (t) - s t e u 0 f (τ ) dτ g 2 (u) 2σ u a x (u) ⊤ ϵ θ (x u , z, u) ∂x u du a x (s) = e t s f (τ ) dτ a x (t) - s t e u s f (τ ) dτ g 2 (u) 2σ u a x (u) ⊤ ϵ θ (x u , z, u) ∂x u du (A.4)To simplify Equation (A.4), recall that f (t) is defined asf (t) = d log α t dt , (A.5)for VP type SDEs. Furthermore, let λ t := log(α t /σ t ) be one half of the log-SNR. Then the diffusion coefficient can be simplified using the log-derivative trick such thatg 2 (t) = dσ 2 t dt -2 d log α t dt σ 2 t = 2σ 2 t d log σ t dt - d log α t dt = -2σ 2 t dλ t dt . (A.6)Using this updated expression of g 2 (t) along with computing the integrating factor in closed form enables us to express Equation (A.4) asa x (s) = α t α s a x (t) + 1 α s s t α u σ u dλ u du a x (u) ⊤ ϵ θ (x u , z, u) ∂x u du. (A.7)Lastly, by rewriting the integral in terms of an exponentially weighted integral α u σ u = α 2 u σ u /α u = α 2 u e -a x (s) = α t α s a x (t) + 1 α s λs λt α 2 λ e -λ a x (λ) ⊤ ϵ θ (x λ , z, λ) ∂x λ dλ, (A.11) a z (s) = a z (t) + λs λt α λ e -λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂z dλ, (A.12) a θ (s) = a θ (t) + λs λt α λ e -λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂θ dλ. (A.13)Then to find the AdjointDEIS solvers we take a k-th order Taylor expansion about λ t and integrate in the log-SNR domain.For k ≥ 1, the (k -1)-th Taylor expansion at λ t of the inner term of the exponentially weighted integral in Equation (A.11) isa x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂x λ = k-1 n=0 (λ -λ t ) n n! d n dλ n α 2 λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂x λ λ=λt + O((λ -λ t ) k ). (A.14)Then plugging this into Equation (A.11) yieldsa x (s) = α t α s a x (t) + 1 α s λs λt e -λ k-1 n=0 (λ -λ t ) n n! d n dλ n α 2 λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂x λ λ=λt dλ +O(h k+1 ) = α t α s a x (t) + 1 α s k-1 n=0 d n dλ n α 2 λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂x λ λ=λt estimated λs λt (λ -λ t ) n n! e -λ dλ analytically computed + O(h k+1 ) omitted , (A.15)where h = λ s -λ t .The exponentially weighted integral λs λt(λ-λt) n n!e -λ dλ can be solved analytically by applying n times integration by parts [16, 28] such thatλs λt e -λ (λ -λ t ) n n! dλ = σ s α s h n+1 φ n+1 (h), (A.16)with the special φ-functions [26] . These functions are defined asφ n+1 (h) := 1 0 e (1-u)h u n n! du, φ 0 (h) = e h , (A.17)which satisfy the recurrence relation φ k+1 (h) = (φ n (h) -φ n (0))/h and have closed forms for k = 1, 2:φ 1 (h) = e h -1 h , (A.18) φ 2 (h) = e h -h -1 h 2 . (A.19)Likewise, the Taylor expansion of the exponentially weighted integral in Equation (A.12) yieldsa z (s) = a z (t) + λs λt e -λ k-1 n=0 (λ -λ t ) n n! d n dλ n α λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂z λ=λt dλ + O(h k+1 ) = a z (t) + k-1 n=0 d n dλ n α λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂z λ=λt estimated λs λt (λ -λ t ) n n! e -λ dλ analytically computed + O(h k+1 ) omitted .(A.20)For k = 1 and omitting the higher-order error term, Equation (A.15) becomes:a x (s) = α t α s a x (t) + 1 α s α 2 t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x t λs λt (λ -λ t ) 0 0! e -λ dλ = α t α s a x (t) + σ s (e h -1) α 2 t α 2 s a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x tBy Equation (A.16).(A.21)Likewise, the continuous adjoint equation for z, Equation (A.20), becomes when k = 1 by omitting the higher-order error term:a z (s) = a z (t) + α t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂z λs λt (λ -λ t ) 0 0! e -λ dλ = a z (t) + σ s (e h -1) α t α s a x (t) ⊤ ∂ϵ θ (x t , z, tBy Equation (A.16).And the first-order solver for a θ (t) can be found in a similar fashion, thus we have derived the AdjointDEIS-1 solvers.For notational convenience let V(x; t) denote the scaled vector-Jacobian product of the adjoint state a x (t) and the gradient of the model w.r.t. x t , i.e.,V(x; t) = α 2 t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x t . (A.23)Consider the following definition of the limit in the log-SNR domaind dλ α 2 λ a x (λ) ⊤ ∂ϵ θ (x λ , z, λ) ∂x λ = lim λr→λt V(x; λ t ) -V(x; λ r ) ρh , (A.24)where ρ = λt-λr h with h = λ s -λ t and where r is some previous step r < t < s. Again V(x; λ t ) is overloaded to mean V(x; t λ (λ t )). Then by omitting higher-order error O(h k+1 ), Equation (A.15) becomes:a x (s) = α t α s a x (t) + 1 α s V(x; λ t ) λs λt (λ -λ t ) 0 0! dλ + V (1) (x; λ t ) λs λt (λ -λ t ) 1 1! dλ = α t α s a x (t) + 1 α s σ s α s (e h -1)V(x; λ t ) + σ s α s (e h -h -1)V (1) (x; λ t ) . (A.25)By applying the same approximation used in Lu et al. [16] ofe h -h -1 h ≈ e h -1 2 , (A.26)then we can rewrite the second term of the Taylor expansion as along with a corresponding definition of first-derivative w.r.t. λ as defined in Equation (A.24). As such Equation (A.20), when k = 2, becomes the following when omitting the higher-order error term:σ s α s (e h -h -1)V (1) (x; λ t ) ≈ σ s α s (e h -h -1) V(x; λ t ) -V(x; λ r ) ρh By Equation (A.24) ≈ σ s α s e h -1 2ρ V(x; λ t ) -V(x; λ r ) By Equation (A.26) = σ s α s e h -1 2ρ α 2 t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂x t -α 2 r a x (r) ⊤ ∂ϵ θ (x r , z,a z (s) = a z (t) + V(z; λ t ) λs λt (λ -λ t ) 0 0! dλ + V (1) (z; λ t ) λs λt (λ -λ t ) 1 1! dλ = a z (t) + σ s α s (e h -1)V(z; λ t ) + σ s α s (e h -h -1)V (1) (z; λ t ). (A.30)The second term of the Taylor expansion can be rewritten asσ s α s (e h -h -1)V (1) (z; λ t ) ≈ σ s α s (e h -h -1) V(z; λ t ) -V(z; λ r ) ρh ≈ σ s α s e h -1 2ρ V(z; λ t ) -V(z; λ r ) By Equation (A.26) = σ s α s e h -1 2ρ α t a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂z -α r a x (r) ⊤ ∂ϵ θ (x r , z, r) ∂z .(A.31)Then Equation (A.30) becomesa z (s) = a z (t) + σ s (e h -1) α t α s a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂z + σ s e h -1 2ρ α t α s a x (t) ⊤ ∂ϵ θ (x t , z, t) ∂z - α r α s a x (r) ⊤ ∂ϵ θ (x r , z, r) ∂z , (A.32)and the corresponding second-order solver for a θ (t) can be found in a similar manner.For notational brevity we denote the scaled vector-Jacobian products of the solution trajectory of AdjointDEIS asṼ(x; t) = α 2 t ãx (t) ⊤ ∂ϵ θ (x t , z, t) ∂ xt , (B.1) Ṽ(z; t) = α t ãx (t) ⊤ ∂ϵ θ (x t , z, t) ∂z . (B.2)For the AdjointDEIS solvers, we make similar assumptions to Lu et al. [16] .Assumption B.1. The total derivatives of the vector-Jacobian products V (k) ({x λ , z, θ}, λ) as a function of λ exist and are continuous for 0 ≤ j ≤ k + 1 (and hence bounded with initial condition x 0 . Then the following map(a, t) → -a ⊤ ∂f θ (x t , z, t) ∂[x t , z, θ]is Lipschitz in a. Moreover, the Lipschitz constant L > 0 is given byL = sup t∈[0,T ] ∂f θ (x t , z, t) ∂x t . (B.3)Proof. Now as x t is continuous andf θ is continuously differentiable in x, so t → ∂f θ ∂[xt,z,θ] (x t , z, t) is a continuous function on the compact set [0, T ], so it is bounded by some L > 0. Likewise, for a ∈ R d the map (a, t) → -a ⊤ ∂f θ (xt,z,t) ∂[xt,z,θ]is Lipschitz in a with Lipschitz constant L and this constant is independent of t.Proof. First we consider the case of the adjoint state a x (t). Recall that the AdjointDEIS-1 solver for a x with higher-order error terms is given bya x (t i+1 ) = α ti α ti+1 a x (t i ) + σ ti+1 (e hi -1) α 2 ti α 2 ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂x ti + O(h 2 i ), (B.4)where we let t i = t, t i+1 = s, h i = λ ti+1 -λ ti from Equation (A.21). By Theorem B.1 and Equation (A.21) it holds thatãx (t i+1 ) = α ti α ti+1 ãx (t i ) + σ ti+1 (e hi -1) α 2 ti α 2 ti+1 ãx (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂ xti = α ti α ti+1 ãx (t i ) + σ ti+1 (e hi -1) α 2 ti α 2 ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂x ti + O(ã x (t i ) -a x (t i )) = α ti α ti+1 a x (t i ) + σ ti+1 (e hi -1) α 2 ti α 2 ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂x ti + O(ã x (t i ) -a x (t i )) = a x (t i+1 ) + O(h 2 max ) + O(ã x (t i ) -a x (t i )). (B.5)Repeat, this argument, from ãx (t 0 ) = a x (0) then we findãx (t M ) = a x (T ) + O(M h 2 max ) = a x (T ) + O(h max ). (B.6)Although the argument for the adjoint state a z (t) follows an analogous form to the one above we explicitly state it for completeness. Recall that the AdjointDEIS-1 solver for a z with higher-order error terms is given bya z (t i+1 ) = a z (t i ) + σ ti+1 (e hi -1) α ti α ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂z + O(h 2 i ). (B.7)By Theorem B.1 and Equation (A.22) it holds that ãz (t i+1 ) = ãz (t i ) + σ ti+1 (e hi -1)α ti α ti+1 ãx (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂z = ãz (t i ) + σ ti+1 (e hi -1) α ti α ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂z + O(ã x (t i ) -a x (t i )) = a z (t i ) + σ ti+1 (e hi -1) α ti α ti+1 a x (t i ) ⊤ ∂ϵ θ (x ti , z, t) ∂z + O(ã x (t i ) -a x (t i ) = a z (t i+1 ) + O(h 2 max ) + O(ã x (t i ) -a x (t i ). (B.8)Repeat, this argument, from ãz (t 0 ) = 0 then we findãz (t M ) = a z (T ) + O(M h 2 max ) = a z (T ) + O(h max ). (B.9)An identical argument can be constructed for a θ thereby finishing the proof.We prove the discretization error of the AdjointDEIS-2M solver. Note for the AdjointDEIS-2M solver we haveh i = λ ti+1 -λ ti-1 and ρ i = λt i -λt i-1 hi . Furthermore, let ∆ i = ∥ã x (t i ) -a x (t i )∥.Without loss of generality, we will prove this only for a x ; the derivation for a z and a θ is analogous.Proof. First we consider the case of the adjoint state a x (t). Recall that the AdjointDEIS-2, see Equation (A.25), solver for a x with higher-order error terms is given bya x (t i+1 ) = α ti α ti+1 a x (t i )+ 1 α ti+1 σ ti+1 α ti+1 (e hi -1)V(x; t i )+ σ ti+1 α ti+1 (e hi -h i -1)V (1) (x; t i ) +O(h 3 i ). (B.10)Taylor's expansion yieldsa x (t i+1 )- α ti α ti+1 a x (t i )+ 1 α ti+1 σ ti+1 α ti+1 (e hi -1)V(x; t i )+ σ ti+1 α ti+1 (e hi -h i -1)V (1) (x; t i ) ≤ Ch 3 i , (B.11)where C is a constant that depends on V (2) (x t , t). Also note thatV (1) (x; t i ) - 1 ρ i h i V(x; t i ) -V(x; t i-1 ) ≤ Ch i . (B.12)Since ρ i is bounded away from zero, and e -hi = 1 -EQUATION)where M i = 1 + sup ti≤t≤ti+1 ∥V (1) (x; t)∥ and L are the Lipschitz constants of V(x; t) by Theorem B.1. Then, ∆ i+1 can be estimated as∆ i+1 ≤ α ti α ti+1 ∆ i + σ ti+1 α 2 ti+1 L∆ i + σ ti+1 α 2 ti+1 (CM i h 3 i + CLh i (∆ i + ∆ i+1 )) + Ch 3 i ≤ α ti α ti+1 ∆ i + Ch i (∆ i + ∆ i+1 + h 2 i ). (B.14)Thus, ∆ i+1 = O(h 2 max ) as long as h max is sufficiently small and ∆ 0 + ∆ 1 = O(h 2 max ), which can be verified via Taylor expansion, thereby finishing the proof.For additional clarity, we let x(t) ≡ x t and likewise, z(t) ≡ z t .Proof. As z(t) is continuously differentiable w.r.t. t there exists some function z ′ (t) such thatdz dt (t) = z ′ (t). (C.1)Consider the augmented state defined byd dt x z (t) = f aug = f θ (x t , z t , t) z ′ (t) , (C.2)and the associated augmented adjoint statea aug (t) := a x a z (t). (C.3)The Jacobian of f aug has the form∂f aug ∂[x, z] = ∂f θ (x,z,t) ∂x ∂f θ (x,z,t) ∂z 0 0 . (C.4)Recall that a x (t) evolves withda x dt (t) = -a x (t) ⊤ ∂f θ (x(t), z(t), t) ∂x(t) . (C.5)Using Equations (C.3) and (C.4) we can define the evolution of the adjoint augmented state asda aug dt (t) = -[a x a z ] (t) ∂f aug ∂[x, z] (t) = -a x ∂f θ (x,z,t) ∂x a x ∂f θ (x,z,t) ∂z (t). (C.6)Therefore, a z (t) evolves with the ODEa z (0) = 0, da z dt (t) = -a x (t) ⊤ ∂f θ (x(t), z(t), t) ∂z(t) . (C.7)We have thus shown the evolution of a z (t) for some continuously differentiable function z(t).Now we prove the solution is unique and exists. As x(t) is continuous and f θ is continuously differentiable in x, it follows that the map t → ∂f θ ∂x (x(t), z(t), t) is a continuous function on the compact set [0, T ], and therefore it is bounded by some L > 0. Correspondingly, for a ∈ R d it follows that the map (a, t) → -a ⊤ ∂f θ ∂[x,z] (x(t), z(t), t) is Lipschitz in a with Lipschitz constant L and this constant is independent of t. Therefore, by the Picard-Lindelölf theorem [49, Theorem 110C] the solution a z (t) exists and is unique.In this section, we provide further details on the continuous adjoint equations for diffusion SDEs that we omitted from the main paper due to their technical nature and for the purpose of brevity.Consider the Itô integral given byx T = T 0 x t dw t , (D.1)where x t is a continuous semi-martingale adapted to the filtration generated by the Wiener process{w t } t∈[0,T ] , {F t } t∈[0,T ] .The following quantity, however, is not defined0 T x t dw t . (D.2)This is because x t and w t are adapted to {F t } t∈[0,T ] which is defined in forwards time. This means x t does not anticipate future events only depends on past events. While this is generally sufficient when we wish to integrate backwards in time we want future events to inform past events.Clearly, we need a different tool to model this backwards SDE. As such, taking inspiration from the work on neural SDEs [50] , we follow the treatment of Kunita [51] for the forward and backward Fisk-Stratonovich integrals using two-sided filtration. Let {F s,t } s≤t;s,t∈[0,T ] be a two-sided filtration, where F s,t is the σ-algebra generated by {w v -w u : s ≤ u ≤ v ≤ t} for s, t ∈ [0, T ] such that s ≤ t.Forward time. For a continuous semi-martingale {x t } t∈[0,T ] adapted to the forward filtration {F 0,t } t∈[0,t] , the Stratonovich stochastic integral is given asT 0 x t • dw t = lim |Π|→0 N k=1 x t k + x t k-1 2 (w t k -w t k-1 ) (D.3)whereΠ = {0 = t 0 < • • • < t N = T } is a partition of the interval [0, T ] and |Π| = max k t k -t k-1 .The forward filtration {F 0,t } t∈[0,t] is analogous to the filtration defined in the prior section; therefore, any continuous semi-martingale adapted to it only considers past events and does not anticipate future events.Reverse time. Consider the backwards Wiener process q w t = w t -w T that is adapted to the backward filtration {F s,T } s∈[0,T ] , then for a continuous semi-martingale {q x t } t∈[0,T ] adapted to the backward filtration, the backward Stratonovich integral isT 0 q x t • d q w t = lim |Π|→0 N k=1 q x t k + q x t k-1 2 ( q w t k-1 -q w t k ) (D.4)The backward filtration {F s,T } s∈[0,T ] is the opposite of the forward filtration in the sense that continuous semi-martingales adapted to it only depend on future events and do not anticipate past events. As such, time is effectively reversed.Remark D.1. While the Stratonovich symmetric integrals give us a powerful tool for integrating forwards and backwards in time with stochastic integrals, it is important that we use the same realization of the Wiener process.Consider the Stratonovich SDE defined asx T = x 0 + T 0 f (x t , t) dt + T 0 g(x t , t) • dw t , (D.5)where f , g ∈ C ∞,1 b , i.e., they belong to the class of functions with infinitely many bounded derivatives w.r.t. the state and bounded first derivatives w.r.t. time. Thus, the SDE has a unique strong solution. Given a realization of the Wiener process, there exists a smooth mapping Φ called the stochastic flow such that Φ s,t (x s ) is the solution at time t of the process described in Equation (D.5) started at x s at time s ≤ t. This then defines a collection of continuous maps S = {Φ s,t } s≤t;s,t∈[0,T ] from R d to itself.Kunita [51, Theorem 3.7.1] shows that with probability 1 this collection S satisfies the flow propertyΦ s,t (x s ) = Φ u,t (Φ s,u (x s )) s ≤ u ≤ t, x s ∈ R d , (D.6)and that each Φ s,t is a smooth diffeomorphism from R d to itself. Hence, S is the stochastic flow of diffeomorphisms generated by Equation (D.5). Moreover, the backward flow q Ψ s,t := Φ -1 s,t satisfies the backwards SDE:EQUATIONfor all s, t ∈ [0, T ] such that s ≤ t. This formulation makes intuitive sense as the backwards SDE differs only from the forwards SDE by a negative sign.Now consider the adjoint flow A s,t (x s ) = ∂L(Φ s,t (x s ))/∂x s , then q A s,t (x t ) = A s,t ( q Ψ s,t (x t )). Li et al. [50] show that q A s,t (x t ) satisfies the backward SDE:q A s,t (x t ) = ∂L ∂x t + t s q A u,t (x t ) ∂f ∂x u ( q Ψ u,t (x t ), u) du + t s q A u,t (x t ) ∂g ∂x u ( q Ψ u,t (x t ), u) • d q w u .(D.8) As the drift and diffusion coefficient of this SDE are in C ∞,1 b , the system has a unique strong solution.We are now ready to put all of this together to prove the result from the main paper.Proof.dx t = f (x t , t) dt + g(t) • dw t .(D.9) By Equation (D.8) the adjoint state admitted by the flow of diffeomorphisms generated by Equation (D.9) evolves with the SDEq A s,t (x t ) = ∂L ∂x t + t s q A u,t (x t ) ∂f ∂x u ( q Ψ u,t (x t ), u) du + t s q A u,t (x t ) ∂g ∂x u (u) • d q w u =0 = ∂L ∂x t + t s q A u,t (x t ) ∂f ∂x u ( q Ψ u,t (x t ), u) du. (D.10)Clearly, the adjoint state evolves with an ODE revolving around only the drift coefficient, i.e., f . Therefore, we can rewrite the evolution of the adjoint state as da x (t) = -a x (In this section, we include some additional experiments that did not fit within the main paper. One of the advantages of AdjointDEIS is that the solver for the diffusion ODE and continuous adjoint equations are distinct. This means that we don't have to force N = M enabling greater flexibility when using AdjointDEIS. As such we explore the impact of using fewer steps to estimate the gradient while keeping the number of sampling steps N = 20 fixed. In Figure 4 we illustrate the impact of the change in the number of discretization steps when estimating the gradients. Unsurprisingly, the fewer steps we take the less accurate the gradients are. This matches the empirical data presented in Table 2 which measures the impact of face morphing performance measured in MMPMR. We measure the impact of the learning rate on guided generation with AdjointDEIS in Table 2 .Unsurprisingly, large learning rates lower performance, especially for less accurate gradients. I.e., when M is small. We illustrate an example of the impact in Figure 5 . Clearly, the learning rate of η = 1 starts to distort the images even if it still fools the FR system. As alluded to in the main paper, one of the drawbacks of diffusion SDEs is that they require small step sizes to work properly. We observe that the missing high frequency content is added back in when the step size is increased, see Figure 6 .F.1 AdjointDEIS-2M AlgorithmFor completeness we have the full AdjointDEIS-2M solver implemented in Algorithm 1.Algorithm 1 AdjointDEIS-2M.Require: Initial values a x (0), monotonically increasing time steps {t i } M i=0 , and noise prediction model ϵ θ (x t , z, t). 1: Denote h i := λ ti+1 -λ ti , for i = 0, . . . , M -1.2: ãx (t 0 ) ← a x (0)▷ Initialize an empty buffer Q.3: ãz (t 0 ) ← 0, ãθ (t 0 ) ← 0. 4: Q buffer ←-ãaug (t 0 ) 5: ãx (t 1 ) ← αt i αt 1 ãx (t 0 ) + σ t1 (e h0 -1) α 2 t 0 α 2 t 1 ãx (t 0 ) ⊤ ∂ϵ θ (xt 0 ,z,t) ∂ xt 0 6: ãz (t 1 ) ← ãz (t 0 ) + σ t1 (e h0 -1) αt 0 αt 1 ãx (t 0 ) ⊤ ∂ϵ θ (xt 0 ,z,t) ∂z 7: ãθ (t 1 ) ← ãθ (t 0 ) + σ t1 (e h0 -1) αt 0 αt 1 ãx (t 0 ) ⊤ ∂ϵ θ (xt 0 ,z,t) ∂θ 8: Q buffer ←-ãaug (t 1 ) 9: for i ← 1, 2, . . . , M -1 do 10: ρ i ← hi-1 hi 11: D i ← 1 + 1 2ρi Ṽ(x; t i ) -1 2ρi Ṽ(x; t i-1 )12:E i ← 1 + 1 2ρi Ṽ(z; t i ) -1 2ρi Ṽ(z; t i-1 )13: if i < M -1 then 18:F i ← 1 + 1 2ρi Ṽ(θ; t i ) -Q buffer ←-ãaug (t i+1 ) 19:end if 20: end for 21: return ãx (t M ), ãz (t M ), ãθ (t M ).Our code for AdjointDEIS will soon be available here at https://github.com/zblasingame/ AdjointDEIS.For reproducibility purposes, we provide a list of links to the official repositories of other works used in this paper.3. The ElasticFace model can be found at https://github.com/fdbtrs/ElasticFace. 4. The AdaFace model can be found at https://github.com/mk-minchul/AdaFace. 5. The official Diffusion Autoencoders repository can be found at https://github.com/ phizaz/diffae. 6. The official MIPGAN repository can be found at https://github.com/ ZHYYYYYYYYYYYY/MIPGAN-face-morphing-algorithm.In this section, we outline the details for the experiments run in Section 5.For completeness, we provide the DiM algorithm from [34] following the notation used in [52] . The original bona fide images are denoted x 0 . The conditional encoder is E : X → Z, Φ is the numerical diffusion ODE solver, Φ + is the numerical diffusion ODE solver as time runs forwards from 0 to T . The algorithm is presented in Algorithm 2.Algorithm 2 DiM Framework.Require: Blend parameter w = 0.5. Time schedule{t i } N i=1 ⊆ [0, T ], t i < t i+1 . 1: z a ← E(x (a) 0 )▷ Encoding bona fides into conditionals.2:z b ← E(x (b) 0 ) 3: for i ← 1, 2, . . . , N -1 do 4: x (a) ti+1 ← Φ + (x (a) ti , ϵ θ (x (a)ti , z a , t i ), t i ) ▷ Solving the probability flow ODE as time runs from 0 to T . 5: x (b) ti+1 ← Φ + (x (b) ti , ϵ θ (x (b) ti , z b , t i ), t iIn our reporting of the NFE we record the number of times the diffusion noise prediction U-Net is evaluated both during the encoding phase, N E , and solving of the PF-ODE or diffusion SDE, N . We chose to report N + N E over N + 2N E as even though two bona fide images are encoded resulting in 2N E NFE during encoding, this process can simply be batched together, reducing the NFE down to N E . When reporting the NFE for the Morph-PIPE model, we report N E + BN where B is the number of blends. While a similar argument can be made that the morphed candidates could be generated in a large batch of size B, reducing the NFE of the sampling process down to N , we chose to report BN as the number of blends, B = 21, used in the Morph-PIPE is quite large, potentially resulting in Out Of Memory (OOM) errors, especially if trying to process a mini-batch of morphs.Using N E + N reporting over N E + BN , the NFE of Morph-PIPE is 350, which is comparable to DiM. The reporting of NFE for AdjointDEIS was calculated as N E + n opt (N + M ) where n o pt is the number of optimization steps and M is the number of discretization steps for the continuous adjoint equations.All of the main experiments were done on a single NVIDIA Tesla V100 32GB GPU. On average the guided generation experiments for our approach took between 6 -8 hours for the whole dataset of face morphs with a batch size of 8. Some additional follow-up work for the camera-ready version used an NVIDIA H100 Tensor Core 80GB GPU with a batch size of 16.The SYN-MAD 2022 dataset is derived from the Face Research Lab London (FRLL) dataset [45] . FRLL is a dataset of high-quality captures of 102 different individuals with frontal images and neutral lighting. There are two images per subject, an image of a "neutral" expression and one of a "smiling" expression. The ElasticFace [46] FR system was used to select the top 250 most similar pairs, in terms of cosine similarity, of bona fide images for both genders, resulting in a total of 489 bona fide image pairs for face morphing [44] , as some pairs did not generate good morphs on the reference set; we follow this minimal subset.All three FR systems use the Improved ResNet (IResNet-100) architecture [53] as the neural net backbone for the FR system. The ArcFace model is a widely used FR system [34, 36, 38, 42] . It employs an additive angular margin loss to enforce intra-class compactness and inter-class distance, which can enhance the discriminative ability of the feature embeddings [41] . ElasticFace builds upon the ArcFace model by using an elastic penalty margin over the fixed penalty margin used by ArcFace. This change results in an FR system with state-of-the-art performance [46] . Lastly, the AdaFace model employs an adaptive margin loss by weighting the loss relative to an approximation of the image quality [47] . The image quality is approximated via feature norms and is used to give less weight to misclassified images, reducing the impact of "low" quality images on training. This improvement allows the AdaFace model to achieve state-of-the-art performance in FR tasks.The AdaFace and ElasticFace models are trained on the MS1M-ArcFace dataset, whereas the ArcFace model is trained on the MS1M-RetinaFace dataset. N.B., the ArcFace model used in the identity loss is not the same ArcFace model used during evaluation. The model used in the identity loss is an IResNet-100 trained on the Glint360k dataset [54] with the ArcFace loss. We use the cosine distance to measure the distance between embeddings from the FR models. All three FR systems require images of 112 × 112 pixels. We resize every image, post alignment from dlib which ensures the images are square, to 112 × 112 using bilinear down-sampling. The image tensors are then normalized such that they take values in [-1, 1]. Lastly, the AdaFace FR system was trained on BGR images so the image tensor is shuffled from the RGB format to the BGR format.In this section we compare several recent methods for training-free guided generation. We broadly classify these techniques into two categories:1. Techniques which directly optimize the solution trajectory during sampling [52, 55, 56] 2. Techniques which search for the optimal latents x T and or z (this can include optimizing the solution trajectory as well) [57, 58] .In Table 4 we compare several different techniques for training-free guided diffusion along this category along with whether the formulation is for diffusion ODEs, SDEs, or both. where u is the control function. This optimization objective learns to alter the flow, f θ , by u. In practice, this amounts to injecting a control step governed by u(t) for a discretized schedule. This technique does not allow for learning an optimal x T , z, or θ.FreeDoM [55] The technique would work for either diffusion ODEs or SDEs.DOODL [57] looks at gradient calculation based on the invertibility of EDICT [60] . This method can find the gradient w.r.t. x T ; however, it cannot for the other quantities. DOODL additionally has further overhead due to the dual diffusion process of EDICT. Further analysis of DOODL compared to continuous adjoint equations for diffusion models can be found in [58] . More closely related to our work is the recent AdjointDPM [58] who also explores the use of adjoint sensitivity methods for backpropagation through the probability flow ODE. While they also propose to use the continuous adjoint equations to find gradients for diffusion models, our work differs in several ways which we enumerate in Table 5 . For clarity, we use orange to denote their notation. They reparameterize Equation (2.4) as dy dρ = εθ (e γ -1 (ρ) 0 f (τ ) dτ y, γ -1 (ρ), c) (H.12)where c denotes the conditional information, ρ = γ(t), and dγ dt = e -t 0 f (τ ) dτ g 2 (t) 2σt . which gives them the following ODE for calculating the adjoint. In our approach we integrate over λ t whereas they integrate over ρ. Moreover, we provide custom solvers designed specifically for diffusion ODEs instead of using a black box ODE solver. Our approach is also interoperable with other forward ODE solvers meaning our AdjointDEIS solver is agnostic to the ODE solver used to generate the output; however, the AdjointDPM model is tightly coupled to its forward solver. Lastly and most importantly our method is more general and supports diffusion SDEs, not just ODEs.Although the original paper omitted a closed form expression for γ -1 (ρ), we provide to give a comparison between both methods, and to ensure AdjointDPM can be fully implemented. In the VP SDE scheme with a linear noise schedule log α t is found to belog α t = - β 1 -β 0 4 t 2 - β 0 2 t (H.14)on t ∈ [0, 1] with β 0 = 0.1, β 1 = 20, following Song et al. [15] . Then γ -1 (ρ) is found to be γ -1 (ρ) = β 0 -β 2 0 + 4 log (β 0 -β 1 )β 0 -β 1 (H.15)Concurrent work to ours by Marion et al. [59] has also explored the method of adjoint sensitivity for guidance of diffusion models. They, however, focus on an efficient scheme to parallelize the solution to the adjoint ODE from the perspective of bi-level optimization rather than the adjoint technique itself. So while we focused on the details of the continuous adjoint equations, they focused on an efficient implementation of the optimization problem from the perspective of bi-level optimization.For completeness, we show how to analytically compute the drift and diffusion coefficients for a linear noise schedule Ho et al. [1] in the VP scenario Song et al. [15] . With a linear noise schedule log α t is found to belog α t = - β 1 -β 0 4 t 2 - β 0 2 t (I.1)on t ∈ [0, 1] with β 0 = 0.1, β 1 = 20, following Song et al. [15] . The drift coefficient becomesf (t) = - β 1 -β 0 2 t - β 0 2 (I.2)and as σ t = 1 -α Importantly, dσt dt does not exist at time t = 0, as σ t is discontinuous at that point, and so an approximation is needed when starting from this initial step. In practice, adding a small ϵ ≪ 1 to t = 0 should suffice.A more precise and technical name would be the empirical adjoint probability flow ODE; however, for brevity's sake we prefer the adjoint diffusion ODE and leave the reader to infer from context that this describes an empirical model with learned score function.The notation '• dwt' denotes Stratonovich integration which is different from Itô. More details on this are found in Appendix D.The SYN-MAD 2022 dataset used in this paper can be found at https://github.com/ marcohuber/SYN-MAD-2022.The ArcFace models, MS1M-RetinaFace dataset, and MS1M-ArcFace dataset can be found at https://github.com/deepinsight/insightface. 


Decision: This work leverages adjoint sensitivity analysis to efficiently compute gradients of the probability flow ODE implementation of diffusion model, in order to optimize latent variables and hyperparameters. Reviewers and AC agreed about the merits of this work, and acceptance is recommended.

Reviews: review: paper_summary: paper_summary: AdjointDEIS uses the method of adjoint sensitivity to compute gradients of diffusion models, which is more efficient and less memory intensive, and robust to the injected noise. This work proposes efficient solvers for both the adjoint probability flow ODE and the adjoint diffusion SDE. Experiments demonstrate the efficacy of the solvers on guided generation for a face morphing problem.
strength_weakness: strength_weakness: - The number of experimental applications of the method seem lacking (just one application to face morphing attacks is included). I’d encourage the authors to think about applications to data assimilation or inverse problems for instance.
- Since the primary contribution of the paper is solver, having a proof of convergence rates or empirical studies of convergence rate would be valuable additions to the paper.
questions: questions: - It would be interesting to visually plot example forward and backward trajectories of the diffusion processes. What is the typical magnitude of the numerical errors that accumulate during the backward ODE solve?
limitations: limitations: N/A
paper_summary: paper_summary: The paper proposes an adjoint sensitivity method -- AdjointDEIS -- for efficiently calculating gradients of diffusion SDE models. Current methods for naive backpropagation rely on discrete adjoints which are memory intensive. The authors introduce an approach based on the stochastic adjoint sensitivity method to solve the gradients with respect to initial noise, conditional information, and model parameters. The authors develop custom solvers for the adjoint problems along with the proposed sensitivity method. The methods are validated on a face morphing problem.
strength_weakness: strength_weakness: * **Limited Experimental Scope**: While the overall analysis and formulations are sound, the overall experimental validation for the method is limited to just a single experiment.
questions: questions: 1. Point 2 in contributions: Aren't the continuous adjoint sensitivity methods used in Neural SDEs already general-purpose enough to handle diffusion problems?
limitations: limitations: N/A
paper_summary: paper_summary: In this paper, the authors proposed an accelerated  method for differentials of pretrained diffusion models with respect to its latent valuables or parameters by making use of 
1. the Taylor expansion of the log-SNR parameter, and
2. the exact integral formula of the derivatives related to the probability flow ODE, which they call Adjoint Diffusion ODE.

They applied this method to morphing of two generated images. As they summarized in Table 1, their results achieved good performance.

For more details:
- The setting
    - Data diffusion SDE $d\mathbf{x}_t = f(t) \mathbf{x}_t + g(t) d\mathbf{w}_t$
        - Its integral $q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t|\alpha_t \mathbf{x}_0, \sigma_t^2)$ where $f(t), g(t)$ are related to $\alpha_t, \sigma_t$ by Eq.(2.5).
    - A pretrained diffusion model $\epsilon_\theta(\mathbf{x},\mathbf{z}, t)$
        - Its generating process: the probability flow ODE $\frac{d\mathbf{x}\_t}{dt} = f(t)\mathbf{x}\_t + \frac{g(t)^2}{2 \sigma_t}\epsilon_\theta(\mathbf{x}, \mathbf{z}, t)$
        - They name the velocity (the RHS of the ODE) as $\mathbf{h}_\theta(\mathbf{x}, \mathbf{z}, t)$
- The purpose
    - Calculation of gradients of $\mathcal{L}(\mathbf{x}_0(\mathbf{x}_T, \mathbf{z}, \theta))$ with respect to $\mathbf{z}$ or $\theta$
    - which are defined by ODEs in Eq.(3.3) and (3.4).
- The method
    - One can recover the derivatives based on a vector $\hat{\mathbf{a}}_t = \frac{\partial \mathcal{L}}{\partial \mathbf{y}_t}$ where $\mathbf{y}_t = \alpha_t^{-1}\mathbf{x}_t$.
    - For $t<s$, one can get $\hat{\mathbf{a}}_s = \hat{\mathbf{a}}_t + \int\_{\lambda_t}^{\lambda_s} ... d\lambda$ (Proposition 3.1).
    - Applying the Taylor expansion on ... with $\lambda$ around $\lambda_t$, and drop $O((\lambda_s - \lambda_t)^{k+1})$ term, we get a concrete relation between $\hat{\mathbf{a}}_s$ and $\hat{\mathbf{a}}_t$ (Eq.(3.11)),
    - which provides explicit (approximated) formulas for the derivatives.
- The empirical result
    - They apply the proposed method to morphing (Section 5 Experiments).
strength_weakness: strength_weakness: - **significance(-)** It is difficult to understand in which respect the proposed method is superior.
- **clarity(-)** I am almost convinced by the explanations, however there is a concern as explained in Questions.
questions: questions: - for significance improvements:
    - **Q1(significance)** In the theoretical sections, the authors apply omitting the higher order terms of the Taylor expansion, so the resultant method should be just an approximation. Why is it good performance in Table 1 compared to other methods nonetheless the method is just an approximation?
- for clarity improvements:
    - **Q1(clarity)** In the Adjoint Empirical Diffusion SDE cases, the integral results in Eq.(4.12) and (4.13) look very similar to ODE cases in Eq.(3.13) and (3.14). The only difference seem to be (doubled) 2nd term. Why such similarity happens?
    - **Q2(clarity)** In the beginning of section 4, the authors wrote their motivation to analyze SDE cases as it possibly improve the result, however, it turns to be that the SDE version works worse than the ODE cases. Why does it happens?
limitations: limitations: They leave comments on
1. Lack of conditional generation in numerical experiments
2. No application of parameter differentiation in numerical experiments

in the final page.
paper_summary: paper_summary: This paper presents a way to optimize the initial seed or control parameters to a diffusion network to optimize a differentiable loss applied to the final samples drawn from a discretization of the diffusion SDE / probability flow ODE. Put differently the idea is to apply the white box integration trick from DPM / DPM++ paper, and the method from the neural ODE paper to differentiate through an ODE solver.

In addition, the paper presents methods for backpropagating through an SDE. 

In terms of experiments the method is demonstrated on a face morphing application.
strength_weakness: strength_weakness: The experimental results seem to be weak.

This paper presents a way to control the final outcome from a diffusion sampler using a differentiable loss, but the only application demonstrated in the paper is on face-morphing where the NFE are not controlled across the baselines. Some of the baselines have fewer NFE and it's not possible to say conclusively whether they underperform because they are worse methods or just because they had less computational budget.

AdjointDEIS-2M is not even implemented and compared. 

Also, a lot more compelling demonstrations could be shown of these new algorithms. For example, a number of "controlnet" type applications could have been tried to create an image that matches a sketch of an output.
questions: questions: Mainly I would like to understand a little bit more about the rationale for selecting the face morph application selection, and whether the authors think the empirical comparisons in Table 1 are apples-to-apples given the different NFE for the methods.


Some more minor comments
- eq (3.1) has a factor of 2 missing.
- around eq(2.5) it might be good to explicitly clarify that while VP typically uses alpha^2 + sigma^2 = 1 ibut you are not making that assumption here.
limitations: limitations: Yes limitations have been addressed.

Comments: title: Response to authors
comment: Thank you for addressing my comment. I have updated my score. I still encourage the authors to do more extensive evaluation for the method in the future.
Thank you for addressing this.

While the overall paper is well written and theoretically sound, I agree with reviewer vPC5 that the experimental scope of this paper is quite limited (with just a single experiment). As such I will keep my recommendation of Accept (score 7)
Thanks for the clear explanations.

On **Q1(significance)**, I understand a little about the reason for good performance that is related to  the capability of direct gradient estimates based on proposed methods in this paper. I think it makes sense.

On **Q1(clarity)**, I understood the reason for appearance of 2. This seems to be in the same reasons of 2 in denominator of the probability flow ODE.

On **Q2(clarity)**, I think it would be making known contributions more clear as the author responded, and I believe it will be improved in the camera-ready version. 

Overall, I think my concerns are solved by the response, and I would like to raise my score a little.