# Randomized stable SDIRK sampler

- ID: 1129
- Canonical URL: https://synthcore.org/idea/1129/randomized-stable-sdirk-sampler
- API JSON: https://synthcore.org/api/idea/1129.json
- API Markdown: https://synthcore.org/api/idea/1129.md
- Verification status: unverified
- Source: [arXiv:2607.18928](https://arxiv.org/abs/2607.18928)
- Category: dynamics
- Solves: stability, accuracy, sample-efficiency
- ML areas: diffusion-sampling, ssm, training-dynamics
- Math tags: numerical-analysis, dynamical-systems, stochastic-processes, linear-algebra
- Ratings: usefulness 6/10; difficulty 7/10; novelty 6/10

## Idea description

Replace the explicit Euler, Heun, or fixed-step midpoint update used for a neural ODE or diffusion probability-flow trajectory with a two-stage randomized SDIRK step. Draw one random scalar per time step, use it in both implicit stage equations, and solve each stage with Newton or damped fixed-point iteration. The randomness targets quadrature error caused by nonsmooth score networks, while the singly diagonal structure permits reuse of the same Jacobian preconditioner for both stage solves.

## Mathematical statement

The paper studies the initial-value problem $u'(t)=f(t,u(t))\in\mathbb{R}^d$, $u(0)=u_0$, on the grid $t_n=nh$ with $h=T/N$. Its randomized SDIRK stages use a common diagonal coefficient, so each stage linearization contains the resolvent $(I-hx f_u(\theta^2,\ell^2))^{-1}$, where $x$ is the randomized diagonal coefficient, $f_u$ is the Jacobian with respect to the state, and $(\theta^2,\ell^2)$ is an intermediate time-state pair. Use the two-stage singly diagonal template $A(x)=\begin{bmatrix}x&0\\1-2x&x\end{bmatrix}$, $b=(1/2,1/2)^T$, and $c(x)=(x,1-x)^T$: $k_1=f(t_n+xh,u_n+hxk_1)$ and $k_2=f(t_n+(1-x)h,u_n+h[(1-2x)k_1+xk_2])$. The update is $u_{n+1}=u_n+\frac h2(k_1+k_2)$. Algebraic stability is checked through $b_i\geq0$ and $M=BA+A^TB-bb^T\succeq0$, where $B=\operatorname{diag}(b_1,b_2)$. The paper's convergence result is root-mean-square order $2.5$, meaning $(\mathbb{E}\|u(T)-u_N\|^2)^{1/2}=O(h^{2.5})$ under its low-regularity assumptions.

## Key formulas

- $$u^{\prime}(t)=f(t,u(t))\in\mathbb{R}^{d},\qquad u(0)=u_{0},\qquad t_n=nh,\quad h=T/N.$$
- $$k_1=f(t_n+xh,u_n+hxk_1),\qquad k_2=f\bigl(t_n+(1-x)h,u_n+h[(1-2x)k_1+xk_2]\bigr).$$
- $$u_{n+1}=u_n+\frac{h}{2}(k_1+k_2),\qquad (I-hx f_u(\theta^2,\ell^2))^{-1}\ \text{is the shared diagonal-stage resolvent}.$$
- $$M=BA+A^TB-bb^T\succeq0,\qquad A=\begin{bmatrix}x&0\\1-2x&x\end{bmatrix},\quad b=\begin{bmatrix}1/2\\1/2\end{bmatrix},\quad B=\operatorname{diag}(b).$$

## Implementation notes

Integrate this at the trajectory-step operation, not inside the score or vector-field network. For a state z_n at time t_n, draw x once per step from the coefficient distribution used by the paper; for an MVP, use x uniformly on [0.25,0.5] and sweep this range to verify algebraic stability. Compute k1 by solving k1=f_theta(t_n+xh,z_n+hx k1). Compute k2 by solving k2=f_theta(t_n+(1-x)h,z_n+h[(1-2x)k1+xk2]), then return z_{n+1}=z_n+h(k1+k2)/2. Use two or three Newton iterations with Jacobian-vector products, or Anderson-accelerated fixed-point iteration; do not form a dense Jacobian. A Newton linear solve should use conjugate gradients, GMRES, or a low-rank preconditioner for the operator I-hx J_theta, matching the paper's resolvent structure. For diffusion models, apply the method to the reverse-time probability-flow ODE; for neural ODE training, use it in the forward solver with checkpointing or differentiation through the unrolled iterations. Compute the paper-derived quantities exactly: A(x), b, c(x), and M. Estimate convergence by comparing against a very small-step reference trajectory over 32 independent random x-sequences, reporting RMS endpoint and path error. The first cheap experiment should use a small score model on CIFAR-10 or a synthetic stiff vector field containing a ReLU kink, comparing Euler, Heun, deterministic SDIRK, and randomized SDIRK at equal network evaluations. Success means lower trajectory error or better FID at fixed evaluations, together with fewer divergences in regimes where explicit Euler is unstable.

## Disclaimer

AI-generated research hypothesis, automatically tested. Not peer-reviewed.