Skip to content

Commit

Permalink
Better preconditioning (#53)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>

* 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 <[email protected]>

* 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>
  • Loading branch information
@nathanaelbosch @github-actions
nathanaelbosch and github-actions[bot] authored Aug 23, 2021
1 parent 5b836ee commit 1616635
Show file tree
Hide file tree
Showing 8 changed files with 94 additions and 33 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
3 changes: 1 addition & 2 deletions src/caches.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
37 changes: 35 additions & 2 deletions src/preconditioning.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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
40 changes: 38 additions & 2 deletions src/priors.jl
View file
Original file line number Diff line number Diff line change
@@ -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))
Expand Down Expand Up @@ -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)

Expand Down
3 changes: 1 addition & 2 deletions src/solution.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
31 changes: 13 additions & 18 deletions test/convergence.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
3 changes: 1 addition & 2 deletions test/preconditioning.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down
8 changes: 4 additions & 4 deletions test/priors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

2 comments on commit 1616635

@nathanaelbosch
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/43376

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.8 -m "<description of version>" 1616635ea29cbc0b48ac40fe404daabf8b848910
git push origin v0.1.8

Please sign in to comment.