Diffusive optical tomography in the Bayesian framework

TL;DR

This paper investigates the relationship between inverse problems for the radiative transfer equation and the diffusion equation in optical tomography, proving convergence of their Bayesian posteriors in the strong scattering regime.

Contribution

It develops a rigorous framework to compare inverse problems for RTE and DE within the Bayesian setting, demonstrating posterior convergence in the asymptotic regime.

Findings

Proves convergence of RTE inverse problem to DE inverse problem in Bayesian framework.

Uses Hellinger metric and Kullback-Leibler divergence to quantify posterior convergence.

Establishes theoretical foundation for multiscale approximation in optical tomography.

Abstract

Many naturally-occuring models in the sciences are well-approximated by simplified models, using multiscale techniques. In such settings it is natural to ask about the relationship between inverse problems defined by the original problem and by the multiscale approximation. We develop an approach to this problem and exemplify it in the context of optical tomographic imaging. Optical tomographic imaging is a technique for infering the properties of biological tissue via measurements of the incoming and outgoing light intensity; it may be used as a medical imaging methodology. Mathematically, light propagation is modeled by the radiative transfer equation (RTE), and optical tomography amounts to reconstructing the scattering and the absorption coefficients in the RTE from boundary measurements. We study this problem in the Bayesian framework, focussing on the strong scattering regime. In this regime the forward RTE is close to the diffusion equation (DE). We study the RTE in the asymptotic regime where the forward problem approaches the DE, and prove convergence of the inverse RTE to the inverse DE in both nonlinear and linear settings. Convergence is proved by studying the distance between the two posterior distributions using the Hellinger metric, and using Kullback-Leibler divergence.

Equations

$$v\cdot \nabla f(x,v)={\int }_{{\mathbb{S}}^{d-1}}k(x,v,v^{\prime })f(x,v^{\prime })\text{d}{v}^{\prime }-\alpha (x,v)f(x,v).$$

$$-\nabla \cdot \left(a(x)\nabla \rho (x)\right)=0.$$

$$y=\mathcal{G}(\sigma )+\eta ,$$

$${\mu }^y(d\sigma )=\frac{1}{Z}{\mu }^{\sigma }(y){\mu }_0(d\sigma ),$$

$$Z={\int }_X{\mu }^{\sigma }(y)\text{d}{\mu }_0(d\sigma )$$

$${\mu }^y(\sigma )=\frac{1}{Z}\exp \left(-{\left(\sigma -m_0\right)}^{\top }{\mathcal{C}}_{\text{prior}}^{-1}\left(\sigma -m_0\right)-{\left(y-\mathcal{G}(\sigma )\right)}^{\top }{\mathcal{C}}_{\text{error}}^{-1}\left(y-\mathcal{G}(\sigma )\right)\right).$$

$${\mathrm{d}}_{\text{KL}}(\mu ,{\mu }^{\prime })={\int }_X\log \frac{d\mu }{d{\mu }^{\prime }}(\sigma )\text{d}\mu (d\sigma );$$

$$d_{\text{Hell}}(\mu ,{\mu }^{\prime }{)}^2=\frac{1}{2}{\int }_X{\left(\sqrt{\frac{d\mu }{d\lambda }}-\sqrt{\frac{d{\mu }^{\prime }}{d\lambda }}\right)}^2\text{d}\lambda .$$

$$y=\mathcal{F}(u)+\eta .$$

$$y\approx \mathcal{F}(u_0)+D\mathcal{F}(u_0)v+\eta .$$

$$z=Gv+\eta$$

$${\mathcal{C}}_{\text{post}}^{-1}={\mathcal{C}}_{\text{prior}}^{-1}+G^T{\mathcal{C}}_{\text{error}}^{-1}G,\text{and}{m}_{\text{post}}=m_{\text{prior}}+{\mathcal{C}}_{\text{post}}{G}^T{\mathcal{C}}_{\text{error}}^{-1}\left(z-G{m}_{\text{prior}}\right).$$

$$\alpha (x,v)=k(x,v,v^{\prime })={ϵ}^{-1}\sigma (x),$$

$${⟨1⟩}_v={\int }_{{\mathbb{S}}^{d-1}}1\text{d}v=1,$$

$$v\cdot \nabla f=\frac{1}{ϵ}\sigma \mathcal{L}f,$$

$$\mathcal{L}f={\int }_{{\mathbb{S}}^{d-1}}f(x,v^{\prime })\text{d}{v}^{\prime }-f={⟨f⟩}_v-f.$$

$${\Gamma }_{\pm }=\{(x,v):x\in \partial \Omega ,\pm v\cdot n_x>0\}$$

$${\Gamma }_{\pm }(y)=\{(x,v):x=y,\pm v\cdot n_y>0\},$$

$$\Gamma ={\Gamma }_{+}\cup {\Gamma }_{-},\text{and}\Gamma (y)={\Gamma }_{+}(y)\cup {\Gamma }_{-}(y).$$

$$f{|}_{{\Gamma }_{-}}=\varphi (x,v).$$

$$\{\begin{array}{l}v\cdot \nabla f=\frac{1}{ϵ}\sigma \mathcal{L}f,(x,v)\in \Omega \times {\mathbb{S}}^{d-1}\\ f{|}_{{\Gamma }_{-}}=\varphi (x,v)\end{array},$$

$$\begin{cases}-\nabla\cdot\left(\frac{1}{\sigma}\nabla\rho\right)=0,\quad x\in\Omega\subset\mathbb{R}^{d}\\ \rho\big{|}_{\partial\Omega}=\xi(x)\end{cases}\,,$$

$$\max \{\Vert \sigma {\Vert }_{L_{\infty }(\Omega )},\Vert {\sigma }^{-1}{\Vert }_{L_{\infty }(\Omega )},\Vert \nabla \left({\sigma }^{-1}\right){\Vert }_{L_{\infty }(\Omega )}\}<C_1;$$

$$\max \{\Vert \xi {\Vert }_{L_{\infty }(\partial \Omega )},\Vert \varphi {\Vert }_{L_{\infty }(\Gamma )}\}<C_1.$$

$$\mathcal{A}=\{\sigma \in C^3(\Omega ):\max \{\Vert \sigma {\Vert }_{L_{\infty }(\Omega )},\Vert {\sigma }^{-1}{\Vert }_{L_{\infty }(\Omega )},\Vert \nabla \left({\sigma }^{-1}\right){\Vert }_{L_{\infty }(\Omega )}\}<C_1\}.$$

$$\underset{\sigma \in \mathcal{A}}{\sup }\Vert \frac{1}{\sigma }{\partial }_n\rho {\Vert }_{\infty }<C.$$

$$\Vert f-\rho {\Vert }_{L_{\infty }(\Omega \times {\mathbb{S}}^{d-1})}<C_{\mathcal{A}}ϵ;$$

$$\Vert f-\rho +\frac{ϵ}{\sigma }v\cdot \nabla \rho {\Vert }_{L_{\infty }(\Omega \times {\mathbb{S}}^{d-1})}<C_{\mathcal{A}}{ϵ}^2.$$

$$f_{\text{in}}=f_0+\text{ϵ}{f}_1+{\text{ϵ}}^2{f}_2+\cdots .$$

$$v\cdot \nabla f_0+\text{ϵ}v\cdot \nabla f_1+{\text{ϵ}}^{2}v\cdot \nabla f_2+\cdots =\frac{1}{\text{ϵ}}\sigma \mathcal{L}[f_0+\text{ϵ}{f}_1+{\text{ϵ}}^2{f}_2+\cdots ].$$

Full text

Optical tomographic imaging is a technique for infering the properties of biological tissue via measurements of the incoming and outgoing light intensity; it may be used as a medical imaging methodology. Mathematically, light propagation is modeled by the radiative transfer equation (RTE), and optical tomography amounts to reconstructing the scattering and the absorption coefficients in the RTE from boundary measurements. We study this problem in the Bayesian framework, focussing on the strong scattering regime. In this regime the forward RTE is close to the diffusion equation (DE). We study the RTE in the asymptotic regime where the forward problem approaches the DE, and prove convergence of the inverse RTE to the inverse DE in both nonlinear and linear settings. Convergence is proved by studying the distance between the two posterior distributions using the Hellinger metric, and using Kullback-Leibler divergence.

1 Introduction

1.1 Background

Optical imaging is one of many medical imaging techniques that uses light to probe tissue structure [29, 3]. Near infrared light is sent into biological tissue, and the outgoing photon current at the surface of the tissue is then measured. Using these measurements, it is possible to infer properties of the tissue. While traditional imaging methods such as X-ray imaging provide good reconstructions of the tissue’s properties, they are typically more expensive and more invasive than optical imaging. Optical imaging can be used for brain, breast, and joint imaging, as well as monitoring blood oxygenation [19, 20, 27].

To study optical imaging mathematically, one may use the radiative transfer equation (RTE). The forward RTE describes the dynamics of photons in materials with given optical properties. We denote the distribution of particles at location $x$ with velocity $v$ by $f(x,v)$, where $x\in\Omega\subset\mathbb{R}^{d}$, $d=2,3$, and $v\in\mathbb{S}^{d-1}$, the unit sphere in $\mathbb{R}^{d}$. The model enforces particle motion with constant unit speed, and the velocity affects the problem only through the direction of travel of the particle. The optical properties are characterized by two parameters, the scattering coefficient and the absorption coefficient. The scattering coefficient, denoted $k(x,v,v^{\prime})$, determined by the probability of a photon, currently moving in direction $v$ at position $x$, scattering off a particle in the material and changing direction to $v^{\prime}$. The total absorption coefficient, denoted by $\alpha(x,v)$, accounts for photon absorbtion into the material where they are lost. With this notation established, the RTE is, in its most general form,

[TABLE]

Here the gradient operation is with respect to $x$, as are related contractions of the gradient to a divergence. The forward problem is to determine the particle distribution function $f$, given the optical properties of the medium as characterized by $k$ and $\alpha.$ Optical imaging amounts to solving a related inverse problem: the map from incoming data (light intensity injected into the tissue) to the measured outgoing data (light intensity collected outside the tissue) is termed the albedo operator, and the absorption and scattering coefficients in RTE are reconstructed from knowledge of the albedo operator. There are a number of theoretical results concerning the inverse RTE, primarily focussed on the setting in which the entire albedo map is known: it was shown in [11] that the medium is uniquely recoverable in dimension $d=3$, and then in [33] that the reconstruction is Lipschitz stable. For further literature surveys see the reviews in [4, 5].

Another model for photon dynamics is the diffusion equation (DE). The diffusion equation typically governs lower energy photons than the RTE, leading to a larger scattering effect and less absorption. Let $\rho(x)$ denote the light intensity at location $x$ where, as before, $x\in\Omega\subset\mathbb{R}^{d}$, and let $a(x)$ denote the diffusion coefficient. Then the DE is

[TABLE]

In this setup, the map from the Dirichlet data (light intensity injected into the tissue) to the Neumann data (light propagating out) is termed the Dirichlet-to-Neumann (DtN) map, and is used to reconstruct the diffusion coefficient $a(x)$. Using the DtN map to reconstruct the medium for the elliptic equation is the famous Calderón problem, widely studied from a theoretical perspective. Two foundational papers are [32], where the uniqueness was shown, and [1] in which logarithmic stability of the inversion was demonstrated. The review [8] contains further citations to literature in this area.

It is natural to examine the relation between the two forward models, and to understand, from both a physical as well as a mathematial perspective, why they give distinct stability performances in the inverse problem. When $a(x)$ and $(k(x,v,v^{\prime}),\alpha(x,v))$ satisfy certain relations, the two forward models are asymptotically “close” when the laser beam is composed of low-energy photons. In the forward setting, physically, high-energy photons experience little scattering before leaving the domain, while in comparison, low-energy photons are scattered frequently by the tissue before being released and measured at the boundary. As a consequence, high-energy photons present a crisp resolution, and the images from low-energy photons are rather blurred. Mathematically, the RTE is taken as the correct forward model, and we can use the Knudsen number to present the number of times that an average photon scatters. In the low-energy regime, the number of times the photon scatters increases, the Knudsen number shrinks to zero, and the RTE converges to the DE in the forward setting. Correspondingly, the inverse RTE is expected to converge to the inverse DE: meaning the information carried in the albedo operator is almost the same as that in the DtN map, and the reconstruction should also converge. This has been numerically observed in [2, 17, 10], and proved rigorously in [21].

The literature refered to thus far focusses on settings in which the entire albedo or DtN map is known, and this leads to a deep mathematical theory. However it is arguably far from the practical setting in which partial and noisy information about these mappings is all that is available. The Bayesian formulation of the inverse problem is useful in this setting as it allows for incorporation of prior information, partial obeservation and noise level in a natural fashion. The practicality of this approach was demonstrated in the monograph [18], and recent work has led to a mathematical framework [22, 23, 31] suitable for well-posedness analyses [15] and computations which blend state-of-the-art computational PDE and computatioal statistical approaches [9, 13]. In the Bayesian approach to the solution of the inverse problem, all quantities are viewed as random variables, and the solution is the probability distribution of the unknown quantity conditioned on the data [15]. Bayes’ theorem allows determination of this conditional distribution (the posterior) from the prior distribution on the unknown and from the likelihood, the distribution of the data conditioned on fixing the unknown. Our work is focussed on understanding the relationship between the two inverse problems in this Bayesian setting.

1.2 Our Contribution

The goal of this paper is to connect the two inverse problems in optical imaging, and specifically to prove convergence of the inverse RTE to the inverse DE in the Bayesian framework. Multiscale techniques provide the desired estimates on the forward problem and we show how these may be transfered to the Bayesian inverse problems. To this end we make the following contributions:

• •

we provide multiscale-based error estimates which relate the solution of the forward problems for the RTE and DE – see Theorem 2 for which we provide a formal asymptotic justification in the main body of the paper, and a rigorous proof in the appendix;

• •

we compare the two posterior distributions for the RTE and DE, measuring distance between them in the Kullback-Leibler divergence and the Hellinger distance in the optically thick regime (zero limit of the Knudsen number) when the scattering coefficient is large – see Theorems 2 and 3;

• •

we extend the convergence result linking posterior distributions to a setting in which the albedo operator’s dependence on the medium is approximated by linearization – see Theorem 3 and Corollary 2.

The approach we adopt will apply to other Bayesian inverse problems whose forward models are linked through multiscale analyses. The paper is organized as follows. In section 2 we provide the mathematical setting for our work, including discussion of the Bayesian formulation of inverse problems, and the diffusion limit of the radiative transfer equation. In section 3 we estimate the distance between the Bayesian solution of the RTE and DE inverse problems, and in section 4 we address the same question in the linearized setting. We conclude in section 5.

2 The Setting

In this section we establish the mathematical framework within which all our results are derived. In subsection 2.1 we describe Bayesian inverse problems in general. We then discuss the setting of linear inverse problems with Gaussian priors and Gaussian additive noise, in which the posterior is also Gaussian; and we discuss linearization of the forward operator to obtain an approximate Gaussian posterior. Subsection 2.2 describes the forward problems from the RTE and for the DE, providing error estimates linking their solutions in the small Knudsen number regime. In subsection 2.3 we formulate the Bayesian inverse problem for the RTE and DE. Subsection 2.4 is denoted to linearization of the forward operator for RTE and DE, and hence forms the basis for defining an approximate Gaussian posterior distributions for their respective inverse problems.

2.1 Bayesian Inversion

Consider the inverse problem of finding $\sigma$ from $y$ where

[TABLE]

$\mathcal{G}$ is a known forward map that takes parameter to the data space, and $\eta$ a noise pollution. In the Bayesian formulation of inversion $y,\sigma$ and $\eta$ are viewed as random variables, linked by equation (2), and assumption that the prior distribution on $\sigma$, denoted by $\mu_{0}$, and the distribution of the noise $\eta$, denoted by $\mu_{\text{error}}$, are known. The objective is to find the conditional distribution on $\sigma$ given $y$ (denoted $\sigma|y$.) In this paper we will assume that $\eta$ is independent of $\sigma$ a priori and denote the distribution of $y$ given $\sigma$, which is then a translation by $\mathcal{G}(\sigma)$ of $\mu_{\text{error}}(\eta)$, by $\mu^{\sigma}(y).$ We will concentrate on the commonly occuring case in which $\eta$ is in a function space and the data $y$ is finite dimensional; then $\mu_{\text{error}}(\eta)$ may be identified with its Lebsegue density, whilst $\mu_{0}$ and $\mu^{y}$ are measures on a separable Banach space.

If we denote by $\mu^{y}$ the posterior distribution on $\sigma$ given observation $y$ then Bayes’ theorem gives

[TABLE]

where

[TABLE]

and $X$ is a subset of a separable Banach space which contains the support of $\mu_{0}$; then $Z$ normalizes $\mu^{y}$ to a probability density. If we make the additional assumption that both $\mu_{0}$ and $\mu_{\text{error}}(\eta)$ are Gaussian and finite dimensional so that $\mu_{0}=\mathcal{N}(m_{0},\mathcal{C}_{\text{prior}})$ and $\mu_{\text{error}}(\eta)=N(0,\mathcal{C}_{\text{prior}})$ then we may write a formula for the Lebesgue density of the posterior:

[TABLE]

We note that analogous formulae are also available in the infinite dimensional case; see [31, 15] and the references therein.

In optical tomography, one has two fundamental models for describing light propagation: the radiative transfer equation (RTE), and the diffusion equation (DE). We will denote the solution of the respective inverse problems by $\mu^{y}_{\textrm{DE}}(\sigma)$ and $\mu^{y}_{\textrm{RTE}}(\sigma)$. This paper is primarily concerned with showing that these two distributions are close in the small Knudsen number regime, and quantifying the difference. There are multiple ways to quantify the distance between two probability distributions $\mu$ and $\mu^{\prime}$. We will use the Kullback-Leibler (KL) divergence and the Hellinger distance. If $\mu$ has density with respect to $\mu^{\prime}$ and $\mu$ has support in $X$ defined as above, then the KL divergence is given by

[TABLE]

if $\mu$ and $\mu^{\prime}$ have density with respect to common reference measure $\lambda$, with support in $X$ defined as above, then the Hellinger distance is given by

[TABLE]

These formulae have interpretations in the infinite dimensional setting; see the appendix of [15]. The KL divergence has an information theoretic interpretation which makes it attractive. However, the Hellinger metric is particularly useful because, for square integrable test functions, it translates directly into bounds of differences of expectations of test functions, see Lemma 7.14 in [15]. The square root of the KL divergence bounds the Hellinger metric, but often sharper bounds on differences of expectations of test functions are obtained by studying the Hellinger distance directly. KL divergence was used to quantify the error incurred when approximating posterior distributions in [26] in finite dimensions, and the Hellinger metric was used in [12] in the infinite dimensional setting required in this paper.

In some contexts the unknown $\sigma$ is naturally a positive random variable and so we seek instead $u$ where $\sigma=\exp(u)$. If we define $\mathcal{F}=\mathcal{G}\circ\exp(\cdot)$ then the inverse problem (2) becomes

[TABLE]

Often we have an approximate solution $u_{0}$ to (5) and it is natural to seek a solution which deviates slightly from this. In this situation we write $u=u_{0}+v$ and linearize (5) to obtain

[TABLE]

This suggests studying the linear inverse problem

[TABLE]

where $z=y-\mathcal{F}(u_{0})$ and $G=D\mathcal{F}(u_{0}).$ If we put Gaussian prior $\mathcal{N}(m_{\text{prior}},\mathcal{C}_{\text{prior}})$ on $v$ then the posterior on $v|z$ is also Gaussian $\mathcal{N}(m_{\text{post}},\mathcal{C}_{\text{post}})$ determined by

[TABLE]

These formulae can also be interpreted in the infinite dimensional setting; see [24] and further citations in [31, 15, 28].

When Bayesian inversion is based on a nonlinear forward model, characterization of the resulting non-Gaussian posterior distribution can be quite complicated, requiring MCMC or SMC techniques [30]. One possible approach to deal with this is to perform the linearization described above and work with Gaussian priors and posterior distributions, leading to closed form solutions. These can be augmented with constraints by means of rejection sampling based on independent sampling from the Gaussian posterior.

2.2 Diffusion Limit Of The RTE

We consider the RTE (1) in the setting where the the absorption coefficient $\alpha$ and the scattering coefficient $k(x,v,v^{\prime})$ are set to

[TABLE]

where $\epsilon$ is the Knudsen number. The thickness of the material physically corresponds to the number of times a photon scatters between being injected in a medium and escaping. The physical quantity is termed the Knudsen number, which stands for the ratio of mean free path and the domain length. The mean free path is the average distance a particle travels before being scattered. When the Knudsen number is small, photons, on average, scatter many times before they are emitted, and the material is thus regarded as optically thick.

In equation (1) $\textrm{d}v$ denotes normalized unit measure, meaning

[TABLE]

where we have used the notation $\braket{\cdot}_{v}$ to denote normalized integration over $v$. Thus equation (1) may be written as

[TABLE]

where the collision operator is

[TABLE]

To ensure a unique solution we impose an incoming boundary condition, the analogue of a Dirichlet boundary condition for equations lacking velocity dependence. To this end define

[TABLE]

which denotes the collection of coordinates on the boundary $x\in\partial\Omega$ on which the velocity $v$ points in/out of the domain, respectively where $\pm v\cdot n_{x}>0$. Here $n_{x}$ is the normal vector at point $x$ pointing out of $\Omega$. The incoming boundary condition is imposed on $\Gamma_{-}$. We also define, for any fixed $y\in\partial\Omega$,

[TABLE]

and set

[TABLE]

For a unique solution to (1), boundary conditions must be imposed on $\Gamma_{-}$ as follows:

[TABLE]

Combining the foregoing considerations we obtain

[TABLE]

with $\mathcal{L}$ as defined in (8). The domain $\Omega$ has a smooth $C^{1}$ boundary $\partial\Omega$. In the small $\epsilon$ regime, it was conjectured in [7] and then proved in [6] that the equation is asymptotically close to the following diffusion equation:

[TABLE]

We make this convergence explicit under the following assumptions.

Assumption 2.1**.**

The functions $\sigma,\phi$ and $\xi$ characterizing the medium and the boundary conditions are smooth functions, bounded in the following sense:

• •

the admissible medium is bounded, meaning there is a constant $C_{1}$ so that:

[TABLE]

• •

and the boundary conditions are smooth and bounded, meaning:

[TABLE]

We also term the set of admissible media:

[TABLE]

Here $C^{3}$ is the collection of third-order differentiable function set.

With this assumption, we first have the uniform boundedness of the Neumann data over $\mathcal{A}$.

Proposition 2.1** ([16]).**

Suppose $\rho$ solves (10) with medium $\sigma$ and the smooth boundary condition $\xi$ satisfying Assumption 2.1, then there is a constant $C$ that only depends on $C_{1}$ and $\Omega$, so that

[TABLE]

Note that we assume only that the medium is smooth enough so the Neumann data is bounded. The regularity of the medium could certainly be relaxed but we do not pursue that direction in this paper. The key point here is to have the uniform bound over the set $\mathcal{A}$.

Theorem 2**.**

Suppose $f(x,v)$ satisfies equation (9) with smooth boundary conditions and $\rho(x)$ solves (10). Then, as $\epsilon\to 0$, $f(x,v)\to\rho(x),$ assuming suitable compatibility relationships between the boundary data $\phi$ and $\xi$ of the two equations. In particular, with compatible boundary conditions at different orders, one approximates $f$ through different forms:

• •

if $\phi(x,v)=\xi(x)$:

[TABLE]

• •

if $\phi(x,v)=\xi(x)-\epsilon\frac{1}{\sigma(x)}v(x)\cdot\nabla\xi(x)$:

[TABLE]

Here the constant $C_{\mathcal{A}}$ depends on $C_{1}$, the upper bound in Assumption 2.1 for the admissible set.

We leave the rigorous proof to the appendix, and present here the formal perturbation expansion derivation; the latter is useful in building intuition.

Proof 2.2** (Sketch Proof: Perturbation Expansion).**

We will use the standard asymptotic expansion technique in $\epsilon$ away from the boundary. Set

[TABLE]

Plugging the expansion into (9), we obtain

[TABLE]

Multiplying by $\epsilon$ and equating in powers of $\epsilon$ gives

[TABLE]

The zeroth order expansion indicates that $f_{0}$ is in the null space of $\mathcal{L}$. From equation (8) we deduce that $f_{0}$ must be velocity independent, and thus we write $f_{0}(x,v)=\rho(x)$. With this expression, considering the equation at $\mathcal{O}(\epsilon)$, we have

[TABLE]

Here we have used the fact that $\mathcal{L}$ is one-to-one on the domain ($\textrm{Null}\,\mathcal{L})^{\perp}$ and that $v\cdot\nabla\rho\perp(\textrm{Null}\,\mathcal{L})^{\perp}$, since $v$ integrates to zero on the unit sphere. To close the system we consider the equation at $\mathcal{O}(\epsilon^{2})$, substituting $f_{0}=\rho$ and $f_{1}=-\frac{1}{\sigma}v\cdot\nabla\rho$ to obtain

[TABLE]

Integrating the equation on both sides with respect to $v$ and taking into account the fact that $\int_{\mathbb{S}^{d-1}}\mathcal{L}[f_{2}]\textrm{d}v=0$, we have, using the summation convention,

[TABLE]

implying that

[TABLE]

Here we have used that

[TABLE]

Note that $C_{d}$ depends on dimension. For example, in $\mathbb{S}^{2}$, $C_{d}=1/3$. Thus, we have shown that the radiative transfer equation in the diffusion limit becomes the diffusion equation, which concludes the sketch proof of the theorem.

Remark 2.3**.**

In the preceding formal derivation we have ignored boundary conditions. In practice, unless these are chosen carefully, there will be a mismatch between the DE and the small $\epsilon$ solution of the RTE near the boundary. The boundary conditions stated in the theorem give different levels of consistency between the two equations, and hence lead to differing error estimates. See the proof in the appendix for details. When the boundary conditions are incompatible the analysis is considerably more subtle – see [7, 34, 25] for details.

2.3 Inverse Problems for the RTE and DE

We now define Bayesian inverse problems for the RTE and DE, relating to partial and noisy observations of the albedo and DtN operators respectively. The first ingredient is definition of the forward map $\mathcal{G}$, which we now do for the RTE and DE equations. We conclude this subsection with a discussion of the prior distribution, which we choose in common between the RTE and DE settings.

In optical tomography, high energy light with a known intensity is injected into the material, and detectors are placed on the tissue boundary to collect the light current emitted. For the RTE, the albedo operator is defined by $\mathcal{H}^{\textrm{RTE}}$ which is a $\sigma-$dependent linear transformation of boundary data $\phi$ into the measurement space, defined by

[TABLE]

where

[TABLE]

and $f$ satisfies (9). It is important to note that, whilst $\mathcal{H}^{\textrm{RTE}}$ is linear in its action on $\phi$, it depends nonlinearly on the unknown medium $\sigma.$ The inverse problem of reconstructing $\sigma$ from measurements of $h^{\textrm{RTE}}$ is thus a nonlinear inverse problem.

In practice, finitely many smooth incoming data $\phi_{k}$ are injected and finitely many measurements are made at the boundary for each $\phi_{k}$; we assume that these measurements may be expressed as linear functionals $l_{j}$ of $h^{\textrm{RTE}}$. We thus define the forward map to be inverted by

[TABLE]

where $(j,k)\in\{1,\cdots,J\}\otimes\{1,\cdots,K\}.$ We assume additive Gaussian noise $\eta$ to obtain the compact representation of the inverse problem

[TABLE]

where $\eta\in\mathbb{R}^{JK}$ is drawn from a Gaussian distribution which we assume to have the form

[TABLE]

meaning

[TABLE]

For the DE model the situation is analogous. The map that takes the Dirichlet data to the Neumann outflow is termed the DtN map and is defined by

[TABLE]

where

[TABLE]

and $\rho$ satisfies (10). In practice, finitely many incoming data $\xi_{k}$ are injected and finitely many linear functionals $l_{j}$ of $h^{\textrm{DE}}$ are measured, noisily, leading to an inverse problem of the form

[TABLE]

where $\eta\in\mathbb{R}^{JK}$ denotes observational noise and where the forward map is defined by

[TABLE]

where $(j,k)\in\{1,\cdots,J\}\otimes\{1,\cdots,K\}.$ For simplicity we assume the same noise model (18) for $\eta.$

Together with (18) and the assumption that $\eta$ and $\sigma$ are a priori independent, (17) and (22) define the likelihood for a Bayesian formulation of the inverse problem of determining $\sigma$ from $\mathsf{y}$ from RTE and DE respectively. We now define the prior on $\sigma$, which we will choose in common between the two inverse problems. To this end recall the set (11) and define prior distribution $\mu_{0}(d\sigma)$ to be a probability measure supported on $\mathcal{A}:$

Assumption 3**.**

The prior measure $\mu_{0}$ is supported on an infinite dimensional separable Banach space, and the support is contained in the admissible set $\mathcal{A}$ given by (11):

[TABLE]

In our case $\mathcal{A}$ is a subset of $C^{3}$. If we further relax the regularity assumptions, to let $\mathcal{A}$ be a subset of $W^{1,\infty}$, for example, then $W^{1,\infty}$ is not separable. But it is possible to construct useful measures with support in $W^{1,\infty}$ which are separable, for example through the closure of sets of random series expansions; see [15] for a related example in $L^{\infty}.$

Bayes’ theorem (3) for both models is then given by

[TABLE]

where

[TABLE]

The functions $\mathcal{G}^{\textrm{RTE}}$ and $\mathcal{G}^{\textrm{DE}}$ are here both viewed as mappings from $\mathcal{A}$ into $\mathbb{R}^{JK}.$ The normalization factors are given by

[TABLE]

Note also that the likelihoods $\mu^{\sigma}_{\textrm{RTE}}(\mathsf{y})$ and $\mu^{\sigma}_{\textrm{DE}}(\mathsf{y})$ are, for fixed $\sigma$, proportional to densities on $\mathbb{R}^{JK}$; hence we write them as functions of $y$. On the other hand $\mu^{\mathsf{y}}_{\textrm{RTE}}(d\sigma),\mu^{\mathsf{y}}_{\textrm{DE}}(d\sigma)$ and $\mu_{0}(d\sigma)$ are measures with support in $\mathcal{A}$, a subset of an infinite dimensional separable Banach space.

Theorem 2 shows that given compatible $\phi$ and $\xi$, $h^{\textrm{RTE}}$ and $h^{\textrm{DE}}$ are close when the Knudsen number $\epsilon$ is small, so that $\mathcal{G}^{\textrm{RTE}}$ and $\mathcal{G}^{\textrm{DE}}$ are close for every fixed $\sigma$, when $\epsilon$ is small. In section 3 we use these facts to demonstrate the convergence of $\mu^{\sigma}_{\textrm{RTE}}$ to $\mu^{\sigma}_{\textrm{DE}}$ as $\epsilon\to 0.$

2.4 Linearized Albedo Operator And DtN Map

We derive linearized versions of the albedo operator and the DtN map by assuming that the unknown medium $\sigma$ is close to a known background $\sigma_{0}$. In order to enforce positivity, we assume $\sigma=e^{u}$, define $u_{0}$ so that $\sigma_{0}=e^{u_{0}}$, and find equations satisfied by the perturbation $w=u-u_{0}$. The corresponding inverse problem amounts to reconstructing $w$ using the measurements and some known information computed using the background medium $\sigma_{0}$. We express the admissible set $\mathcal{A}$ given in (11) on the log-scale and write

[TABLE]

This implies

[TABLE]

To start, we recall equation (9),

[TABLE]

We assume that there is a background scattering coefficient $u_{0}\in\mathcal{A}_{u}$, and that $w(x)\in C^{3}(\Omega)$ is a small fluctuation of $u$ around the background $u_{0}$: $w=u-u_{0}$. Then

[TABLE]

We define a new function $f_{\textrm{lin}}$ which solves the RTE with the background scattering coefficient and the same boundary condition,

[TABLE]

The difference between $f$ and $f_{\textrm{lin}}$, termed $\mathsf{f}=f-f_{\textrm{lin}}$, then satisfies, neglecting terms of $\mathcal{O}(w^{2})$, the following error equation:

[TABLE]

To extract boundary data from (32), we define the adjoint equation, with a delta function on the boundary,

[TABLE]

Here we have used the fact that $\mathcal{L}$ is self-adjoint, and for the adjoint equation, the incoming boundary condition should be imposed on $\Gamma_{+}$. We have also imposed a delta function concentrated at $y\in\partial\Omega$. Multiplying (32) by $g$ and (33) by $\mathsf{f}$ and subtracting the two obtained equations, integrated over $x$ and $v$, we obtain, upon using Green’s identity

[TABLE]

Noting that $\mathsf{f}|_{\Gamma_{-}}=0$ and $g|_{\Gamma_{+}}=\delta_{y}$, we may further simplify the left hand side, obtaining

[TABLE]

As in the nonlinear case, we have finitely many measurements and experiments conducted. In the $K$ experiments, we use incoming data $\phi_{k}$, and for each experiment we measure data using the measurement-operator $l_{j}$:

[TABLE]

Letting $g_{j}$ denote the solution to (33) with $\delta_{x_{j}}(x)$ on the boundary, and $f_{k}$, and $f_{\textrm{lin},k}$ denote the solutions to (29) and (31) with $\phi_{k}$ as boundary data, we define

[TABLE]

and define

[TABLE]

with $\langle\cdot,\cdot\rangle_{x}$ denoting the inner-product defined by integration over $\Omega$ in the $x$ variable alone. We note that, by  (34),

[TABLE]

because

[TABLE]

Here $\mathsf{f}_{k}=f_{k}-f_{\textrm{lin},k}.$

Deriving the linearized DtN map for the diffusion equation is similar. For ease of notation we start with the form of the diffusion equation as in (10), where the scattering coefficient is shown in the denominator. For positivity, we use $\sigma=e^{u}$, and $\sigma_{0}=e^{u_{0}}$ as before. We now derive an equation for $\tau=\rho-\rho_{lin}$, first noting that

[TABLE]

where $\rho_{\textrm{lin}}$ solves:

[TABLE]

Subtracting the two equations and neglecting higher order terms, we have

[TABLE]

with $w=u-u_{0}$ as in the RTE case.

We also define $\rho_{g}$ that solves the adjoint equation:

[TABLE]

Multiplying (40) by $\rho_{g}$, (41) by $\tau$ and integrating over $\Omega$, we obtain:

[TABLE]

Similarly to the nonlinear case, we conduct finitely many experiments and make finitely many measurements as in (35). Defining

[TABLE]

where $\rho_{g,j}$ satisfies (41) with $\delta_{y_{j}}$ boundary condition and $\rho_{\textrm{lin},k}$ satisfies (39) with $\xi_{k}$ with as the boundary condition. Using the approximation that $\tau$ satisfies (40), we write (42) as

[TABLE]

where again $G^{\textrm{DE}}_{jk}$ is the linearized approximation and we have used the estimate

[TABLE]

Here $\tau_{k}=\rho_{k}-\rho_{lin,k}.$ This defines the (linear) action of $G^{DE}_{jk}$ on $w$. Notice that the linearized albedo operator (37) and the linearized DtN map (44) have the same format: they are both Fredholm first type integrals, determined by the $\gamma_{jk}^{\textrm{RTE}}$ and $\gamma_{jk}^{\textrm{DE}}$ respectively defined in (36) and (43). To show the convergence of the two problems amounts to showing, in the small $\epsilon$ regime, the convergence of the two forward maps $\gamma_{jk}^{\textrm{RTE}}\sim\gamma_{jk}^{\textrm{DE}}$ for all $j$ and $k$, and the convergence of the data $G^{\textrm{RTE}}_{jk}(w)$ to $G^{\textrm{DE}}_{jk}(w)$ for reasonably small $w$.

3 Nonlinear Inverse Problems

In this section we analyze the distance between the posterior distributions of the nonlinear inverse problems defined by the RTE and DE, namely $\mu^{\mathsf{y}}_{\textrm{RTE}}(\sigma)$ and $\mu^{\mathsf{y}}_{\textrm{DE}}(\sigma)$ respectively. We consider the setting in which the Knudsen number $\epsilon$ is small. We show that the two distributions converge in the Kullback-Leibler divergence and the Hellinger distance as $\epsilon\to 0$. The three subsections concern, in turn, the following convergence results as $\epsilon\to 0$:

convergence of the forward map $\mathcal{G}^{\textrm{RTE}}(\sigma)$ to $\mathcal{G}^{\textrm{DE}}(\sigma)$ for a fixed list of $(\phi_{k},l_{j})$;

• 2.

convergence of the KL divergence between $\mu^{\mathsf{y}}_{\textrm{RTE}}(\sigma)$ and $\mu^{\mathsf{y}}_{\textrm{DE}}(\sigma)$ to zero;

• 3.

convergence of the Hellinger metric $\mu^{\mathsf{y}}_{\textrm{RTE}}(\sigma)$ and $\mu^{\mathsf{y}}_{\textrm{DE}}(\sigma)$ to zero.

Before these three pieces of analysis, recall that the forward measurement operators for the RTE and DE are defined in (16) and (23) respectively and that Bayes’ theorem (3) delivers the formulae for the posterior distributions in (24)–(26).

3.1 Convergence Of The Forward Map

For simplicity we assume that $l_{j}$ is the linear funtional that takes corresponding to evaluation at point $x_{j}\in\partial\Omega$:

[TABLE]

Other linear functionals can be handled with similar analysis. Then

[TABLE]

and

[TABLE]

where $f$ and $\rho$ are the solutions to the RTE and the DE with $\phi_{k}$ and $\xi_{k}$ as incoming conditions, respectively.

We now have the following proposition:

Proposition 1**.**

Assume that $\phi_{k}(x,v)=\xi(x)-\epsilon\frac{1}{\sigma(x)}v(x)\cdot\nabla\xi_{k}(x)$. Then, under Assumption 2.1, the forward maps $\mathcal{G}^{\textrm{RTE}}$ and $\mathcal{G}^{\textrm{DE}}$ satisfy

[TABLE]

Furthermore, there is a constant $C$ that only depends on $C_{1}$ and $\Omega$ so that

[TABLE]

Proof 3.1**.**

To show the first item it is enough to prove that for every $j$ and $k$,

[TABLE]

Note that, for any $y$, ${\mathbb{S}}^{d-1}\backslash\Gamma(y)$ is the set on which $(v\cdot n)=0$. Hence, employing (13), (46) and (47), and defining

[TABLE]

we then have

[TABLE]

Here we used Theorem 2 which states

[TABLE]

which concludes the proof of equation (48). Equation (49) is a direct consequence of Proposition 2.1 and by combining the conclusion of Proposition 2.1 with equation (48).

3.2 Convergence In Kullback-Leibler Divergence

We use the convergence of the forward map to show the convergence in the posterior distribution using the Kullback-Leibler divergence.

Theorem 2**.**

Let the assumptions of Proposition 1, together with Assumption 3, hold. Then

[TABLE]

Proof 3.2**.**

We first note that, over the set $\mathcal{A}$ which contains the support of the common prior measure $\mu_{0}$, the likelihoods $\mu^{\sigma}_{\textrm{RTE}}(\mathsf{y})$ and $\mu^{\sigma}_{\textrm{DE}}(\mathsf{y})$ are bounded uniformly from above and below. Hence the measures $\mu^{\mathsf{y}}_{\textrm{RTE}}$ and $\mu^{\mathsf{y}}_{\textrm{DE}}$ are mutually absolutely continuous and have densities with respect to one another. In particular we may define

[TABLE]

where $\sigma\in\mathcal{A}$. Clearly $\mu^{\mathsf{y}}_{\textrm{DE}}$ has no $\epsilon$ dependence, and so it suffices to show that $\log\frac{d\mu^{\mathsf{y}}_{\textrm{RTE}}}{d\mu^{\mathsf{y}}_{\textrm{DE}}}$ is $\mathcal{O}(\epsilon)$, uniformly on $\mathcal{A}$. Using (24)–(26), we find

[TABLE]

We will show that both terms are $\mathcal{O}(\epsilon)$. Recalling (25) we see that

[TABLE]

where $c<\infty$ is the Lipschitz constant for $\exp(-|x|/2\gamma^{2}).$ Now note that

[TABLE]

and according to Proposition 1,

[TABLE]

we deduce that

[TABLE]

Using the definition of the normalization factor and noting that $\int\textrm{d}\mu_{0}(\mathcal{A})=1$, we also have

[TABLE]

Noting that $Z^{\textrm{DE}}$ and $\mu^{\sigma}_{\textrm{DE}}(\mathsf{y})$ are bounded from below uniformly with respect to $\sigma\in\mathcal{A}$, we deduce from the two preceding displays that

[TABLE]

which completes the proof.

3.3 Convergence In Hellinger Metric

Convergence in the Hellinger metric has a very similar proof to that used in KL divergence.

Theorem 3**.**

Let the assumptions of Proposition 1, together with Assumption 3, hold. Then

[TABLE]

Proof 3.3**.**

We first recall the definition of the Hellinger distance between two distributions in section 2.1, using $\lambda=\mu_{0}$ as the reference measure:

[TABLE]

Following [31] we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

With the same argument, we have

[TABLE]

where we have used

[TABLE]

and the Lipschitz argument as in (52). To deal with $I_{2}$, we notice that

[TABLE]

using the fact that

[TABLE]

According to the definition of $Z^{\textrm{RTE},\textrm{DE}}$, we have

[TABLE]

where we used (52). Plugging these back in (55), we have

[TABLE]

Together with the boundedness of $I_{1}$ and the inequality (53), we conclude

[TABLE]

4 Linearized Inverse Problems

In this section we study approximations of the two Bayesian inverse problems in the linearized setting. We show asymptotic closeness of the posterior distributions in the small Knudsen number regime $\epsilon\ll 1.$ Equations (36)–(38) and (44) give rise to the following approximate inverse problems:

[TABLE]

where

[TABLE]

is a vector of length $JK$ and can be regarded as the linearized data. It can be obtained by subtracting $\mathsf{y}$, the collected measurements with $\{\phi_{k},\,k=1,\dots,K\}$ being the input data and $\{l_{j},\,j=1,\dots,J\}$ being the pointwise evaluation operator, as defined in (45), and $\mathcal{G}^{\textrm{RTE},\textrm{DE}}(\log(u_{0}))$, the background data that is precomputed using (31) or (39) with the same input and measurement operator, and the background medium $u_{0}$.

Assuming $\eta\sim\mathcal{N}(0,\gamma^{2}\mathbb{I})$ as always, now we have the likelihood functions:

[TABLE]

The two models use the same prior distribution function $\nu_{0}(w)$, satisfying

[TABLE]

The posterior distributions are then

[TABLE]

where $Z^{\textrm{RTE}}$ and $Z^{\textrm{DE}}$ are the normalization factors.

4.1 Convergence Of Linearized Forward Map

We first show the convergence of the interpreters $\gamma^{\textrm{RTE},\textrm{DE}}$, which will allow us to show the convergence of the two forward maps.

Proposition 1**.**

Assume $u_{0}\in\mathcal{A}_{u}$, then for $\epsilon$ sufficiently small, $\gamma^{\textrm{RTE}}\to\gamma^{\textrm{DE}}$. More specifically, for every $j$ and $k$,

[TABLE]

Here the constant $C$ only depends on $C_{\mathcal{A}}$ and $C_{1}$.

Proof 4.1**.**

Recall the definition of $\gamma_{jk}$ in (36)

[TABLE]

where $g_{j}$ and $f_{\textrm{lin},k}$ solve (33) and (31) with $\delta_{y_{j}}$ and $\phi_{k}$ as boundary conditions. We further recall Theorem 2, so that we have

[TABLE]

where $\rho_{g,j}$ and $\rho_{\textrm{lin},k}$ solve

[TABLE]

with boundary condition $\delta_{y_{j}}$ and $\xi_{k}$ respectively. Recalling $\mathcal{L}\rho=0$ for all $\rho$, and that $\int_{\mathbb{S}^{d-1}}v\textrm{d}{v}=0$, then

[TABLE]

We conclude the proof since this holds for every $j$ and $k$.

We emphasize that $\gamma^{\textrm{RTE}}$ is uniquely determined by $g$ and $f_{\textrm{lin}}$ that solve (33) and (31), and that the two equations depend merely on $u_{0}$. So the convergence holds true as long as $u_{0}\in\mathcal{A}_{u}$, and there is no requirement on $w$.

Corollary 2**.**

For any fixed $u_{0}\in\mathcal{A}_{u}$, assume $w\in C^{3}$, if $\epsilon$ significantly small, then $G^{\textrm{RTE}}\to G^{\textrm{DE}}$. More specifically,

[TABLE]

where $G^{\textrm{RTE},\textrm{DE}}$ are two vectors of length $JK$, and $C$ only depends on $C_{\mathcal{A}}$, $J$ and $K$.

Proof 4.2**.**

According to the definition of $G^{\textrm{RTE},\textrm{DE}}$,

[TABLE]

We conclude using the result from Proposition 1, and that $\|\gamma^{\textrm{RTE}}_{jk}-\gamma^{\textrm{DE}}_{jk}\|_{L_{2}(\Omega)}\lesssim\|\gamma^{\textrm{RTE}}_{jk}-\gamma^{\textrm{DE}}_{jk}\|_{L_{\infty}(\Omega)}$ for:

[TABLE]

4.2 Convergence In Hellinger Metric

The proof of the following result is a straightforward extension of Theorem 3 and hence we only sketch it.

Theorem 3**.**

Consider the linearized setting with ${u}_{0}\in\mathcal{A}_{u}$ and assume that $\nu_{0}$ is a centred Guassian measure supported on the space $C^{3}(\Omega).$ Then the Hellinger distance between the posterior distribution $\nu^{\mathsf{y}}_{\textrm{RTE}}$ and $\nu^{\mathsf{y}}_{\textrm{DE}}$ is bounded by $\mathcal{O}(\epsilon)$:

[TABLE]

Proof 4.3** (Sketch Proof).**

The primary difference of this theorem with Theorem 3 is that the data $\mathsf{y}$ is subtracted by $\mathcal{G}^{\textrm{RTE}}(\log(u_{0}))$ and $\mathcal{G}^{\textrm{DE}}(\log(u_{0}))$, and that the linear operator can be made explicit: $G^{\textrm{RTE},\textrm{DE}}=\langle\gamma^{\textrm{RTE},\textrm{DE}},w\rangle$. As a result, the estimates for $I_{1}$ and $I_{2}$ change accordingly. The proof is rather similar to that for Theorem 3, so we omit the details and only estimate $I_{1}$ here:

[TABLE]

For the second inequality to hold true, we first use the conclusion from Corollary 2, and to bound the second term, we simply use:

[TABLE]

to obtain

[TABLE]

Application of the Fernique theorem [14] shows that we have $\int_{C^{3}}\|w\|^{p}_{L_{2}(\Omega)}\textrm{d}\nu_{0}$ bounded by a constant (independent of $\epsilon$) for any $p$ and that

[TABLE]

The estimate for $I_{2}$ is very similar, and therefore

[TABLE]

which leads to the conclusion of Theorem 3.

Comparing the preceding theorem with Theorem 3, a very interesting phenomenon we immediately see is that the convergence in the linearized setting has a higher rate. This higher rate is a direct consequence of the convergence in $\gamma$ in which the $\mathcal{O}(\epsilon)$ cancels due to the symmetry between the forward model and the adjoint.

4.3 Implications For Posterior Convergence

In the linear setup, if the prior distribution and the likelihood functions are both Gaussian functions, the posterior distribution is also Gaussian, and is thus completely determined by its mean and covariance, or the leading two moments. In our case, $u_{0}\in\mathcal{A}_{u}$, and $w\in C^{3}$, and the prior is supported in $C^{3}(\Omega)$ for $w$. Thus distances between the posterior distributions computed using the RTE and DE can be estimated from distances between means and covariances. To this end, consider the following lemma.

Lemma 4** (Lemma 7.14 from [15]).**

Let $\nu$ and $\nu^{\prime}$ be two probability measures on a separable Banach space $X$. Assume also that $f:X\to E$, where $(E,\|\cdot\|)$ is a separable Banach space, is measurable and has second moments with respect to both $\nu$ and $\nu^{\prime}$. Then

[TABLE]

Furthermore, if $E$ is a separable Hilbert space and $f:X\to E$ also has fourth moments, then

[TABLE]

where the operator norm on $E$ is employed.

When applied in our case, we obtain:

Corollary 5**.**

Let $m^{\textrm{RTE},\textrm{DE}}_{\text{post}}$ and $\mathcal{C}^{\textrm{RTE},\textrm{DE}}_{\text{post}}$ denote the mean function and the covariance operator computed from the posterior distribution of the radiative transfer and diffusion model in the linearized setting. Then

[TABLE]

Here the norm for the mean is the standard norm on $C^{3}(\Omega)$ and the norm for the covariance is the operator norm on $H^{3}(\Omega).$

Proof 4.4**.**

Let $f$ as in the statement of Lemma 4 be the identity map, and the spaces $X$ and $E$ be $C^{3}(\Omega)$ equipped with $L_{2}$ norm, then $f(w)=w$. Let $\nu$ and $\nu^{\prime}$ be the posterior distributions $\nu^{\mathsf{y}}_{\textrm{RTE}}$ and $\nu^{\mathsf{y}}_{\textrm{DE}}$ respectively. Then:

[TABLE]

Since

[TABLE]

and that

[TABLE]

where we have again used $\int_{\mathcal{C}^{3}(\Omega)}\|w\|^{2}_{L_{2}(\Omega)}\textrm{d}\nu^{\mathsf{y}}_{\textrm{RTE}}<C$ using the Fernique theorem, and that $\textrm{d}\nu^{\mathsf{y}}_{\textrm{RTE},\textrm{DE}}\lesssim\textrm{d}\nu_{0}$ [15], we achieve the convergence of the mean function by applying Theorem 3. The same analysis is applied to analyze the covariance.

5 Conclusion

In this paper, we study the inverse problem of diffuse optical tomography to reconstruct the scattering coefficient. Partial and noisy data is assumed, and hence a Bayesian formulation of inversion is natural. Two forward models are described, one employing the radiative transfer equation and the other employing the diffusion equation respectively. Multiscale analysis demonstrates that solutions of the two forward models are close in the optically thick (small Knudsen number) regime, and this allows us to quantify the convergence of the two Bayesian inverse problems. In particular, we show that $\mu^{\mathsf{y}}_{\textrm{RTE}}$ and $\mu^{\mathsf{y}}_{\textrm{DE}}$, the two posterior distribution functions, are $\mathcal{O}(\epsilon)$ away from each other in both Kullback-Leibler divergence sense, and the Hellinger sense, for both nonlinear and linear cases. Forward solution of the diffusion equation is computationally less burdensome than the radiative transfer equation, and the theory justifies using it within Bayesian inversion algorithms where multiple forward model evaluations are required. We have employed a setting in which compatible boundary conditions are used for the two forward models. It would also be of interest to study extensions of this. However the primary stumbling block here is the analysis of the forward problem itself. The approach developed in this paper will apply to other Bayesian inverse problems whose forward problems are close, and can be used to justify inversion algorithms which employ an averaged or (as in this case) homogenized approximate forward model, in order to speed-up computation.

Acknowledgements AMS is supported by US AFOSR grant FA9550-17-1-0185. KN and QL are supported by NSF DMS 1619778, 1750488 and NSF TRIPODS 1740707.

6 Appendix

We give the rigorous proof here for Theorem 2. The two statements are proved in the same way, and for generality, we will only prove the second one, and the proof for the first statement, or even for higher order expansions, are easy extensions.

Proof 6.1**.**

We repeat the RTE with a specially designed boundary condition,

[TABLE]

where $\rho$ satisfies

[TABLE]

Now we decompose the solution to the RTE as

[TABLE]

where $f_{0}=\rho$, $f_{1}=-\frac{1}{\sigma}v\cdot\nabla\rho$, and $f_{2}=\frac{1}{\sigma}\mathcal{L}^{-1}\left[(v\cdot\nabla)\frac{1}{\sigma}(v\cdot\nabla)\rho\right]$. Note that for the definition of $f_{2}$ to hold true, it is necessary that

[TABLE]

which in turn requires

[TABLE]

Since $\rho$ is $\frac{1}{\sigma}$-harmonic with smooth boundary condition $\|\xi\|_{L_{\infty}(\partial\Omega)}<C_{1}$, then by the maximum principle [16],

[TABLE]

Then since $\mathcal{L}^{-1}$ is a bounded operator on $\textrm{Null}\,\mathcal{L}^{\perp}$, we know that both $f_{1}$ and $f_{2}$ are bounded, meaning there is a constant $C_{4}$ that depends on $C_{1}$, $C_{2}$ and $C_{3}$:

[TABLE]

and

[TABLE]

where we used the boundedness of $\|\sigma^{-1}\|_{L_{\infty}(\Omega)}<C_{1}$, $\|\nabla\left(\frac{1}{\sigma}\right)\|_{L_{\infty}(\Omega)}<C_{1}$, and the boundedness of the harmonic function and its derivatives.

Plugging it back into the equation, we have

[TABLE]

Since $\rho$ is a constant in $v$, and is thus in $\textrm{Null}\,\mathcal{L}$, then $\mathcal{L}\rho=0$. Using the definition of $\mathcal{L}$, we also have $\mathcal{L}(v\cdot\nabla\rho)=-v\cdot\nabla\rho$. Now we cancel the terms and obtain the following equation for $f_{r}$,

[TABLE]

It is immediate that $f_{r}$ satisfies RTE and is equipped with a source term of $\mathcal{O}(\epsilon^{2})$. The boundary condition for $f_{r}$ is of the same order,

[TABLE]

By maximum principle, the solution is bounded in $L_{\infty}$ by the boundary condition and the source term, so we have

[TABLE]

where $C_{5}\leq e^{C_{1}l}$ and $l$ is the longest radius of the domain. This leads to the fact that:

[TABLE]

We note again the constant merely depends on the boundedness of $C_{1}$, the upper bound of the infinite norm of $\frac{1}{\sigma}$, its derivative and $\xi$.