From 1616635ea29cbc0b48ac40fe404daabf8b848910 Mon Sep 17 00:00:00 2001 From: Nathanael Bosch Date: Mon, 23 Aug 2021 15:21:20 +0200 Subject: [PATCH] Better preconditioning (#53) * Better smoothing stability * CompatHelper: add new compat entry for "Tullio" at version "0.3" (#51) Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * CompatHelper: add new compat entry for "Octavian" at version "0.3" (#50) * CompatHelper: add new compat entry for "Octavian" at version "0.3" * Update Project.toml Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Nathanael Bosch * CompatHelper: bump compat for "ModelingToolkit" to "6" (#46) * CompatHelper: bump compat for "ModelingToolkit" to "6" * Update Project.toml Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Nathanael Bosch * Add a function to return the empty preconditioning matrices * Changed the preconditioning and the priors to use SotA * Fix a test and improve code * Updated the convergence tests * Bump version Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- Project.toml | 2 +- src/caches.jl | 3 +-- src/preconditioning.jl | 37 +++++++++++++++++++++++++++++++++++-- src/priors.jl | 40 ++++++++++++++++++++++++++++++++++++++-- src/solution.jl | 3 +-- test/convergence.jl | 31 +++++++++++++------------------ test/preconditioning.jl | 3 +-- test/priors.jl | 8 ++++---- 8 files changed, 94 insertions(+), 33 deletions(-) diff --git a/Project.toml b/Project.toml index 3e632c322..8f8761bc7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ProbNumDiffEq" uuid = "bf3e78b0-7d74-48a5-b855-9609533b56a5" authors = ["Nathanael Bosch"] -version = "0.1.7" +version = "0.1.8" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/src/caches.jl b/src/caches.jl index 271291c92..6b87e3a2b 100644 --- a/src/caches.jl +++ b/src/caches.jl @@ -87,8 +87,7 @@ function OrdinaryDiffEq.alg_cache( # Prior dynamics @assert alg.prior == :ibm "Only the ibm prior is implemented so far" - P = Diagonal(ones(uElType, D)) - PI = Diagonal(ones(uElType, D)) + P, PI = init_preconditioner(d, q, uElType) A, Q = ibm(d, q, uElType) diff --git a/src/preconditioning.jl b/src/preconditioning.jl index e96d682de..2a017dd87 100644 --- a/src/preconditioning.jl +++ b/src/preconditioning.jl @@ -1,4 +1,11 @@ -@fastmath @inbounds function make_preconditioner!(P, h, d, q) +function init_preconditioner(d, q, elType=typeof(1.0)) + D = d*(q+1) + P = Diagonal(ones(elType, D)) + PI = Diagonal(ones(elType, D)) + return P, PI +end + +@fastmath @inbounds function make_preconditioner_old!(P, h, d, q) val = h^(-q-1/2) @simd for j in 0:q @simd for i in 1:d @@ -10,7 +17,7 @@ end -@fastmath @inbounds function make_preconditioner_inv!(PI, h, d, q) +@fastmath @inbounds function make_preconditioner_old_inv!(PI, h, d, q) val = h^(q+1/2) @simd for j in 0:q @simd for i in 1:d @@ -32,3 +39,29 @@ function make_preconditioners!(post::GaussianODEFilterPosterior, dt) make_preconditioner!(P, dt, d, q) make_preconditioner_inv!(PI, dt, d, q) end + + + + +@fastmath @inbounds function make_preconditioner!(P, h, d, q) + val = factorial(q) / h^(q+1/2) + @simd for j in 0:q + @simd for i in 1:d + P[j*d + i,j*d + i] = val + end + val /= (q-j) / h + end + return P +end + + +@fastmath @inbounds function make_preconditioner_inv!(PI, h, d, q) + val = h^(q+1/2) / factorial(q) + @simd for j in 0:q + @simd for i in 1:d + PI[j*d + i,j*d + i] = val + end + val *= (q-j) / h + end + return PI +end diff --git a/src/priors.jl b/src/priors.jl index 354acf773..3525e3bf0 100644 --- a/src/priors.jl +++ b/src/priors.jl @@ -1,10 +1,13 @@ ######################################################################################## # Integrated Brownian Motion ######################################################################################## -"""Generate the discrete dynamics for a q-IBM model. INCLUDES AUTOMATIC PRECONDITIONING! +""" + ibm(d::Integer, q::Integer, elType=typeof(1.0)) + +Generate the discrete dynamics for a q-IBM model. INCLUDES AUTOMATIC PRECONDITIONING! Careful: Dimensions are ordered differently than in `probnum`!""" -function ibm(d::Integer, q::Integer, elType=typeof(1.0)) +function ibm_old(d::Integer, q::Integer, elType=typeof(1.0)) A_base = diagm(0=>ones(elType, d*(q+1))) Q_base = zeros(elType, d*(q+1), d*(q+1)) @@ -56,6 +59,39 @@ function ibm(d::Integer, q::Integer, elType=typeof(1.0)) end +function ibm(d::Integer, q::Integer, elType=typeof(1.0)) + D = d*(q+1) + + # Make A + A_breve = zeros(elType, q+1, q+1) + @simd ivdep for j in 1:q+1 + @simd ivdep for i in 1:j + @inbounds A_breve[i, j] = binomial(q-i+1, q-j+1) + end + end + A = kron(A_breve, I(d)) + @assert istriu(A) + A = UpperTriangular(A) + + # Make Q + Q_breve = zeros(elType, q+1, q+1) + @fastmath _transdiff_ibm_element(row::Int, col::Int) = + one(elType) / (2 * q + 1 - row - col) + @simd for col in 0:q + @simd for row in col:q + val = _transdiff_ibm_element(row, col) + @inbounds Q_breve[1 + col, 1 + row] = val + @inbounds Q_breve[1 + row, 1 + col] = val + end + end + QL_breve = cholesky(Q_breve).L + QL = kron(QL_breve, I(d)) + Q = SRMatrix(QL) + + return A, Q +end + + """Same as above, but without the automatic preconditioning""" function vanilla_ibm(d::Integer, q::Integer) diff --git a/src/solution.jl b/src/solution.jl index 4efade096..53ab4f11a 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -161,8 +161,7 @@ function GaussianODEFilterPosterior(alg, u0) SolProj = u0 isa ArrayPartition ? [Proj(1); Proj(0)] : Proj(0) A, Q = ibm(d, q, uElType) - P = Diagonal(ones(uElType, D)) - PI = Diagonal(ones(uElType, D)) + P, PI = init_preconditioner(d, q, uElType) GaussianODEFilterPosterior( d, q, SolProj, A, Q, P, PI, false) end diff --git a/test/convergence.jl b/test/convergence.jl index 0f6ea3a54..3ac9b5a18 100644 --- a/test/convergence.jl +++ b/test/convergence.jl @@ -9,31 +9,26 @@ TESTTOL = 0.2 # Simple linear problem prob = ODEProblem( ODEFunction( - (u, p, t) -> [1.01].*u, - jac=(u, p, t) -> [1.01]; - analytic=(u0, p, t) -> u0.*[exp(1.01*t)], - ), [big(1/2)], (0.0, 1.0)) + (u, p, t) -> 1.01 .* u, + jac=(u, p, t) -> 1.01; + analytic=(u0, p, t) -> u0 .* exp(1.01*t), + ), [big(1/2)], big.((0.0, 1.0))) -@testset "EK0(order=$q)" for q in 1:3 - dts = 1 .//2 .^(9:-1:2) - sim = test_convergence(dts, prob, EK0(order=q)) - @test sim.𝒪est[:final] ≈ q+1 atol=TESTTOL - @test sim.𝒪est[:l2] ≈ q+1 atol=TESTTOL - @test sim.𝒪est[:l∞] ≈ q+1 atol=TESTTOL -end -@testset "EK0(order=$q)" for q in 4:5 +@testset "EK0(order=$q) convergence" for q in 1:5 dts = 1 .//2 .^(8:-1:4) sim = test_convergence(dts, prob, EK0(order=q)) - @test sim.𝒪est[:final] ≈ q+1 atol=TESTTOL - @test sim.𝒪est[:l2] ≈ q+1 atol=TESTTOL+0.1 - @test sim.𝒪est[:l∞] ≈ q+1 atol=TESTTOL+0.1 + @test sim.𝒪est[:final] ≈ q+1 atol=0.2 + @test sim.𝒪est[:l2] ≈ q+1 atol=0.25 + @test sim.𝒪est[:l∞] ≈ q+1 atol=0.2 end -@testset "EK1(order=$q)" for q in [1, 3, 4, 5] - dts = 1 .//2 .^(8:-1:3) +@testset "EK1(order=$q) convergence" for q in 1:5 + dts = 1 .//2 .^(8:-1:4) sim = test_convergence(dts, prob, EK1(order=q)) - @test sim.𝒪est[:l2] ≈ q+1 atol=TESTTOL+0.1 + @test sim.𝒪est[:final] ≈ q+1 atol=0.55 + @test sim.𝒪est[:l2] ≈ q+1 atol=0.5 + @test sim.𝒪est[:l∞] ≈ q+1 atol=0.5 end diff --git a/test/preconditioning.jl b/test/preconditioning.jl index f867f262c..0b5bfe845 100644 --- a/test/preconditioning.jl +++ b/test/preconditioning.jl @@ -28,8 +28,7 @@ prob = ProbNumDiffEq.remake_prob_with_jac(prob) # First test that they're both equivalent D = d*(q+1) - P = Diagonal(ones(D)) - PI = Diagonal(ones(D)) + P, PI = ProbNumDiffEq.init_preconditioner(d, q) ProbNumDiffEq.make_preconditioner!(P, h, d, q) ProbNumDiffEq.make_preconditioner_inv!(PI, h, d, q) @test Qh_p ≈ P * Qh * P' diff --git a/test/priors.jl b/test/priors.jl index 758854a79..8e478e52f 100644 --- a/test/priors.jl +++ b/test/priors.jl @@ -47,13 +47,13 @@ end A, Q = ProbNumDiffEq.ibm(d, q) Qh = Q * σ^2 - AH_21_PRE = [1 1 0.5 + AH_21_PRE = [1 2 1 0 1 1 0 0 1] - QH_21_PRE = σ^2 * [1/20 1/8 1/6 - 1/8 1/3 1/2 - 1/6 1/2 1] + QH_21_PRE = σ^2 * [1/5 1/4 1/3 + 1/4 1/3 1/2 + 1/3 1/2 1] @test AH_21_PRE ≈ A @test QH_21_PRE ≈ Qh