Fix SROCK2 NoiseWrapper and non-square noise (#3188, #3170) - #3468
Fix SROCK2 NoiseWrapper and non-square noise (#3188, #3170)#3468singhharsh1708 wants to merge 1 commit into
Conversation
4134cbe to
8480565
Compare
|
@oscardssmith any reviews on this fix ? |
| # Similarly tᵢ₋₂ = tₛ₋₂, tᵢ₋₁ = tₛ₋₁, tᵢ = tₛ | ||
|
|
||
| if (W.dW isa Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob) | ||
| if (W.dW isa Number) || is_diagonal_noise(integrator.sol.prob) |
| function init_χ!(vec_χ, W) | ||
| r = rng(W) | ||
| for i in eachindex(vec_χ) | ||
| vec_χ[i] = 2 * floor(rand(r) + 0.5) - 1 |
There was a problem hiding this comment.
this assumes Float64, also isn't gpu compatible.
|
Needs convergence tests, I don't think the algorithm looks right. |
243293e to
157e080
Compare
|
The length(W.dW) == 1 check was removed because it caused incorrect behavior for non-diagonal noise with a single noise source (e.g. noise_rate_prototype = ones(n, 1)). In that case length(W.dW) == 1 is true but the noise is non-diagonal, so taking the scalar path u += Gₛ .* W.dW fails — Gₛ is n×1 and adding it to u (length-n vector) causes a DimensionMismatch. The mul! path is always correct for non-diagonal noise regardless of m. The is_diagonal_noise check already handles the scalar path correctly. |
|
DimensionMismatch for non-square noise (noise_rate_prototype = zeros(n, m)) |
|
Is the algorithm correct for non diagonal though? Or does it need levy area calculations? What does the paper say? |
|
The key paper is Abdulle, Vilmart, Zygalakis (2013) SIAM J. Sci. Comput. 35(4). Eq. (17) defines J_{q,r} = (ξ_q ξ_r ± χ_q) · h/2 for q ≠ r using Rademacher χ variables — those are exactly the terms commented out at line 461. For our PR's use case (noise_rate_prototype = zeros(n, 1), m=1), the q ≠ r cross terms never appear (loop runs once, i=j=1 only), so the Lévy area is trivially zero and the algorithm is correct as-is — confirmed by the weak order ~2 convergence test. The commented-out Lévy area is a pre-existing gap for general m>1 non-diagonal noise; implementing it would be a separate PR. |
|
Does the paper say it should converge order 2 on general non diagonal? |
|
Weak convergence tests should probably be in the special sets for the special runners since they are long running |
|
Yes — Abdulle, Vilmart, Zygalakis (2013) SIAM J. Sci. Comput. 35(4) Theorem 3.4 proves weak order 2 for general non-commutative non-diagonal noise of arbitrary dimension m. The full formula (eq. 17) requires the cross terms J_{q,r} = (ξ_q ξ_r ± χ_q) · h/2 for q ≠ r using Rademacher variables χ. These were previously commented out in both the OOP and IIP SROCK2 paths — now activated. Verified: m=2 non-diagonal gives order ~2.03, m=1 tests still pass. |
ff59b7f to
0f21faf
Compare
|
@ChrisRackauckas can you check once if it still needs changes.suddenly so many tests started failing ? |
|
The big v7 branch merged. The repo won't be in a normal development state until that's all out. Hopefully Friday, trying to get it complete but also it's very delicate and shouldn't be rushed. |
73f4b33 to
8219443
Compare
8219443 to
573a646
Compare
…ML#3170) Non-diagonal problems with a single noise channel were routed into the diagonal code paths by length(W.dW) == 1 shortcuts and crashed with MethodError (oop) or DimensionMismatch (iip). Routing now depends only on is_diagonal_noise and the W.dW type. The general-noise paths drew Rademacher variables through W.rng, which does not exist on NoiseWrapper. They now go through an init_chi! helper that finds the generator on the wrapped source process, and the commented J_qr cross terms from Abdulle, Vilmart and Zygalakis (2013), eq. (17), are activated so SROCK2 keeps weak order 2 for m > 1. The KomBurSROCK2 stage corrections from the original PR SciML#3468 already landed in SciML#3883 and are unchanged here.
ebfadf5 to
056919b
Compare
…ML#3170) Non-diagonal problems with a single noise channel were routed into the diagonal code paths by length(W.dW) == 1 shortcuts and crashed with MethodError (oop) or DimensionMismatch (iip). Routing now depends only on is_diagonal_noise and the W.dW type. The general-noise paths drew Rademacher variables through W.rng, which does not exist on NoiseWrapper. They now go through an init_chi! helper whose accessor walks the source chain of a wrapped process and falls back to the global generator for processes that replay a recorded path and carry none, such as NoiseGrid. KomBurSROCK2 draws from the same helper, so its trajectories for general noise differ from the previous global rand draws. The commented J_qr cross terms of Abdulle, Vilmart and Zygalakis (2013), eq. (17), are activated. They cancel for commutative diffusions, so the new coverage for them is a second moment check on a problem whose two noise channels do not commute. The KomBurSROCK2 stage corrections from the original submission already landed in SciML#3883 and are unchanged here.
056919b to
951f707
Compare
2 participants
Fixes #3188
Fixes #3170
Rebuilt on current master. The KomBurSROCK2 stage math that was part of the original submission landed separately in #3883, so what is left here are the two crashes plus the eq. (17) cross terms.
The first crash is one-channel non-square noise being routed into the diagonal code paths by
length(W.dW) == 1shortcuts. The MWE from #3170:On master that throws MethodError out of place and DimensionMismatch in place for SROCK2, SROCKEM, SROCKC2, KomBurSROCK2 and TangXiaoSROCK2. Dropping the length shortcut from the cache constructors and the perform_step! branches leaves routing to
is_diagonal_noiseand theW.dWtype, and all of those cases pass.The second crash is a non-diagonal problem with a supplied noise process, because the general-noise paths read
W.rngdirectly.NoiseWrapperkeeps the generator on its wrapped source, and that source can itself be a wrapper.NoiseGridreplays a recorded path and has no generator at all. The Rademacher draw now goes through aninit_χ!helper whose accessor walks the source chain and falls back toRandom.default_rng()when there is nothing left to walk to, so a wrapper, a wrapper around a wrapper, and a NoiseGrid all work. SROCK2 with a NoiseGrid crashed on master too and is fixed by the same change.KomBurSROCK2 now takes its Rademacher variables from that helper instead of a bare
rand(), so its trajectories for general noise change. For a NoiseGrid the helper resolves to the global generator and the solve is identical to master value for value; for a problem carrying its own noise process the draw comes from that process. The variables are symmetric either way, so this is a trajectory change only.Finally, the J_{q,r} cross terms of Abdulle, Vilmart and Zygalakis (2013), eq. (17), which were commented out in both SROCK2 general-noise paths, are activated.
On tests, a Core-group file covers m = 1 non-diagonal noise for five algorithms out of place and in place, plus NoiseWrapper, nested NoiseWrapper and NoiseGrid for SROCK2 and KomBurSROCK2. A local-only SROCK2NonDiagonalConvergence group checks weak order 2 against analytic solutions with 5e4 trajectories: 2.02 for m = 1 out of place, 2.02 for m = 1 in place, 2.03 for m = 2 out of place, 2.05 for m = 2 in place.
Those four problems do not cover eq. (17). Every B_k in them is diagonal, so their diffusions commute and the cross terms cancel; reverting the eq. (17) hunks leaves all four order estimates unchanged to 15 digits. Coverage for the cross terms is therefore a separate second moment check on g(u) = [0 1.2*u2; u1 0], whose channels B1 = [0 0; 1 0] and B2 = [0 1.2; 0 0] do not commute. For a linear Ito SDE the second moment solves P' = A P + P A' + sum_k B_k P B_k', so vec(P)(T) = exp(M T) vec(u0 u0') with M = kron(I, A) + kron(A, I) + sum_k kron(B_k, B_k), and here the cross entry is exactly Eu1 u2 = exp(-T). With 2e5 trajectories at dt = 1/2 the measured error is 5.8e-3 with eq. (17) active and about 6e-2 with the terms commented out, so the tolerance is a regression guard for them.
A weak order fit on that problem is not practical: the order-2 bias drops under the Monte Carlo error by dt = 1/4 at any trajectory count that fits the runtime budget, so the check is an accuracy assertion at one step size rather than an order estimate. The whole group runs in about 84 seconds.
StochasticDiffEqROCK is bumped to 2.0.3 and Random is added to its dependencies.
AI Disclosure
Claude assisted with this work.