Third-order DPM++ multistep coefficients fail a quadratic forcing test #15056

Describe the bug

I found a coefficient problem in the third-order DPM++ multistep update. With a scalar data predictor m(lambda) = lambda**2, the deterministic update returns half the exact integral, and the SDE mean update has the wrong sign.

The history is at lambda -2h, -h, 0, with target h. It gives D0=D1=0 and D2=h*h. The exact forcing is alpha_t * k * integral(exp(-k*(h-u))*u**2, u=0..h), with k=1 for DPM++ and k=2 for SDE-DPM++, from equations (8) and (14) of the DPM++ paper. Zero noise isolates the SDE drift term; the predictor here is independent of the state.

The multistep D2 is the quadratic coefficient, but the deterministic multiplier is half the required value. The SDE multiplier also subtracts 0.5 where the integral requires adding it. A local correction matches constant, linear, and quadratic controls, including unequal history spacing. On a cubic forcing problem, it restores observed order 3; the current updates approach orders 2 and 1 respectively. These are scalar consistency tests, not image-quality measurements or a general SDE convergence claim.

I checked the earlier third-order work in #9104. The expressions above remain in the current multistep scheduler. The deterministic expression also appears in the original DPM-Solver implementation. The single-step Heun SDE path has the same sign expression on inspection, though its D2 normalization differs.

I can attach the derivation and convergence test. Would it make sense to handle the related third-order coefficient cases together?

Reproduction

This does not need a model download:

import math
import torch
from diffusers import DPMSolverMultistepScheduler

h = 0.2
for algorithm, k in [("dpmsolver++", 1), ("sde-dpmsolver++", 2)]:
    scheduler = DPMSolverMultistepScheduler(
        algorithm_type=algorithm, solver_order=3
    )
    scheduler.sigmas = torch.tensor(
        [math.exp(2*h), math.exp(h), 1., math.exp(-h)],
        dtype=torch.float64,
    )
    scheduler._step_index = 2
    outputs = [torch.tensor([v], dtype=torch.float64)
               for v in [4*h*h, h*h, 0.]]
    actual = scheduler.multistep_dpm_solver_third_order_update(
        outputs,
        sample=torch.zeros(1, dtype=torch.float64),
        noise=torch.zeros(1, dtype=torch.float64),
    ).item()
    alpha = 1 / math.sqrt(1 + math.exp(-2*h))
    exact = alpha * (h*h - 2*h/k - 2*math.expm1(-k*h)/(k*k))
    print(algorithm, actual, exact)

Logs

System Info

CPU, float64; Python 3.12.13, PyTorch 2.13.0+cu130, Diffusers 0.35.2.

Reproduced with the installed scheduler and with the scheduler source from commit efe46a5b177af67c64b3666da8ddddac7ee12163, using the installed Diffusers runtime dependencies. The scheduler file is unchanged at 7ecb1e1b4fe500be7fb7381a7bb79a24d7fa3941 on October 15, 2026.