Consistent Inversion
of Noisy Non-Abelian X-Ray Transforms
FRANÇOIS MONARD
University of California, Santa Cruz
RICHARD NICKL
University of Cambridge
AND
GABRIEL P. PATERNAIN
University of Cambridge
Abstract
For M a simple surface, the nonlinear statistical inverse problem of recovering a
matrix field  W M ! so.n/ from discrete, noisy measurements of the SO.n/-
valued scattering data C of a solution of a matrix ODE is considered (n  2).
Injectivity of the map 7! C was established by Paternain, Salo, and Uhlmann
in 2012.
A statistical algorithm for the solution of this inverse problem based on Gauss-
ian process priors is proposed, and it is shown how it can be implemented by
infinite-dimensional MCMC methods. It is further shown that as the number
N of measurements of point evaluations of C increases, the statistical error in
the recovery of  converges to 0 in L2.M/-distance at a rate that is algebraic
in 1=N and approaches 1=
p
N for smooth matrix fields . The proof relies,
among other things, on a new stability estimate for the inverse map C ! .
Key applications of our results are discussed in the case n D 3 to polarimet-
ric neutron tomography. © 2020 The Authors. Communications on Pure and
Applied Mathematics published by Wiley Periodicals LLC
Contents
1. Introduction 2
2. Theoretical Results for the Deterministic Inverse Problem 7
3. Bayesian Inversion of Non-Abelian X-Ray Transforms 9
4. Implementation of the Algorithm 14
5. Proofs 18
Bibliography 52
Communications on Pure and Applied Mathematics, 0001–0055 (PREPRINT)
© 2020 The Authors. Communications on Pure and Applied Mathematics published by Wiley Peri-
odicals LLC
This is an open access article under the terms of the Creative Commons Attribution License, which
permits use, distribution and reproduction in any medium, provided the original work is properly
cited.
2 F. MONARD, R. NICKL, AND G. P. PATERNAIN
1 Introduction
1.1 Non-Abelian X-Ray Transforms
Our object of study is the non-Abelian X-ray transform, a mapping from a
matrix-valued field defined on a Riemannian surface with boundary .M; g; @M/,
to its scattering data C, defined at the influx boundary @CSM of M , given by
@CSM D f.x; v/ 2 TM; x 2 @M; gx.v; v/ D 1; hv; xig  0g;
where TM is the tangent bundle of M , and x denotes the outward unit normal at
x 2 @M .
We will assume that the surfaceM is simple in the sense that it is (topologically)
a disk, it has no conjugate points, and a strictly convex boundary. Strictly convex
domains in the plane (and small perturbations of them) are examples of simple
surfaces. In this context, all unit-speed geodesics1 inM exitM in finite time. This
fact allows us to identify @CSM with the space of geodesics on M , by associating
to any .x; v/ 2 @CSM the unique geodesic  passing through .x; v/.
Let WM ! Cnn be a smooth map. Given a unit-speed geodesic  W 0; T !
M with endpoints .0/; .T / 2 @M , we may define the scattering data of  on 
to be C./ WD U.0/, where U W 0; T  ! Cnn satisfies the linear system of
ODEs
PU C..t//U D 0; U.T / D id:
This problem, backward in time for convention here, is well-posed and leads to
a unique definition of U.0/, containing cumulated information about  along the
geodesic  . Note that when  is scalar, we obtain logU.0/ D R T0 ..t//dt ,
which is the classical X-ray/Radon transform of  along the curve  . Considering
the collection of all such data makes up the scattering data (or non-Abelian X-ray
transform) of , viewed here as a map
CW @CSM ! Cnn;
and we are concerned with the problem of recovering  from C. Inverting
Abelian and non-Abelian X-ray transforms are examples of inverse problems in
integral geometry, an active field permeating several tomographic imaging meth-
ods; see, e.g., the recent topical review [23].
The problem of inverting the nonlinear mapping  7! C in this generality has
been recently solved in [35]. Previous injectivity results were obtained, either by
adding curvature conditions on the manifold or by fixing a Lie group G (realised
as matrices, for simplicity) and its Lie algebra g, in turn asking whether a g-valued
field  can be recovered from its G-valued scattering data C. In this paper, we
will mainly use the Lie groups SO.n/ D fU 2 Rnn; U TU D id; detU D 1g,
U.n/ D fU 2 Cnn; U U D idg, and SU.n/ D U.n/ \ fdet D 1g, and their Lie
1 Unit-speed geodesics are locally defined dynamically through the equation r P P D 0 with r
the Levi-Civita connection and satisfying g.t/. P.t/; P.t// D 1 for all t where .t/ is defined.
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 3
algebras so.n/ D fA 2 Rnn; AT C A D 0g, u.n/ D fA 2 Cnn; A C A D 0g,
and su.n/ D u.n/ \ ftr D 0g. Above, T, , det, and tr refer to matrix ‘transpose’,
‘conjugate transpose’, ‘determinant’, and ‘trace’, respectively. Note the inclusions
SO.n/  SU.n/  U.n/:(1.1)
The state of the art on this question can be written as follows:
THEOREM 1.1. Let .M; g/ be a simple surface. The map  7! C is injective in
the following cases:
(a) G D U.n/ [34],
(b) G D GL.n;C/ [35].
The proof of (b) consists of a reduction to the unitary case in (a) via a factor-
ization theorem in loop groups. Earlier injectivity results have been obtained by
several authors; cf. [14, 32, 33] and references therein, particularly when .M; g/ is
a domain in the Euclidean plane.
The absence of concrete reconstruction formulas for the inverse map C ! 
when n  2, and the challenge of dealing with physical experiments such as those
arising in polarimetric neutron tomography, where N discrete and noisy measure-
ments DN  PN of C 2 SO.3/ are made (see Section 1.2 and Section 1.3 for
details), motivate the main contribution of this article, which is to present a statis-
tical algorithm x.DN / that allows to recover . The implementation of x.DN /
is detailed in Section 4, and our main theoretical result is the statistical analogue
of the injectivity result Theorem 1.1, namely the frequentist consistency of recon-
struction in the large sample limit, which somewhat informally can be stated as
follows:
THEOREM 1.2. Suppose the dataDN is generated from the probability distribution
PN0 where 0 WM ! so.n/ is any smooth matrix field 0. Then we have that, as
sample size N !1, and in PN0-probability,x.DN /  0L2.M/ ! 0:
See Theorem 3.2 in Section 3 for a fully rigorous statement of this result, which
in fact requires significantly weaker hypotheses on 0, and also specifies an ex-
plicit ‘algebraic’ rate of convergence N  in the last limit.
The proof of the previous theorem relies on ideas from Bayesian nonparamet-
ric statistics [16, 43] and on new ‘quantitative versions’ of the injectivity result in
Theorem 1.1, which are of independent interest and stated in Section 2.
1.2 Polarimetric Neutron Tomography (PNT)
The basic problem in PNT consists in finding a magnetic field from spin mea-
surements of neutrons [11, 22, 25, 37]. In this case the explicit relation is
 D
24 0 B3  B2 B3 0 B1
B2  B1 0
35
4 F. MONARD, R. NICKL, AND G. P. PATERNAIN
where B D .B1; B2; B3/ is the magnetic field. In the case of PNT one assumes
that the underlying surfaceM is just the disc in the plane (by slicing with 2D discs
one can solve the 3D problem).
The details of the experiment of polarimetric neutron tomography may be found,
e.g., in [37]. Here we give a description that is suitable for our purposes. The data
produced by the experiment is the orthogonal matrix C 1 .x; v/ D C T.x; v/ 2
SO.3/, where C.x; v/ is the scattering data described above. The significance of
this in terms of spin, is as follows: if a neutron traveling along the ray determined
by .x; v/ enters the magnetic field with a spin sin 2 S2 (S2 denotes the Euclidean
unit sphere in R3), it exits the field with spin sout D C 1 .x; v/sin 2 S2 (for an
ensemble of polarized neutrons in a magnetic field it can be shown that they be-
have like a particle with a classical magnetic moment). The magnetic field B is
defined in 3D space, but the experiment makes measurements on a 2D plane and
produces a global reconstruction by slicing. The geometry of the experiment is
thus a 2D parallel beam geometry that is easily converted into fan-beam geometry
as considered above. The question is then how to manipulate the spin to produce
the orthogonal matrix. This is done with an ingenious sequence of spin flippers
and rotators placed before and after the magnetic field being measured. The ma-
terial containing the magnetic field can also be rotated so as to produce parallel
beams from different angles. After the spin has been manipulated it goes through
an analyser; this device is essentially a spin filter that only lets those neutrons with
vertically aligned spin go through. The neutron count is then measured with a de-
tector that produces an intensity reading. The spin of the entering beam is perfectly
aligned with the spin of the analyser, so that the intensity measurement is actually
a measurement of the angle of rotation of the spin due to the magnetic field. The
key relation is given by [25, eq. 1]
(1.2) I D I0A
1
2
.1C cos'/;
where A is the attenuation of the medium, I0 is the intensity of the incoming beam,
and ' is the angle by which the spin has rotated.
The use of the spin flipper allows the measurement of
I 0 D I0A
1
2
.1   cos'/;
and from this one deduces that
cos' D I   I
0
I C I 0 ;
which then becomes an entry of our matrix C T.x; v/. By rotating by =2 and
flipping (rotation by ) one can thus produce the entire orthogonal matrix as data.
In other words, if fe1; e2; e3g is the canonical basis of 3-space, cos' gives, for all
i; j , C T.x; v/ei  ej and hence all the matrix entries. In some situations, where the
attenuation of the medium is known, the use of spin flippers is not necessary and
can be calibrated out. Assuming an additive Gaussian noise in the intensities I ,
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 5
equation (1.2) approximately produces an additive Gaussian noise in the entries of
the matrix C, which is precisely the noise model we adopt below.
As in the articles [12,37], our approach reconstructs 3D magnetic fields of arbi-
trary direction and distribution. This provides a method able to investigate samples
without imposing any a priori knowledge of the magnetic field orientation and re-
quires understanding of the full nonlinear inverse problem. The recent preprint [12]
introduces a modified Newton-Kantorovich type algorithm for the solution of the
nonlinear problem, a Newton-type algorithm where the inversion of the Jacobian
at each iteration only uses the differential of the map  7! C at the base point
0  0.
As pointed out in [12], the algorithm appears to work well for small enough
fields (or large enough velocities of neutrons), but may fail due to ‘phase wrapping’
when the field is large enough. Our approach does not exhibit this problem.
1.3 The Statistical Observation Scheme
Consider a simple surface M as above with influx boundary @CSM , and a
matrix-valued map
 WM ! g
and scattering data
C W @CSM ! G:
Here we take G D SO.n/ for some n  2, with corresponding Lie algebra g D
so.n/, the set of skew-symmetric matrices. Recall that in the key application to
PNT from the previous subsection, M is the flat disk and n D 3. We could take
G D SU.n/ and g D su.n/ just as well, but for the sake of conciseness prefer to
avoid a complex-valued statistical noise model in what follows.
To describe the statistical observation setting, let  be the uniform distribution
(volume element) on @CSM (see (1.5) below for a precise definition), and consider
‘design’ random variables
.Xi ; Vi /
N
iD1 i:i:d:  on @CSM:
These draws represent a randomised choice of the geodesics for which experiments
are performed—they have to be equally spaced in a statistical sense throughout
geodesic space @CSM . For each resulting measurement of C..Xi ; Vi // the sta-
tistical observational error arising in the experiment is modeled by independent
Gaussian matrix noise. More precisely, let
."i;j;k W 1  j; k  n/NiD1
be i.i.d. N.0; ff2/; ff > 0, random variables that are independent of the .Xi ; Vi /’s,
and let Ei D ."i;j;k/ be the random n  n noise matrix that adds a Gaussian noise
variable in each matrix entry to C..Xi ; Vi //. Our observations then consist of the
sequence of N random n  n matrices
(1.3)
Yi D .Yi;j;k/; Yi;j;k D C..Xi ; Vi //j;k C "i;j;k;
i D 1; : : : ; N; 1  j; k  n:
6 F. MONARD, R. NICKL, AND G. P. PATERNAIN
The variables Yi;j;k are all independent, and even i.i.d. for j; k fixed. Conditionally
on .Xi ; Vi / D .xi ; vi / they are multivariate normal random variables with diagonal
covariance and (vectorised) mean Cffi.xi ; vi /j;k . Note that while C.x; y/ takes
values in SO.n/, the Yi are not in SO.n/ as we have not constrained Ei at all—this
is in line with the physical experiments for PNT described in Section 1.2 where
statistical errors arise from noisy measurements of each matrix entry of C.x; v/.
For the theory we will assume that the noise variance ff2 > 0 is fixed and known—
in practice it can be replaced by the estimated sample variance of the Yi;j;k’s.
To fix notation: The joint law of the random variables .Yi ; .Xi ; Vi //NiD1 in (1.3)
on .Rnn  @CSM/N will be denoted by PN D NiD1P i, where we note P i D
P 1 for all i . We also write P
N
" for the law of the .Ei /NiD1’s, N for the law of the
.Xi ; Vi /
N
iD1, and
(1.4) DN D fY1; : : : ; YN ; .X1; V1/; : : : ; .XN ; VN /g
for the full data vector. The corresponding expectation operators are obtained by
replacing ‘P ’ by ‘E’ in the preceding expressions. The dependence on ff2 will be
suppressed in the notation.
1.4 Some Geometric Background and Basic Notation
We conclude this section by introducing some more basic notation that will be
used throughout.
Our background geometry is a simple surface with boundary .M; g; @M/. By
‘simple’, we mean (i) M is nontrapping (in the sense that every maximal geodesic
in M has finite length), (ii) M has no conjugate points, and (iii) @M is strictly
convex (i.e., @M has positive definite second fundamental form). We denote by
SM the unit tangent bundle of M , namely
SM D f.x; v/ 2 TM; gx.v; v/ D 1g:
Its boundary @SM WD f.x; v/ 2 SM W x 2 @M g can be split into ‘influx’
and ‘outflux’ boundary, depending on whether the tangent vector points inside or
outside, namely we define, for x is the outer unit normal at x 2 @M ,
@SM WD f.x; v/ 2 @SM W hv; xig  0g:
The manifolds M , @M , SM , and @CSM all carry natural volume elements,
allowing us to define L2 spaces below. Specifically, the Riemannian metric g
induces an area form dx on M and restricts to a metric on @M . The unit sphere
bundle SM carries the volume element d3 D dx dv where dv is the length
element in the unit circle Sx  TxM . Finally, the boundary @SM of SM carries
the area form d2 D ds dv where dv is as above and ds is the arclength (w.r.t.
the metric g) along the boundary. Its restriction to @CSM will be denoted by
(1.5)   1
Area.@CSM/
d2j@CSM :
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 7
The spaces Cn and Cnn will be equipped with the canonical Hermitian inner
product h ; i and induced norm j  j. For elements in Cnn, this corresponds to the
Frobenius norm jAj2F WD tr.AA/ D
Pn
i;jD1 jAi;j j2, which is U.n/-invariant in
the sense that for any U 2 U.n/ and A arbitrary, jAU jF D jUAjF D jAjF .
Given .N; h/ a d -dimensional Riemannian manifold (eitherM , @M , SM , @SM ,
or @CSM as explained above), one may adapt the usual function spaces to Cn- or
C
nn-valued functions as follows: L2.N;Cnn/, L1.N;Cnn/ with norms
kU k2
L2
WD
Z
N
jU j2F dVolh; kU kL1 WD sup
y2N
jU.y/jF :
One may differentiate functions using partial derivatives f@yj gdjD1 in coordinate
charts, or equivalently, using fTj gdjD1 a global basis of smooth vector fields on
N that pairwise commutes (it will be useful to adopt the latter viewpoint in later
sections). Given a d -index  D .1; : : : ; d /, one may define jj D 1C  Cd
and T  D T 11   T dd . The metric h equips N with a distance function dh.x; y/,
and for   0, we can thus define Hölder spaces C  .N;Cnn/ with norm
kU kC D
X
jjbc
sup
y2N
T U.y/
F
C
X
jjDbc
sup
x¤y2N
jT U.x/   T U.y/jF
dh.x; y/
 bc
;
with the second term removed when  is an integer. We will also use L2-based
Sobolev spaces H s.N;Cnn/ with norm
kU k2H s D
X
jjs
kT U k2
L2
for s 2 N, and defined by interpolation otherwise (see, e.g., [41, chap. 4]).
As above, when clear from the context, the domain and/or codomain will be
dropped from the notation. In the following sections, spaces of functions with
codomain SO.n/, SU.n/, or their Lie algebras will make use of the same topology
of the corresponding spaces of Cnn-valued functions. The c-subscript attached
to a space of maps defined on M denotes the linear subspace of those maps that
vanish identically outside of a compact subset of the interior M int of M .
2 Theoretical Results for the Deterministic Inverse Problem
When discrete measurements of the forward data C are corrupted by statistical
noise, the injectivity result Theorem 1.1 is not useful for reconstructing  from
the observations, and we will discuss in the next section how to develop statistical
methods that consistently solve this statistical inverse problem. The proofs that
substantiate these methods are based on quantitative versions of Theorem 1.1—
stability estimates—as well as continuity properties of the forward map, and we
describe in this section the analytical results we obtain.
8 F. MONARD, R. NICKL, AND G. P. PATERNAIN
The results to follow hold when the codomain of the matrix fields is the largest of
the three compact Lie groups introduced before Theorem 1.1, namely U.n/ (with
Lie algebra u.n/); see equation (1.1).
THEOREM 2.1. Let .M; g/ be a simple surface. Given two matrix fields  and 
in C 1.M; u.n//, there exists a constant c.;/ such that
k  kL2.M/  c.;/
CC 1   idH1.@CSM/;
where c.;/ is a continuous function of kkC1 _ kkC1 , explicitly
c.;/ D C1.1C .kkC1 _ kkC1// eC2.kkC1_kkC1 /;(2.1)
and where the constants C1; C2 only depend on .M; g/.
The proof of Theorem 2.1 initially follows the approach for obtainingL2 ! H 1
stability estimates for the geodesic X-ray transform I as presented in [39, theorem
3.4.3]. Our starting point is the pseudo-linearisation formula
CC
 1
 D idC I.;/.  /
where I.;/ is a geodesic X-ray transform with suitable weights; see Lemma
5.5. To prove Theorem 2.1 it suffices to show that
k  kL2.M/  c.;/kI.;/.  /kH1.@CSM/:
To this end, we use the energy identity (Pestov identity) developed in [34] for ma-
trix weights arising for connections and matrix fields. The presence of the weights
produces additional terms in the identity that need to be controlled to obtain the es-
timate above, and this is where most of the work lies. The main idea for controlling
them comes from [34], where a connection with the right curvature is artificially in-
troduced to control these terms. The connection is later removed by using (scalar)
holomorphic integrating factors whose existence is guaranteed by the microlocal
properties of the normal operator associated to the geodesic X-ray transform act-
ing on functions. Taming these integrating factors has a cost that is reflected in the
constant c.;/ given in (2.1).
For the proof of Theorem 3.2 below, we also require ‘forward’ estimates in
Sobolev and Hölder scales. These are less sophisticated in nature than the stability
estimate above, and hold under less restrictive assumptions. Recall that .M; g/
is said to be nontrapping if there is no geodesic with infinite length (any simple
manifold is nontrapping).
THEOREM 2.2. Let .M; g/ be a nontrapping surface with strictly convex boundary.
For any integer k  0 and for every; 2 C k.M; u.n//, the following continuity
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 9
estimates hold:
kC   CkHk.@CSM;Cnn/(2.2)
. .1C kkCk C kkCk /2kk  kHk.M;Cnn/;
kC   CkCk.@CSM;Cnn/(2.3)
. .1C kkCk C kkCk /kk  kCk.M;Cnn/;
where by. we mean that the inequality holds with some constant that only depends
on .M; g/ and k.
In fact, in the proof of Theorem 3.2 we shall use instead of Theorem 2.1 the
following corollary of the previous two results:
COROLLARY 2.3. Under the same hypotheses as in Theorem 2.1 and c.;/ as
in (2.1), then
(2.4) k  kL2.M/  C 0c.;/.1C kkC1/kC   CkH1.@CSM/;
where C 0 is independent of  or .
3 Bayesian Inversion of Non-Abelian X-Ray Transforms
3.1 Main Results
The main goal of this section is to introduce a method to infer the matrix field
 2 C.M; so.n// from discrete observations DN of the scattering data C de-
scribed in Section 1.3. We follow the general paradigm of Bayesian inverse prob-
lems advocated by A. Stuart [10, 40], which is also related to the paradigm of
Bayesian numerical analysis [4, 13] in the noiseless case (ff D 0). The idea is
to start from a Gaussian process prior  for the parameter  and to use Bayes’
theorem to infer the best posterior guess for  given data DN .
Below we will state a theorem to the effect that the posterior mean fields xN D
EjDN  corresponding to a flexible class of Lie-algebra-valued Gaussian pro-
cess priors  for  consistently recover the ‘true’ 0 in the frequentist large sam-
ple limit as N !1, when noisy experiments have been performed under PN0 in
the model (1.3). In fact, we will provide a stochastic convergence rate to 0 of the
recovery error that is algebraic in inverse sample size 1=N .
The proof of Theorem 3.2 below provides a template to establish rigorous sta-
tistical guarantees for the Bayesian approach to other nonlinear inverse problems
as well. See Section 5.4 and Remark 3.6 for more discussion.
We emphasise that obtaining probabilistic consistency under PN0 entails ap-
proximate uniformity of the design .Xi ; Vi / and rules out ‘adversarial’ designs.
Fixed (nonrandom) design .xi ; vi / that is sufficiently ‘equally spaced’ throughout
@CSM could be considered as well in the theory that follows, either via appealing
to asymptotic statistical equivalence results in nonparametric regression [36] or by
10 F. MONARD, R. NICKL, AND G. P. PATERNAIN
tracking the numerical discretisation error explicitly through all the proofs that fol-
low. For the purposes of the present paper we opt for the random design setting, as
it allows for a cleaner, unified probabilistic treatment of the measurement process.
To introduce the Bayesian approach more concisely, consider a prior  for a
vector field .B1; : : : ; Bxn/ by prescribing a Borel probability measure on the space
xnjD1C.M/ where
xn D n.n   1/
2
D dim.so.n//:
The natural isomorphism between xnjD1C.M/ and the space C.M; so.n// of con-
tinuous functions from M to so.n/ in turn generates a prior  for  by forming a
so.n/-valued field from the Bi . For instance, in the case n D 3 so that also xn D 3,
relevant in PNT, we construct  from
.x/ D
24 0 B3.x/  B2.x/ B3.x/ 0 B1.x/
B2.x/  B1.x/ 0
35; x 2M:(3.1)
Then we make the Bayesian model assumption that
.Yi ; .Xi ; Vi //
N
iD1j  PN on .Rnn  @CSM/N ;
which by Bayes’ rule generates on C.M; so.n// a conditional posterior distribution
ofj.Yi ; .Xi ; Vi //NiD1—it will be denoted by.  j.Yi ; .Xi ; Vi //NiD1/  .  jDN /.
The posterior distribution arises from a dominated family of probability measures
(see (5.43) below) and is hence given by
.AjDN /  .AjY1; : : : ; YN ; .X1; V1/; : : : ; .XN ; VN //
D
R
A e
`N ./d./R
e`N ./d./
;
(3.2)
for any Borel set A in C.M; so.n//. Here
`N ./ D
X
iN
`i ./; where
`i ./ D  
1
2ff2
X
1j;kn

Yi;j;k   C..Xi ; Vi //j;k
2
;
(3.3)
is, up to additive constants, the log-likelihood function of the observations.
While what precedes was not specific to the choice of a particular prior, the main
theorem to follow will hold for priors arising from certain so.n/-valued Gaussian
processes. These will be constructed from a Gaussian base prior 0 from which
the coordinates Bj of xnjD1C.M/ will be drawn independently. In fact, we will
require draws from 0 to have -Hölder-continuous sample paths on M almost
surely. We refer, e.g., to [18, secs. 2.1 and 2.6] for the basic definitions of Gaussian
measures and processes and their reproducing kernel Hilbert spaces (RKHS).
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 11
Condition 3.1. For  > 0 and  >  C 1, let 0 be a centred Gaussian Borel
probability measure on the Banach space C.M/ that is supported in a separable
(measurable) linear subspace of C  .M/, and assume its RKHS .H; kkH/ is con-
tinuously imbedded into the Sobolev space H.M/.
See Remark 3.4 for concrete examples and constructions of such Gaussian process
priors with ‘maximal choice’ H D H.M/ and arbitrary  >  C 1.
Now given a random draw f 0  0 we define a new random function
(3.4) B.x/ D BN .x/ D
f 0.x/p
N 1=.C1/
; x 2M; f 0  0;
and denote its law in C.M/ by B D B;N . Then let B1; : : : ; Bxn be random
functions on M drawn as i.i.d. copies from B , and let the prior  D xnjD1B
for  be the resulting centred Gaussian product probability measure in the space
C.M; so.n// ' xnjD1C.M/ (see (3.1) for n D 3). Shrinking the prior towards the
origin in a N -dependent way as in (3.4) is crucial in our proofs; see Remark 3.5
for discussion.
The following theorem gives a bound for the convergence rate of the posterior
mean:
(3.5) xN D x
 
.Yi ; .Xi ; Vi //
N
iD1
 D Ej.Yi ; .Xi ; Vi //NiD1
towards the true field 0 in L2.M/-loss, under the law PN0 of the observations.
Note that this mean (expected value) is understood in the usual sense of Bochner in-
tegrals, and hence x takes values inC.M; so.n//. For fixed data vector Yi ; .Xi ; Vi /
and since for C 2 SO.n/ the norms kCkL1 are bounded by a fixed constant,
this expected value exists almost surely by (3.2) and a basic application of Fer-
nique’s theorem (see [18, exer. 2.1.5]). Let us say  2 H if all matrix entries of 
are contained in H.
THEOREM 3.2. Suppose the Gaussian prior  for  arises as after (3.4) with
base prior 0 satisfying Condition 3.1 for  > C 1,  > 2. Let xN be the mean
(3.5) of the posterior distribution.  j.Yi ; .Xi ; Vi //NiD1/ arising from observations
(1.3). Assume 0 2 C .M; so.n// \H. Then we have, for some  > 0,
PN0
 xN  0L2.M/ > N ! 0 as N !1:
The proof is given in Section 5.4. We note that the constraint  > 2 (and
hence  > 3) could be relaxed to  > 1 (and hence  > 2) at the expense of
more technical proofs (see Remark 5.20). We further remark that in the proof we
establish in particular that the random posterior measure .  j.Yi ; .Xi ; Vi //NiD1/
on C.M; so.n// concentrates with probability approaching 1 in a N -diameter
L2.M/-ball centred at 0; see Theorem 5.19.
12 F. MONARD, R. NICKL, AND G. P. PATERNAIN
3.2 Remarks and Discussion
Remark 3.3 (The exponent ). In the proof (see (5.66)) we show that
 <

.2 C 2/
x   1
x for any integer
x s.t. 1 < x < ;
is permitted in the previous theorem. If 0 2 C1.M/ D
T
>0H
.M/ and if
we take priors that verify Condition 3.1 for large enough ;  andH D H.M/
(possible by Remark 3.4), then we can make  as close to 1=2 as desired, and it is
easy to show that  D 1=2 cannot be improved upon by any algorithm. So at least
for smooth 0 the recovery guarantee from Theorem 3.2 is (near-)optimal. In the
‘low regularity case’ where  is not large, our bound for  may not be optimal. A
conjecture for the optimal value for  can be obtained from the much simpler linear
and Abelian case (n D 1) corresponding to the classical Radon transform, which is
treated in [31, exam. 2.5], where the exponent  D =.2 C 3/ is attained, which
can be shown to be optimal in this special case.
Remark 3.4 (Construction of Gaussian priors). We describe here some Gaussian
process priors verifying Condition 3.1 with H D H.M/.
As a first basic example consider the case where M equals the unit disk D D
f.x1; x2/ 2 R2 W x21 C x22  1g in R2 with ‘flat’ (Euclidean) geometry, relevant in
PNT. For arbitrary  > 0we can then take for0 the restriction toD of a stationary
Gaussian process onR2 with appropriate (Whittle-)Matérn covariance function k
(see [16, p. 313] and Section 4 below). This gives a Gaussian prior on C.D/ with
RKHS H equal to the space of restrictions to D of elements of H.R2/ (using
exercise 2.6.5 in [18]). This space is well-known (e.g., [41, chap. 4]) to coincide
with H.D/, and the sample paths of this process lie in the separable subspace
C 0.D/ of C  .D/ for any  < 0 <    1; see [16, p. 575f] for a proof.
The preceding construction works for any smooth bounded domainD in R2. In
particular, a simple surface M is diffeomorphic to a disc and the Sobolev spaces
H.D/ andH.M/ coincide with equivalent norms—the Matérn prior can thus be
used even when M equals D equipped with a different Riemannian metric. Alter-
natively, one can embed M isometrically into a larger closed compact (boundary-
less) manifold S and use the orthonormal basis of eigenfunctions fekg of the
Laplace-Beltrami operator on S to generate Gaussian random series fS .x/ DP
k ffkgkek.x/, gk i:i:d: N.0; 1/, x 2 S; which after restriction to M and for
suitable choice of ffk > 0, generate Gaussian priorswith any prescribed Sobolev
space H.M/ as RKHS.
Remark 3.5 (Rescaled Gaussian Priors). While the use of Gaussian process tech-
niques [3, 15, 26] in the proof of Theorem 3.2 is inspired by previous work in
[42, 43] and also [17] for ‘direct’ problems, the inverse setting poses several chal-
lenges, particularly in the nonlinear case. In our proofs we show how these chal-
lenges can be overcome by shrinking common Gaussian process priors towards the
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 13
origin as in (3.4)—the shrinkage enforces the necessary additional ‘a priori’ regu-
larisation of the posterior distribution to permit the use of our stability estimates.
While similar rescaled priors have been shown to work in some ‘direct’ settings
before (they appear as special cases of the rescaled priors studied in [42], see their
theorem 3.2), in our setting they play a crucial role: Without rescaling the exponen-
tial growth in the C 1-norms of  of the constant (2.1) would render our stability
estimate useless in the proofs.
Remark 3.6 (Related literature on Bayesian nonlinear inverse problems). The study
of statistical guarantees for the Bayesian approach to nonlinear inverse problems
has seen a recent surge of interest. In the references [29, 30, 44] nonlinear inverse
problems of elliptic and parabolic type are studied. The results therein however
only hold for specific ‘uniformly bounded wavelet’ type priors—while these are
useful to develop a first theoretical understanding of Bayesian inversion algorithms,
they posit very strong a priori assumptions on the parameter of interest and the
efficient computability of the resulting posterior distribution is also unclear.
The recent reference [31] obtains convergence rate results for optimisation based
MAP estimates (see Section 4.2 for a brief discussion of those) in a general class
of nonlinear inverse problems. For nonlinear forward maps as the ones relevant
here, these MAP estimates can be difficult to compute, and at any rate may behave
quite differently from the posterior mean: The algorithm Ej.Yi ; .Xi ; Vi //NiD1
studied here is a Bochner integral with respect to an infinite-dimensional and non-
Gaussian posterior distribution, and variational ideas from optimisation cannot be
used directly in its analysis. In the proof of Theorem 3.2 we develop new tech-
niques that allow to prove convergence rates for such algorithms—see Section 5.4
for a discussion of the key ideas that are relevant in other settings, too. Indeed,
the very recent references [1, 19] have already succeeded in adapting our proof
template to other nonlinear inverse problems. For instance, [1] studies statistical
versions of a conceptually related boundary value problem arising with electrical
impedance tomography (‘Calderón problems’). Our results imply that statistical
inversion of non-Abelian X-ray transforms (for ‘smooth parameters’ ) admits
better (i.e., polynomial) convergence rates than the necessarily logarithmic (in in-
verse noise level) recovery guarantees derived in [1] for the Calderón problem (with
smooth conductivities).
Remark 3.7 (Towards uncertainty quantification). Theorem 3.2 also serves as a
starting point to prove more refined Bernstein–von Mises theorems entailing that
the posterior distribution is approximated in a suitable infinite-dimensional space
by a canonical Gaussian measure (cf. [5, 6]). For a nonlinear elliptic inverse prob-
lem a first result of this kind was recently proved in [29], and for the linearisation
of the nonlinear problem considered here, such results were obtained in [28]. In
principle, joining the ideas of [28,29] with the techniques of the present paper, one
can conjecture that Bernstein–von Mises theorems should also hold true for the
case of non-Abelian X-ray transforms—this is the subject of ongoing research.
14 F. MONARD, R. NICKL, AND G. P. PATERNAIN
FIGURE 4.1. Left to right: An example of mesh with Nv D 886 ver-
tices. Some geodesics for the metric we use in the examples that follow.
A contour plot of the ‘sound speed’ c D e x is superimposed.
4 Implementation of the Algorithm
In this section, we present some numerical reconstructions of an su.2/-valued
matrix field  from its noisy scattering data C 2 SU.2/. In this case,  is
generated by three real-valued components B1; B2; B3 through the relation  D
B1ff1 C B2ff2 C B3ff3, where we have used the basis of su(2):
ff1 D
1
2

i 0
0  i

; ff2 D
1
2

0 1
 1 0

; ff3 D
1
2

0 i
i 0

;
with structure equations ff1; ff2 D ff3, ff2; ff3 D ff1, and ff3; ff1 D ff2. The
approach presented easily adapts to any so.n/-, su.n/-, or u.n/-valued field (in-
cluding the so.3/-valued case of polarimetric neutron tomography, a close cousin
of the present case), with some minor Lie-group-specific modifications to be made
for an accurate computation of forward data.
4.1 Numerical Domain and Forward Operator
The computational domain is an unstructured triangular mesh discretising the
unit disk M D fx2 C y2  1g made of Nv vertices, and functions on it are
piecewise linear, uniquely determined by their values at the vertices; see Figure 4.1.
In particular,  is regarded as an element of R3Nv .
The metric is isotropic, written as g D e2x.x;y/id, with scalar function x given
by
x.x; y/ D 0:3.e ..xC0:3/2Cy2/=2fi2   e ..x 0:3/2Cy2/=2fi2/; fi D 0:25:
Such an example can be seen to be nontrapping and have no conjugate points and
a strictly convex boundary; i.e., .M; g/ is simple. The case of Euclidean geometry
would correspond to x  0. Geodesic (data) space, modeled as @CSM , is param-
eterised in fan-beam coordinates .; / 2 .0; 2/  . =2; =2/ (with uniform
probability measure d D d d=.22/).
Below we will draw N geodesics uniformly at random, characterised by N ini-
tial conditions .i ; i / 2 @CSM , 1  i  N , and our statistical algorithm will
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 15
require numerical evaluation of the forward data C.i ; i /, which we now de-
scribe: Out of each data point .i ; i /, we first compute a geodesic using a forward
scheme with step size h to solve a discretisation of the system
Px.t/ D e x cos ; Py.t/ D e x sin ; P.t/ D e x.  sin @xxC cos @yx/;
with initial condition x.0/ D cosi , y.0/ D sini , and .0/ D i CCi , until
the geodesic exits the domain. This produces a discretised geodesic
i D f.xi .tj /; yi .tj //; tj D jh; 0  j  Jig:
Once such a geodesic is computed, we must then discretise the matrix ODE
PU.i .t/; Pi .t//C.i .t//U.i .t/; Pi .t// D 0; U.i .0// D id:
(The problem here is forward in time unlike that given in the introduction, though
since  is u.n/-valued, this amounts to computing the conjugate transpose of C,
which leads to the same problem.)
To discretise the above ODE, we denote U .i;j / WD U.i .tj /; Pi .tj // and imple-
ment the scheme
U .i;j / D exp. h.i;j 1//  U .i;j 1/; 1  j  Ji ;(4.1)
where we have defined .i;j 1/ D .xi .tj 1/; yi .tj 1//. In fact, the code im-
plements a predictor-corrector variant of this scheme for improved accuracy on the
computation of the exponentials.
The use of matrix exponentials in (4.1) (compared to standard forward-marching
schemes) ensures that the matrix solutionU numerically remains in SU.2/, and the
computation of these exponentials can be done via an explicit formula; namely, for
A D aff1 C bff2 C cff3 and denoting jaj WD
p
a2 C b2 C c2, we have for l 2 R
exp.lA/ D cos

l jaj
2

idC sinc

l jaj
2

lA .sinc x WD .sin x/=x/:
(Note that the formula above would need to be adapted if a Lie algebra g dif-
ferent from su.2/ is of interest.) The evaluation of .i;j 1/ is done by barycen-
tric combination of the values of  at the three vertices of the triangle containing
.xi .tj 1/; yi .tj 1//.
After implementing (4.1), the scattering data C.i / is nothing but U .i;Ji / (in
fact, the other values U .i;j / for j < Ji are not kept in memory after computation).
The magnetic field  we will use in the experiments below as well as its noiseless
scattering data C are visualised in Figure 4.2.
As we will use Monte Carlo Markov chains (MCMCs) in the following section,
let us mention that once the mesh is fixed, some computations are done prior to
the MCMC, namely, all geodesics as well as the triangle indices and barycentric
weights along them.
16 F. MONARD, R. NICKL, AND G. P. PATERNAIN
FIGURE 4.2. Top: the three components .B1; B2; B3/ of the magnetic
field realised as 0 D B1ff1 C B2ff2 C B3ff3. Bottom: real (left 2  2
block) and imaginary (right 22 block) parts of the scattering data C0 W
@CSM ! SU.2/ for the magnetic field 0 visualised on top.
4.2 Statistical Estimation through MCMC
Given data as in (1.3), a common approach to inverse problems would be to
compute a Tikhonov regulariser that minimises a penalised least squares fit func-
tional (with, e.g., Sobolev-norm penalty)
(4.2) QN ./ D
1
2ff2
NX
iD1
jYi   C.Xi ; Vi /j2F C
1
2
kk2H
over the space of all matrix fields WM ! g where g is the Lie algebra describing
the constraint on the codomain of . The map QN is not convex, and efficient
computation of the global minimiser may be challenging. One approach would be
to use a gradient-based iterative scheme [24], but the algorithmic stability of these
(or other variational) methods is unclear in the setting considered here.
The optimiser of the functional (4.2) can be shown to correspond to a posterior
mode, or ‘maximum a posteriori estimate (MAP)’, of a Gaussian process prior 
on C.M; g/ with RKHS equal to H (see [9] for a general result of this kind). In-
stead of computing that maximiser, one may compute other posterior characteris-
tics such as the posterior mean (average) EjDN  D Ej.Yi ; .Xi ; Vi //NiD1,
which in our nonlinear setting is different from the MAP estimate.
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 17
For Gaussian priors, MCMC algorithms such as the preconditioned Crank-Nicol-
son (pCN) method (see [7]) are available to sample from the posterior distribution.
To introduce the algorithm, note that as in (3.3), the log-likelihood function given
the data .Yi ; .Xi ; Vi //NiD1 equals, up to additive constants,
`./ WD   1
2ff2
NX
iD1
jYi   C.Xi ; Vi /j2F :
One then approximates the posterior mean Ej.Yi ; .Xi ; Vi //NiD1 by a Monte
Carlo average b D 1
Ns
PNs
nD0n of a Markov chain .n/ of lengthNs as follows:
Let  be a Gaussian prior for ; initialise n D 0 for n D 0, then repeat:
(1) Draw    and for  > 0 define the proposal pn WD
p
1   2n Cp
2.
(2) Set
nC1 D
(
pn with probability 1 ^ exp.`.pn/   `.n//;
n otherwise:
The algorithm is terminated at n D Ns and requires evaluation of `.n/ and
thus of the scattering data Cn.Xi ; Vi / for every n and .Xi ; Vi /. For g D su.2/
relevant in the simulations that follow, this can be done as described in Section 4.1.
The invariant measure of the Markov chain fng equals the posterior distri-
bution .  jDN /, and under certain conditions that are compatible with our set-
ting, [21] derived dimension-free spectral gaps that imply that the distribution
of n mixes rapidly towards .  jDN /. The approximation of EjDN  byb D 1
Ns
PNs
nD0n can thus be expected to compare to the one of the standard
central limit theorem, with corresponding nonasymptotic error guarantees; see sec-
tion 4 in [21].
To perform numerical simulations, we discretise DP3iD1Biffi WM ! su.2/
as in Section 4.1 and for each Bi choose an independent Matérn prior (cf. Remark
3.4) with parameters .; `/, which on functions on the mesh (i.e., vectors in RNv )
uses the covariance matrix Ci;j D k;`.jxi   xj j/ for 1  i; j  Nv, with positive
definite kernel
k;`.r/ WD
21 
./
 p
2r
`
!
K
 p
2r
`
!
;
with K the modified Bessel function of the second kind. The constant  con-
trols the Sobolev regularity while ` controls the characteristic length scale of the
samples; see Figure 4.3 for an illustration.
We draw N geodesics at random according to the uniform law for .; / (some
samples on @CSM of size N D 200; 400; 800 are visualised in Figure 4.4), and
then generate synthetic data .Yi ; .Xi ; Vi //NiD1 as explained in Section 4.1 for the
18 F. MONARD, R. NICKL, AND G. P. PATERNAIN
FIGURE 4.3. Left to right: two Matérn prior samples with ` D 0:1, 0:2
and 0:3, respectively. In all samples,  D 3.
FIGURE 4.4. Left to right: Examples of sample draws on @CSM for
N D 200; 400; 800.
magnetic field0 displayed in Figure 4.2, adding Gaussian noiseN.0; ff2/ to each
matrix entry of C0 .
We then implement the pCN algorithm to approximately compute the posterior
mean xN D Ej.Yi ; .Xi ; Vi //NiD1 from Theorem 3.2. The step size  is ad-
justed so that after ‘burn-in’, the acceptance rate of proposals stabilises around
25%. Once the chain is computed we visualise b D 1
Ns
PNs
nD0n—examples of
outcomes corresponding to increasing data set are given in Figure 4.5, illustrat-
ing the improvement in ‘reconstructions’ as the number N of measurement points
increases.
5 Proofs
5.1 Geometric Preliminaries
Let .M; g/ be a compact, oriented, two-dimensional Riemannian manifold with
smooth boundary @M . As before, SM will denote the unit circle bundle that is a
compact 3-manifold with boundary given by
@.SM/ D f.x; v/ 2 SM W x 2 @M g:
We let X be the geodesic vector field, i.e., the infinitesimal generator of the geo-
desic flow of M . Since M is assumed to be oriented, there is a circle action on
the fibers of SM with infinitesimal generator V called the vertical vector field. It
is possible to complete the pair X; V to a global frame of T .SM/ by considering
the vector field X? WD X; V . There are two additional structure equations given
by X D V; X? and X;X? D  V where  is the Gaussian curvature of the
surface. Using this frame, we can define a Riemannian metric on SM by declaring
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 19
FIGURE 4.5. Top to bottom: The posterior mean field y for sample
sizes N D 200; 400; 800, to be compared with the true field 0 from
Figure 4.2. The number of Monte-Carlo iterations equals Ns D 100000.
Other parameters: ff D 0:05,  D 3, ` D 0:2,  D 0:000025, Nv D 886.
fX;X?; V g to be an orthonormal basis, and the volume form of this metric will
be denoted by d3. The fact that fX;X?; V g are orthonormal together with the
commutator formulas implies that the Lie derivative of d3 along the three vector
fields vanishes.
Given functions u; v W SM ! Cn, we consider the inner product
.u; v/ D
Z
SM
hu; viCn d3:(5.1)
Upon defining .x; v/ WD  gx.v; x/ for .x; v/ 2 @SM , the following formula
(known as Santaló’s formula) holds for any f 2 L1.SM/:Z
SM
f .x; v/d3 D
Z
@CSM
Z fi.x;v/
0
f .'t .x; v//dt .x; v/d
2;(5.2)
where 't is the geodesic flow.
We now discuss the manifold @CSM and its geometry. One may define a natural
frame on @CSM , given by
V WD V j@CSM ; T WD .?X C X?/j@CSM where ? WD V(5.3)
20 F. MONARD, R. NICKL, AND G. P. PATERNAIN
(T represents horizontal differentiation along the tangent direction). It is easily
seen that V; T  D 0 and that these two vector fields are orthonormal for the met-
ric on @SM induced by the metric defined on SM . In particular, .T; V / is an
orthonormal frame for @CSM , and we may define H s.@CSM I / with respect to
that frame. We now prove a useful lemma that will simplify later calculations.
LEMMA 5.1. Let .M; g/ be a nontrapping surface with strictly convex boundary.
Then the vector field X can be completed into a global, pairwise commuting frame
fX;PT ; PV g of T .SM/. This frame is smooth on SMnS@M and continuous on
SM , and satisfies PT j@CSM D T and PV j@CSM D V .
PROOF OF LEMMA 5.1. For .x; v/ 2 @CSM n S@M and t 2 .0; fi.x; v//, we
define two vector fields on SM int,
.PT /'t .x;v/ WD d't j.x;v/T.x;v/; .PV /'t .x;v/ WD d't j.x;v/V.x;v/:
Since for .x; v/ 2 @CSM n S@M the map .x; v; t/ 7! 't .x; v/ is smooth and
injective and t 2 .0; fi.x; v//, this defines global, smooth sections of T .SM int/ so
thatX;PT ; PV pairwise commute. Via direct computation of the differential of the
flow (see, e.g., [27, sec. 4.2]), one may obtain the following expressions on SM int
PT D .?/ X C  .aX?   .Xa/V /; PV D  bX? C .Xb/V;
where a;bWSM ! R satisfy
X2aC a D X2bC b D 0 .SM/;  a b
Xa Xb
j@CSM D id;
and where for h W @CSM ! C, one defines h W SM ! C through the relation
h .'t .x; v// D h.x; v/; .x; v/ 2 @CSM; t 2 0; fi.x; v/:
One further notices that the definition of PV ; PT extends by continuity to @.SM/,
with the appropriate restrictions claimed in the statement of the lemma. 
5.2 Forward Estimates and Proof of Theorem 2.2
In this section, we derive various continuity estimates for the forward map 7!
C. Recall that if the boundary @M is strictly convex, by [38, lemma 4.1.2, p. 113]
there is a constant C0.M; g/ > 0 such that
fi.x; v/  C0.x; v/ 8.x; v/ 2 @CSM:(5.4)
We start with the following basic estimates.
LEMMA 5.2 (Workhorse lemma). Let .M; g/ be a nontrapping surface with strictly
convex boundary and 2 C.M; u.n//. Suppose F 2 C.SM;Cnn/ and consider
the unique continuous solution G W SM ! Cnn to XGCG D F on SM with
Gj@ SM D 0. Then there exists a constant C1.M; g/ such that
kGj@CSMkL1.@CSM;Cnn/  kGkL1.SM;Cnn/  C1kF kL1.SM;Cnn/;(5.5)
kGkL2.SM;Cnn/  C1kF kL2.SM;Cnn/;(5.6)
kGj@CSMkL2.@CSM;Cnn/  C1kF kL2.SM;Cnn/:(5.7)
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 21
The constant C1 can be chosen as C1 D max.fi1;
p
C0/, with fi1 the diameter
of M and C0 the constant given in (5.4).
PROOF. It is easy to check that
G.x; v/ D  U.x; v/
Z fi.x;v/
0
.U 1 F /.'t .x; v// dt; .x; v/ 2 SM;
where U is the unique solution U toXU CU D 0 on SM with U j@CSM D id.
Taking the Frobenius norm, using U.n/-invariance and the fact that U is unitary,
we get
jGjF .x; v/ D
Z fi.x;v/
0
.U 1 F /.'t .x; v//dt

F

Z fi.x;v/
0
U 1 F F .'t .x; v//dt  Z fi.x;v/
0
jF jF .'t .x; v//dt:
Upon bounding the right-hand side crudely by fi1kF kL1 , this immediately im-
plies (5.5). On to theL2 estimates, applying Cauchy-Schwarz yields for all .x; v/ 2
SM , we obtain
jGjF .x; v/2  fi.x; v/
Z fi.x;v/
0
jF j2F .'t .x; v//dt  fi1
Z
x;v
jF j2F ;(5.8)
where x;v is the maximal geodesic passing through .x; v/. Now fix .x; v/ 2
@CSM and integrate the inequality above along the geodesic flow 't .x; v/ to arrive
atZ fi.x;v/
0
jGjF .'t .x; v//2 dt  fi21
Z fi.x;v/
0
jF jF .'t .x; v//2 dt; .x; v/ 2 @CSM:
Multiplying both sides by , integrating w.r.t. d2 and using Santaló’s formula
yields (5.6).
For the estimate on L2.@CSM/, looking at (5.8) for .x; v/ 2 @CSM and using
(5.4), we arrive at
jGjF .x; v/2  C0
Z fi.x;v/
0
jF j2F .'t .x; v//dt .x; v/; .x; v/ 2 @CSM:
Integrating w.r.t. d2 and using Santaló’s formula (5.2) on the right-hand side
immediately gives (5.7). Lemma 5.2 is proved. 
We now prove the main result on forward estimates, Theorem 2.2. We shall
follow the model proof of [38, theorem 4.2.1], which shows that the standard X-
ray transform I maps H s to H s . We do this in two stages: first we explain in
Section 5.2 the proof in the simpler case in which the matrix fields have support
contained in the interior of M , and then we explain in Section 5.2 how to derive
the general case.
22 F. MONARD, R. NICKL, AND G. P. PATERNAIN
Proof of Theorem 2.2 Assuming  and  with Support in the Interior of M
As a preliminary identity, given  and  two skew Hermitian matrix fields,
consider the two U.n/-valued solutions U; U such that XU CU D 0 with
boundary condition Uj@ SM D id. It is immediate to find that the relation
X.U   U/C.U   U/ D  .  /U
holds pointwise on SM , and that .U   U/j@ SM D 0. Using the fact that
.U   U/j@CSM D C   C with estimate (5.5) yields
kC   CkL1  C1k.  /UkL1 D C1k  kL1 :
Similarly, combining the observation with (5.7) yields (2.2), and we can also ob-
tain, using (5.6),
kU   UkL2  C1k  kL2 :(5.9)
To prove theC 1 continuity estimate, consider the functionW WD PV .U U/,
such that W j@CSM D V.C   C/ and for brevity set P D PV . The following
identity is immediate:
XW CW D  .P.U   U/C .P.  //U C .  /PU/:
In addition, since  and  are compactly supported in M int, the functions U,
U equal the identity matrix in a neighbourhood of @ SM and in particular,
W j@ SM D 0.
Using estimates (5.5)–(5.7) and U.n/-invariance of Frobenius norms gives
kV.C   C/kL2
 C1kP.U   U/C .P.  //U C .  /PUkL2
 C1.kPk1kU   UkL2 C kP.  /kL2 C k  kL2kPUkL1/:
We also haveX.PU/CPU D  .P/U with PUj@ SM D 0, so by (5.5),
we get kPUkL1  C1kPkL1 . Combining this fact with (5.9), we arrive at
kV.C   C/kL2
 C1
 
C1.kPkL1 C kPkL1/k  kL2 C kP.  /kL2

;
and a similar bound for kPV .U   U/kL2 . Obtaining a similar estimate for
T .C   C/, we arrive at
kC   CkH1 . .1C kkC1 C kkC1/k  kH1 :
Similar arguments using sup norms everywhere yield
kC   CkC1 . .1C kkC1 C kkC1/k  kC1 :
To proceed to higher-order derivatives, if P D P 1V P 2T is a derivative of
order jj, setting W D P.U  U/, we have W j@CSM D V 1T 2.C  C/,
W j@ SM D 0, and
XW CW D  P; .U   U/;
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 23
where the right-hand side involves derivatives of  of order at most jj, and
derivatives of U   U of order at most jj   1. Combining the estimates of
Lemma 5.2 and an induction on k (whose formulation also involves control on
kP 1V P 2T .U   U/kL2.SM/ for all 1 C 2  k, and where the commuting
frame fX;PV ; PT g avoids the proliferation of terms due to nontrivial commuta-
tors) proves the theorem for higher-order derivatives. 
Proof of Theorem 2.2  and  Supported Up to @M
Consider a compact nontrapping surface .M; g/ with strictly convex bound-
ary, and let  2 C.M;Cnn/ be a matrix-valued field. We shall call R 2
C.SM;GL.n;C// an integrating factor for  if R is differentiable along the
geodesic vector field X and XRCR D 0. Let U denote the unique integrat-
ing factor with Uj@ SM D id. Recall that C D Uj@CSM . First note that the
work of the previous section also proves for every k  0 that if  and  are C k
matrix fields compactly supported inside of M int, we also have
kU   UkCk.SM/ . .1C kkCk C kkCk /kk  kCk.M/;
kU   UkHk.SM/ . .1C kkCk C kkCk /kk  kHk.M/:
(5.10)
Let  W @SM ! @SM denote the scattering relation of the metric (i.e., the
map that takes initial conditions of a geodesic at the moment of entry to final con-
ditions at the moment of exit). If R denotes any other integrating factor for ,
then it must have the form UF ], where F ] is the first integral (i.e., XF ] D 0)
determined by F 2 C.@CSM;GL.n;C//. Thus R D UF ], and from this we
deduce
(5.11) C D R.R 1  /:
In particular, given two continuous matrix fields , , Equation (5.11) and
X.R 1 R/ D R 1 .  /R
imply the identity on @CSM :
C   C D .R  R/R 1   CR.R 1  R 1 /  (5.12)
D RI.R 1 .  /R.R 1  /;(5.13)
where I is the standard X-ray transform acting on functions in SM . To complete
the proof of Theorem 2.2 for ; supported up to the boundary, we then need to
construct integrating factors with good regularity on SM (i.e, at @0SM included)
and which behave continuously in terms of  and . To this end, we consider
.M; g/ isometrically embedded in a closed manifold .S; g/. The Seeley extension
theorem asserts that for any k  0 there is a continuous extension map
Ek W C k.M/! C k.S/; Ek WHk.M/! Hk.S/:
(It also works for C1.) We consider a slightly larger compact manifold with
boundary M  S engulfing M so that . M;g/ stays nontrapping and with strictly
24 F. MONARD, R. NICKL, AND G. P. PATERNAIN
convex boundary. We fix once and for all a smooth cutoff function  so that it has
compact support in M int and it equals 1 near M . Thus given  2 C k.M; u.n//,
z WD Ek./ 2 C kc . M; u.n//;
and since Ek is continuous,
(5.14) kzkCk . kkCk ; kzkHk . kkHk :
Now by virtue of the work in Section 5.2 applied to z on M , we can deduce
estimates of the form
kUz   idkCk . kzkkCk ; kU 1z   idkCk . kzk
k
Ck
:(5.15)
We then take as smooth integrating factors R WD UzjSM and R WD UzjSM .
Combining (5.14) and 5.15 we derive
(5.16) kR   idkCk . kkkCk ; kR 1   idkCk . kkkCk :
Combining (5.14) and (5.10) applied to Uz and Uz, we obtain
kR  RkCk  kUz   UzkCk
. .1C kzkCk C kzkCk /kkz   zkCk.M/
. .1C kkCk C kkCk /kk  kCk.M/;
(5.17)
and similarly for kR 1  R 1 kCk and forHk norms. Then the proof for Theorem
2.2 for ; supported up to the boundary consists in applying the product rule to
(5.12) (for C k norms) and (5.13) (for Hk norms) and using estimates (5.16) and
(5.17) together with the boundedness of the standard X-ray I between Sobolev
spaces [38, theorem 4.2.1].
5.3 Stability Estimate—Proof of Theorem 2.1
Setting, Main Results, and Proofs of Theorem 2.1 and Corollary 2.3
Before considering the nonlinear inverse problem, we must establish a stabil-
ity estimate for a linear inverse problem, that of reconstructing a function f 2
C1.M;Cn/ from its attenuated X-ray transform, where the attenuation is matrix-
valued. Namely, given a smooth skew-Hermitian matrix inM , we define If WD
uf

@CSM
, where u D uf W SM ! Cn is the unique solution to the problem
XuCu D  f .SM/; uj@ SM D 0:
The injectivity of such a transform was proved in [34], and we now provide a
stability estimate for it.
THEOREM 5.3. Let .M; g/ be a simple Riemannian surface with boundary and 
a smooth, skew-Hermitian matrix field in M . Then for any f 2 C1.M/, we have
the following stability estimate:
kf kL2.M;Cn/  C1.1C kkC1/eC2kkC1kIf kH1.@CSM;Cn/:(5.18)
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 25
Remark 5.4 (Dependence of C1; C2). The constants C1; C2 only depend on the
geometry of .M; g/. The constant C1 blows up like .   1/ 1, where  is the
terminator constant of .M; g/. This is one of the ways that this stability estimate
ceases to hold as one approaches nonsimplicity. The main other quantity appear-
ing in C1; C2 is w1, the sup norm of the integrating factor defined below. The
behavior of such a quantity, while finite on any simple surface, remains to be better
understood.
On to the nonlinear stability estimate, injectivity of the operator  ! C re-
stricted to u.n/-valued fields was initially proved in [34], and Theorem 2.1 up-
grades this result with a stability estimate. While the remaining sections will focus
on the proof of Theorem 5.3, we now explain how this result implies Theorem 2.1.
The main additional ingredient needed is a pseudo-linearization identity, relating
scattering data to attenuated X-ray transforms:
LEMMA 5.5 (Pseudo-linearization). Let .M; g/ be a nontrapping surface with
strictly convex boundary. For any ; 2 C.M;Cnn/, the following relation
holds:
(5.19) CC 1 D idC I.;/.  /;
where I.;/WL2.M;Cnn/! L2.@CSM;Cnn/ is an attenuated X-ray trans-
form with matrix field .;/, an endomorphism of Cnn with pointwise action
.;/  U D U   U; U 2 Cnn:
PROOF OF LEMMA 5.5. With U; U the fundamental solutions of XU C
U D 0 with Uj@ SM D id and Uj@CSM D C (similarly for ), denote
W WD UU 1   id. A direct computation shows that
XW CW  W D  .  / .SM/; W j@ SM D 0I
and thus by the definition of the attenuated X-ray transform,
W j@CSM D I.;/.  /:
Since we also have by construction W j@CSM D CC 1   id, identity (5.19)
follows. 
PROOF OF THEOREM 2.1. Appealing to the pseudo-linearization (5.19), one
may notice that if ,  are skew-Hermitian, then the field .;/ is skew-
Hermitian when viewed as an endomorphism ofCnn. Moreover, since the entries
of .;/ are linear in the entries of  and , we directly have that
k.;/kC1  C.kkC1 _ kkC1/;
26 F. MONARD, R. NICKL, AND G. P. PATERNAIN
with C a universal constant. Then relation (5.19), together with Theorem 5.3 im-
mediately implies
k  kL2.M;Cnn/
 C1.1C k.;/kC1/ eC2k.;/kC1kI.;/.  /kH1.@CSM/
 C 01.1C kkC1 _ kkC1/ eC
0
2
.kk
C1
_kk
C1
/
CC 1   idH1.@CSM/:
This shows Theorem 2.1 when; 2 C1.M; u.n//. Since all quantities involved
above do not depend on derivatives of ; of order higher than 1, and C1; C2
are independent of ;, approximating ; 2 C 1.M; u.n// by sequences in
C1.M; u.n// (and using Theorem 2.2) will yield the same stability estimate for
C 1 matrix fields. 
We also cover the proof of Corollary 2.3, based on the previous result and the
forward estimate Theorem 2.2.
PROOF OF COROLLARY 2.3. It is enough to show thatCC 1   idH1 . .1C kkC1/kC   CkH1 :(5.20)
To show this, we write at the pointwise levelCC 1   idF D .C   C/C 1 F D jC   CjF ;
hence kC CkL2 D kCC 1   idkL2 . To control first derivatives, take P D V
or T ; we haveP  CC 1   idF D P.C   C/C  id   CC 1 PCF
 jP.C   C/jF C jPC jF
id   CC 1 F
using triangle inequality and submultiplicativity. Squaring, taking the sup norm of
jPCjF , and integrating on @CSM , we obtainP  CC 1   id2L2  2 kP.C   C/k2L2 C kPCk21kC   Ck2L2:
Combining the estimates for P D V and P D T we arrive at
kCC 1   idk2H1  .1C 2kVCk2L1 C 2kTCk2L1/kC   Ck2L2
C 2kV.C   C/k2L2 C 2kT .C   C/k2L2 :
Now using the forward estimate (2.3) with k D 1 and   0 (thus C D id), we
deduce that
1C 2kVCk2L1 C 2kTCk2L1 . 1C kk2C1 :
This yields the estimate kCC 1   idk2H1 . .1 C kk2C1/kC   Ck2H1 , and
taking square roots yields (5.20) (using that
p
1C x2=.1Cx/ is uniformly bounded
for x 2 0;1/). 
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 27
Proof of Theorem 5.3—Main Outline
As in [34], the main method of proof involves an energy identity (or Pestov
identity), based on integrations by parts on SM . To do this, let us recall that with
the inner product .u; v/ defined in (5.1), and upon also denoting
.u; v/@SM WD
Z
@SM
uv d2;
the following integrations-by-parts formulas holds for u; v 2 C1.SM;Cn/:
(5.21)
.V u; v/ D  .u; V v/; .V u; v/@SM D .u; V v/@SM ;
.Xu; v/ D  .u;Xv/C .u; v/@SM ; .x; v/ WD  hv; xi:
We will also use extensively the harmonic decomposition on the fibers of SM .
Namely, the space L2.SM;Cn/ decomposes orthogonally as a direct sum
L2.SM;Cn/ D
M
k2Z
Hk
where Hk is the eigenspace of  iV corresponding to the eigenvalue k. A function
u 2 L2.SM;Cn/ has a Fourier series expansion
u D
1X
kD 1
uk;
where uk 2 Hk . Let k D C1.SM;Cn/ \ Hk . Of special interest are the
operators
 WD
1
2
.X  iX?/;(5.22)
with the property that .k/  k1 for all k 2 Z. For more details on the
operators  and the Fourier expansion, we refer the reader to [20] where these
tools were first introduced.
DEFINITION 5.6. A function u W SM ! Cn is said to be holomorphic if uk D 0
for all k < 0. Similarly, u is said to be antiholomorphic if uk D 0 for all k > 0.
To control the terms involving the matrix field, one must introduce an artificial
connection as we will see below. This first requires that we derive a Pestov identity
for X-ray transforms with connection A and matrix field.2 Namely, given a skew
Hermitian pair .A;/ on the bundle M  Cn and f 2 C1.M;Cn/, we define
IA;f D uj@CSM , where u is the unique solution to the problem
Gu D  f .SM/; uj@ SM D 0; .G WD X C AC/:
While previous Pestov identities have been derived in [34], the present one ac-
counts for nonzero boundary terms, and in particular reflects more precisely how
the stability constant degrades as .M; g/ approaches nonsimplicity. This is cap-
tured by the concept of terminator constant Ter: given a simple surface .M; g/,
2 The matrix field  is also referred to as a ‘Higgs’ field in the literature.
28 F. MONARD, R. NICKL, AND G. P. PATERNAIN
there exists a number Ter > 1 such that for any  2 .1; Ter, there exists a
smooth function r D r W SM ! R, solution to the Riccati-type equation
Xr C r2 C  D 0.
THEOREM 5.7. Let .M; g/ a simple surface with boundary, with terminator con-
stant Ter > 1 and .A;/ a skew-Hermitian pair on the bundle M Cn. Then for
any u 2 C1.SM;Cn/ and  2 .1; Ter, the following identity holds:
1

kGV u   rV uk2 C
   1

kGV uk2 C kGuk2   kVGuk2
  .?FAu; V u/  <.u;Gu/  <..?dA/u; V u/
D <.rT;Au; V u/@SM C<.hv?; iV u; u/@SM
  1

. rV u; V u/@SM :
(5.23)
In the identity above,
?dA WD X?C ; V .A/; rT;Au D T uC A.x; ?/u;(5.24)
and r is a smooth function on SM that only depends on the surface. The quantity
?FA is the curvature of the connection A, which upon a judicious choice of con-
nection, can have a controlled sign. To achieve this, consider the scalar Hermitian
connection a WD i'id, where ' is a smooth 1-form such that d' D !g (the area
form of the metric g). We choose a specific ' of the form ' D ?dh for h a real-
valued function satisfying ?d ?dh D 1 with Neumann condition dh./ D 0 at the
boundary. The latter condition implies that rT;sau D T u for any real s. Then we
have
a D i.X?h/id; a1 D Ch; a 1 D   h D  a1;
with  defined in (5.22) and i ? Fa D  1.
By [34], we can construct a holomorphic scalar function w 2 C1.SM/ sat-
isfying Xw D  iX?h. Without loss of generality, w can be chosen even. The
condition on w0 reads  .w0   h/ D 0, for which it is sufficient to use w0 D h.
With this choice of a and s 2 R, in what follows, we will denoteGs WD XCsaC
and G D G0. With w as above, we have Gsu D eswG.e swu/. Moreover, Sw (the
complex conjugate of w) is antiholomorphic and solves XSw D CiX?h, so also
Gsu D e sSwG.esSwu/.
Lastly, we will denote by the projection onto positive and negative harmon-
ics. Namely, u D
P
k>0 uk . We have the following commutator formulas
for any u 2 C1.SM/:
 ; X C saCu D .  C sa 1/u0   .C C sa1/u 1;
C; X C saCu D .C C sa1/u0   .  C sa 1/u1:
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 29
The following lemma will help us control u by versions of u that are conjugated
by special integrating factors.
LEMMA 5.8. With the holomorphic function esw and antiholomorphic function
e s
0 xw and any s; s0 2 R, we have
 u D  .e sw .eswu///; Cu D C.es0 xwC.e s0 xwu//I
in particular, we get the equality
u D u0 C .e sw .eswu//CC.es0 xwC.e s0 xwu//:(5.25)
PROOF. We only prove  u D  .e sw .eswu///, and the rest is similar.
It is enough to notice that for any holomorphic function f , the equality .f u/ D
 .f  u/ holds, as this amounts to saying that the negative harmonics of f u
do not depend on the nonnegative harmonics of u. This is immediate since
.f u/k D
X
p0
fpuk p:
Then we compute immediately
 .e
 sw .e
swu/// D   .e sweswu/ D  u;
hence the result. 
OUTLINE OF PROOF OF THEOREM 5.3. At first we are going to assume that
the solution u to the transport problem XuC u D  f , uj@ SM D 0 is C1. If
f is supported all the way to the boundary, this may not be the case, as u may fail
to be smooth at the glancing @0SM because fi is not smooth at @0SM . However,
there is a standard way to fix this issue and we shall do this at the very end. For
now we will proceed as if u were smooth in SM .
The initial transport equation, projected onto the harmonic term of degree 0,
reads
 f D Cu 1 C  u1 Cu0 D .Cu 1 Cu0=2/C . u1 Cu0=2/;
so that, in particular,
kf k2  2 kCu 1 Cu0=2k2 C k u1 Cu0=2k2:(5.26)
The crux is then to find how to bound the quantities on the right by the boundary
values of u. Using a Pestov identity with a special connection sa defined as above
(and its holomorphic integrating factor esw ), we show how to control the first term
using control over  .eswu// for s > 0. Similar work can be done to control the
second term using control over C.e s
0 xwu/ for s0 < 0.
30 F. MONARD, R. NICKL, AND G. P. PATERNAIN
We first derive in Section 5.3 the following identity:
Cu 1 C
1
2
u0
D ..eswu/.C   sa1/.e sw//0 C
1
2
e sw0.eswu/0
C i
2
.e swGsV .e
swu//0 C
i
2
e sw0.GsV .e
swu//0:
(5.27)
Since .C   sa1/.e sw/ only has strictly positive harmonic terms, the first term
in the right-hand side of (5.27) only depends on  .eswu/. Upon defining vs WD
 .e
swu/, the identity (5.27) reads
Cu 1 C
1
2
u0 D .vs.C   sa1/.e sw//0 C
1
2
.es.w w0/u/0
C i
2
.e swGsV vs/0 C
i
2
e sw0.GsV vs/0:
(5.28)
Denoting w1 D supSM jwj, we straightforwardly obtain the estimate
(5.29)
Cu 1 C 12u0
2
 C0
 jwj2
C1
se2sw1kvsk2
C jj2
C0
k.es.w w0/u/0k2 C e2sw1kGsV vsk2

;
and control on kCu 1 C 12u0k2 will be obtained after controlling each term in
the last right-hand side. We first control k.es.w w0/u/0k2 by kvsk2 C kGsV vsk2
via the estimate
(5.30)
k.es.w w0/u/0kL2.M/
 C 0e2sw1 kGsV vsk2L2.M/ C jj2C0kvsk2L2.M/ C kIf kL2.@SM/:
We then control kvsk2 and kGsV vsk2 by boundary terms via the Pestov identity
and setting up an appropriate threshold on s. To do this, we consider the transport
problem for vs , written as
Gsvs D Gs. .eswu// D Gs; .eswu/
D .C C sa1/.eswu/ 1   .  C sa 1/..eswu/0/
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 31
We then use the Pestov identity (5.23) for vs , with ?Fsa D is id and ?dsa D
X?:
1

kGsV vs   rV vsk2 C
   1

kGsV vsk2 C k.C C sa1/.eswu/ 1k2
C s.vs; iV vs/  <.vs; Gsvs/  <..X?/vs; V vs/
D <.T vs; V vs/@SM C<.hv?; iVvs; vs/@SM
  1

. rV vs; V vs/@SM :
(5.31)
Before choosing s appropriately, we need additional work (tedious as in [34]) on
the term <.vs; Gsvs/. Taking into account boundary terms, and upon defining
B1 WD , we prove in Section 5.3 that
<.vs; Gsvs/ D
1X
kD1
. 1/k j.vs/ kj2  <.B 1.vs/ k; .vs/ k 1/
C<.ex.v/.vs/ k; .vs/ k 1/@SM

;
(5.32)
with ex.v/ defined in (5.40). The last term in the sum will move to the right-hand
side of (5.31), while the other two need to be controlled with a large s. To achieve
this, we prove in Section 5.3 the following:
LEMMA 5.9. There exists a universal constantC > 0 such that for all s  C jjC1
s.vs; iV vs/  
1X
kD1
. 1/k j.vs/ kj2  <.B 1.vs/ k; .vs/ k 1/
 <..?dsa/vs; V vs/  0:
(5.33)
In particular, for s D C jjC1 C 1, identity (5.31) becomes
1

kGsV vs   rV vsk2 C
   1

kGsV vsk2
C k.C C sa1/.eswu/ 1k2 C
1X
kD1
kjv kj2
 <.T vs; V vs/@SM C<.? Vvs; vs/@SM
  1

. rV vs; V vs/@SM
 
1X
kD1
. 1/k<.ex.v/.vs/ k; .vs/ k 1/@SM
(5.34)
We now explain how to bound the right-hand side in terms of kIf k2H1.@CSM/.
Recall that vs D  .eswu/. The first claim is that  ; V  D  ; T  D 0. The
32 F. MONARD, R. NICKL, AND G. P. PATERNAIN
first one is obvious because both operators are diagonal of the fiberwise Fourier
decomposition C1.@SM/ D Lk2Z ker.id   ikV /. That T is also diagonal on
this decomposition follows from the fact that T; V  D 0. With this in mind, we
have, on @SM :
V vs D  V.eswu/ D  esw.s.V w/uC V u/;
T vs D  esw.s.T w/uC T u/;
and since uj@ SM D 0, V u and T u will be controlled by kIf kH1.@CSM/. The
right-hand side of (5.34) is thus bounded by
C 0.s2 C sjjC0 C 1/e2sw1kIf k2H1.@CSM/;
where the constant C 0 does not depend on .
Using this bound and throwing out the first and third terms of the left-hand side
of (5.34), we obtain
   1

kGsV vsk2 C
1X
kD1
kjv kj2
 C 0.s2 C sjjC0 C 1/e2sw1kIf k2H1.@CSM/:
The second term in the left-hand side controls kvskL2 directly, and we can write
(5.35)
.   1/kGsV vsk2 C kvsk2
 C 0.s2 C sjjC0 C 1/e2sw1kIf k2H1.@CSM/;
with C 0 some constant independent of . Recalling that s D C jjC1 C 1 and
combining estimates (5.26), (5.29), (5.30), and (5.35), we arrive at estimate (5.18),
completing the proof of Theorem 5.3. 
Pestov Identity with Boundary Term for Ray Transforms
with Skew-Hermitian Pairs
Let A and  be a skew-Hermitian pair, and define
G WD X C AC; G? WD X?   AV where AV WD V.A/:
We have the following structure equations:
(5.36)
G; V  D G?; V;G? D G  ;
G;G? D  V   ?FA   ?dA;
where ?dA D X?CAV  AV, or when the connectionA is scalar, ?dA D
X?, where .x/ is the Gaussian curvature. In what follows, we will need to
integrate by parts with boundary terms, and using (5.21), we obtain for G:
.Gu; v/ D  .u;Gv/C .u; v/@SM :
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 33
PROOF OF THEOREM 5.7. We first write a differential identity using the struc-
ture equations (5.36):
GV VG   VGGV D GV V;GC G; V GV
D  GVG? CG?VG
D  GV;G?C G?; GV
D  G2 CGC V 2 C ?FAV C .?dA/V;
where Gf WD G.f /. We record this here as
GV; VG D  G2 CGC V 2 C ?FAV C .?dA/V:(5.37)
Now, considering u smooth and supported up the boundary, we write
kVGuk2   kGV uk2
D .VGu; VGu/   .GV u;GV u/ D
D  .V VGu;Gu/C .GGV u; V u/   .GV u; V u/@SM
D .GV; VGu; u/   .V VGu; u/@SM   .GV u; V u/@SM
D kGuk2   .Gu; u/@SM   .u;Gu/C .u; u/@SM
C .V 2u; u/C .?FAV u; u/C ..?dA/V u; u/
  .V VGu; u/@SM   .GV u; V u/@SM
We now arrange the four boundary terms using integration by parts in V and the
formulas
V D hv?; i D ?; V 2 D V? D  :
First notice that
.V VGu; u/@SM D  .VGu; .V/u/@SM   .VGu; V u/@SM
D  .VGu;?u/@SM   .VGu; V u/@SM
D  .Gu; u/@SM C .Gu;?V u/@SM   .VGu; V u/@SM :
We then obtain
.Gu; u/@SM C .V VGu; u/@SM C .GV u; V u/@SM   . u; u/@SM
D .Gu;? V u/@SM   .VGu; V u/@SM
C .GV u; V u/@SM   . u; u/@SM
D .? GuC  G?u; V u/@SM   . u; u/@SM :
We now simplify, using that V.A/.x; v/ D A.x; v?/ and ?X C X? D T ,
? GuC  G?u D T uC A.x; ?/uC ?u DW rT;AuC ?u:
34 F. MONARD, R. NICKL, AND G. P. PATERNAIN
The boundary terms then simplify into
.? GuC  G?u; V u/@SM   .u; u/@SM
D .rT;Au; V u/@SM C .? u; V u/@SM   . u; u/@SM
D .rT;Au; V u/@SM C .?V u; u/@SM :
With this notation, the full Pestov identity takes the form
kGV uk2   .V u; V u/C kGuk2   kVGuk2   .u;Gu/
C .?FAV u; u/C ..?dA/V u; u/
D .rT;Au; V u/@SM C .?V u; u/:
(5.38)
To recover [34, eq. (8)], we take the real part of the equality above and notice that
.?FAV u; u/ D  .?FAu; V u/ because V.?FA/ D 0; then
..?dA/V u; u/ D .V ..?dA/u/; u/   .V .?dA/u; u/
D  ..?dA/u; V u/   ..dA/u; u/:
Since the last term is purely imaginary, the real parts of the other terms agree, and
upon taking the real part of (5.38), we obtain
kGV uk2   .V u; V u/C kGuk2   kVGuk2  <.u;Gu/
  .?FAu; V u/  <..?dA/u; V u/
D <.rT;Au; V u/@SM C<.?V u; u/@SM :
(5.39)
(Note that the second boundary term is purely real so the < is just ornamental.)
We finally explain how the index form term kGV uk2 .V u; V u/ can be rewrit-
ten as the sum of a nonnegative term and a boundary term. With Ter as in the
statement, and the function r D r WSM ! R solving Xr C r2 C  D 0, we
now compute, for any  2 C1.SM;Cn/
kG   r k2 D kG k2   .G ; r /   .r ;G /C kr k2:
We simplify
.G ; r /C .r ;G / D .X ; r /C .r ;X /
D
Z
SM
.X /r x C r .X x /
D
Z
SM
X. r x /   .Xr/ x 
D . r ; /@SM C
Z
SM
.r2 C / x :
We arrive at
kG   r k2 D kG k2   .r ; /@SM   . ; /;
and we may rearrange this as
.kG k2   . ; // D kG   r k2 C .   1/kG k2 C . r ; /@SM :
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 35
Plugging this last relation into (5.39) with  D V u yields (5.23). 
Remaining Estimates and Lemmata
PROOF OF EQUALITY (5.27). We write, using Lemma 5.8
Cu 1 D C.e sw .eswu// 1
D C
hX
kD0
.e sw/2k.e
swu/ 1 2k
i
D
1X
kD0
 
..C   sa1/.e sw/2k/.eswu/ 1 2k
C .e sw/2k.C C sa1/.eswu/ 1 2k

D ..eswu/.C   sa1/.e sw//0 C
1X
kD0
.e sw/2k.C C sa1/.eswu/ 1 2k :
To rewrite the last term, from the equation Gs.eswu/ D  eswf , note the relation
.C C sa1/.eswu/ 1 2k C .  C sa 1/.eswu/1 2k C.eswu/ 2k D 0:
Then we have, for k > 0,
.GsV .e
swu// 2k D V.Gs .eswu// 2k  
D0
C.Gs; V  .eswu// 2k
D  i.C C sa1/.eswu/ 1 2k C i.  C sa 1/.eswu/1 2k
D  2i.C C sa1/.eswu/ 1 2k   i.eswu/ 2k;
where we used the transport equation in the last line. For k D 0,
.GsV .e
swu//0 D  i.C C sa1/.eswu/ 1:
Plugging this back into the equation for Cu 1, we get
Cu 1 D ..eswu/.C   sa1/.e sw//0 C .e sw/0i.GsV .eswu//0
C
X
k>0
.e sw/2k

i
2
.GsV .e
swu// 2k  
1
2
.eswu/ 2k

:
We now writeX
k>0
.e sw/2k.e
swu/ 2k D
1X
kD0
.e sw/2k.e
swu/ 2k   e sw0.eswu/0
D u0   e sw0.eswu/0;
and similarlyX
k>0
.e sw/2k.GsV .e
swu// 2k
D .e swGsV .eswu//0   e sw0.GsV .eswu//0:
36 F. MONARD, R. NICKL, AND G. P. PATERNAIN
Using the last two computations, we arrive at (5.27). 
PROOF OF ESTIMATE (5.30). The transport equation for eswu projected onto
the harmonic term of degree  1 reads
.  C sa 1/.eswu/0 D  .C C sa1/.vs/ 2  .vs/ 1:
For our choice of connection, a 1 D   w0, so the left side can be rewritten as
.  C sa 1/.eswu/0 D esw0 .e sw0.eswu/0/ D esw0 .es.w w0/u/0I
hence we obtain
 .e
s.w w0/u/0 D  e sw0.C C sa1/.vs/ 2   e sw0.vs/ 1:
We then rewrite the latter right-hand side in terms of GsV vs . Notice that
.GsV vs/ 1 D .C C sa1/.V vs/ 2 C.V vs/ 1
D  2i.C C sa1/.vs/ 2   i.vs/ 1;
so
 .C C sa1/.vs/ 2 D  
i
2
.GsV vs/ 1 C
1
2
.vs/ 1;
and thus
 .e
s.w w0/u/0 D  
e sw0
2
..GsV vs/ 1 C.vs/ 1/:
Upon deriving an estimate of the form
kf kL2.M/  C.k f kL2.M/ C kf j@MkL2.@M//;
we can write
k.es.w w0/u/0kL2.M/ . k .es.w w0/u/0kL2.M/ C k.es.w w0/u/0j@MkL2.@M/
.
1
2
ke sw0..iGsV vs/ 1 C.vs/ 1kL2.M/
C k.es.w w0/u/0/j@MkL2.@M/;
and (5.30) follows. 
PROOF OF (5.32). We first need to write an integration by parts for  defined
in (5.22). Using integrations by parts (5.21) we first derive an integration by parts
for X? D XV   VX : for any u;w smooth on SM ,
.X?u;w/C .u;X?w/ D .XV u;w/   .VXu;w/C .u;XVw/   .u; VXw/
D .XV u;w/C .V u;Xw/C .Xu;w/C .u;XVw/
D .V u;w/@SM C .u; V w/@SM
D  ..V/u;w/@SM D  .?u;w/@SM :
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 37
We now compute, using that C D   
.u; Cw/C . u;w/ D
1
2
..u; .X C iX?/w/C .X   iX?u;w//
D 1
2
..C i?/u;w/@SM D .ex.v/u;w/@SM
where we define
ex.v/ WD
1
2
..x; v/C i?.x; v//:(5.40)
Similarly, for the skew-Hermitian connection considered,
.u; .C C sa1/w/C ..  C sa 1/u;w/ D .ex.v/u;w/@SM :
Now, using the fact that
.Gsvs/ 1 D .  C sa 1/.eswu/0 D  .C C sa1/.vs/ 2  .vs/ 1;
we compute
<.vs; Gsvs/ D <..vs/ 1; .Gsvs/ 1/
D <..vs/ 1; .C C sa1/.vs/ 2/   j.vs/ 1j2
D <..  C sa 1/..vs/ 1/; .vs/ 2/
 <.ex.v/.vs/ 1; .vs/ 2/@SM   j.vs/ 1j2
D  j.vs/ 1j2 C<.b 1.vs/ 1; .vs/ 2/
 <.ex.v/.vs/ 1; .vs/ 2/@SM C p1;
where p1 WD <..  C sa 1/.vs/ 1; .vs/ 2/. Upon defining
pn WD <..  C sa 1/.vs/ n; .vs/ n 1/; n  1;(5.41)
we now prove by induction the following claim:
<.vs; Gsvs/
D
nX
kD1
. 1/k j.vs/ kj2  <.b 1.vs/ k; .vs/ k 1/
C<.ex.v/.vs/ k; .vs/ k 1/@SM
C . 1/nC1pn:
(5.42)
38 F. MONARD, R. NICKL, AND G. P. PATERNAIN
The case n D 1 is proved above, and the induction step .n H) nC 1/ follows
from the calculation
pn D <..  C sa 1/.vs/ n; .vs/ n 1/
D  <..vs/ n 1; .  C sa 1/.vs/ n/
D <..vs/ n 1; .C C sa1/.vs/ n 2 C.vs/ n 1/
D j.vs/ n 1j C <..vs/ n 1; .C C sa1/.vs/ n 2/@SM
D j.vs/ n 1j   <..  C sa 1/..vs/ n 1/; .vs/ n 2/
C<.ex.v/.vs/ n 1; .vs/ n 2/@SM
D j.vs/ n 1j   <.b 1..vs/ n 1/; .vs/ n 2/
C<.ex.v/.vs/ n 1; .vs/ n 2/@SM   pnC1:
Putting this equality back into (5.42) proves the induction. Since vs 2 H 1.SM/,
we have that limn!1 pn D 0, and thus (5.32) follows. 
PROOF OF LEMMA 5.9. The term that ultimately controls everything is
s.vs; iV vs/ D s
X
k<0
jkjj.vs/kj2:
The infinite sum in (5.33) can then be controlled by
1X
kD1
. 1/k j.vs/ kj2  <.B 1.vs/ k; .vs/ k 1/  C1jjC1 X
k<0
j.vs/kj2;
with C1 a universal constant. As for the last term of the left-hand side of (5.33),
we write
..X?/vs; V vs/ D .. iB1 C iB 1/vs; V vs/
D .B1vs   B 1vs; iV vs/
D
X
k<0
k.B1.vs/k 1   B 1.vs/kC1; .vs/k/
j..X?/vs; V vs/j  C2jjC1
X
k<0
jkjj.vs/kj2;
where C2 is a universal constant. Lemma 5.9 follows upon taking C D C1 C
C2. 
Conclusion: Dealing with the Glancing
Consider a function  2 C1.M/ such that it coincides with M 3 x 7!
d.x; @M/ in a neighbourhood of @M and such that   0 and @M D  1.0/.
Clearly r.x/ D  .x/ for x 2 @M . Using , we extend  to the interior of M
as .x/ D  r.x/ for x 2M . We let .x; v/ WD hv; .x/i and
T WD V./X C X?:
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 39
Note that T is now defined on all SM and agrees with the vector field T defined
previously on @SM . In fact, T and V are tangent to every @SM" D f.x; v/ 2
SM W x 2  1."/g, where M" D  1. 1; ". The next lemma for fi is the key
input to deal with the glancing; cf. [38, lemma 4.1.3], [39, lemma 3.2.3], and [8,
lemma 5.1].
LEMMA 5.10. The functions V fi and T fi are bounded on SM n @0SM .
To substantiate the previous claim that the behaviour of u D uf is the same
as that of fi we proceed as follows. We consider a smooth integrating factor
R W SM ! GL.n;C/ such that XR C R D 0. These always exist for any
nontrapping manifold with strictly convex boundary. A simple calculation shows
that we may write u in terms of R as
u.x; v/ D R.x; v/
Z fi.x;v/
0
.R 1f /.'t .x; v//dt for .x; v/ 2 SM;
where 't is the geodesic flow of .M; g/. Thus directly from Lemma 5.10 we obtain
the following:
LEMMA 5.11. The functions V u and T u are bounded on SM n @0SM .
Next we note that all the previous work that we have done assuming u smooth
may be summarized as follows:
THEOREM 5.12. Let .M; g/ be a simple Riemannian surface with boundary and
 a smooth, skew-Hermitian matrix field in M . Then for any f 2 C1.M/, we
have the following stability estimate:
kf kL2.M;Cn/  C1.1C kkC1/eC2kkC1kvkH1.@SM;Cn/;
where v is any smooth solution of Xv Cv D  f .
PROOF OF THEOREM 5.3 IN FULL GENERALITY. LetM" for small "be the sur-
face considered above. We let u W SM ! Cn be the unique solution to the problem
XuCu D  f .SM/; uj@ SM D 0:
The function v WD ujSM" is smooth in SM" and solves Xv C v D  f since u
does. Hence we may apply Theorem 5.12 in M" to obtain
kf kL2.M";Cn/  C1.1C kkC1/eC2kkC1kukH1.@SM";Cn/;
where we might as well use the constants for M that bound those for M". We now
let "! 0; we clearly have
kf kL2.M";Cn/ ! kf kL2.M;Cn/;
and using Lemma 5.11 we see that
kukH1.@SM";Cn/ ! kukH1.@SM;Cn/:
40 F. MONARD, R. NICKL, AND G. P. PATERNAIN
Since
u.x; v/ D
(
I.f /.x; v/; .x; v/ 2 @CSM;
0; .x; v/ 2 @ SM;
the theorem is proved. 
5.4 Consistency of the Posterior Mean: Proof of Theorem 3.2
We assume ff2 D 1, the general case 0 < ff2 < 1 requires only notational
changes.
The overall strategy we pursue here, which has also been used in some form
in [29–31,44], is to show first that the Bayesian algorithm recovers the ‘regression
function’ C consistently in a natural statistical distance function, and to combine
this with quantitative stability estimates for the inverse map C 7!  in appropri-
ate metrics. This exploits crucially that the estimated Bayesian regression outputs
lie in the (nonlinearly constrained) range of the forward map C, so that the sta-
bility estimate applies to them. To make this approach work with ‘unbounded’
Gaussian priors is challenging, and our proofs proceed as follows: We first es-
tablish the posterior contraction Theorem 5.13 under general conditions, borrow-
ing from Bayesian nonparametric theory (e.g., [16, theorem 8.19] or [18, theorem
7.3.3]), slightly strengthening the usual statement of such theorems to give explicit
exponential bounds for the convergence rate to 0 of certain posterior probabili-
ties. Since our regression functions C take values in SO.n/, they are uniformly
bounded and the usual Hellinger distance occurring in such contraction theorems
is then Lipschitz-equivalent to the standard L2-distance (see Lemma 5.14). Then
Lemma 5.16 uses results of [26] to show that the key small ball condition in Theo-
rem 5.13 can be verified for the Gaussian priors from Condition 3.1 even after they
have been shrunk towards 0, if the true matrix field 0 belongs to the RKHS H.
Next, Lemma 5.17 exploits fine properties [3,15] of infinite-dimensional Gauss-
ian measures to show that such ‘shrunk’ priors charge ‘sufficiently regular’ matrix
fields (effectively C  -balls) with probability close enough to 1 that the posterior
distributions inherit these regularity properties. This is crucial to apply the ‘for-
ward’ estimate Theorem 2.2 and the ‘stability’ estimate (2.4) in the proof of The-
orem 5.19—effectively the specific structure of our inverse problem enters only in
this theorem and only through these two estimates.
Finally, the exponential convergence to 0 of the order e .CC3/N
2
N obtained in
Theorem 5.19 permits a ‘quantitative uniform integrability argument’ in Section
5.4 to deduce convergence of the whole posterior (Bochner-) mean towards the
true matrix field 0.
Let us mention that in the recent contributions [1, 19] (written after the first
version of this manuscript was completed), the general proof template developed
here has already been used effectively in two different nonlinear inverse problems
(arising with elliptic PDEs); see also Remark 3.6.
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 41
A General Contraction Theorem
Consider a collection P of probability density functions on some measurable
space .X ;A/ with respect to a dominating measure ; specifically in our measure-
ment model (1.3) we take
P D

p D
dP 1
d
W  2 C.M; so.n//
ff
; X D Rnn  @CSM;
where X is equipped with its natural product Borel-ff algebra A, where d D
dy  d with dy equal to Lebesgue measure on Rnn and  given in (1.5). By the
Gaussianity of the "1;j;k these probability densities are of the form
(5.43) p.y; .x; v// D
1
.2/n
2=2
exp

 1
2
X
1j;kn

yj;k   .C..x; v//j;k/
2ff
;
where .y; .x; v// 2 X . Since the map .; y; .x; v// 7! p.y; .x; v// is jointly
Borel-measurable from C.M/  X to R (using (2.3) and that point evaluation is
k  k1-continuous), the posterior distribution (3.2) exists by standard arguments (
[16], p.7) and has the desired form. In the proof of the following theorem, we show
in particular that the marginal density
R QN
iD1 p.Yi ; .Xi ; Vi //d./ is positive
on events of PN0-probability approaching 1, so that (3.2) is well-defined also in
the frequentist setting where DN  PN0 . We also define the Hellinger distance h
on such densities by
(5.44) h2.p; p/ D
Z
X
 p
p  pp
2
d; ; 2 C.M; so.n//:
Denote by N.F; h; / the minimal number of Hellinger-balls of radius  required
to cover a set F of -densities on X . We then have the following:
THEOREM 5.13. Consider a prior for  arising from a sequence  D N of
Borel probability measures on F  C.M; so.n// and let .  j.Yi ; .Xi ; Vi //NiD1/
be the posterior distribution arising from i.i.d. observations .Yi ; .Xi ; Vi //NiD1j 
PN . Let 0 2 F , let N ! 0 be a sequence such that
p
NN !1 as N !1,
and define sets
(5.45)
BN D

 2 F W E10

log
p0
p
..Y; .X; V //

 2N ;
E10

log
p
p0
..Y; .X; V //
2
 2N
ff
:
Suppose for some constant C > 0 the prior  satisfies for all N large enough
(5.46) .BN /  e CN
2
N ;
and that for some sequence FN  F of approximating sets for which
(5.47) 
 F n FN   Le .2CC6/N2N for some 0 < L <1;
42 F. MONARD, R. NICKL, AND G. P. PATERNAIN
we have the complexity bound
(5.48) logN.FN ; h; N /  cN2N
for some fixed constant c > 0. Then for some large enough constant m D
m.C; c/ > 0
PN0


 FN \ f W h.p; p0/  mN gj.Yi ; .Xi ; Vi //NiD1
 1   e .CC3/N2N

!N!1 0:
(5.49)
PROOF. Recall from (1.4) that we write DN D .Yi ; .Xi ; Vi //NiD1. The proof
proceeds as in the proof of [18, theorems 7.3.1 and 7.3.3]: We first use [18, lemma
7.3.2] and the hypothesis (5.46) to deduce that the events
(5.50) AN D
Z
F
NY
iD1
p
p0
.Yi ; .Xi ; Vi // d./  e .2CC/N
2
N
ff
satisfy PN0.AN /! 1 asN !1. Moreover, using (5.48) and [18, theorem 7.1.4]
with choices "0 D m0N , any m0 < m and logN."/ D cN2N constant in " > "0,
we deduce that for every k > 1 there exists m0; m large enough such that we can
find ‘tests’ (random indicator functions) N D N .DN / for which
(5.51)
PN0.N D 1/!N!1 0 and
sup
2FN Wh.p;p0 />mN
EN .1  N /  e kN
2
N :
Now let us write
xFN D FN \ fh.p; p0/  mN g
for the event whose posterior probability we want to bound. Then by (3.2) and as
N !1,
PN0
 

  xF cN jDN   e .CC3/N2N 
D PN0
 R
xF c
N
QN
iD1
p
p0
.Yi ; .Xi ; Vi //d./R
F
QN
iD1
p
p0
.Yi ; .Xi ; Vi //d./
 e .CC3/N2N ; N D 0; AN
!
C o.1/
 PN0
Z
xF c
N
NY
iD1
p
p0
.Yi ; .Xi ; Vi //d./.1  N /  e .2CC5/N
2
N

C o.1/:
By Markov’s inequality, decomposing
xF cN D FcN [ fh.p; p0/ > mN g
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 43
and using Fubini’s theorem as well as
(5.52) EN0
NY
iD1
p
p0
.Yi ; .Xi ; Vi //.1  N / D EN .1  N /  1;
we further bound the last probability as
e.2CC5/N
2
N
Z
xF c
N
EN .1  N /d./
 e.2CC5/N2N

2.FcN /C
Z
2FN Wh.p;p0 />mN
EN .1  N /d./

 2Le N2N C e.2CC5 k/N2N !N!1 0;
where we have used (5.47) and (5.51) with k and then m large enough. 
The ‘information-theoretic distance’ h arises naturally in such posterior con-
traction theorems; see [16]. The following lemma, which adapts a result due to
Birgé [2] to the setting of SO.n/-valued functions, shows that the Hellinger dis-
tance is Lipschitz-equivalent to the standard L2-metric
kC   CkL2 D
s X
1j;kn
kC;j;k   C;j;kk2L2 :
LEMMA 5.14. For 2 C.M; so.n//, letCW @CSM ! SO.n/ be its non-Abelian
X-ray transform. Then there exist positive constants c0 D c0.n/; c1 D c1.n/ such
that
(5.53)
1
c0
kC   Ck2L2  h2.p; p/  c1kC   Ck2L2 ;
8; 2 C.M; so.n//:
PROOF. Write
(5.54) .p; p/ 
Z
X
p
ppd D 1  
1
2
h2.p; p/
44 F. MONARD, R. NICKL, AND G. P. PATERNAIN
for the ‘Hellinger affinity’. By (5.43) and using the standard formula for the mo-
ment generating function ofN.0; 1/-variables with probability density ffi, the quan-
tity .p; p/ equals
D 1
.2/n
2=2
Z
X
exp

1
4
X
j;k
h
 yj;k   .C..x; v//j;k/2
  yj;k   .C..x; v//j;k/2iff
D
Z
@CSM
exp

 1
4
X
j;k

C 2..x; v//j;k C C 2..x; v//j;k
ff

Z
Rnn
e
1
2
P
j;k yj;k.C..x;v//j;kCC..x;v//j;k/
Y
j;k
ffi.yjk/dy d.x; v/
D
Z
@CSM
exp

 2
8
X
j;k

C 2..x; v//j;k C C 2..x; v//j;k
ff
 exp

1
8
X
j;k

C..x; v//j;k C C..x; v//j;k
2ff
d.x; v/
D
Z
@CSM
exp

 1
8
jC.x; v/   C.x; v/j2F
ff
d.x; v/:
By Jensen’s inequality the last integral is greater than or equal to expf kC  
Ck2L2=8g and using standard inequalities for 1 e ´; ´ > 0; the right-hand side of
(5.53) follows. Next we notice that since C.x; v/; C.x; v/ 2 SO.n/, their ma-
trix entries are all bounded by 1, and we hence have jC.x; v/   C.x; v/j2F =8 
B2 for some constant B D B.n/. We can thus proceed exactly as in the proof
of [2, prop. 1] (or see lemma 21 in [19]) to also deduce the left-hand side inequal-
ity in (5.53). 
Verification of the Prior Mass Condition
We now verify condition (5.46) in the last theorem for an explicit constantC > 0
and the Gaussian prior from Theorem 3.2. To do this we first show that one can
reduce to checking small ball conditions for k  kL2.M/-norms on the level of the
original matrix parameter .
LEMMA 5.15. For 0 2 C.M; so.n// and  > 0 define
BN ./ D f 2 C.M; so.n// W k  0kL2.M/  N =g;
and let BN ;; N ; be as in Theorem 5.13. Then for some  D .M; n/ large
enough, we have BN ./  BN , and thus in particular, for every N 2 N,
.BN /  .BN .//:
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 45
PROOF. From (1.3) with  D 0 and (5.43) we have
logp.Y1; .x; v//   logp0.Y1; .x; v//
D  
X
1j;kn

1
2
.C..x; v//j;k   C0..x; v//j;k/2
C "1;j;k.C..x; v//j;k   C0..x; v//j;k/

:
Therefore, since E1" "1;j;k D 0 and  is the unit volume measure on @CSM ,
E10

log
p0
p
..Y; .X; V //

D 1
2
kC   C0k2L2.@CSM/
 C
2
1
2
k  0k2L2.M/;
where we have used the forward estimate (2.2) with Lipschitz constant C1 D
C1.M; n/. Thus if   2=C 21 the first inequality defining BN is verified for
 2 BN ./. To verify the second, note that all C.x; v/ 2 SO.n/ are bounded in
the kkL1.@CSM/-norm by some fixed constant B D B.n/. Thus
E10

log
p
p0
.Y; .X; V //
2
 2E1
X
j;k
1
2
.C..X; V //j;k   C0..X; V //j;k/2
2
C 2E1E1"
X
j;k
"j;k.C..X; V //j;k   C0..X; V //j;k/
2
 c0.B; n/kC   C0k2L2  c.n/C1k  0k2L2
for some constant c.n/ > 0, where we have also used that "j;k i:i:d: N.0; 1/
implies, for .x; v/ 2 @CSM fixed,X
j;k
"j;k
 
C..x; v//j;k   C0..x; v//j;k
  N.0; jC.x; v/   C0.x; v/j2F /;
and again (2.2), so that the overall result follows from the appropriate choice of
 > 0. 
We now turn to lower bound the small ball probabilities.BN .// for the prior
 featured in Theorem 3.2 where for the given  we will choose
(5.55) N D N =.2C2/ so that
p
NN D N 1=.2C2/:
Note that
p
NN precisely equals the rescaling of the prior in (3.4). Let us recall
the base RKHS H from Condition 3.1.
46 F. MONARD, R. NICKL, AND G. P. PATERNAIN
LEMMA 5.16. Let  D xnjD1B be the prior for  from Theorem 3.2 with  >
 C 1,  > 0, assume 0 2 H, and choose N as in (5.55). Let BN ./ be as in
Lemma 5.15. Then for every  > 0 there exists a constant C 0 D C 0.; ; k0kH;
n;M/ such that for every N 2 N,
.BN .//  exp
 C 0N2N 	:
In particular, for BN as in (5.45) in Theorem 5.13, there exists a finite constant
C D C.; k0kH; n;M/ > 0
such that for every N 2 N,
(5.56) .BN /  exp
 CN2N 	:
PROOF. Since k 0kL2.M/  xnmaxj kBj  B0;j kL2.M/, to prove the first
inequality it suffices to lower bound, by independence of the Bj ’s,
xnY
jD1
B
 
B W kB   B0;j kL2.M/  N =.xn/

; xn D dim.so.n//:
The sets fb W kbkL2.M/  cg, c > 0, are convex and symmetric, hence by [18,
cor. 2.6.18] we have for every j fixed,
B.kB   B0;j kL2.M/  N =.xn//
 e kB0;j k2RKHS.B/=2B.kBkL2.M/  N =.xn//
D e N2N kB0;j k2H=2B.kBkL2.M/  N =.xn//
where we have used that
kB0;j k2RKHS.B/ D N
2
N kB0;j k2H <1
in view of (3.4), (5.55), (and where we refer to [18, exer. 2.6.5] or [16, lemma I.16]
for standard preservation properties of RKHS under linear transformations).
We next bound the centred probability that by (3.4), (5.55) equals
B.kBkL2.M/  N =.xn// D 0.kf 0kL2.M/ 
p
N2N =.xn//:
By Condition 3.1 the RKHS of the Gaussian law of f 0 in C.M/ is continuously
imbedded into H.M/. The unit ball U of this space satisfies the bound
(5.57) logN.U; k  kL2.M/; /  .A=/2=; 0 <  < A for some A > 0;
for its L2.M/-covering numbers: Indeed, since the simple surface M is diffeo-
morphic to a disk, we can extend all functions f in H.M/ to elements fe of
the Sobolev space H.I2/ on the 2-torus I2 D .0; 12 ff M , with Sobolev-norm
increased by at most a fixed multiplicative constant (chap. 4 in [41]). An appro-
priate bound for the L2.I2/-covering numbers of ffe W f 2 U g is then provided
in [18, (4.184)], which in turn (since kf   f 0kL2.M/  kfe   f 0ekL2.I2/ for all
f; f 0 2 L2.M/) also bounds the L2.M/-covering numbers of U as required.
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 47
To proceed, we can now use (5.55) and [26, theorem 1.2] (with the value of 
there equal to our 2=) to lower bound the last small ball probability by
(5.58)
exp
 cpN2N =.xn/  4=2 .2=/ 	  exp c0N2N 	;
noting
p
N2N D N . 1/=.2C2/
for constants c D c./; c0 D c0.; n; / and since  > 1. Combining what
precedes proves the first inequality of the lemma with
(5.59) C 0 D 1
2
xnX
jD1
kB0;j k2H C xnc0:
The second inequality (5.56) now follows from the first and Lemma 5.15. 
We note that the proof in fact shows that the constant C depends only on upper
bounds for k0kH.
Excess Mass and Complexity Conditions
Having determined the constant C in (5.46) for the Gaussian prior in Theo-
rem 3.2, we now turn to verifying the remaining conditions (5.47) and (5.48) in
Theorem 5.13 for a suitable choice of FN that will provide sufficient regularity
of the posterior distribution to combine it with our stability estimates for the map
 7! C.
LEMMA 5.17. Let be the prior from Theorem 3.2 with  > C1,  > 0, let N
be as in (5.55) and assume N2N  1. Form > 0 define subsets of C.M; so.n// as
FN D

 W  D 1 C2; k1kL2  N ; k2kH  m; kkC  m
	
(a) Then for every K > 0 we can choose m large enough such that
.FN /  1   e KN2N :
(b) Moreover, for some c D .m; ; n; vol.M// we have
logN.FN ; h; N /  cN2N :
PROOF.
(a) Recalling (3.4), (5.55), we can identify a prior draw  with the vector field
.B1; : : : ; Bxn/ D
1p
NN
 
f 01; : : : ; f
0
xn

; f 0j i:i:d: 0:
We denote by 0xn the product measure describing the law of the centred Gaussian
random variable .f 01; : : : ; f
0
xn/ in the Banach space xnjD1C.M/.
48 F. MONARD, R. NICKL, AND G. P. PATERNAIN
Write next FN D FN;1 \ FN;2 where, with f 0i; corresponding to i , i D 1; 2;
FN;1 D
 
f 0j D f 01;j C f 02;j
xn
jD1
W
xnX
jD1
kf 01;j k2L2  N4N ;
xnX
jD1
kf 02;j k2H.M/  m2N2N
ff
;
FN;2 D

.f 01; : : : ; f
0
xn/ W
xnX
jD1
kf 0j kC.M/  m
p
NN
ff
;
so that it suffices to bound the prior probabilities of the complements ofFN;1;FN;2.
We first turn to FN;2. By Condition 3.1 the vector field .f 01; : : : ; f 0xn/ defines a
Gaussian Borel random variable in a separable linear subspace S of xnjD1C  .M/.
By the Hahn-Banach theorem its xnjD1C  .M/-norm can then be represented as a
countable supremum
 f 01; : : : ; f 0xnxn
jD1
C.M/
D sup
t2T
t f 01; : : : ; f 0xn
of bounded linear functionals T D .tm W m 2 N/ defined on .S; kkxn
jD1
C.M//.
We then apply Fernique’s theorem [15], concretely [18, theorem 2.1.20], to the
centred Gaussian process .X.t/ WD t .f 01; : : : ; f 0xn/ W t 2 T / to deduce that for some
fixed constant D > 0,
E
xnX
jD1
f 0j C.M/  D <1;
and then also, for m D m.D/ large enough and since N2N  1,

 FcN;2  xxn xnX
jD1
kf 0j kC.M/  E
xnX
jD1
kf 0j kC.M/  m
p
NN =2

 2e km2N2N
for k a fixed constant, which can be made less than e KN
2
N =2 for anyK provided
m D m.K; k;D/ is chosen large enough.
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 49
It remains to show that.FN;1/  1  12 expf KN2N g form large enough. Us-
ing the continuous imbedding H  H.M/ with imbedding constant c0 (cf. Con-
dition 3.1), it suffices to lower bound
0xn
 
f 0j D f 01;j C f 02;j
xn
jD1
W
xnX
jD1
f 01;j2L2  N4N ; xnX
jD1
kf 02;j k2H
1=2
 m
c0
p
NN

D 0xn
  xAN CmNOH
where OH is the unit ball in xnjD1H and where we define
xAN 

! 2 xnjD1C.M/ W k!kL2 
p
N2N
	
; mN 
m
p
NN
c0
:
By Borell’s [3] isoperimetric inequality (see [18, theorem 2.6.12]) the last proba-
bility is bounded below by
(5.60) 
 
 1
 
0xn.
xAN /
CmN 
where  D Pr.Z  / is the cumulative distribution function of a N.0; 1/ random
variable Z. By the same arguments as those leading to (5.58) above, we have
0xn
  xAN   exp c22N2N 	 for c2 D c2.n; / > 0;
and using the basic inequality  1.u/   p2 log  u, 0 < u < 1 (see [16, lemma
K.6]) and monotonicity of  we can further lower bound (5.60) by


  c2
p
2C m
c0
p
NN

:
Now given K, define
m0N    1

exp
  KN2N =2;
which by the previous inequality for  1 can be made to be less than or equal to
.m
c0
  c2
p
2/
p
NN whenever m D m.K; c2; c0/ is large enough. Conclude that
the penultimate display is lower bounded by

   1exp. KN2N /=2 D 1    1exp  KN2N =2
D 1   1
2
exp
 KN2N 	;
completing the proof of Part (a).
(b) To prove Part (b), note first that to construct a N -covering ofFN in kkL2.M/-
distance it suffices, by definition ofFN , to construct such a covering for aH.M/-
ball of radius m, so that (5.57) and the definition of N give (with A0 > 0)
(5.61)
logN.FN ; kkL2.M/; N /  .A0=N /2=  bN2N
for some b D b.m; ; xn/ > 0:
50 F. MONARD, R. NICKL, AND G. P. PATERNAIN
Lemma 5.14 and (2.2) imply that such a covering induces a .C1
p
c1/N -covering
ofFN in the Hellinger distance h of log-cardinality at most bN2N . Since kkL2.M/
is a norm and hence homogeneous, we can increase the constant from b to c D
c.b; c1; C1; / in (5.61) and obtain a n-covering for h. The desired inequality in
Part (b) follows. 
Remark 5.18. We note that the introduction of the set FN;1 and the use of Borell’s
inequality in the previous lemma can be avoided if one wishes to prove Theorem
3.2 only for any  > 0 (in this case a minor adaptation of Theorem 5.13 and
of (5.61) can be shown to give a slightly worse rate 0N D N =.2C2/ in (5.62)
below). We give this argument however to obtain our sharper bound for  in (5.66).
Final Contraction Theorem
We now put everything together to establish a posterior contraction theorem for
 and subsequently deduce Theorem 3.2 from it.
THEOREM 5.19. Under the hypotheses of Theorem 3.2, with  >  C 1,  > 0,
N D N =.2C2/, and C from (5.56), we have for all m0 large enough that
PN0


 
 W kC   C0kL2.@CSM/  m0N ;
kkC.M/  m0jDN /  1   e .CC3/N
2
N

! 1
(5.62)
as N ! 1. Moreover, if  > 2, then we have for every integer x such that
1 < x <  and all m00 large enough,
PN0
 

 
 W k  0kL2.M/  m00.
x 1/= x
N jDN
  e .CC3/N2N !N!1 0:
Remark 5.20. The constraint  > 2 in the second limit in Theorem 5.19 is only
required to allow space for an integer x 2 .1; / in the following proof, when
combining the interpolation inequality (5.64) with Theorem 2.2 for k D x 2 N. If
a version of Theorem 2.2 were established for noninteger k, then  > 1 and real
x 2 .1; / would be permitted in Theorem 5.19 (and then also in Theorem 3.2).
PROOF. From Lemmata 5.16 and 5.17 withK D 2C C 6 and Theorem 5.13 we
deduce for m D m.C/ large enough, and as N !1
PN0


 
 W f W h.p; p0/  mN g \ fkkC.M/  mgjDN

 1   e .CC3/N2N

! 0:
Applying Lemma 5.14 gives the first limit (5.62) with m0 D .1Cpc0/m.
To prove the second limit we will apply the stability estimate Theorem 2.1 in the
form (2.4) with  D 0. By hypothesis we have k0kC1.M/ . k0kC.M/ <
1; as a consequence for all contained in the event in (5.62) with  > 2, the con-
stants c.;0/ from (2.1) are uniformly bounded by a fixed constant that depends
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 51
on m0;M; k0kC1.M/ and hence for those ’s
(5.63) k  0kL2.M/  D.k0kC1.M/;M;m0/kC   C0kH1.@CSM/:
To proceed we will need a standard interpolation result for Sobolev spaces on the
manifold @CSM to the effect that
(5.64) kW kH1.@CSM/ . kW k
.k 1/=k
L2.@CSM/
kW k1=k
Hk.@CSM/
for all W 2 Hk.@CSM/ and any k > 1. [For real-valued functions this can be
proved using standard arguments from chapter 4 in [41], and these results extend to
matrix fields in a straightforward way.] Moreover, we will use the basic inequality
(5.65) kk
H x.M/
. kk
C x.M/
 kkC.M/
for all 2 C  .M/. Now Theorem 2.2 implies that for all’s in the event in (5.62)
the corresponding kCkH x.@CSM/’s are uniformly bounded by a fixed constant
that depends on m0;M; ; x only. Likewise
kC0kH x.@CSM/  kC0kC.@CSM/ . .1C k0kC / <1
in view of Theorem 2.2 and since 0 2 C  for  > x by hypothesis. Hence for
such ’s the combination of (5.63) and (5.64) with W D C   C0 , k D x gives
k  0kL2.M/ . kC   C0k.
x 1/= x
L2.@CSM/
kC   C0k1=
x
H x.@CSM/
. 
. x 1/= x
N :
The second conclusion of Theorem 5.19 now follows from the preceding inequali-
ties and (5.62).
Completion of the Proof of Theorem 3.2
The last step is to show that the posterior contraction rate in the second limit
of Theorem 5.19 carries over to the posterior mean EjDN . For any integerx 2 .1; / and every
(5.66) 0 <  <

2 C 2 
x   1
x ;
we have as N !1
N WD m00.
x 1/= x
N ' N
  
2C2
x 1
x D o.N /:
Then by the inequalities of Jensen and Cauchy-Schwarz
kEjDN   0kL2.M/
 Ek  0kL2.M/jDN 
 N CEk  0kL2.M/1fk  0kL2.M/  N gjDN 
 N C Ek  0k2L2.M/jDN 1=2.k  0kL2.M/  N jDN /1=2;
52 F. MONARD, R. NICKL, AND G. P. PATERNAIN
and it suffices to show that the second summand is stochastically O.N / as N !
1.
Arguing as in the proof of Theorem 5.13 and using Lemma 5.16 implies that the
sets AN from (5.50) with C from (5.56) satisfy PN0.AN /! 1 as N !1. Now
Theorem 5.19, (3.2), and Markov’s inequality imply
PN0
 
Ek  0k2L2.M/jDN  .k  0kL2.M/  N jDN / > 2N

 PN0
 
Ek  0k2L2.M/jDN e .CC3/N
2
N > 2N
C o.1/
 PN0
 
e .CC3/N
2
N
R k  0k2L2.M/QNiD1 pp0 .Yi ; .Xi ; Vi //d./R QN
iD1
p
p0
.Yi ; .Xi ; Vi //d./
> 2N ; AN
!
C o.1/
 e N2N  2N EN0
Z
k  0k2L2.M/
NY
iD1
p
p0
.Yi ; .Xi ; Vi //d./
 e N2N  2N
Z
k  0k2L2.M/ d./ . e N
2
N  2N !N!1 0
where we have also used Fubini’s theorem, (5.52), and that the Gaussian measure
 is supported in L2.M/ and hence integrates kk2
L2
to a finite constant (see,
e.g., [18, exer. 2.1.5]). 
Acknowledgment. We would like to thank the referee for helpful remarks and
suggestions. We are further very grateful to Bill Lionheart for having introduced
us to polarimetric neutron tomography and its connection to the non-Abelian X-
ray transform. FM was supported by National Science Foundation Grant DMS-
1814104 and a UC Hellman Fellowship. RN was supported by the European Re-
search Council under ERC Grant No. 647812 (UQMSI). GPP thanks the University
of California at Santa Cruz and the University of Washington for hospitality while
this work was in progress. GPP was supported by the Leverhulme trust and EPSRC
Grant EP/R001898/1.
Bibliography
[1] Abraham, K.; Nickl, R. On statistical Calderón problems. Math. Stat. Learn. 2 (2019), no. 2,
165–216. doi:10.4171/MSL/14
[2] Birgé, L. Model selection in Gaussian regression with random design. Bernoulli 10 (2004),
no. 6, 1039–1051.
[3] Borell, C. The Brunn-Minkowski inequality in Gauss space. Invent. Math. 30 (1975), no. 2,
207–216. doi:10.1007/BF01425510
[4] Briol, F.-X.; Oates, C. J.; Girolami, M.; Osborne, M. A.; Sejdinovic, D. Probabilistic integra-
tion: A role in statistical computation? Statist. Sci. 34 (2019), no. 1, 1–22. doi:10.1214/18-
STS660
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 53
[5] Castillo, I.; Nickl, R. Nonparametric Bernstein-von Mises theorems in Gaussian white noise.
Ann. Statist. 41 (2013), no. 4, 1999–2028. doi:10.1214/13-AOS1133
[6] Castillo, I.; Nickl, R. On the Bernstein-von Mises phenomenon for nonparametric Bayes pro-
cedures. Ann. Statist. 42 (2014), no. 5, 1941–1969. doi:10.1214/14-AOS1246
[7] Cotter, S. L.; Roberts, G. O.; Stuart, A. M.; White, D. MCMC methods for functions: modifying
old algorithms to make them faster. Statist. Sci. 28 (2013), no. 3, 424–446. doi:10.1214/13-
STS421
[8] Dairbekov, N. S.; Paternain, G. P.; Stefanov, P; Uhlmann, G. The boundary rigidity
problem in the presence of a magnetic field. Adv. Math. 216 (2007), no. 2, 535–609.
doi:10.1016/j.aim.2007.05.014
[9] Dashti, M.; Law, K. J. H.; Stuart, A. M.; Voss, J. MAP estimators and their consistency in
Bayesian nonparametric inverse problems. Inverse Problems 29 (2013), no. 9, 095017, 27 pp.
doi:10.1088/0266-5611/29/9/095017
[10] Dashti, M.; Stuart, A. M. The Bayesian approach to inverse problems. Handbook of uncertainty
quantification. Vol. 1, 2, 3, 311–428. Springer, Cham, 2017.
[11] Dawson, M.; Manke, I.: Kardjilov, N.: Hilger, A.; Strobl, M.; Banhart, J. Imaging with polarized
neutrons. New Journal of Physics 11 (2009), 043013. doi:10.1088/1367-2630/11/4/043013
[12] Desai, N. M.; Lionheart, W. R. B.; Sales, M.; Strobl, S.; Schmidt, S. Polarimetric neutron to-
mography of magnetic fields: uniqueness of solution and reconstruction. Inverse Problems 36
(2020), no. 4, 045001. doi:10.1088/1361-6420/ab44e0
[13] Diaconis, P. Bayesian numerical analysis. Statistical decision theory and related topics, IV,
Vol. 1 (West Lafayette, Ind., 1986), 163–175. Springer, New York, 1988.
[14] Eskin, G. On non-abelian Radon transform. Russ. J. Math. Phys. 11 (2004), no. 4, 391–408.
[15] Fernique, X. Regularité des trajectoires des fonctions aléatoires gaussiennes. (French) École
d’Été de Probabilités de Saint-Flour, IV-1974, 1–96. Lecture Notes in Mathematics, 480.
Springer, Berlin, 1975.
[16] Ghosal, S.; van der Vaart, A. Fundamentals of nonparametric Bayesian inference. Cambridge
Series in Statistical and Probabilistic Mathematics, 44. Cambridge University Press, Cam-
bridge, 2017. doi:10.1017/9781139029834
[17] Giné, E.; Nickl, R. Rates of contraction for posterior distributions in Lr -metrics, 1  r  1.
Ann. Statist. 39 (2011), no. 6, 2883–2911. doi:10.1214/11-AOS924
[18] Giné, E.; Nickl, R. Mathematical foundations of infinite-dimensional statistical models. Cam-
bridge Series in Statistical and Probabilistic Mathematics, 40. Cambridge University Press,
New York, 2016. doi:10.1017/CBO9781107337862
[19] Giordano, M.; Nickl, R. Consistency of Bayesian inference with Gaussian process priors in
an elliptic inverse problem. Inverse Problems 36 (2020), no. 4, 085001. doi:10.1088/1361-
6420/ab7d2a
[20] Guillemin, V.; Kazhdan, D. Some inverse spectral results for negatively curved 2-manifolds.
Topology 19 (1980), no. 3, 301–312. doi:10.1016/0040-9383(80)90015-4
[21] Hairer, M.; Stuart, A. M.; Vollmer, S. J. Spectral gaps for a Metropolis-Hastings algorithm in
infinite dimensions. Ann. Appl. Probab. 24 (2014), no. 6, 2455–2490. doi:10.1214/13-AAP982
[22] Hilger, A.; Manke, I.; Kardiljov, N.; Osernberg, M.; Markotter, H.; Banhart, J. Tensorial neutron
tomography of threedimensional magnetic vector fields in bulk materials. Nature Communica-
tions 9 (2018), 1–7. doi:10.1038/s41467-018-06593-4
[23] Ilmavirta, J.; Monard, F. Integral geometry on manifolds with boundary and applications. The
radon transform: the first 100 years and beyond, 43–114. de Gruyter, Berlin-Boston, 2019.
doi:10.1515/9783110560855
[24] Kaltenbacher, B.; Neubauer, A.; Scherzer, O. Iterative regularization methods for nonlinear
ill-posed problems. Radon Series on Computational and Applied Mathematics, 6. de Gruyter,
Berlin, 2008. doi:10.1515/9783110208276
54 F. MONARD, R. NICKL, AND G. P. PATERNAIN
[25] Kardjilov, N.; Manke, I.; Strobl, M.; Hilger, A.; Treimer, W.; Meissner, M.; Krist T.; Banhart, J.
Three-dimensional imaging of magnetic field with polarized neutrons. Nature Physics 4 (2008)
399–403. doi:10.1038/nphys912
[26] Li, W. V.; Linde, W. Approximation, metric entropy and small ball estimates for Gaussian
measures. Ann. Probab. 27 (1999), no. 3, 1556–1578. doi:10.1214/aop/1022677459
[27] Merry, W; Paternain, G. P. Lecture notes on inverse problems and dynamics. Unpublished notes,
2011.
[28] Monard, F.; Nickl, R.; Paternain, G. P. Efficient nonparametric Bayesian inference for X-ray
transforms. Ann. Statist. 47 (2019), no. 2, 1113–1147. doi:10.1214/18-AOS1708
[29] Nickl, R. Bernstein-von Mises theorems for statistical inverse problems I: Schrödinger equa-
tion. J. Eur. Math. Soc. (JEMS) 22 (2020), no. 8, 2697–2750. doi:10.4171/JEMS/975
[30] Nickl, R.; Söhl, J. Nonparametric Bayesian posterior contraction rates for discretely observed
scalar diffusions. Ann. Statist. 45 (2017), no. 4, 1664–1693. doi:10.1214/16-AOS1504
[31] Nickl, R.; van de Geer, S.; Wang, S. Convergence rates for penalized least squares estimators
in PDE constrained regression problems. SIAM/ASA J. Uncertain. Quantif. 8 (2020), no. 1,
374–413. doi:10.1137/18M1236137
[32] Novikov, R. G. Non-abelian Radon transform and its applications. The radon transform: the first
100 years and beyond, 115–128. de Gruyter, Berlin-Boston, 2019. doi:10.1515/9783110560855
[33] Novikov, R. G. On determination of a gauge field on Rd from its non-abelian Radon
transform along oriented straight lines. J. Inst. Math. Jussieu 1 (2002), no. 4, 559–629.
doi:10.1017/S1474748002000166
[34] Paternain, G. P.; Salo, M.; Uhlmann, G. The attenuated ray transform for connections and Higgs
fields. Geom. Funct. Anal. 22 (2012), no. 5, 1460–1489.
[35] Paternain, G. P.; Salo, M. The non-Abelian X-ray transform on surfaces. Preprint, 2020.
arXiv:2006.02257 [math.DG]
[36] Reiß, M. Asymptotic equivalence for nonparametric regression with multivariate and random
design. Ann. Statist. 36 (2008), no. 4, 1957–1982. doi:10.1214/07-AOS525
[37] Sales, M.; Strobl, M.; Shinohara, T.; Tremsin, A.; Kuhn, L. T.; Lionheart, W. R.; Desai, N. M.;
Dahl, A. B.; Schmidt, S. Three dimensional polarimetric neutron tomography of magnetic
fields. Scientific Reports 8 (2018), no. 1, 2214.
[38] Sharafutdinov, V. A. Integral geometry of tensor fields. Inverse and ill-posed problems series.
VSP, Utrecht, 1994. doi:10.1515/9783110900095
[39] Sharafutdinov, V. A. Ray transform on Riemannian manifolds. Eight Lectures on Integral Ge-
ometry. Lecture notes, 1999.
[40] Stuart, A. M. Inverse problems: a Bayesian perspective. Acta Numer. 19 (2010) 451–559.
doi:10.1017/S0962492910000061
[41] Taylor, M. E. Partial differential equations I. Basic theory. Second edition. Applied Mathemat-
ical Sciences, 115. Springer, New York, 2011. doi:10.1007/978-1-4419-7055-8
[42] van der Vaart, A. W.; van Zanten, H. Bayesian inference with rescaled Gaussian process priors.
Electron. J. Stat. 1 (2007), 433–448. doi:10.1214/07-EJS098
[43] van der Vaart, A. W.; van Zanten, J. H. Rates of contraction of posterior distribu-
tions based on Gaussian process priors. Ann. Statist. 36 (2008), no. 3, 1435–1463.
doi:10.1214/009053607000000613
[44] Vollmer, S. J. Posterior consistency for Bayesian inverse problems through stability and
regression results. Inverse Problems 29 (2013), no. 12, 125011, 32 pp. doi:10.1088/0266-
5611/29/12/125011
CONSISTENT INVERSION OF NOISY NON-ABELIAN X-RAY TRANSFORMS 55
FRANÇOIS MONARD
Department of Mathematics
University of California, Santa Cruz
Santa Cruz, CA 95064
USA
E-mail: fmonard@ucsc.edu
GABRIEL P. PATERNAIN
Department of Pure Mathematics
and Mathematical Statistics
University of Cambridge
Cambridge CB3 0WB
UNITED KINGDOM
E-mail: g.p.paternain@
dpmms.cam.ac.uk
RICHARD NICKL
Department of Pure Mathematics
and Mathematical Statistics
University of Cambridge
Cambridge CB3 0WB
UNITED KINGDOM
E-mail: r.nickl@
statslab.cam.ac.uk
Received May 2019.
Revised May 2020.