IS  //  Input Systems

Notes · debugging

An unbounded recursion in torchsde: why Stable Audio Open's default scheduler crashes on CPU, MPS, and AMD ROCm

Stable Audio Open's diffusers pipeline dies with a RecursionError deep inside torchsde on CPU, MPS, and AMD ROCm. Here's the real mechanism, why it isn't a CUDA dependency, and the one-line call-site fix.

Jason · June 2026

Versions used throughout: diffusers 0.38.0 · torchsde 0.2.6 · torch 2.12.0

If you try to run Stable Audio Open's diffusers pipeline on Apple Silicon (MPS or CPU) or AMD ROCm, there's a good chance it never produces a single sample. Instead, after the model loads and denoising begins, the process dies with a RecursionError deep inside torchsde, a library you may not even realize the pipeline depends on.

This post walks through what actually happens, why the common "it only works on CUDA" framing is misleading, and a small call-site patch that unblocks the pipeline on these backends.

1. The problem: running the pipeline off CUDA

StableAudioPipeline("stabilityai/stable-audio-open-1.0") works in the officially published example, which runs on NVIDIA CUDA (pipe = pipe.to("cuda"), torch.Generator("cuda")). The natural assumption when it crashes elsewhere is that there's a CUDA dependency somewhere: a custom kernel, a fused op, an fp16 path that only NVIDIA supports.

That assumption is wrong, and it's worth being precise about why, because it changes where you look for the fix.

The crash lives entirely in torchsde, which is pure Python. There is no kernel involved, no device-specific code path, nothing that NVIDIA hardware does differently at the level where the failure occurs. The mechanism is device-independent in principle.

What we can say empirically:

What we cannot honestly claim: that NVIDIA CUDA is immune. We did not run the official CUDA example on NVIDIA hardware, and no source states that CUDA works while other backends fail. The most defensible explanation for the working official example is simply that CUDA is the only path the authors published and tested, not that the bug requires non-CUDA hardware to manifest. Per-device floating-point rounding could land queries differently and dodge the recursion, but that is unverified. The fp16-vs-fp32 story isn't clean either: the AMD repro used fp16 and crashed, ours used fp32 and crashed.

So the accurate framing is this: a scheduler-bounds bug that surfaces reliably on these backends, not a hardware limitation. Treat "CUDA works" as an observation about test coverage, not a property of the bug.

2. The symptom and how we traced it

The visible symptom is a RecursionError: maximum recursion depth exceeded, with a traceback that bottoms out repeating through a handful of torchsde/_brownian/brownian_interval.py frames (_split → _split_exact → _Interval.__init__ → _split …). Just before it diverges, torchsde emits a telling warning:

UserWarning: Should have ta>=t0 but got ta=0.0 and t0=0.3.

Two things about this are diagnostic:

1. Raising sys.setrecursionlimit does nothing. This is the key clue. If the recursion were merely deep, a higher limit would let it finish. It doesn't, because the recursion is unbounded: there is no base case the descent can ever reach. A larger limit just postpones the same crash.

2. The warning names the culprit directly. ta=0.0 and t0=0.3 says a Brownian query is being made at time 0.0 against a tree whose domain starts at 0.3. The query is out of range. Everything downstream follows from that.

We traced it by running the actually-installed diffusers 0.38.0 / torchsde 0.2.6 source step by step against a real cosine exponential schedule, reproducing both the warning and the non-terminating traceback. The whole thing reproduces on CPU in well under a second with no model weights, no UNet, and no network, using just the scheduler's noise sampler and a fabricated final-step query. A self-contained script is in section 4.

3. Technical deep-dive

To see why a query at sigma=0 is fatal, you need three pieces: which scheduler is actually running, what its sigma schedule looks like at the tail, and how torchsde's halfway_tree split behaves when handed a target that rounds onto its own boundary.

3.1 The scheduler and its noise sampler

Despite StableAudioPipeline type-hinting EDMDPMSolverMultistepScheduler, stable-audio-open-1.0 loads CosineDPMSolverMultistepScheduler at runtime per its config, so trace the bug through that scheduler, not the pipeline's signature.

This is an SDE-style multistep solver. On each step(), it injects a Brownian noise increment between consecutive sigma levels. It builds the noise sampler lazily, once, and then queries it on every step with the current and next sigma:

# scheduling_cosine_dpmsolver_multistep.py, inside step(), ~lines 656-664
if self.noise_sampler is None:
    ...
    self.noise_sampler = BrownianTreeNoiseSampler(
        model_output,
        sigma_min=self.config.sigma_min,   # 0.3 for Stable Audio Open
        sigma_max=self.config.sigma_max,   # 500
        seed=seed,
    )
noise = self.noise_sampler(self.sigmas[self.step_index],
                           self.sigmas[self.step_index + 1])

BrownianTreeNoiseSampler (in scheduling_dpmsolver_sde.py) wraps a torchsde.BrownianInterval. For the cosine scheduler, no transform is passed, so it defaults to the identity. The tree's domain is therefore literally [sigma_min, sigma_max] = [0.3, 500], and the values queried are the sigmas themselves. The interval is constructed with tol=1e-6, pool_size=24, and crucially halfway_tree=True:

# scheduling_dpmsolver_sde.py, BatchedBrownianTree.__init__
torchsde.BrownianInterval(t0=0.3, t1=500, tol=1e-6,
                          pool_size=24, halfway_tree=True)

tol=1e-6 makes BrownianInterval round all query times to 6 decimal places (self._round = lambda x: round(x, 6)). halfway_tree=True forces every interval split down the dyadic-bisection branch. Hold onto both facts.

3.2 The sigma schedule, and where the zero comes from

set_timesteps builds the schedule. For Stable Audio Open: sigma_schedule="exponential", sigma_min=0.3, sigma_max=500. The exponential branch computes:

sigmas = torch.linspace(log(sigma_min), log(sigma_max), num).exp().flip(0)

in float32, descending from ~500 down to sigma_min. The last nonzero entry is exp(log(0.3)), which in float32 does not round-trip to 0.3: it evaluates to 0.30000001192092896, i.e. a few ULPs above 0.3.

Then the scheduler's default final_sigmas_type="zero" appends a literal zero:

# set_timesteps, final_sigmas_type == "zero"
sigma_last = 0
self.sigmas = torch.cat([sigmas, sigmas.new_tensor([0.0])])

So the tail of the schedule is:

[..., 0.3485, 0.3233, 0.30000001192092896, 0.0]
                       └ sigmas[-2]          └ sigmas[-1]

On the final denoising step, step() therefore queries the sampler with (sigmas[-2], sigmas[-1]) = (0.30000001192…, 0.0). That 0.0 is the appended final-zero sigma, and it sits below the tree's lower bound of 0.3. This is the whole bug in one line.

(Switching to final_sigmas_type="sigma_min" does not rescue you. The final query becomes (0.30000001192…, 0.3); after tol-rounding both collapse to 0.3, the interval has ~zero width, and BrownianTreeNoiseSampler.__call__ divides the increment by (t1 - t0).abs().sqrt() ≈ 0, yielding NaN/degenerate noise. You trade an obvious crash for silent garbage audio.)

3.3 Inside torchsde: why the bisection never terminates

Now follow the out-of-range query into BrownianInterval.

BatchedBrownianTree.__call__ first sorts the pair, so (0.30000001192…, 0.0) becomes t0=0.0, t1=0.30000001192…, and calls the tree with (0.0, 0.30000001192…).

In BrownianInterval.__call__, ta=0.0 is less than self._start=0.3. So it emits the Should have ta>=t0 warning and clamps ta up to 0.3. But tb=0.30000001192… stays as-is, sitting above _start and below _end. Critically, after clamping, ta (0.3) != tb (0.30000001192…), so the code does not take the ta == tb zero-return fast path. It proceeds into _loc(0.3, 0.30000001192…).

_loc_inner sees ta == self._start with no midway recorded yet, so it calls self._split(tb) with tb = 0.30000001192….

Here is the trap, in the halfway_tree branch of _split:

# brownian_interval.py, _split (halfway_tree branch), ~lines 321-330
self._split_exact(0.5 * (self._start + self._end))   # midpoint = 250.15
...
if midway < self._midway:
    self._left_child._split(midway)                  # recurse left

The first split creates children [0.3, 250.15] and [250.15, 500]. The target 0.30000001192… is less than 250.15, so it recurses into the left child [0.3, 250.15]. Its midpoint is 125.2; the target is below that, so left again. Then [0.3, 62.7], [0.3, 31.5], and so on:

[0.3, 500]
   └ [0.3, 250.15]
        └ [0.3, 125.2]
             └ [0.3, 62.7]
                  └ [0.3, 31.5]
                       └ ...  forever

The left edge is pinned at 0.3 and the target 0.30000001192… sits below every halfway point, so the descent always takes the left branch. It can never exit, because:

That's why this is unbounded recursion, not deep recursion. There is no base case to reach. sys.setrecursionlimit cannot help by construction. The same shape explains the AMD report's upper-bound variant: fp32 endpoint reconstruction can drift a few ULPs past 500, putting tb=500.00006… outside t1=500.0 and driving the symmetric runaway split at the top of the tree.

The root cause, stated plainly: the caller built the Brownian tree over [config.sigma_min, config.sigma_max] but then queries it outside that domain. BrownianTreeNoiseSampler is only valid for queries inside the interval it was constructed over; here the contract is violated by the scheduler, not by torchsde.

4. The fix

The correct place to fix this is the call site in CosineDPMSolverMultistepScheduler.step(), not inside BrownianTreeNoiseSampler and not inside torchsde. Derive the tree's bounds from the actual sigma schedule the scheduler will query, rather than from config. This is exactly what the open upstream PR (#13754) does:

--- a/src/diffusers/schedulers/scheduling_cosine_dpmsolver_multistep.py
+++ b/src/diffusers/schedulers/scheduling_cosine_dpmsolver_multistep.py
@@ -653,10 +653,16 @@ def step(
                 seed = (
                     [g.initial_seed() for g in generator] if isinstance(generator, list) else generator.initial_seed()
                 )
+            # Use the actual extrema of `self.sigmas` rather than the config bounds:
+            # the Karras/exponential reconstruction of the endpoints in fp32 can drift
+            # by a few ULPs, and `final_sigmas_type="zero"` makes `sigmas[-1] == 0`,
+            # both of which fall outside `[config.sigma_min, config.sigma_max]`. An
+            # out-of-range query drives `torchsde` into unbounded recursive splitting
+            # of its Brownian interval and eventually raises `RecursionError` (#13274).
             self.noise_sampler = BrownianTreeNoiseSampler(
                 model_output,
-                sigma_min=self.config.sigma_min,
-                sigma_max=self.config.sigma_max,
+                sigma_min=self.sigmas.min().item(),
+                sigma_max=self.sigmas.max().item(),
                 seed=seed,
             )

With the identity transform, the queried values are exactly self.sigmas, whose minimum is 0 under final_sigmas_type="zero" and whose maximum is the true schedule max (covering the fp32 drift past 500). Building the tree on [0.0, 500.00006…] makes the final-step query (0.3, 0.0) land inside the domain by construction, which is precisely the invariant torchsde requires.

Why this is minimal and correct

Why not the obvious alternatives

To be clear about credit: this is not a novel discovery. The maintainers' issue tracker already has it. PR #13754 (open, from an external contributor, awaiting review as of this writing) makes exactly this change and adds a regression test, test_step_does_not_recurse_with_zero_final_sigma; PR #13515 was an earlier attempt at the same idea. What we add is an independent reproduction on a different backend and precision (Apple Silicon CPU/MPS, fp32) than the original AMD/fp16 report, which helps confirm the root cause is the schedule bounds rather than anything device-specific.

A self-contained reproduction

No model weights, no UNet, no network. Verified on diffusers 0.38.0 / torch 2.12.0 / torchsde 0.2.6, CPU.

"""
Reproduces the RecursionError from CosineDPMSolverMultistepScheduler's
Brownian noise sampler when final_sigmas_type="zero" (its default), then
shows the call-site fix. Runs on CPU in well under a second.
"""
import sys
import torch
from diffusers.schedulers.scheduling_dpmsolver_sde import BrownianTreeNoiseSampler

SIGMA_MIN = 0.3     # CosineDPMSolverMultistepScheduler default (Stable Audio Open)
SIGMA_MAX = 500.0
FINAL_SIGMA = 0.0   # final_sigmas_type="zero" -> last sigma is 0


def fresh_latent():
    return torch.zeros(1, 64, 1024)  # only shape/dtype/device matter


def reproduce_crash():
    # Exactly how step() builds the sampler today.
    sampler = BrownianTreeNoiseSampler(
        fresh_latent(), sigma_min=SIGMA_MIN, sigma_max=SIGMA_MAX, seed=0
    )
    start = float(sampler.tree.trees[0]._start)
    print(f"  tree domain = [{start}, {float(sampler.tree.trees[0]._end)}]")
    print(f"  final-step query = ({SIGMA_MIN}, {FINAL_SIGMA})  "
          f"({FINAL_SIGMA} is below {start})")
    sampler(torch.tensor(SIGMA_MIN), torch.tensor(FINAL_SIGMA))  # diverges


def build_fixed_sampler():
    # Mirrors the call-site fix: bounds from the real schedule, whose min is 0.
    sigmas = torch.tensor([SIGMA_MAX, 50.0, 5.0, SIGMA_MIN, FINAL_SIGMA])
    return BrownianTreeNoiseSampler(
        fresh_latent(),
        sigma_min=float(sigmas.min()),  # 0.0
        sigma_max=float(sigmas.max()),  # 500.0
        seed=0,
    )


def show_fix():
    sampler = build_fixed_sampler()
    print(f"  tree domain = [{float(sampler.tree.trees[0]._start)}, "
          f"{float(sampler.tree.trees[0]._end)}]")
    w = sampler(torch.tensor(SIGMA_MIN), torch.tensor(FINAL_SIGMA))
    print(f"  final-step query succeeded: shape={tuple(w.shape)}, "
          f"std={float(w.std()):.4f}, finite={bool(torch.isfinite(w).all())}")


def main():
    print(f"sys.recursionlimit = {sys.getrecursionlimit()}\n")
    print("[1] Default sampler (current diffusers behavior):")
    try:
        reproduce_crash()
    except RecursionError as e:
        print(f"  REPRODUCED -> RecursionError: {e}\n")
    else:
        print("  did NOT crash (bug not reproduced in this environment)\n")
        return 1
    print("[2] Sampler built from the real sigma schedule (the fix):")
    show_fix()
    print("\nOK: the patch eliminates the RecursionError and yields finite noise.")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Observed output (CPU):

sys.recursionlimit = 1000

[1] Default sampler (current diffusers behavior):
  tree domain = [0.3, 500.0]
  final-step query = (0.3, 0.0)  (0.0 is below 0.3)
  (UserWarning: Should have ta>=t0 but got ta=0.0 and t0=0.3.)
  REPRODUCED -> RecursionError: maximum recursion depth exceeded

[2] Sampler built from the real sigma schedule (the fix):
  tree domain = [0.0, 500.0]
  final-step query succeeded: shape=(1, 64, 1024), std=0.9997, finite=True

5. Impact

This unblocks Stable Audio Open's diffusers pipeline on Apple Silicon (CPU and MPS) and AMD ROCm, where the default scheduler currently crashes before producing any audio. Until the upstream fix lands and ships in a release, the same effect can be obtained by deriving the Brownian tree's bounds from self.sigmas as above (or, for the cosine scheduler's identity transform, by widening the tree's lower bound to 0).

Two honest caveats worth keeping attached to any claim here:

1. The "works" output is non-degenerate, not perceptually validated. With the fix, the final-step Brownian increment is finite with sensible variance: the self-contained script above reports std ≈ 1.0. That confirms the recursion is gone and the increment isn't NaN. It is not a statement about end-to-end perceptual audio quality, which depends on the full pipeline and weights.

2. We make no claim about CUDA. We did not run this on NVIDIA hardware. The most we'll say is that CUDA appears to be the published-and-tested path, which is a statement about coverage, not immunity. Do not read this as "CUDA is required" or "CUDA cannot hit this." The recursion is pure-Python and device-independent in principle, and the only public third-party repro is on AMD ROCm.

The broader lesson is mundane but recurring in diffusion stacks: a Brownian sampler is only valid over the domain it was built for, and schedules that legitimately integrate down to sigma = 0 will query right at (or just past) their declared bounds. Floating-point endpoint reconstruction makes "just past" the common case. Derive sampler bounds from the schedule you will actually traverse, not from the config you declared.

Reproductions were run against the as-installed diffusers 0.38.0, torchsde 0.2.6, and torch 2.12.0 source on CPU and Apple Silicon MPS. The fix shown is diffusers PR #13754; credit for the upstream patch is theirs.