A list of my books, preprints and articles.

All the preprints and published articles are available on arXiv .

Submitted Manuscripts

Linear inverse problems with nonnegativity constraints through the β-divergences: sparsity of optimisers
with Pouchol, C.


We pass to continuum in optimisation problems associated to linear inverse problems $y = Ax$ where $x$ is constrained to the cone $\{x ≥ 0\}$. We focus on the case where the noise model leads to maximum likelihood estimation through the so-called $β$-divergences, which cover several of the most common noise statistics such as Gaussian, Poisson and multiplicative Gamma. Considering~$x$ as a Radon measure over the domain on which the reconstruction is taking place, we show a general sparsity result. In the high noise regime corresponding to the data not in the image cone, optimisers are typically sparse in the form of sums of Dirac measures. We hence provide an explanation as to why any possible algorithm successfully solving the optimisation problem will lead to undesirably spiky-looking images when the image resolution gets finer, a phenomenon well documented in the literature. We illustrate these results with several numerical examples inspired by medical imaging.

ML-EM algorithm with known continuous movement model
with Pouchol, C.


In Positron Emission Tomography, movement leads to blurry reconstructions when not accounted for. Whether known a priori or estimated jointly to reconstruction, motion models are increasingly defined in continuum rather that in discrete, for example by means of diffeomorphisms. The present work provides both a statistical and functional analytic framework suitable for handling such models. It is based on time-space Poisson point processes as well as regarding images as measures, and allows to compute the maximum likelihood problem for line-of-response data with a known movement model. Solving the resulting optimisation problem, we derive an Maximum Likelihood Expectation Maximisation (ML-EM) type algorithm which recovers the classical ML-EM algorithm as a particular case for a static phantom. The algorithm is proved to be monotone and convergent in the low-noise regime. Simulations confirm that it correctly removes the blur that would have occurred if movement were neglected.

Task adapted reconstruction for inverse problems
with Adler, J.; Lunz, S.; Schönlieb, C.-B. & Öktem, O.


It is notoriously difficult to perform a task defined on a model parameter that is only observed indirectly through noisy data in an ill-posed inverse problem. We describe the steps of reconstruction and task as appropriate estimators (non-randomized decision rules). We use deep neural networks to provide a differentiable parametrization of the family of estimators for both steps. These networks are combined and jointly trained against suitable supervised training data in order to minimize a joint differentiable loss function. This results in an end-to-end task adapted reconstruction method. The suggested framework is generic, yet adaptable, with a plug-and-play structure for adjusting both the inverse problem and the task at hand.

Journal/Conference Publications


What makes nonholonomic integrators work?
with Modin, K.


Numerische Mathematik


A nonholonomic system is a mechanical system with velocity constraints. It has been observed numerically that many integrators exhibit excellent long-time behaviour when applied to some nonholonomic problems. We give a rigorous explanation of this behaviour by introducing a family of nonholonomic problems that are foliated over reversible integrable systems. We give numerical examples on unbiased test problems; all but one of the tested methods fail to produce good long-time behaviour.

Invariant connections, Lie algebra actions, and foundations of numerical integration on manifolds
with Munthe-Kaas, H. & Stern, A.


SIAM Journal on Applied Algebra and Geometry

Vol. 4, Issue 1 (2020) doi:10.1137/19M1252879

Motivated by numerical integration on manifolds, we relate the algebraic properties of invariant connections to their geometric properties. Using this perspective, we generalize some classical results of Cartan and Nomizu to invariant connections on algebroids. This has fundamental consequences for the theory of numerical integrators, giving a characterization of the spaces on which Butcher and Lie–Butcher series methods, which generalize Runge–Kutta methods, may be applied.

The ML-EM algorithm in continuum: sparse measure solutions
with Pouchol, C.


Inverse Problems

Vol. 36, Issue 3 (2020) doi:10.1088/1361-6420/ab6d55

Linear inverse problems $A μ = δ$ with Poisson noise and non-negative unknown $\mu \geq 0$ are ubiquitous in applications, such as Positron Emission Tomography (PET) in medical imaging. The associated maximum likelihood problem is most often solved by the ML-EM algorithm, which gives spiky looking results. This work provides an explanation for this phenomenon. We first consider $μ$ as a measure. We prove that if the measurements $δ$ are not in the cone $\{A \mu, \mu \geq 0\}$, which is typical of short exposure times, likelihood maximisers as well as ML-EM cluster points must be sparse, i.e., typically a sum of point masses. In the long exposure regime, we prove that cluster points of ML-EM will instead be measures without singular part.


Spatiotemporal PET reconstruction using ML-EM with learned diffeomorphic deformation
with Öktem, O. & Pouchol, C.


Machine Learning for Medical Image Reconstruction (MLMIR 2019)


Patient movement in emission tomography deteriorates reconstruction quality because of motion blur. Motion correction algorithms have been implemented to take into account all the gated data, but they do not scale well in computation time, especially not in 3D. We propose a novel motion correction algorithm which addresses the scalability issue. Our approach is to combine an enhanced ML-EM algorithm with deep learning based movement registration. The training is unsupervised, and with artificial data.

Currents and finite elements as tools for shape space
with Benn, J.; Marsland, S.; McLachlan, R. I. & Modin, K.


Journal of Mathematical Imaging and Vision

Vol. 61, Issue 8, pp 1197–1220 (2019) doi:10.1007/s10851-019-00896-x

The nonlinear spaces of shapes (unparameterized immersed curves or submanifolds) are of interest for many applications in image analysis. In this paper we study a general representation of shapes that is based on currents. This representation is robust to noise. We develop the theory of currents for shape spaces by considering both the analytic and numerical aspects of the problem. In particular, we study the analytical properties of the current map and the $H^{-s}$ norm that it induces on shapes. We determine the conditions under which the current determines the shape. We then provide a finite element discretization of the currents that is a practical computational tool for shapes. Finally, we demonstrate this approach on a variety of examples.


Full affine equivariance and weak natural transformations in numerical analysis - the case of B-Series.


Discrete Mechanics, Geometric Integration and Lie–Butcher Series

vol 267 doi:10.1007/978-3-030-01397-4_11

Many algorithms in numerical analysis are affine equivariant: they are immune to changes of affine coordinates. This is because those algorithms are defined using affine invariant constructions. There is, however, a crucial ingredient missing: most algorithms are in fact defined regardless of the underlying dimension. As a result, they are also invariant with respect to non-invertible affine transformation from spaces of different dimensions. We formulate this property precisely: these algorithms fall short of being natural transformations between affine functors. We give a precise definition of what we call a weak natural transformation between functors, and illustrate the point using examples coming from numerical analysis, in particular B-Series.

A Numerical Algorithm for $C^2$-splines on Symmetric Spaces
with Bogfjellmo, G. & Modin, K.


SIAM Journal on Numerical Analysis

Vol. 56, Issue 4 (2018) doi:10.1137/17M1123353

Cubic spline interpolation on Euclidean space has countless applications in science and technology. In several emerging fields (computer vision, quantum control), there is a growing need for spline interpolation on curved spaces. We provide an iterative algorithm for $C^2$-splines on symmetric spaces, and we prove convergence under a verifiable set of conditions. We demonstrate the algorithm for three geometries of interest: $n$-spheres, complex projective spaces, and Grassmannians. It is thereby readily applicable, for example, in computer vision tasks of shape recognition and in quantum control of qubits.

Extension of Matrix Pencil Reduction to Abelian Categories


Journal of Algebra and its Applications

Vol. 16, Number 11 (2018) doi:10.1142/S0219498818500627

Matrix pencils, or pairs of matrices, are used in a variety of applications. By the Kronecker decomposition theorem, they admit a normal form. This normal form consists of four parts, one part based on Jordan decomposition, one part made of nilpotent matrices, and two other dual parts, which we call the observation and control part. The goal of this paper is to show that large parts of that decomposition are still valid for pairs of morphisms in modules or in abelian groups, and more generally in any abelian category.


Butcher Series - A Story of Rooted Trees and Numerical Methods for Evolution Equations
with McLachlan, R. I.; Modin, K. & Munthe-Kaas, H.


Asia Pacific Mathematics Newsletter

Vol. 7, Number 1 (2017) doi:10.5281/zenodo.1297011

Runge–Kutta methods are amongst the most used ODE integrators in engineering. John Butcher discovered a surprising relation with rooted trees. There are further connections to perturbative quantum field theory, and rough paths. The understanding of rooted trees and their connections to ODE integrators recently culminated with the paper of the same authors “B-series methods are exactly the affine equivariant methods”. This paper gives a historical background to rooted trees and numerical methods for evolution equations.

Numerical Study of Nonlinear Dispersive Wave Models with SpecTraVVave
with Kalisch, H. & Moldabayev, D.


Electronic Journal of Differential Equations

Vol. 2017, Number 62, pp. 1–23 (2017) doi:10.5281/zenodo.398855

In nonlinear dispersive evolution equations, the competing effects of nonlinearity and dispersion make a number of interesting phenomena possible. We describe a Python package which computes traveling-wave solutions of equations of the form $u_t + [f(u)]_{x} + \mathcal{L} u_x = 0$, where $\mathcal{L}$ is a self-adjoint operator, and $f$ is a real-valued function with $f(0) = 0$. We use a continuation method coupled with a spectral projection. For the Whitham equation, we obtain numerical evidence that the main bifurcation branch features three distinct points of interest: a turning point, a point of stability inversion, and a terminal point corresponding to a cusped wave. We also found that two solitary waves may interact in such a way that the smaller wave is annihilated. We also observed that bifurcation curves of periodic traveling-wave of the Benjamin solutions may cross and connect high up on the branch in the nonlinear regime.

A minimal-variable symplectic integrator on spheres
with McLachlan, R. I. & Modin, K.


Mathematics of Computation

Vol. 86, No. 307, pp. 2325–2344 (2017) doi:10.1090/mcom/3153

Spin systems model spinning black holes binaries, ferromagnetic materials, or hurricanes on gas planets. Their phase spaces, products of spheres, have a natural symplectic structure for which those systems are Hamiltonian. Is there a method with no extra variable which preserves symplecticity? The answer is yes, and it is a remarkable modification of the midpoint method of symplectic vector spaces. The proof of symplecticity is unexpectedly difficult.


Symmetry reduction for central force problems
with McLachlan, R. I. & Modin, K.


European Journal of Physics

Vol. 37, Number 5 (2016) doi:10.1088/0143-0807/37/5/055003

The Kepler problem, describing the trajectory of the earth around the sun, is a particular instance of a central force problem. Finding explicit solutions of those problems is related to the modern technique of symplectic reduction. We give a detailed mathematical and historical account of that reduction process for the central force problem. It also gives a striking example of a dual pair, namely the pair $\mathfrak{so}(3)-\mathfrak{sl}(2)$.

Geometry of discrete-time spin systems
with McLachlan, R. I. & Modin, K.


Journal of Nonlinear Science

Vol. 25, Issue 5, pp. 1507–1523 (2016) doi:10.1007/s00332-016-9311-z

The spherical midpoint method, developed by the same authors, is related to many other techniques in symplectic and Riemannian geometry. We show that it can be regarded as a midpoint collective integrators on the one hand, and as a Riemannian midpoint method on the other hand.

Multisymplectic discretization of wave map equations
with Cohen, D.


SIAM Journal on Scientific Computing

Vol. 38, Issue 2, pp. A953–A972 (2016) doi:10.1137/15M1014322

How to integrate wave map equations both efficiently and without energy dissipation? We present an explicit integrator for wave maps, which is also a multisymplectic integrator when applied to constrained Hamiltonian PDEs.

A classification of volume preserving generating forms in R³
with Xue, H. & Zanna, A.


Discrete and Continuous Dynamical Systems - Series A

Vol. 36, Issue 4, pp. 2285–2303 (2016) doi:10.3934/dcds.2016.36.2285

Volume preserving integrator are of first importance in statistical physics and weather simulations. One can produce such integrators using generating forms. This generalises the generating function methods for symplectic maps. Remarkably, there turns out to be exactly five distinct cases: (a) a “discrete Lagrangian“ case, (b) a “symplectic Euler” case, (c) a mixture of those two, and (d), (e) two extra mysterious cases which are completely different.

B-series methods are exactly the affine equivariant methods
with McLachlan, R. I.; Modin, K. & Munthe-Kaas, H.


Numerische Matematik

Vol. 133, Issue 3, pp 599–622 (2016) doi:10.1007/s00211-015-0753-2

A class of ODE integrators has been identified as having a particular expansion called B-Series. The question remains: which methods are B-Series? It turns out that they are exactly the methods which respect affine transformation. This gives a concrete and simple way to check that a method is a B-Series.

Aromatic Butcher Series
with Munthe-Kaas, H.


Foundations of Computational Mathematics

Vol. 16, Issue 1, pp 183–215 (2016) doi:10.1007/s10208-015-9245-0

The affine equivariance of an ODE integrator is a very natural property: it means that the integrator is insensitive to affine changes of coordinates. The notion of locality is also natural: it means that the integrator computes fixed points exactly. What are the local and affine equivariant methods? We show here that an integrator is local and affine equivariant if and only if it has an expansion in aromatic Butcher series which are a generalisation of Butcher Series.

High order semi-Lagrangian methods for the incompressible Navier–Stokes equations
with Celledoni, E. & Kometa, B. K.


Journal of Scientific Computing

Vol. 66, Issue 1, pp 91–115 (2016) doi:10.1007/s10915-015-0015-6

We seek high order time integrators for the Navier–Stokes equations. The Navier–Stokes equations are a differential equation with a constraints (DAE). Standard DAE Runge–Kutta (RK) methods are known to have maximal order much lower than their RK counterpart. The idea is to use a standard, high order RK method with a suitable projection, in order to obtain high order integrators for Navier–Stokes.

Integrators on homogeneous spaces: Isotropy choice and connections
with Munthe-Kaas, H.


Foundations of Computational Mathematics

Vol. 16, Issue 4, pp 899–939 (2016) doi:10.1007/s10208-015-9267-7

Runge–Kutta methods are typically defined on affine spaces. In that case they are known to be affine equivariant: they are immune to affine changes of coordinates—a fundamental property in practice. What about Runge–Kutta methods on other homogeneous spaces? It turns out that this depends on the choice of isotropy, which is a way to relate the original vector field with an infinitesimal displacements. We show that if such a choice is equivariant, then the various RK methods also are equivariant. We also show that the existence of such equivariant choices of isotropy is related to the property of reductivity of the underlying homogeneous space.


Collective Lie–Poisson integrators on R³
with McLachlan, R. I. & Modin, K.


IMA Journal of Numerical Analysis

Vol. 35, Issue 2, pp. 546-560 (2015) doi:10.1093/imanum/dru013

We describe a Lie–Poisson integrator constructed using the Hopf momentum map, and the related dual pair $\mathfrak{so}(3)$–$\mathfrak{u}(1)$. We give several proofs that the method is symplectic, and preserves the coadjoint orbits.


Symplectic integrators for spin systems
with McLachlan, R. I. & Modin, K.


Physical Review E

Vol. 89, Issue 6 (2014) doi:10.1103/PhysRevE.89.061301

Overview of a new generating function for phase spaces which are products of spheres. In particular, this yields new integrable discretizations of the reduced motion of a free rigid body.

Collective symplectic integrators
with McLachlan, R. I. & Modin, K.



Vol. 27, Issue 6, p. 1525 (2014) doi:10.1088/0951-7715/27/6/1525

We describe a remarkably general recipe to construct symplectic integrators on Lie–Poisson spaces. These are Lie–Poisson integrators which preserve symplectic leaves. The idea is to consider the Hamiltonian action on a symplectic vector space of the Lie group $G$ associated to the Lie–Poisson space. The corresponding momentum map is a Poisson map. If the orbits of $G$ are quadratic, all symplectic Runge–Kutta methods descend to a symplectic integrator on the Lie–Poisson space. We investigate various ways to weaken that assumption.

Geometric Generalisations of Shake and Rattle
with McLachlan, R. I.; Modin, K. & Wilkins, M.


Foundations of Computational Mathematics

Vol. 14, Issue 2, pp 339–370 (2014) doi:10.1007/s10208-013-9163-y

The Shake and Rattle algorithms are widely used to integrate mechanical systems with holonomic constraints. In this paper we show exactly when and how these methods work. In particular, we show that they work no matter the choice of the auxiliary symplectic solver in the ambient space. We also show that the methods work if and only if the constraint manifold is a coisotropic submanifold.

Integrability of Nonholonomically Coupled Oscillators
with Modin, K.


Discrete and Continuous Dynamical Systems - Series A

Vol. 34, Issue 3, pp 1121–1130 (2014) doi:10.3934/dcds.2014.34.1121

We study the non-holonomic particle system. It can be seen as a model for a continuously varying transmission. We show that this system is reversible integrable. This explains the excellent properties of any reversible integrators, in particular of some purported non-holonomic integrators studied in the literature.

Reductions of Operator Pencils


Mathematics of Computation

Vol. 83, pp 189-214 (2014) doi:10.1090/S0025-5718-2013-02740-8

The Stokes equation, a linear model of fluid dynamics, is but one example of a linear PDE with constraints. Other examples are saddle point problems and mixed formulation. In finite dimensions, such systems are called (linear) pencils. We show that there can be no notion of index for such linear pencils in infinite dimensions. We also show the close relationship between the celebrated inf-sup condition and the property to be a regular pencil.


Symplectic integrators for index one constraints
with McLachlan, R. I.; Modin, K. & Wilkins, M.


SIAM Journal on Scientific Computing

Vol. 35, Issue 5, pp 2150–2162 (2013) doi:10.1137/120885085

Problems arising in control theory can be of sub-Riemannian type. One has to integrate a differential equation with contraints on the velocity. We show that symplectic Runge–Kutta method can be good integrators for such problems.


Simplified apriori Estimate for the Time Periodic Burgers Equation
with Fontes, M.


Proceedings of the Estonian Academy of Sciences

Vol. 59, Issue 1, pp 34–41 (2010) doi:10.3176/proc.2010.1.06

A simplified proof of the existence and uniqueness of time-periodic solutions of the Burgers equations for data with regularity -1 in space.


Time-Periodic Solutions of the Burgers Equation
with Fontes, M.


Journal of Mathematical Fluid Mechanics

Vol. 11, Issue 2, pp 303–323 (2009) doi:10.1007/s00021-007-0260-z

We show that the viscous, nonlinear, time periodic Burgers operator is a diffeomorphism between a function space with fractional regularity and its dual.


book cover
with Führer, C. & Solem, J.-E.
Packt Publishing, 2016. ISBN: 978-1-78646-351-7

Introduction to Python for mathematicians and engineers.