Initialize without discretizing to save some matrix exponentials

nathanaelbosch · May 14, 2023 · 3556779 · 3556779
1 parent a904f54
commit 3556779
Show file tree

Hide file tree

Showing 2 changed files with 4 additions and 3 deletions.
diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl
@@ -5,11 +5,10 @@ function calc_H!(H, integ, cache)
     if integ.alg isa EK0
         calc_H_EK0!(H, integ, cache)
         if integ.alg.prior isa IOUP
-            H .-= cache.prior.rate_parameter * E0
+            ProbNumDiffEq._matmul!(H, cache.prior.rate_parameter, E0, -1.0, 1.0)
         end
     elseif integ.alg isa EK1
         calc_H_EK0!(H, integ, cache)
-        # @assert integ.u == @view x_pred.μ[1:(q+1):end]
         if integ.f isa SplitFunction
             ddu .= integ.f.f1.f
         else

diff --git a/src/priors/common.jl b/src/priors/common.jl
@@ -32,7 +32,9 @@ See also: [`make_transition_matrices`](@ref).
 [1] N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
 """
 function initialize_transition_matrices(p::AbstractODEFilterPrior{T}, dt) where {T}
-    Ah, Qh = discretize(p, dt)
+    d, q = p.wiener_process_dimension, p.num_derivatives
+    D = d * (q+1)
+    Ah, Qh = zeros(T, D, D), PSDMatrix(zeros(T, D, D))
     P, PI = initialize_preconditioner(p, dt)
     A = P * Ah * PI
     Q = X_A_Xt(Qh, P)