A Hybrid Inversion Scheme Combining Markov Chain Monte Carlo and Iterative Methods for Determining Optical Properties of Random Media

Abstract

Near-infrared spectroscopy (NIRS) including diffuse optical tomography is an imaging modality which makes use of diffuse light propagation in random media. When optical properties of a random medium are investigated from boundary measurements of reflected or transmitted light, iterative inversion schemes such as the Levenberg–Marquardt algorithm are known to fail when initial guesses are not close enough to the true value of the coefficient to be reconstructed. In this paper, we investigate how this weakness of iterative schemes is overcome using Markov chain Monte Carlo. Using time-resolved measurements performed against a polyurethane-based phantom, we present a case that the Levenberg–Marquardt algorithm fails to work but the proposed hybrid method works well. Then, with a toy model of diffuse optical tomography we illustrate that the Levenberg–Marquardt method fails when it is trapped by a local minimum but the hybrid method can escape from local minima by using the Metropolis–Hastings Markov chain Monte Carlo algorithm until it reaches the valley of the global minimum. The proposed hybrid scheme can be applied to different inverse problems in NIRS which are solved iteratively. We find that for both numerical and phantom experiments, optical properties such as the absorption and reduced scattering coefficients can be retrieved without being trapped by a local minimum when Monte Carlo simulation is run only about 100 steps before switching to an iterative method. The hybrid method is compared with simulated annealing. Although the Metropolis–Hastings MCMC arrives at the steady state at about 10,000 Monte Carlo steps, in the hybrid method the Monte Carlo simulation can be stopped way before the burn-in time.

1. Introduction

In near-infrared spectroscopy (NIRS), we estimate optical properties of biological tissue by solving inverse diffuse problems [1,2]. Such inverse problems are commonly solved by means of iterative methods. In the case of a homogeneous medium, absorption and reduced scattering coefficients of the medium can be obtained with iterative methods such as the Levenberg–Marquardt algorithm [3,4] from time-resolved measurements (for example, see the review article [5]). Neuroimaging for the human brain by NIRS through the neurovascular coupling has been developed and is called functional NIRS (fNIRS) [6,7,8,9]. Since the region of interest on the head can be small, it is possible to assume a simple geometry of the half space. Quantitative measurements of inter-regional differences in neuronal activity requires accurate estimates of optical properties.

In Ref. [10], the Nelder–Mead simplex method was used to retrieve optical parameters in layered tissue. Heterogeneity of optical properties can be obtained by diffuse optical tomography [1,11,12]. Iterative numerical schemes are used to minimize the cost function when solving these inverse problems. A gradient-based approach [13] was used to detect breast cancer [14]. The brain activity of a newborn infant was investigated [15] with diffuse optical tomography by TOAST (temporal optical absorption and scattering tomography) [16,17], in which iterative algorithms such as the nonlinear conjugate gradients, damped Gauss–Newton method, and Levenberg–Marquardt method are implemented. Diffuse optical tomography was performed on human lower legs and a forearm with the algebraic reconstruction algorithm in the framework of the modified generalized pulse spectrum technique [18]. See Refs. [19,20] for numerical techniques for diffuse optical tomography. For these iterative numerical schemes to work, it is important to choose a good initial guess for the initial value of the iteration.

Solving inverse problems by the Bayesian approach has been sought as an alternative way. In Ref. [21], a novel use of the Bayesian approach was considered to take modeling error into account. The Bayesian inversion with the Metropolis–Hastings Markov chain Monte Carlo was used in Refs. [22,23]. The Bayesian approach was used to determine optical parameters of the human head [24]. Although the Markov chain Monte Carlo (MCMC) approach is in principle able to escape from local minima, it is computationally time consuming. Hence, despite the above-mentioned efforts, the use of Monte Carlo for inverse problems in NIRS has been extremely limited.

In this paper, we shed light on Markov chain Monte Carlo once again by combining it with an iterative scheme, and investigate the use of it in NIRS. In particular, we test a hybrid numerical scheme of Markov chain Monte Carlo and an iterative method. In the proposed hybrid scheme, the Markov chain Monte Carlo algorithm is first used to provide an initial guess using jumps in the landscape of the cost function, and then an iterative method is used after the initial large fluctuation. Thus, the proposed method realizes a fast reconstruction while, in the beginning, obtained values at Monte Carlo steps jump around in the landscape of the cost function. The MCMC simulation is necessary only for the first 100 steps. Then the hybrid method starts to use an iterative scheme and reconstructs optical properties by searching the global minimum. The computation time of the iterative scheme is negligible compared with that of the Monte Carlo simulation. Since the Monte Carlo simulation can be stopped way before the burn-in time for the hybrid method, the proposed scheme is at least ten times faster than simulated annealing, which is implemented by the naive use of the Metropolis–Hastings Markov chain Monte Carlo.

The rest of the paper is organized as follows. We develop diffusion theory in Section 2. A polyurethane-based phantom and numerical phantom are also described in Section 2. Section 3 is devoted to experimental and numerical results. Finally, discussion is given in Section 4.

2. Materials and Methods

2.1. Diffusion Theory

2.1.1. Diffuse Light in Three Dimensions

Let us suppose that a random medium occupies the three-dimensional half space. We assume that the random medium is characterized by the diffusion coefficient D and absorption coefficient ${\mu }_a$. We have $D=1/(3{\mu }_s^{\prime })$, where ${\mu }_s^{\prime }$ is the reduced scattering coefficient. Position in the half space ($-\infty <x<\infty$, $-\infty <y<\infty$, $0<z<\infty$) is denoted by $\mathbf{r}=(\mathit{\rho },z)$, $\mathit{\rho }=(x,y)$. Let t denote time. Let c be the speed of light in the medium. We assume an incident beam on the boundary at the origin in the x-y plane. The energy density u of light in the medium obeys the following diffusion equation.

$$\left(\frac{1}{c}\frac{\partial }{\partial t}-D\Delta +{\mu }_a\right)u(\mathbf{r},t)=0,$$

(1)

with the boundary condition

$$-\ell \frac{\partial }{\partial z}u(\mathbf{r},t)+u(\mathbf{r},t)=\delta (x)\delta (y)q(t),$$

(2)

and the initial condition $u(\mathbf{r},0)=0$. The right-hand side of the boundary condition is the incident laser beam which illuminates the medium at the origin $(0,0,0)$ with the temporal profile $q(t)$. In the phantom experiment described below, the incident light enters the phantom in the positive z direction. The extrapolation distance ℓ is a nonnegative constant. If we consider the diffuse surface reflection, we have [25]

$$\ell =2D\frac{1+r_d}{1-r_d},$$

(3)

where the internal reflection $r_d$ is given by [26]$r_d=-1.4399n^{-2}+0.7099n^{-1}+0.6681+0.0636n$. Let us consider the corresponding surface Green’s function $G_s(\mathbf{r},t;{\mathit{\rho }}^{\prime },s)$, which satisfies Equation (1) and the boundary condition

$$\left(-\ell \frac{\partial }{\partial z}+1\right)G_s(\mathbf{r},t;{\mathit{\rho }}^{\prime },s)=\delta (x-x^{\prime })\delta (y-y^{\prime })\delta (t-s),$$

(4)

together with the initial condition $G_s(\mathbf{r},0;{\mathit{\rho }}^{\prime },s)=0$. We obtain [27,28,29,30]

$$\begin{array}{cc}\hfill G_s(\mathbf{r},t;{\mathit{\rho }}^{\prime },s)& =\frac{cD}{\ell }[\frac{2e^{-{\mu }_{a}c(t-s)}}{{\left(4\pi Dc(t-s)\right)}^{3/2}}{e}^{-\frac{{(x-x^{\prime })}^2+{(y-y^{\prime })}^2}{4Dc(t-s)}}{e}^{-\frac{z^2}{4Dc(t-s)}}\hfill \\ & -\frac{e^{-{\mu }_{a}c(t-s)}}{4\pi \ell Dc(t-s)}{e}^{-\frac{{(x-x^{\prime })}^2+{(y-y^{\prime })}^2}{4Dc(t-s)}}{e}^{\frac{z}{\ell }}{e}^{\frac{Dc(t-s)}{{\ell }^2}}\mathrm{erfc}\left(\frac{z+2Dc(t-s)/\ell }{\sqrt{4Dc(t-s)}}\right)],t>s,\hfill \end{array}$$

(5)

with $G_s(\mathbf{r},t;{\mathit{\rho }}^{\prime },s)=0$ for $t\le s$. The complementary error function is defined by $\mathrm{erfc}(x)=(2/\sqrt{\pi }){\int }_x^{\infty }exp(-t^2)dt$. Let ${\mathbf{r}}_d$ be a point in the x-y plane with $|{\mathbf{r}}_d|>0$. We obtain

$$\begin{array}{c}u({\mathbf{r}}_d,t)={\int }_0^t{G}_s({\mathbf{r}}_d,t;0,s)q(s)ds\hfill \\ ={\int }_0^t{e}^{-{\mu }_{a}c(t-s)}{e}^{-\frac{|{\mathbf{r}}_d{|}^2}{4Dc(t-s)}}\frac{2q(s)}{{\left(4\pi Dc(t-s)\right)}^{3/2}}\left[1-\frac{\sqrt{4\pi Dc(t-s)}}{2\ell }{e}^{\frac{Dc(t-s)}{{\ell }^2}}\mathrm{erfc}\left(\frac{\sqrt{Dc(t-s)}}{\ell }\right)\right]ds.\hfill \end{array}$$

(6)

This $u({\mathbf{r}}_d,t)$ in Equation (6) is used for the calculation in Section 3.1.

2.1.2. Diffuse Light in Two Dimensions

For the later purpose of a numerical experiment, we here consider light propagation in two dimensions. Let us suppose that the two-dimensional half space of positive y is occupied with a random medium in the x-y plane. The energy density u of light at position $\mathit{\rho }=(x,y)$ in the medium due to the incident beam on the boundary (the x-axis) at ${\mathit{\rho }}_s^i$ ($i=1,\dots ,M_s$) obeys the following diffusion equation.

$$\left(\frac{1}{c}\frac{\partial }{\partial t}-D\Delta +{\mu }_a(\mathit{\rho })\right)u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=0,$$

(7)

with the boundary condition

$$-\ell \frac{\partial }{\partial y}u(\mathit{\rho },t;{\mathit{\rho }}_s^i)+u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=\delta (x-x_s^i)\delta (t),$$

(8)

and the initial condition $u(\mathit{\rho },0;{\mathit{\rho }}_s^i)=0$. The incident laser beam on the right-hand side of the boundary condition was assumed to be a pulse at the position ${\mathit{\rho }}_s^i=(x_s^i,0)$. We suppose that the diffusion coefficient D is a positive constant but ${\mu }_a$ varies in space.

Below we develop the Rytov approximation. Let us write ${\mu }_a(\mathit{\rho })\ge 0$ as

$${\mu }_a(\mathit{\rho })={\mu }_{a0}+\delta {\mu }_a(\mathit{\rho }),$$

(9)

where ${\mu }_{a0}$ is a constant and the perturbation $\delta {\mu }_a(\mathit{\rho })$ spatially varies. We note the relation,

$$u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)-{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds,$$

(10)

where $u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)$ is the solution to the diffusion Equation (7) with the zeroth-order coefficient replaced by $b_0$. We introduce the Green’s function $G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)$, which satisfies

$$\left(\frac{1}{c}\frac{\partial }{\partial t}-D\Delta +{\mu }_{a0}\right)G=\delta (\mathit{\rho }-{\mathit{\rho }}^{\prime })\delta (t-s),$$

(11)

with the boundary condition $-\ell \frac{\partial G}{\partial y}+G=0$ at $y=0$, and the initial condition $G=0$ at $t=0$. The above relation (10) can be directly verified. By recursive substitution, we obtain u as

$$u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)-{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u_0({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds+O({(\delta {\mu }_a)}^2).$$

(12)

By neglecting higher-order terms assuming that $\delta {\mu }_a$ is small, we arrive at the (first) Born approximation [31], in which u is given by

$$u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)-c{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u_0({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds.$$

(13)

We note that the Green’s function is obtained as [27,28,29,30]

$$\begin{array}{cc}\hfill G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)& =\frac{e^{-{\mu }_{a0}c(t-s)}}{4\pi D(t-s)}{e}^{-\frac{{(x-x^{\prime })}^2}{4Dc(t-s)}}[e^{-\frac{{(y-y^{\prime })}^2}{4Dc(t-s)}}+e^{-\frac{{(y+y^{\prime })}^2}{4Dc(t-s)}}\hfill \\ & -\frac{\sqrt{4\pi Dc(t-s)}}{\ell }{e}^{-\frac{{(y+y^{\prime })}^2}{4Dc(t-s)}}{e}^{{\left(\frac{y+y^{\prime }+2Dc(t-s)/\ell }{\sqrt{4Dc(t-s)}}\right)}^2}\mathrm{erfc}\left(\frac{y+y^{\prime }+2Dc(t-s)/\ell }{\sqrt{4Dc(t-s)}}\right)],\hfill \end{array}$$

(14)

for $t>s$, and otherwise $G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)=0$. Moreover, we obtain

$$u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)=\frac{e^{-{\mu }_{a0}ct}}{2\pi Dt}{e}^{-\frac{{(x-x_s^i)}^2+y^2}{4Dct}}\left[1-\frac{\sqrt{\pi Dct}}{\ell }{e}^{{\left(\frac{y+2Dct/\ell }{\sqrt{4Dct}}\right)}^2}\mathrm{erfc}\left(\frac{y+2Dct/\ell }{\sqrt{4Dct}}\right)\right].$$

(15)

The above expression of $u_0$ is similar to the formula in Equation (6) but does not have a time integral because in this case we assumed the delta function $\delta (t)$ for the incident beam.

To obtain the expression of the Rytov approximation, we introduce ${\psi }_0,{\psi }_1$ as [31]

$$u_0=e^{{\psi }_0},u=e^{{\psi }_0+{\psi }_1}.$$

(16)

By plugging the above expressions of $u,u_0$ into the Born approximation and neglecting terms of $O({(\delta {\mu }_a)}^2)$, we obtain

$$\begin{array}{cc}\hfill e^{{\psi }_1}& =1-e^{-{\psi }_0}{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u_0({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds\hfill \\ & =exp\left[-e^{-{\psi }_0}{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u_0({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds\right].\hfill \end{array}$$

(17)

The (first) Rytov approximation is thus derived as

$$u(\mathit{\rho },t;{\mathit{\rho }}_s^i)=u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)exp\left[-\frac{1}{u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)}{\int }_0^t{\int }_0^{\infty }{\int }_{-\infty }^{\infty }G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)\delta {\mu }_a({\mathit{\rho }}^{\prime })u_0({\mathit{\rho }}^{\prime },s;{\mathit{\rho }}_s^i)d{x}^{\prime }d{y}^{\prime }ds\right].$$

(18)

Let us define

$$\begin{array}{cc}\hfill g(y,t;y^{\prime },s)& =\frac{1}{4\pi D(t-s)}{e}^{-\frac{{(y+y^{\prime })}^2}{4Dc(t-s)}}[1+e^{\frac{{(y+y^{\prime })}^2-{(y-y^{\prime })}^2}{4Dc(t-s)}}-\frac{\sqrt{4\pi Dc(t-s)}}{\ell }\hfill \\ & \times e^{{\left(\frac{y+y^{\prime }}{2\sqrt{Dc(t-s)}}+\frac{\sqrt{Dc(t-s)}}{\ell }\right)}^2}\mathrm{erfc}\left(\frac{y+y^{\prime }}{2\sqrt{Dc(t-s)}}+\frac{\sqrt{Dc(t-s)}}{\ell }\right)].\hfill \end{array}$$

(19)

Then we have

$${\int }_0^{t}G(\mathit{\rho },t;{\mathit{\rho }}^{\prime },s)u_0({\mathit{\rho }}^{\prime },s)ds=e^{-{\mu }_{a0}ct}{\int }_0^t{e}^{-\frac{{(x-x^{\prime })}^2}{4Dc(t-s)}}{e}^{-\frac{{(x^{\prime }-x_s^i)}^2}{4Dcs}}g(y,t;y^{\prime },s)g(y^{\prime },s;0,0)ds.$$

(20)

Therefore, Equation (18) can be rewritten as

$$\begin{array}{cc}\hfill u(\mathit{\rho },t;{\mathit{\rho }}_s^i)& =u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)exp[-\frac{e^{-{\mu }_{a0}ct}}{u_0(\mathit{\rho },t;{\mathit{\rho }}_s^i)}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\delta {\mu }_a({\mathit{\rho }}^{\prime })\hfill \\ & \times \left({\int }_0^t{e}^{-\frac{{(x-x^{\prime })}^2}{4Dc(t-s)}}{e}^{-\frac{{(x^{\prime }-x_s^i)}^2}{4Dcs}}g(y,t;y^{\prime },s)g(y^{\prime },s;0,0)ds\right)d{x}^{\prime }d{y}^{\prime }].\hfill \end{array}$$

(21)

This expression (21) is used to compute the forward data in Section 3.2.

2.2. Inverse Problems by an Iterative Scheme

We suppose that on the surface of a two- or three-dimensional random medium, for each source at ${\mathit{\rho }}_s^i$ or ${\mathbf{r}}_s^i$ ($i=1,\dots ,M_s$) exiting light is detected at ${\mathit{\rho }}_d^j$ or ${\mathbf{r}}_d^j$ ($j=1,\dots ,M_d$) and is measured at times $t^k$ ($k=1,\dots ,M_t$). Let $y_{ijk}$ be measured data whereas u is a solution to the diffusion equation. Let us suppose that u depends on a vector $\mathbf{a}$ which contains unknown parameters. We wish to reconstruct $\mathbf{a}$. In Section 3.1, $\mathbf{a}=({\mu }_a,D)$, and a is a scalar (a one-dimensional vector) in Section 3.2. For example, in the former case the solution u depends on $\mathbf{a}$ since the calculated value of u depends on ${\mu }_a,D$.

Let us introduce vectors $\mathbf{U}$ and $\mathbf{F}(\mathbf{a})$ of dimension $M_s{M}_d{M}_t$ as

$${\left\{\mathbf{U}\right\}}_{ijk}=y_{ijk},{\left\{\mathbf{F}(\mathbf{a})\right\}}_{ijk}=u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i;\mathbf{a})oru({\mathbf{r}}_d^j,t^k;{\mathbf{r}}_s^i;\mathbf{a}),$$

(22)

where we wrote $u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i)=u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i;\mathbf{a})$ and $u({\mathbf{r}}_d^j,t^k;{\mathbf{r}}_s^i)=u({\mathbf{r}}_d^j,t^k;{\mathbf{r}}_s^i;\mathbf{a})$ emphasizing that u depends on $\mathbf{a}$. We find optimal $\mathbf{a}$ by minimizing ${\parallel \mathbf{U}-\mathbf{F}(\mathbf{a})\parallel }_2^2$. Here we particularly consider the Levenberg–Marquardt method [3,4], i.e., the reconstructed value ${\mathbf{a}}^{*}={lim}_{k\to \infty }{\mathbf{a}}_k$ is computed by the iteration given by

$${\mathbf{a}}_{k+1}={\mathbf{a}}_k+{\left[F^{\prime }{({\mathbf{a}}_k)}^T{F}^{\prime }({\mathbf{a}}_k)+\lambda I\right]}^{-1}{F}^{\prime }{({\mathbf{a}}_k)}^T\left(\mathbf{U}-\mathbf{F}({\mathbf{a}}_k)\right),$$

(23)

where $F^{\prime }(\mathbf{a})$ is the Jacobian matrix, which contains derivatives of $\mathbf{F}(\mathbf{a})$ with respect to $\mathbf{a}$, and the parameter $\lambda$ is nonnegative. By modifying the original algorithm according to Ref. [32], our algorithm of the Levenberg–Marquardt method, which we call Algorithm 1, is described below.

Algorithm 1: Levenberg–Marquardt (LM)

Set $k=0$ and $\lambda =1$. Calculate $\mathbf{F}({\mathbf{a}}_k)$ and $F^{\prime }({\mathbf{a}}_k)$. Calculate $S({\mathbf{a}}_k)={\mathbf{d}}^T\mathbf{d}$, where $\mathbf{d}=\mathbf{U}-\mathbf{F}({\mathbf{a}}_k)$. Prepare $A=F^{\prime }{({\mathbf{a}}_k)}^T{F}^{\prime }({\mathbf{a}}_k)$ and $\mathbf{v}=F^{\prime }{({\mathbf{a}}_k)}^T\mathbf{d}$. Find $\delta$ from $(A+\lambda I)\delta =-\mathbf{v}$. Obtain $S({\mathbf{a}}_k+\delta )$ and $R=\frac{S({\mathbf{a}}_k)-S({\mathbf{a}}_k+\delta )}{-{\delta }^T(2\mathbf{v}+A\delta )}$. If $R<0.25$, then set $\nu =10$ (${\alpha }_c<0.1$), $\nu =1/{\alpha }_c$$(0.1\le {\alpha }_c\le 0.5)$, or $\nu =2$ (${\alpha }_c>0.5$), where ${\alpha }_c={\left[2-\left(S({\mathbf{a}}_k+\delta )-S({\mathbf{a}}_k)\right)/{\delta }^T\mathbf{v}\right]}^{-1}$. If $R<0.25$ and $\lambda =0$, set $\lambda =1/\parallel A^{-1}\parallel$ and $\nu =\nu /2$. In the case of $R<0.25$, we set $\lambda =\nu \lambda$. If $R>0.75$, then we set $\lambda =\lambda /2$. If $R>0.75$ and $\lambda <1/\parallel A^{-1}\parallel$, set $\lambda =0$. Otherwise when $0.25\le R\le 0.75$, no update for $\lambda$. If $S({\mathbf{a}}_k+\delta )\ge S({\mathbf{a}}_k)$, then return to Step 5. If $S({\mathbf{a}}_k+\delta )<S({\mathbf{a}}_k)$, set ${\mathbf{a}}_{k+1}={\mathbf{a}}_k+\delta$. Then put $k+1\to k$ and go back to Step 2. Repeat the above procedure until one of the stopping criteria $\parallel \delta \parallel <{10}^{-4}$ and $S<{10}^{-14}$ is fulfilled.

2.3. Inverse Problems by Markov Chain Monte Carlo

For simplicity in this section, we describe the Metropolis–Hastings Markov chain Monte Carlo algorithm (MH-MCMC) assuming a is the scalar appearing in Section 3.2. The extension of the derivation to the general case of vector $\mathbf{a}$ is straightforward.

Suppose that the coefficient ${\mu }_a$ is unknown and ${\mu }_a$ depends on a parameter a. We will reconstruct a within the framework of the Bayesian inversion algorithm [33,34], i.e., we will find the probability distribution $\pi (a|\mathbf{U})$ of a for measured data $\mathbf{U}$. Let $f_{\mathrm{prior}}(a)$ be the prior probability density and $\pi (\mathbf{U}|a)$ be the likelihood density or the conditional probability density of $\mathbf{U}$ for a. Then the joint probability density $\pi (a,\mathbf{U})$ of $a,\mathbf{U}$ is given by $\pi (a,\mathbf{U})=\pi (\mathbf{U}|a)f_{\mathrm{prior}}(a)$. According to the Bayes formula, the conditional probability density $\pi (a|\mathbf{U})$ is given by

$$\pi (a|\mathbf{U})=\frac{L(\mathbf{U}|a)f_{\mathrm{prior}}(a)}{{\int }_{-\infty }^{\infty }L(\mathbf{U}|a)f_{\mathrm{prior}}(a)da},$$

(24)

where $L(\mathbf{U}|a)$ is a function proportional to $\pi (\mathbf{U}|a)$. Assuming Gaussian noise, we put

$$L(\mathbf{U}|a)=e^{-\frac{1}{2{\sigma }_e^2}{\parallel \mathbf{U}-\mathbf{F}(a)\parallel }_2^2}.$$

(25)

In this paper, we simply set $f_{\mathrm{prior}}(a)=1$, i.e., we can say $f_{\mathrm{prior}}(a)\propto {\mathbf{1}}_{[a_{\min },a_{\max }]}(a)$ and the interval $[a_{\min },a_{\max }]$ is large enough so that all a’s appearing in the Markov chain fall into this interval. Here, ${\mathbf{1}}_A(a)$ is the indicator function defined as ${\mathbf{1}}_A(a)=1$ for $a\in A$ and $=0$ for $a\notin A$. General uniform distributions can be used for $f_{\mathrm{prior}}(a)$ if we use the prior-reversible proposal that satisfies $f_{\mathrm{prior}}(a)q(a^{\prime }|a)=f_{\mathrm{prior}}(a^{\prime })q(a|a^{\prime })$ (see below for the proposal distribution $q(a^{\prime }|a)$) [35]. Another possible choice of $f_{\mathrm{prior}}(a)$ is a Gaussian distribution, which turns out to be the Tikhonov regularization term in the cost function.

Using the Metropolis–Hastings algorithm, we can evaluate $\pi (a|\mathbf{U})$ even when the normalization factor is not known and only the relation $\pi (a|\mathbf{U})\propto L(\mathbf{U}|a)f_{\mathrm{prior}}(a)$ is available [33,34]. We can find $\pi (a|\mathbf{U})$ using a sequence $p_1,p_2,\dots$ as

$$\pi (a|\mathbf{U})=\underset{k\to \infty }{lim}{p}_k(a),$$

(26)

where $p_{k+1}(a)$ is obtained from $p_k(a)$ (Markov chain) as

$$p_{k+1}(a^{\prime })={\int }_{-\infty }^{\infty }K(a^{\prime },a)p_k(a)da.$$

(27)

For all $a,a^{\prime }$, the transition kernel satisfies

$$K(a^{\prime },a)\ge 0,{\int }_{-\infty }^{\infty }K(a^{\prime },a)d{a}^{\prime }=1.$$

(28)

Let us write $K(a^{\prime },a)$ as

$$K(a^{\prime },a)=\alpha (a^{\prime },a)q(a^{\prime }|a)+r(a)\delta (a^{\prime }-a),$$

(29)

where $q(a^{\prime }|a)$ is the proposal distribution and

$$r(a)=1-{\int }_{-\infty }^{\infty }\alpha (a^{\prime },a)q(a^{\prime }|a)d{a}^{\prime }.$$

(30)

In the Metropolis–Hastings algorithm we set $\alpha (a^{\prime },a)=min\left\{1,\pi (a^{\prime }|\mathbf{U})q(a|a^{\prime })/\left[\pi (a|\mathbf{U})q(a^{\prime }|a)\right]\right\}$. With this choice of $\alpha (a^{\prime },a)$, the detailed balance is satisfied and $K(a,a^{\prime })\pi (a^{\prime }|\mathbf{U})=K(a^{\prime },a)\pi (a|\mathbf{U})$. A common choice of $q(a^{\prime }|a)$ is the normal distribution, i.e., $q(·|a)=\mathcal{N}(a,{\epsilon }^2)$ with the mean a and the standard deviation $\epsilon >0$. We note that $q(a^{\prime }|a)=q(a|a^{\prime })$ in this case. We have

$${\int }_{-\infty }^{\infty }h(a)\pi (a|\mathbf{U})da=\underset{k_{\max }\to \infty }{lim}\frac{1}{k_{\max }}\sum _{k=1}^{k_{\max }}h(a_k),$$

(31)

where $a_k\sim p_k(·)$ and h is an arbitrary continuous bounded function.

Simulated annealing [36] is a type of the Metropolis–Hastings MCMC algorithm in which the temperature ${\sigma }_e$ decreases during the simulation. The algorithm is summarized below as Algorithm 2. In this paper we consider two temperatures.

Algorithm 2: Two-temperature simulated annealing (SA)

Set large ${\sigma }_e$ as a high temperature. Generate $a^{\prime }\sim q(·|a_k)=\mathcal{N}(a_k,{\epsilon }^2)$ with $\epsilon >0$ for given $a_k$. Calculate $\alpha (a^{\prime },a_k)=min\left\{1,\pi (a^{\prime }|\mathbf{U})/\pi (a_k|\mathbf{U})\right\}$. Update $a_k$ as $a_{k+1}=a^{\prime }$ with probability $\alpha (a^{\prime },a_k)$ but otherwise set $a_{k+1}=a_k$. Continue while $k\le k_b$. Then decrease ${\sigma }_e$ to a smaller value as a low temperature, and continue to update $a_k$.

Now we propose an MCMC-iterative hybrid method by combining Algorithms 1 with the Metropolis–Hastings Markov chain Monte Carlo. We first run the Markov chain Monte Carlo. Although the Metropolis–Hastings MCMC eventually gives the correct solution after about 10,000 steps when the chain reaches the steady state, we stop the Monte Carlo simulation way before the burn-in time. At the $k_b$th Monte Carlo step after initial rapid changes ceases, we record the obtained reconstructed value $a_{k_b}$ and switch to Algorithm 1 with the initial value the recorded $a_{k_b}$. It is important that $k_b$ is chosen after the initial rapid changes although $k_b$ can be a significant way from the burn-in time. Otherwise, the final reconstructed results are quite robust against the choice of $k_b$. We refer to the following hybrid algorithm as Algorithm 3.

Algorithm 3: Hybrid

Choose an initial guess $a_0$, which is not necessarily close to the global minimum. Generate $a^{\prime }\sim q(·|a_k)=\mathcal{N}(a_k,{\epsilon }^2)$ with $\epsilon >0$ for given $a_k$. Calculate $\alpha (a^{\prime },a_k)=min\left\{1,\pi (a^{\prime }|U)/\pi (a_k|U)\right\}$. Update $a_k$ as $a_{k+1}=a^{\prime }$ with probability $\alpha (a^{\prime },a_k)$ but otherwise set $a_{k+1}=a_k$. Obtain a reconstructed $a_{k_b}$. Switch to Algorithm 1 with the initial guess $a_{k_b}$.

We close this subsection by the Gelman-Rubin convergence diagnostic, which uses the intra-chain variance and inter-chain variance [37,38]. We run $M_{\mathrm{MC}}$ different chains with different initial values. Let $a_k^m$ denote the kth value in the mth chain ($k=1,\dots ,k_b$, $m=1,\dots ,M_{\mathrm{MC}}$). We discard the first $k_a-1$ steps before the initial rapid changes cease. Then we compute the following intra-chain average and variance.

$${\overline{a}}_m=\frac{1}{k_b-k_a+1}\sum _{k=k_a}^{k_b}{a}_k^m,{\sigma }_m^2=\frac{1}{k_b-k_a}\sum _{k=k_a}^{k_b}{(a_k^m-{\overline{a}}_m)}^2.$$

(32)

Next we compute the inter-chain mean and variance.

$$\overline{a}=\frac{1}{M_{\mathrm{MC}}}\sum _{m=1}^{M_{\mathrm{MC}}}{\overline{a}}_m,B=\frac{k_b-k_a+1}{M_{\mathrm{MC}}-1}\sum _{m=1}^{M_{\mathrm{MC}}}{({\overline{a}}_m-\overline{a})}^2.$$

(33)

Now we introduce

$$\widehat{V}=\frac{k_b-k_a}{k_b-k_a+1}W+\frac{1}{k_b-k_a+1}B,$$

(34)

where $W=\frac{1}{M_{\mathrm{MC}}}{\sum }_{m=1}^{M_{\mathrm{MC}}}{\sigma }_m^2$. This is an unbiased estimator of the true variance. But W is also an unbiased estimate of the true variance if the chains converge. Therefore, we have $\sqrt{\widehat{V}/W}\approx 1$ if converged. Below, we will see that we can choose $k_b$ for which $\sqrt{\widehat{V}/W}$ is not close to 1. In this paper, we set $M_{\mathrm{MC}}=10$.

2.4. TRS Measurements of a Polyurethane-Based Phantom

In this section, we consider time-resolved measurements for a phantom. The solid phantom (biomimic optical phantom) is made of polyurethane to simulate biological tissues (INO, Quebec, QC, Canada). The absorption coefficient and reduced scattering coefficient are ${\mu }_a=0.0209{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.853{\mathrm{mm}}^{-1}$ at wavelength 800 nm. The refractive index of the phantom is $n=1.51$. The measurements were performed by the TRS (time-resolved spectroscopy) instrument (TRS-80, Hamamatsu Photonics K.K., Hamamatsu, Japan). Two optical fibers from TRS-80 are attached on the top of the phantom with the separation $13\mathrm{mm}$ ($M_s=M_d=1$). The phantom is illuminated by near-infrared light of the wavelength $760\mathrm{nm}$ through one optical fiber and the reflected light is detected by the other optical fiber. The time interval of the measured data is $10\mathrm{ps}$ and we used the data from $t_1=2.00\mathrm{ns}$ to $t_{M_t}=8.00\mathrm{ns}$ ($M_t=601$). Measured counts, i.e., the number of photons, are shown in Figure 1. The upper panel of Figure 1 shows the instrument response function (IRF), which is given as a property of the experimental device. The measured reflected light is shown in the lower panel of Figure 1. The function $q(t^k)$ is the instrument response function divided by the maximum value of the measured reflected light. The peak of q is at $2.56\mathrm{ns}$. ${\left\{U\right\}}_{11k}=y_{11k}$ is the measured reflected light normalized by its maximum value. In the lower panel of Figure 1, $t_1$ and $t_{M_t}$ are marked by vertical dashed lines. The size of the phantom is large enough that we can assume the three-dimensional half space. Then the energy density of the detected light is computed by Equation (6). The parameter ℓ is obtained from Equation (3).

2.5. Numerical Phantom

To examine when the iterative scheme fails by being trapped by a local minimum and how the Markov chain Monte Carlo is capable of escaping from such local minima, a toy model is devised which is simple enough to explicitly understand the structure of the cost function but is complicated enough that the cost function has one local minimum and one global minimum.

We consider diffuse optical tomography in the two-dimensional space. In our toy model we suppose that the diffusion coefficient D is constant everywhere in the medium but there is absorption inhomogeneity and the absorption coefficient $\delta {\mu }_a$ is unknown. Moreover, we assume that $\delta {\mu }_a(\mathit{\rho })$ is given by

$$\delta {\mu }_a(\mathit{\rho })=\eta f_a(x)\delta (y-y_0),$$

(35)

where $\eta ,y_0$ are given positive constants. Here, $f_a(x)$ is given by

$$f_a(x)=\left[a^3+3\left(1+\frac{tanh{x}^2}{10}\right)a^2\right]\left(1-tanh{x}^2\right),$$

(36)

where $a>1.1$ is a constant. Thus, a is the parameter to be reconstructed in our toy inverse problem. As is shown in Figure 2, the function $f_a(x)$ has one peak at $x=0$ and the maximum value is $f_a(0)=a^2(a+3)$, $f_a$ monotonically decays for $|x|>0$, and $f_a\to 0$ as $|x|\to \infty$.

Now we describe how the forward data is computed. The unit of length and unit of time are taken to be $\mathrm{mm}$ and $\mathrm{ps}$, respectively. On the x-axis, we place two sources ($M_s=2$) at ${\mathit{\rho }}_s^1=(-20,0)$, ${\mathit{\rho }}_s^2=(20,0)$ and three detectors ($M_d=3$) at ${\mathit{\rho }}_d^1=(-40,0)$, ${\mathit{\rho }}_d^2=(0,0)$, ${\mathit{\rho }}_d^3=(40,0)$. We set ${\mu }_s^{\prime }=1$, ${\mu }_a=0.02$, $n=1.37$. Suppose that there is absorption inhomogeneity at depth 5. For $\delta {\mu }_a$, we put $\eta =300/c$, $y_0=5$, and

$$a=1.5.$$

(37)

To distinguish, hereafter the true value of a is denoted by $\overline{a}$. We assume 3% Gaussian noise and give the measured data as

$$y_{ijk}=u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i;\overline{a})(1+e_{ijk}),$$

(38)

where $e_{ijk}\sim \mathcal{N}(0,{0.03}^2)$. In this numerical experiment, $M_t=500$ and $t^{k+1}-t^k=t^1=5$ ($k=1,\dots ,M_t-1$).

By substituting the form (35) of $\delta b$ in (21), we can express the energy density as

$$\begin{array}{cc}\hfill u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i;a)& =u({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i)\hfill \\ & =u_0({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i)exp[-\frac{\eta e^{-{\mu }_{a0}c{t}^k}}{u_0({\mathit{\rho }}_d^j,t^k;{\mathit{\rho }}_s^i)}{\int }_0^{t^k}g(0,t^k;y_0,s)\hfill \\ & \times g(y_0,s;0,0)\left({\int }_{-\infty }^{\infty }{f}_a(x^{\prime })e^{-\frac{{(x_d^j-x^{\prime })}^2}{4Dc(t^k-s)}}{e}^{-\frac{{(x^{\prime }-x_s^i)}^2}{4Dcs}}d{x}^{\prime }\right)ds],\hfill \end{array}$$

(39)

where

$$u_0=\frac{e^{-{\mu }_{a0}c{t}^k}}{2\pi D{t}^k}{e}^{-\frac{{(x_d^j-x_s^i)}^2}{4Dc{t}^k}}\left[1-\frac{\sqrt{\pi Dc{t}^k}}{\ell }{e}^{{\left(\frac{\sqrt{Dc{t}^k}}{\ell }\right)}^2}\mathrm{erfc}\left(\frac{\sqrt{Dc{t}^k}}{\ell }\right)\right].$$

(40)

3. Results

3.1. Determination of Optical Properties

Here we consider reconstructions from measured data in the phantom experiment. Let us first consider initial guesses below.

$$\mathrm{Case}1:{\mu }_a=0.01{\mathrm{mm}}^{-1},{\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}.$$

(41)

In this Case 1, the following ${\mu }_a,{\mu }_s^{\prime }$ are obtained both by Algorithm 1 (LM) and Algorithm 3 (hybrid).

$${\mu }_a=0.016{\mathrm{mm}}^{-1},{\mu }_s^{\prime }=0.63{\mathrm{mm}}^{-1}.$$

(42)

The results for ${\mu }_a$ and ${\mu }_s^{\prime }$ are shown in Figure 3 and Figure 4, respectively. In Figure 3, reconstructed values of ${\mu }_a$ are plotted against the iteration number k. In the top panel of Figure 3, we see that Algorithm 1 (LM) quickly converges. In Algorithms 2 and 3, we set $k_b=99$. As is seen in the middle panel of Figure 3, the convergence of Algorithm 2 (SA) is slow. The temperature is decreased from ${\sigma }_e={10}^{-6}$ to ${\sigma }_e={10}^{-7}$ at the $k_b$th Monte Carlo step. Similarly, $\epsilon$ is changed from $0.1$ to $0.001$. Algorithm 2 returns the correct values after many Monte Carlo steps; we have $\sqrt{\widehat{V}/W}=1.02,1.06$ for ${\mu }_a,{\mu }_s^{\prime }$, respectively, when $k_a=9000$ and $k_b=10000$. Finally, in the bottom panel of Figure 3, we see that the iteration of Algorithm 3 (hybrid) immediately converges after we switch from the Metropolis–Hastings MCMC (${\sigma }_e={10}^{-6}$, $\epsilon =0.1$) to Algorithm 1 (LM). No convergence of the Monte Carlo chain is required, and we found $\sqrt{\widehat{V}/W}=8.6$ for ${\mu }_a$ and $=11.0$ for ${\mu }_s^{\prime }$ ($k_a=49$, $k_b=99$). A similar behavior is observed for the reconstruction of ${\mu }_s^{\prime }$ in Figure 4. In the top panel of Figure 4, we find that Algorithm 1 (LM) works best and converges after a few iterations, whereas Algorithm 2 (SA) in the middle panel of Figure 4 has slow convergence and reconstructed values around at $k_b=99$ are still away from ${\mu }_s^{\prime }$ in (42). By switching from MH-MCMC to Algorithm 1 (LM) using Algorithm 3 (hybrid), convergence is easily obtained as shown in the bottom panel of Figure 4.

Now we start the simulation by setting the following initial values.

$$\mathrm{Case}2:{\mu }_a=0.5{\mathrm{mm}}^{-1},{\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}.$$

(43)

In this Case 2, Algorithm 1 (LM) fails and returns ${\mu }_a=0.068{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.75{\mathrm{mm}}^{-1}$ whereas Algorithm 3 (hybrid) still gives the correct values. The reconstructed values of ${\mu }_a$ and ${\mu }_s^{\prime }$ at each iteration are shown in Figure 5 and Figure 6, respectively. In Figure 5, the vertical axes of the left three panels are from 0 to $0.6$ whereas the vertical axes of the right three panels are from 0 to $0.08$. In Figure 6, the vertical axes of the left three panels are from 0 to 2 while the vertical axes of the right three panels are from 0 to $1.3$. In the top panels of Figure 5, we see that the reconstruction of ${\mu }_a$ is unsuccessful by Algorithm 1 (LM). Algorithm 2 (SA) approaches ${\mu }_a$ in Equation (42) but still deviates from that value in the middle panels of Figure 5. As shown in the bottom panels of Figure 5, Algorithm 3 (hybrid) converges in a few iterations after switching to Algorithm 1 (LM). In the top panels of Figure 6, Algorithm 1 (LM) fails to reconstruct ${\mu }_s^{\prime }$. In the middle panels of Figure 6, reconstructed values by Algorithm 2 (SA) come close to ${\mu }_s^{\prime }$ in Equation (42) but suffer from slow convergence. In the bottom panels of Figure 6, we see that Algorithm 3 (hybrid) arrives at ${\mu }_s^{\prime }$ in Equation (42) without any problem.

The initial guesses for Case 1 are reasonably close to the values found in Equation (42). However, we started with initial guesses which are quite different from the above-mentioned values in Case 2. It is not surprising that Algorithm 1 (LM) does not work for Case 2, but Algorithm 3 (hybrid) can arrive at the correct values. The numerical calculations were performed on Matlab (i7-8700 CPU, 16 GB memory). In the hybrid method, the Metropolis–Hastings MCMC does not reach the steady state at about 100 steps but can be switched to the Levenberg–Marquardt method. For Figure 5 and Figure 6, the computation time was $5\mathbf{s}$. The simulated annealing method returns the correct value after about 10,000 steps, but it takes $8\min$. The computation for Algorithm 1 (LM) stopped in $0.5\mathbf{s}$.

Below we summarize the reconstructed values on

Table 1. Reconstructed values of ${\mu }_a,{\mu }_s^{\prime }$ are shown for Case 1 and Case 2. The units of ${\mu }_a,{\mu }_s^{\prime }$ are both ${\mathrm{mm}}^{-1}$. For Algorithm 2, values at 120th Monte Carlo step are shown.

	Case 1 (${\mathit{\mu }}_{\mathit{a}}$, ${\mathit{\mu }}_{\mathit{s}}^{\prime }$)	Case 2 (${\mathit{\mu }}_{\mathit{a}}$, ${\mathit{\mu }}_{\mathit{s}}^{\prime }$)
initial values	($0.01$, $1.0$)	($0.5$, $1.0$)
Algorithm 1 (LM)	($0.016$, $0.63$)	($0.068$, $1.75$)
Algorithm 2 (SA)	($0.032$, $1.02$)	($0.028$, $0.92$)
Algorithm 3 (hybrid)	($0.016$, $0.63$)	($0.016$, $0.63$)

. Although Algorithm 2 (SA) and Algorithm 3 (hybrid) return the same results after a long time, the hybrid scheme is about ten times faster, and moreover there is no need for choosing the low temperature for the latter algorithm. For Algorithm 3, it is not necessary to wait until the burn-in time, but it is enough if the initial rapid change ceases. In Section 3.2, we illustrate that Algorithm 3 works once the Monte Carlo simulation escapes from a local minimum and the algorithm does not require that the calculation reaches the steady state.

3.2. Determination of Absorption Inhomogeneity

We perform a numerical experiment of diffuse optical tomography using the toy model. Let us consider when the iterative scheme fails. We see $F^{\prime }(0)={\left.\frac{\partial }{\partial a}u({\mathit{\rho }}^j,t^k;{\mathit{\rho }}^i;a)\right|}_{a=0}=0$. For sufficiently small $ϵ>0$, we have $F^{\prime }(ϵ)>0$, $F^{\prime }(-ϵ)<0$, $|F^{\prime }(\pm ϵ)|\simeq 0$, and $U-F(\pm ϵ)<0$. Therefore, if we start the iteration from $a_0=ϵ$, we obtain

$$|a_1|<ϵ,|a_2|<|a_1|,\dots .$$

(44)

Thus, the sequence $a_k$ approaches 0 and can never arrive at $\overline{a}$ ($=1.5$).

Let us consider how the difference $U-F(a)$ depends on a. We introduce

$$h(t)=\frac{1}{t}exp\left(-\frac{y_0^2}{4Dct}\right)\left[1-\frac{\sqrt{\pi Dct}}{\ell }{e}^{{\left(\frac{y_0}{2\sqrt{Dct}}+\frac{\sqrt{Dct}}{\ell }\right)}^2}\mathrm{erfc}\left(\frac{y_0}{2\sqrt{Dct}}+\frac{\sqrt{Dct}}{\ell }\right)\right].$$

(45)

The following form is obtained using Equation (39). By neglecting noise, we have

$$\begin{array}{c}u({\mathit{\rho }}^j,t^k;{\mathit{\rho }}^i;\overline{a})-u({\mathit{\rho }}^j,t^k;{\mathit{\rho }}^i;a)\hfill \\ =\frac{\eta e^{-{\mu }_{a0}c{t}^k}}{{(2\pi D)}^2}{\int }_0^{t^k}h(t^k-s)h(s)\left[{\int }_{-\infty }^{\infty }{d}_{\overline{a}}(a,x^{\prime })\left(1-tanh{x^{\prime }}^2\right)e^{-\frac{{(x_d^j-x^{\prime })}^2}{4Dc(t^k-s)}}{e}^{-\frac{{(x^{\prime }-x_s^i)}^2}{4Dcs}}d{x}^{\prime }\right]ds,\hfill \end{array}$$

(46)

where $d_{\overline{a}}(a,x^{\prime })=\xi (\overline{a},x^{\prime })-\xi (a,x^{\prime })$ with

$$\xi (a,x^{\prime })=a^2\left[a+3\left(1+\frac{tanh{x^{\prime }}^2}{10}\right)\right].$$

(47)

For a given $x^{\prime }$, the function $|d_{\overline{a}}(a;x^{\prime }){|}^2$ has one global minimum at $a=\overline{a}$, one local minimum at $a=-2\left(1+\frac{tanh{x^{\prime }}^2}{10}\right)$, and one local maximum at $a=0$. Therefore, the above expression implies that Algorithm 1 (LM) fails when the initial value $a_0$ is negative and the correct value $\overline{a}=1.5$ is reconstructed only for a positive initial guess. Indeed, the computation ends up with $-2.05$ if we start the iteration from $a_0=-0.1$ (see below) or $-0.01$, and the value $1.68$ is obtained when $a_0=0.01$. Figure 7 shows $|d_{\overline{a}}{(a,x)|}^2$ for $\overline{a}=1.5$ and $tanh(x^2)=0.5$.

In Figure 8, we plot reconstructed values of a against iteration numbers. The initial value is set to $a_0=-0.1$. The reconstruction by Algorithm 1 (LM) is shown in the top panel of Figure 8. As we can predict from Figure 7, Algorithm 1 (LM) converges to the local minimum and can never arrive at the global minimum. Monte Carlo simulation can jump over the local maximum and approach the global minimum, but keeps fluctuating as shown in the middle panel of Figure 8. We initially set $({\sigma }_e,\epsilon )=({10}^{-6},0.5)$ for Algorithms 2 and 3. In Algorithm 2, we set $({\sigma }_e,\epsilon )=({10}^{-7},0.005)$ after the temperature decreases at the $k_b$th Monte Carlo step. In the bottom panel of Figure 8, Algorithm 3 (hybrid) successfully arrives at the global minimum. Using Matlab, the computation time was $17\min$ while we needed $3\mathrm{hr}$ for Algorithm 2 (SA) (1000 steps).

After the initial time with large jumps, it is found that we can set $k_b=99$ in our simulation of the MCMC-iterative hybrid method. The correct value of $\overline{a}$ is reconstructed by Algorithm 3 (hybrid) even when the simulation experiences a local minimum. Starting from $a_{k_b}=1.8257$, the calculation stops at $a^{*}=1.6841$ ($a^{*}=a_{k_b+3}$). We note that the reconstructed $a^{*}$ is not exactly $1.5$ due to noise.

When the initial guess $a_0=-0.1$ lies in the valley of the local minimum (see Figure 7), the reconstructed a by Algorithm 3 (hybrid) moves to the valley of the global minimum with the help of Monte Carlo simulation as shown in Figure 8, while the reconstructed a by Algorithm 1 (LM) falls to the local minimum as the nature of Newton-type iterative methods. For Algorithm 3 (hybrid), we obtained $a^{*}=1.6841$. There is a possibility that negative reconstructed values are obtained by Algorithm 2 (SA) and Algorithm 3 (hybrid) for the first a few iterations if different pseudo-random numbers are used. These negative reconstructed values, however, will turn positive and the behavior of reconstructed a by Algorithm 2 (SA) and Algorithm 3 (hybrid) is always more or less similar to the middle and bottom panels of Figure 8. For the initial guesses of $a_0=-0.01$ and $a_0=0.01$, Algorithm 3 (hybrid) works without any problem and Algorithm 2 (SA) also returns a reasonable result after a sufficiently large number of iterations (results not shown).

4. Discussion

In this paper, we have proposed a hybrid numerical scheme which uses Markov chain Monte Carlo in the first step and then uses an iterative method in the second step. We switch from MH-MCMC to LM by observing proposed parameter values. For the typical jump length of parameters in MH-MCMC, $\epsilon =0.1$ was used in Section 3.1 and $\epsilon =0.5$ was chosen in Section 3.2. Although these values were set so that the MH-MCMC calculation was efficiently performed, other choices of $\epsilon$ are also possible. The proposed scheme is quite general and can be applied to different inverse problems in NIRS which are solved by iterative methods even when the forward problem must be solved fully numerically with finite difference method or finite element method. It is an interesting future study how the hybrid scheme can be extended to diffuse optical tomography, which has many unknowns.

More sophisticated algorithms than the Metropolis–Hastings Markov chain Monte Carlo used in this paper have been proposed to overcome slow convergence. The delayed rejection scheme reduces the net rejection rate [39]. In the adaptive Metropolis algorithm, parameters in the proposal distribution are adjusted during Monte Carlo steps [40]. The DRAM algorithm, which combines the above-mentioned two schemes, was also proposed [41]. Two-level MCMC algorithms [42,43] and a multi-level MCMC [23] have been investigated to improve the MCMC algorithm. Such Monte Carlo schemes might improve the first step of our hybrid method by finding the valley of the global minimum more easily.

Related to the Metropolis–Hastings Markov chain Monte Carlo algorithm, quantum annealing [44] has been developed in addition to simulated annealing. Aiming at escaping from local minima, brute-force search and genetic algorithm [45] are also well-known optimization algorithms. The introduction of these methods in NIRS may be found useful in the future.

For the clinical use of NIRS, it has been recognized from early days that finding absolute values of the absorption and scattering coefficients is important [2]. In Ref. [15], it is emphasized that the obtained absolute values highly depend on the starting values of the initial estimate for the study of the infant brain with their measurement system. By performing Markov chain Monte Carlo before using standard iterative schemes, we may automatically acquire good initial values. Such clinical application is a natural next step of our research.