From 7c3db97366ead019adaa8979d97efdeb02797ceb Mon Sep 17 00:00:00 2001 From: Filip Tronarp Date: Thu, 12 Oct 2023 16:43:01 +0200 Subject: [PATCH 1/3] preconditioning --- src/iwp.jl | 15 +++++++++++++++ test/runtests.jl | 8 ++++++++ 2 files changed, 23 insertions(+) diff --git a/src/iwp.jl b/src/iwp.jl index ac3c567..817caa9 100644 --- a/src/iwp.jl +++ b/src/iwp.jl @@ -110,6 +110,21 @@ function transition_cov_cholf_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where { return L end + +function transition_cov_cholf_precond_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T} + L = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff)) + precond = @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff) + @simd ivdep for n = 0:ndiff + @simd ivdep for m = 0:ndiff + if m <= n + @inbounds L[n+1, m+1] = + sqrt(2 * m + 1) * factorial(n)^2 / factorial(n - m) / factorial(n + m + 1) + end + end + end + return precond, L +end + """ transition_cov(model::IWP{T}, dt::T, method::AbstractRealizationMethod) where {T} diff --git a/test/runtests.jl b/test/runtests.jl index e55a898..cf7e475 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,5 +31,13 @@ using Test end + @testset "preconditioning" begin + for ndiff = 0:9 + L = IntegratedWienerProcesses.transition_cov_cholf_1d(ndiff, dt, ReverseTaylor()) + precond, Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, dt, ReverseTaylor()) + @test L ≈ Diagonal(precond) * Lbreve + end + end + end From 5ff2c046044356e997b428bf487708d3257d6348 Mon Sep 17 00:00:00 2001 From: Filip Tronarp Date: Fri, 13 Oct 2023 16:23:56 +0200 Subject: [PATCH 2/3] preconditioned transition matrix --- src/iwp.jl | 23 ++++++++++++++++++----- test/runtests.jl | 9 +++++++-- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/src/iwp.jl b/src/iwp.jl index 817caa9..59d9d39 100644 --- a/src/iwp.jl +++ b/src/iwp.jl @@ -64,10 +64,10 @@ function transition_matrix( method::AbstractRealizationMethod, ) where {T} V = promote_type(real(T), typeof(dt)) - return kron(_transition_matrix_1d(ndiff(model), V(dt), method), one(T)I(dim(model))) + return kron(transition_matrix_1d(ndiff(model), V(dt), method), one(T)I(dim(model))) end -function _transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<:Real} +function transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<:Real} irange = collect(0:ndiff) #g = @. δt^irange / factorial(irange) g = @. exp(irange * log(dt) - logfactorial(irange)) @@ -82,6 +82,20 @@ function _transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T< end +precond(ndiff::Integer, dt::Real, ::ReverseTaylor) = @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff) + +function transition_matrix_precond_1d(ndiff::Integer, T, ::ReverseTaylor) + Φ = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff)) + @simd ivdep for col = 0:ndiff + @simd ivdep for row = col:ndiff + @inbounds Φ[row+1, col+1] = binomial(row, col) + end + end + return Φ +end + + + """ transition_cov_cholf(model::IWP{T}, dt::T, method::AbstractRealizationMethod) where {T} @@ -111,9 +125,8 @@ function transition_cov_cholf_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where { end -function transition_cov_cholf_precond_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T} +function transition_cov_cholf_precond_1d(ndiff::Integer, T, ::ReverseTaylor) L = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff)) - precond = @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff) @simd ivdep for n = 0:ndiff @simd ivdep for m = 0:ndiff if m <= n @@ -122,7 +135,7 @@ function transition_cov_cholf_precond_1d(ndiff::Integer, dt::T, ::ReverseTaylor) end end end - return precond, L + return L end """ diff --git a/test/runtests.jl b/test/runtests.jl index cf7e475..b19d3f3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,9 +33,14 @@ using Test @testset "preconditioning" begin for ndiff = 0:9 + Φ = IntegratedWienerProcesses.transition_matrix_1d(ndiff, dt, ReverseTaylor()) L = IntegratedWienerProcesses.transition_cov_cholf_1d(ndiff, dt, ReverseTaylor()) - precond, Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, dt, ReverseTaylor()) - @test L ≈ Diagonal(precond) * Lbreve + prec = IntegratedWienerProcesses.precond(ndiff, dt, ReverseTaylor()) + Φbreve = IntegratedWienerProcesses.transition_matrix_precond_1d(ndiff, Float64, ReverseTaylor()) + Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, Float64, ReverseTaylor()) + + @test Φ ≈ Diagonal(prec) * Φbreve / Diagonal(prec) + @test L ≈ Diagonal(prec) * Lbreve end end From 3b3354263f80ca7c081620ed279d9099d5321da9 Mon Sep 17 00:00:00 2001 From: Filip Tronarp Date: Wed, 18 Oct 2023 16:10:03 +0200 Subject: [PATCH 3/3] wrap preconditioner in Diagonal --- src/iwp.jl | 2 +- test/runtests.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/iwp.jl b/src/iwp.jl index 59d9d39..4786291 100644 --- a/src/iwp.jl +++ b/src/iwp.jl @@ -82,7 +82,7 @@ function transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<: end -precond(ndiff::Integer, dt::Real, ::ReverseTaylor) = @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff) +precond(ndiff::Integer, dt::Real, ::ReverseTaylor) = Diagonal(@. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff)) function transition_matrix_precond_1d(ndiff::Integer, T, ::ReverseTaylor) Φ = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff)) diff --git a/test/runtests.jl b/test/runtests.jl index b19d3f3..7527e31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,8 +39,8 @@ using Test Φbreve = IntegratedWienerProcesses.transition_matrix_precond_1d(ndiff, Float64, ReverseTaylor()) Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, Float64, ReverseTaylor()) - @test Φ ≈ Diagonal(prec) * Φbreve / Diagonal(prec) - @test L ≈ Diagonal(prec) * Lbreve + @test Φ ≈ prec * Φbreve / prec + @test L ≈ prec * Lbreve end end