# A Homotopy Bayesian Approach for Inverse Problems

Authors: Xiao-Mei Yang; Zhi-Liang Deng
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In solving Bayesian inverse problems, it is often desirable to use a common density parameterization to characterize the prior and posterior. Typically, we seek an optimal approximation of the true posterior from the same distribution family as the prior. As one of the most important classes of distributions in statistics, the exponential family is considered as the parameterization. The optimal parameter values for representing the approximated posterior are achieved by minimizing the deviation between the parameterized density and a homotopy that deforms the prior density into the posterior density. Instead of the usual moment parameters, we introduce a new parameter system in the process, by which we reduce the dimension of the parameter and therefore the computational cost. The parameters are ruled by a system of explicit first-order differential equations. Solving this system over finite “time” interval yields the desired optimal density parameters. This method is proven to be effective using some numerical examples.


---

## Body

## 1. Introduction

Inverse problems occur widely in mathematical and engineering fields. The related theories and algorithms have been developed by many authors [1]. To analyze the uncertainty of the unknowns, the Bayesian approach is widely applied to find the solutions to inverse problems [2–4]. It provides quantitative estimations of the probability distribution of the solution which can be used to quantify the confidence. The prior distribution is transited to the posterior distribution of the unknowns. This transition is important and causes numerical challenge in finding the posterior distribution. The Markov chain Monte Carlo (MCMC) methods, e.g., the Metropolis–Hastings (MH) sampling algorithm, are essential tools for analyzing it. These algorithms allow us to evaluate arbitrary probability distributions; however, they are inherently sequential in nature due to the Markov property, which causes the autocorrelation of the samples.It is inevitable to increase the uncertainty of the estimation of posterior quantities of interest, such as the means and the variance.We seek an optimal approximation of the posterior distribution rather than draw the posterior samples by the MCMC method. We assume that the approximated posterior and the prior are of the same distribution type and share the common parameterization density. The transition from the prior to the approximated posterior involves the change of the distribution parameters. For such problems, the homotopy method is a promising approach, which characterizes solution spaces by smoothly tracking solution from one formulation (typically an “easy” problem) to another (typically a “hard” problem) [5]. Hanebeck et al. [6] introduced a general framework to perform the measurement update of filter problems with a homotopy. Furthermore, Hagmar et al. [7] developed this framework in filter problems and derived a so-called homotopy differential equation (HDE) to rule the progressive processing of the distribution parameters. In [6, 7], the Gaussian and mixed Gaussian are adopted in the approximated process. To the authors’ knowledge, there is no related work using the homotopy idea to solve the Bayesian inverse problems. In this paper, we extend this homotopy idea and use it to treat the Bayesian inverse problems. As an extension of [6, 7], the approximated posterior density family is chosen as the exponential and the mixed exponential family, and the usual moment parameters are replaced by natural parameters.Specially for Gaussian, the dimension of the parameters can be reduced.Compared with the MCMC method, the proposed method is possible to avoid to draw plenty of autocorrelated posterior samples.We apply the proposed method to solve two inverse problems in partial differential equations: inverse heat conduction problem and inverse scattering problem. The finite element method (FEM) and boundary integral equation method are used to solve the forward problem [8, 9]. For some other numerical techniques for differential equations, we refer to [10–18]. The Bayesian perspective to inverse problems is developed rapidly in the past two decades. There are two main points about the algorithms in the Bayesian inversion: MCMC-based sampling methods [19–27] and variational methods [28–30]. For the MCMC-based sampling algorithms, some fast numerical techniques are studied by many authors. In [4–33], the polynomial chaos expansion is used to the surrogate of the forward model or the likelihood function. The Gaussian process regression method is also used to the acceleration of the Bayesian inversion [34, 35]. Despite the large body of works on theories and applications of Bayesian approaches in inverse problems, the homotopy idea is still quite an unexplored terrain. Our work develops it and is expected to extend to more fields of applications.The rest of this paper is organized as follows: in Section2, we give the basic framework of the Bayesian inversion using homotopy. In Section 3, we introduce the exponential and the mixed exponential families. In Section 4, the approximation of the HDE is derived. In Section 5, some examples are given to verify the effectiveness of the proposed algorithm. Finally, we give the conclusion.


## 2. Bayesian Inversion Using Homotopy

Inverse problems are concerned with the subject of converting observational data into information which are not observed directly in systems. In mathematics, an inverse problem takes the abstract form:(1)y˜=ςk+Ξ,in which Ξ is an additive noise, the unknown κ∈U is to be determined, given the data y˜∈Y, where U and Y are Banach spaces. We apply a probabilistic viewpoint, the Bayesian approach, to give the solution to (1), in which all quantities including the unknown κ, the noise Ξ, and the observations y˜ are regarded as random variables. In the Bayesian framework, the information about the unknown is updated by blending prior beliefs with observed data. Typically, the prior and posterior are coded in the corresponding probability measures, which are linked by the Bayes’ formula:(2)dμy˜dμ0κ∝Lκ;y˜≔exp−Φκ;y˜,where Lκ;y˜ is the likelihood function and Φκ;y˜ is the negative log likelihood.For simplicity, we letU and Y be finite dimensional spaces. The posterior density corresponding to μy˜κ is given as(3)py˜κ=exp−Φκ;y˜qκZ,where qκ represents the prior density and Z is the normalization constant. The main target is the numerator, which is expressed in the following equation:(4)exp−Φκ;y˜qκ≔pκ.We try to find an optimal approximation ofpκ using homotopy.The homotopy Bayesian approach is the one to perform progressive processing for the posterior distributionpκ. In this method, instead of directly approximating the true density pκ, it starts with a tractable density and continuously approaches the true density via intermediate densities. Choose the homotopy as follows:(5)pκ,t=Ltκ;y˜qκ=exp−tΦκ;y˜qκ,where t∈0,1 is a parameter. It is obvious that pκ,1=pκ and pκ,0=qκ. For each t∈0,1, an optimal approximation of pκ,t is to be sought from some distribution family, which has the density gκ;η of parameters η. We minimize a deviation function Gη,t between pκ,t and gκ;η and obtain the optimal approximation of pκ,t. Obviously, when t=1, we get an optimal approximation of the posterior pκ.We usually use Kullback–Leibler (KL) divergence and the squared Hellinger distance to measure the deviation between two probability distributions over the same variableκ. In this paper, we assume that the random variables are continuous. For random variables X and Y, we denote their probability densities by xκ and yκ, respectively. The KL divergence and squared Hellinger distance between X and Y are given, respectively.(6)DKLX‖Y=∫xκlogxκyκdκ,H2X,Y=12∫xκ−yκ2dκ.Taking the deviationG with x=pκ,t and y=gκ;η in the sense of (5), i.e., G=DKL or G=H2, we seek to find the solution(7)ηt=argminηGη,t.The minimization necessary condition allows us to get(8)Gηηηt,tη′t+Gηtηt,t=0.Fort=0, the density gκ;η0 is an optimal approximation of qκ. If the approximated distributions gκ;η are chosen from the same family as the prior distribution qκ, the initial condition is set according to gκ;η0=qκ. Then, let qκ be identified as gκ;η0. Thereby, we set the initial condition η0=η0. Solving for η′t in (8), we arrive at the initial value problem of the HDE, which is expressed as follows:(9)η′t=−Gηη−1ηt,tGηtηt,t,η0=η0.


## 3. Exponential Family and Mixed Exponential Family

The exponential family (EF) and mixed exponential family serve as the approximation of the posterior distribution in this paper. For a numeric random variableK, the parametric EF probability density can be written as follows:(10)qκ;θ=hκexpTκ,θ−Aθ,where θ is called the natural (canonical) parameter, Tκ is the sufficient statistic, and Aθ is the log normalizer given by(11)Aθ=log∫hκexpTκ,θdκ.It is easy to know that(12)dlogqκ;θdθ=Tκ−A′θ.We can computeA′θ by A′θ=ETκ.For ad-dimensional Gaussian random variable K, K∼Nϰ,Σ, the probability density is given as follows:(13)qκ|ϰ,Σ=12πd/2exp−12trΣ−1κκ⊤+ϰ⊤Σ−1κ−12ϰ⊤Σ−1ϰ−12logΣ.In this case, we can reduce the dimension of parameters. Denote the precision matrix byP=Σ−1. By the symmetric positive definiteness of P, we have the Cholesky factorization P=ℛ⊤ℛ with ℛ being a lower triangular matrix. Then, by introducing the new parameter(14)θ≔ϰvechℛ≕θ1θ2,we have(15)∂logqκ|θ∂θ1=−ℛ⊤ℛκ−ϰ,∂logqκ|θ∂θ2=vechdiagℛ−1−ℛκ−ϰκ−ϰ⊤.Remark 1.
It is obvious that the dimension of the new parameter system decreases dramatically compared with the moment parameter or the natural parameter system. Parameterθ1 is the mean, and θ2 has a lower dimension than the covariance matrix.
Subsequently, for the MEF, we letgκ;η=∑i=1Mwiqκ;θi, i.e.,(16)gκ;η=∑i=1MwihκexpΦκ,θi−Aiθi,where ∑i=1Mwi=1 and wi≥0 for i=1:M. To remove the weight constraint, let wi=π/2+arctanλi/∑j=1Mπ/2+arctanλj and λM=0. Then, it follows by simple calculations that(17)∂loggκ;η∂θi=wiqκ;θigκ;η∂logqκ;θi∂θi,i=1:M,∂loggκ;η∂λi=1/1+λi2qκ;θi/gκ;η−1∑j=1Mπ/2+arctanλj,i=1:M−1.Remark 2.
The parameter dimension number of the mixed Gaussian exponential family isM−1+Md+dd+1/2.


## 4. Approximation of the HDE

In HDE (9), the homotopy density is involved in the integration term. The numerical evaluations are a challenge. Therefore, it is necessary to give an approximation version of (9) in a real numerical simulation process to avoid the direct computations for the homotopy densities. Next, we derive an approximated HDE in inverse problems.For a partition0=t0<t1<⋯<tn=1, we have from the Taylor expansion of ηt(18)ηi+1=ηi+ηi′Δti,where Δti=ti+1−ti, ηi, and ηi′ are the approximations of ηti and η′ti, respectively. We then can take(19)ηi′=−G˜ηη−1ηi,tiG˜ηtηi,ti,with G˜ηηηi,ti and G˜ηtηi,ti being some approximations of Gηηηti,ti and Gηtηti,ti in (9), respectively. Next, let us give these approximations using the deviations of measures. For the KL divergence, we have(20)Gη,t=∫pκ,tlogpκ,tgκ;ηdκ.When we implement the iteration (18) with (19), the density pκ,t at ti has been approximated by gκ;ηi. Therefore, we take the approximations of Gηηηti,ti and Gηtηti,ti as follows:(21)G˜ηηηi,ti=∫∂logg∂η∂logg∂η⊤g|ηidκ≔Iηi,G˜ηtηi,ti=∫Φκ;y˜∂gκ;η∂η|ηidκ,where the term Iη is the Fisher information matrix when gκ;η is a probability density. For the squared Hellinger metric, we have the corresponding expressions as follows:(22)G˜ηηηi,ti=14Iηi,G˜ηtηi,ti=14∫Φκ;y˜∂gκ;η∂η|ηidκ.For both cases, the approximation (19) can be given by(23)ηi′=−Iηi−1∫Φκ;y˜∂gκ;η∂η|ηidκ.Thereby, we have the approximated HDE(24)dη˜dt=−I−1η˜∫Φκ;y˜∂g∂η˜dκ.With these computations, we can implement the homotopy Bayesian approach for inverse problems.


## 5. Numerical Examples

We list the homotopy Bayes process in Algorithm 1.Algorithm 1: Homotopy Bayes algorithm.
(1)
Initial condition: Take the initial parameter from the prior distribution.(2)
Whilet=ti, compute the integration in (23) using the Monte Carlo algorithm and get ηi′. And set ηi+1=ηi+ηi′Δti and set ti+1=ti+Δti.(3)
Stop untilt=1 and use ηn as the approximated posterior parameter.


### 5.1. Inverse Heat Conduction Problem

We consider the reconstruction of the thermal conductivitiesκ using the measurements of solution u governed by the system, which is as follows:(25)∇⋅κ∇u=0in0,1×0,0.6,∂u∂x10,x2=∂u∂x11,x2=0,ux1,0.6=T=200,−κ0∂u∂x2x1,0=2000.This example is originally presented in [4]. We solve the problem with the homotopy Bayesian approach and investigate the performance of the algorithm. As displayed in Figure 1(a), the background thermal conductivity is denoted as κ0 that is known, while the conductivities of the material inclusions are termed as κ1 and κ2, respectively.Figure 1
Heat conduction setups. (a) 2D IHCP. (b) 6D IHCP.
(a)(b)We consider the inverse heat conduction problem (IHCP) that is posed when the thermal conductivitiesκ=κ1,κ2⊤ are unknown, and their inference is intended.We suppose that a number ofN measurements y˜=ux of the temperature field at the measurement locations x1,x2,…,xN⊤ are available. The forward model G:κ⟶y˜ is established by the finite element discretization of (25) [9]. With the discretization model, the measured temperatures y˜ are generated by adding ansatz noise to the numerical solutions as follows:(26)y˜=y+Ξ=Gκ+Ξ,where Ξ∼N0,δ2I. The prior is set to a multivariate lognormal distribution qκ=ɛ∏i=12qκi with independent marginals qκi=LNκi|ϰ0,σ0 with ϰ0=30 and σ0=6. Parameters ϰ0 and σ0 describe the mean and standard deviation of the lognormal prior. They are related to the parameters of the associated normal distribution Nlogκi|λ0,ζ02 via ϰ0=expλ0+ζ02/2 and σ02=expζ02−1exp2λ0+ζ02. The unknown parameters are represented as κi=expλ0+ζ0ξi in terms of the standardized variables ξi∼N0,1.In this computation process, we take a uniform partition for homotopy parameter interval0,1. We set Δti=0.01 in (18). These integrations in (23) are computed by the Monte Carlo numerical integration. We take the true κ1=32 and κ2=28. For noise δ=0.25, we display the numerical results in Figure 2. The estimated posterior density for ξ is given in Figure 2(c). We transform it to κ by κ=expλ0+ζ0ξ and show the plot in Figure 2(a). Accordingly, we plot the posterior for ξ and κ using the FEM evaluation for the forward problem in Figures 2(b) and 2(d), respectively.Figure 2
The numerical results for the 2D IHCP: the reconstruction forκ=expλ0+ζ0ξ ((a) and (b)) and the posterior density ((c) and (d)).
(a)(b)(c)(d)We also consider the IHCP with six unknown conductivities. The setup is displayed in Figure1(b).The unknown conductivitiesκ=κ1,…,κ6⊤ are inferred with N=20 noisy measurements y˜=ux1,…,ux20⊤. Taking δ=0.05 and the true parameters κ=20,24,…,40, we display the reconstructed probability in Figure 3 using the proposed algorithm.Figure 3
The reconstructed posterior probability density for 6D IHCP.We see from these numerical results that the homotopy method achieves satisfactory reconstructions for 2D and 6D IHCP.


### 5.2. Inverse Acoustic Obstacle Scattering

We apply the proposed algorithm to an inverse acoustic scattering problem with a sound-soft obstacle. We consider the scattering by long cylindrical obstacles with cross-sectionsΩ⊂ℝ2 that are star-like with respect to the origin. In mathematics, we assume that Ω⊂ℝ2 is a bounded, simply connected domain with C2 boundary ∂Ω. Then, ∂Ω can be uniquely represented by a periodic function r:0,2π⟶ℝ+:(27)∂Ω≔rscoss,sins=expqscoss,sins,where s∈0,2π,qs=logrs, 0<rs<rmax.For a given incident plane waveuix≔expikx⋅d,x∈ℝ2,d∈S≔x^∈ℝ2|x^=1, where k>0 is the wavenumber, i=−1 and d≔cosς,sinς is the direction, the scattering problem is to find the scattered field us, or the total field u=ui+us, such that(28a)Δu+k2u=0inℝ2Ω¯,u=0on∂Ω,(28b)limr⟶∞r∂us∂r−ikus=0.It is well-known from the Sommerfeld radiation condition (28b) that the scattered field admits the following asymptotic expansion:(29)usx,d=expiπ/48kπexpikrru∞x^,d+O1r,as r≔|x|⟶∞ uniformly for all directions x^=x/|x|. The function u∞x^,d defined on the unit circle S⊂ℝ2 is called the far field pattern of us. We can write the forward problem by [36](30)u∞x^,d=GΩ,where G is the shape-to-measurement operator. Using the parameterization (27) and taking the noise in measurements into account, the inverse model is given by(31)y˜=Gqx^,d+Ξ,x^,d∈Γo×Γi⊂S×S,where Γo is the aperture of the observation and Γi is the aperture of the incident wave.The priorq is taken as the truncated Fourier series [36](32)qNs=κ02π+∑n=1Nκnnυcosnsπ+κ˜nnυsinnsπ,where κn and κ˜n are i.i.d. (independent and identically distributed) with κn,κ˜n∼N0,1 and υ is a positive constant.The ansatz data are used to the numerical reconstruction. The synthetic Gaussian noise is added to the true forward model, i.e.,(33)y˜=Gqx^,d+δGqx^,dξ1+ξ2i,ξ1,ξ2∼N0,1.In the numerical implementation, we fix wavenumberk=1, δ=0.01 and take υ=2.2,N=5 in (32). For two incident waves Γi=cosπ/2,sinπ/2∪cos3π/2,sin3π/2, we collect the full aperture data, Γo=coss,sins, s∈0,2π. Some shapes are chosen as the test examples. We draw 100 random samples from the approximated posterior distribution, which is displayed in Figure 4. It can be seen that the proposed algorithm provides satisfactory reconstructions for the obstacle scattering problem.(a)
Threelobes:rs=0.5+0.25e−sin3s−0.1sins;(b)
Pear:rs=5+sin3s/6;(c)
Acorn:rs=3/517/4+2cos3s;(d)
Roundedtriangle:rs=2+0.5coss;(e)
Roundrect:rs=cos4s+2/3sins4−1/4;(f)
Kite:x1=coss+0.65cos2s−0.65, x2=1.5sins.Figure 4
The numerical reconstructions for inverse scattering problem using full aperture data with two incident waves: (a)–(c) with incident anglesς=0 and ς=π/2 and (d)–(f) with ς=0 and ς=π. The blue curves are the posterior samples drawn from the approximated posterior distribution. (a) Three lobes. (b) Pear. (c) Acorn. (d) Rounded triangle. (e) Round rect. (f) Kite.
(a)(b)(c)(d)(e)(f)


## 5.1. Inverse Heat Conduction Problem

We consider the reconstruction of the thermal conductivitiesκ using the measurements of solution u governed by the system, which is as follows:(25)∇⋅κ∇u=0in0,1×0,0.6,∂u∂x10,x2=∂u∂x11,x2=0,ux1,0.6=T=200,−κ0∂u∂x2x1,0=2000.This example is originally presented in [4]. We solve the problem with the homotopy Bayesian approach and investigate the performance of the algorithm. As displayed in Figure 1(a), the background thermal conductivity is denoted as κ0 that is known, while the conductivities of the material inclusions are termed as κ1 and κ2, respectively.Figure 1
Heat conduction setups. (a) 2D IHCP. (b) 6D IHCP.
(a)(b)We consider the inverse heat conduction problem (IHCP) that is posed when the thermal conductivitiesκ=κ1,κ2⊤ are unknown, and their inference is intended.We suppose that a number ofN measurements y˜=ux of the temperature field at the measurement locations x1,x2,…,xN⊤ are available. The forward model G:κ⟶y˜ is established by the finite element discretization of (25) [9]. With the discretization model, the measured temperatures y˜ are generated by adding ansatz noise to the numerical solutions as follows:(26)y˜=y+Ξ=Gκ+Ξ,where Ξ∼N0,δ2I. The prior is set to a multivariate lognormal distribution qκ=ɛ∏i=12qκi with independent marginals qκi=LNκi|ϰ0,σ0 with ϰ0=30 and σ0=6. Parameters ϰ0 and σ0 describe the mean and standard deviation of the lognormal prior. They are related to the parameters of the associated normal distribution Nlogκi|λ0,ζ02 via ϰ0=expλ0+ζ02/2 and σ02=expζ02−1exp2λ0+ζ02. The unknown parameters are represented as κi=expλ0+ζ0ξi in terms of the standardized variables ξi∼N0,1.In this computation process, we take a uniform partition for homotopy parameter interval0,1. We set Δti=0.01 in (18). These integrations in (23) are computed by the Monte Carlo numerical integration. We take the true κ1=32 and κ2=28. For noise δ=0.25, we display the numerical results in Figure 2. The estimated posterior density for ξ is given in Figure 2(c). We transform it to κ by κ=expλ0+ζ0ξ and show the plot in Figure 2(a). Accordingly, we plot the posterior for ξ and κ using the FEM evaluation for the forward problem in Figures 2(b) and 2(d), respectively.Figure 2
The numerical results for the 2D IHCP: the reconstruction forκ=expλ0+ζ0ξ ((a) and (b)) and the posterior density ((c) and (d)).
(a)(b)(c)(d)We also consider the IHCP with six unknown conductivities. The setup is displayed in Figure1(b).The unknown conductivitiesκ=κ1,…,κ6⊤ are inferred with N=20 noisy measurements y˜=ux1,…,ux20⊤. Taking δ=0.05 and the true parameters κ=20,24,…,40, we display the reconstructed probability in Figure 3 using the proposed algorithm.Figure 3
The reconstructed posterior probability density for 6D IHCP.We see from these numerical results that the homotopy method achieves satisfactory reconstructions for 2D and 6D IHCP.


## 5.2. Inverse Acoustic Obstacle Scattering

We apply the proposed algorithm to an inverse acoustic scattering problem with a sound-soft obstacle. We consider the scattering by long cylindrical obstacles with cross-sectionsΩ⊂ℝ2 that are star-like with respect to the origin. In mathematics, we assume that Ω⊂ℝ2 is a bounded, simply connected domain with C2 boundary ∂Ω. Then, ∂Ω can be uniquely represented by a periodic function r:0,2π⟶ℝ+:(27)∂Ω≔rscoss,sins=expqscoss,sins,where s∈0,2π,qs=logrs, 0<rs<rmax.For a given incident plane waveuix≔expikx⋅d,x∈ℝ2,d∈S≔x^∈ℝ2|x^=1, where k>0 is the wavenumber, i=−1 and d≔cosς,sinς is the direction, the scattering problem is to find the scattered field us, or the total field u=ui+us, such that(28a)Δu+k2u=0inℝ2Ω¯,u=0on∂Ω,(28b)limr⟶∞r∂us∂r−ikus=0.It is well-known from the Sommerfeld radiation condition (28b) that the scattered field admits the following asymptotic expansion:(29)usx,d=expiπ/48kπexpikrru∞x^,d+O1r,as r≔|x|⟶∞ uniformly for all directions x^=x/|x|. The function u∞x^,d defined on the unit circle S⊂ℝ2 is called the far field pattern of us. We can write the forward problem by [36](30)u∞x^,d=GΩ,where G is the shape-to-measurement operator. Using the parameterization (27) and taking the noise in measurements into account, the inverse model is given by(31)y˜=Gqx^,d+Ξ,x^,d∈Γo×Γi⊂S×S,where Γo is the aperture of the observation and Γi is the aperture of the incident wave.The priorq is taken as the truncated Fourier series [36](32)qNs=κ02π+∑n=1Nκnnυcosnsπ+κ˜nnυsinnsπ,where κn and κ˜n are i.i.d. (independent and identically distributed) with κn,κ˜n∼N0,1 and υ is a positive constant.The ansatz data are used to the numerical reconstruction. The synthetic Gaussian noise is added to the true forward model, i.e.,(33)y˜=Gqx^,d+δGqx^,dξ1+ξ2i,ξ1,ξ2∼N0,1.In the numerical implementation, we fix wavenumberk=1, δ=0.01 and take υ=2.2,N=5 in (32). For two incident waves Γi=cosπ/2,sinπ/2∪cos3π/2,sin3π/2, we collect the full aperture data, Γo=coss,sins, s∈0,2π. Some shapes are chosen as the test examples. We draw 100 random samples from the approximated posterior distribution, which is displayed in Figure 4. It can be seen that the proposed algorithm provides satisfactory reconstructions for the obstacle scattering problem.(a)
Threelobes:rs=0.5+0.25e−sin3s−0.1sins;(b)
Pear:rs=5+sin3s/6;(c)
Acorn:rs=3/517/4+2cos3s;(d)
Roundedtriangle:rs=2+0.5coss;(e)
Roundrect:rs=cos4s+2/3sins4−1/4;(f)
Kite:x1=coss+0.65cos2s−0.65, x2=1.5sins.Figure 4
The numerical reconstructions for inverse scattering problem using full aperture data with two incident waves: (a)–(c) with incident anglesς=0 and ς=π/2 and (d)–(f) with ς=0 and ς=π. The blue curves are the posterior samples drawn from the approximated posterior distribution. (a) Three lobes. (b) Pear. (c) Acorn. (d) Rounded triangle. (e) Round rect. (f) Kite.
(a)(b)(c)(d)(e)(f)


## 6. Conclusions

In this paper, we apply the homotopy to Bayesian inverse problems. This method provides an optimal approximation for the posterior density. We use the EF and MEF to approximate the posterior in this process. For the EF and MEF, we get the unified derivative formula in the numerical implementation. Especially for Gaussian EF, the dimension of parameters is reduced. We take the KL divergence and Hellinger distance to measure the deviation between distributions and give the HDE of the distribution parameters, which yield the approximated posterior distribution. The approximation of HDE is derived and used to the numerical implementation. Some numerical results show that the proposed algorithm is effective. In the following research studies, some models with multisolutions are to be explored. For the present method, the high-dimension problem is still a challenge.



---

*Source: 9680613-2022-09-10.xml*