STORK: Faster Diffusion and Flow Matching Sampling by Resolving both Stiffness and Structure-Dependence
Authors: Zheng Tan
Weizhen Wang
Andrea L. Bertozzi
Ernest K. Ryu
Venue: ICLR 2026 Poster (2026)
Keywords: diffusion model
fast sampling method
stabilized Runge--Kutta
training-free
Primary area: generative models
TL;DR: We propose a new training-free fast sampler for diffusion models by resolving both stiffness and structure-dependence.

Abstract: Diffusion models (DMs) and flow-matching models have demonstrated remarkable performance in image and video generation. However, such models require a significant number of function evaluations (NFEs) during sampling, leading to costly inference. Consequently, quality-preserving fast sampling methods that require fewer NFEs have been an active area of research. However, prior training-free sampling methods fail to simultaneously address two key challenges: the stiffness of the ODE (i.e., the non-straightness of the velocity field) and dependence on the semi-linear structure of the DM ODE (which limits their direct applicability to flow-matching models). In this work, we introduce the Stabilized Taylor Orthogonal Runge–Kutta (STORK) method, addressing both design concerns. We demonstrate that STORK consistently improves the quality of diffusion and flow-matching sampling for image and video generation.

Published as a conference paper at ICLR 2026
STORK: FASTER DIFFUSION AND FLOW MATCH- ING SAMPLING BY RESOLVING BOTH STIFFNESS AND BRUUCTRRE-DEPNDENNE
Zheng Tan' Weizhen Wang* Andrea L. Bertozzi' Ernest K. Ryu'
I Deparrmment of Mathematics, University of California, Los Angeles
DDpaarment of Computer Science, University of California, Los Angeles 1zhengtan bertozzi, eryu @math. aclla.ddi rzwang1210@cs udua.edu
ABSTRACT
Diffusion models (DMs) and flow-matching models have demonstrated remarkable performance in image and video generation. However, such models require a significant number of function evaluations (NFEs) during sampling, leading to costly inference. Consequently, quality-preserving fast sampling methods that require fewer NFEs have been an active area research. However, prior training- free sampling methods fail to simultaneously address two key challenges: the stiffness of the ODE (i.e., the non-straightness of the velocity field) and depen- dence the semi-linear structure of the DM ODE (which limits their direct applicability to flow-matching models). In this work, we introduce the Stabilized Taylor Orthogonal Runge-Kutte (STORK) method, addressing both design con- cerns. We demonstrate that STORK consistently improves the quality of diffusion and flow-matching sampling for image and video generation. Code is available at htttps://gtthub com/ TTTT22222001S000
INTRODUCTION
Diffusion models (DMs) (Sohl-Dickstein et al. 2015 Ho et al.
2020 Song et al. 2021 have
demonstrated remarkable achievements during the past several years on various tasks, including image generation (Dhariwal & Nichol 2021 Meng et al. 2022). text-to-image generation (Ramesh
ct al. 2022, Rombach et al. [022 Gu et al. 2022 Podell et al. 2023 Mo et al. 2023 Zhang et al.
2023). video generation (Ho et al. 2022 Kong et al. 2024] Jin al. 2025 Liu et al. (025), and diffusion policy (Chi et al. 2024). As shown in Song et al. (2021), the sampling of DMs can be treated solving SDE number of function evaluations (NFEs) of the learned model when performing sampling, leading to costly inference.
equivalent ODE backward in time. However, DMs require a significant
Consequently, quality-preserving fast sampling methods have become an active area of research, and many methods that achieve effective sampling with fewer NFES have been proposed (Song et al. 2022. Liu et al. [022]] Zhhng & Chen 2023 Lu et al. 2022: 2025 Xie et al. 2025 Zhao et al. 2023). As noted in (Liu et al.) 2024a), the "straightness" of the velocity lield is key consideration in such fast sampling methods, and this observation corresponds to the notion of sriffness in classical numerical analysis (Burden & Faires [0TT). Building on this connection, exponential integrator techniques for stiff ODEs (Hochbruck & Ostermann) 2010) have been successfully adapted into the DPM-Solver (Lu et al. (022) and DEIS (Zhang & Chen, (023) samplers. However, these approaches rely on the semi-linear structure of the ODE formulation in noise-based diffusion models, and thus cannot be directly applied to flw-mmtcchin-bee
models. such cases, additional approximations,
such the data prediction step, are required (Lu et al. 2025 Xie al. 2025
Contribution. In this work, we introduce fast training- frren sampler, the Stabilized Taylor Orthogonal Runge-Kunnn Kutta (STORK) method, built upon stabilized Runge- Kutta (SRK) methods (Ab dulle 2002 Meyer al. 2014 Skaras al. 2021). As emphasized in Figure 4 STORK is

---

Published as a conference paper at ICLR 2026
stiff solver that does not rely the semi-linear structure of noise- based diffusion model's ODE formulation, making it applicable to both boise-based and flow-base how-based models. Our experiments show that STORK consistently improves the quality of diffusion and flow-matche mawching sampling for image and video generation. Flow-Fuler
Hlo--PPM-Solvertt
Flow-UniPC
STORK (Ours)
stork
strakr
stok
Stork
"A STORK holding plaque with the text "STORK" it. Cyberpunk, 8K resolution, highly detailed, ReX-On, super resolution" have the correct word with better details.
"A young black cowboy country music singer performing with a guitar on his horse in field outside the city to a crowd of fans, 8k, full HD" We have more vivid background and better realism.
"Joker logo with the initials the lstter Q" Only our method both the Joker logo and the letter Q.
NIEGE
"Create an image that evokes nostalgia, us using imagery of video games, movies, and other objects that remind us of our past" Our image is further detailed and has diverse objects matching the prompt.
Xie et al. 2025), Flow-UniPC (Zhao al. (202),, and STORK. All images are generated using the SANA 1.6B model (Xie et al. 2025) at 1024 x 1024 resolution with only 8 NFEs. Prompts are displayed beneath each image pair, accompanied by our commentary explaining why STORK's generations are superior. STORK achieves much better visual fidelity at the extremely low NFE case, showing its effectiveness fast sampling method. Zoom in for better visual details.
Figure 1: Comparison between the Flow-Fuler, Nlo--DPM-Solver+++ (Lu et al. 2025
2

---

Published as a conference paper at ICLR 2026
(a) Flow-UniPC video gemeration, 8 NFEs
(b) STORK video generation, NFEs
Figure 2: Video generation on Hunyuan model (Kong et al. 2024 with prompt:
"Iron Man is walking towards the camera in the rain at night, with lot of fog behind him. Science fiction movie, close-up", Our video portrays Iron Man more clearly and has rain in the background.
BACKGROUND AND RELATED WORKS
Diffusion models.
Diffusion models (DMs) generate samples by numerically solving the reverse-
time ODE
dx di f()xx(t) + g(t)²
20. (t))
[0,7],
(I)
where the "initial" condition æт sampled as Gaussian, f(t), g()), ar are known functions, and (s(x(t), t) is the noise prediction model, a trained neural network. For further information how diffusion models are trained and how they are able to generate high-quality images, we refer the
readers to standard references (Sohl-Dickster) et al. 2015 Ho al. 2020; Song et al. 2022 2021
Notice that equation I exhibits a semi-linear structure with linear term f(t)x(tt) and non-linear term involving eg- The DPM-Solver (Lu et al.) 2022 numerical ODE solver that exploits the semi-linear structure of equation] [ and it is one of the most widely used fast sampling methods for diffusion models.
Flow matching objectives. Flow matching is a class of generative models closely related to diffusion models and has recently attracted significant attention. Unlike diffusion models, which learn a noise model (score function), flow matching models directly learn a vector field vy and generate samples by numerically solving the forward-time ODE
dx
dt v()(t) t [0,7],
(2)
where the initial condition 2e is sampled as a Gaussian. For further information on how flow matching models trained and how they able to generate high-quality images, refer the readers to standard references (Lipman et al. 2023 Podell et al. 2023
Xie et al. 2025).
Because the flow matching ODE equation lacks semi-linear structure, the DPM-Solver (Lu et al. 2022) cannot be directly applied. Although follow-up works (L.u et al. 2025 Xie et al. [025) Zheng al. 2023) cireumvent this limitation using so-called data prediction step, this approach introduces errоrs at each step.
Classifier-free guidance for conditional sampling. In conditional sampling, generative model generates samples relevant to a conditional variable and classifier-free guidance is the most widely used technique for conditional sampling with diffusion models (Dhariwal & Nichol 2021]] Ho Salimans, 2022). Specifically, given parametrized model c(2,, e). the conditional noise is chosen be
E(#;, f, c) = sc(#,, f,
s)c( f)
where 0 is the classifier-f free guidance-scale, <(x,, t, c) is the noise model conditioned c, and c(x,, t, 0) represents the unconditional noise with @ serving as a place holder. Then, one simply replaces the noise e(x,,() !) in equation [ with e e( 1, 0).

---

Published as a conference paper at ICLR 2026
Main phase
Euler step
t-0
T-4h
T-3h
T-2h
T-h
NFE
Virtual NFE
Derivative dependency
x(0)
x(T)
Denoise via Flow ODE (noise to data)
Figure 3: Illustration of NFE evaluations for STORK 4 with & = 4 and 1st order Taylor approximation. For presentation clarity. we uniform timesteps with size h. "NFE" denotes actual NFEs, while the "virtual NFE" denotes NFEs approximated with the Taylor expansion. The arrows indicate that the previously computed velocity is used for first derivative approximations. Euler's method is used for the first step since there is no previous velocity.
For flow matching models such as Stable- Saabbl--Difuuiio--33 (Esser et al. 2024). FLUX. 1-dev Forest 2024), and SANA (Xie et al. 2025), similar approach is used, replacing the velocity model vr (( (1), of equation with similarly defined (tt
Current fast sampling methods. Existing fast sampling methods generally fall into two categories:
needs additional training, such knowledge distillation (Salimans & Ho 2022: Starodubcev
al. 2025), ) consistency models (Song et al.] 22033]
M023,La Lu & Song 2025 Chen et al. 2025). and one-step
diffusion models (Liu al., [2024a] Yin et al. [024] Chen al.
2025); the other training-free
which uses various numerical solvers for Stochastic Differential Equations (SDEs) and Ordinary
Differential Equations (ODEs) (Song et al., 20022 Karras al. 2022 Jolicoeur-Martineau et al. 2021
Liu et al. 2022 Zhang & Chen, 2023 Lu et al., Zheng et al. 2023 Zhao et al. 2023 2024). It
was first shown Song et al. (2021) that the sampling in DMs equivalent to a backward SDE or
ODE. Gotta Go Fast (Jolicour--Mariii Martineau et al. 2021) tried to accelerate the backward SDE sampling.
EDM (Karras et al. 2022) then tried to solve the issue that the trajectory between image and noise distribution is not straight enough with better training of DMs, and used Heun's 's method for solving the ODE. PNDM (Liu et al. 2022) demonstrated that direct utilization of classical Runge-Kuttu Kutta (RK) methods cannot result in fast sampling, therefore proposed pseudo-numerical method the noise manifold, using 4-step Adams- Bashforth and 4-step RK as initial steps. DEIS (Zhang & Chen 2023)
and DPM-Solver (Lu et al. 2022) were concurrent works that tried to utilize the semi-linear structure in the ODE in noise-based DMs, and used the exponential integrator (Hochbruck & Ostermann 2010)
technique solve it. Following-up works of the DPM-Solver (Lu et al.) [025 |025| al. 2023 Xie
et al. 2025) used data prediction to adopt the DPM-Solver into flowbbase models, and introduced
multi-step version of DPM-Solver. UniPC and DC-Solver (Zhao al.) [2033 2024) proposed
predictor-correct formulation that can be plugged into any fast sampling method in order to boost the convergence order for sampling with extremely small NFEs.
3 STABILIZED TAYLOR ORTHOGONAL RUNGE-KUTTH (STORK)
We now present our main method Stabilized Taylor Orthogonal Runge- Kutta (STORK), which based on stabilized Runge-Kungg Kutta (SRK) methods (Verwer et al. 1990: Abdulle & Medovikov
2001 Abdulle 2002 Meyer al. 2014: Sulllvvn 2019: Skaras et al. 2021 Tan al. 2025).
SRK methods in classical numerical analysis use orthogonal polynomials such as the Chebyshev polynomial /Veweer al. 1990 Abdulle & Medovikov 2001 Abdulle 2002) mitigate the effect of
stiffness (Burden & Faires (2011) and to allow for larger timesteps. However, standard SRK methods require an excessive number of function evaluations (NFEs), we approximate them and obtain our main methods STORK-2 and STORK-4, described fully as Algorithms and in Appendix
3.1 STABILIZED RUNGE- KUTTA (SRK) METHODS
Stabilized Runge- Kutta (SRK) methods is class of novel explicit solvers in numerical analysis. To motivates the necessity of applying SRK to noise-based and matching- based diffusion models, two important aspects named stiffness and structure-dependencc need to be considered.
4

---

Published as a conference paper at ICLR 2026
Stiff ODEs and their solvers. Prior works (Wizadwongsa & Suwajanakom 2023 Lu et al. 2022
Karras et al., 2022) have shown that the trajectory between the target distribution and the noise distribution may be "not straight enough" and that this causes challenges for fast sampling and high- quality image generation. This, in fact, corresponds to the classical notion of stiffness in numerical analysis (Burden & Faires (01TT) Loosely speaking, ODE is stiff if the local change in the solution slope is excessively fast so that large timesteps produce inaccurate or even unstable solutions that blow up, causing failure of typical explicit numerical methods. One modern technique is the exponential integrator (Hochbruck & Ostermann 2010), which is the motivation of the DPM-Solver and similar methods (Zhang & Chen [023] Lu et al. 2022 2025 Xie et al. 2025). More detailed discussion of the stiffness concept can be found in Appendix B
Structure-dependenncy Although the methods derived from the exponential integrator technique resolves the stiffness issue, they largely depend on the semi-linear structure of the ODE of the noise- based diffusion model. More precisely, this class of methods are designed for solving stiff ODEs of the form de La NN (x(t)) where L is non-zero linear operator and N is the non-linear part, and therefore called structure-dd -eeetuddd-deeeeett solvers. Methods without requiring any special form for the ODE are, therefore, called structure -nddeeenddeet solvers. Note that the ODEs for flo--matching models do not have the semi-linear structure. Therefore, seucture-deppndee solvers (Lu et al.) [2022]
2023 Xie et al., (025) cannot be directly applied, and additional techniques such as data prediction need to be conducted.
SRK methods. We now introduce the SRK method in numerical analysis under the flow matching setting; similar derivation holds by replacing v by Eg for noise-based DMs. Con- sider the ODE
Stiff solvers
Structure independent
Euler
DPM-Solver
SRK
Heun
dx dt
virtual
JUL
v(x(t),t), e [0,
DPM-Solver++
Runge-Kutta
STORK
PNDM
where x E Rd. The SRK methods is a class of single-step, multi-stage methods that are primar- ily designed allow one to numerically solve stiff ODEs with larger timesteps. In practice, the most widely used SRK methods are the second- order and fourth-order SRK methods (Abdulle
DEIS
UniPC
Figure 4: Derived from the SRK method, STORK both stiff problem solver and a structure- independent solver.
Medovikov 2001 Abdulle 2002 Verwer
et al.
1990 Meyer et al. 2014 0 Sullivan
2019
Skaras et al. 2021).
Let h denote the timesteps used for solving the flow matching ODE step from ==fo) 2((o h); in general, à can be non-uniform across different timesteps. An example of the second-order SRK (SRK2) method is called the Kunge-Kutt--Gegenbuurr (RKG2) method proposed by [Skaras et al. 202 1) which takes the form of
Yo xt(o); Y Yo- hijn Yo, to):
Y, = 1 + v Y, (1 Hj - 1/j Yn
(3)
P (Y)-,,t-,,)) hYo to). 2....,
x (to h) = Y..
Here, the SRK2 coefficients satisfy
2j + bj M4 =
[35 = 11yW1 = j+1 bj
./5 = 115
j b b
by-2
where
6
4(j - 1)(j +
(j (j+1)(1 + 2)
W1
bj
3j(j + 1)(1 + 2)( + +3) 3)
0==0 =1
by-
(s + ()(s 1)
2
5

---

SRK vs. RK methods. Classical Runge- Kutta (RK) methods and SRK methods crucially differ in their ability handle stiffness: RK methods are not stiff equation solvers, while the SRK methods are designed primarily for stiff equations. With more sub-steps 8, RK methods converge with higher arbitrarily Published as a conference paper at ICLR 2026
Algorithm 1 STORK-4 (Fourth order STORK)
Require: initial value ær, timesteps {4 velocity network o(; )), intermediate step number s, Taylor order 1,2,3. h; f1-1 fgp
a(t,-1) = xxt,) h (
for = to M n do
x=t;-1) xx4,) + (z5,,v(x(tt))))) 2,) 1( end for
for M
to to do Yo - æ(1;), Y1 ar4)) + bppp(((())
for j to s do
U Uapprox (Y; (1,11,t/ 1,)1) = TaylorExpansion(n. ,, (((),,
if Y, = hully Dapprox 1 - - j-
4 then
else
Y,
Y + VappyVappos (Y--1,
end if
end for
x x(t,-1) = Y
end for
return x(fo)
A fourth-order SRK method (SRK4) is the orthogonal Leage-Kutta-Chhbysse method developed
by Abdulle (002). The s-stage SRK4 method for solving this ODE from xtoo) to æ(to h) is
Yo (toa): Y = Yo h111 to),
Y; -1 5- 1- (-11) u, Kig Y j-2;
== 2,,, 2_.... ...,8 j == 3,
(4)
Y, = Y,-4-
a=to h) Y.
All parameters My, Vy: Ky are ODE-independent, pre-computed constants; the intermediate times E [to h, fo] are only dependent on tn and to h. The exact values of the parameters and more details of the derivation can be found in [[bduulle (2002) and Appendix(C] The SRK4 method allows one to choose h (ss²) maintaining stability (Abdulle) 2002), therefore handles the stiffness. call an update from Yj-1 1 to Y, as one sub-step, and an update from æ(4)) to æ(1o h) as super-step. Notice that in all the derivations above. never assume any special structures of v (x) (()))))) t). With the stiffness analysis in Appendix [ SRK methods uniquely belongs to both stiff' solvers and sructure-independent solvers, shown in Figure 4
order, while SRK methods handle stiffness better. Naturally, SRK methods large number with a general formula, while RK methods cannot. Detailed comparisons can be found in Appendix D
3.2 STORK: SRK WITH VIRTUAL NFE
Despite the success of SRK methods in classical numerical analysis, fatal issue when applying them on diffusion and fow-matching model sampling that an 8-stage SRK method requires & NFEs for updating super-step. In practice, s usually needs to be chosen to be around 10 to 50, leading to inordinate NFE count. As shown in Table naive application SRK4 to the CIFAR-10 (Krizhevsky et al. 2009) dataset results in very poor sampling results, especially in the range of small NFEs. Therefore, a mechanism to reduce the NEE count must be applied for the SRK methods to become practical.
natural choice to reduce NFEs in each super-step is to approximate v (t,(t,),t) using the Taylor expansion in the time variable vi Yo, fo). By treating the velocity purely (-dependent function,

---

Published as a conference paper at ICLR 2026
Table Ablation study comparison to the fourth-order Kutte-Kutta method and vanilla fourth-order stabilized Linge-Kuttal on CIFAR-10 1 et al. (2009) dataset. The NFEs in parentheses indicate the NFEs used for the RK4 method, since they must be multiples of for RK4.
Method NEF
10(12)
30(32)
50(52)
RK4
121.411
33,662
4.504
5.059
5.091
SRK4
443.812
40.828
6.225
6.324
6.167
STORK-4 (Ours)
5.497
4.167
3.888
3.809
3,789
Taylor expansion yields
= v PYo. to) + (t, - to)oooYo. (t, to)?
v (ty(tt)
v" to)
(t) - lo) v"' (Yo to) + O(tt, - t)
Since the exact evaluation of the higher-order derivatives also incurs costs no less than an NFE, we further approximate the derivatives of the velocity field vYY,, tp) using finite-difference approxi- mation. Except for the several initial super-steps (depending on the order of Taylor expansion used), we store the previously computed velocities and use the forward finiteddifference method to get an approximation for the derivatives of the noise or velocity up to the desired order. Regarding the initial steps, we use one Euler's method step followed by 2-step Adams-Bashforth method steps until there
enough points for Taylor expansion approximation. The intermediate Uuppeox (Y, ty) are called virtual NFEs, since they are approximated using the Taylor expansion so that no additional NFEs are used for those points opposed to the actual NFEs.
Plugging in the Taylor expansion and the velocity derivative approximations above into SRK4 de- scribed by equation suitable modification, we get the fourth-order Stabilized Taylor Orthogonal Runge-Kutta (STORK) method, denoted by STORK-4. The overall STORK method pipeline is illus- trated in Figure Similarly, we denote the first and second order STORK methods as STORK-1 I and STORK-2. A similar derivation works for noise-based DMs by using Taylor expansion on the noise. Algorithm summarizes the complete STORK-4 algorithm, with
TaylorExpanssonn(rder Yy(t)), Yo, to)
v vYY,, to) + (t, to Vapprox v(Yo, to) + (t, to Uapeeos (Y, to) + (t,- vi
(Yo to);
order 1.
2 apppoxx (Yo, to),
order 2.
o o(Yo, to) + (t, to ) Vapprox to) + (to=
Upppoox (Yo, (o)
+ (t, to ) "pprox (Y, to);
order 3.
More the details of the algorithms are in the Appendix E
We note that only first, second, and fourth order SRK methods exist for reasons related to the roots of the so-called stability polynomial (Abdulle, 2002). Therefore, third and higher-order SRK methods do not exist, and STORK-k can exist only for k = 1,2, 4. We find that STORK-4 consistently outperforms STORK-1 and STORK-2, so we conduct all the experiments using the STORK-4 method. The Taylor expansion order n is empirically chosen to be n with both unconditional and conditional moise-based generation, and n = 1 with conditional latent space flow-matching generation. Even though the choice of s does not affect the NFEs in the STORK method, the number of sub-steps s also requires tuning to maximize image quality. Intuitively, excessively large & causes errors of the Taylor expansion to accumulate; excessively small does not provide large enough region of absolute stability (Burden & Faires 201T). More studies on order and s are in Appendix G.T
Finally, via classical numerical analysis arguments in Abdulle (2002); Meyer et al. (2014); Skaras cet al. (202T) and the error term the Taylor expansion, the order of convergence can be summarized as shown in the following theorem and proved in Appendix [
Theorem 1. Assume e ca(0,, t) or v(æ,, t) satisfies the assumptions in Appendix. Let ( LM be
the sequence computed by STORK-k k with timesteps Mo For k = 2, 4, if Taylor expansion is not used for the virtual NFEs, the STORK-k solver converges to the expected solution with order k.
7

---

Published as a conference paper at ICLR 2026
20
20
30
DPM WWW-Solee44
--- DPM Salvere
Unift
10
DPM Solvere
STORK (Ours)
STORK (Ours)
11
UNIPE
F12
STORK (Ours)
9 10
12
15
10
NFE
(a) CIFAR-IO (32px, DDIM, Uncond)
(b) LSUN--eedoon
(c) MJHQ-30K (512px,
(256px, DDIM, Uncond)
Pixart-o, CFG=4.5)
Figure 5: Sample quality measured by FID for unconditional (Uncond) and classifier-fre--uiddee greeguudede (CFG) generation with noise- preddction models. As shown, STORK constantly outperforms other methods across datasets and image scales.
i.e. Eto To O(h²) where h max (t: 14-1). The Taylor expansion STORK-k solver in ISIS
Algorithm 2 and converges to the non Taylor STORK-K solver with order O
EXPERIMENTS
We compare the image and video generation quality of STORK against the most state-of-the-ard (SOTA) training-frree sampling methods, DPM-Solver++ (Lu et al. 2025) and UniPC (Zhao et al.) 2023). For image generation, we consider unconditional generation and conditional generation using both pixel and latent-space noise-predicting models (Song al., 2022, Chen al. 2023) and latent-space flow-matching models (Xie et al. 2025 Esser et al. 20244 BlackForest
To demonstrate the broad applicability of STORK, we benehmark on video generation using the Hunyuan Video (Kng et al., 2024) model. We also conduct an inference time and memory consumption analysis across different sampling methods to show that STORK-4 requires similar wall-clock running time and memory compared to other methods when using same amount of NFEs.
unless explicitly mentioned, all the noise-based generations use STORK- 4 in Algorithm | with second-order Taylor expansion, and flw-matcching-or based generations use STORK-4 with first-order Taylor expansion. Experiment details, additional metries, and supplemen- tary visualizations are provided in Appendices[G G and H
Following from Section 3
As shown in Figure 6 and Table 2 STORK consistently outperforms SOTA methods across tasks, model scales, and generation scales. We believe empirical evidence strongly favors STORK as better training-free sampling method.
4.1 IMAGE GENERATION: UNCONDITIONAL AND CONDITIONAL NOISE-PREDICTING MODELS
Experimental setup. We conduct both unconditional and conditional generations using both pixel- space and latent-space models. For unconditional generation, we experiment on CIFAR-10 (32 x 32) (Krizhevsky et al. 2009) and LSUNNBBdroom (256 256) (Yu al., 2016 using the pre-trained
DDIM (Ho et al. 2020) checkpoints provided by PNDM (L.iu et al. 2022). For each examined method, we generate 50K images to calculate Fréchet Inception Distance (FID) (Heusel et al. 2018) with respect to the Inception activation statistics provided by the PNDM codebase. For conditional generation, we use the Pixart-a on the MJHQ-30K (Li et al. 2024) dataset with 30,000 samples to calculate FID, with classifier-free guidance (CFG) scale 4.5. By empirical evidence, 14 is used on the CIFAR-10 10 dataset, and s = 24 is used for LSUN-Bedrowdd Bedroom and MJHQ-30K. To avoid the singularity at L O, denoise up a small value 0.
Results. As shown in Figure 5a bottom until all methods converge to a similar value. These results show the robustness of STORK across datasets and NFEs.
5b and 5c the FID curve for STORK consistently remains the
8

---

Published as a conference paper at ICLR 2026
-a- Flow-DPM-Solveet
11
HooNDDMM-Solve+++
10-
FIoN-UniPC.
STORK (Ours)
9
STORK (Ours)
8-
NFE
(a) MJHQ-30K (512px, SANA-0.6B, CFG=4.5)
(b) MJHQ-30K (1024px. SANA-1.6B, CFG=4.5)
20
Flow:
Flow-DPM DPM Solver+
Flow-DPM Www.DDMM-Sollccce Solver++
14-
Flow-UniPC
Flow-UniPC
STORK (Ours)
18
STORK (Ours)
17-
0
11
15-
10
12
14
12
MUHQ-30K (512px, FLUUX. 1-dev, CFG=3.5)
(d) MS-COCO (512px, SD-3.5-Large, CFG=3.5)
Figure 6: Sample quality measured by FID for conditional generation with latent space flow- matching models. As shown, STORK constantly outperforms other methods across various text-to- image flow-matching based generative models large margin.
4.2 IMAGE GENERATION: CONDITIONAL FLOW-MATCHIN MODELS
Experimental setup. We benchmark STORK on classifier-fre free text-to text-to-image generation using SANA (Xie al. 2025), FLUX. I-dev Forssk 2024), and Stable-Diffusion- Diffusion--35--Lagg (SD-3.5- L) (Esser et al.) 2024) as our latent-space flow-matching models. To demonstrate the scalability of STORK in terms of model size and generation scale, we use SANA-0.6B to benchmark at 512px resolution and SANA-L.6B to benchmark at 1024px resolution, using prompts from the MJHQ- 30K et al. (024) datasets. To further demonstrate STORK's synergy with more popular and 512px using FLUX. 1 1-dev and SD-3.5-L.
also benchmark generation Moreover, we prompt SD-3.5-L on the validation split of MS-COCO Lin et al. (2015) dataset to validate the robustness of STORK in terms of data distribution. All reported FIDs are calculated using 30k samples, with reference images resized to the corresponding generation resolution for Inception statistics calculation. Following (Esser al. 2024), we randomly sampled 30k validation image-promp pairs from MS-COCO as reference for FID calculations. Finally, the CFG scales for each benchmark are to the models' defaults, which are 3.5 for SD-3.5-L and FLUX.I-des and 4.5 for SANA.
larger-scale models,
Results. As shown in Figure ba to Figure 6d STORK demonstrates superior FIDs when using 7-15 NFEs. The non-trivial performance gap from STORK other sampling methods adds quantitative evidence for the superiority of STORK, besides the qualitative samples from Figure
4.3 VIDEO GENERATION: CONDITIONAL LATENT FLOW MATCHING
Experimental setup. We finally test our STORK method on the text-to-video generation task. To the best of our knowledge, we are the first training-frr free fast sampling work with experiments on video generation. We use the Hunyuan video diffusion model (Kong et al. (024).. with each frame generated at 512 x 320 resolution. For each video, 129 frames are generated at 15 frames per second (fps). Classifier-fra free guidance scale is set to Hunyuan's default value of 6. We benchmark on the EvalCrafter (Liu et al.) 2024b) evaluation suite, which consists of 700 prompts. Four different sub-metrict caleulated and aggregated final score.

---

Published as a conference paper at ICLR 2026
Table 2: EvalCrafter I evaluation of Hunyuan video diffusion model, with different sampling methods. Four sub-metrics and the final score are recorded. As shown in the table, STORK method methods in the final score.
consistently outperforms the Flow-DPM-Solven DPM-Solver++ and the
Method NFE
4
7
45.02
46.42
48.03
49.51
50.74
Visual Quality
Flow-UniPC
45.51
47,46
49.23
50.72
51.88
STORK (Ours)
50.00
51.91
52.11
52.72
52.68
Too--PPM-Solverttt
40.90
46.10
47.54
48.78
47.83
Text-Videe Alignment
Flow-UniPC
41.38
46.35
47.43
48.60
46.59
STORK (Ours)
43.18
46.11
46.64
50.09
46.92
Flow-DPM +
55.35
54.49
54.26
53.74
53.89
Motion Quality
Plow-nnPPC
55.24
54.56
54.06
54.12
53.37
STORK (Ours)
55.02
54.50
54.26
54.18
54.06
MDw-DPM--oolve+++
64.17
63.80
63.48
63.25
63.07
Temporal Consistency
Flow-UniPC
63.94
63.43
63.20
62.92
62.76
STORK (Ours)
61.41
61.67
62.08
62.14
62.26
205
211
213
215
216
Final Score
Flow-UniPC
206
212
214
216
215
STORK (Ours)
210
214
215
219
216
Table 3: Inference wall-clock profiling across different Now-matching sampling methods, all with 10 NFEs, over 10 trials.
Table 4: GPU memory usage during inference for different flow-matching sampling methods, with 10 NFEs.
Method
Avg (8)
Min (s)
Max (8)
Method
Max alloc. (GB) Max reserv. (GB)
Flow-DPM-Solveret
1.100
1.042
1,583
32.069
32.580
Flow-UniPC
1.223
1.167
1.699
Flow-UniPC
32.069
32.580
STORK-4 (Ours)
1.224
1.167
1.687
STORK-4 (Ours)
32.071
32.584
Results. Table 2 showcases the text-to-video generation results, As shown in the table, STORK constantly outperforms the -ow-DPM-Solvert+ (Lu et al.] 2025] Xie et al. (025) and the flow- UniPC (Zhao et al., 2023) methods in terms of the final score. For the sub-metrics, STORK constantly outperforms the other methods in terms of the visual quality, especially when NFEs are extremely small, and achieves comparable results in other metrics. This better video generation further support the general applicability of STORK to various generation tasks.
4.4 INFERENCE SPEED AND MEMORY CONSUMPTION ANALYSIS
We finally test the wall-clock time inference and the memory consumptions across different sampling methods. All experiments are conducted on a single H100-NVML GPU and generating images using batch size 1 with 11.UX.1-Des. loaded in torch. b oat generating at (512x512) resolution. For the time efficiency measurement, 10 NFES are used for each method with 10 trials. As shown in T the time usage similar for our proposed STORK-4 compared to other methods. For the memory usage, shown in Table 4 both the maximum allocated memory and the maximum reserved memory are similar across different methods.
CONCLUSION
We propose a novel training-free fast sampling method called STORK for unconditional and condi- tional diffusion and flow matching models. Our method is ODE-structure independent; therefore, it can be applied to noise and flow matching-bart based DMs without any further adjustments. Nevertheless, sampling using
smaller NFEs with STORK and its possible future extensions remains to be explored. Therefore, further investigation into the potentials of STORK for effective sampling of diffusion and flow matching models presents promising direction for future research.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
Zheng Tan and Andrea L. Bertozzi are supported by NSF grant DMS-2152717, ONR grant DOD N00014-23-1-2565, and DOE grant DE-SC0025589, Ernest K. Ryu is supported by NSF grant CCF-2504627.
The first author would like to thank Sicheng Mo for valuable discussions on diffusion models.
REFERENCES
A. Abdulle. Fourth order Chebyshev methods with recurrence relation. SIAM Journal on Scientific Computing, 23(6): 2002. doi: 10. 10.1137///10448770033555554 URL https: /doi org/10. 113./1 0tt.11117/5111422503733333333330
A. Abdulle and A. A. Medovikov. Second order Chebyshev methods based on orthogonal polynomiass Numerische Mathematik, 90(1):1-18, 2001. doi: 10. 00.10077/00221111002222
BlackForest. Black forest labs; frontier ai lab, 2024. URL /ttp/tttpss/bbbaccfreetttllaaa.. URL https://lackfooretlabsssa
R. Burden and D. Faires. Numerical Analysis. Brooks/Cole, Boston, 9th edition, 2011.
J. Chen, Yu, C. Ge, L. Yao, Xie, Y. Wu, Z. Wang. J. Kwok, Luo, H. Lu, and Z. Li. PixArt-a: Fast training of diffusion transformer for photorealistic text-to-image synthesis. 2023,
J. Chen, S. Xue, Y. Zhao, J. Yu, S/ Paul, J. Chen, H. Cai, E. Xie, and S. Han. SANA-Sprint: One-step diffusion with continuous-time consistency distillation. arXiv preprint arXiv:2503.0964/ 2025.
C. Chi, Xu, Feng, E. Cousineau, Y. Du, B. Burchfiel, R. Tedrake, and Song. Diffusion policy: Visuomotor policy learning via action diffusion. Robotics: Science and Systems, 2024.
P. Dhariwal and A. Nichol. Diffusion models beat GANs on image synthesis. Advances in Neural Information Processing Systems, 2021.
Esser, S. Kulal, A. Blattmann. R. Entezari, Müller, H. Saini, Levi. Lorenz, A. Sauer, F. Boesel, D. Podell, T. Dockhorn, Z. English, K. Lacey, A. Goodwin, Y. Marek, and R. Rombach. Scaling rectified flow transformers for high-resolution image synthesis. International Conference on Machine Learning. 2024.
S. Gu, D. Chen, J. Bao, F. Wen, B. Zhang, D. Chen, L. Yuan, and B. Guo. Vector quantized diffusion model for text-to-image synthesis. Computer Vision and Pattern Recognition, 2022.
M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter. GANs trained by two time-scale update rule converge to a local Nash equilibrium. Advances in Neural Information Processing Systems, 2018.
J. Ho and T. Salimans. Classifier-free diffusion guidance. Advances in Neural Information Processing Systems, 2022.
J. Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 2020.
Ho, T. Salimans, A. Gritsenko, W. Chan, M. Norouzi, and D. J. Fleet. Video diffusion models. International Conference on Learning Representations, 2022.
M. Hochbruck and A. Ostermann. Exponential integrators. Acta Numerica, 19:209-2888 2010. doi: 10.101775096249991000000
Y. Jin, Z. Sun, N. Li, K. Xu, K. Xu, H. Jiang, N. Zhuang, Q. Huang, Y. Song. Y. Mu. and Z. Lin. Pyramidal flow matching for efficient video generative modeling, 2025. URL. htttss://axxi org/abs/ 2410. 05954
A. Jolicoeur-Marte Martineau, K. Li. R Taillefe,,, T. Kachman, and I. Mitliagkas, Gotta go fast when generating data with score-based models, 2021. URL https:stt..........ttt //t/iv..vvvvv. orq/ org/abs /2005 14080
11

---

Published as a conference paper at ICLR 2026
T. Karras, Miika. Aittala, T. Aila, and S. Laine. Elucidating the design space of diffusion-bian -ased generative models. Advances in Neural Information Processing Systems, 2022.
W. Kong, Q. Tian, Z. Zhang, R. Min, Z. Dai, J. Zhou, J. Xiong, X. Li, B. Wu, J. Zhang, al. Hunyuanvideo: A systematic framework for large video generative models. arXiv preprins arXiv:2412.03603, 2024.
A. Krizhevsky, G. Hinton, et al. Learning multiple layers of features from tiny images. Master's thesis, University of Toronto, 2009.
D. Li, A. Kamko, E. Akhgari, A. Sabet, L. Xu, and S. Doshi. Playground v2.5: Three insights towards enhancing aesthetic quality in text-to-image generation, 2024.
T.-Y. Lin, M. Maire, S. Belongie, L. Bourdev, R. Girshick, J. Hays, P. Perona, D. Ramanan, L.C. Zitnick, and P. Dollár. Microsoft COCO: Common objects in context, 2015.
Y. Lipman, R. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le. Flow matching for generative modeling. International Conference Learning Representations, 2023.
D. Liu, S. Li, Y. Liu, Z. Li, K. Wang, X. Li, Q. Qin, Y. Liu, Y. Xin, Z. Li, et al. Lumina-videw Efficient and flexible video generation with multi-scale next-dit, arXiv preprint 0rXiv22502.06782, 2025.
L. Liu, Y. Ren, Z. Lin, and Z. Zhao. Pseudo numerical methods for diffusion models manifolds. International Conference on Learning Representations, 2022.
X. Liu, X. Zhang, J. Ma, Peng. and Q. Liu. InstaFlow: One step is enough for high-quality diffusion-based text-to-image generation. International Conference on Learning Representations, 2024a.
Y. Liu, X. Cun, X. Liu, X. Wang, Y. Zhang, H. Chen, Y. Liu, T. Zeng, R. Chan, and Y. Shan. Evalcrafter: Benchmarking and evaluating large video generation models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 22139-22149, 2024b.
C. Lu and Song. Simplifying, stabilizing and scaling continuous-time consistency models, 2025. URL
11081
C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, and Zhu. DPM-Solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps. Advances in Neural Information Processing Systems, 2022.
C. Lu, Y. Zhou, F. Bao, J. Chen, Li, and J. Zhu. DPM-Solvert+: Fast solver for guided sampling of diffusion probabilistic models. Machine Intelligence Research, 22:730-751, August 2025. doi: 10. 1007/s11633-000 133--33----2---- htLpttps:: / /oi. org wttp/////................. 1007 001.:00//1111130------1155111111
C. Meng, Y. He, Song, J. Song, J. Wu, Zhu, and Ermon. SDEdit: Guided image syn- thesis and editing with stochastic differential equations. International Conference on Learning Representations, 2022.
C. D. Meyer, D. S. Balsara, and T. D. Aslam. A stabilized Runge-Kutte Kutta- method for explicit super-time-stepping of parabolic and mixed equations. Journal of Computational Physics, 257: 594-626, 2014. ISSN 0021-9991. doi: ott.a//ccc........acc.aooaaaaaa.
021. URL https
www. www.:ccieeccee.c rect com/
S. Mo. F. Mu, K. H. Lin, Y. Liu, B. Guan, Y. Li, and B. Zhou. FreeControl: Shaining-free spatial control of any text-to-image diffusion model with any condition. Computer Vision and Pattern Recognition, 2023.
S. O' Sullivv.. Runge- Kutta Gegenbauer explicit methods for adveecioondiiffssiiii probless Journal of Computational Physics, 388:209 223, 2019, ISSN 0021-9991. doi: http://tbb.org/10.111//// 2019.03.001 URL https: www. www.:scciieeeecceeeeeeee https://w...ceeeeddreee.... com httpp//ewe.e.ceeeeceeeeeeeeaaaaa S0021 199199991111 9301706
12

---

Published as a conference paper at ICLR 2026
D. Podell, Z. English, K. Lacey, A. Blattmann, T. Dockhor,, J. Müller, J. Penna, and R. Rombach. SDXL: Improving latent diffusion models for high-resolution image synthesis. International Conference Learning Representations, 2023.
A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, and M. Chen. Hierarchical textcncnioal image generation with CLIP latents, 2022. URL
R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer. High-resolution image synthesis with latent diffusion models. Computer Vision and Pattern Recognition, 2022.
T. Salimans and J. Ho. Progressive distillation for fast sampling of diffusion models. International Conference on Learning Representations, 2022.
Skaras, T. Saxton, C. Meyer, and T. D. Aslam. Super-time-steppinn schemes for parabolic equations with boundary conditions. Journal of Computational Physics, 425: 2021. ISSN 0021-9991. doi:
109879. URL https: http:://www.cccc. wtwp://www.wei.eeeeeeeee.e.eeeee com/
J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. International Conference on Machine Learning, 2015.
J. Song, C. Meng, and S. Ermon. Denoising diffusion implicit models. International Conference on Learning Representations, 2022.
Y. Song, J. ative modeling through stochastic differential equations. International Conference on Learning Reprexentations, 2021.
D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based gener-
Y. Song, P. Dhariwal, M. Chen, and I. Sutskever. Consistency models. International Conference Machine Learning, 2023.
N. Starodubcev, D. Kuznedelev, A. Babenko, and D. Baranchuk. Scale-wise distillation of diffusion models, 2025. URL
G. Steinebach.
ROCK4.
www htt.:///..///////..occccoooooo..
Fileexchange/ (itpaaannne22299-roaaee rock4
2006. Version 6.2006.
Z. Tan, D. Aslam, and A. L. Bertozzi. Explicit monotone stable super-time-steppinn methods for finite time singularities. arXiv preprint arXiv:2507.17062, 2025.
J. G. Verwer, W. H. Hundsdorfer, and B. P. Sommeijer. Convergence properties of the Runge- Kutta- Chebyshev method. Numerische Mathematik, 57:157-178, 1990. doi: 10. 10.10077/F001386405 URL
https:///oi. org/ 101.:000110000011111100000111100 3F4448
P. Platen, S. Patil, A. Lozhkov, P. Cuenca, N. Lambert, K. Rasul, M. Davaadorj. D. Nair, S. Paul, W. Berman, Y. Xu, S. Liu, and T. Wolf. Diffusers: State-of-the-art diffusion models. https:://githubee com
diffusers) 2022.
S. Wizadwongsa and S. Suwajanakom. Accelerating guided diffusion sampling with splitting numerical methods. International Conference on Learning Representations, 2023.
X. Wu, Y. Hao, K. Sun, Chen, F. Zhu, R. Zhao, and H. Li. Human preference score v2: A solid benchmark for evaluating human preferences of text-to-image synthesis. arXiv preprint 2XXx:2300600941 2023.
E. Xie, J. Chen, J. Chen, H. Cai, H. Tang. Y. Lin, Z. Zhang, M. Li, L. Zhu, Y. Lu, and S. Han. SANA: Efficient high-resolution image synthesis with linear diffusion transformers. International Conference on Learning Representations, 2025.
T. Yin, M. Gharbi, R. Zhang, E. Shechtman, Durand, W. Freeman, and T. Park. One-step diffusion with distribution matching distillation. Computer Vision and Pattern Recognition, 2024.
F. Yu, A. Seff, Y. Zhang, S. Song, T. Funkhouser, and J. Xiao. LSUN: Construction of a large-scale image dataset using deep learning with humans in the loop, 2016. ARL/atpps: abs 2bs/1506 03365
13

---

Published as a conference paper at ICLR 2026
L. Zhang, A. Rao, and M. Agrawala. Adding conditional control to text-to-image diffusion models. International Conference on Computer Vision, 2023.
Q Zhang and Y. Chen. Fast sampling of diffusion models with exponential integrator. International Conference Learning Representations, 2023.
W. Zhao, L. Bai, Rao, Zhou, and J. Lu. UniPC: unified predictor-corret correccoo framework for fast sampling of diffusion models. Advances in Neural Information Processing Systems, 2023.
W. Zhao, H. Wang, J. Zhou, and J. Lu. DC-Solver: Improving predictor- diffusion sampler via dynamic compensation. European Conference on Computer Vision, 2024.
K. Zheng, C. Lu, Chen, and J. Zhu. DPM-Solver-v3: Improved diffusion ODE solver with empirical model statistics. Advances in Neural Information Processing Systems, 2023.
14

---

Published as a conference paper at ICLR 2026
APPENDIX
A LARGE LANGUAGE MODELS (LLMS) USAGE
We did not LLMs for the writing of this paper, and LLMs did not help to the extent that they could be regarded as contributor during the research process.
B BRIEF ILLUSTRATION OF STIFFNESS
this section, we briefly illustrate the classical notion of stiffness in numerical analysis. Consider the ODE,
dæ 20xr xx0) 1, E [0, 1].
dt
It can be easily shown that s=t) -20 However, when the ODE above is solved numerically using various methods, some methods exhibit spurious oscillations, as illustrated in Figure While the numerical solution converges the exact solution in the limit the step size h significant errors can arise for moderate values of h, particularly when the analytical solution trajectory is "not straight enough."
20r,
1.00
0.75
0.50
0.25
0.00
-0.25
Euler (RK1)
-0.50
Heun (RK2)
-0.75
RKG2 (STORK-2 without Taylor)
Analytical solution
-1.00
0.0
0.2
0.4
0.6
0.8
1.0
Figure 7: Example of a stiff ODE. The plot shows the analytical and numerical solutions of the ODE equation using Euler, Heun, and the second-order Runge-Kutee Kutta- Gegenbauee (RKG2) method with 4 sub-steps. Euler's method and Heun's method are two of the most popular flow matching sampling methods, implemented in diffuser (von Platen et al. 2022). The RKG2 method is our STORK-2 method with exact sub-steps. Ten timesteps are used for each of the numerical solutions.
As shown the plot, Euler and Heun's method has significant errors while RKG2 is close to exact solution.
The region of absolute stability characterizes the notion of stability for numerical methods. Consider the test problem
dx
Ar.
di
where À E C can be arbitrary complex number. The region of absolute stability of numerical method is defined as the set of hy values such that the numerical method with step size h applied to equation yields solution that remains bounded. Since the theoretical exact solution to equation
15

---

Published as a conference paper at ICLR 2026
Region of Absolute Stability for Different Methods
Euler
Heun
4
RK4
RKG2, 5=3
RKG2, s=5
2
RKG2, s=7
Im(2)
-2
-4
-6
-25
-0
-15
-10 Re(z)
5
Figure 8: Region of absolute stability for different numerical methods. Euler and Heun's 's methods are widely used in the diffusion model sampling, especially in the flow matching setting. The RK4 method is the most widely used explicit Runge-Kutta method in numerical analysis. The RKG2 method the basis for our STORK-2 method; this plot, and sub-steps are used. As shown in the plot, the RKG2 method has much larger region of absolute stability than the other methods, so that it is much more stable, therefore more suitable for solving stiff problems. Moreover, the region of absolute stability becomes significantly larger as increases.
remains bounded for all in the left half of the complex plane, a numerical method is considered more stable if its region of absolute stability encompasses larger portion of that half-plane.
As an example, when applying Euler's method to the test problem equation 6 the numerical method becomes
(tt h) xtt) + htx((t) (1 + ht)x(t).
Therefore, the stability polynomial is R(=) = + the region of absolute stability is therefore +|| < 1, which corresponds to a circle centered around with radius in the complex plane.
The region of absolute stability for Heun's method, fourth-order Kunge-Kutta (RK4) method, and the second-order
(RKG2) with sub-steps are plotted Figure g The RKG2 method, for which the STORK-2 method is derived based on, has region of absolute stability that contains a much larger portion of the left half complex plane than the Euler, Heun, and RK4 methods. Analytical derivation in Appendix C further confirms that s-stage RKG2 method has large region of absolute stability that grows with order O(s2).
C STABILIZED RUNGE-KUTTA METHODS
present the derivation of the second-order Kutae-KKtta-Geeebbauue (RKG2) method (Skaras)
et al. 202 and the fourth order orthogonal Runge-Kutter Kuttta-hhhbbbhee (ROCK4) method (Abdulle 2002).
The derivations closely follow the original derivations Skaras et al. (202 and Abdulle
(2002). Notice RKG2 corresponds to the SRK2 method and ROCK4 corresponds to the SRK4 method in the main content.
Consider the ODE
da Mx, x0) = Zo! M X
DE R" .nnnn N4.
dt
16

---

Published as a conference paper at ICLR 2026
The analytical solution for the above linear equation is defined by
(tM) ym m! Zo.
x(1) = (x())
Since would like derive single-step method, the method takes the form of x(t - h) = ((hM)x(()).
(8)
where the timestep size, and R(hM) is called the stability polynomial. order to enlarge the size of the region of absolute stability, the shifted Gegenbauer polynomial with parameter to be the stability polynomial; this leads to the Runge-Kurte BeageeKKtta-GGegenbe Gegenbauer (RKG) method (Skaras
$
2 is chosen
et al. 2021: O' Suullivan 2019). More precisely, the stability polynomial of an s-stage RKG method Rg (M) satisfies
Ra(2) a. -b, (1 wi=),
(3/2)
where 0gy bs w1 are parameters to chosen, and C////2 (x) the s-degree Gegenbauer polynomial with parameter 3 (Skaras et al. 2021). In order for the method to be convergent, we need to equate equation 7 and equation [ to the highest order possible. Since there are three unknowns, the best possible order is to equate the first three terms and get second-order (h²) convergence. Therefore, the three equations that are needed to be satisfied are
R,(0)= R2 (0) 1.
Equating and solving the differential equation, we get the result for RKG2 method:
4(j - 4)j + + 2)
0j = 1 (j + 1) + +2)
W1 4)s by - +
by.
1)
2
To ensure numerical stability in Appendix B
one only needs ensure that 1 + w.hM has eigenvalue within [- so that by the boundedness of the Gegenbauer polynomial we get numerical stability for free Skaaas et al. 2021). This requires that h hepplicit (e+1))e-1) O(s²). where hexpliiit is the maximum timesteps that ensures Euler's method to be numerically stable.
Finally, would like to turn the RKG2 method into the Runge-Kutta Kutta form as shown in equation 3 so that it can be implemented easily in practice. Using the well-known Gegenbauer polynomial inductive relationship
cla) (2) _ 2z(s + a 1)c (s 2a - 2)) 8
can be derived that the stability polynomial satisfies the relationship
aj b, c (3/2) (1+ = +
(3/2) (3/2)) (1 w1z))
(3/2) (1 + (1- - U2)
+ b
- + 7tr
where
2j + bj
+ bj
Mi j
ß = 1 Vj j b
93z -1z
by-1"
Therefore, the RKG2 method in the Lunge-Kutte. formulation is
Yo = xxt),
Y Yo ii
Y: = Hy M1 + 12 Y22 + M4 v1
2
ztt h) Y..
Replace each MY, in equation 9 by v (Y, t, exactly recovers equation of our STORK-2 2 method follows as Section
Then the derivation
The fourth-order orthogonal Runge-Kutta Cheitt-hv (ROCK4) method can be derived similarly, with an extra composition method term. The stability polynomial designed to be
(2)=)) wa Pa
17

---

Published as a conference paper at ICLR 2026
which satisfies [R)(2)] I 1 for [-1,, 0), with as large possible. The first polynomial
wa(2)= (1 ; 21 )) (1-
1
1
22
Z2
is a fourth order polynomial that serves composition method (Burden & Faires 2011 1) to boost the entire method to fourth order, and the second polynomial Pa-4(z) is an orthogonal polynomial
that is orthogonal with respect the weight function (a)² and normalized such that Po)4(0) 1. vT-27
Detailed derivation of the polynomial (2) is similar to the RKG2 method and can be found in Abdulle (2002). For the expression of (=) that serves
composition method, it evokes the concept of Butcher's 's Tableau (Burden & Faires (011) in classical numerical analysis, and the coefficients are derived by solving system equations with 10 unknown variables, that two degrees of freedom exist and the original ROCK4 paper Abdulle (2002) chose the optimal parameters by experiments. The coefficients material. Both our implementation and the coefficients in the STORK method are heavily based on the implementation in Steinebach (2006).
precomputed and provided along with the supplementary
D COMPARISON WITH SIMILAR METHODS
In this section, compare STORK with existing sampling methods and highlight their relationships. The following discussion motivates the effectiveness of STORK, and the ablation study in Table shows the necessity of modifications from stabilized Runge-Kutta STORK methods.
D.I COMPARISON BETWEEN STORK AND CLASSICAL RUNGE-KUTTA METHODS
When = 1, STORK-2 and STORK-4 -4 both reduce to Euler's ' Method, which is a 1 1-step Runge- Kutta (RK) method. However, for 1, STORK-2 and STORK-4 are not classical RK methods, and they exhibit some fundamental qualitative differences. The classical RK methods, such as Hcun's s method or RK4, designed for higher-order accuracy the number of sub-steps increases (Burden & Faires 2011), while the STORK methods (or SRK methods) are designed to address the stiffness in ODES and PDEs (Verwer al. 1990 Abbuulle 20022) with regard to the step size. Regardless of the choice of sub-step number the STORK-K method in Algorithm and [ with a fixed A: always has the same order as shown in Theorem
Another essential difference is that the number of sub-steps s in STORK is an easily tunable hyperparameter. For stiffer problems, one can straightforwardly increase With s-step RK methods, the parameter 8 is usually not considered a hyperparameter, since changing s completely changes the method. By tuning s and using many sub-steps between the super-step from i, to f;-1, STORK can effectively address stiffness, while direct applications of traditional RK methods have not lead to significant improvements in the sampling quality (Liu et al. 2022 Lu et al. 2022) due to the stepsize constraint that stiffness imposes on RK methods.
D.2 COMPARISON BETWEEN STORK AND VANILLA STABILIZED RUNGE KUTTA METHODS
Although the SRK method single-step method, i.e. one only needs the value of xxt) in order to get the value of æ=( h), the Taylor expansion makes STORK a multi-step method (Burden & Faires 201T). However, STORK is not linear multi-step method. linear multi-step method, such as the 2-step Clams-Bashffrth method
att h) x(t) 2/2 (((()))) f(x(t 1
reuses each previous function evaluation once per step. Other pseudo multi-step numerical methods, such as PNDM (Liu et al.) 20022), also have the same mechanism. STORK, on the other hand, repeatedly uses previous function evaluations in the Taylor expansion approximation. Intuitively speaking, more information about the local velocity fields in between t and h is extracted from those evaluations, so that the solution would follow the learned velocity field more closely, leading to better approximation. Experiments in Section demonstrate 4
the effectiveness of this mechanism.
18

---

Published as a conference paper at ICLR 2026
Algorithm 2 STORK-2 (Second order STORK)
Require: initial value ær, timesteps { velocity network o(; )), intermediate step number s, 1,2,3.
Taylor order
h; f1-1 44
a(t--1) = ztt,) + h
for = to M do
x=t;-1) =ff,) + 1.5h,v E.) -05h)--(((((-))
end for
for M
to 0
Yo = æ(1;), Y1 ar44) + h.p)
for j to s do
U Vapprox (Y,-1, 1 = TaylorExpansion(t. 1, Yo. ts).
Y, My y- + v Y y-
Vgy Yo + Kiy Vuppoux (Y, y Hy,,,
end for x ((,--1) Y.
end for
return zffo)
E STORK-2 AND STORK-4 ALGORITHMS
We now present the details of STORK-2 2 and STORK-4 algorithms described Section 3 For both methods, one Euler's step is used the beginning. If Taylor expansion order is then one more step of 2-step Adams-Bass Bashforth method used, and if Taylor expansion order is 2. two more step of 2-step Adams-Bashforth method is used. The later steps follow the derivation of SRK2 and SRK4 methods in Section 3 Details of the algorithms are shown in Algorithm and Algorithm for flow-based models. Similar algorithms can be written down using exactly the same procedure as in Section 3 on noise-based models, by Taylor expansion on the noise. We later noticed that we essentially swapped Yo and Y1 in our official STORK-4 implementation, so the same change needs be done if one wants to reproduce the results in the paper.
F CONVERGENCE PROOF FOR THEOREM
In this section, we present the proof of Theorem
Proof of Theorem 7
By derivation in C the RKG2 method converges with order O(h?) and the ROCK4 method converges with order Ohh) expected.
Now we'll show that the STORK-2 method and STORK-4 method converge to RKG2 and ROCK4, respectively, with order O(h2) for both cases. By Taylor expansion, know that
tt, v(Yo,to) to) + (t, - ta) (Yo to) + O(t, to)
Since the three point forward velocity approximation has order ((h) by classical numerical analysis, and 1, to & h in both STORK-2 and STORK-4, we know that
o (t,t)),tt))) vYo. to) + (ty to) to Upppoos (Yo, to) + (t, to) 2 approx (Yo, to) O(h).
Plugging into Algorithm 2 and Algorithm RKG2 equation 3 and ROCK4 equation respectively, with order O(h2) since both STORK-2 and STORK-4, there exists an additional h in front of each virtual NFE a (Y (t,tt),tz).
know that the STORK-2 and STORK-4 converge to
G EXPERIMENTS
In this section, discuss experiment setup in greater detail. We present all available data tables for clearer demonstration. For UniPC (Zhao et al. 2023). DPM-Solver++ (Lu et al. 2025), DEIS (Zhang & Chen 2023). Flow-Fuler, and DDIM (Song et al. 2022) schedulers,
19

---

Published as a conference paper at ICLR 2026
the implementation provided by the Diffusers (von Platen et al. 2022) package version 0.35. D devO.
Note that in order for the compatibility of FLUX.I-dev and Hunyuan Video (Kong et al., 2024) pipelines with DPM-Solver++ and UniPC, we made two lines of modification on the original Diffusers source code. Essentially, we enforce sigmas None if the scheduler's config has use_flow flow sigmas = True. We can provide the updated implementation upon request for reproducibility.
G.1 STUDIES ON THE PARAMETERS
In this subsection, we conduct ablation studies on the number of sub-steps s used, the order of the STORK method, and the order of the Taylor expansion.
G.1.1 EFFECT OF SUB-STEPS
We investigate the effect of sub-steps the generation fidelity of diffusion models. For these experiments, we use the SANA (Xie et al. 2025) 0.6B variant and generate images at 512 x 512 resolution for FID calculation using 30000 samples. As shown in Table while all choices of lead to decent results, the choice better stiffness handling ability, overly large results in too much Taylor series approximation errors. We found that 8 = 9 is optimal for the conditional flow-matching generation, and recommend users experiment with choices of < 100 for the downstream models and tasks.
cannot be excessively large or small. However, larger & enables
Table 5: Effects of sub-steps sample fidelity. As shown, sub-step leads to decent performance. and we empirically find achieves the best performance with SANA-0.6B by (Xie et al. 2025) generating at 512 x 512 resolution and in general for flow-matching models.
s NFE
8
9
10
s=5
8.000
7.999
8.537
9.168
8=9
7.659
6.667
6.526
6.270
s=14
8.351
7.242
7.026
6.671
8=24
8.686
7.526
7.286
6.884
s=54
8.836
7.649
7.384
6.992
s= 104
8.866
7.669
7.412
7.016
G.1.2 EFFECT OF SOLVER ORDER AND TAYLOR EXPANSTON ORDER
We further conducted ablation studies for both Taylor expansion order and solver order. As mentioned in Section 3 STORK-I, STORK-2, and STORK-4 are the only possible configurations for SRK methods. We now examine their corresponding performance using the SANA-0.6B model for 512x512 resolution on the MJHQ-30K datasets, using 9 for all trials. First, second, and third order Taylor expansions for each solver order case is tested. shown in Table the best combination is STORK-4 with first order Taylor expansion.
G.2 IMAGE GENERATION METRICS
In this section, we provide the detailed benchmarks on Fréchet Inception Distance (FID) as reported in the main paper. We further benchmark human preference of the generated samples using HPSv2 (Wu ct al. 2023), and we report the averaged metrics across categorics. Unless otherwise specified, the STORK's hyperparameter is set to 9.
As shown in this section. STORK generally outperforms other fast sampling methods, from noise- predicting flow-mate matching models across image and video generation tasks and datasets. Moreover, STORK demonstrates scalable performance terms of model size and generation resolution. There- fore, experimental results strongly support the superiority of STORK fast sampling method.
20

---

Published as a conference paper at ICLR 2026
Table 6: Effects of order. STORK-I, STORK-2, STORK-3 denote the order of SRK method as the base for STORK, and "Ist", "2nd", and "3rd" denote the Taylor expansion order. FID J for MJHQ-30K using SANA 0.6B (Xie al. 2025) tested. The gray numbers denote the best result in the current order, and the bold numbers denote the absolute best numbers.
Method NFE
9
10
STORK-1-Ist
23.015
21.372
20.799
19.994
STORK-1-2nd
33.856
25,460
22.351
20.913
-13RR
22.768
21.066
20.316
19.593
STORK-2-1st
8.485
7.585
7.417
7.052
STORK-2-2nd
23.602
12.646
8.561
7.227
STORK-2-3rd
14.420
11.790
10,479
9.442
STORK-4-14 -s-1st
7.659
6.667
6.526
6.270
STORK-4-2nd
54.179
30.168
15.254
9.619
STORK-4-3r 20.484
14.976
12.556
10.712
G.2.1 FRÉCHET INCEPTION DISTANCE
Unconditional generation. The benchmark CIFAR-10 (Krizhessky et al., 2009) is presented in table The benchmark on LSUN Beer-Bmmroom (Yu et al., (016) is presented in table For these datasets, the FID is calculated using 50,000 samples, and we use the Inception statistics provided by PNDM (Liu et al. 2022) as reference.
Conditional generation. The benchmark for MJHQ-30K (Li al. 2024) dataset using Pixart-
(Chen et al. 2023) is shown in table The benchmark for MJHQ-30 HJO using SANA-0.6B (Xic
al. 2025) shown in table The benchmark for MJHQ-30K at 1024 resolution using SANA-
1.6B is shown in table II The benchmark for MJHQ-30K using FLUX.I-der (BlackForesl 2024)
shown in table 12 The benehmark for MS-COCS (Lin al. 2015) using
Large (Esser et al. 2024) is presented in table 13 All FIDs are calculated using 30,000 samples
to expedite experiments. For MS-COCO, we randomly subsampled 30,000 images from the entire
validation split. Additionally, the weights for Stable-Diffusion-3355-
and FLUX..1dee are loaded
in torch.bfloatl..
Table 7: Unconditional generation on CIFAR-10 dataset, (DDIM, 32px)
Method \ NFE
9
10
12
15
20
30
50
100
DDIM (Song et al.
23.260
20.390
18.500
15.477
12.848
10.900
8.657
6,990
5,520
DPM-
2025
8.669
7.250
6.471
5.525
4.745
4,015
3.947
3.859
3.851
UniPC Zhao
19.732
17.879
16.666
14.593
12.623
10.843
8.828
7.085
5.685
STORK-4, c=e-2 s=14
5.743
5.497
4.964
4.592
4.168
3.888
3.789
ST0RK-4,<=le3,.5
7.816
7.505
6.077
4.831
3.879
3.337
3.204
3.484
Table 8: Unconditional generation on the L-SUN-Bedrood dataset. (DDIM, 256px)
Method NFF
10
12
20
30
40
50
DDIM
22.164
18.730
16.355
13.228
10.566
8.324
6.768
6.088
5.880
DPM-Solvert+
17.679
14.894
13.101
11.160
8.441
7.488
6.934
6.665
6.491
UniPC
13.665
12.858
12.146
11.213
9.592
8.208
6.481
6.043
5.772
STORK-4, e=le-2, s=24
16.353
11.918
9.594
7.647
6.546
6.107
6.202
6.457
6.658
21

---

Published as a conference paper at ICLR 2026
Table Conditional generation MJHQ-30K. Our method more quickly converges plateau value around 5.5. (Pixart-a, 512px, CFG=4.5)
Method
8
9
10
12
15
20
30
40
50
DEIS Zhang & Chen 2023
9.046
8.162
7.645
6.835
5.535
5.528
5.525
5,540
5.544
DoM-Solver++
8.734
7.859
7.481
6.799
6,485
6.189
5.952
5.846
5.764
UniPC
9.361
8.588
7.912
7.012
6.594
6.260
5,923
5.813
5.738
STORK-4, e=le-2, =424
8.835
6.712
6.035
5.626
5.591
5.535
5.505
5,495
5.508
Table 10: Conditional generation on MJHQ-30K (SANA-0.6B, 512px, CFG=4.5)
Method
7
8
10
12
15
20
30
40
50
Flow-Euler
15.856
13.320
11,851
10.782
9,434
8,430
7.551
6,932
6.637
6.494
Now-DPM-Solvers
9,629
8.390
7.628
7.278
6.774
6,424
6.282
6.095
6.086
6.097
Flow- Flow-UniPC
9.406
8.493
7.801
7.443
6.878
6,452
6.265
6.090
6.085
STORK-4
7.659
6.667
6.526
6.270
5.899
5,851
5.897
5,928
5.989
5,995
Table 11: Conditional generation on MJHQ-30K (SANA-1.68 1024px, CFG=4.5)
Method
8
10
12
20
Flow-Euler 18.348
14.849
12.603
10.980
8.828
7.142
6.069
How-DPMM-SoNee
11.005
9.299
8.288
7.532
6.609
5.929
5,428
Flow-UniPC
10.304
8.909
7.981
7.223
6.412
5.856
5,448
STORK-4
9.804
7.910
6.901
6.258
5.407
5.029
4.879
Table 12: Conditional generation on MJHQ-30K (FL.UX11dde) 512px, CFG=3.5)
Method
7
9
10
12
15
20
Flow-Euler
16,046
14.609
13.968
13.493
12,804
12.588
12.467
Plow-DPM-Solver+t
12.922
12.337
11.925
11.659
11.515
11.569
11.779
Flow-UniPC
12.248
11.874
11.810
11,637
11.595
11.585
12.044
STORK-4
11.456
11.021
10.952
10.978
10.948
11.316
11.705
Table 13: Conditional generation on MS-COCO (SD-3.5-Large, 512px, CFG=3.5)
Method
7
8
10
12
15
20
Flow-Euler
20,748
19.037
18.408
17.975
17.268
17.078
16.858
Hlw-DPMM-Solvers
18.553
18.001
17.633
17.322
17.053
16.889
16.633
Flow-UniPC
19.319
18.220
17.759
17.439
16.954
16.784
16.628
STORK-4
18.722
15.964
15.017
14.661
14.767
15.061
15.291
G.2.2 HPSV2
We further evaluate STORK in terms of human preferences using HPSv2 (Wu al. 2023) benchmark. and we report the averaged score over categories in this section. As shown in table 14fand table [5 generated samples by STORK is consistently preferred over other sampling methods.
22

---

Published as a conference paper at ICLR 2026
Table 14: HPSv2 Benchmark. (Pixart-a: 512px, CFG=4.5)
Method NNFE
8
10
12
15
20
30
40
50
DEIS
27.28
28.09
28.56
28.85
29.27
29.80
29.89
29.97
29.95
29.99
DPM-Solver++
27.60
28.38
28.86
29.12
29.48
29.75
29.90
29.98
29.99
30.02
UniPC
27.33
28.17
28.66
29.04
29.51
29.78
29.94
30.02
30.04
30.05
STORK-4, iFle,,, ==24
27.21
28.00
28.75
29.17
29.66
29.92
30.03
30.07
30.09
30.08
Table 15: HPSv2 Benchmark. (FLUX.1-dex, 512px, CFG=3.5)
Method NFE
7
10
12
20
30
Flow-Euler
28.60
29.09
29.52
29.69
30.13
30.33
30.48
30.64
Flow-DPM-SS0 +
29.34
29.58
29.93
30.09
30.30
30.51
30.60
30.66
Flow-UniPC
29.61
29.77
30.07
30.21
30.39
30.56
30.60
30.67
STORK-4
30.27
30.63
30.84
30.84
31.00
30.95
30.90
30.71
G.3 VIDEO GENERATION METRICS
In this section, we provide more comprehensive results the text-to-video generation using Hunyuan model (Kong al. 2024), up NFE=30. As can be seen, the final score converges for all the methods after 10 NFEs, and STORK generates the best results in terms of the final score. For the sub-metries, our STORK method wins by majority in visual quality and motion quality, and get comparable results in other metrics.
Table 16: EvalCrafter evaluation of Hunyuan Model, with different sampling methods. Scores of the four sub-metrics and the total score are recorded. As shown the table, STORK method consistently outperforms the flow-DPM-Soo DPM-Solver++ and the Hlow-UnIPC methods in terms of the final score for small NFEs, all get similar final score for larger NFEs.
Method NFE
8
9
10
12
20
30
Moo-DDMM-oolvere
45.02
46.42
48.03
49.51
50.74
51.48
52.40
53.11
54.82
54.90
Visual Quality
Flow-UniPC
45.51
47.46
49.23
50.72
51.88
52.66
53.32
54.15
55.08
55.17
STORK (Ours)
50.00
51.91
52.11
52.72
52.68
53.17
53.60
53.74
53.75
53.14
Flow DPM Solvert+
40.90
46.10
47.54
48.79 78
47.83
47.93
50.13
48.97
50.07
49.85
Rext-Video Alignment
Flow-UnIPC
41.38
47.43
46.59
48.76
49.95
50.46
$1.16
50.37
STORK (Ours)
43.18
46.11
46.64
50.09
46.92
48.94
48.01
49.37
50.40
52.44
Hlos-DMM-Solverr
55.35
54.49
54.26
53. 74
53.89 89
53.55
53.
53.54
53.22
53.62
Motion Quality
Flow-UniPC
55.24
$4.56
54.06
54.12 12
53.37
53.51
53.40
53.43
53.32
53.63
STORK (Ours)
55.02
$4.50
54.26
$4.18
$4.06
53.47
$3.79
53.81
$3.95
$3.80
Mlow-DMM.S PPM--Solver++
64.17
63.80
63.48
63.25
63.07
62.92
62.78
62.58
62.
61.85
Temporal Consistency
Flow-UniPC
63,94
63.43
63.20
62.92
62.76
62.61
62.47
62.30
61.97
61.65
STORK (Ours)
61.41
61.67
62.08
62.14
62.26
62.35
62.08
62.12
61.66
61.25
Flow DPM Solvert+
205
211
213
21.5
216
216
219
218
220
220
Final Score
Uni-U
206
212
214
216
215
218
219
220
222
221
STORK (Ours)
210
214
215
219
216
218
217
219
220
221
H ADDITIONAL VISUALIZATIONS
We provide additional visualizations. Generated samples for STORK on CIFAR-10 Krizhevsky et af. (209)) are in Figure 9 and LSUN-Bedroom (Yu et al. 2016) are Figure 10 Visualizations
for MS-COCO-2014 (Lin al.,
2015) are in Figure [] Visualizations for MJHQ-30K (Li al.
2024) are in Figure 12 and Figure T3
Visualizations for Hunyuan (Kong et al.
2024 model video
generation are Figure 14 and Figure
23

---

Published as a conference paper at ICLR 2026
NFE=10
NFE=15
NFE=20
DDIM (ong et al. 2022
DPM-Solver++ Lu et al. 2025
UniPC++ ZZhao et al. 2023
STORK-4
Figure 9: Unconditional generation on CIFAR-10
Krizhevsky et al. 2009 Generated using DDIM
model.
24

---

Published as a conference paper at ICLR 2026
NFE=10
NFE=15
NFE=20
2022
DDIM (Song et al.
OPM-Sollve++ ( et al. 2025)
UniPC+ (Zhao et aL. 2023
STORK
Figure 10: Unconditional generation on L.SUN-Bedroom model.
Yu et al. 2016). Generated using DDIM
25

---

Published as a conference paper at ICLR 2026
Flow-Fuler
Flow- #low-DPM-Solver++
Flow-UniPC
STORK (Ours)
woman dressed as pilgrim in old kitchen." "
"A man taking a photograph skate boarders in the rain."
"Four people gathered in the kitchen looking at a refrigerator.
"A baseball player prepares to swing at a pitch at home plate."
"A black dog has just captured Frisbee using its mouth." "
Figure 11: Comparison between the FloweEuler, Bow-DPM-Soovvet+ Lu et al. 2025 Xie et al. 2025), FloW-UniPC (Zhao et al. 2023), and STORK. Generated using Diffuu--i 3.5-Large (Esser et al. 2024 with 15 NFEs.
26

---

Published as a conference paper at ICLR 2026
Flow-Fuler
Flow-DPM-Soover++
Flow-UniPC
STORK (Ours)
"Darth Vader eating birthday cake in star destroyer interior, cinematic, 50mm". Our image has clearer and lighter background, while other methods have darker backgrounds.
"1920s fantasy interior architecture rebel art, hammer 40K universe, gothic painting". Taogg.Thhe The and Now-DPM-Soffe Solver++ fail to generate a decent image, while the flow-UniPC method has worse visualization quality than the STORK method.
"Watercolor illustration of cute pink and while calf portrait surrounded by small while flowers". Our method has better image quality and more details than other methods.
"3D, photorealistie, photo shot, psychedelic, neon colors, neon lights, palm trees, tall, rectangular hotel building, 1960s cars out front, dranatic, dark sky, rain, sinematie, dramatic Lighting, hyperrealistic, 8k". Our method clearly shows the rain in the background.
Figure 12: More comparisons between the Flow-Fuler, Plow-DPM-Solver+++ (Lu et al. 2025 Xie
et al. 2025). Flow-UniPC (Zhao al.) 2023), and STORK. All images generated using the SANA 1.6 model (Xie et al., (025) at 1024 x 1024 resolution with only 8 NFEs, using prompts from MJHQ-30K. Prompts are displayed beneath each image pair, accompanied by our commentary explaining why STORK's generations are superior.
27

---

Published as a conference paper at ICLR 2026
Flow-Euler
Flov-DPM-Solvert+
Flow-UniPC
STORK (Ours)
"the cockpit of the Roeing 787 Dreaml iner, on night flight, the night sky, clouded, dynamic style, depth Field, F2.8, high Cont 8K, Cinematic 1. lighting, ethereal light, intricate details, ext detailed, incredible details, full colored, complex details, by Weta Digital, Photography, Photoshoot, Shot on 70mm, real photo, einematic, high detailed, HDR, hyper realistic,
"Mangrove emerging from the seed in glass cup.". ".
"French A Alps, photorealisti dramatic 1 lighting, mystique scene, ethereal, majestic, aesthetics, cinematic Lighting, deep focus, super adobe, detailed texture, neuro cognitive art, photoshop, octane render, Pinterest art, awardwinring landscape photography, world renowned, high resolution, color grading, high art, no blurs, ultra wideangle lenses, photo realism, 300 dpi, ultra Quality, 32k "
"Paradise sunner flowers, clipart, sticker, in the style of aquarellist, made of flowers, fernando amorsolo, graphic design elements, vector illustration, in the style of dreamy watercolor florals, floral explosions, thomas w schaller, sticker art, detailed foliage, cute and tropical set of floral pattern, in the style of Luminous watercolors and ink, art deco sensibilities, barbiecore, timeless elegance, floralpunk, colorful arrangements "
Figure 13: More comparisons between the Flow-Euler, al. 2025).. Flow-UniPC (Zhao 2023 ), and STORK on MJHQ-30K, using NFEs.
(Lu et al. 2025
Xie
28

---

Published as a conference paper at ICLR 2026
Flow-UniPC video generation
Hunyuan model KKng al.
2024), 8 NFE
(b) STORK video generation on Hunyuan model
Kong al. 2024 NFE
Figure 14: Video generation comparison with prompt "doctors are constructing a
robot". "
Our video has more than one doctor, and has more realistic scenario.
(a) Flow-UniPC video generation
Hunyuan model
Kong et al. 2024), 8 NFE
(b) STORK video generation on Hunyuan model Kong et al. 2024 NFE
Figure 15: Video generation comparison with prompt "A young woman in her mid-twenties with long blond hair walked briskly down the busy city street, wearing red coat and black boots". Our video has higher quality with a better background.
29


Reviews:
[{"review_id": "YwCGj9JDOF", "reviewer": "Reviewer_Hn4U", "summary": "The paper proposes an advanced ODE sampler, STORK, to handle the stiff dynamics of ODE, while remaining applicable to both score-matching models and flow-matching models. Authors tackle the limitation of prior samplers such as the famous DPM-Solvers, which rely on a semi-linear assumption, thus not natural to be directly applied to flow-matching models. STORK is a variant of stabilized Runge–Kutta (SRK), which is known to handle stiffness well, but it requires expensive intermediate model evaluations. So authors suggest Taylor approximation-based virtual NFEs to approximate internal stages. In experiments, they show strong empirical results across various tasks, including unconditional/conditional image generation and T2V generations, showing consistent improvements under low NFE setups.", "questions": "This is a minor question, but Fig5a looks somewhat unusual, where it shows DPM-Solver++ outperforming UniPC, which contradicts to already reported results from other papers - why is this?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. Clear motivation and message: the authors clearly present what problem they aim to solve and how they approach it.\n2. Conceptually sound method: Leveraging SRK to handle stiffness, and using Taylor-based virtual NFEs to reduce computational cost, finally making a viable and efficient sampler for flow-matching models\n3. Comprehensive empirical evaluation: covering unconditional/conditional image generation and text-to-video tasks, with consistently strong results.\n4. Well-written manuscript: the paper is clear and reader-friendly, e.g., additional explanation on Appendix B about stiffness and numerical analysis, parts like this make the paper accessible to broad readers", "weaknesses": "1. Runtime analysis: authors are only reporting NFE here, but for a more thorough analysis/comparison, they need to report wall clock time and GPU(VRAM) usage. So that readers could better understand in detail, e.g., how much of that time it takes to virtual NFE calculation?\n2. On Table 1,\n    - Clarification on NFE report: since they're comparing samplers, with different count evaluations such as inner NFEs for higher-order/intermediate calculation, or virtual NFEs, and so on... Table 1 is a bit hard to read thoroughly and feels unclear.\n    The authors should consider explicitly separating NFEs for the super-step/sub-step for clearer comparison.\n    - SRK vs STORK: so SRK is super-step (real NFE) + sub-step(real NFE), while STORK is super-step (real NFE) + sub-step(virtual NFE).\n    It would be helpful to show the comparison not just in \"budget-fair\" regime, but in \"algorithm-fair\" regime, e.g., evaluate both methods with the same super/sub-step structure where SRK operates properly (i.e., with sufficient NFE budget).\n    It could strongly back the message “SRK works well but is too expensive, hence we propose STORK to replace those sub-step NFEs with virtual ones.”\n3. Stiffness analysis\n    - Since the \"stiffness\" is one of the main motivations, the actual empirical stiffness analysis on pretrained diffusion/flow matching models would strengthen the argument. While calculating Jacobian could be expensive, even a toy-level experiment or approximate analysis would be highly informative and appreciated.\n\n[minor]\n\n1. I understand the limited page for the main part, but it would be beneficial if the ablation table for the main hyperparameter (Table 3 in the supplement), or pseudocode/algorithm (algo1, algo2) could be moved or at least mentioned in the main paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "UZMiFiPOju", "reviewer": "Reviewer_1FHX", "summary": "This paper introduces STORK, a new type of training-free numerical ODE solver for diffusion and flow models. STORK is based on SRK methods for solving stiff ODEs. To reduce NFEs, STORK further introduces virtual NFEs, which uses finite difference to approximate the derivatives of the velocity field to extrapolate intermediate velocities without extra NFEs. Experiments across various dataset and pretrained models show that STORK outperforms state-of-the-art ODE solvers such as UniPC and DPM++.", "questions": "Please see weaknesses.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The results presented in this paper are solid. Quantitative evaluations demonstrate a clear advantage of STORK over UniPC and DPM-Solver in few-step settings. Qualitative results also look promising: STORK seems to generate more details than other solvers under few-step settings, especially on video generation.\n- Despite the dense math, the presentation quality of this paper is high, making it easily readable. Notably, this paper makes a connection with the notion of stiffness in classical numerical analysis, and then adapts the widely used SRK methods for flow models with a clear motivation.", "weaknesses": "- As mentioned in L300, it is stated that \"naive application of SRK4 to the CIFAR-10 dataset results in very poor sampling results\", so the authors propose to plug Taylor approximation to SRK4. This is an interesting observation, but the reason is not analyzed in depth. If SRK is considered a common method for solving stiff ODEs, why would it not work in this case? Using Taylor expansion and Adams-Bashforth approximation is common in previous flow ODE solvers, so plugging this into SRK seems to make it less unique.\n- I wonder how the methods are categorized in Fig. 4. In particular, UniPC, which works very well in practice, is categorized as not stiff, whereas DPM-Solver is considered stiff. Yet it is clear that UniPC can be considered as an extension of DPM-Solver with an additional corrector. Please check this carefully.\n- Another claimed advantage of the proposed method is structure-independency. However, in most cases, structure-dependency of ODE solvers is not really an issue in practice. It is very trivial to adopt diffusion solvers using exponential integrators for flow models by rescaling the noise schedule and reparametrizing the prediction format.\n- Why is FID getting worse with more NFEs? (fig. 6, Table 6, 10, 11). This behavior is clearly different from DPM-Solver and UniPC, which often converge monotonically. Presenting images generated under higher NFEs could help analyzing this behavior.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DtkhaXvGUK", "reviewer": "Reviewer_jYva", "summary": "The paper introduces STORK, a method for faster sampling from diffusion models by using the Taylor expansion version of a solver (SRK) that is more accurate for stiff ODEs.", "questions": "- Can you clarify why other methods, like the DPM solver, cannot apply to flow matching? What do you mean by the error introduced in each step when you talk about the data prediction step?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The paper introduces a novel method for faster sampling from diffusion models.\n2. The performance of STORK-4 shows improvements when NFE is small.", "weaknesses": "1. The paper motivates the proposed method by claiming diffusion models exhibit stiff dynamics, but this is not properly justified empirically or theoretically. There is insufficient evidence presented that diffusion models are actually stiff, making the motivation unclear. Table 1 shows that the STORK-4 is significantly better than SRK4. More intriguingly, even with 50 NFE of SRK4, its performance is around 6.167 for FID score which is lower than its STORK-4 at 10 NFE. This raises the question of whether SRK4 is a good solver for diffusion models and whether the stiffness assumption is valid.\n2. Theorem 1 establishes the theoretical guarantee of STORK-4 based on its approximation quality to SRK4. However, it is not clear whether the theory actually supports the empirical finding that 10 NFE of STORK-4 can be better than 50 NFE of SRK4. Additionally, since the goal is to achieve fewer NFE in sampling, does the asymptotic rate (provided in Theorem 1) even matter in practice? The constant in front of the rate may matter more for the low NFE regime that this paper targets. Theorem 1 does not help understand the performance of few-step sampling, leaving the empirical gains unexplained by the theoretical analysis.\n3. STORK-4 has many hyperparameters, which makes it unclear whether the observed performance gains come from genuine improvements or from more extensive hyperparameter search. Maybe fix a parameter like s=9 suggested in the paper and report the results for all experiments could help understand the performance gains.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TW6mtq3wYJ", "reviewer": "Reviewer_VACL", "summary": "This paper addresses the challenges of non-straight velocity fields and the dependence on the semi-linear structure in diffusion model ODEs, which make fast but quality-preserving sampling difficult. STROKE tackles these issues by designing a solver that explicitly handles stiffness and structural dependency problems inherent in diffusion and flow-matching models.", "questions": "- In Appendix B, the authors illustrate stiffness through a toy example, but can the same stiffness analysis be quantitatively applied to real diffusion or flow-matching models? A visualization or empirical stiffness profile for these models would help support the theoretical motivation.\n- In Tables 6 and 7, the baseline outperforms STROKE at smaller NFEs, raising the question: how does STROKE behave in 512×512 image generation when NFE ≤ 5? It would be interesting to see whether SRK’s stability advantages appear under such extreme sampling constraints.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The authors redesign the Stable Runge–Kutta (SRK) method to align with the structural characteristics of diffusion and flow dynamics.\n- They propose a general framework applicable to both noise-predicting (diffusion-based) and flow-based generative models, extending beyond the semi-linear assumption used in most prior solvers.\n- The paper combines rigorous theoretical analysis (including convergence and stability proofs) with comprehensive quantitative and qualitative experiments, demonstrating consistent improvements across different datasets and models.", "weaknesses": "- Despite claiming “Stable Runge–Kutta,” the paper does not explore low-NFE (≤5) regimes, where stability advantages would be most pronounced. It remains unclear whether STROKE can outperform existing methods when the step count is drastically reduced.\n- The paper lacks comparisons on larger benchmarks such as ImageNet-256 and ImageNet-512, which are standard for evaluating scalability and visual fidelity at higher resolutions.\n- Although NFEs are matched across methods, a fair comparison in terms of actual computational cost is missing. The paper introduces the concept of *virtual NFEs*—intermediate approximations computed via Taylor expansion and finite-difference—but does not show how much wall-clock time, memory usage, or throughput efficiency is affected. Internal computation overhead (e.g., Taylor expansion, finite-difference updates, velocity caching) incurs additional cost even without extra NFEs. The paper abstracts this under “virtual NFEs,” but a detailed analysis of real runtime and resource usage would strengthen the claim of efficiency.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)