The DDIM Update Rule
Problem Statement
Implement the general DDIM update, with the interpolation parameter Ī· that sweeps continuously from a deterministic ODE solver to full DDPM ancestral sampling.
Background
DDIM (Song et al., 2021) observes that the DDPM objective constrains only the marginals q(xtāā£x0ā), so you are free to pick a non-Markovian reverse process with the same marginals. That yields a family of samplers indexed by Ī·ā„0:
Ļtā=Ī·1āαĖtā1āαĖtā²āāā1āαĖtā²āαĖtāāā
x_{t'} = \underbrace{\sqrt{\bar{\alpha}_{t'}}\,\hat{x}_0}_{\text{predicted } x_0} + \underbrace{\sqrt{1-\bar{\alpha}_{t'}-\sigma_t^2}\;\hat{\varepsilon}}_{\text{direction pointing to } x_t} + \sigma_t z$$ where $t'$ is the *next* (smaller) timestep, not necessarily $t-1$. Read the middle term carefully: the amount of noise pointed back along $\hat{\varepsilon}$ is reduced by exactly the amount of fresh noise $\sigma_t$ that is about to be injected, so the total variance stays $1-\bar{\alpha}_{t'}$ whatever $\eta$ is. At $\eta = 0$ nothing random enters and the sampler becomes a deterministic map from the initial latent to the image -- which is what makes DDIM invertible and makes latent interpolation meaningful. At $\eta = 1$ with $t' = t-1$ it reduces exactly to DDPM ancestral sampling. ## Your Task Implement: ```python def ddim_step(x_t, eps, alpha_bar_t, alpha_bar_prev, eta=0.0, z=None): ``` Return $x_{t'}$. When **eta == 0** the **z** argument is unused and may be **None**. ## Input Format - **x_t**, **eps**: NumPy arrays of matching shape. - **alpha_bar_t** (float): $\bar{\alpha}$ at the current step, in **(0, 1)**. - **alpha_bar_prev** (float): $\bar{\alpha}$ at the target step, in **(0, 1]**. Use **1.0** for the final step to a clean sample. - **eta** (float): **0.0** is deterministic DDIM, **1.0** is DDPM. - **z**: array shaped like **x_t**, or **None**. ## Output Format A NumPy array shaped like **x_t**. ## Sample ```python x_t = np.array([0.4, -0.6]) eps = np.array([0.1, 0.2]) print(np.round(ddim_step(x_t, eps, 0.5, 0.8, eta=0.0), 4).tolist()) ``` **x0_hat = (x_t - sqrt(0.5)*eps)/sqrt(0.5)**, **sigma = 0** so the direction term is **sqrt(1-0.8)*eps**, and the output is **sqrt(0.8)*x0_hat + sqrt(0.2)*eps**.Example:
x_t = np.array([0.4, -0.6]) eps = np.array([0.1, 0.2]) print(np.round(ddim_step(x_t, eps, 0.5, 0.8, eta=0.0), 4).tolist())
[0.4612, -0.8484]
With eta = 0 the sigma term vanishes entirely. x0_hat = (x_t - sqrt(0.5)*eps)/sqrt(0.5) = [0.4657, -1.0485]. The direction term is sqrt(1 - 0.8 - 0) * eps = sqrt(0.2)*eps. Adding sqrt(0.8)*x0_hat to it gives the result. Note that the schedule moves up here (0.5 to 0.8), which is fine: DDIM only needs the two alpha-bar values.
Constraints:
- Clamp the argument of the direction-term square root at zero: use
max(1 - alpha_bar_prev - sigma**2, 0.0)to absorb floating-point negatives. - Add the
sigma * zterm only wheneta > 0. alpha_bar_prev = 1.0is legal and must yield exactly x^0ā (sigma and the direction term both vanish).- Do not draw randomness;
zis supplied. - Do not round inside the function.
1. Background Knowledge
Diffusion models generate data by reversing a gradual noising process. The original DDPM (Denoising Diffusion Probabilistic Models) defines a Markovian reverse process where each step xtā1ā depends only on xtā and a predicted noise ε^. However, DDIM (Denoising Diffusion Implicit Models) reveals that the training objective only constrains the marginal distribution q(xtāā£x0ā), not the conditional transition q(xtā1āā£xtā). This freedom allows us to define a non-Markovian reverse process that shares the same marginals but behaves differently.
The key innovation in DDIM is the introduction of an interpolation parameter Ī·ā„0. When Ī·=0, the reverse process becomes deterministic, effectively solving an ordinary differential equation (ODE). This determinism makes the process invertible, enabling applications like latent interpolation and image editing. When Ī·=1 and steps are unit-sized (tā²=tā1), the update rule reduces exactly to the standard DDPM ancestral sampling. For 0<Ī·<1, the sampler interpolates between these two extremes, allowing for faster sampling with fewer steps while maintaining high sample quality.
The update rule relies on predicting the clean data x^0ā from the noisy input xtā and the predicted noise ε^. The next state xtā²ā is constructed as a linear combination of x^0ā and ε^, plus a controlled amount of fresh noise Ļtāz. The coefficient Ļtā is carefully designed so that the total variance of the resulting distribution matches the required marginal variance 1āαĖtā²ā, regardless of the value of Ī·. This ensures that the generated samples remain consistent with the learned data distribution.
2. Algorithm Approach
The core task is to implement the vectorized update equation provided in the problem statement. The approach involves three main computational stages:
- Predict Clean Data: Calculate x^0ā using the current noisy sample xtā, the predicted noise ε^, and the current noise schedule value αĖtā.
- Calculate Noise Scale: Compute Ļtā based on Ī·, αĖtā, and αĖtā²ā. This term controls how much randomness is injected into the next step.
- Construct Next State: Combine the predicted clean data, the direction towards the current noisy sample (scaled by the remaining variance), and the fresh noise term.
The implementation must handle the special case where Ī·=0 (deterministic), in which case the noise term z is irrelevant and Ļtā becomes zero. The function should be robust to None inputs for z when Ī·=0.
3. Step-by-Step Strategy
- Compute x^0ā: Use the formula x^0ā=αĖtāāxtāā1āαĖtāāε^ā.
- Calculate 1āαĖtāā and αĖtāā.
- Perform element-wise operations to get x^0ā.
- Compute Ļtā: Use the formula Ļtā=Ī·1āαĖtā1āαĖtā²āāā1āαĖtā²āαĖtāāā.
- Note: If Ī·=0, Ļtā is simply 0.
- Be careful with the order of operations and square roots.
- Handle the case where αĖtā=αĖtā²ā (though typically tā²<t implies αĖtā²ā>αĖtā in standard schedules, the math holds).
- Compute the Direction Term: The term 1āαĖtā²āāĻt2āāε^ represents the direction pointing back towards xtā adjusted for the new noise.
- Calculate the coefficient 1āαĖtā²āāĻt2āā.
- Multiply this coefficient by ε^ (element-wise).
- Compute the Noise Term:
- If Ī·>0, ensure z is provided. Multiply Ļtā by z.
- If Ī·=0, this term is zero.
- Combine Terms: Calculate xtā²ā=αĖtā²āāx^0ā+DirectionĀ Term+NoiseĀ Term.
- Multiply x^0ā by αĖtā²āā.
- Add the direction term.
- Add the noise term.
- Return Result: Return the resulting array xtā²ā.
4. Common Pitfalls
- Incorrect Variance Calculation: A common mistake is misinterpreting the term 1āαĖtā²āāĻt2āā. This term ensures that the sum of variances from the predicted mean and the noise equals the target marginal variance 1āαĖtā²ā. Do not confuse this with 1āαĖtā²āā.
- Handling Ī·=0: When Ī·=0, Ļtā becomes 0. The code should not crash if z is None. Ensure you check for eta == 0 or handle None gracefully before attempting to multiply z.
- Square Root Domains: Ensure that the arguments to np.sqrt are non-negative. Due to floating-point precision, values like 1āαĖtā might become slightly negative if αĖtā is very close to 1. Using np.maximum(0, value) inside sqrt can prevent NaNs.
- Broadcasting Issues: Ensure that scalar values like α_bar_t are broadcast correctly against the array shapes of xtā and ε. NumPy usually handles this, but explicit scalar multiplication is safer.
- Confusing t and tā²: Remember that tā² is the previous (smaller) timestep in the reverse process, meaning αĖtā²ā>αĖtā (since Ī±Ė decreases as t increases). The formula uses αĖtā²ā for the target state's scale.
5. Time & Space Complexity
- Time Complexity: O(N), where N is the total number of elements in the input array xtā. The operations involve element-wise arithmetic (addition, multiplication, square roots) across the entire array. There are no nested loops or complex dependencies between elements.
- Space Complexity: O(N), to store the intermediate arrays x^0ā, the direction term, and the final result xtā²ā. If implemented in-place or with careful memory management, it could be reduced, but typically temporary arrays of size N are created.