Publications

with Führer, C. & Solem, J.-E.

Packt Publishing, 2016. ISBN: 978-1-78646-351-7

Introduction to Python for mathematicians and engineers.

with Cece, E.; Meyrat, P.; Torino, E. & Colarieti-Tosti, M.

Journal of Imaging

doi:10.3390/jimaging9100231

The detection of cancer lesions of a comparable size to that of the typical system resolution of modern scanners is a long-standing problem in Positron Emission Tomography. In this paper, the effect of composing an image-registering convolutional neural network with the modeling of the static data acquisition (i.e., the forward model) is investigated. Two algorithms for Positron Emission Tomography reconstruction with motion and attenuation correction are proposed and their performance is evaluated in the detectability of small pulmonary lesions. The evaluation is performed on synthetic data with respect to chosen figures of merit, visual inspection, and an ideal observer. The commonly used figures of merit (Peak Signal-to-Noise Ratio, Recovery Coefficient, and Signal Difference-to-Noise Ratio) give inconclusive responses, whereas visual inspection and the Channelised Hotelling Observer suggest that the proposed algorithms outperform current clinical practice.

with Sagemüller, J.

arXiv:2209.12459

xAI 2023: Explainable Artificial Intelligence

Vol. 1901 (2023) doi:10.1007/978-3-031-44064-9_19

Various types of saliency methods have been proposed for explaining black-box classification. In image applications, this means highlighting the part of the image that is most relevant for the current decision. We observe that several of these methods can be seen as edge cases of a single, more general procedure based on finding a particular ablation path through the classifier's domain. This gives additional geometric insight to the existing methods. We also demonstrate that this ablation path method can be used as a technique in its own right, the higher computational cost being traded against additional information given by the path.

with Laurent, A.; Munthe-Kaas, H. & McLachlan, R. I.

arXiv:2301.10998

Forum of Mathematics, Sigma

Vol. 11 (2023) doi:10.1017/fms.2023.63

Aromatic B-series were introduced as an extension of standard Butcher series for the study of volume preserving integrators. It was proven with their help that the only volume preserving B-series method is the exact flow of the differential equation. The question was raised whether there exists a volume preserving integrator that can be expanded as an aromatic B-series. In this work, we introduce a new algebraic tool, called the aromatic bicomplex, similar to the variational bicomplex in variational calculus. We prove the exactness of this bicomplex and use it to describe explicitly the key object in the study of volume-preserving integrators: the aromatic forms of vanishing divergence. The analysis provides us with a handful of new tools to study aromatic B-series, gives insights on the process of integration by parts of trees, and allows to describe explicitly the aromatic B-series of a volume-preserving integrator. In particular, we conclude that an aromatic Runge–Kutta method cannot preserve volume.

with Pouchol, C.

arXiv:2006.15845

Inverse Problems and Imaging

Vol. 18, Issue 1 (2024) doi:10.3934/ipi.2023026

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.

with Adler, J.; Lunz, S.; Schönlieb, C.-B. & Öktem, O.

arXiv:1809.00948

Inverse Problems

Vol. 38, Issue 7 (2022) doi:10.1088/1361-6420/ac28ec

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.

with Pouchol, C.

arXiv:2004.02444

Journal of Mathematical Imaging and Vision

Vol. 63, pp. 650–663 (2021) doi:10.1007/s10851-021-01021-7

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.

with Modin, K.

arXiv:1707.04514

Numerische Mathematik

Vol. 145, pp. 405–435, (2020) doi:10.1007/s00211-020-01126-y

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.

with Munthe-Kaas, H. & Stern, A.

arXiv:1903.10056

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.

with Pouchol, C.

arXiv:1909.01966

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.

with Öktem, O. & Pouchol, C.

arXiv:1908.09515

Machine Learning for Medical Image Reconstruction (MLMIR 2019)

pp. 151–162 doi:10.1007/978-3-030-33843-5_14

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.

with Benn, J.; Marsland, S.; McLachlan, R. I. & Modin, K.

arXiv:1702.02780

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.

arXiv:1605.07344

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.

with Bogfjellmo, G. & Modin, K.

arXiv:1703.09589

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.

arXiv:1604.00272

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.

with McLachlan, R. I.; Modin, K. & Munthe-Kaas, H.

arXiv:1512.00906

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.

with Kalisch, H. & Moldabayev, D.

arXiv:1606.01465

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.

with McLachlan, R. I. & Modin, K.

arXiv:1402.3334

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.

with McLachlan, R. I. & Modin, K.

arXiv:1512.04631

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)$.

with McLachlan, R. I. & Modin, K.

arXiv:1505.04035

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.

with Cohen, D.

arXiv:1411.6516

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.

with Xue, H. & Zanna, A.

arXiv:1412.6345

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.

with McLachlan, R. I.; Modin, K. & Munthe-Kaas, H.

arXiv:1409.1019

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.

with Munthe-Kaas, H.

arXiv:1308.5824

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.

with Celledoni, E. & Kometa, B. K.

arXiv:1207.5147

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.

with Munthe-Kaas, H.

arXiv:1402.6981

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.

with McLachlan, R. I. & Modin, K.

arXiv:1307.2387

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.

with McLachlan, R. I. & Modin, K.

arXiv:1402.4114

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.

with McLachlan, R. I. & Modin, K.

arXiv:1308.6620

Nonlinearity

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.

with McLachlan, R. I.; Modin, K. & Wilkins, M.

arXiv:1207.3367

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.

with Modin, K.

arXiv:1212.6618

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.

arXiv:1208.3154

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.

with McLachlan, R. I.; Modin, K. & Wilkins, M.

arXiv:1207.4250

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.

with Fontes, M.

arXiv:1205.0783

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.

with Fontes, M.

arXiv:1205.0780

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.