A list of my books, preprints and articles.

All the preprints and published articles are available on arXiv .

A list of my books, preprints and articles.

All the preprints and published articles are available on arXiv .

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.

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.

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.

`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.

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.

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.

`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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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*.

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.

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*.

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.

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.

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.

Introduction to Python for mathematicians and engineers.