Merge pull request #166 from pnkraemer/new_precond

Preconditioner

Author: pnkraemer

src/probnum/diffeq/odefiltsmooth/prior.py

@@ -41,7 +41,7 @@ class ODEPrior(LTISDEModel):
     By default, we choose :math:`P` to be the
     matrix that maps to filtering iteration to the Nordsieck vector,
 
-    .. math:: P = \\text{diag }(1, h, h^2, ..., h^q).
+    .. math:: P = \\text{diag }(h^{-q}, h^{-q+1}, ..., 1).
 
     Here, :math:`h` is some expected average step size. Note that we
     ignored the factorials in this matrix. Our setting makes it easy to
@@ -70,6 +70,10 @@ class ODEPrior(LTISDEModel):
     steps and higher order methods which especially makes smoothing
     algorithms unstable.
 
+    Another advantage of this preconditioning is that the smallest value
+    that appears inside the algorithm is :math:`h^{q}` (with
+    preconditioning) instead of :math:`h^{2q+1}` (without preconditioning).
+
     The matrices :math:`F, u, L` are the usual matrices for
     IBM(:math:`q`), IOUP(:math:`q`) or Matern(:math:`q+1/2`) processes.
     As always, :math:`B(t)` is
@@ -149,7 +153,7 @@ def precond2nordsieck(self, step):
             )
             warnings.warn(message=warnmsg, category=RuntimeWarning)
             step = 1e-15 ** (1 / self.ordint)
-        powers = np.arange(self.ordint + 1)
+        powers = np.arange(start=-self.ordint, stop=1)
         diags = step ** powers
         precond = np.kron(np.eye(self.spatialdim), np.diag(diags))
         invprecond = np.kron(np.eye(self.spatialdim), np.diag(1.0 / diags))

tests/test_diffeq/test_odefiltsmooth/test_odefiltsmooth.py

@@ -469,3 +469,23 @@ def test_filter_ivp_h_mat52_ukf(self):
 
     def test_filter_ivp_h_mat72_kf(self):
         probsolve_ivp(self.ivp, step=self.step, which_prior="matern72", method="eks0")
+
+
+class TestPreconditioning(unittest.TestCase):
+    """
+    Solver with high order and small stepsize should work up to a point where
+    step**order is below machine precision.
+    """
+
+    def setUp(self):
+        initdist = RandomVariable(distribution=Dirac(20 * np.ones(2)))
+        self.ivp = ode.lotkavolterra([0.0, 1e-4], initdist)
+        self.step = 1e-5
+        self.prior = "ibm3"
+
+    def test_small_step_feasible(self):
+        """
+        With the 'old' preconditioner, this is impossible because step**(2*order + 1) is too small.
+        With the 'new' preconditioner, the smallest value that appears in the solver code is step**order
+        """
+        probsolve_ivp(self.ivp, step=self.step, which_prior=self.prior, method="eks0")

tests/test_diffeq/test_odefiltsmooth/test_prior.py

@@ -46,7 +46,7 @@
 
 QH_21_PRE = (
     DIFFCONST ** 2
-    * STEP ** 5
+    * STEP
     * np.array([[1 / 20, 1 / 8, 1 / 6], [1 / 8, 1 / 3, 1 / 2], [1 / 6, 1 / 2, 1]])
 )