Inverse radiative transfer with goal-oriented hp-adaptive mesh refinement: adaptive-mesh inversion

The inverse problem for radiative transfer is important in many applications, such as optical tomography and remote sensing. Major challenges include large memory requirements and computational expense, which arise from high-dimensionality and the need for iterations in solving the inverse problem. Here, to alleviate these issues, we propose adaptive-mesh inversion: a goal-oriented hp-adaptive mesh refinement method for solving inverse radiative transfer problems. One novel aspect here is that the two optimizations (one for inversion, and one for mesh adaptivity) are treated simultaneously and blended together. By exploiting the connection between duality-based mesh adaptivity and adjoint-based inversion techniques, we propose a goal-oriented error estimator, which is cheap to compute, and can efficiently guide the mesh-refinement to numerically solve the inverse problem. We use discontinuous Galerkin spectral element methods to discretize the forward and the adjoint problems. Then, based on the goal-oriented error estimator, we propose an hp-adaptive algorithm to refine the meshes. Numerical experiments are presented at the end and show convergence speed-up and reduced memory occupation by the goal-oriented mesh adaptive method.


Introduction
This paper concerns the numerical solution of inverse radiative transfer equation, which serves as the mathematical foundation for applications such as optical tomography [2,4,32], remote sensing [13,14,29,34], and neutron transport [27,28,30].Despite its wide applications, devising numerical methods for inverse radiative transfer is notoriously challenging because of the high-dimensionality of the forward problem, for which standard discretizations with sufficient accuracy would usually require large memory occupation, and can render the solver too slow to be applicable.Therefore, it is of great interest to devise numerical methods which use fewer degrees of freedoms (DOFs) or less memory, while still maintain the accuracy requirement.
Among the various types of memory reduction methods, we are interested in hp-adaptive mesh refinement (hp-AMR) methods.This is for several reasons.Firstly, the method is versatile in the sense that it can efficiently represent the solution in regions where it is smooth, while also capturing local features.This advantage is especially beneficial for applications such as optical tomography and remote sensing, for which it is common to observe a pattern with a few local features embedded in a smooth background distribution (clouds in the sky for remote sensing; narrow, dirac-delta like inflow laser beam for optical tomography; etc).Another reason concerns the well-developed theoretical understanding of AMR, and its success in the field of solid and fluid mechanics.Indeed, the concept of hp-AMR traces back to at least the 1980s by the pioneering work of Babuška and his collaborators [18]; see also the review paper [5] for hp-finite element methods (FEM).It was shown that hp-adaptive FEM can achieve exponential convergence even when the solution presents singularities.This suggests the great potential of using hp-AMR to reduce the DOFs for radiative transfer.Our recent work [15] demonstrates this potential for the forward problem, while here we consider solving the inverse radiative transfer problem.
Before we proceed to the introduction of our proposed method, we first review some related existing work on inverse radiative transfer, or more generally on inverse transport problems.The well-posedness of inverse transport problems was studied in [6,12,37] by exploiting the singular decomposition of the albedo operators.Then, the stability of the inverse problem in different scaling was studied in [10,11,25,26].For the numerical discretization of the inverse problem, the time-independent inverse radiative transfer was considered in [1,[22][23][24], and the frequency domain problem was considered in [33].Finally, [32] gives a review on the numerical techniques for the inverse transport problems in medical imaging.In [3], the effect of the numerical error on the quality of the reconstruction was discussed, where it was shown that the numerical error on the forward/adjoint solver can lead to significant errors in the reconstruction of the optical images.
Considering the efficiency of AMR in reducing numerical error, it is not surprising to see that AMR has also been applied to inverse problems.In [19,20], the authors considered the diffuse optical imaging problem, where the connection between the forward/adjoint solver error and the reconstruction error was studied theoretically, and an AMR method was proposed based on their theoretical findings.In [35,36], finite volume (FV) based AMR methods for fluorescence/luminescence imaging were studied numerically.In [7,8], a FEM based AMR method was proposed for fluorescence-enhanced optical tomography, and a large savings in DOFs was observed in numerical tests.
In this paper, we propose a goal oriented hp-AMR method which is distinctive to the existing AMR approaches of solving the inverse radiative transfer equation in the following aspects.First, we devise a novel goal-oriented error estimator, which provides more efficient meshrefinement strategies compared to energy or L 2 -norm based error estimators; see [9,39].In addition, we use hp-AMR instead of only h-AMR with fixed p.The extra p-adaptivity allows more efficient representation of the solution where it is smooth.Finally, we consider the full radiative transfer equation instead of its diffusion-type approximation.The full radiative transfer equation is a more computationally challenging task because of the interplay of both the advection and the scattering terms.Since we aim for maximal generality of our method for the different regimes of the radiative transfer equation, we use discontinuous Galerkin (DG) methods to approximate the forward and the adjoint equations.The DG methods can handle well both the advection-dominated and the diffusive regimes, by their flexible choices of the numerical traces.For instance, upwind-type fluxes allow the schemes to suppress unphysical oscillations in the advection-dominated regime, while in the diffusive regime, HDG-type fluxes enable the techniques of static condensation to save more DOFs.Finally, the DG setting allows a straightforward implementation of hp-adaptivity, which lies at the center of our memory reduction method.
The development of our method is based on two key observations.The first is that there exists a goal function, and by minimizing this goal function, the discretized error landscape Φ h (σ) will be sandwiched by the original landscape Φ(σ) in the PDE setting (see ( 4) and (13) for their definitions).This observation motivates the use of this goal function for devising error estimators for AMR.Another observation is that when the numerical error is large, we can use the adjoint solution, which was originally calculated for updating the optical parameters, to also devise the goal-oriented error estimators for doing AMR.This observation suggests that the goal-oriented AMR is especially suitable for the inverse problem setting, since the adjoint solution only needs to be solved for one time but can serve two purposes, namely, the gradient calculation and the error estimator calculation.
The rest of the paper is organized as follows.In section 2, we consider the forward and the inverse problems of the radiative transfer in the PDE setting, and propose an algorithm of reconstructing the optical coefficients.This algorithm will serve as the 'PDE map' for the numerical discretization of the inverse equation.In section 3, we first introduce the DG discretization for the forward and the adjoint problems, and the gradient calculation.Then, we propose the goal function and devise the goal-oriented error estimators.At the end, we combine the discretization and the error estimator and propose an algorithm of reconstructing the optical coefficient with hp-mesh adaptivity.Finally, in section 4, we present numerical experiments to test the performance of our proposed method.

Forward and inverse problems
In this section, we introduce the forward and the inverse problems of the radiative transfer equation.Then, an algorithm (algorithm 1) is proposed to reconstruct the optical coefficients in the PDE setting.

Forward problem
Let Ω ⊂ R d be a bounded Lipschitz domain, and S represent the unit sphere in R d .Let ∂Ω be the boundary of Ω, and n the outward-pointing unit normal vector on the boundary.We denote by the outflow and the inflow boundary of Ω × S, respectively.We consider the following (timeindependent) radiative transfer equation: find u ∈ V such that In the above equation, u(x, s) is the radiative intensity at spatial location x and along the direction s ∈ S, while σ a (x) 0, σ s (x) 0, f(x, s), and g(x, s) are the absorption coefficient, the scattering coefficient, the source term, and the inflow radiation, respectively.We assume the scattering phase function P(s, s ′ ) to have the form of the Henyey-Greenstein function: where g is the asymmetric parameter, and c is chosen such that ´S P(s, s ′ )ds ′ ≡ 1.Here V is the solution space where V g is the affine space of V satisfying the inflow condition (1b).Note that V 0 is a subspace of V.
To render the presentation more concise, we rewrite (1) into a more compact form.Let F = (f, g) T , and L(σ) = (D(σ), γ − ) T , where γ − is the trace operator to L 2 (Γ − ), and σ represents the collection of both σ e and σ s , namely σ(x) := (σ e (x), σ s (x)).Then the equation (1) transforms into In the equation (1) (or (3)), we use the square bracket [•] to indicate that D(σ) (or L(σ)) depends linearly on u; the round bracket suggests that D (or L) can depend non-linearly on σ.In the setting of the forward problem, the optical property represented by σ and the source/boundary term F are given as input data.Then the task it to solve the equation (1) (or (3)) for the radiative intensity u.We refer to [16] and the references therein for the well-posedness of the forward problem.

Inverse problem
In the setting of the inverse problem, we are given certain measurements of the solution, here denoted as Mu, while the task is to infer the distribution of the optical property σ or the source term f.The former task fits into the field of optical tomography, while the latter one is usually referred to as inverse source problems.We refer to [2,6] for a review on these topics.
In this paper, we focus on the optical tomography problem.To be more specific, we consider the task of reconstructing the optical property of the medium σ, based on a set of tests determined by the source terms F i with i = 1, . . ., N t , and a set of measurements on the outflow boundary Γ + , here denoted as M j with j = 1, . . ., N m .This problem can be formulated as an optimization problem.Suppose σ * is the target optical coefficient we would like to reconstruct.
We define y ij = M j u i = M j (L(σ * )) −1 [F i ] as the measurement data.For many applications, the data are polluted by noise.In such a case we define y ij := M j u i + δ ij where δ ij represent the noise.The optical tomography problem becomes where w ij are the weights associated to each test-measurement signal y ij .A typical choice is w ij = 1 |y ij | 2 ; see [32].Here, αR(σ) is the regularization term and α is the regularization parameter.Some widely used examples of regularization include the L 2 -regularization, for which we choose R(σ) = ∥σ − σ∥ L 2 (Ω) (where σ is a spatial average), or H 1 -regularization, for which we choose R(σ) = ∥∇σ∥ L 2 (Ω) .For more on regularization techniques and the choices of regularization parameters, we refer to the monograph [17] and the references therein.
Note that Φ(σ) is nonlinear and there is no explicit expression for Φ(σ) in general cases.To proceed, we transform (4) into a PDE-constrained optimization problem: subject to L (σ) To solve the problem, we introduce the Lagrangian where v i is the Lagrange multiplier.By the Karush-Kuhn-Tucker (KKT) necessity condition [38], the solution of (5a), if exists, is also the first-order stationary point of Ψ in (6).By taking Fréchet derivatives of (6), we obtain where ⟨•, •⟩ is the trial-test bracket, manifested as the summation of the L 2 -inner product on Ω × S and Γ − , and L(σ) * is the adjoint operator of L(σ).Now, the optimization problem reduces to finding the stationary point (σ * , u * , v * ) such that equations (7) equal 0. Based on (7), we propose the following algorithm 1 of finding σ * .In the algorithm, the forward and the adjoint problems associated with (7c) and (7b) are first solved to obtain u i and v i , which are then used to update σ by an updating direction δσ k and the step length γ > 0. The step length γ can be determined by a line search [31].The updating direction can be simply chosen as the gradient given by (7a), or by Quasi-Newton method such as the BFGS method.We refer to [31] more details.
1: Set the maximal iteration count N max_iter and the error tolerance ϵ > 0 2: Initialize σ 0 3: Set k = 0 4: while k N max_iter do 5: Using σ k , solve (7c) = 0 to obtain u i 6: Using σ k and u i , solve (7b) = 0 to obtain v i 7: Using σ k , u i , v i , and (7a) to obtain the updating direction δσ k 8: Update the optical coefficients by σ k+1 = σ k − γ δσ k , where γ is determined by a line search 9: Update the iteration count: Break the loop 12: end if 13: end while

Numerical discretization
In this section, we introduce the numerical approximations to the inverse problem (4), manifested by using the discontinuous Galerkin spectral element (DGSE) methods to discretize the forward problem (7c) = 0, the adjoint problem (7b) = 0, and the gradient calculation (7a).Then, we propose a goal-oriented error estimator which can be easily calculated based on the forward and the adjoint numerical solution.Finally, based on the error estimator, we introduce an hp-AMR algorithm to solve the inverse radiative transfer problem numerically.We begin by considering the approximation spaces.For simplicity, we restrict ourselves to consider only a 2D domain Ω discretized by rectangular meshes T h .The generalization to 3D and polyhedral meshes is straightforward under the DG setting.Since any K ∈ T h is rectangular, we can write K = K x × K y .We require that the mesh T h satisfies the one-irregular condition, i.e. for any element K ∈ T h , there are at most two neighbour elements connecting to K through each edge.Here h is both an index for the triangulations T h and the mesh-size, which is defined as h := max K∈T h h K .For each K ∈ T h , we denote by T a,K h the triangulation of the unit sphere S. We also let F K be the collection of the faces of K. Now we introduce the approximation space: are a set of polynomial bases on K ⋆ , where p K ⋆ represents the polynomial degree and ⋆ ∈ {x, y, a}.Here we choose φ is the Lagrange polynomial basis associated to the ith Gauss-Legendre-Lobatto quadrature points, and + a is the push-forward map.This choice of basis gives us spectral element methods.Note that the space V h is completely determined by the mesh discretization T h and the polynomial degrees p K a , p K y , and p K a which are associated to each spatial-angular element K × K a .We will sometimes write T hp to denote the mesh T h encoded with also the polynomial degree information.Finally, to incorporate the boundary condition, we introduce As a special case, V 0 h is a subspace of V h with zero inflow radiation.For notation conciseness, we define as the L 2 inner product on K × K a , and we denote the associated norm as K×K a .Also, when there are two elements K + and K − which share a common face F, we introduce the jump notation 3.1.2.DG for the forward and adjoint problems.Now we are ready to introduce the DG numerical discretization for the radiative transfer equations.We first define the following bilinear form: where we take the upwind numerical trace u h : and u nbr h is the restriction of u h from the neighbour elements of K on F. The DG approximation for the ith test of the forward problem (7c) reads as On the other hand, the DG approximation for the ith test of the adjoint problem (7b) reads as where are the measurements.For different tests, we choose different underlying hp-meshes.We denote by T i hp and T ⋆−i hp as the underlying hp-mesh for the ith test of the forward and the adjoint problem, respectively.We shall assume that T ⋆−i hp is obtained by refining T i hp .For instance, T ⋆−i hp can be a one-level h or p-refinement of T i hp .This extra refinement for the adjoint meshes will be essential for devising the error estimators for the inverse problem, and will be discussed in more detail in section 3.2.2.
3.1.3.Gradient calculation in the discrete setting.For the rest of this subsection, we consider how to use the forward and the adjoint numerical solutions, namely u i h and v i h , to calculate the gradient and update the optical parameter in the DG setting.Namely, we consider the DG discretization for the gradient calculation (7a).
We begin by introducing the approximation spaces for the optical parameter σ.Since we use the DG spectral element method as our radiative transfer equations solver, it is natural to choose the discretized optical parameter σ h to live in a DG spectral element space.We shall denote by T σ hp as the underlying mesh for this space, and by Σ h , or Σ h (T σ hp ) as the space.Here we adopt the parametric reconstruction method [32] as a regularization to relieve the ill-posedness of the original inverse problem (4).This means that we shall set T σ hp to live in a coarse mesh with large mesh-size.In the extreme case, we choose T σ hp to be the one-element mesh covering the whole domain Ω.Then, the polynomial degree p can be adjusted in the sense that a lower degree p represents a stronger regularization effect (restriction to the lowfrequency modes) and vice versa.Now, the discretized version of (7a) becomes for any test function δσ h ∈ Σ h , where u i,s h := ´S P(s, s ′ )u i h (s ′ )ds ′ , and u i h and v i h are the solution of ( 9) and (10), respectively.By letting the test function δσ h go over all the basis polynomials of Σ h , we obtain a vector as the updating gradient for σ h .
We conclude section 3.1 by the following algorithm 2.
Break the loop 14: end if 15: end while 3.2.Goal-oriented error estimator 3.2.1.The goal function.The numerical solver ( 9) and ( 10) can be rewritten in a compact form as follows: where DGSE and DGSE * represent the forward and the adjoint solver, respectively.Using this notation, the discretized inverse problem (4) becomes min where I := {1, . . ., N t } is the index set for all tests.Note that the discretized error landscape functional Φ h depends both on the optical property σ and the set of the hp-meshes , then the error landscape Φ(σ) can be related to its discretized version Φ h (σ, T hp ) as follows: As a result, when Err h (σ) → 0, we have The above estimate suggests that, when the error Err h (σ) goes to zero, the minimizer of Φ(σ) also minimizes Φ(σ, T I hp ) and vice versa.This observation suggests that we should aim for devising numerical methods such that Err h (σ) is reduced in the most efficient way.We next explain how this can be achieved in more detail.
Toward this aim, for simplicity, write u i (σ Recall that we denote by σ * the minimizer of Φ(σ) so we have y ij = M j u i (σ * ) by definition.Now, for each test i, we introduce the functional J i (σ): Then, it holds that The above identity suggests that we can use the functional J i (σ) to devise error estimators to guide the mesh-refinement of T i hp .This approach fits into the well-known framework of goaloriented error control and AMR; see [9].To be more specific, the error estimators are obtained by solving the dual equation (10) except the right-hand-side term is different.Now note that Therefore, when |y ij − M j u i (σ)| ≪ |y ij − M j u i h (σ)|, or namely, when the numerical error dominates the error caused by optical parameters, we have . Thus, in this case, we can use M i (u i h (σ)) as an alternative of J i (σ) to devise error estimators to guide the mesh refinement of T i hp .As a result, we have z i h ≈ v i h .This suggests that we do not need to solve another adjoint equation for z i h but can simply use v i h to devise the error estimators.
3.2.2.The error estimators.For the rest of this subsection, we show how a goal-oriented error estimator can be devised based on a given functional J in the DG setting.Namely, we would like to control |J(u − u h )| by a posteriori error estimators.To make the presentation concise, we shall consider the forward problem (1), and its DG discretization formulated as follows: find u h ∈ V g h such that Before we present the main result, we first introduce two lemmas on the consistency and the Galerkin orthogonality of the DG bilinear form (8a).
Lemma 3.1 (consistency of the DG bilinear form).Suppose u solves (1).Then Proof.By multiplying (1a) with v * and integrating on K × K a , we obtain Since u is the exact solution, we have (s • n)u = (s • n) u by the consistency of the numerical trace, see (8b).Then we take summation for all K ∈ T h and K a ∈ T a,K h .This completes the proof.
Proof.We simply choose v * ∈ V 0 h in (15) and then take the difference between ( 14) and ( 15).This completes the proof.
Consider the following dual equation: find z ∈ V 0 such that Then we define where In the above equations, ϕ h can be any function in V 0 h .Proposition 3.1.Let u ∈ V g be the solution of (1), u h ∈ V g h solves (14), and z ∈ V 0 be the solution of the dual problem (17).Then, we have Proof.Since u h ∈ V g h and u ∈ V g , we know (u − u h ) ∈ V 0 .Thus, we can let v = u − u h in (17), which gives us , where the Galerkin orthogonality identity ( 16) is used for the second equal sign.Then, by the definition of the bilinear form (8a) and its consistency property (15), we can proceed as follows: This completes the proof.
Proposition 3.1 suggests that we need to solve a dual equation (17) to calculate the error estimators.Since a direct solve for z is not feasible, we calculate an approximate solution to z by z h .Here z h should live in a more refined space than the space for u h .Otherwise, we can choose ϕ h = z h and this will render the error estimators vanishing.To be more specific, we seek z h ∈ Ṽ0 h such that Here Ṽh is a more refined space than V h .For instance, we can choose Ṽh to be an h-refined or prefined version of V h .In this paper, we shall use the p-refined version since it is straightforward to implement.Once z h is calculated, we choose ϕ h = Π V z h where Π V is the L 2 projection onto the space V h .Namely, the weighting terms become For the rest of this section, we use the error estimators derived here to design an hp-AMR method to solve the inverse radiative transfer problem (4).

Goal-oriented mesh adaptation for the inverse problem
Up to this point we have introduced the forward solver ( 9), the adjoint solver (10), the gradient calculation (12), and the error estimators (18).Next, we show how they can be combined to solver the inverse problem ( 13) with mesh adaptation.We shall first go over the basic elements of the hp-AMR and then propose the full algorithm.

hp-AMR.
Given a mesh T hp , we can refine it based on the error estimators introduced in (18).For h-refinement only, ( 18) is sufficient for an AMR algorithm.For hp-AMR, we will need an additional smoothness estimators to complete the algorithm.Since we use a more refined mesh for the duality solution v h , the dominating error comes from u h .Hence, here we use u h for devising the smoothness estimator.We consider the smoothness estimators proposed in [21].These estimators are obtained by examining the decaying pattern of the coefficients obtained from a Legendre expansion of the numerical solution.To be more specific, given the forward solution u h , we first calculate its mean intensity u mI h := ´S u h .Then, for any element K ∈ T h , we calculate the Legendre expansion coefficients of u mI h on K, here denoted as a K i,j with i = 0, . . ., p x and j = 0, . . ., p y , where p x and p y are the polynomial degrees of the tensor-product basis function on K in the x and the y directions, respectively.Then, the smoothness estimator is defined as follows: With the error estimator η K (defined in (18)) and the smoothness estimator η s K (defined in ( 19)), we propose the following algorithm 3 to refine a given mesh T hp .Algorithm 3. hp-adaptive mesh refinement.
1: Sort all element K ∈ T h in increasing order according to the error estimator η K defined in (18).2: Mark the largest r ref percentage of the elements for refinement-T mark h .3: For all K ∈ T mark h , calculate the smoothness estimators η s K defined in (19).4: For all K ∈ T mark h , we perform the following refinement 5: if p K + 1 < η s K − 1 2 then 6: Perform p-refinement for K. 7: else 8: Perform h-refinement for K. 9: end if Note that for each test i, the error estimators η i K are computed based on the forward solution u i h and the adjoint solution v i h ; see equation (18), where z is replaced by v i h .
3.3.2.The full algorithm.We conclude this section by proposing the algorithm 4 of solving the discretized inverse problem ( 13) with mesh adaptation.
1: Set T σ hp , the mesh for the optical parameter 2: Set the maximal iteration count N max_iter , the maximal DOFs count N dofs , and the error tolerance ϵ > 0. 3: For each test i, initialize T i hp , the mesh for solving the forward and the adjoint problems 4: Initialize the optical parameter σ 0 h .5: Set k = 0 6: while k N max_iter and DOFs N dofs do 7: Using σ k h , solve (9) 13)) then 13: For each test i, refine the mesh T i hp by algorithm 3 14: Reinitialize the optical parameter end if 16: end while Note that previously in algorithm 2, we stop the iteration when the relative difference Here in algorithm 4, we apply mesh-refinement instead of stopping the iteration.This procedure is repeatedly applied until the maximal DOFs count N dofs is reached.This can prevent the algorithm from refining indefinitely.The constant N dofs can be chosen according to the memory-occupation restriction or computational time restriction.The parameter λ determines how much information will be inherited from the optical parameter recovered in the previous mesh setting.If the previous mesh setting can provide a relatively good approximation of σ, then we can increase λ for a better initial guess.Otherwise, we choose λ = 0 to avoid potential bias from an inaccurate numerical solution.In the numerical tests we carried out in this paper, it only takes a few steps for σ k h to converge, so we always choose λ = 0.

Numerical experiments
In this section, we carry out numerical experiments to test out the performance of our proposed numerical methods.In the first subsection, we test the goal-oriented AMR method for the forward problem.Then, in the second subsection, we test the performance of the method for the inverse problem.

Goal-oriented AMR for the forward problem
In this subsection, we consider the numerical approximation to the forward problem (1).
Before testing the goal-oriented AMR methods, we first test the convergence of the DGSE method without adaptivity.We consider a rectangular spatial domain [0, L x ] × [0, L y ] with the inflow radiation coming from the top and the left boundary: For the other sides of the boundary, we assume zero inflow radiation and zero source term f.A scatterer is placed at the center of the domain, so that the extinction coefficient σ e is We choose the single scattering albedo to be ω = 10  11 , so the scattering coefficient is σ s = ωσ e and the absorption coefficient is σ a = σ e − σ s .We choose k 0 = 2.The asymmetric parameter in the scattering phase function in ( 1) is chosen as g = 0.8, which is a typical value for water clouds in the atmosphere.
For this test we consider the DG method ( 9) with the polynomial degree p x = p y = 1 and also with p x = p y = 2 for the spatial discretization.To compare, we also consider an FV spatial discretization method, which is equivalent to the DG method with p x = p y = 0.For all these methods, we shall uniformly partition the domain [0, L x ] × [0, L y ] into a collection of elements and slowly decrease the element size.For the angular discretization, we apply a uniform partition into 32 pieces by P 0 element (FV angular discretization).Figure 1 shows that the DG method has errors that converge as O(DOFs −(p+1)/d ) where the dimension d = 2.This is the expected optimal convergence order for the DG method and the FV method (upwind).
We next consider the forward problem with the goal-oriented AMR algorithm.We consider a rectangular spatial domain [0, L x ] × [0, L y ] with the inflow radiation coming from the top boundary: For the other sides of the boundary, we assume zero inflow radiation and zero source term f.For this test, we also consider a scatterer placed at the center of the domain following the form of equation (20).To better test the AMR algorithm, here we consider the scatterer with a sharper boundary transition by taking k 0 = 10.The other settings are the same as the previous test; namely, we choose ω = 10/11 and g = 0.8.See figure 2 for a visualization of σ e .Our aim is to find an efficient numerical approximation to J(u), where u is the solution of (1).Here we set the goal function J to be (a smoothed version of) a point measurement on the outflow radiation, located at the right side of the domain:   Goal-oriented AMR versus standard AMR.The plots show the mean intensity of the radiant intensity u h and the meshes.Left: using the goal-oriented error estimator η K introduced in (18).Right: using the standard jump estimator η jump K defined in (21).The standard jump estimator induces refinement near sharp gradients, whereas the goaloriented error estimator also induces refinement along the path between the scatterer and the measurement location on the boundary at (x, y) = (L x, (15/21)Ly).
Specifically, the (smoothed) point measurement is put at (L x , 15  21 L y ).In the first test, we use the error estimator η K defined in (18) to guide mesh refinement.The mesh refinement algorithm 3 is repeatably applied with the refinement ratio r ref = 0.2, until the DOFs for u h reach the threshold N max_dof = 2 × 10 5 .To better visualize where the meshes are refined, here we use a fixed polynomial degree p = 3.
To compare, we also consider a standard jump estimator Since the exact solution u is not known, we calculate |J(u h ) − J(ũ)| as an approximation to the true functional error |J(u h ) − J(u)|, where ũ is calculated on a more refined mesh than u.Here ũ is calculated with an extra p-refinement for each element K ∈ T h .Figure 3 shows the numerical solutions and the meshes calculated based on the goaloriented estimator η K and the standard jump estimator η jump K .We observe that the goal-oriented estimator can successfully guide the mesh to refine at the place where the solution has a sharp gradient, namely, top of the scatter.In addition, it refines along the path connecting the place of the sharp gradient and the place of the point measurement (located at the right side boundary at y = 15  21 L y ).In contrast, the standard estimator η jump only refines where there are sharp gradients, without taking into consideration of the effect of the measurement J.
In figure 4, we plot the error |J(u) − J(u h )|.We observe that the goal-oriented estimator reduces the error efficiently while the standard estimator fails to decrease the error.This is consistent with our observation in figure 3, where the goal-oriented estimator refines both the places of the sharp gradient and the path connecting to the measurement, while the standard estimator fails to refine those elements concerning the measurement J.The experiment shows that the goal-oriented estimator can be much more efficient than the standard one, when the goal function is not an L 2 or energy-based norm, but instead a measurement on the domain boundary.Note that this is exactly the case for the applications such as optical tomography, which we shall consider in the next subsection.

Goal-oriented AMR for the inverse problem
In this subsection, we consider the numerical tests for the inverse problem (13).We aim to test the performance of algorithm 4 and compare it with other mesh-refinement methods.We consider a rectangular domain [0, L x ] × [0, L y ] where on each side, we put 2 inflow radiation laser beams and 20 measurements collecting the angularly-averaged outflow radiation.As a result, we have N tst = 8 and N msm = 80 in the formulation (4) or (13).The absorption is fixed as σ a = 0.1 and the asymmetric parameter is chosen as g = 0.1.The scattering coefficient σ s is decomposed into a summation σ s = σ 0 s + σs where σ 0 s is a background state and σ represents a perturbation.We assume the background state σ 0 s and the absorption σ a are known and we aim to recover the scattering coefficient σ s .We shall consider two test cases with the same background state σ 0 s = 1 but with different perturbations σs .See the upper-left sub-figures of figures 6 and 9 for the true scattering coefficients σ s for the first and the second test case, respectively.
To generate the data y ij = M j ũi h , we solve the forward problem (9) on a very refined mesh, namely, a 8 by 8 mesh with the polynomial degree chosen as p = 9, so we can regard ũi h as the 'true' solution.The optical parameter σ h is set to live on the one-element mesh with p = 19.
To solve the inverse problem, we apply an H 1 -regularization, namely R(σ) = ∥∇σ∥ L 2 (Ω) , with the regularization parameter chosen as α = 10 −1 ; see (4) where α was first introduced.We implement algorithm 4, with ϵ = 10 −3 , and r ref = 0.2 for the mesh refinement algorithm 3.For the iteration method, we use the limited memory BFGS algorithm to calculate the updating direction δσ k h .The parameter γ is determined by a line search method, for which we use the backtracking algorithm which starts with µ = 1 and then followed by the candidate step length µ = 2 −1 , 2 −2 , 2 −3 , . . ., until the Armijo condition is satisfied, where c 1 = 10 −4 .For more details on the BFGS method and the line search algorithm, we refer to [31].
We test the performance of algorithm 4 by comparing it with other mesh-refinement methods.We set the initial mesh T i hp (i = 1, . . ., N tst ) by a 4 by 4 mesh with p = 2.To be more specific, we shall compare algorithm 4 with other methods, for which, instead of applying the AMR algorithm 3 in the algorithm 4, we apply (1) h-refinement, for which we partition the domain into more elements (uniformly divided), ( 2) p-refinement, for which we increase the polynomial degree on each element, (3) h-AMR-goal, for which we only apply h-refinement in the algorithm 3, and (4) hp-AMR-standard, for which we use the standard jump error estimator (18) to guide the mesh-refinement in the algorithm 3. We remark that the h-refinement method with p = 2 is similar to a FV method as mesh resolution is increased, while the prefinement method with an increased polynomial order p is similar to a spectral method in the sense that an exponential (spectral) convergence can be achieved when the solution is sufficiently smooth.Now we show the results for the first test case.Figure 5 shows that the hp-AMR-goal method is the most efficient method in reducing the error ∥σ h − σ * h ∥ L 2 (Ω) .Namely, it uses the fewest DOFs to achieve the best recovery.The h-AMR-goal method is overall less efficient than the hp-AMR-goal method but their performances are similar.For the hp-AMR-standard, h-Ref, and p-Ref methods, we observe their errors first decrease and then increase, and in all DOFs levels, their recovered optical property is poorer than the one recovered by the hp-AMR-goal method.The error for the hp-AMR-standard method decreases fast at the beginning but then the error increases to back to the original value, but even in its decreasing phase its performance is less efficient than the hp-AMR-goal method.This consolidates our claim that the goal-oriented estimator is necessary for inverse problems.For h-Ref and p-Ref methods, we observe that they decrease the errors much more slowly than hp-AMR-goal.In summary, the hp-AMRgoal method behaves better than other methods in almost any measure, and the h-AMR-goal method behaves less well but similar to hp-AMR-goal.
To see this more clearly, we plot in figure 6, the recovered scattering coefficients by the different refinement methods.To make a fair comparison, for all the refinement methods, we take the snapshot at the refinement step such that the DOFs are closest to 10 6 .The results agree with what we observe in figure 5. Namely, the hpAMR-goal method gives the best approximation to the true scattering coefficient, while using the fewest number of DOFs.
As additional detail of what is seen in figure 6, we observe that the hRef method and hpAMR-standard method both fail to provide a satisfactory recovery to the true distribution (left-top sub-figure).All of the pRef, hAMR-goal, and hp-AMR methods can successfully recover the correct location and the approximate shape of the perturbation σs .In further detail, though, the h-AMR-goal and hp-AMR-goal methods provide a more accurate recovery for the amplitude of perturbation in comparison to the p-Ref method.Note that here we chose a fixed regularization parameter α for all refinement methods for a fair comparison.Despite this, we remark that if it is allowed for the regularization parameter α to be adapted properly to the different refinement algorithms, the recovery for pRef, hAMR-goal, and hpAMR-goal might further improve.
To visualize how the meshes are refined by the different methods, we plot the forward solution and the corresponding meshes in figure 7.For goal-oriented hp-AMR method, we observe that the mesh is refined both at the place where the forward solution has a sharp gradient, and also at the left boundary where the measurement is located.In contrast, the standard hp-AMR method only refines the mesh where the sharp gradients are located, but fail to refine also at the measurements.The goal-oriented h-AMR method refines both at the sharp gradients and at the boundary where the measurements are located.However, it does not use a high order approximation as in the hp-AMR cases.This could explain why its performance is overall less great than hp-AMR-goal.Therefore, the goal-oriented hp-AMR is the only method that supports both (1) a balanced refinement strategy between resolving the forward and the adjoint solution, and (2) an automatic high-order approximation by p-adaptivity.Considering this, it is not surprising to see that the hp-AMR-goal method recovers the optical coefficient with the best quality in figure 6.
To test the reliability of our proposed method, we conduct another experiment (the second test case) with the same setting used for figure 5 but with a different scattering coefficient σ s to be recovered; see the top-left figure of figure 9 for a visualization of σ s .From figures 8 and 9, we again observe that our proposed goal-oriented hpAMR method has the best performance compared to all the other refinement approaches.Namely, it uses the fewest DOFs but provides the best quality of the recovered scattering coefficient σ s .
To test the robustness of the proposed goal-oriented hp-adaptive algorithm, we conduct additional numerical tests with noise added to the signals y ij .Namely, instead of using y ij , we use ỹij := y ij (1 + δX ij ) where X ij are independent, identically distributed random variables from a uniform distribution on [−1, 1], and δ represents the noise level.We consider two cases of δ = 1% and δ = 10%.
Figure 10 shows that when the noise level δ = 1%, the performance of the different refinement methods are similar to the noise-free case (figure 5).We again observe that the hp-AMRgoal methods gives the smallest error while using the fewest DOFs.When the noise level δ increases to 10%, we observe an increase of the errors for most refinement methods.Despite this, the hp-AMR-goal method still provides the overall smallest error compared to the other refinement approaches, especially in the last refinement level.In figures 11 and 12, we plot the recovered scattering coefficient by the different refinement methods with the noise level 1% and 10%, respectively.The results are consistent with what we observe in figure 10.We repeat the noise test for the second test case (with the true scatterer shown in the upperleft sub-figure of figure 9).In figure 13, we observe that in both cases of the noise level of δ = 1% and δ = 10%, the hp-AMR-goal method provides the overall smallest error.We plot in figures 14 and 15 the recovered scattering coefficients.The results are consistent with what      we observe in figure 13 in the sense that the hp-AMR-goal methods provides the best quality of the recovery in the different noise levels.

Conclusion
In this paper we propose a goal-oriented hp-AMR method to solve the inverse radiative transfer equation.The method is based on the development of a novel goal-oriented error estimator, which is achieved by connecting two kinds of duality arguments in different fields, namely, (1) the duality-based based mesh adaptivity for goal-oriented error minimization, and (2) the adjoint-based inversion techniques for solving inverse problems.The numerical tests suggest that the proposed method solves the inverse problem with the best quality of the recovered optical coefficients, while using the fewest DOFs.While our method is proposed here for inverse radiative transfer, the general principles of devising the error estimators and designing the refinement algorithms should be able to be naturally extended to enable adaptivemesh inversion for other types of inverse problems, such as the Calderón problem or inverse scattering.
Potential future work includes further study of combinations of more involved regularization strategies, and applications to remote sensing problems.

Algorithm 2 . 2 : 5 :
Reconstructing optical coefficients-DG discretization.1: Set the mesh for the optical parameter-T σ hp Set the mesh for each test i-T i hp for the forward and T ⋆−i hp for the adjoint problem 3: Set the maximal iteration count N max_iter and the error tolerance ϵ > 0 4: Initialize σ 0 h Set k = 0 6: while k

Figure 1 .
Figure 1.Convergence test for the forward problem, without AMR.Left: mean intensity of the radiative solution calculated by the DG-P 2 method.Right: L 2 error of the mean intensity calculated by the finite volume method, DG with P 1 , and DG with P 2 polynomial approximation.All methods converge at the expected rates of O(DOFs −(p+1)/d ) with dimension d = 2, and the DG methods with larger p converge faster.

Figure 2 .
Figure 2. Extinction coefficient σe, to represent an idealized cloud located at the center of the domain.

Figure 3 .
Figure 3. Goal-oriented AMR versus standard AMR.The plots show the mean intensity of the radiant intensity u h and the meshes.Left: using the goal-oriented error estimator η K introduced in (18).Right: using the standard jump estimator η jump

Figure 4 .
Figure 4. Evaluation of the error |J(u h ) − J(u)| by the goal-oriented and the standard error estimator.The goal-oriented estimator leads to a much smaller error.

Figure 5 .
Figure 5.Total DOFs of the forward and the adjoint solutions (the summation of all tests) versus the L 2 error of the optical property ∥σ h − σ * h ∥ L 2 (Ω) for the first test case.

Figure 6 .
Figure 6.Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the first test case.L 2 errors are printed at the upper-right corners for each sub-figures.

Figure 7 .
Figure 7. Top row: forward solution u 1 h for the first test case, by goal-oriented hpAMR (left), goal-oriented hAMR(middle), and standard hpAMR(right) at the last step of the refinement.Bottom row: the corresponding meshes T 1 hp , where the colorbar represents the polynomial degrees.

Figure 8 .
Figure 8.Total DOFs of the forward and the adjoint solutions (the summation of all tests) versus the L 2 error of the optical property ∥σ h − σ * h ∥ L 2 (Ω) for the second test case.

Figure 9 .
Figure 9. Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the second test case.L 2 errors are printed at the upper-right corners for each sub-figures.

Figure 10 .
Figure 10.Total DOFs of the forward and the adjoint solutions (the summation of all tests) versus the L 2 error of the optical property ∥σ h − σ * h ∥ L 2 (Ω) for the first test case.Left: 1% noise is added to y ij .Right: 10% noise is added to y ij .

Figure 11 .
Figure 11.Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the first test case.L 2 errors are printed at the upper-right corners for each sub-figures.1% noise is added to y ij .

Figure 12 .
Figure 12.Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the first test case.L 2 errors are printed at the upper-right corners for each sub-figures.10% noise is added to y ij .

Figure 13 .
Figure 13.Total DOFs of the forward and the adjoint solutions (the summation of all tests) versus the L 2 error of the optical property ∥σ h − σ * h ∥ L 2 (Ω) for the second test case.Left: 1% noise is added to y ij .Right: 10% noise is added to y ij .

Figure 14 .
Figure 14.Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the second test case.L 2 errors are printed at the upper-right corners for each sub-figures.1% noise is added to y ij .

Figure 15 .
Figure 15.Recovered scattering coefficients σs for the different mesh refinement methods, with the total DOFs at around 10 6 for the second test case.L 2 errors are printed at the upper-right corners for each sub-figures.10% noise is added to y ij .