[rk-mlu]

First slice of milestone M1 Core and Euler–Maruyama (see the roadmap in docs/development/architecture.md): everything from the problem definition to the solution, tested end to end on geometric Brownian motion and a linear system with multidimensional noise.

What is added

Area	Objects
Problems	ItoSDE (base SDEProblem), NoiseStructure (declared and verified flags), Calculus
Noise	StepNoise, NoiseSource; BrownianPath is now a noise source (steps(requires))
Schemes	Scheme, OneStepScheme, EulerMaruyama, matvec
Solver	Solver (checks scheme, noise and coefficient shapes before the loop; save="all" | "final"), Solution
Benchmarks	GeometricBrownianMotion, LinearSystem, registry numsde.benchmarks.PROBLEMS, get_problem
Backend	asarray_like (public), asarray_in, device_of
Docs	theory pages schemes/euler_maruyama, problems/geometric_brownian_motion, problems/linear_system; API reference; changelog

The part of the API sketch up to solve now runs as tests/unit/test_api_sketch.py.

Decisions (recorded in the decision log)

• A noise source iterates over its own grid (steps(requires) instead of steps(grid)); other step sizes come from BrownianPath.coarsen.
• A problem carries its initial time t0; the noise grid must start there, and exact_solution(t, W) receives W(t) − W(t0). drift and diffusion receive t as a Python float.
• Given a noise source, the solver runs in its namespace, dtype and device; an array x0 must agree with it. Without a source, solve(times=..., n_paths=..., seed=...) samples a Wiener path following x0 or explicit backend, dtype, device.
• LinearSystem is built from a common eigenbasis, so its matrices commute (commutativity recorded as verified) and the exact solution needs no matrix exponential and no SciPy.
• Ruff accepts the en dash in names such as Euler–Maruyama.

Tests

• Unit tests: hand-computed EM steps (scalar and multidimensional noise), all solver checks and error messages, benchmark exact solutions against a truncated exponential series.
• Property tests: batch independence, reduction to explicit Euler without noise, exactness for dX = dW, coarsened paths as the same realization.
• Statistical tests (slow): fitted strong order 0.5 ± 0.1 on both problems (error uniformly on the grid, p = 2), weak error of E[X_h^N] for GBM against its closed form.
• PyTorch: NumPy and PyTorch solutions agree; tensor x0 selects the backend; a CUDA test with the gpu marker.
• Coverage 100 % lines and branches, with and without PyTorch installed.

Please check

• Section references on the EM theory page (8441325): Kloeden–Platen Theorem 10.2.2 covers the final-time error in L^1, the L^2 final-time error now cites Milstein–Tretyakov (2021), Theorem 1.1.1; Section 14.1 (weak order) confirmed.
• Strong order 1/2 uniformly on the grid for p = 2: Mao (2008), Chapter 2, Theorem 7.3. Divergence for superlinear coefficients: Hutzenthaler–Jentzen–Kloeden (2011), Theorem 2.1, checked against the article.
• Kanagawa (1988), Yokohama Math. J. 36, 79–86, for p ≥ 2: not yet consulted (journal repository down for maintenance); marked † and kept out of references.bib. Deferred to Check Kanagawa (1988) for the strong order of Euler–Maruyama in L^p, p ≥ 2 #14; not a merge blocker.
• Mao2008: 2nd ed., Horwood Publishing, Chichester, 2008, ISBN 978-1-904275-34-3, confirmed by the maintainer (8661bb9).
• HutzenthalerJentzenKloeden2011: Proc. R. Soc. A 467(2130), 1563–1576, 2011 (printed issue, as in MR2795791; online first December 2010), confirmed by the maintainer (69e7d5a).
• Default parameters of geometric Brownian motion (μ = 0.05, σ = 0.2, as in the API sketch), confirmed by the maintainer.
• Default parameters of the linear system (see its theory page), confirmed by the maintainer.
• The GPU test tests/unit/test_solver_torch.py::test_cuda_solution_matches_numpy on your machine (nox -s gpu): passed on 2026-11-03 together with test_brownian_torch.py::test_cuda_path_matches_numpy (NVIDIA RTX A4000, driver 595.99.02, PyTorch 2.15.0+cu132, Python 3.14).

Next in M1

Streaming noise source, StrongConvergenceStudy with pathwise bootstrap, Brownian motion on the circle, first tutorial, ports of the GBM legacy scripts.

https://claude.ai/code/session_018KaTyMUfivRVAGYQR6kGnp

---

Generated by Claude Code

[rk-mlu]

Review of commit 690cbc0

Recommendation: address the following two findings before merging.

[P1] Internally sampled noise bypasses problem compatibility checks

In solver.py, lines 154–158, _check_noise is called only when the caller supplies a noise source. Internally sampled paths therefore bypass the initial-time and x0 namespace/dtype/device checks.

Reproduced example:

import numsde as ns

problem = ns.benchmarks.GeometricBrownianMotion(mu=1.0, sigma=0.0, t0=2.0)
solver = ns.Solver(problem, ns.schemes.EulerMaruyama())
grid = ns.uniform_grid(1.0, 4)
solution = solver.solve(times=grid, seed=0)
# Succeeds on [0, 1], although the initial condition is prescribed at t0=2.
# solution.final == [[2.44140625]]

path = ns.WienerProcess(1).sample_path(grid, seed=0)
solver.solve(path)
# Correctly raises ValueError: grid starts at 0.0, but t0=2.0.

A separate reproduction with a custom ItoSDE and a float64 array x0 confirmed that solve(times=..., dtype=np.float32) silently downcasts the initial condition, contrary to the documented requirement that an array x0 agree with the noise source.

Suggested fix: Separate validation of mutually exclusive sampling arguments from noise/problem compatibility validation. Run the latter for both supplied and internally sampled noise. Add regression tests for mismatched initial times and array x0 dtype/backend/device.

[P2] The stated assumptions do not guarantee strong order 1 for additive noise

The EM theory page, lines 48–66 assumes only 1/2-Hölder continuity in time. Under that assumption, additive noise and smooth drift do not suffice for strong order 1.

For example, on [0,1], consider

$$ \mathrm{d}X(t)=\sqrt{t},\mathrm{d}W(t),\qquad X(0)=0. $$

This has additive noise, zero (hence smooth) drift, and a diffusion coefficient satisfying the stated time regularity. By the Itô isometry, on a uniform grid with h=1/N the squared terminal error is

$$ \mathbb{E}|X(1)-X_h^N|^2 =\sum_{n=0}^{N-1}\int_{nh}^{(n+1)h} (\sqrt{t}-\sqrt{nh})^2,\mathrm{d}t \asymp h^2\log(1/h). $$

Thus even the terminal L2 error is not O(h). The spatial Milstein correction vanishes, but that does not eliminate the error from freezing a time-dependent diffusion coefficient.

Suggested fix: Explicitly add sufficient time regularity of g (for example, Lipschitz continuity in time), alongside suitable drift regularity, to the additive-noise order-1 claim. Document this qualification for expected_strong_order as well; additive=True alone is not a regularity certificate.

Validation

Check	Result
Regular suite including doctests	287 passed; 100% line and branch coverage
Statistical tests	7 passed
Ruff lint and formatting	Passed
Codespell	Passed
Offline Sphinx build with warnings treated as errors	Passed
Mypy, rerun with --no-incremental	Passed, 24 source files

The initial cached mypy invocation encountered an internal tool error; the clean non-incremental run succeeded.

PyTorch and CUDA tests were not run because PyTorch was not installed in the review environment. The bibliography/section references marked † remain unverified. No source changes were made during this review.

[rk-mlu]

Thanks for the review. Both findings are accepted and addressed in 9057fec.

P1: sampled noise bypassed the compatibility checks

Confirmed with both reproductions. solve now separates the two concerns: _check_no_sampling_options rejects sampling options combined with a supplied source, and _check_compatibility (noise dimension, grid starting at t0, namespace, dtype and device of an array x0) runs for supplied and sampled noise. Explicit backend, dtype or device must now agree with an array x0; the earlier test test_sampling_follows_an_array_initial_value had encoded the silent downcast and was changed accordingly.

Regression tests (tests/unit/test_solver.py, tests/unit/test_solver_torch.py): grid not starting at t0 (sampled and supplied), dtype mismatch in both directions, plain NumPy vs. array_api_compat.numpy, NumPy noise for a tensor x0. The solve docstring has a doctest for the t0 case.

P2: assumptions for strong order 1 with additive noise

Confirmed; the closed form of your example gives E|X(1) − X_h^N|²/h² = 0.74, 0.97, 1.20, 1.43, 1.66 for N = 2^4, 2^8, …, 2^20, growing like log(1/h). The EM theory page now requires g Lipschitz continuous in t and f twice continuously differentiable in x with bounded derivatives and Lipschitz in t for order 1, explains why (the freezing error of g), and includes your example. The docstrings of EulerMaruyama and expected_strong_order state that the noise-structure flags certify no time regularity. The reference for the order-1 statement is marked † until the literature check.

Both points are recorded in the decision log. No errata entry, since the code was never merged.

Validation

Check	Result
Regular suite with PyTorch (CPU)	306 passed, 100 % lines and branches
Regular suite without PyTorch	293 passed, 4 skipped, 100 %
Statistical tests	7 passed
Ruff, mypy, codespell, Sphinx -W	passed
CI on 9057fec	all jobs green

Still open: the CUDA test (maintainer's machine) and the references marked †.

---

Generated by Claude Code

[rk-mlu]

Follow-up review of commit 9057fec

Both findings from the previous review are substantively addressed:

• P1: Supplied and internally sampled noise now pass the same compatibility check. Regression tests cover mismatched initial times, dtype mismatches, and backend conflicts.
• P2: The theory page now explicitly states sufficient time regularity for strong order 1 with additive noise. The qualification is also documented for expected_strong_order.

Remaining documentation fix: literal control characters in LaTeX

In docs/theory/schemes/euler_maruyama.md, lines 68–81, some LaTeX commands have been stored as control characters:

Intended command	Actual contents
\big	U+0008 (backspace) followed by ig
\vert	U+000B (vertical tab) followed by ert
\asymp	U+0007 (bell) followed by symp

Affected lines: 68 and 80 (big), 79 (vert, twice), and 81 (asymp).

Please replace these sequences with the literal LaTeX commands so that the formulas render correctly. The intended displayed equation is:

\mathbb{E}\,\vert X(1) - X_h^N \vert^2
= \sum_{n=0}^{N-1} \int_{nh}^{(n+1)h} \big(\sqrt{t} - \sqrt{nh}\big)^2\,\mathrm{d}t
\asymp h^2 \log(1/h).

Validation and conclusion

All 12 GitHub checks for the reviewed commit are successful, including the CPU PyTorch tests, documentation, lint/type checks, wheel tests, and the aggregate CI passed check. Tests were not rerun locally during this follow-up because the execution environment was unavailable; this assessment is based on inspection of the changes, regression tests, raw documentation content, and GitHub check results.

After the small documentation correction above, this follow-up review identifies no further merge blockers. The previously outstanding bibliography/section verifications remain unverified. No source changes were made.

[rk-mlu]

Thanks, confirmed and fixed in 2132d7f.

The three commands had passed through a non-raw string in an editing script, where \b, \v and \a are escape sequences. They are now the literal \big, \vert and \asymp; the displayed equation reads exactly as in your comment. A scan of all tracked text files found no other control characters, and the built HTML contains the commands.

Since sphinx-build -W did not notice the problem, tests/unit/test_repository_text.py now rejects control characters other than tab, line feed and carriage return in all text files of the repository. Run against the previous commit, it reports exactly the seven occurrences of your table (lines 68, 79, 80 and 81). It splits lines at \n only, because str.splitlines() treats U+000B as a line break and would have hidden the \vert cases. The rule is recorded in the decision log.

All 12 checks pass on 2132d7f. #13, stacked on this branch, has been rebased and is green as well.

---

Generated by Claude Code