A Statistical Perspective on Inverse and Inverse Regression
Problems
Debashis Chatterjeea , Sourabh Bhattacharyaa,∗

arXiv:1707.06852v1 [stat.ME] 21 Jul 2017

a Interdisciplinary

Statistical Research Unit, Indian Statistical Institute,
203 B. T. Road, Kolkata - 700108, India

Abstract
Inverse problems, where in broad sense the task is to learn from the noisy response
about some unknown function, usually represented as the argument of some known
functional form, has received wide attention in the general scientific disciplines. However, in mainstream statistics such inverse problem paradigm does not seem to be as
popular. In this article we provide a brief overview of such problems from a statistical,
particularly Bayesian, perspective.
We also compare and contrast the above class of problems with the perhaps more
statistically familiar inverse regression problems, arguing that this class of problems
contains the traditional class of inverse problems. In course of our review we point out
that the statistical literature is very scarce with respect to both the inverse paradigms,
and substantial research work is still necessary to develop the fields.
Keywords: Bayesian analysis, Inverse problems, Inverse regression problems,
Regularization, Reproducing Kernel Hilbert Space (RKHS), Palaeoclimate
reconstruction

1. Introduction
The similarities and dissimilarities between inverse problems and the more traditional forward problems are usually not clearly explained in the literature, and often
“ill-posed” is the term used to loosely characterize inverse problems. We point out that
these two problems may have the same goal or different goal, while both consider the
same model given the data. We first elucidate using the traditional case of deterministic
differential equations, that the goals of the two problems may be the same. Consider a
dynamical system
dxt
= G(t, xt , θ),
(1.1)
dt
where G is a known function and θ is a parameter. In the forward problem the goal
is to obtain the solution xt ≡ xt (θ), given θ and the initial conditions, whereas, in
∗ Corresponding

author.
Email addresses: debashis1chatterjee@gmail.com (Debashis Chatterjee),
bhsourabh@gmail.com, sourabh@isical.ac.in (Sourabh Bhattacharya)

Preprint submitted to XYZ

July 24, 2017

the inverse problem, the aim is to obtain θ given the solution process xt . Realistically,
the differential equation would be perturbed by noise, and so, one observes the data
y = (y1 , . . . , yT )T , where
yt = xt (θ) + t ,
(1.2)
for noise variables t having some suitable independent and identical (iid) error distribution q, which we assume to be known for simplicity of illustration. A typical
method of estimating θ, employed by the scientific community, is the method of calibration, where the solution of (1.1) would be obtained for each θ-value on a proposed
grid of plausible values, and a set ỹ(θ) = (ỹ1 (θ), . . . , ỹT (θ))T is generated from the
iid

model (1.2) for every such θ after simulating, for i = 1, . . . , T , ˜t ∼ q; then forming
ỹt (θ) = xt (θ) + ˜t , and finally reporting that value θ in the grid as an estimate of the
true values for which ky − ỹ(θ)k is minimized, given some distance measure k·k; maximization of the correlation between y and ỹ(θ) is also considered. In other words, the
calibration method makes use of the forward technique to estimate the desired quantities of the model. On the other hand, the inverse problem paradigm attempts to directly
estimate θ from the observed data y usually by minimizing some discrepancy measure
between y and x(θ), where x(θ) = (x1 (θ), . . . , xT (θ))T . Hence, from this perspective the goals of both forward and inverse approaches are the same, that is, estimation
of θ. However, the forward approach is well-posed, whereas, the inverse approach is
often ill-posed. To clarify, note that within a grid, there always exists some θ̂ that minimizes ky − ỹ(θ)k among all the grid-values. In this sense the forward problem may be
thought of as well-posed. However, direct minimization of the discrepancy between y
and x(θ) with respect to θ is usually difficult and for high-dimensional θ, the solution
to the minimization problem is usually not unique, and small perturbations of the data
causes large changes in the possible set of solutions, so that the inverse approach is
usually ill-posed. Of course, if the minimization is sought over a set of grid values of θ
only, then the inverse problem becomes well-posed.
From the statistical perspective, the unknown parameter θ of the model needs to be
learned, in either classical or Bayesian way, and hence, in this sense there is no real
distinction between forward and inverse problems. Indeed, statistically, since the data
are modeled conditionally on the parameters, all problems where learning the model
parameter given the data is the goal, are inverse problems. We remark that the literature
usually considers learning unknown functions from the data in the realm of inverse
problems, but a function is nothing but an infinite-dimensional parameter, which is a
very common learning problem in statistics.
We now explain when forward and inverse problems can differ in their aims, and are
significantly different even from the statistical perspective. To give an example, consider the palaeoclimate reconstruction problem discussed in Haslett et al. [18] where
the reconstruction of prehistoric climate at Glendalough in Ireland from fossil pollen
is of interest. The model is built on the realistic assumption that pollen abundance depends upon climate, not the other way around. The compositional pollen data with the
modern climates are available at many modern sites but the climate values associated
with the fossil pollen data are missing. The inverse nature of the problem is associated
with the fact that it is of interest to predict the fossil climate values, given the pollen
assemblages. The forward problem would result, if given the fossil climate values (if
2

known), the fossil pollen abundances (if unknown), were to be predicted.
Technically, given a data set y that depends upon covariates x, with a probability
distribution f (y|x, θ) where θ is the model parameter, we call the problem ‘inverse’
if it is of interest to predict the corresponding unknown x̃ given a new observed ỹ
(see Bhattacharya and Haslett [9]), after eliminating θ. On the other hand, the more
conventional forward problem considers the prediction of ỹ for given x̃ with the same
probability distribution, again, after eliminating the unknown parameter θ. This perspective clearly distinguishes the forward and inverse problems, as opposed to the other
parameter-learning perspective discussed above, which is much more widely considered in the literature. In fact, with respect to predicting unknown covariates from the
responses, mostly inverse linear regression, particularly in the classical set-up, has been
considered in the literature. To distinguish the traditional inverse problems from the
covariate-prediction perspective, we use the phrase ‘inverse regression’ to refer to the
latter. Other examples of inverse regression are given in Section 7.
Our discussion shows that statistically, there is nothing special about the existing
literature on inverse problems that considers estimation of unknown (perhaps, infinitedimensional) parameters, and the only class of problems that can be truly regarded
as inverse problems as distinguished from forward problems are those which consider
prediction of unknown covariates from the dependent response data. However, for the
sake of completeness, the traditional inverse problems related to learning of unknown
functions shall occupy a significant portion of our review.
The rest of the paper is structured as follows. In Section 2 we discuss the general
inverse model, providing several examples. In Section 3 we focus on linear inverse
problems, which constitute the most popular class of inverse problems, and review the
links between the Bayesian approach based on simple finite difference priors and the
deterministic Tikhonov regularization. Connections between Gaussian process based
Bayesian inverse problems and deterministic regularizations are reviewed in Section
4. In Section 5 we provide an overview of the connections between the Gaussian process based Bayesian approach and regularization using differential operators, which
generalizes the discussion of Section 3 on the connection between finite difference priors and the Tikhonov regularization. The Bayesian approach to inverse problems in
Hilbert spaces is discussed in Section 6. We then turn attention to inverse regression
problems, providing an overview of such problems and discussing the links with traditional inverse problems in Section 7. Finally, we make concluding remarks in Section
8.
2. Traditional inverse problem
Suppose that one is interested in learning about the function θ given the noisy observed responses y n = (y1 , . . . , yn )T , where the relationship between θ and y n is
governed by following equation (2.1) :
yi = G(xi , θ) + i ,

(2.1)

for i = 1, . . . , n, where xi are known covariates or design points, i are errors associated with the i-th observation and G is a forward operator defined appropriately, which
is usually allowed to be non-injective.
3

Note that since n = (1 , . . . , n )T is unknown, the noisy observation vector y n
itself may not be in the image set of G. If θ is a p-dimensional parameter, then there
will often be situations when the number of equations is smaller than the number of
unknowns, in the sense that p > n (see, for example, Dashti and Stuart [14]). Modern statistical research is increasingly coming across such inverse problems termed
as “ill-posed” which are not in the exact domain of statistical estimation procedures
(O’Sullivan [27]) where the maximum likelihood solution or classical least squares
may not be uniquely defined and with very bad perturbation sensitivity of the classical
solution. However, although such problematic issues are said to characterize inverse
problems, the problems in fact fall in the so-called “large p small n” paradigm and
has received wide attention in statistics; see, for example, Bühlmann and van de Geer
[10], Giraud [17]. A key concept involved in handling such problems is inclusion of
some appropriate penalty term in the discrepancy to be minimized with respect to θ.
Such regularization methods are initiated by Tikhonov [34] and Tikhonov and Arsenin
[35]. Under this method, usually a criterion of the following form is chosen for the
minimization purpose:
n

1X
2
[yi − G(xi , θ)] + λJ(θ), λ > 0.
n i=1

(2.2)

The functional J is chosen such that highly implausible or irregular values of θ has
large values (O’Sullivan [27]). Thus, depending on the problem at hand, J(θ) can be
used to induce “sparsity” in an appropriate sense so that the minimization problem may
be well-defined. We next present several examples of classical inverse problems based
on Aster et al. [3].
2.1. Examples of inverse problems
2.1.1. Vertical seismic profiling
In this scientific field, one wishes to learn about the vertical seismic velocity of the
material surrounding a borehole. A source generates downward-propagating seismic
wavefront at the surface, and in the borehole, a string of seismometers sense these seismic waves. The arrival times of the seismic wavefront at each instrument are measured
from the recorded seismograms. These times provide information on the seismic velocity for vertically traveling waves as a function of depth. The problem is nonlinear if
it is expressed in terms of seismic velocities. However, we can linearize this problem
via a simple change of variables, as follows. Letting z denote the depth, it is possible to
parameterize the seismic structure in terms of slowness, s(z), which is the reciprocal
of the velocity v(z). The observed travel time at depth z can then be expressed as:
Z z
Z ∞
t(z) =
s(u)du =
s(u)H(z − u)du,
(2.3)
0

0

where H is the Heaviside step function. The interest is to learn about s(z) given observed t(z). Theoretically, s(z) = dt(z)
dz , but in practice, simply differentiating the
observations need not lead to useful solutions because noise is generally present in the
observed times t(z), and naive differentiation may lead to unrealistic features of the
solution.
4

2.1.2. Estimation of buried line mass density from vertical gravity anomaly
Here the problem is to estimate an unknown buried line mass density m(x) from
data on vertical gravity anomaly, d(x), observed at some height, h. The mathematical
relationship between d(x) and m(x) is given by
Z ∞
h
d(x) =
3 m(u)du.
−∞ [(u − x)2 + h2 ] 2
As before, noise in the data renders the above linear inverse problem difficult. Variations of the above example has been considered in Aster et al. [3].
2.1.3. Estimation of incident light intensity from diffracted light intensity
Consider an experiment in which an angular distribution of illumination passes
through a thin slit and produces a diffraction pattern, for which the intensity is observed. The data, d(s), are measurements of diffracted light intensity as a function
of the outgoing angle −π/2 ≤ s ≤ π/2. The goal here is to obtain the intensity of
incident light on the slit, m(θ), as a function of the incoming angle −π/2 ≤ θ ≤ π/2,
using the following mathematical relationship:
Z

π/2

d(s) =

(cos(s) + cos(θ))

2

−π/2



sin(π (sin(s) + sin(θ)))
π (sin(s) + sin(θ))

2
m(θ)dθ.

2.1.4. Groundwater pollution source history reconstruction problem
Consider the problem of recovering the history of groundwater pollution at a source
site from later measurements of the contamination at downstream wells to which the
contaminant plume has been transported by advection and diffusion. The mathematical model for contamination transport is given by the following advection-diffusion
equation with respect to t and transported site x:
∂C
∂2C
∂C
=D 2 −ν
∂t
∂x
∂x
C(0, t) = Cin (t)
C(x, t) → 0 as x → ∞.
In the above, D is the diffusion coefficient, ν is the velocity of the groundwater flow,
and Cin (t) is the time history of contaminant injection at x = 0. The solution to the
above advection-diffusion equation is given by
Z

T

Cin (t)f (x, T − t)dt,

C(x, T ) =
0

where

"

(x − ν(T − t))
f (x, T − t) = p
exp
3
4D(T − t)
2 πD(T − t)
x

2

It is of interest to learn about Cin (t) from data observed on C(x, T ).

5

#
.

2.1.5. Transmission tomography
The most basic physical model for tomography assumes that wave energy traveling
between a source and receiver can be considered to be propagating along infinitesimally
narrow ray paths. In seismic tomography, if the slowness at a point x is s(x), and the
ray path is known, then the travel time for seismic energy transiting along that ray path
is given by the line integral along `:
Z
t = s(x(l))dl.
(2.4)
`

Learning of s(x) from t is required. Note that (2.4) is a high-dimensional generalization of (2.3). In reality, seismic ray paths will be bent due to refraction and/or reflection,
resulting in nonlinear inverse problem.
The above examples demonstrate the ubiquity of linear inverse problems. As a
result, in the next section we take up the case of linear inverse problems and illustrate
the Bayesian approach in details, also investigating connections with the deterministic
approach employed by the general scientific community.
3. Linear inverse problem
The motivating examples and discussions in this section are based on Bui-Thanh
[11].
Let us consider the following one-dimensional integral equation on a finite interval
as in equation (3.1):
Z
G(x, θ) =

K(x, t) θ(t) dt,

(3.1)

where K(x, ·) is some appropriate, known, real-valued function given x Now, let the
dataset be y n = (y1 , y2 , . . . , yn )T . Then for a known system response K(xi , t) for the
dataset, the equation can be written as follows:
Z
yi = G(xi , θ) + i ; i ∈ {1, 2, . . . , n}
(3.2)
R1
As a particular
example, let G(x, θ) = 0 K(x, t) θ(t) dt, where K(x, t) =

√ 1 2 exp −(x − t)2 /2ψ 2 is the Gaussian kernel and θ : [0, 1] 7→ R is to be
2πψ

learned given the data y n and xn = (x1 , . . . , xn )T . We first illustrate the Bayesian
approach and draw connections with the traditional approach of Tikhonov’s regularization when the integral in G is discretized. In this regard, let xi = (i − 1)/n, for
i = 1, . . . , n. Letting θ = (θ(x1 ), . . . , θ(xn ))T and K be the n × n matrix with the
(i, j)-th element K(xi , xj )/n, and n = (1 , . . . , n )T , the discretized version of (3.2)
can be represented as
y n = Kθ + n .
(3.3)

We assume that n ∼ Nn 0n , σ 2 I n , that is, an n-variate normal with mean 0n , an
n-dimensional vector with all components zero, and covariance σ 2 I n , where I n is the
n-th order identity matrix.
6

3.1. Smooth prior on θ
To reflect the belief that the function θ is smooth, one may presume that
θ(xi−1 ) + θ(xi+1 )
+ ˜i ,
(3.4)
2

iid
where, for i = 1, . . . , n, ˜i ∼ N 0, σ̃ 2 . Thus, a priori, θ(xi ) is assumed to be
an average of its nearest neighbors to quantify smoothness, with an additive random
perturbation term. Letting


−1 2 −1 0 · · · · · ·
 0 −1 2 −1 0 · · · 


..
..
..
..
..
.. 
1

.
.
.
.
. 
(3.5)
L=  .
,
2 .

.
.
.
.
.
..
..
..
..
.. 
 ..
0
0 · · · −1 2 −1
θ(xi ) =

and ˜ = (˜
1 , . . . , ˜n )T , it follows from (3.4) that
Lθ = ˜,

(3.6)

Now, noting that the Laplacian of a twice-differentiable real-valued function f with
Pk
2
independent arguments z1 , . . . , zk is given by ∆f = i=1 ∂∂zf2 , we have
i

∆θ(xj ) ≈ n2 (Lθ)j ,

(3.7)

where (Lθ)j is the j-th element of Lθ.
However, the rank of L is n−1, and boundary conditions on the Laplacian operator
is necessary to ensure positive definiteness of the operator. In our case, we assume
that θ ≡ 0 outside [0, 1], so that we now assume θ(0) = θ(x2 1 ) + ˜0 and θ(xn ) =

θ(xn−1 )
+ ˜n , where ˜0 and ˜n are iid N 0, σ̃ 2 . With this modification, the prior on
2
θ is given by


1
2
π(θ) ∝ exp − 2 kL̃θk ,
(3.8)
2σ̃
where k · k is the Euclidean norm and

2 −1
 −1 2

 0 −1

..
..
1
.
.
L̃ = 

2 .
.
..
 ..

 0
0
0
0

0
−1
2
..
.
..
.
···
···

7

0
0
−1
..
.
..
.
−1
0

···
···
0
..
.
..
.
2
−1

···
···
···
..
.
..
.
−1
2







.






(3.9)

Rather than assuming zero boundary conditions,

 moregenerally
 one may assume
σ̃ 2
σ̃ 2
that θ(0) and θ(xn ) are distributed as N 0, δ2 and N 0, δ2 , respectively. The
n
0
resulting modified matrix is then given by


2δ0 0
0
0 ··· ···
 −1 2 −1 0 · · · · · · 


 0 −1 2 −1 0 · · · 


..
..
..
..
.. 
1  ..
.
.
.
.
.
.
.
(3.10)
L̂ = 

2

 .
.
.
.
.
.
..
..
..
..
.. 
 ..


 0
0 · · · −1 2 −1 
0
0 ··· 0
0 2δn
To choose δ0 and δn , one may assume that
V ar [θ(0)] =

 T −1


σ̃ 2
σ̃ 2
2 T
[n/2] ,
=
V
ar
[θ(x
)]
=
=
V
ar
θ(x
)
=
σ̃

L̂ L̂
n
[n/2]
[n/2]
δ02
δn2

where [n/2] is the largest integer not exceeding n/2, and [n/2] is the [n/2]-th canonical
basis vector in Rn+1 . It follows that
δ02 = δn2 =

1

.
−1
[n/2]
T[n/2] L̂ L̂


T

Since this requires solving a non-linear equation (since L̂ contains δ0 and δn ), for
avoiding computational complexity one may simply employ the approximation
δ02 = δn2 =

1

,
−1
T[n/2] L̃ L̃
[n/2]


T

where L̃ is given by (3.9).
3.2. Non-smooth prior on θ
To begin with, let us assume that θ has several points of discontinuities on the grid
of points {x0 , . . . , xn }. To reflect this information in the prior, one may assume that
θ(0) = 0 and for i = 1, . . . , n, θ(xi ) = θ(xi−1 ) + ˜i , where, as before, ˜i are iid
N 0, σ̃ 2 . Then, with


1
0
0
0 ··· ···
 −1 1
0
0 ··· ··· 


 0 −1 1
0
0 ··· 


..
..
..
..
.. 
1  ..
,
.
.
.
.
.
.
L∗ = 
(3.11)

2
 .

.
.
.
.
.
..
..
..
..
.. 
 ..


 0
0 · · · −1 1
0 
0
0 ··· 0
−
1
8

the prior is given by



1
∗
2
π(θ) ∝ exp − 2 kL θk .
2σ̃

(3.12)

One may also flexibly account for any particular big jump. For instance, if for some
` ∈ {0, . . . , n}, the jump θ(x` ) − θ(x`−1 ) is particularly large compared
 to the
 other
2

jumps, then it can be assumed that θ(x` ) = θ(x`−1 ) + ∗` , with ∗` ∼ N 0, σ̃ξ2 , where

ξ < 1. Letting D ` be the diagonal matrix with ξ 2 being the `-th diagonal element and
1 being the other diagonal elements, the prior is then given by


1
π(θ) ∝ exp − 2 kD ` L∗ θk2 .
(3.13)
2σ̃
A more general prior can be envisaged where the number and location of the
jump discontinuities are unknown. Then we may consider a diagonal matrix D =
diag{ξ1 , . . . , ξn }, so that conditionally on the hyperparameters ξ1 , . . . , ξn , the prior
on θ is given by


1
∗
2
(3.14)
π(θ|ξ1 , . . . , ξn ) ∝ exp − 2 kDL θk .
2σ̃
Prior on ξ1 , . . . , ξn may be considered to complete the specification. These may also
be estimated by maximizing the marginal likelihood obtained by integrating out θ,
which is known as the ML-II method; see Berger [6]. Calvetti and Somersalo [12] also
advocate likelihood based methods.
3.3. Posterior distribution
For convenience, let us generically denote the matrices L, L̃, L̂, L∗ , D ` L∗ , DL∗ ,
1
by Γ− 2 . Then it can be easily verified that the posterior of θ admits the following
generic form:
 

1
1
− 21
2
2
π (θ|y n , xn ) ∝ exp −
ky
−
Kθk
+
kΓ
θk
.
(3.15)
2σ 2 n
2σ̃ 2
Note that the exponent of the posterior is of the form of the Tikhonov functional, which
we denote by T (θ). The maximizer of the posterior, commonly known as the maximum
a posteriori (MAP) estimator, is given by
θ̂ M AP = arg max π (θ|y n , xn ) = arg min T (θ).
θ

(3.16)

θ

In other words, the deterministic solution to the inverse problem obtained by Tikhonov’s
regularization is nothing but the Bayesian MAP estimator in our context.
Writing H = σ12 K T K + σ̃12 Γ−1 , which is the Hessian of the Tikhonov functional
1
(regularized misfit), and writing k·kH = kH 2 ·k, it is clear that (3.15) can be simplified
to the Gaussian form, given by
(
)
2
1 −1 −1
π (θ|y n , xn ) ∝ exp − θ − 2 H K y n
.
(3.17)
σ
H
9

It follows from (3.17) that the inverse of the Hessian of the regularized misfit is the
posterior covariance itself. From the above posterior it also trivially follows that
θ̂ M AP

1
1
= 2 H −1 K −1 y n = 2
σ
σ



1 T
1
K K + 2 Γ1
σ2
σ̃

−1

K T Y n,

(3.18)

which coincides with the Tikhonov solution for linear inverse problems. The connection between the traditional deterministic Tikhonov regularization approach with
Bayesian analysis continues to hold even if the likelihood is non-Gaussian.
3.4. Exploration of the smoothness conditions
For deeper investigation of the smoothness conditions, let us write


1
1
1
2
2
2
2
θ̂ M AP = arg min T (θ) = σ
ky − ỹ n k + %kΓ̃ θk ,
2 n
2
θ
1

(3.19)

1

where ỹ n = Kθ, % = σ 2 /σ̃ 2 and Γ̃ 2 = Γ− 2 . Now, from (3.7) it follows that for the
smooth priors with the zero boundary conditions, our Tikhonov functional discretizes
T∞ (θ) =

1
1
ky − ỹ n k2 + %k∆θk2L2 (0,1) ,
2 n
2

(3.20)

R1
where k · k2L2 (0,1) = 0 (·)2 dt.
On the other hand, for the non-smooth prior (3.12), rather than discretizing ∆θ,
∇θ, that is, the gradient of θ, is discretized. In other words, for non-smooth priors, our
Tikhonov functional discretizes
T∞ (θ) =

1
1
ky − ỹ n k2 + %k∇θk2L2 (0,1) .
2 n
2

(3.21)

Hence, realizations of prior (3.12) is less smooth compared to those of our smooth
priors. However, the realizations (3.12) must be continuous. The priors given by
(3.13) and (3.14) also support continuous functions as long as the hyperparameters
are bounded away from zero. These facts, although clear, can be rigorously justified
by functional analysis arguments, in particular, using the Sobolev imbedding theorem
(see, for example, Arbogast and Bona [1]).
4. Links between Bayesian inverse problems based on Gaussian process prior and
deterministic regularizations
In this section, based on Rasmussen and Williams [31], we illustrate the connections between deterministic regularizations such as those obtained from differential
operators as above, and Bayesian inverse problems based on the very popular Gaussian
process prior on the unknown function. A key tool for investigating such relationship
is the reproducing kernel Hilbert space (RKHS).

10

4.1. RKHS
We adopt the following definition of RKHS provided in Rasmussen and Williams
[31]:
Definition 4.1 (RKHS). Let H be a Hilbert space of real functions θ defined on an
index set X . Then H is called an RKHS endowed with an inner product h·, ·iH (and
norm kθkH = hθ, θiH ) if there exists a function K : X × X 7→ R with the following
properties:
(a) for every x, K(·, x) ∈ H, and
(b) K has the reproducing property hθ(·), K(·, x)iH = θ(x).
Observe that since K(·, x), K(·, x0 ) ∈ H, it follows that hK(·, x), K(·, x0 )iH =
K(x, x0 ). The Moore-Aronszajn theorem asserts that the RKHS uniquely determines
K, and vice versa. Formally,
Theorem 1 (Aronszajn [2]). . Let X be an index set. Then for every positive definite
function K(·, ·) on X × X there exists a unique RKHS, and vice versa.
R
Here, by positive definite function K(·, ·) on X×X, we mean K(x, x0 )g(x)g(x0 )dν(x)dν(x0 ) >
0 for all non-zero functions g ∈ L2 (X, ν), where L2 (X, ν) denotes the space of functions square-integrable on X with respect to the measure ν.
Indeed, the subspace H0 of H spanned by the functions {K(·, xi ); i = 1, 2, . . .}
is dense in H in the sense that every function in H is a pointwise limit of a Cauchy
sequence from H0 .
To proceed, we require the concepts of eigenvalues and eigenfunctions associated
with kernels. In the following section we provide a briefing on these.
4.2. Eigenvalues and eigenfunctions of kernels
We borrow the statements of the following definition of eigenvalue and eigenfunction, and the subsequent statement of Mercer’s theorem from Rasmussen and Williams
[31].
Definition 4.2. A function ψ(·) that obeys the integral equation
Z
C(x, x0 )ψ(x)dν(x) = λψ(x0 ),

(4.1)

X

is called an eigenfunction of the kernel C with eigenvalue λ with respect to the measure
ν.
We assume that the ordering is chosen such that λ1 ≥ λ2 ≥ · · · . The eigenfunctions
are orthogonal with respect to ν and can be chosen to be normalized so that
R
ψ
(x)ψ
i
j (x)dν(x) = δij , where δij = 1 if i = j and 0 otherwise.
X
The following well-known theorem (see, for example, König [21]) expresses the
positive definite kernel C in terms of its eigenvalues and eigenfunctions.

11

Theorem 2 (Mercer’s theorem). Let (X, ν) be a finite measure space and C ∈ L∞ X2 , ν 2
be a positive definite kernel. By L∞ X2 , ν 2 we mean the set of all measurable
functions C : X2 7→ R which are essentially bounded, that is, bounded up to a set
of ν 2 -measure zero. For any function C in this set, its essential supremum, given
by inf {C ≥ 0 : |C(x1 , x2 )| < C, for almost all (x1 , x2 ) ∈ X × X} serves as the norm
kCk.
Let ψj ∈ L2 (X, ν) be the normalized eigenfunctions of C associated with the eigenvalues λj (C) > 0. Then
∞

(a) the eigenvalues {λj (C)}j=1 are absolutely summable.
P∞
(b) C(x, x0 ) = j=1 λj (C)ψj (x)ψ¯j (x0 ) holds ν 2 -almost everywhere, where the series converges absolutely and uniformly ν 2 -almost everywhere. In the above, ψ¯j
denotes the complex conjugate of ψj .
It is important to note the difference between the eigenvalue λj (C) associated with
the kernel C and λj (Σn ) where Σn denotes the n × n Gram matrix with (i, j)-th
element C(xi , xj ). Observe that (see Rasmussen and Williams [31]):
0

Z

λj (C)ψj (x ) =

n

C(x, x0 )ψj (x)dν(x) ≈

X

1X
C(xi , x0 )ψj (xi , x0 ),
n i=1

(4.2)

where, for i = 1, . . . , n, xi ∼ ν, assuming that ν is a probability measure. Now
substituting x0 = xi ; i = 1, . . . , n in (4.2) yields the following approximate eigen
system for the matrix Σn :
Σn uj ≈ nλj (C)uj ,
(4.3)
where the i-th component of uj is given by
uij =

ψj (xi )
√ .
n

(4.4)

Since ψj are normalized to have unit norm, it holds that
n

uTj uj =

1X 2
ψ (xi ) ≈
n i=1 j

Z

ψ 2 (x)dν(x) = 1.

(4.5)

X

From (4.5) it follows that
λj (Σn ) ≈ nλj (C).
−1

(4.6)

Indeed, Theorem 3.4 of Baker [5] shows that n λj (Σn ) → λj (C), as n → ∞.
For our purposes the main usefulness of the RKHS framework is that kθk2H can
be perceived as a generalization of θ T K −1 θ, where θ = (θ(x1 ), . . . , θ(xn ))T and
K = (K(xi , xj ))i,j=1,...,n , is the n × n matrix with (i, j)-th element K(xi , xj ).

12



4.3. Inner product
Consider a real positive semidefinite kernel K(x, x0 ) with an eigenfunction exPN
pansion K(x, x0 ) = i=1 λi φi (x)φi (x0 ) relative to a measure µ. Mercer’s theorem
ensures
that the eigenfunctions are orthonormal with respect to µ, that is, we have
R
φi (x)φj (x)dµ(x) = δij . Consider a Hilbert space of linear combinations of the
PN
PN θi2
eigenfunctions, that is, θ(x) =
i=1 θi φi (x) with
i=1 λi < ∞. Then the inner
PN
PN
product hθ1 , θ2 iH between θ1 = i=1 θ1i φi (x), and θ2 = i=1 θ2i φi (x) is of the
form
N
X
θ1i θ2i
hθ1 , θ2 iH =
.
(4.7)
λi
i=1
PN θ2
This induces the norm k · kH , where kθk2H = i=1 λii . A smoothness condition on
the space is immediately imposed by requiring the norm to be finite – the eigenvalues
must decay sufficiently fast.
The Hilbert space defined above is a unique RKHS with respect to K, in that it
satisfies the following reproducing property:
hθ, K(·, x)i =

N
X
θi λi φi (x)
i=1

λi

= θ(x).

(4.8)

= K(x, x0 ).

(4.9)

Further, the kernel satisfies the following:
hK(x, ·), K(x0 , ·)i =

N
X
λ2 φi (x)
i

i=1

λi

PN
Now, with reference to (4.6), observe that the square norm kθk2H = i=1 θi2 /λi
and the quadratic form θ T Kθ have the same form if the latter is expressed in terms
of the eigenvectors of K, albeit the latter has n terms, while the square norm has N
terms.
4.4. Regularization
The ill-posed-ness of inverse problems can be understood from the fact that for
any given data set y n , all functions that pass through the data set minimize any given
measure of discrepancy D(y n , θ) between the data y n and θ. To combat this, one
considers minimization of the following regularized functional:
R(θ) = D(y n , θ) +

τ
kθk2H ,
2

(4.10)

where the second term, which is the regularizer, controls smoothness of the function
and τ is the appropriate Lagrange multiplier.
The well-known representer theorem (see, for example, Kimeldorf and Wahba [20],
O’Sullivan et al. [28], Wahba [38], Schölkopf and
Pn Smola [32]) guarantees that each
minimizer θ ∈ H can be represented as θ(x) = i=1 ci K (x, xi ), where K is the corresponding reproducing kernel. If D (y n , θ) is convex, then there is a unique minimizer
θ̂.
13

4.5. Gaussian process modeling of the unknown function θ
For simplicity, let us consider the model
yi = θ(xi ) + i ,

(4.11)

iid

for i = 1, . . . , n, where i ∼ N (0, σ 2 ), where we assume σ to be known for simplicity
of illustration. Let θ(x) be modeled by a Gaussian process with mean function µ(x)
and covariance kernel K associated with the RKHS. In other words, for any x ∈ X,
E [θ(x)] = µ(x) and for any x1 , x2 ∈ X, Cov (θ(x1 ), θ(x2 )) = K(x1 , x2 ).
Assuming for convenience that µ(x) = 0 for all x ∈ X, it follows that the posterior
distribution of θ(x∗ ) for any x∗ ∈ X is given by

π(θ(x∗ )|y n , xn ) ≡ N µ̂(x∗ ), σ̂ 2 (x∗ ) ,
(4.12)
where, for any x∗ ∈ X,
µ̂(x∗ ) = sT (x∗ ) K + σ 2 In

−1

yn ;

σ̂ 2 (x∗ ) = K(x∗ , x∗ ) − sT (x∗ ) bK + σ 2 In

(4.13)
−1

s(x∗ ),

(4.14)

T

with s(x∗ ) = (K(x∗ , x1 ), . . . , K(x∗ , xn )) .
Observe that the posterior mean admits the following representation:
µ̂(x∗ ) =

n
X

c̃i K(x∗ , xi ),

(4.15)

i=1

−1
where c̃i is the i-th element of K + σ 2 In
yn .
In other words, the posterior mean of the Gaussian process based model is consistent with the representer theorem.
5. Regularization using differential operators and connection with Gaussian process
For x = (x1 , . . . , xd )T ∈ Rd , let
m

2

Z

kL θk =

X

∂ m θ(x)

j1 +···+jd =m

∂xj11 · · · ∂xjdd

and
kPθk2 =

M
X

bm kLm θk2 ,

!2
,

(5.1)

(5.2)

m=0

for some M > 0, where the co-efficients bm ≥ 0. In particular, we assume for our
purpose that b0 > 0. It is clear that kPθk2 is translation and rotation invariant. This
norm penalizes θ in terms of its derivatives up to order M .

14

5.1. Relation to RKHS
It can be shown, using the fact that the complex exponentials exp(2πisT x) are
eigen functions of the differential operator, that
2

kPθk =

Z X
M

2
m
bm 4π 2 sT s
θ̃(s) ds,

(5.3)

m=0

where θ̃(s) is the Fourier transform of θ(s). Comparison of (5.3) with (4.7) yields the
hP
m i−1
M
2 T
power spectrum of the form
b
4π
s
s
which yields the following
m=0 m
kernel by Fourier inversion:
Z
exp(2πisT (x − x0 ))
0
0
ds.
(5.4)
K(x, x ) = K(x − x ) =
PM
2 T m
m=0 bm (4π s s)
Calculus of variations can also be used to minimize R(θ) with respect to θ, which
yields (using the Euler-Lagrange equation)
θ(x) =

n
X

bi G(x − xi ),

(5.5)

i=1

with

m
X

(−1)m bm ∇m G = δx−x0 ,

(5.6)

i=1

where G is known as the Green’s function. Using Fourier transform on (5.6) it can be
shown that the Green’s function
Pm is nothing but the kernel K given by (5.4). Moreover,
it follows from (5.6) that i=1 (−1)m bm ∇m and K are inverses of each other.
Examples of kernels derived from differential operators are as follows. For d = 1,
setting b0 = b2 , b1 = 1 and bm = 0 for m ≥ 2, one obtains K(x, x0 ) = K(x − x0 ) =
1
0
2b exp (−b|x − x |), which is the covariance of the Ornstein-Uhlenbeck process. For
m
0
0
general d dimension,
bm = b2m
 1 setting
 /(m!2 ), yields K(x, x ) = K(x − x ) =
1
0 T
0
exp − 2b2 (x − x ) (x − x ) .
(2πb2 )d/2
Considering a grid xn , note that
!
M
M
X
X
T
T
kPθk2 ≈
bm (Dm θ) (Dm θ) = θ T
Dm
Dm θ,
(5.7)
m=0

m=0

where Dm is a suitable finite-difference approximation of the differential operator.
Note that such finite-difference approximation has been explored in Section 3, which
we now investigate in a rigorous setting. Also, since (5.7) is quadratic in θ, assuming a prior for θ, the logarithm of which has this form, and further assuming that
log [D(y n , θ)] is a log-likelihood quadratic in θ, a Gaussian posterior results.

15

5.2. Spline models and connection with Gaussian process
Let us consider the penalty function to be kLm θk2 . Then polynomials up to degree
m − 1 are not penalized and so, are in the null space of the regularization operator. In
this case, it can be shown that a minimizer of R(θ) is of the form
θ(x) =

k
X

dj ψj (x) +

j=1

n
X

ci G(x, xi ),

(5.8)

i=1

where {ψ1 , . . . , ψk } are polynomials that span the null space and the Green’s function
G is given by (see Duchon [15], Meinguet [23])

cm,d |x − x0 |2m−d log |x − x0 | if 2m > d and d even
G(x, x0 ) = G(x − x0 ) =
,
cm,d |x − x0 |2m−d
otherwise.
(5.9)
where cm,D are constants (see Wahba [38] for the explicit form).
We now specialize the above arguments to the spline set-up. As before, let us con
iid
sider the model yi = θ(xi ) + i , where, for i = 1, . . . , n, i ∼ N 0, σ 2 . For simplicity, we consider the one-dimensional set-up, and consider the cubic spline smoothing
problem that minimizes
Z
n
X
2
R(θ) =
(yi − θ(xi )) + τ

1

2

[θ00 (x)] dx,

(5.10)

0

i=1

where 0 < x1 < · · · < xn < 1. The solution to this minimization problem is given by
θ(x) =

1
X

dj xj +

n
X

ci (x − xi )3+ ,

(5.11)

i=1

j=0

where, for any x, (x)+ = x if x > 0 and zero otherwise.
Following Wahba [37], let us consider
f (x) =

1
X

βj xj + θ(x),

(5.12)

j=0



where β = (β0 , β1 )T ∼ N 0, σβ2 I2 , and θ is a zero mean Gaussian process with
covariance


Z 1
v3
|x − x0 |v 2
+
,
(5.13)
σθ2 K(x, x0 ) =
(x − u)+ (x0 − u)+ du = σθ2
2
3
0
where v = min{x, x0 }.
Taking σβ2 → ∞ makes the prior of β vague, so that penalty on the polynomial
terms in the null space is effectively washed out. It follows that


−1
E [θ(x∗ )|y n , xn ] = h(x∗ )T β̂ + s(x∗ )T K̂
y n − H T β̂ ,
(5.14)
16

where, for any x, h(x) = (1, x)T , H = (h(x1 ), . . . , h(xn )), K̂ is the covariance

−1
−1
−1
matrix corresponding to σθ2 K(xi , xj ) + σ 2 δij , and β̂ = H K̂ H
H K̂ y n .
Since the elements of s(x∗ ) are piecewise cubic polynomials, it is easy to see that
the posterior mean (5.14) is also a piecewise cubic polynomial. It is also clear that
(5.14) is a first order polynomial on [0, x1 ] and [xn , 1].
5.2.1. Connection with the `-fold integrated Wiener process
Shepp [33] considered the `-fold integrated Wiener process, for ` = 0, 1, 2 . . ., as
follows:
Z 1
(x − u)`+
Z(u)du,
(5.15)
W` (x) =
`!
0
where Z is a Gaussian white noise process with covariance δ(u − u0 ). As a special
case, note that W0 is the standard Wiener process. In our case, note that
K(x, x0 ) = Cov (W1 (x), W1 (x0 )) .

(5.16)

2
R
The above ideas can be easily extended to the case of the regularizer f (m) (x) dx,
for m ≥ 1 by replacing (x − u)+ with (x − u)m−1
/(m − 1)! and letting h(x) =
+

m−1 T
1, x, . . . , x
.
6. The Bayesian approach to inverse problems in Hilbert spaces
We assume the following model
y = G(θ) + ,

(6.1)

where y, θ and  are in Banach or Hilbert spaces.
6.1. Bayes theorem for general inverse problems
We will consider the model stated by equation (6.1). Let Y and Θ denote the
sample spaces for y and θ, respectively. Let us first assume that both are separable
Banach spaces. Assume µ0 to be the prior measure for θ. Assuming well-defined joint
distribution for (y, θ), let us denote the posterior of θ given y as µy . Let  ∼ Q0 where
Q0 such that  and θ are independent. Let Q0 be the distribution of . Let us denote the
conditional distribution of y given θ by Qθ , obtained from a translation of Q0 by G(θ).
Assume that Qθ  Q0 . Thus, for some potential Φ : Θ × Y 7→ R,
dQθ
= exp (−Φ(θ, y)) .
dQ0

(6.2)

Thus, for fixed θ, Φ(θ, ·) : Y 7→ R is measurable and EQ0 [exp (−Φ(θ, y))] = 1. Note
that −Φ(·, y) is nothing but the log-likelihood.
Let ν0 denote the product measure
ν0 (dθ, dy) = µ0 (dθ)Q0 (dy),
17

(6.3)

and let us assume that Φ is ν0 -measurable. Then (θ, y) ∈ Θ×Y is distributed according
to the measure ν(dθ, dy) = µ0 (dθ)Qθ (dy). It then also follows that ν  ν0 , with
dνθ
(θ, y) = exp (−Φ(θ, y)) .
dν0

(6.4)

Then we have the following statement of Bayes’ theorem for general inverse problems:
Theorem 3 (Bayes theorem for general inverse problems). Assume that Φ : Θ × Y 7→
R is ν0 -measurable and
Z
C=
exp (−Φ(θ, y)) µ0 (dy) > 0,
(6.5)
Θ

for Q0 -almost surely all y. Then the posterior of θ given y, which we denote by µy ,
exists under ν. Also, µy  µ0 and for all y ν0 -almost surely,
dµyθ
1
(θ) = exp (−Φ(θ, y)) .
dµ0
C

(6.6)

Now assume that Θ and Y are Hilbert spaces. Suppose  ∼ N(0, Γ). Then the
following theorem holds:
Theorem 4 (Vollmer [36]).


dµy
1
∝ exp − kG(θ)k2Γ + hy, G(θ)iΓ ,
dµ0
2

(6.7)

where h·, ·iΓ = hΓ−1 ·, ·i, and k · kΓ is the norm induced by h·, ·iΓ .

iid
For the model yi = θ(xi ) + i for i = 1, . . . , n, with i ∼ N 0, σ 2 , the posterior
is of the form
!
n
2
X
dµy
(yi − θ(xi ))
.
(6.8)
∝ exp −
dµ0
2σ 2
i=1
6.2. Connection with regularization methods
It is not immediately clear if the Bayesian approach in the Hilbert space setting
has connection with the deterministic regularization methods, but Vollmer [36] prove
consistency of the posterior assuming certain stability results which are used to prove
convergence of regularization methods; see Engl et al. [16].
We next turn to inverse regression.
7. Inverse regression
We first provide some examples of inverse regression, mostly based on Avenhaus
et al. [4].

18

7.1. Examples of inverse regression
7.1.1. Example 1: Measurement of nuclear materials
Measurement of the amount of nuclear materials such as plutonium by direct chemical means is an extremely difficult exercise. This motivates model-based methods.
For instance, there are physical laws relating heat production or the number of neutrons emitted (the dependent response variable y) to the amount of material present,
the latter being the independent variable x. But any measurement instrument based on
the physical laws first needs to be calibrated. In other words, the unknown parameters
of the model needs to be learned, using known inputs and outputs. However, the independent variables are usually subject to measurement errors, motivating a statistical
model. Thus, conditionally on x and parameter(s) θ, y ∼ P (·|x, θ), where P (·|x, θ)
denotes some appropriate probability model. Given y n and xn , and some specific ỹ,
the corresponding x̃ needs to be predicted.
7.1.2. Example 2: Estimation of family incomes
Suppose that it is of interest to estimate the family incomes in a certain city through
public opinion poll. Most of the population, however, will be unwilling to provide
reliable answers to the questionnaires. One way to extract relatively reliable figures is
to consider some dependent variable, say, housing expenses (y), which is supposed to
strongly depend on family income (x); see Muth [26], and such that the population is
less reluctant to divulge the correct figures related to y. From past survey data on xn
and y n , and using current data from families who may provide reliable answers related
to both x and y, a statistical model may be built, using which the unknown family
incomes may be predicted, given their household incomes.
7.1.3. Example 3: Missing variables
In regression problems where some of the covariate values xi are missing, they may
be estimated from the remaining data and the model. In this context, Press and Scott
[29] considered a simple linear regression problem in a Bayesian framework. Under
special assumptions about the error and prior distributions, they showed that an optimal
procedure for estimating the linear parameters is to first estimate the missing xi from
an inverse regression based only on the complete data pairs.
7.1.4. Example 4: Bioassay
It is usual to investigate the effects of substances (y) given in several dosages on
organisms (x) using bioassay methods. In this context it may be of interest to determine the dosage necessary to obtain some interesting effect, making inverse regression
relevant (see, for example, Rasch et al. [30]).
7.1.5. Example 5: Learning the Milky Way
The modelling of the Milky Way galaxy is an integral step in the study of galactic
dynamics; this is because knowledge of model parameters that define the Milky Way
directly influences our understanding of the evolution of our galaxy. Since the nature
of the Galaxy’s phase space, in the neighbourhood of the Sun, is affected by distinct
Milky Way features, measurements of phase space coordinates of individual stars that

19

live in this neighbourhood of the Sun, will bear information about the influence of such
features. Then, inversion of such measurements can help us learn the parameters that
describe such Milky Way features. In this regard, learning about the location of the Sun
with respect to the center of the galaxy, given the two-component velocities of the stars
in the vicinity of the Sun, is an important problem. For k such stars, Chakrabarty et al.
[13] model the k × 2-dimensional velocity matrix V as a function of the galactocentric
location (S) of the Sun, denoted by V = ξ(S). For a given observed value V ∗ of V ,
it is then of interest to obtain the corresponding S ∗ . Since ξ is unknown, Chakrabarty
et al. [13] model ξ as a matrix-variate Gaussian process, and consider the Bayesian
approach to learning about S ∗ , given data {(S i , V i ) : i = 1, . . . , n} simulated from
established astrophysical models, and the observed velocity matrix V ∗ .
We now provide a brief overview of of the methods of inverse linear regression,
which is the most popular among inverse regression problems. Our discussion is generally based on Hoadley [19] and Avenhaus et al. [4].
7.2. Inverse linear regression
Let us consider the following simple linear regression model: for i = 1, . . . , n,
yi = α + βxi + σi ,

(7.1)

iid

where i ∼ N (0, 1).
For simplicity, let us consider a single unknown x̃, associated with a further set of
m responses {ỹ1 , . . . , ỹm }, related by
ỹi = α + β x̃ + τ ˜i ,

(7.2)

iid

for i = 1, . . . , m, where ˜i ∼ N (0, 1) and are independent of the i ’s associated with
(7.1).
The interest in the above problem is inference regarding the unknown x. Based on
(7.1), first least squares estimates of α and β are obtained as
Pn
(y − ȳ)(xi − x̄)
Pn i
;
(7.3)
β̂ = i=1
2
i=1 (xi − x̄)
α̂ = ȳ − β̂ x̄,
(7.4)
Pn
Pn
where ȳ = i=1 yi /n and ȳ = i=1 xi /n. Then, letting ỹ¯ = i=1 ỹi /n, a ‘classical’ estimator of x is given by
ỹ¯ − α̂
,
(7.5)
x̂C =
β̂
which is also the maximum likelihood estimator for the likelihood associated with (7.1)
and (7.2), assuming known σ and τ . However,
h
i
2
E (x̂C − x) |α, β, σ, τ, x = ∞,
(7.6)
Pn

which prompted Krutchkoff [22] to propose the following ‘inverse’ estimator:
¯
x̂I = γ̂ + δ̂ ỹ,
20

(7.7)

where
Pn
(y − ȳ)(xi − x̄)
Pn i
;
δ̂ = i=1
2
i=1 (yi − ȳ)

(7.8)

γ̂ = x̄ − δ̂ ȳ,

(7.9)

are the least squares estimators of the slope and intercept when the xi are regressed on
the yi . It can be shown that the mean square error of this inverse estimator is finite.
However, Williams [39] showed that if σ 2 = τ 2 and if the sign of β is known, then
the unique unbiased estimator of x has infinite variance. Williams advocated the use of
confidence limits instead of point estimators.
Hoadley [19] derive
Pn confidence limits setting σ = τ and assuming without loss
of generality that i=1 xi = 0. Under these assumptions, the maximum likelihood
estimators of σ 2 with xn and y n only, ỹ n = (ỹ1 , . . . , ỹn )T only, and with the entire
available data set are, respectively,
n

σ̂12 =

2
1 X
yi − α̂ − β̂xi ;
n − 2 i=1

σ̂22 =

1 X
¯ 2;
(ỹi − ỹ)
m − 1 i=1

(7.11)

σ̂ 2 =



1
(n − 2)σ12 + (m − 1)σ22 .
n−2+m−1

(7.12)

(7.10)

n

2

Now consider the F -statistic F = nσ̂β̂2 for testing the hypothesis β = 0. Note that
under the null hypothesis this statistic has the F distribution with 1 and n + m degrees
of freedom. For m = 1,
r
n
β̂ (x̂C − x)
σ 2 (n + 1 + x2 )
has a t distribution with n − 2 degrees of freedom. Letting Fα;1,ν denote the upper α
point of the F distribution with 1 and ν degrees of freedom, a confidence set S can be
derived as follows:

{x : xL ≤ x ≤ xU }
if
F > Fα;1,n−2 ;


n+1
{x
:
x
≤
x
}
∪
{x
≥
x
}
if
F
≤ F < Fα;1,n−2 ;
2
L
U
α;1,n−2
S=
n+1+x̂C

n+1

(−∞, ∞)
if
F < n+1+x̂2 Fα;1,n−2 ,
C
(7.13)
where xL and xU are given by


 1
Fα;1,n−2 (n + 1) (F − Fα;1,n−2 ) + F x̂2C 2
F x̂C
±
.
F − Fα;1,n−1
F − Fα;1,n−2
n+1
Hence, if F < n+1+x̂
2 Fα;1,n−2 , then the associated confidence interval is S =
C
(−∞, ∞), which is of course useless.
Hoadley [19] present a Bayesian analysis of this problem, presented below in the
form of the following two theorems.

21

Theorem 5 (Hoadley [19]). Assume that σ = τ , and let x be independent of (α, β, σ 2 )
a priori. With any prior π(x) on x and the prior
π(α, β, σ 2 ) ∝

1
σ2

on (α, β, σ 2 ), the posterior density of x given by
π(x|y n , xn , ỹ n ) ∝ π(x)L(x),
where
L(x) = h
1+

n
m

 m+n−3
n
2
+ x2
1+ m
,


i m+n−2
2
2
F
2
+ Rx̂C + m+n−3 + 1 (x − Rx̂C )

where
R=

F
.
F +m+n−3

For m = 1, Hoadley [19] present the following result characterizing the inverse
estimator x̂I :
Theorem 6 (Hoadley [19]). Consider the following informative prior on x:
x = tn−3

n+1
,
n−3

where tν denotes the t distribution with ν degrees of freedom. Then the posterior
distribution of x given y n , xn and ỹ n has the same distribution as
s
x̂2
n + 1 + RI
.
x̂I + tn−2
F +n−2
In particular, it follows from Theorem 6 that the posterior mean of x is x̂I when
m = 1. In other words, the inverse estimator x̂I is Bayes with respect to the squared
error loss and a particular informative prior distribution for x.
Since the goal of Hoadley [19] was to provide a theoretical justification of the
inverse estimator, he had to choose a somewhat unusual prior so that it leads to x̂I as
the posterior mean. In general it is not necessary to confine ourselves to any specific
prior for Bayesian analysis of inverse regression. It is also clear that the Bayesian
framework is appropriate for any inverse regression problem, not just linear inverse
regression; indeed, the palaeoclimate reconstruction problem (Haslett et al. [18]) and
the Milky Way problem (Chakrabarty et al. [13]) are examples of very highly non-linear
inverse regression problems.
7.3. Connection between inverse regression problems and traditional inverse problems
Note that the class of inverse regression problems includes the class of traditional
inverse problems. The Milky Way problem is an example where learning the unknown,
22

matrix-variate function ξ (inverse problem) was required, even though learning about
S, the galactocentric location of the sun (inverse regression problem) was the primary
goal. The Bayesian approach allowed learning both S and ξ simultaneously and coherently.
In the palaeoclimate models proposed in Haslett et al. [18], Bhattacharya [7] and
Mukhopadhyay and Bhattacharya [25], although species assemblages are modeled conditionally on climate variables, the functional relationship between species and climate
are not even approximately known. In all these works, it is of interest to learn about the
functional relationship as well as to predict the unobserved climate values, the latter
being the main aim. Again, the Bayesian approach facilitated appropriate learning of
both the unknown quantities.
7.4. Consistency of inverse regression problems
In the above linear inverse regression, notice that if τ > 0, then the variance of the
estimator of x can not tend to zero, even as the data size tends to infinity. This shows
that no estimator of x can be consistent. The same argument applies even to Bayesian
approaches; for any sensible prior on x that does not give point mass to the true value
of x, the posterior of x will not converge to the point mass at the true value of x as the
data size increases indefinitely. The arguments remain valid for any inverse regression
problem where the response variable y probabilistically depends upon the independent
variable x. Not only in inverse regression problems, even in forward regression problems where the interest is in prediction of y given x, any estimate of y or any posterior
predictive distribution y will be inconsistent.
To give an example of inconsistency in non-linear and non-normal inverse problem,
iid
consider the following set-up: yi ∼ Poisson (θxi ), for i = 1, . . . , n, where θ > 0 and
xi > 0 for each i. Let us consider the prior π(θ) ≡ 1 for all θ > 0. For some
i∗ ∈ {1, . . . , n} let us assume the leave-one-out cross-validation set-up in that we wish
to learn x = xi∗ assuming it is unknown, from the rest of the data. Putting the prior
π(x) ≡ 1 for x > 0, the posterior of x is given by (see Bhattacharya and Haslett [9],
Bhattacharya [8])
π(x|xn \xi , y n ) ∝

xyi
(x +

P

j6=i

xj )(

Pn

j=1

yj +1)

.

(7.14)

Figure 7.1 displays the posterior of x when i∗ = 10, for increasing sample size. Observe that the variance of the posterior does not decrease even with sample size as large
as 100, 000, clearly demonstrating inconsistency. Hence, special, innovative priors are
necessary for consistency in such cases.
8. Conclusion
In this review article, we have clarified the similarities and dissimilarities between
the traditional inverse problems and the inverse regression problems. In particular, we
have argued that only the latter class of problems qualify as authentic inverse problems
in they have significantly different goals compared to the corresponding forward problems. Moreover, they include the traditional inverse problems on learning unknown
23

0.30

Posterior of x

n=10

0.25

n=100

n=10000

0.15

n=100000

0.00

0.05

0.10

posterior density

0.20

n=1000

0

5

10

15

20

x

Figure 7.1: Demonstration of posterior inconsistency in inverse regression problems. The vertical line denotes the true value.

functions as a special case, as exemplified by our palaeoclimate and Milky Way examples. We advocate the Bayesian paradigm for both classes of problems, not only
because of its inherent flexibility, coherency and posterior uncertainty quantification,
but also because the prior acts as a natural penalty which is very important to regularize
the so-called ill-posed inverse problems. The well-known Tikhonov regularizer is just
a special case from this perspective.
It is important to remark that the literature on inverse function learning problems
and inverse regression problems is still very young and a lot of research is necessary to
develop the fields. Specifically, there is hardly any well-developed, consistent model
adequacy test or model comparison methodology in either of the two fields, although
Mohammad-Djafari [24] deal with some specific inverse problems in this context, and
Bhattacharya [8] propose a test for model adequacy in the case of inverse regression
problems. Moreover, as we have demonstrated, inverse regression problems are inconsistent in general. The general development in these respects will be provided in the
PhD thesis of the first author.

24

References
References
[1] Arbogast, T., and Bona, J. L. [2008], “Methods of Applied Mathematics,”. University of Texas at Austin.
[2] Aronszajn, N. [1950], “Theory of Reproducing Kernels,” Transactions of the
American Mathematical Society, 68, 337–404.
[3] Aster, R. C., Borchers, B., and Thurber, C. H. [2013], Parameter Estimation and
Inverse Problems, Oxford, UK: Academic Press.
[4] Avenhaus, R., Höpfinger, E., and Jewell, W. S. [1980], “Approaches
to Inverse Linear Regression,”.
Technical Report. Available at
https://publikationen.bibliothek.kit.edu/270015256/3812158.
[5] Baker, C. T. H. [1977], The Numerical Treatment of Integral Equations, Oxford:
Clarendon Press.
[6] Berger, J. O. [1985], Statistical Decision Theory and Bayesian Analysis, New
York: Springer-Verlag.
[7] Bhattacharya, S. [2006], “A Bayesian Semiparametric Model for Organism Based
Environmental Reconstruction,” Environmetrics, 17(7), 763–776.
[8] Bhattacharya, S. [2013], “A Fully Bayesian Approach to Assessment of Model
Adequacy in Inverse Problems,” Statistical Methodology, 12, 71–83.
[9] Bhattacharya, S., and Haslett, J. [2007], “Importance Resampling MCMC for
Cross-Validation in Inverse Problems,” Bayesian Analysis, 2, 385–408.
[10] Bühlmann, P., and van de Geer, S. [2011], Statistics for High-Dimensional Data,
New York: Springer.
[11] Bui-Thanh, T. [2012], “A Gentle Tutorial on Statistical Inversion
Using the Bayesian Paradigm,”.
ICES Report 12-18. Available at
http://users.ices.utexas.edu/ tanbui/PublishedPapers/BayesianTutorial.pdf.
[12] Calvetti, D., and Somersalo, E. [2007], Introduction to Bayesian Scientific Computing: Ten Lectures on Subjective Computing, New York: Springer.
[13] Chakrabarty, D., Biswas, M., and Bhattacharya, S. [2015], “Bayesian Nonparametric Estimation of Milky Way Parameters Using Matrix-Variate Data, in a New
Gaussian Process Based Method,” Electronic Journal of Statistics, 9, 1378–1403.
[14] Dashti, M., and Stuart, A. M. [2015], “The Bayesian Approach to Inverse Problems,”. eprint: arXiv:1302.6989.
[15] Duchon, J. [1977], Splines Minimizing Rotation-Invariant Semi-norms in
Sobolev Spaces,, in Constructive Theory of Functions of Several Variables, eds.
W. Schempp, and K. Zellner, Springer-Verlag, New York, pp. 85–100.
25

[16] Engl, H. W., Hanke, M., and Neubauer, A. [1996], Regularization of Inverse
Problems, Dordrecht: Kluwer Academic Publishers Group. Volume 375 of Mathematics and its Applications.
[17] Giraud, C. [2015], Introduction to High-Dimensional Statistics, New York: Chapman and Hall.
[18] Haslett, J., Whiley, M., Bhattacharya, S., Salter-Townshend, M., Wilson, S. P.,
Allen, J. R. M., Huntley, B., and Mitchell, F. J. G. [2006], “Bayesian Palaeoclimate Reconstruction (with discussion),” Journal of the Royal Statistical Society:
Series A (Statistics in Society), 169, 395–438.
[19] Hoadley, B. [1970], “A Bayesian Look at Inverse Linear Regression,” Journal of
the American Statistical Association, 65, 356–369.
[20] Kimeldorf, G., and Wahba, G. [1971], “Some Results on Tchebycheffian Spline
Functions,” Journal of Mathematical Analysis and Applications, 33, 82–95.
[21] König, H. [1986], Eigenvalue Distribution of Compact Operators, : Birkhäuser.
[22] Krutchkoff, R. G. [1967], “Classical and Inverse Regression Methods of Calibration,” Technometrics, 9, 425–435.
[23] Meinguet, J. [1979], “Multivariate Interpolation at Arbitrary Points Made Simple,” Journal of the Applied Mathematics and Physics, 30, 292–304.
[24] Mohammad-Djafari, A. [2000], “Model Selection for Inverse Problems: Best
Choice of Basis Function and Model Order Selection,”.
Available at
https://arxiv.org/abs/math-ph/0008026.
[25] Mukhopadhyay, S., and Bhattacharya, S. [2013], “Cross-Validation Based Assessment of a New Bayesian Palaeoclimate Model,” Environmetrics, 24, 550–
568.
[26] Muth, R. F. [1960], The Demand for Non-Farm Housing,, in The Demand for
Durable Goods, ed. A. C. Harberger. The University of Chicago.
[27] O’Sullivan, F. [1986], “A Statistical Perspective on Ill-Posed Inverse Problems,”
Statistical Science, 1, 502–512.
[28] O’Sullivan, F., Yandell, B. S., and Raynor, W. J. [1986], “Automatic Smoothing
of Regression Functions in Generalized Linear Models,” Journal of the American
Statistical Association, 81, 96–103.
[29] Press, S. J., and Scott, A. [1975], Missing Variables in Bayesian Regression,, in
Studies in Bayesian Econometrics and Statistics, eds. S. E. Fienberg, and A. Zellner, North-Holland, Amsterdam.
[30] Rasch, D., Enderlein, G., and Herrendörfer, G. [1973], “Biometrie,”. Deutscher
Landwirtschaftsverlag, Berlin.

26

[31] Rasmussen, C. E., and Williams, C. K. I. [2006], Gaussian Processes for Machine
Learning, Cambridge, Massachusetts: The MIT Press.
[32] Schölkopf, B., and Smola, A. J. [2002], Learning with Kernels, USA: MIT Press.
[33] Shepp, L. A. [1966], “Radon-Nikodym Derivatives of Gaussian Measures,” Annals of Mathematical Statistics, 37, 321–354.
[34] Tikhonov, A. [1963], “Solution of Incorrectly Formulated Problems and the
Reguarization Method,” Soviet Math. Dokl., 5, 1035–1038.
[35] Tikhonov, A., and Arsenin, V. [1977], Solution of Ill-Posed Problems, New York:
Wiley.
[36] Vollmer, S. [2013], “Posterior Consistency for Bayesian Inverse Problems
Through Stability and Regression Results,” Inverse Problems, 29. Article number
125011.
[37] Wahba, G. [1978], “Improper Priors, Spline Smoothing and the Problem of
Guarding Against Model Errors in Regression,” Journal of the Royal Statistical
Society B, 40, 364–372.
[38] Wahba, G. [1990], “Spline Functions for Observational Data,”. CBMS-NSF Regional Conference series, SIAM. Philadelphia.
[39] Williams, E. J. [1969], “A Note on Regression Methods in Calibration,” Technometrics, 11, 189–192.

27

