Analytically compute the square root of the transition covariance

src/priors/iwp.jl

@@ -44,13 +44,37 @@ function to_1d_sde(p::IWP)
     return LTISDE(F_breve, L_breve)
 end
 
+function make_preconditioned_transition_cov_cholf_1d(q::Integer, ::Type{elType}) where {elType}
+    if q >= 10 && !(q isa BigInt)
+        return make_preconditioned_transition_cov_cholf_1d(big(q), elType)
+    end
+
+    L = if q isa BigInt
+        zeros(elType, q + 1, q + 1)
+    else
+        Array{elType}(calloc, q + 1, q + 1)
+    end
+
+    @simd ivdep for n in 0:q
+        @simd ivdep for m in 0:q
+            if m <= n
+                @inbounds L[q+1-m, q+1-n] =
+                    sqrt(2 * m + 1) * factorial(n)^2 / factorial(n - m) /
+                    factorial(n + m + 1)
+            end
+        end
+    end
+    return L
+end
+
 function preconditioned_discretize_1d(iwp::IWP{elType}) where {elType}
     q = iwp.num_derivatives
 
     A_breve = binomial.(q:-1:0, (q:-1:0)')
-    Q_breve = Cauchy(collect(q:-1.0:0.0), collect((q+1):-1.0:1.0)) |> Matrix  # for Julia1.6
+    # Q_breve = Cauchy(collect(q:-1.0:0.0), collect((q+1):-1.0:1.0)) |> Matrix  # for Julia1.6
 
-    QR_breve = cholesky(Q_breve).L' |> collect
+    # QR_breve = cholesky(Q_breve).L' |> collect
+    QR_breve = make_preconditioned_transition_cov_cholf_1d(q, elType)
     A_breve, QR_breve = elType.(A_breve), elType.(QR_breve)
     Q_breve = PSDMatrix(QR_breve)