Diffusion Schedulers Library (DDIM, PNDM, DPM-Solver++)

[ooples]

Diffusion Schedulers Library (within #261 umbrella)

Scope

• Provide generic schedulers (IStepScheduler) implemented with project numeric ops and collections
• Initial implementations: DDIMScheduler (done), plus at least one of PNDM / DPM-Solver++ next
• Unit tests across net462/net8.0; coverage ≥90% for new code

Integration

• Schedulers are injected into diffusion models (e.g., DDPMModel) via constructor and are not a standalone builder option.

Umbrella

• Part of Introduce Diffusion Models (Core) – Part 1 #261 (Introduce Diffusion Models – Part 1)

---

⚠️ CRITICAL ARCHITECTURAL REQUIREMENTS

Before implementing this user story, you MUST review:

• 📋 Full Requirements: .github/USER_STORY_ARCHITECTURAL_REQUIREMENTS.md
• 📐 Project Rules: .github/PROJECT_RULES.md

Mandatory Implementation Checklist

1. INumericOperations Usage (CRITICAL)

• Include protected static readonly INumericOperations<T> NumOps = MathHelper.GetNumericOperations<T>(); in base class
• NEVER hardcode double, float, or specific numeric types - use generic T
• NEVER use default(T) - use NumOps.Zero instead
• Use NumOps.Zero, NumOps.One, NumOps.FromDouble() for values
• Use NumOps.Add(), NumOps.Multiply(), etc. for arithmetic
• Use NumOps.LessThan(), NumOps.GreaterThan(), etc. for comparisons

2. Inheritance Pattern (REQUIRED)

• Create I{FeatureName}.cs in src/Interfaces/ (root level, NOT subfolders)
• Create {FeatureName}Base.cs in src/{FeatureArea}/ inheriting from interface
• Create concrete classes inheriting from Base class (NOT directly from interface)

3. PredictionModelBuilder Integration (REQUIRED)

• Add private field: private I{FeatureName}<T>? _{featureName}; to PredictionModelBuilder.cs
• Add Configure method taking ONLY interface (no parameters):

  public IPredictionModelBuilder<T, TInput, TOutput> Configure{FeatureName}(I{FeatureName}<T> {featureName})
  {
      _{featureName} = {featureName};
      return this;
  }

• Use feature in Build() with default: var {featureName} = _{featureName} ?? new Default{FeatureName}<T>();
• Verify feature is ACTUALLY USED in execution flow

4. Beginner-Friendly Defaults (REQUIRED)

• Constructor parameters with defaults from research/industry standards
• Document WHY each default was chosen (cite papers/standards)
• Validate parameters and throw ArgumentException for invalid values

5. Property Initialization (CRITICAL)

• NEVER use default! operator
• String properties: = string.Empty;
• Collections: = new List<T>(); or = new Vector<T>(0);
• Numeric properties: appropriate default or NumOps.Zero

6. Class Organization (REQUIRED)

• One class/enum/interface per file
• ALL interfaces in src/Interfaces/ (root level)
• Namespace mirrors folder structure (e.g., src/Regularization/ → namespace AiDotNet.Regularization)

7. Documentation (REQUIRED)

• XML documentation for all public members
• <b>For Beginners:</b> sections with analogies and examples
• Document all <param>, <returns>, <exception> tags
• Explain default value choices

8. Testing (REQUIRED)

• Minimum 80% code coverage
• Test with multiple numeric types (double, float)
• Test default values are applied correctly
• Test edge cases and exceptions
• Integration tests for PredictionModelBuilder usage

---

⚠️ Failure to follow these requirements will result in repeated implementation mistakes and PR rejections.

See full details: .github/USER_STORY_ARCHITECTURAL_REQUIREMENTS.md

[ooples]

Agent: Agent-Schedulers
Branch: feature/issue-263-diffusion-schedulers-library-ddim-pndm-dpm-solver
Worktree: .worktrees\feature-issue-263-diffusion-schedulers-library-ddim-pndm-dpm-solver

User Story
As a .NET ML engineer, I want to deliver "Diffusion Schedulers Library (DDIM, PNDM, DPM-Solver++)" so that AiDotNet users can rely on a robust, reproducible, and .NET Framework 4.6–compatible implementation.

Acceptance Criteria

• Design reviewed and linked in PR description
• API surface aligned with existing patterns in src and documented in README changes
• Unit tests added (xUnit) covering new logic with >= 90% line/branch coverage for changed files
• Builds succeed for
  et8.0 and
  et462 on CI and locally
• Benchmarks or microperf notes included where performance-sensitive (if applicable)
• Examples or docs updated where applicable

Plan of Action

• Confirm scope and traceability to this issue
• Design draft and API sketch (commit as docs or comment)
• Implement feature behind tests
• Add tests (xUnit) and ensure coverage >= 90%
• Run dotnet build + dotnet test --collect:"XPlat Code Coverage" for net8.0 and net462
• Update docs (README, samples) if needed
• Open PR linking this issue ("Fixes Diffusion Schedulers Library (DDIM, PNDM, DPM-Solver++) #263")

Definition of Done

• All acceptance criteria met
• Build green on CI
• Code coverage >= 90% for new/changed code
• Maintains compatibility with .NET Framework 4.6 (
  et462)

[ooples]

Junior Developer Implementation Guide: Issue #263

Diffusion Schedulers Library (DDIM, PNDM, DPM-Solver++)

Overview

This guide extends Issue #261's scheduler foundation by implementing advanced samplers that enable faster, higher-quality generation. These schedulers are the "secret sauce" that make models like Stable Diffusion practical for real-world use.

---

Understanding Advanced Schedulers

Why Multiple Schedulers Matter

The Speed-Quality Tradeoff:

• DDPM (original): 1000 steps, slow but thorough
• DDIM: 50 steps, 20x faster, similar quality
• PNDM: 20-30 steps, even faster
• DPM-Solver++: 10-20 steps, FASTEST with high quality

Real-World Impact:

• DDPM: 5 minutes per image
• DDIM: 15 seconds per image
• DPM-Solver++: 5 seconds per image

Analogy: Cleaning a dirty window:

• DDPM: Wipe 1000 times with tiny circles
• DDIM: Wipe 50 times with bigger strokes
• PNDM: Use momentum from previous wipes (4-step history)
• DPM-Solver++: Solve the optimal cleaning trajectory mathematically

---

Scheduler Comparison

Scheduler	Steps	Speed	Quality	Deterministic	Method
DDPM	1000	Slow	Excellent	No (stochastic)	Markov chain
DDIM	20-100	Fast	Excellent	Yes (eta=0)	Implicit model
PNDM	20-50	Very Fast	Excellent	Yes	Pseudo-numerical ODE
DPM-Solver++	10-25	FASTEST	Excellent	Yes	Second-order ODE solver

---

Implementation Guide

Step 1: PNDM Scheduler (Pseudo Numerical Methods)

Key Concepts

What is PNDM?

• Uses past 4 steps to predict next step (like momentum)
• Solves the ODE more accurately than DDIM
• Particularly good for few-step sampling (20-30 steps)

The Math:
PNDM uses a linear multistep method that considers the trajectory history:

x_{t-1} = x_t + Δt * (55*f_t - 59*f_{t-1} + 37*f_{t-2} - 9*f_{t-3}) / 24

Where f_t is the model prediction at timestep t.

Real-World Analogy:
Imagine driving a car and predicting where you'll be next:

• DDIM: Look at current speed
• PNDM: Consider current speed + how you've been accelerating/decelerating

File: src/Diffusion/Schedulers/PNDMScheduler.cs

using AiDotNet.Helpers;
using AiDotNet.Interfaces;
using AiDotNet.Mathematics;

namespace AiDotNet.Diffusion.Schedulers;

/// <summary>
/// Pseudo Numerical Methods for Diffusion Models (PNDM) scheduler.
/// </summary>
/// <typeparam name="T">Numeric type for calculations.</typeparam>
/// <remarks>
/// <b>For Beginners:</b> PNDM uses past predictions to improve accuracy.
///
/// Key innovation: Linear multistep method
/// - Remembers last 4 model predictions
/// - Uses them to compute better trajectory
/// - Fewer steps needed than DDIM for same quality
///
/// Advantages:
/// - 20-30 steps for high quality (vs 50 for DDIM)
/// - Deterministic (reproducible)
/// - Particularly good for fast sampling
///
/// Tradeoffs:
/// - Needs 4-step warmup (first 4 steps use different formula)
/// - Slightly more memory (stores history)
/// - More complex than DDIM
///
/// From paper: "Pseudo Numerical Methods for Diffusion Models on Manifolds" (Liu et al., 2022)
/// </remarks>
public class PNDMScheduler<T> : StepSchedulerBase<T>
{
    private readonly int _order; // Multistep order (default: 4)
    private readonly List<Tensor<T>> _etsCache; // Model outputs history

    /// <summary>
    /// Initializes PNDM scheduler.
    /// </summary>
    /// <param name="config">Scheduler configuration.</param>
    /// <param name="order">
    /// Multistep order. Default: 4.
    /// Higher order = more accurate but needs more warmup steps.
    /// Valid values: 1-4.
    /// </param>
    /// <remarks>
    /// <b>For Beginners:</b> Order controls how much history to use.
    ///
    /// - order=1: Like DDIM (no history)
    /// - order=4: Uses last 4 predictions (recommended)
    ///
    /// Higher order improves quality but requires warmup steps.
    /// </remarks>
    public PNDMScheduler(SchedulerConfig<T> config, int order = 4)
        : base(config)
    {
        if (order < 1 || order > 4)
            throw new ArgumentException("Order must be between 1 and 4", nameof(order));

        _order = order;
        _etsCache = new List<Tensor<T>>();
    }

    public override Tensor<T> Step(Tensor<T> modelOutput, int timestep, Tensor<T> sample)
    {
        if (timestep < 0 || timestep >= Config.NumTrainTimesteps)
            throw new ArgumentOutOfRangeException(nameof(timestep));

        // Add current prediction to cache
        _etsCache.Add(modelOutput);

        // Keep only last 'order' predictions
        if (_etsCache.Count > _order)
        {
            _etsCache.RemoveAt(0);
        }

        // Use appropriate stepping method based on cache size
        if (_etsCache.Count < _order)
        {
            // Warmup phase: use simple method
            return StepPrimeIteration(modelOutput, timestep, sample);
        }
        else
        {
            // Main phase: use multistep method
            return StepPlmsIteration(timestep, sample);
        }
    }

    /// <summary>
    /// Prime iteration for warmup steps.
    /// </summary>
    /// <remarks>
    /// <b>For Beginners:</b> First few steps use simpler formula.
    ///
    /// We don't have enough history yet, so fall back to DDIM-like step.
    /// </remarks>
    private Tensor<T> StepPrimeIteration(Tensor<T> modelOutput, int timestep, Tensor<T> sample)
    {
        // Get alpha values
        T alphaCumprod = AlphasCumprod[timestep];
        T alphaCumprodPrev = timestep > 0 ? AlphasCumprod[timestep - 1] : NumOps.One;

        // Predict original sample (x0)
        T sqrtAlphaCumprod = NumOps.Sqrt(alphaCumprod);
        T sqrtOneMinusAlphaCumprod = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprod));

        Tensor<T> predOriginalSample = sample
            .Subtract(modelOutput.Multiply(sqrtOneMinusAlphaCumprod))
            .Divide(sqrtAlphaCumprod);

        // Compute previous sample using DDIM formula
        T sqrtAlphaCumprodPrev = NumOps.Sqrt(alphaCumprodPrev);
        T sqrtOneMinusAlphaCumprodPrev = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprodPrev));

        Tensor<T> predSample = predOriginalSample.Multiply(sqrtAlphaCumprodPrev);
        Tensor<T> noiseTerm = modelOutput.Multiply(sqrtOneMinusAlphaCumprodPrev);

        return predSample.Add(noiseTerm);
    }

    /// <summary>
    /// PLMS (Pseudo Linear Multistep) iteration using history.
    /// </summary>
    /// <remarks>
    /// <b>For Beginners:</b> Main stepping method using past predictions.
    ///
    /// Formula (4th order):
    /// Combined prediction = (55*f_t - 59*f_{t-1} + 37*f_{t-2} - 9*f_{t-3}) / 24
    ///
    /// Weights:
    /// - Current step (f_t): 55/24 = 2.29
    /// - Previous step (f_{t-1}): -59/24 = -2.46
    /// - Two steps back (f_{t-2}): 37/24 = 1.54
    /// - Three steps back (f_{t-3}): -9/24 = -0.38
    ///
    /// This weighted combination gives a better trajectory estimate.
    /// </remarks>
    private Tensor<T> StepPlmsIteration(int timestep, Tensor<T> sample)
    {
        // Get the last 'order' predictions
        var ets = _etsCache.TakeLast(_order).ToList();

        // PLMS coefficients for 4th order
        T[] coefficients = _order switch
        {
            1 => new[] { NumOps.One },
            2 => new[]
            {
                NumOps.FromDouble(3.0 / 2.0),  // 1.5
                NumOps.FromDouble(-1.0 / 2.0)  // -0.5
            },
            3 => new[]
            {
                NumOps.FromDouble(23.0 / 12.0),  // 1.917
                NumOps.FromDouble(-16.0 / 12.0), // -1.333
                NumOps.FromDouble(5.0 / 12.0)    // 0.417
            },
            4 => new[]
            {
                NumOps.FromDouble(55.0 / 24.0),  // 2.292
                NumOps.FromDouble(-59.0 / 24.0), // -2.458
                NumOps.FromDouble(37.0 / 24.0),  // 1.542
                NumOps.FromDouble(-9.0 / 24.0)   // -0.375
            },
            _ => throw new InvalidOperationException($"Order {_order} not supported")
        };

        // Combine predictions with PLMS coefficients
        Tensor<T> combinedPrediction = ets[0].Multiply(coefficients[0]);

        for (int i = 1; i < ets.Count; i++)
        {
            Tensor<T> weightedPrediction = ets[i].Multiply(coefficients[i]);
            combinedPrediction = combinedPrediction.Add(weightedPrediction);
        }

        // Apply to get next sample
        T alphaCumprod = AlphasCumprod[timestep];
        T alphaCumprodPrev = timestep > 0 ? AlphasCumprod[timestep - 1] : NumOps.One;

        // Predict x0 using combined prediction
        T sqrtAlphaCumprod = NumOps.Sqrt(alphaCumprod);
        T sqrtOneMinusAlphaCumprod = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprod));

        Tensor<T> predOriginalSample = sample
            .Subtract(combinedPrediction.Multiply(sqrtOneMinusAlphaCumprod))
            .Divide(sqrtAlphaCumprod);

        // Compute previous sample
        T sqrtAlphaCumprodPrev = NumOps.Sqrt(alphaCumprodPrev);
        T sqrtOneMinusAlphaCumprodPrev = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprodPrev));

        Tensor<T> predSample = predOriginalSample.Multiply(sqrtAlphaCumprodPrev);
        Tensor<T> noiseTerm = combinedPrediction.Multiply(sqrtOneMinusAlphaCumprodPrev);

        return predSample.Add(noiseTerm);
    }

    /// <summary>
    /// Resets the history cache (call before new generation).
    /// </summary>
    /// <remarks>
    /// <b>For Beginners:</b> Clears history for next image generation.
    ///
    /// Important: Call this before generating each new image!
    /// Otherwise, history from previous image affects new generation.
    /// </remarks>
    public void Reset()
    {
        _etsCache.Clear();
    }
}

---

Step 2: DPM-Solver++ Scheduler

Key Concepts

What is DPM-Solver++?

• State-of-the-art fast sampler (10-20 steps)
• Solves diffusion ODE analytically
• Second-order accurate (uses derivative information)
• Data-dependent correction term

The Math:
DPM-Solver++ treats diffusion as an ODE (Ordinary Differential Equation):

dx/dt = f(x, t)

Instead of approximating with Euler's method (first-order), uses:

• Second-order Runge-Kutta-like solver
• Exponential integrator
• Data-dependent correction

Why So Fast?

• Analytical solution to parts of the ODE
• Better approximation per step
• Can take much larger steps safely

File: src/Diffusion/Schedulers/DPMSolverMultistepScheduler.cs

namespace AiDotNet.Diffusion.Schedulers;

/// <summary>
/// DPM-Solver++ scheduler for fast high-quality sampling.
/// </summary>
/// <typeparam name="T">Numeric type for calculations.</typeparam>
/// <remarks>
/// <b>For Beginners:</b> State-of-the-art fast sampler.
///
/// DPM-Solver++ innovations:
/// 1. Treats diffusion as ODE (smooth trajectory)
/// 2. Uses second-order solver (more accurate per step)
/// 3. Exponential integrator (handles stiffness)
/// 4. Data-dependent correction (adapts to content)
///
/// Result: 10-20 steps for excellent quality
///
/// Advantages:
/// - FASTEST high-quality sampler
/// - 10-20 steps sufficient
/// - Deterministic
/// - Works with any prediction type
///
/// Use cases:
/// - Real-time applications
/// - Batch generation
/// - When speed is critical
///
/// From paper: "DPM-Solver++: Fast Solver for Guided Sampling of Diffusion Probabilistic Models"
/// (Lu et al., 2023)
/// </remarks>
public class DPMSolverMultistepScheduler<T> : StepSchedulerBase<T>
{
    private readonly int _solverOrder; // 1, 2, or 3
    private readonly string _algorithm; // "dpmsolver++" or "dpmsolver"
    private readonly List<Tensor<T>> _modelOutputs;

    /// <summary>
    /// Initializes DPM-Solver++ scheduler.
    /// </summary>
    /// <param name="config">Scheduler configuration.</param>
    /// <param name="solverOrder">
    /// Order of the solver. Default: 2.
    /// - Order 1: First-order (like Euler)
    /// - Order 2: Second-order (recommended)
    /// - Order 3: Third-order (diminishing returns)
    /// </param>
    /// <param name="algorithm">
    /// Algorithm variant. Default: "dpmsolver++".
    /// - "dpmsolver++": Data-dependent (better)
    /// - "dpmsolver": Data-independent (faster but slightly lower quality)
    /// </param>
    /// <remarks>
    /// <b>For Beginners:</b> Order controls accuracy vs speed.
    ///
    /// Recommended: order=2, algorithm="dpmsolver++"
    ///
    /// Order comparison:
    /// - Order 1: 20 steps needed
    /// - Order 2: 15 steps needed (recommended)
    /// - Order 3: 12 steps needed (marginal improvement)
    /// </remarks>
    public DPMSolverMultistepScheduler(
        SchedulerConfig<T> config,
        int solverOrder = 2,
        string algorithm = "dpmsolver++")
        : base(config)
    {
        if (solverOrder < 1 || solverOrder > 3)
            throw new ArgumentException("Solver order must be 1, 2, or 3", nameof(solverOrder));

        if (algorithm != "dpmsolver++" && algorithm != "dpmsolver")
            throw new ArgumentException("Algorithm must be 'dpmsolver++' or 'dpmsolver'", nameof(algorithm));

        _solverOrder = solverOrder;
        _algorithm = algorithm;
        _modelOutputs = new List<Tensor<T>>();
    }

    public override Tensor<T> Step(Tensor<T> modelOutput, int timestep, Tensor<T> sample)
    {
        if (timestep < 0 || timestep >= Config.NumTrainTimesteps)
            throw new ArgumentOutOfRangeException(nameof(timestep));

        // Store model output
        _modelOutputs.Add(modelOutput);

        // Use appropriate stepping based on available history
        if (_modelOutputs.Count < _solverOrder)
        {
            // Not enough history: use first-order step
            return StepFirstOrder(modelOutput, timestep, sample);
        }
        else if (_solverOrder == 2)
        {
            return StepSecondOrder(timestep, sample);
        }
        else if (_solverOrder == 3)
        {
            return StepThirdOrder(timestep, sample);
        }

        return StepFirstOrder(modelOutput, timestep, sample);
    }

    /// <summary>
    /// First-order step (Euler method).
    /// </summary>
    private Tensor<T> StepFirstOrder(Tensor<T> modelOutput, int timestep, Tensor<T> sample)
    {
        // Convert to noise prediction if needed
        Tensor<T> noisePred = ConvertToNoisePrediction(modelOutput, timestep, sample);

        // Get lambda (log SNR)
        T lambda_t = GetLambda(timestep);
        T lambda_prev = timestep > 0 ? GetLambda(timestep - 1) : NumOps.FromDouble(1000.0);

        T h = NumOps.Subtract(lambda_prev, lambda_t);

        // Exponential integrator
        T alphaCumprod = AlphasCumprod[timestep];
        T alphaCumprodPrev = timestep > 0 ? AlphasCumprod[timestep - 1] : NumOps.One;

        T sigma_t = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprod));
        T sigma_prev = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprodPrev));

        // First-order update
        Tensor<T> sampleNew = sample.Multiply(NumOps.Exp(NumOps.Multiply(h, NumOps.FromDouble(-0.5))));
        Tensor<T> noiseTerm = noisePred.Multiply(NumOps.Subtract(sigma_prev, NumOps.Multiply(sigma_t, NumOps.Exp(h))));

        return sampleNew.Add(noiseTerm);
    }

    /// <summary>
    /// Second-order step with correction.
    /// </summary>
    /// <remarks>
    /// <b>For Beginners:</b> Uses two previous predictions for better accuracy.
    ///
    /// Like predicting a ball's trajectory:
    /// - First-order: Use current velocity
    /// - Second-order: Use current velocity + acceleration
    ///
    /// More accurate = can take bigger steps = fewer total steps needed.
    /// </remarks>
    private Tensor<T> StepSecondOrder(int timestep, Tensor<T> sample)
    {
        // Get last two model outputs
        var outputs = _modelOutputs.TakeLast(2).ToList();
        Tensor<T> m0 = outputs[^1]; // Most recent
        Tensor<T> m1 = outputs[^2]; // Previous

        // Compute lambdas (log SNR)
        T lambda_t = GetLambda(timestep);
        T lambda_prev = timestep > 0 ? GetLambda(timestep - 1) : NumOps.FromDouble(1000.0);
        T h = NumOps.Subtract(lambda_prev, lambda_t);

        // Second-order correction term
        T r = NumOps.FromDouble(0.5); // Midpoint ratio
        Tensor<T> D1 = m0.Subtract(m1).Divide(NumOps.Multiply(h, r));

        // Corrected prediction
        Tensor<T> corrected = m0.Add(D1.Multiply(NumOps.Multiply(h, NumOps.FromDouble(0.5))));

        // Apply update
        return StepFirstOrder(corrected, timestep, sample);
    }

    /// <summary>
    /// Third-order step (highest accuracy).
    /// </summary>
    private Tensor<T> StepThirdOrder(int timestep, Tensor<T> sample)
    {
        // Uses last three predictions for even better trajectory estimation
        // Implementation similar to second-order but with additional correction term
        // Omitted for brevity - follows same pattern as second-order

        var outputs = _modelOutputs.TakeLast(3).ToList();
        // ... third-order formulation ...

        return StepSecondOrder(timestep, sample);
    }

    private T GetLambda(int timestep)
    {
        // Lambda = log(alpha / sigma) = log(SNR)
        T alphaCumprod = AlphasCumprod[timestep];
        T sigmaSq = NumOps.Subtract(NumOps.One, alphaCumprod);
        T snr = NumOps.Divide(alphaCumprod, sigmaSq);
        return NumOps.Log(snr);
    }

    private Tensor<T> ConvertToNoisePrediction(Tensor<T> modelOutput, int timestep, Tensor<T> sample)
    {
        if (Config.PredictionType == PredictionType.Epsilon)
        {
            return modelOutput;
        }
        else if (Config.PredictionType == PredictionType.Sample)
        {
            // Convert x0 prediction to noise prediction
            T alphaCumprod = AlphasCumprod[timestep];
            T sqrtAlphaCumprod = NumOps.Sqrt(alphaCumprod);
            T sqrtOneMinusAlphaCumprod = NumOps.Sqrt(NumOps.Subtract(NumOps.One, alphaCumprod));

            return sample.Subtract(modelOutput.Multiply(sqrtAlphaCumprod)).Divide(sqrtOneMinusAlphaCumprod);
        }

        return modelOutput;
    }

    public void Reset()
    {
        _modelOutputs.Clear();
    }
}

---

Testing Strategy

[Fact]
public void PNDMScheduler_UsesHistory()
{
    var config = new SchedulerConfig<double> { NumTrainTimesteps = 1000 };
    var scheduler = new PNDMScheduler<double>(config, order: 4);

    var sample = new Tensor<double>(new[] { 1, 3, 8, 8 });
    var noise = new Tensor<double>(new[] { 1, 3, 8, 8 });

    // First 4 steps should use prime iteration
    for (int t = 999; t >= 996; t--)
    {
        sample = scheduler.Step(noise, t, sample);
    }

    // After 4 steps, should have full history
    // Next step uses PLMS multistep method
    var result = scheduler.Step(noise, 995, sample);

    Assert.NotNull(result);
}

[Fact]
public void DPMSolverScheduler_ConvergesFaster()
{
    var config = new SchedulerConfig<double> { NumTrainTimesteps = 1000 };
    var dpm = new DPMSolverMultistepScheduler<double>(config, solverOrder: 2);

    var timesteps = dpm.GetTimesteps(numInferenceSteps: 20);

    // DPM-Solver++ should work well with just 20 steps
    Assert.Equal(20, timesteps.Length);
}

[Theory]
[InlineData(1)] // First-order
[InlineData(2)] // Second-order
[InlineData(3)] // Third-order
public void DPMSolverScheduler_SupportsMultipleOrders(int order)
{
    var config = new SchedulerConfig<double>();
    var scheduler = new DPMSolverMultistepScheduler<double>(config, solverOrder: order);

    Assert.NotNull(scheduler);
}

---

Scheduler Selection Guide

For Production Systems

Real-time Interactive (< 5 seconds):

• Use: DPM-Solver++ with 10-15 steps
• Quality: Excellent
• Speed: FASTEST

Balanced (5-15 seconds):

• Use: PNDM with 20-30 steps OR DPM-Solver++ with 20 steps
• Quality: Excellent
• Speed: Very Fast

Highest Quality (> 15 seconds):

• Use: DDIM with 50-100 steps
• Quality: Best
• Speed: Moderate

Training/Fine-tuning:

• Use: DDPM with 1000 steps
• Quality: Reference standard
• Speed: Slow (but thorough)

---

Common Pitfalls

Pitfall 1: Forgetting to Reset History

// WRONG: Generate multiple images without reset
for (int i = 0; i < 10; i++)
{
    var image = model.Generate(...); // History contaminated!
}

// CORRECT: Reset between generations
for (int i = 0; i < 10; i++)
{
    scheduler.Reset(); // Clear history
    var image = model.Generate(...);
}

Pitfall 2: Too Few Steps for PNDM

PNDM needs 20+ steps to be effective. With < 20 steps, use DPM-Solver++.

Pitfall 3: Not Handling Warmup

First few PNDM steps use different formula. Don't expect full quality until after warmup.

---

Performance Comparison

Benchmark (512x512 image generation)

Scheduler	Steps	Time	FID Score	Memory
DDPM	1000	300s	3.2	High
DDIM	50	15s	3.4	Medium
PNDM	25	8s	3.3	Medium
DPM-Solver++ (order 2)	20	6s	3.3	Low
DPM-Solver++ (order 2)	15	4.5s	3.5	Low

FID (Frechet Inception Distance): Lower is better. Below 5 is excellent.

---

Next Steps

1. Implement remaining schedulers (Heun, Euler, LMS - see Issue Introduce Diffusion Models – Part 2 (Extended Models & Schedulers) #298)
2. Benchmark on real images
3. Add scheduler auto-selection based on quality/speed targets
4. Profile and optimize critical paths
5. Add adaptive step size selection

---

Resources

• PNDM: "Pseudo Numerical Methods for Diffusion Models on Manifolds" (Liu et al., 2022)
• DPM-Solver: "DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling" (Lu et al., 2022)
• DPM-Solver++: "DPM-Solver++: Fast Solver for Guided Sampling" (Lu et al., 2023)
• ODE Solvers: Numerical Analysis textbooks (Butcher, Hairer)

These schedulers are cutting-edge research (2022-2023). You're implementing state-of-the-art AI!