Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization

: In this paper we study an l1-regularized multilevel approach for bioluminescence tomography based on radiative transfer equation with the emphasis on improving imaging resolution and reducing computational time. Simulations are performed to validate that our algorithms are potential for efficient high-resolution imaging. Besides, we study and compare reconstructions with boundary angular-averaged data, boundary angular-resolved data and internal angular-averaged data respectively.


Introduction
Bioluminescence tomography (BLT) [1][2][3] is an emerging molecular imaging modality, which can be used to monitor physiological and pathological activities at molecular levels, specially visualize primary tumor growth and tumor cell metastasis. In BLT [4][5][6][7][8][9][10][11][12][13], one tries to recover the location and intensity of isotropic bioluminescent sources from boundary measurements. However, it remains challenging in this area that it is highly desirable but extremely difficult to improve the image resolution as well as reconstruction stability and efficiency in biomedical and clinical practice.
Regarding the forward model of in vivo photon migration, radiative transfer equation (RTE) [14][15][16][17] is a well-accepted golden standard and needs to be appropriately discretized numerically to handle complex geometry [18]. In practice, due to its computational burden as a multi-dimensional problem involving both angular directions and spatial variables, the most popular forward model is diffusion approximation (DA) [19], a lower-order approximation of RTE. However, the validation of DA breaks down in small animal imaging, for example, in the region close to light sources or with relatively low scattering. Recently, we have developed a fast forward solver to reduce the computational complexity so that numerical RTE solver is accessible in practice. We will present numerical examples to show that BLT based on RTE can achieve satisfactory resolution, which is usually not the case for BLT based on DA.
In this work, we assume that the optical property of the underlying medium is known. So BLT becomes a linear inverse source problem. In practice the main challenge is due to severe ill-posedness and large degrees of freedom for inverse problem. The ill-posedness may come from two sources: (1) insufficient measurements, such as angular-averaged data, that cause an underdetermined system; (2) ill-conditioning of the system matrix whose columns correspond to Green's functions of RTE. The ill-posedness makes inversion or reconstruction sensitive to small perturbation in data, such as measurement error and noise. Moreover, the computation cost for computing the Green's function for RTE increases dramatically and the inversion becomes even more ill-posed when degrees of freedom are increased. This type of problem is typical in many inverse problems. Construction of efficient, accurate and stable reconstruction algorithms is problem-dependent and challenging.
In our application we will explore sparsity structure in the inverse problem. For many BLT applications, the bioluminescent source is sparsely distributed in space. We propose a multilevel imaging method to fully take advantage of the sparsity to improve both efficiency and stability. In this multilevel framework, we first reconstruct the sparse source efficiently on a coarse mesh with good localization, which provides a good initial guess as well as a small region of interest (ROI) for the fine reconstruction. And next reconstruction on fine mesh is carried out within the localized ROI with dramatically reduced degrees of freedom. Through such coarse-to-fine reconstructions, one can stably extract high-resolution images while substantially reducing the computational time. As another consequence, the amount of measurements can also be reduced by utilizing source sparsity in our multilevel reconstruction.
The first crucial step in our multilevel approach is efficient and stable reconstruction on the coarse mesh with good localization property, which is equivalent to finding the sparse solution to the linear system defined on the coarse mesh. How to find the sparse solution efficiently, accurately and stably with as few measurements as possible is at the heart of compressive sensing methods. Compressive sensing has been studied for decades and recently has received a lot of attention due to theoretical justification and many practical implementations, see [20-25] and references therein. The key idea is that, given the sparsity of the signal, an efficient and stable algorithm through l1 regularization can be applied to recover or approximate the original signal with appropriate measurements that have cardinality comparable to the sparsity of the signal. In statistics community, a large body of recent work on least absolute shrinkage and selection operator (LASSO) is also on l1 regularization for sparse model selection, see [26,27] and references therein. However, in order to recover the sparse signal accurately and stably, current theory requires some nice property of the underdetermined system or the measurement matrix. Essentially it requires near-orthogonality or incoherence among columns of the system matrix. For example, the commonly used restricted isometry property (RIP) requires well-conditioning for all submatrices of size comparable to the sparsity or near-orthogonality among columns. RIP is used for both recoverability and stability results in [22][23][24][25]. Although a non-RIP result is established in [28] by using characterization of null space of the measurement matrix, the stability results still depend on the conditioning of the system. In recent work on LASSO from [26], an irrepresentable condition is proposed for consistent model selection by LASSO with the right amount of l1 regularization. However, the irrepresentable condition is again related to angle condition between columns of the measurement matrix that are used in the true model and the rest of the columns that are irrelevant. These requirements cannot be satisfied by our problem as well as most inverse problems due to the ill-posedness. In this work we show that by combining l1 regularization with the multilevel approach we can still take advantage of the sparsity and improve computation efficiency, stability and reconstruction results while requiring fewer measurements.
The inverse source problem for RTE has been studied in several papers in mathematics and engineering communities [29][30][31][32][33][34][35][36][37][38]. Here we approach this inverse source problem directly by taking advantage of the source sparsity with both l1 regularization and multilevel reconstruction. As we finish up this paper, we recognize a similar l1-regularized strategy for BLT was studied recently in [39], in which BLT based on DA with regularization by differentiable approximation of l1 norm instead of l1 norm was studied due to the difficulty of dealing with non-differentiable l1 norm. Also adaptive mesh rather than multi-level mesh has been studied in the context of adaptive finite element method for BLT [57]. However, our l1regularized multilevel approach is entirely different although we both use "multilevel".
Here is the outline of this paper. In section 2, we shall first introduce notations and discretization methods for RTE forward solver, formulate BLT as a linear inverse source problem, describe l1 regularization strategy, and then summarize our algorithms in the multilevel framework; in section 3, we shall validate our algorithms with simulations; in section 4, we shall conclude and comment on our future work.

Forward problem
In this paper, we model light propagation by RTE with the vacuum boundary condition where the quantities in the equation are the photon flux Ψ which depends on both space r and angleŝ , the bioluminescent source q , the Henyey-Greenstein scattering function f with anisotropic factor g , the absorption coefficient a µ , the scattering coefficient s µ , the outer normaln at boundary ∂Ω of domain Ω and the angular space 1 n S − , i.e., the unit circle in 2D or the unit sphere in 3D. The methods and algorithms in this paper can be extended easily to RTE with reflection boundary condition. We shall only describe RTE with vacuum boundary condition for the presentation purpose.
An efficient solver of (1) was proposed in [18] with angular discretization by finite element method and spatial discretization by piecewise-linear Discontinuous Galerkin method. We design multigrid methods in both angular and spatial space with improved source iteration as relaxations to achieve fast convergence in various regimes. We shall briefly introduce notations and formulate forward problem as follow.
Note in BLT, we assume the unknown bioluminescent source q is isotropic, i.e., Then we discretize (1) into the following matrix form . For details, we refer readers to [18]. Note that in practice we do not explicitly form matrices in (2) since our method is local in the sense that (2)

RTE-based BLT
In this section, we describe RTE-based BLT through Green's functions.
Since the light source is assumed to be isotropic, let( , ; ) G r s r′ be the Green's function of Mathematically, we need enough angular-resolved data in (3) to match the dimension of unknowns. However, what is usually available in practice are boundary angular-averaged data, for which we have Due to dimension mismatch between unknowns and measurements in (4), no unique inversion can be defined mathematically. Numerically, insufficient data lead to a underdetermined linear system, especially on the fine mesh. In this work we apply the proposed multilevel reconstruction algorithm to both cases of (3) and (4) and show that computation efficiency, image resolution and reconstruction stability can be achieved when the source is sparse. In particular by exploring sparsity with multilevel approach one can alleviate the ill-posedness, improve the efficiency and reduce the amount of measurements. Now we denote measurements by X , Green's functions by G and write (3) and (4) in the following matrix form for BLT GQ X = .
(5) In particular, the i th column of G corresponds to measurements from the unit source sitting on i th element and the j th row of G corresponds to adjoint solutions of RTE from j th measurement operator.
The main difficulty to invert the linear system (5) with boundary data is due to its illposedness. First (5) is underdetermined with insufficient data, i.e., s M N << , especially in fine reconstructions with limited angular-averaged data. Second, matrix G is ill-conditioned, i.e., columns of G for adjacent sources are almost parallel since adjacent sources, especially those in depth, produce the similar measurements at the boundary. This can also be seen through resolution analysis [32,40]. Moreover, G is even more ill-conditioned on finer mesh.
As a result, in BLT, the reconstructed image is usually blurred with low resolution, and hence does not provide good localization of the source; lots of measurements are essential for a stable reconstruction, which requires costly experiments and intensive computation; the reconstruction is usually sensitive to noise.
To deal with the ill-posedness, proper regularization is crucial, but subtle for inverse problem. A good choice of regularization should take into account a priori knowledge of the inverse problem. Based on sparsity assumption in our case, we choose l1 regularization and combine with a multilevel approach to explore source sparsity. For comparison reason we first present the conventional l2 regularization.
Before doing that, since the analytic Green's function is only available for some simple geometry and it has to be computed numerically in most cases, first we shall start from (2) and derive discrete version of (5) in the next subsection.

Computation of numerical Jacobian
. Since the solution Φ linearly depends on Q , we have , which is so called Jacobian or sensitivity matrix. Next we shall compute Jacobian.
Differentiating both sides of (2) with respect to i q , we have , the adjoint solution of RTE with j P as the source. Thus for adjoint method by (7), we need to compute M forward problems.
In the conventional setting for BLT with boundary angular-averaged data, M is usually much smaller than s N and the adjoint method is preferred in this case. However, there are two cases that s N is smaller or much smaller than M . First is when a large amount of data is available, such as boundary angular-resolved data or internal data. Second occurs in our multilevel approach for sparse source reconstruction since multilevel method allows us to reduce the number of unknowns dramatically on the fine mesh through good localization by coarse reconstruction so that the direct method may be favorable, and thus greatly reduces the number of forward problem computation Regarding the computational time, the most time-consuming part is to compute the Jacobian J . However, it can be done in parallel, which we will implement in the future. On the other hand, it is obvious that the computation cost of (6) is proportional to the number of measurements. So it is highly desirable to have an accurate and stable reconstruction with as few measurements as possible.

l2 regularization
The popular way for BLT is to formulate it as a least-square minimization problem with an extra l2 regularization term, i.e., where 2 λ is the parameter to balance the data fitting and l2 regularization.
In general, a good choice of 2 λ is not known as a priori. Therefore, we choose a nonlinear optimization method through Levenberg-Marquardt method [41,42] to automatically adjust 2 λ adaptively, in which (8) is reduced to a few iterations Note that l2 regularization is embedded in the algorithm with adaptive 2,k Although (8) is relatively easy to solve, it does not take advantage of sparsity of source distribution. As a matter of fact the reconstructed image from BLT is usually blurred with low resolution, and does not localize well for sparse source. With insufficient data and noise, the method becomes more problematic, which we will see in simulations.

l1 regularization
When the source in BLT is sparse, l1 regularization is a natural and viable choice for finding a good approximate sparse solution. More importantly, finding a good sparse solution on a coarse mesh means good localization which is crucial in identifying the region of interest on finer mesh in our multilevel approach. Although the beautiful compressive sensing theory, which says that with both l1 regularization and appropriate measurements whose cardinality is comparable to the sparsity of the unknown one can stably recover the sparse signal, does not apply to our inverse problem due to severe ill-conditioning of the system matrix, we demonstrate that by combining l1 regularization with the multilevel approach we can still take advantage of the sparsity and improve computation efficiency, stability and reconstruction results while requiring fewer measurements.
Therefore, we regularize l1 norm instead of l2 norm as follow where 1 1 || || : is a spatially weighted l1 norm with area/volume i A and takes into account of mesh nonuniformity. The simulations in next sections show the improvement of l1 regularization (10) over l2 regularization (8) especially in terms of localization and stability when the source is sparse and the amount of measurements is small.
The numerical algorithm for problem (10) is not as easy as for problem (8) since l1 norm is non-differentiable. Since the main purpose of this paper is not about efficient solver of (10), we use a standard interior-point method [44-46] for minimizing (10) and point out that there are many recent works [47-52] in developing efficient algorithms for (10). All these methods can be directly used in our algorithms to improve the optimization efficiency further. Based on the interior-point method we choose for (10), around eighty to ninety percents of the computational time are spent on computing Jacobian J .
In this paper, we solve the convex optimization problem (10) by the barrier method [44]. We shall briefly formulate the method as follow.
First, we replace the l1 norm by inequality constraints as with constraints , 1, , , and then add the constraints into the minimization problem by logarithmic barrier penalties as Next we solve a sequence of t -subproblems (12) with increasing t , where the solution t Q of (12) is no more than 2 / s N t -suboptimal which implies t Q converges to the optimal point of (11) as t → ∞ .
In general, the computational time for (11) is longer than (9). To reduce the optimization time by (11), a preconditioned conjugate gradient method [45,46,53] is implemented to compute the search direction for backtracking line search.

Multilevel reconstruction
It is typical for an inverse problem to become significantly more difficult to solve as mesh is refined. This is because not only degrees of freedom are increased as a more underdetermined linear system but also more importantly the system becomes much more ill-conditioned. Even one can afford significant amount of measurements and computation cost, it is often that using a fine mesh does not necessary lead to a high-resolution reconstruction in practice. The key points for our multilevel approach for sparse source reconstruction for BLT are as follow.
First, the inverse problem on the coarse mesh can be solved more efficiently and stably. Second, the solution from the coarse mesh provides not only a good initial guess but more importantly a good localization for the sparse source with l1 regularization. The localization will allow us to solve an inverse problem restricted to a small region on fine mesh with the significantly reduced degrees of freedom and improved stability.
Last, the amount of measurements for this multilevel approach only needs to provide a good localization on coarse mesh and matches the degree of freedom in the confined region on fine mesh, which is determined by both sparsity of the source and resolution of the fine mesh.
In our multilevel method, we assume the spatial meshes can be refined from one to another. First we describe a two-level algorithm including three steps: (1) minimize (10) on the coarse mesh : (2) identify the support c S of c Q on the coarse mesh by simple threshold of c Q , e.g., with respect to the maximal values, i.e., c , and interpolate the support c S from c T onto the fine mesh ≤ ≤ on the fine mesh; (3) minimize (10) within the support f S on f T to obtain For sparse source reconstruction, (10) is preferred and usually we have f . In support detection, ε is between 0.01 and 0.05, and we actually extend the support a little bit by one or two grids according to the distance function between centers of grids.
Algorithm: multilevel approach for BLT otherwise, find i S by threshold of 1 i Q − , and interpolate from Step 2: compute Jacobian with either (6) or (7).
Step 3: reconstruct i Q on i T by solving (8) or (10).

END
Next, we describe our multilevel algorithm in detail, which also incorporates the flexibility from the dual-mesh scheme [54], i.e., the mesh for unknowns is not necessarily the same as the mesh for solving RTE in the reconstruction. The advantage of the dual-mesh scheme is that we not only improve the reconstruction stability by doing coarse reconstruction and also keep the accuracy of forward solver by solving it on fine mesh.
We start to formulate our multilevel method by defining a sequence of meshes {( , , ), }

Simulations
In this section we use various numerical simulations to validate our l1-regularized multilevel imaging algorithm for BLT. Most simulations are performed in two dimensions (2D), which serves the purposes of this paper. Although we give an example (Fig. 11) in three dimensions (3D), we will study 3D in vivo BLT further in more practical setting in our future study. In all figures shown in this section, reconstructed images are plotted as what they originally are from reconstructions without further image processing, such as thresholding. As pointed out at the beginning, radiative transport equation (RTE) is used as the underlying model for both forward and inverse problem. Data for forward problem is generated on a fine mesh that is different from the fine mesh used for inverse problem. The fast forward solver developed in [18] is used to compute both the forward problem and the Green's functions for inverse problem. In particular we design different simulations to demonstrate the following points.
First, for sparse source, l1 regularization localizes the source better than l2 regularization, especially when there is noise and small amount of data. See Fig. 1, 2, 4 and 5.
Second, our multilevel algorithm is more effective and efficient than solving the inverse problem directly on the fine mesh. See Fig. 6.
Third, our l1-regularized multilevel algorithm works well with both angular-resolved data and angular-averaged data. Besides, reconstruction with angular-resolved data has the potential for high-resolution imaging. See Fig. 4 and 5.
Fourth, RTE based model can reconstruct better than DA based model with the same amount of angular-averaged data. See Fig. 7 and 8.
Last, when interior measurements are available, better and more stable reconstruction can be achieved. See Fig. 9 and 10.
The following 2D simulations are performed on a square domain with 20mm side length with homogeneous optical properties In the first example (Fig. 1), a small inclusion (c) with 1mm diameter is centered at (5mm, 5mm). The simulated angular-averaged measurements with 12 M = is generated with 64 angular directions and a spatial mesh (a) of 697 nodes and 1312 elements, while the reconstruction (b) mesh is of 32 directions, 697 nodes and 1312 elements.
The image from l2 reconstruction without noise is displayed in (d) in Fig. 1. While the images from l1 reconstruction without or with noise are displayed in Fig. 2  , where σ is the signal to noise ratio, and N is a Gaussian random variable with unity mean and unity variation. In Fig. 2, we vary σ from 0 percent to 50 percents with 10 percents increment and plot l1-regularized image. Comparing images from l2 and l1 regularization, we conclude that l1 regularization can efficiently localize the sparse source with a few measurements, even in the presence of large noise, while l2 regularization fails to reconstruct the sparse source.     Figures (a), (c) and (e) are reconstructed images by l1 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively. Figures  (b), (d) and (f) are reconstructed images by l2 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively.
In the second example (Fig. 3), three letter inclusions (a) are put in the domain to evaluate our reconstruction algorithms; the forward mesh (b) has 1509 nodes and 2896 elements; the coarse mesh (c) for reconstructions has 365 nodes and 668 elements and the fine one (d) has 1397 nodes and 2672 elements.
In Fig. 4 and Fig. 5, the reconstructions with anisotropic factor g = 0.9 and 0.0 respectively are carried out to demonstrate that l1 regularization ((a), (c) and (e) in Fig. 4 and Fig. 0.5) in general reconstructs better images than l2 regularization ((b), (d) and (f) in Fig. 4 and Fig. 0.5), especially when inverse problem is more severely ill-posed due to fewer available data, e.g., angular-resolved data.
Comparing (a) and (c), we find that angular-resolved data are preferred over angularaveraged data for the similar amount of data, especially in the forward-peaking regime with big anisotropic factor g since anisotropic angular information improves the probing ability. When a rich amount of angular-resolved data is available as in (e), our algorithm with l1 regularization is able to reconstruct high-resolution images.
A justification of this is that boundary angular-resolved data are essentially threedimensional in 2D, which is enough to match the dimensions of unknown sources. As a result, the reconstruction based on RTE has a great potential for high-resolution images even with only boundary data, while one does not expect a good image from DA-based reconstruction due to the dimension mismatch between the boundary data and unknowns since DA does not contain angular information.
In conclusion, one needs RTE model to alleviate the ill-posedness of standard BLT when only boundary data are available; on the other hand, efficient algorithms and proper regularization, e.g., l1 regularization for sparse source, are essential to recover high-resolution images. Besides, one may already notice a noise dot in the middle of reconstructed images in (a), (c) and (e) in Fig. 5. due to the ill-posedness of inverse problem. However, this dot has small contrast on the coarse mesh, thus can be removed through our multilevel approach by choosing proper threshold since its value is much smaller than that of other sources.  Fig. 6, we conclude that, with smaller amount of data, our multilevel approach gives better reconstruction than the single-mesh one. On the other hand, it implies that the multilevel approach is able to reduce the computational time for a given accuracy. Moreover, we may further accelerate the reconstruction by solving much less number of forward problems with the direct method (6) for Jacobian on the fine mesh when the sources are sparsely localized.  In Fig. 7 and 8, we compare the RTE-based reconstructions with DA-based reconstructions in diffusive or non-diffusive regime, for which l1-regularized reconstructions based on either RTE or DA are performed with 120 boundary angular-averaged measurements. To avoid the biased model error, for DA-based reconstructions, the simulation measurements are generated with DA as well.
In Fig. 7, we use the same meshes ((b) and (d) in Fig. 3) respectively for data generation and inverse problem for DA. In DA, we have the reduced scattering coefficient with g = 0.0 and 0.9 respectively. From reconstruction results, we conclude that BLT with DA does not produce as good images as with RTE. Actually with anisotropic g or within the domain of a few mean free paths, BLT with RTE produces much better results than with DA. Besides, RTE may utilize angular-resolved data for high-resolution images while DA is unable to incorporate the angular information.
Clearly the reconstruction is in non-diffusive regime containing both strong absorbing region and void-like region. From reconstruction results as in Fig. 8, we conclude that RTEbased BLT performs better DA-based BLT in the presence of non-diffusive regions. On the other hand, reconstruction with l1 norm tends to give sparse solution, which may have smaller support than the original support as shown in (d). In our next paper, we shall introduce the methods to resolve this issue.
As the last simulations in 2D, we compute two examples in Fig. 9 and 10 when internal angular-averaged data are available. Fig. 9. Comparison of l1-regularized reconstructions with boundary angular-averaged data and internal angular-averaged data when g = 0.9. Figure (a) is the true object. Figure (b) is the mesh to generate the forward data. Figure (c) is the result with 120 boundary angular-averaged data. Figure (d) is the result with 120 random internal angular-averaged data.
In Fig. 9, an L-shaped source inclusion (a) is embedded at the center of the square domain with the forward mesh (b) of 1469 nodes and 2816 elements. The fine mesh in the second example above ((d) in Fig. 3) is used for reconstruction. We randomly pick up 120 internal angular-averaged measurements and plot the reconstruction result in (d) while (c) is with 120 boundary angular-averaged measurements. Clearly, internal angular-averaged data give better reconstruction image than boundary angular-averaged data in reconstructing sources in depth. Fig. 10. l1-regularized reconstructions with internal angular-averaged data when g = 0.9. Figures (a) and (b) are reconstructed images through the multilevel approach with 668 measurements on the coarse mesh and 324 measurements on the fine mesh respectively.
In Fig. 10, we perform the multilevel method when full internal angular-averaged data are available. However, only partial data are necessary on the fine reconstruction. Again it indeed confirms the superiority of l1-regularized multilevel reconstruction: we first efficiently localize the support of unknowns through l1 minimization on coarse mesh, and then resolve the high-resolution image on fine mesh by optimization only in the region of interest while computing much less number of forward problems than the single-mesh reconstructions.
In conclusion, the internal data will better-pose the inverse problem, which motivates us to explore in the future the potential imaging modality to generate internal data, such as one similar to photoacoustic tomography [55]. ; two isotropic sources are embedded into two cubes with side length 2mm centered at (4mm, 18mm, 7mm) and (10mm, 4mm, 13mm) respectively. The forward mesh for data generation has 258 angular directions and a spatial mesh (b) of 3500 nodes and 17683 elements; the reconstruction mesh has 66 angular directions and a spatial mesh (c) of 2837 nodes and 14592 elements; 128 boundary angularaveraged measurements are used for the reconstruction, and thus 128 adjoint RTE solvers are computed in the reconstruction, which takes roughly three hours in total on a desktop with Intel CPU E6850 (3.00GHz). Comparing true source (b) and reconstructed source (c), we conclude that our algorithm is able to reconstruct the source in a 3D non-diffusive regime, e.g., high absorbing and void-like region. As an extension of the proposed methods in this paper, we shall study 3D RTE-based BLT in details in the more practical setting in our future work. Besides, although our algorithm is pretty efficient even in 3D, it can be accelerated further by parallelizing RTE forward solvers in both forward and inverse problem.

Conclusions and discussions
In this paper, we first introduce l1-regularized reconstruction, which is proven to perform in general better than l2-regularized reconstruction, especially for sparse sources with less accessible data and large noise. Second, we propose a multilevel reconstruction approach, which reconstructs better and takes less computational time than the single-mesh reconstruction. With the combination of both, our RTE-based l1-regularized multilevel approach has the potential for high-resolution images with even only boundary data.
On the other hand, we compare reconstructions with boundary angular-average data, boundary angular-resolved data and internal angular-averaged data. With either rich angularresolved data or internal data, we are able to reconstruct almost the exact light sources. As a result, we shall explore on the potential imaging modalities which make those data accessible.
As a similar problem to BLT, linearized optimal tomography based on perturbation of a priori background can be turned into an inverse source problem exactly as BLT, to which the algorithms here can be directly extended.
It is also worthy of mentioning another benefit from multilevel approach: that is coarse reconstruction is able to provide a good initial guess for fine reconstruction, which is crucial for nonlinear inverse problem, such as optical tomography [56].
Last, the l1 regularization can be very efficient for sparse source reconstruction. However, in many interesting cases, the total variation instead of unknowns itself is sparse. Thus we shall consider regularizing total variation instead in our next paper.