preconditioning

[filtron]

add preconditoned version of cholesky computation as per feature request.

[nathanaelbosch]

Great, thanks for adding this! I think it might make sense to also add a preconditioned version of the transition matrix A to complete this PR?

[filtron]

Great, thanks for adding this! I think it might make sense to also add a preconditioned version of the transition matrix A to complete this PR?

Hm.... should they both be using the same preconditioner? and also should both functions compute the preconditioner?
If yes to the former that seems like a lot of superfluous computations for the standard use case...
One alternative is to make functions

transition_matrix_and_cov_cholf_1d 
transition_matrix_and_cov_cholf_precond_1d

[nathanaelbosch]

should they both be using the same preconditioner?

This is what the "Stable Implementation ..." paper does at least.

should both functions compute the preconditioner?

Probably not as that would lead to duplication? So maybe either have one function for A, Q, P each, or as you suggest one function that computes both the transition matrix and covariance. I don't know what the best version is; I think that personally in ProbNumDiffEq.jl I would still use my own preconditioner implementation (here) to update them in-place at each step, so for me it would not matter too much which you choose.

[filtron]

So maybe the way to go is to have functions

precond
transition_matrix_precond_1d
cov_cholf_precond_1d

where the latter two only return preconditioned matrices, the correct matrices can then be obtained by multiplying with

precond(ndiff, dt, ReverseTaylor()) |> Diagonal 

[filtron]

@nathanaelbosch

A slight problem with the above of course is that the routines dont know with eltype to use for the arrays.
I therefore changed the call signatures to

transition_matrix_precond_1d(ndiff::Integer, T, ::ReverseTaylor)
transition_cov_cholf_precond_1d(ndiff::Integer, T, ::ReverseTaylor)

where T is the element type desired for the output arrays.
Another alternative is of course to make the user supply arrays to be written into with the desired eltype. Thoughts?

[nathanaelbosch]

Is there a particular reason for why this doesn't output a Diagonal directly? Then all the functiosns return matrices.

[filtron]

We could certainly return a Diagonal matrix instead. But now, weekend.

[filtron]

@nathanaelbosch fine now?

[nathanaelbosch]

Looks good to me! I think the implementation could be made slightly more efficient, at least from what I recall I replaced some parts to compute products cumulatively instead of calling factorial and things like that, but the result is much less readable, so I'd say this is perfectly fine :) Thanks a lot!