{'fraction_non_alphanumeric': 0.03798230301571059, 'fraction_numerical': 0.013724191897911275, 'mean_word_length': 4.299362041467305, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Iterative Algorithms for Quantitative Photoacoustic Tomography in the Radiative Transport Regime

In this paper, we describe the numerical reconstruction method for quantitative photoacoustic tomography (QPAT) based on the radiative transfer equation (RTE), which models light propagation more accurately than diffusion approximation (DA). We investigate the reconstruction of absorption coefficient and/or scattering coefficient of biological tissues. Given the scattering coefficient, an improved fixed-point iterative method is proposed to retrieve the absorption coefficient for its cheap computational cost. And we prove the convergence. To retrieve two coefficients simultaneously, Barzilai-Borwein (BB) method is applied. Since the reconstruction of optical coefficients involves the solution of original and adjoint RTEs in the framework of optimization, an efficient solver with high accuracy is improved from~\cite{Gao}. Simulation experiments illustrate that the improved fixed-point iterative method and the BB method are the comparative methods for QPAT in two cases.


Introduction
Photoacoustic tomography (PAT) is a developing medical imaging technique using non-ionizing waves. It combines ultrasound imaging and optical imaging [29], so it can achieve high contrast and high resolution simultaneously. Photoacoustic imaging is based on photoacoustic effect. A short pulse of electromagnetic energy propagates through the tissue rapidly at speed of higher magnitude than sound after the tissue is illuminated. Then the electromagnetic energy is absorbed by the tissue and converted to heat. There quickly follows the temperature of tissue increasing which leads to expansion spatially, and the expansion induces a pressure field. Acoustic waves generated by pressure distribution will be measured by detectors placed at the surface. They can be used to extract physical information about tissue inside. In medical imaging science, this technique is mostly applied to tomography imaging of skin, breast cancer and small animal imaging [29].
In the physical view, there are two inverse problems associated with two stages: the first stage of generating pressure after absorbing electromagnetic energy and second stage of initial pressure generating acoustic waves. First problem (aka PAT) is to reconstruct initial pressure distribution from temporal measurements of photoacoustic waves, and the second one, called quantitative photoacoustic tomography (QPAT), is to estimate optical coefficients from initial pressure distribution. Although initial pressure reflects some information inside, it is the product of the absorption of energy relied on properties of the object indirectly. Only the intrinsic optical coefficients provide direct information. Essentially, spatially varying optical coefficients lead to spatially varying initial pressure distribution. Moreover, reference [22] mentions that initial pressure decreases quickly in the tissue area with small absorption coefficient, so that we can not distinguish visually in the finial image of initial pressure. In a word, the original information in the investigated object can not be obtained completely until the optical coefficients are recovered. Generally, estimated optical coefficients comprise of absorption coefficient, scattering coefficient and Grüeneisen coefficient, where absorption coefficient is usually of the major clinic interest [14]. Nowadays, algorithms on reconstruction of initial pressure have been developed maturely [20,29,28,1,21], so a lot of researchers are focusing on QPAT.
In terms of mathematical model of QPAT, light propagation inside the tissue can be modeled by radiative transfer equation (RTE) [25,22,18] or diffusion approximation (DA) [4,12]. However, DA approximates to RTE only in the situation where light behaves diffusely, which is invalid for small object or in regions close to light source, so that the optical coefficients estimated from DA model are not accurate. The reconstruction theories and some algorithms based on DA model are presented in [23,24,5,13,31,7]. The two models are compared numerically in [27] and the derivation from RTE model to DA model can be found in [10].
In terms of measurement, there are mainly three kinds of methodologies: single illumination, multi-source illumination and multi-spectral illumination. It has been proved that single illumination can recover any one of absorption and Grüeneisen coefficient stably given scattering coefficient, and even there is an explicit analytic formula for the non-scattering media [22]. Multi-source can recover two of the three coefficients mentioned above uniquely and simultaneously [22,18,3,26]. Furthermore, multi-spectral illumination can achieve unique reconstruction of three coefficients [2]. We focus on the RTE model with multi-source illumination and aim to investigate more efficient reconstruction method in this setting.
In term of reconstruction algorithm, assuming the knowledge of scattering coefficient, fixed-point iteration is used to recover absorption coefficient in [19]. However, the requirement of close enough initial value and the lack of theoretical results are obstacles of its application. When the scattering coefficient is unknown, QPAT can be formulated as a nonlinear optimization problem [18,14] and furthermore can be simplified to a linear optimization problem by Born linearization [22]. Optimization approaches such as Jacobian-based method and gradient-based method are investigated in reference [25], in which Limited-memory BFGS (LBFGS), as a state of art gradient-based method, has been applied to the problem [14].
In the optimization framework, it is inevitable to solve original RTE in forward problem and adjoint RTE in the process of calculating gradient of objective function. Finite element method combined with streamline diffusion modification is applied to solve stationary RTE [18,25,30]. In this way, a large sparse linear system needs to be solved, which is still difficult. Inspired by [15], we apply Discontinuous Galerkin (DG) method combined with multigrid method to solve the two RTEs.
Our contributions in this paper are in several folds.
(i) We propose an improved fixed-point iterative algorithm and prove its convergence given scattering coefficient. Numerically, the improved fixed-point iteration converges faster than the optimization-based Barzilai-Borwein (BB) method.
(ii) We apply BB method to reconstruct absorption coefficient and scattering coefficient simultaneously. To the best of our knowledge, it is the first time to apply BB method to this problem. (iii) We revise the numerical scheme for RTE proposed in [15], which fails to solve adjoint RTE. We revise it to make it suitable for the adjoint RTE and prove its convergence. Then the RTE can be reduced to a sparse block diagonal linear system finally.
The rest of the paper is organized as follows. We introduce mathematical model of forward problem and inverse problem as well as notation conventions in Section 2. Given scattering coefficient, an improved fixed-point iterative algorithm and relevant proofs for convergence are presented in Section 3.1. The BB method for QPAT and the deduction of the gradient of objective function are detailed in Section 3.2. And we discuss solvers for original RTE and adjoint RTE respectively in Section 4.1. The simulation results are shown in Section 4.2. The conclusion comes in Section 5.

The mathematical description
In the following, we consider the convex bounded object region Ω P R n pn " 2, 3q with Lipschitz boundary BΩ and angular space S n´1 . Describing boundary condition, we divide boundary BΩ into inflow boundary and outflow boundary Γ´:" px, θq P BΩˆS n´1 : νpxq¨θ ă 0 ( , where νpxq is the outward normal vector at the position x of the boundary BΩ.

QPAT model
We use stationary RTE to describe the light propagation in the form of rθ¨∇`µ a pxq`µ s pxqs φpx, θq´µ s pxq where θ and x are the direction and the position of interest respectively, φpx, θq is the density of energy at the position x in the direction θ, µ a pxq and µ s pxq are spatially varying absorption and scattering coefficients respectively, and qpx, θq is the source term. Function f pθ, θ 1 q describing the probability of photon traveling from direction θ 1 to θ is called scattering phase function and usually is characterized by anisotropic factor g of the form which is the well-known Henyey-Greenstein (H-G) scattering function. For the sake of simplicity, we define scattering operator K by pKφqpx, θq :" We assume inflow boundary condition φpx, θq " q b px, θq, px, θq P Γ´, that is no photon getting in Ω except for boundary source q b px, θq. In QPAT, source term qpx, θq is usually regarded as zero. Absorbed energy density is presented by where Φpxq :" pAφqpxq " in which A is the accumulation operator over all directions.
Owing to absorbed energy, initial acoustic pressure distribution is generated of the form p 0 pxq :" γpxqhpxq, where γpxq is spatially varying Grüeneisen coefficient. γpxq is the efficiency of conversion of absorbed energy to pressure and is assumed to be known. Throughout this work, it is rescaled to be one, i.e. γpxq " 1.

PAT model
On account of initial pressure distribution, acoustic waves propagate through object region Ω and a series of temporal acoustic signals are detected by ultrasonic detectors at the surface BΩ. The behavior is described by where cpxq is the acoustic speed, ppx, tq is the acoustic pressure at the position x and time t, and p 0 pxq is the initial pressure distribution. The goal of PAT is to recover optical coefficients µ a and µ s given temporal data ppx, tq ppx, tq P BΩˆp0, T sq. This is often accomplished in two stages: first one recovers p 0 pxq from temporal data ppx, tq ppx, tq P BΩˆp0, T sq, and next one reconstructs µ a and µ s from initial pressure p 0 (or h) given illumination conditions. We assume the first stage has been done, and the second stage QPAT is our concern. We assume tissue is illuminated for M times with known boundary sources q bm pm " 0, 1, . . . , M´1q, then inverse problem QPAT is to determine µ a and µ s from data h m , that is hpx; µ a , µ s , q bm q pm " 0, 1, . . . , M´1q.

Notation conventions
We denote the absorption coefficient and scattering coefficient by µ i a and µ i s in the ith iteration respectively, and corresponding solution of RTE and heat with mth illumination of boundary source q bm pm " 0, 1, . . . , M´1q is denoted by φ i m px, θq and h i m pxq. The exact coefficients and data are denoted by µå , µs , φm, and hm.

Reconstruction for optical coefficients
In this section, we propose an improved fixed-point iterative method to recover absorption coefficient given scattering coefficient and then apply Barzilai-Borwein (BB) method to simultaneously recover scattering coefficient as well. First, it is necessary to make several assumptions on the two optical coefficients.

The case of given scattering coefficient
Under the Assumption 1, we provide an effective fixed-point iterative method.
Since one measurement can recover absorption coefficient, subscripts 'm' of symbols mentioned in Section 2.3 are omitted for brevity. The algorithm is detailed as follows.
Algorithm 1 Improved fixed-point iteration Input: Given initialization µ 0 a mentioned in Assumption 1, data h˚, scattering coefficient µ s , boundary source q b , and tolerance . Solve stationary RTE with absorption µ i a and scattering coefficients µ s to obtain the solution φ i . Then let h i pxq " µ i a pxqpAφ i qpxq;

3:
If h i´h˚ 1 ă , break with µ a " µ i a ; or implement next step;

6: end for
Remark. The Assumption 1 is important because φ is a denominator in the Algorithm 1. In practice, we can also use φ`ε instead with the small ε to avoid the situation where φ " 0. Fixed-point iterative method has been mentioned in many references [19,9,14,11], which is usually unstable and sensitive to initial value. However, the modification of fifth step in Algorithm 1 can guarantee convergence given scattering coefficient.
Carrying out Algorithm 1, we can obtain function sequences tĂ µ a i pxqu 8 i"0 , tµ i a pxqu 8 i"0 , th i pxqu 8 i"0 and tφ i pxqu 8 i"0 . We claim that the sequence h i´h˚ 1 converges to zero. The proof is presented as follows.
Firstly, we state a theorem on stationary RTE.
). For RTE (2) satisfying Assumption 1, if source q ě 0 and boundary source q b ě 0, there exists a unique non-negative and continuous solution φpx, θq.
In the view of physics, the intuitive explanation behind Theorem 1 is straightforward. For non-negative initial energy, the energy will not become negative while it was transferred and absorbed. Lemma 2 reveals the monotonicity of φ with respect to µ a . Lemma 2. Let φ 1 pxq and φ 2 pxq be the solutions of RTEs with µ a being µ 1 a and µ 2 a respectively. Then φ 1 px, θq ě φ 2 px, θq provided µ 1 a pxq ď µ 2 a pxqp@x P Ωq. Note that the superscripts of µ 1 a and µ 2 a are used to distinguish the different absorption coefficients.
Proof. We see that φ :" φ 1´φ2 solves the RTE px, θq P Γ´( 10) with absorption coefficient µ 2 a and scattering coefficient µ s . Since source term pµ 2 a´µ 1 a qφ 1 and boundary source are both non-negative, we can derive using Theorem 1.
Using Lemma 2, we investigate the properties of function sequences from Algorithm 1.
Then it suffices to establish (13) for pi`1q.
According to Algorithm 1, An easy induction completes the proof.
Monotonicity of µ i a can easily deduce the monotonicity of φ i .

Lemma 4.
For any px, θq P ΩˆS n´1 , the sequence tφ i px, θqu 8 i"0 obtained in the Algorithm 1 satisfies Proof. It is obvious using Lemma 2 and Lemma 3.
We analyze the update process for µ a by introducing new definitions. For any integer i, we divide region Ω into two parts in the form of # Ω ì :" tx P Ω : h˚pxq ě h i pxqu, Then we have Actually, it is a process of keeping and updating from µ i a to µ i`1 a , keeping the value of µ i a of Ω í , updating the value of Ω ì . Eventually, we prove the convergence of function sequence th i pxqu 8 i"0 . Theorem 5. For any x P Ω, the sequence th i pxqu 8 i"0 obtained in the Algorithm 1 satisfies lim Proof. We divide our proof into four claims. (i): Function sequences tµ i a pxqu 8 i"0 , tφ i px, θqu 8 i"0 and th i pxqu 8 i"0 converge pointwisely.
Since monotonicity and boundedness of tµ i a pxqu 8 i"0 and tφ i px, θqu 8 i"0 for any px, θq P ΩˆS n´1 , it follows that they converge pointwisely, whose limits are denoted by µ a pxq and φpx, θq respectively. Obviously, th i pxqu 8 i"0 converges pointwisely and its limit is denoted by hpxq. And obviously, hpxq " µ a pxqpAφqpxq.
(ii): For any x P Ω, if there is some I such that if x P Ω Ì , then x P Ω ì pi ą Iq. It is obvious that µ I`1 a pxq " Ă µ a I`1 pxq ě µ I a pxq. According to Lemma 4, we can obtain h I`1 pxq " µ I`1 a pAφ I`1 qpxq ď µ I`1 a pxqpAφ I qpxq " h˚pxq. And an easy induction shows that x P Ω ì pi ą Iq.
Thanks to x P Ω ì pi ą 0q for any x P Ω 0 , it follows that where φ 0 " sup xPΩ φ 0 pxq. According to Lemma 4, gives last inequality in (16). Then Owing to 0 ď µ a´µ 0 a ă 8, we can obtain that that is, h˚pxq´h i pxq converges to 0 pointwisely. So hpxq " h˚pxq px P Ω 0 q.
Integrating (16) with respect x over Ω 0 gives We divide region Ω 0 into two parts of form # Ω 01 : " tx P Ω 0 : Exist some i such that x P Ω ì u, According the proof of (iii) and the definition of Ω í , it can easily be seen that and From (19), it is clear that Considering RTEs with absorption coefficient µ i a and µå , it is easy to obtain Since phase function f pθ, θ 1 q satisfies Integrating (22) with respect to θ gives Integrating (23) with respect to x over Ω gives Applying Green's formula to the second term of (24), we have Apparently, we have Naturally, when i is sufficiently large, we have Hence, when i tends to infinity, we have This completes the proof of (iv). Therefore, lim iÑ8 }h˚´h i } 1 " 0.
Given scattering coefficient, by Algorithm 1, we can get monotonically increasing sequence µ i a pxq and monotonically decreasing sequence φ i px, θq, which converge pointwisely toμ a ď µ˚andφ ě φ˚respectively. Sequence h i pxq converges to exact function h˚pxq in the L 1 -norm.

Reconstruction of µ a and µ s simultaneously
In practice, the scattering coefficient is also unknown, so the improved fixed-point iteration method is not applicable for the two unknown coefficients case. It is necessary to establish a more general method to recover two coefficients simultaneously. From now on, we follow optimization approach to estimate µ a and µ s . First, we define error functional to measure the distance between measurement data hm and estimated data h m , where h m px; µ a , µ s q equals µ a pAφ m qpx; µ a , µ s q for estimated µ a and µ s . Then the reconstruction of QPAT can be reformulated as min µa,µs Fpµ a , µ s q.
Remark. We replace h m p¨; µ a , µ s q´hm 2 2 {2 by logph m p¨; µ a , µ s qq´logphmq 2 2 {2 as data-fidelity term is due to that the latter boosts the convergence of minimization method according to [27]. We discuss its advantage and properties in Appendix C.
Remark. According to inverse problem theory, regularization is useful for handling ill-conditioned problem. Usually, in image science, ill-conditionedness results in the edges of image blur and noise amplifying. Even though QPAT is a typical nonlinear problem, as pointed out in [2,22], it is relatively stable with the help of multiple measurements. We focus on the minimization of data-fidelity term for the reason that multiple measurements alleviate the ill-posedness numerically. The algorithms on the minimization of objective function with regularization term, such as some a priori information, can be easily derived by our framework. We do not explore this in this paper.
Using the definition of Fréchet derivative ∇F of F in feasible direction ph µa , h µs q P D aˆDs , we can write F 1 pµ a , µ s qph µa , h µs q :" ∇F, ph µa , h µs q L 2 pΩq .
The solution of RTE with coefficients µ a and µ s can be regarded as a function with respect to µ a and µ s , that is, φ " φpx, θ; µ a , µ s q. Fortunately, the directional derivative of φpx, θ; µ a , µ s q with respect to µ a and µ s in any feasible direction exists according to [18]. Furthermore, the directional derivative of data-fidelity with respect to µ a and µ s can be expressed analytically by µ a , φpx, θ; µ a , µ s q and its solution of adjoint RTE. This concludes in Proposition 6. Similar to the proof of Proposition 3.3 in [18], we can express the gradient of objective functional as follows.
Proposition 6. For any pairs pµ a , µ s q P D aˆDs and feasible direction ph µa , h µs q P D aˆDs , we have Notice that A˚is the adjoint operator of A and pA˚f qpx, θq " f pxq, @f P L 2 pΩq.
Proof. Following the proof in [18] (Proposition 3.3), we have (31) Notice that φ 1 m p¨; µ a , µ s qph µa , h µs q is the directional derivative of φ m p¨; µ a , µ s q with respect to pµ a , µ s q in direction ph µa , h µs q.
It completes the proof.
After deducing the gradient of error functional (26), all we need to do is to find appropriate stepsize in the negative gradient direction to decrease the functional. Many optimization methods can achieve this goal, including steepest descent, Quasi-Newton, and so on. However, these gradient-based methods usually involve the linesearch process. A linesearch step needs to solve original RTE or adjoint RTE for several times. Since solving RTE dominates the computational cost, linesearch is computationally intensive. To mitigate the heavy computational cost, the well-known BB gradient method is applied to compute stepsize. It is derived from a two-point approximation to the scant equation underlying Quasi-Newton methods [6,17]. And it is R-superlinearly convergent in the two-dimensional quadratic case [6]. Without loss of generality, we denote µ a or µ s by µ , then the update formula at kth step is Generally, there are two choices about the stepsize α k : where s k " µ k´µk´1 and y k " ∇F k´∇ F k´1 , and F k is the iterative sequence of error functional F, see (26). BB algorithm is detailed in Algorithm 2.

3:
Solve stationary RTEs (2) with absorption and scattering coefficients µ i a and µ i s respectively to obtain the solution φ i m . Then let h i m pxq " µ i a pxqpAφ i m qpxq;

4:
If F i ă 1 , end up with µ a " µ i a and µ s " µ i s ;

6:
Calculate the gradient of F i : ∇ µa F i and ∇ µs F i with respect to µ a and µ s . If ∇ µa F i ď 2 and/or ∇ µs F i ď 3 , let flag a " 0 and/or flag s " 0;

7:
If i " 0 or 1, updating µ a and/or µ s in negative gradient direction with small step size such that F i decreases; otherwise, if flag a " 1 and/or flag s " 1, update µ a and/or µ s by BB stepsize; 8: end for

Numerical simulations
The reconstructions of absorption and scattering coefficients are investigated with simulations in two cases: the reconstruction of µ a given µ s and the reconstruction of µ a and µ s simultaneously.

Solver for RTE
We consider the numerical solver for RTE in 2-D. And it can be extended to 3-D with little effort. Finite element method combined with streamline diffusion modification is applied to solve stationary RTE [30,18,25], where P 1 Lagrangian elements in spatial and angular space are used. In this way, a large sparse linear system needs to be solved which is still difficult. By improving the algorithm proposed in [15], we can solve original RTE as well as adjoint RTE on 2D and 3D unstructured mesh by Discontinuous Galerkin (DG) method combined with multigrid method, which reduces the problem to a sparse block diagonal linear system.
We divide angular space into P equal intervals, and the corresponding directions and angles are denoted by θ 0 , θ 1 , . . . , θ P´1 and β 0 , β 1 , . . . , β P´1 with interval ∆β. The spatial domain is discretized into unstructured triangular mesh. Suppose spatial domain is divided into N triangles and each triangle contains n d nodes, where n d " 3 for 2-D spatial domain. Lagrangian elements and piecewise linear DG elements are used to discretize RTE, that is, numerical discrete scheme of RTE is where L k pθq is piecewise linear basis function in angular space which takes value of 1 in direction θ k and 0 in other directions, ϕ ij is the linear basis function in spatial domain in the direction θ k which takes value of 1 in the jth nodes of ith triangle and 0 in other nodes and triangles, and φ i,j,k is the value of φpx, θq in direction θ k in jth node of ith element. Such spatial basis function can approximate discontinuous solution which is more suitable for actual situation, such as some edges of inclusions in object region. On account of scattering term in RTE, angular Gauss-Seidel iteration is applied in [15]. It iteratively solves RTE in every fixed direction in the form of where ω kk 1 " with ω 0 kk 1 " f pθ k , θ k 1 q. Obviously, equation (33) can be solved by a lot of numerical methods, and we use DG method.
A multigrid scheme is applied to solve (33), and it is reduced to solve a sparse block diagonal system. Multiplying test function ϕ ij in both sides of (33), integrating it over ith triangle τ i with respect to x, we havé from Green's formula, where subscripts i and j of φ i,j,k are omitted. In (34),φ is the value of neighboring element in inflow direction, which is the product of upwind scheme. The specific explanation of why the value inflow direction is used is as follows.
For an intuitive explanation, we use rectangular mesh to illustrate in Figure 1. For rectangular mesh, in direction θ k and node x, RTE can be discretized into where φ k :" φpx, θ k q. Backward difference is used to approximate one order derivative. Then when 0 ď β k ă π 2 , its upwind scheme is Obviously, the scheme converges because the equation (36) is diagonally dominant, and Gauss-Seidel scheme will accelerate the convergence. In fact, from (36), φ i,j,k is updated using Apparently, φ 0,j,k and φ i,0,k are known as inflow boundary condition. Therefore, (37) is an explicit scheme which updates φ i,j,k , φ i,1,k , φ 2,j,k , φ i,2,k ,¨¨¨pi, j, " 1, 2,¨¨¨q successively. We can see that (36) and (34) both use inflow information to update outflow information. This may be a kind of intuitive explanation. Indeed, the convergence of (34) for vacuum boundary condition is proved in [16], see Appendix A. However, the above numerical scheme fails to solve the adjoint RTE with same updating order, that is to use inflow information to update outflow information. We provide a reverse updating order for elements to guarantee the convergence of solver for adjoint RTE. The solver for adjoint RTE is detailed as follows.
For adjoint RTE (30), the corresponding discrete scheme in direction θ k for rectangular mesh iś Differently, forward difference is used to approximate the one order derivative, and when 0 ď β k ă π 2 , its upwind scheme is pa`b`µ a`µs´µs ω kk qφ i,j,k´p aφ i`1,j,k`b φ i,j`1,k q´µ s ÿ k 1 ‰k Then φ i,j,k is updated by The use of forward difference makes (39) diagonally dominant again. Besides, outflow information is used to update inflow information in (40). Heuristically, for unstructured mesh, we have whereφ is the value of neighboring element in outflow direction which is used to update inflow information. Using similar discussion, the convergence of (41) is proved in Appendix B. Furthermore, we apply multigrid method to solve RTE to accelerate convergence.
From above discussion, no matter which mesh is used, the key of convergence is to update outflow information using inflow information for the original RTE (2) and update the latter using the former for the adjoint RTE (41). The convergence in rectangular mesh is obvious. As for triangular mesh, the update scheme results from variation analysis, so the convergence proofs are obtained by discussing corresponding variations, see Appendix A and B.
Therefore, for original RTE (2), there are three layers of loops. We apply the multigrid scheme at the outermost loop. Since it involves two coordinate systems, i.e., angular-coordinate and spatial-coordinate, there are variable updating schemes in terms of the iteration order, such as the angle-prior and space-prior. The second-layer loop is about direction, that is, the corresponding spatial equation (33) is solved in turn for each direction θ k pk " 0, 1, . . . , P´1q. In third-layer loop, for each element in the direction θ k , the value of φ i,j,k is updated iteratively through solving a 3ˆ3 linear system (34). Note that the updating order is in the following order: consider the each element successively from the boundary along the direction of θ k . We refer the interested reader to [15] for its algorithm details and to [16] for its theory. On the contrast, for adjoint RTE, in third-layer loop, inspired by (39) we propose to consider each element successively from the boundary along the direction of´θ k , which just reverse the updating order of original RTE. And linear system (41) needs to be solved. The corresponding pseudo-code is presented in Algorithm 3.

Numerical results
In this subsection, we present some numerical results to demonstrate the numerical performance of improved fixed-point iteration and BB method. For the sake of simplicity, only 2D reconstruction is investigated. The anisotropic factor g equals 0.9. Applying RTE solver described in 4.1, we can get discrete energy density H of for each direction θ k pk " 0, 1, . . . , P´1q do 3: for each element τ i do 4: Updating φ i,j,k pj " 0, 1, 2q by solving a 3ˆ3 linear system (34) for original RTE or (41) for adjoint RTE. For original RTE, the updating order is from boundary along θ k through each elements to other side of region. For adjoint RTE, the updating order is reverse.
These regions are depicted in Figure 2, where the optical coefficients of 1st, 3rd, and 4st templates are piecewise constant, and the second one is smooth in object region Ω 0 . These templates contain discontinuous borders as well as continuous borders, whose corners are both sharp and smooth.  4.2.1. Reconstruction of µ a given µ s Using one measurement with boundary source in (20,0), we apply improved fixed-point iterative method and BB method to retrieve absorption coefficients of templates given scattering coefficients. Our initial guesses are set to be the same as background. Figure 3 shows results after 50 iterations for noiseless data. Then 5% Gaussian noise (i.e. " 5%) is added to data according to (42), and the final reconstructed images are showed in Figure 4. Specific iterative relative errors are showed in Figure 5. From Figure 3 and Figure 4, it seems that both methods can retrieve absorption coefficient with almost the same accurate solutions in the case of noise-free and noisy measurements. However, from Figure 5, it is obvious that improved fixed-point iteration converges more rapidly than BB method. Improved fixed-point iteration achieve the critical point only after a few iterations. To attain the same accuracy, the number of iteration of improved fixed-point iteration is about the half of the one of BB method. Owing to solving the adjoint RTEs in BB method, improved fixed-point iteration can achieve the accuracy as the BB method with about quarter computational time of the latter.

4.2.2.
Reconstruction of µ a and µ s simultaneously We apply BB method to reconstruct absorption and scattering coefficients simultaneously from four measurements. The boundary point sources are placed in the top, bottom, left and right sides respectively. The reconstruction results from noiseless data are showed in Figure 6. Corresponding reconstruction results from data added by 5% Gaussian noise are showed in Figure 7. The relative L 2 errors of reconstructions of optical coefficients of four templates are tabulated in Table 1. From Figure 6 and Figure 7, we can see that BB method retrieves absorption coefficient accurately in piecewise constant and smooth cases, and even is stable to noisy data. In the first, third and fourth templates, some jagged borders of inclusions are reasonable, which are caused by the interpolation of data from fine mesh to coarse mesh. Nevertheless, the reconstruction of scattering coefficient is not satisfactory. In piecewise constant case, the edges between different pieces are blur, which can be seen from the first, third, and fourth template in Figure 6 and Figure 7. This is due to the insensitive scattering coefficient. After all, there is no explicit µ s in the formula (6).

Conclusion
In this paper, we investigate the reconstruction of absorption and scattering coefficients in QPAT in two cases. Given scattering coefficient, we propose an improved fixed-point iterative method to reconstruct the absorption coefficient and prove its convergence. The advantage of this reconstruction algorithm is that its fast convergence and it does not require the initial guess close to the exact solution.
Meanwhile, it does not need to solve adjoint RTE in the optimization approach.
For the simultaneous reconstruction of the two coefficients, we apply a state of art BB method, which dose not need linesearch to compute stepsize. Indeed, linesearch involves solving RTE for several times, which is expensive. BB method only use the values of previous two steps to update current estimation. Numerical results show that improved fixed-point iteration can achieve almost the same accuracy with BB method. It takes about quarter the computational time of the BB method. Moreover, the algorithm is stable to noise. In the unknown scattering coefficient case, numerical results show that BB method can obtain quite accurate absorption estimate. However, the result of the estimation of scattering coefficient is not satisfactory. There is no explicit µ s in (6), so it is insensitive to the objective function. There is obvious blur near the borders in piecewise constant templates. In future, we will design a better error functional to improve the reconstruction of µ s . Although BB method can recover numerically absorption and scattering coefficients simultaneously, its convergence remains to be studied. Besides, for more practical application, we will do some realistic simulations, such as 3-D templates or real physical QPAT data.
The convergence of spatial discrete scheme (34) is discussed next. Assume triangulation T h with h " sup SPT h diampSq, V d h :" tv : v| SPT h P P d pSq, v| ΩzS " 0u, where P d pSq is the space of d-degree polynomials. Let pu, vq S :" ş S uv dx, Γ h :" pY SPT h BSqzBΩ and BΩ ḱ p`q :" tx P BΩ : ν¨θ k ď pąq0u. Then the convergence result is as follows.
From Theorem 11 and Theorem 12, we see if sup f pxq ă 1, taking log not only makes the error functional steeper near the minimizer for some fixed direction, but also the farther between x and x˚, the change from J 1 pxq to J 2 pxq is more significant. Then, for multivariate function, taking log makes the contours not so narrow. Thence, log-type error functional improves the convergence of minimization method. As for the general case of sup f pxq ą 1, we can replace f pxq by f pxq{c for some constant c ą sup f pxq.