Fix issues with Cholesky of ForwardDiff.Dual matrices

Author: nathanaelbosch

src/ProbNumDiffEq.jl

@@ -56,6 +56,14 @@ cov2psdmatrix(cov::Diagonal; d) =
 cov2psdmatrix(cov::AbstractMatrix; d) =
     (@assert size(cov, 1) == size(cov, 2) == d; PSDMatrix(Matrix(cholesky(cov).U)))
 cov2psdmatrix(cov::PSDMatrix; d) = (@assert size(cov, 1) == size(cov, 2) == d; cov)
+make_hermitian_if_fowarddiff(M::AbstractMatrix) = M
+make_hermitian_if_fowarddiff(M::AbstractMatrix{<:ForwardDiff.Dual}) = begin
+    # Since ForwardDiff 1.0 `ishermitian` fails for small fluctuations in the dual
+    if !ishermitian(ForwardDiff.value.(M))
+        error("Matrix is not Hermitian.")
+    end
+    return 0.5 * (M + M')
+end
 
 """
     add!(out, toadd)

src/diffusions/calibration.jl

@@ -8,6 +8,7 @@ Needed for MLE diffusion estimation.
 invquad
 invquad(v, M::Matrix; v_cache, M_cache) = begin
     v_cache .= v
+    M = make_hermitian_if_fowarddiff(M)
     M_chol = cholesky!(copy!(M_cache, M))
     ldiv!(M_chol, v_cache)
     dot(v, v_cache)

src/filtering/update.jl

@@ -104,6 +104,7 @@ function update!(
         _matmul!(K2_cache, P_p, H')
     end
 
+    _S = make_hermitian_if_fowarddiff(_S)
     S_chol = length(_S) == 1 ? _S[1] : cholesky!(_S)
     rdiv!(K, S_chol)
 

src/solve.jl

@@ -2,12 +2,10 @@
 function DiffEqBase.__init(
     prob::DiffEqBase.AbstractODEProblem{uType,tType,false},
     alg::AbstractEK,
-    timeseries_init=(),
-    ts_init=(),
-    ks_init=(),
-    recompile::Type{Val{recompile_flag}}=Val{true};
+    args...;
     kwargs...,
-) where {uType,tType,recompile_flag}
+) where {uType,tType}
+    # @info "Inside ProbNumDiffEq's overloaded DiffEqBase.__init function"
     @warn "The given problem is in out-of-place form. Since the algorithms in this " *
           "package are written for in-place problems, it will be automatically converted."
     if prob.f isa DynamicalODEFunction
@@ -43,13 +41,12 @@ function DiffEqBase.__init(
             prob.kwargs...,
         )
     end
+
+    # @info "Calling DiffEqBase.__init now"
     return DiffEqBase.__init(
         _prob,
         alg,
-        timeseries_init,
-        ts_init,
-        ks_init,
-        recompile;
+        args...;
         kwargs...,
     )
 end