Mega high order integrated Wiener process priors

[pnkraemer]

EDIT: we probably want https://link.springer.com/article/10.1007/BF03167904. Everything written below is outdated.

_This issue is a reminder how to get to crazy high order IBM priors (@nathanaelbosch).

The bottleneck (without damping) is the Cholesky decomposition of the process noise covariance: maximum is ordint=10. With a little damping (see #404), the maximum becomes ordint=15for smoothing.

It is possible to remove the damping entirely by means of high-precision cholesky decompositions (for instance, these could be precomputed for the preconditioned, i.e. step=1 transition) using https://mpmath.org/doc/current/matrices.html#linear-algebra.

The following snippet produces those factors:

from mpmath import *
mp.dps = 75   # gets to ordint=53; higher values give higher ordints; what does dps mean???
for i in range(100):
    try:
        cholesky(hilbert(i))
    except ValueError:  # not posdef anymore
        print(f"Problem occurs at i={i}")
        break

While calling this high-precision arithmetic at any solve sounds crazy to me, one could precompute those matrices, store them as .npy or whatever, put them in statespace (the files should be tiny).

Of course, fixing this Cholesky factor does not enable filtering and smoothing of arbitrarily high order: at some point (order 15? 20?) the smoothing iteration overflows. However, the current bottleneck would be resolved by this trick here :)

See also: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.hilbert.html and https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.invhilbert.html, which may be interesting to us._

[pnkraemer]

While writing this, I stumbled upon https://link.springer.com/article/10.1007/BF03167904 which might be an easier, and even better solution :)

[pnkraemer]

and https://people.eecs.berkeley.edu/~wkahan/MathH110/HilbMats.pdf

[JonathanWenger]

@pnkraemer If this issue is outdated can it be closed or the content updated appropriately? Alternatively we can move this to discussions under research.