A Two-Stage Method for Piecewise-Constant Solution for Fredholm Integral Equations of the First Kind

Abstract: A numerical method is proposed for estimating piecewise-constant solutions for Fredholm integral equations of the first kind. Two functionals, namely the weighted total variation (WTV) functional and the simplified Modica-Mortola (MM) functional, are introduced. The solution procedure consists of two stages. In the first stage, the WTV functional is minimized to obtain an approximate solution f∗ TV. In the second stage, the simplified MM functional is minimized to obtain the final result by using the damped Newton (DN) method with f∗ TV as the initial guess. The numerical implementation is given in detail, and numerical results of two examples are presented to illustrate the efficiency of the proposed approach.


Introduction
In many physical problems, the relation between the quantity observed and the quantity to be measured can be formulated as a Fredholm integral equation of the first kind: where the kernel function k and the right-hand side g are known, while f is the unknown to be determined.The Fredholm integral equation of the first kind is ill-posed; see for instance [1], Chapter 2. In practical applications, there is noise in the right-hand side; therefore, Equation (1) should be revised as: where η represents the error.The above equation can be written in the operator form: where (K f )(x) = 1 0 k(x, t) f (t)dt, 0 ≤ x ≤ 1 and h = g + η.Numerical methods for obtaining a reasonable approximate solution to the Fredholm integral equation of the first kind have attracted many researchers, and many research results have been achieved; see, for instance, ( [2], Chapter 12), and [3][4][5][6][7].Due to the ill-posedness nature of the problem, numerical solutions are extremely sensitive to perturbations caused by observation and rounding errors.Therefore, regularization is required to obtain a reasonable approximate solution.Many regularization methods have been proposed; see, for instance, [8] Chapters 2 and 8.The Tikhonov regularization method [9], the truncated singular value decomposition (TSVD) method [10], the modified TSVD Mathematics 2017, 5, 28; doi:10.3390/math5020028www.mdpi.com/journal/mathematics(MTSVD) method [11], the Chebyshev interpolation method [12], the collocation method [13,14], the projected Tikhonov regularization method [15], and so on, are applied to obtain approximate continuous solutions of Equation (2).The total variation (TV) regularization method [16][17][18][19][20][21][22], adaptive TV methods [23][24][25][26][27], the piecewise-polynomial TSVD (PP-TSVD) method [28], and so on, are applied to obtain approximate piecewise-continuous solutions.
In this paper, we focus on the case where the solution of Equation ( 1) is piecewise-constant and the possible function values are known, that is: where c 1 < c 2 < • • • < c m are given constants.This kind of problem arises in many applications, for instance barcode reading, where m = 2, c 1 = 0, c 2 = 1, and image restoration, where m = 256, c p = p − 1, p = 1, 2, . . ., 256.
This paper is organized as follows.In Section 2, two objective functionals are introduced; one is based on a weighted TV functional to allow the TV regularization to be spatially varying, and the other one is based on a simplified Modica-Mortola (MM) functional to make use of the a priori knowledge of the solution.In Section 3, the implementation of numerical methods for solving Equation ( 2) is presented in detail.In Section 4, numerical examples are presented to illustrate the effectiveness of the proposed approach.Finally, concluding remarks are given in Section 5.

The Objective Functionals
The TV regularization method has been shown to be an effective way to estimate piecewise-constant solutions.This method looks for a numerical solution that has small TV, which is not inclined to a continuous or discontinuous solution.The weighted TV regularization scheme is more efficient [23], which allows the TV regularization to be spatially varying, that is a small weight is used if there is a possible edge and a large weight if there is no edge.
A commonly-used weighted TV functional is defined by: where β > 0 is a parameter and ω(x) ≥ 0 is a weighting function.An example of ω(x) is: where θ > 0 is a parameter; see [23].We can see that 0 < ω(x) ≤ 1, and ω(x) decreases as | is large, the corresponding weight is small.The above weighting function is unsmooth, so we modify it as: One can use the weighted TV functional for Equation (2) to obtain a piecewise-constant solution.In this paper, we consider solving the following minimization problem: where: and α > 0 is a small regularization parameter.
To impose the constraint ∏ m p=1 ( f (x) − c p ) = 0 to the numerical solutions, we consider the following simplified MM functional: where γ > 0 is a regularization parameter and: Note that the MM functional, which is given by: (see for instance [29]) is a useful tool for solving constant constraint problems, e.g., Bogosel and Oudet applied the functional to analyze a spectral problem with the perimeter constraint [30].
It must be pointed out that the MM functional M( f ) is not convex and has many minimal points.In general, we can only obtain a local minimizer for Φ M ( f ).In other words, the numerical solution obtained by minimizing Φ M ( f ) depends on the initial guess.This observation motivates us to consider estimating piecewise-constant solutions for (2) by a two-stage method: compute the minimizer f * TV of Φ TV ( f ) and then use f * TV as the initial guess to obtain a minimizer of Φ M ( f ) to obtain the final result.

Discretization
To solve the minimization problems ( 5) and ( 6) numerically, we need to discretize the relevant functionals.We use the midpoint quadrature and the central divided difference to discretize integrals and first order derivatives, respectively.Let the interval [0, 1] be partitioned uniformly into n subintervals [(i − 1)∆x, i∆x], i = 1, 2, . . ., n, then the quadrature points . Then, the discretization form of (2) is given by: The discretization of the functionals J( f ), TV ω,β ( f ), and M( f ) are as follows.
(1) The discretization of J( f ) is given by: where • denotes the vector two-norm.(2) Let e i be the i-th column of the n × n identity matrix, and let: and: We approximate the weighted TV functional TV ω,β ( f ) (cf. (3)) by: where ω i = ω((i − 1/2)∆x) and ψ(x) = x 2 + β (we replace β(∆x) 2 by β for the sake of simplicity).Obviously, if the weighting factors are set to ω i = 1, i = 1, 2, . . ., n − 1, the weighted TV functional is just the traditional TV functional.In this paper, we choose smooth factors by approximating (4): (3) As for M( f ), we simply approximate it by: Here, we omit the factor ∆x for the sake of simplicity.
Therefore, the discretization of Φ TV ( f ) is: and the discretization of Φ M ( f ) is:

Numerical Implementation
In this section, the implementation of the damped Newton (DN) method for minimization of the functions Φ TV (f) and Φ M (f) (cf.( 8) and ( 9)) is given.We first derive the gradients and the Hessians of J(f), TV ω,β (f) and M(f).Lemma 1.Let β and θ be positive constants, and define: (1) The gradient and the Hessian of J(f) = Kf − h 2 are given by: and: respectively.
Then, the gradient and the Hessian of TV ω,β (f) are given by: and: respectively, where: and: The gradient of M(f) is given by: where the partial derivatives ∂M ∂ f i , i = 1, 2, . . ., n, are given by: The Hessian of M(f) is a diagonal matrix given by: where: Proof.Results (1) and ( 2) of the lemma are well known; see [8], Section 8.2. ( , then we have: where φ(x) and ρ(x) are defined by (10).Let θ = θ( β + θ); it can be easily checked that for any v, where L(f) is given by (18).
To obtain the Hessian of TV ω,β , we consider TV ω,β (f + τv + ξw).From the expression of ζ (x), we have that for any v, w: where L(f) and L (f)f are given by ( 18) and ( 19) respectively.Consequently: (4) It is easy to see that the partial derivatives ∂M ∂ f i , i = 1, 2, . . ., n, are given by: which can be rewritten in the form of (21).For the second order partial derivatives, we have: and: Therefore, the Hessian of M is a diagonal matrix given by: where . ., n, are given by (23).
From Formulas ( 11)-( 23), we can get the gradients and the Hessians of Φ TV (f) and Φ M (f).We introduce the following DN method for minimization of Ψ(f) = Φ TV (f) or Φ M (f).Update the approximate solution:

Remark 1.
(1) Since we carry out exact line search (Step 6) in the iterative process, the DN method converges to a local minimal point if the Hessian of Ψ(f) is invertible.In the case that Ψ(f) is singular, modification is required.In our numerical tests in Section 4, no modification is needed.It must be pointed out that for large-scale systems, the most expensive part of the algorithm is solving Hs = −g to obtain s.If the matrix H has a Toeplitz structure, then the conjugate gradient method with the fast Fourier transform (FFT) can be applied to solve the system efficiently; see [31] for details.
(2) For the minimization of Φ TV (f), if we approximate the Hessian of Φ TV (f) by the positive definite matrix K T K + αL(f) (assume that K1 = 0 where 1 is the vector with all elements equal to one) and set the step-size τ ν to 1, we get the Gauss-Newton method.If we obtain τ ν by using line search method, we get a modified Gauss-Newton (MGN) method.In this paper, we also consider the MGN method for minimization of Φ TV (f).Obviously, the MGN method converges to a local minimal point.
To end this section, we state the process for solving Equation (2).We first obtain the minimizer f * TV of Φ TV (f), and then, we obtain a minimizer f * M of Φ M (f) by using the DN method with f * TV as the initial guess.The approach can be summarized as Algorithm 2.
Obtain the minimizer f * TV of Φ TV (f) by using the MGN method or the DN method; Obtain a minimizer f * M of Φ M (f) by using the DN method with f * TV as the initial guess; Output f * M .

Numerical Examples
In this section, we present numerical results for two examples to illustrate that our approach is indeed capable of estimating piecewise-constant solutions to Equation (2).
Besides the numerical solution obtained by using our approach, approximate solutions obtained by using the TV regularization method and the weighted TV regularization method with weighting factors 7)) are presented.In the following tables and figures, we use the symbols weighted total variation Modica-Mortola (WTVMM), TV and WTV to denote the above three methods.
We set the number of subintervals to n = 128.All tests were carried out by using MATLAB, and the termination parameter is set to = 10 −5 ; see Step 8 of Algorithm 1.In the tables, the column "(α, β, θ, γ)" gives relevant parameters; the column "#it" gives the number of iterations (for the WTVMM method, i 1 + i 2 denotes that the WTV method and the MM method require i 1 and i 2 iterations, respectively).The column "error" gives the relative error defined by: where f exact and f app are the vectors of the exact solution and the approximate solution at the quadrature points, respectively.Here, one iteration refers to Steps 4-8 of Algorithm 1.
There are several regularization parameter choosing methods, e.g., the discrepancy principle, the generalized discrepancy principle, the generalized cross validation, the L-curve and the normalized cumulative periodogram; see, for instance, [32,33], ([1], Chapters 5-6), ([8], Chapter 7), ( [34], Chapters 5, 6 and 9).However, we choose the parameters by experiments since there are two additional parameters in the regularization term TV ω,β (f).We consider the following rules in choosing the parameters: (1) Since we use x 2 + β as an approximation of |x|, we set β to a very small positive number.
We note that if β is significantly larger than x 2 , then: As a result, the regularization term x 2 + β leads to an H 1 -like regularization.(2) If θ is much larger than x 2 , then the weighting function ω(x) is close to one, i.e., It follows that the weighted TV regularization is about the same as the traditional TV regularization.Therefore, we should not choose a large value for θ. (3) We choose the value for α by the help of the discrepancy principle and observe if the numerical solution is near piecewise-constant.After choosing a value for α, we consider adjustment of β and θ: if the numerical solution is over-smooth, we choose a smaller value for β or θ; on the other hand, if the numerical solution is oscillating, we choose a larger value for β or θ.
We have tried the DN method and the MGN method for minimization of the functions corresponding to the TV functional (ω i = 1) and the WTV functional (the weighted factors are given by ( 7)) in our numerical tests.Numerical results show that the DN method is better than the MGN method when ω i = 1, and the MGN method is better than the DN method when ω i is given by (7).In the following, we show the numerical results carried out by the DN method for minimization of the function corresponding to the TV functional and those carried out by the MGN method for minimization of the function corresponding to the WTV functional.
Example 1.The kernel of Equation ( 2) is given by: k where σ = 0.05.In this case, the condition number of the matrix K is 7.8654 × 10 18 .The right-hand side g(x) is obtained by setting the exact solution to: where c 1 = 1, c 2 = 2, c 3 = 3.Two noisy right-hand sides, which contain 1% and 10% of white noise, respectively, are tested.
Since the average of c 1 , c 2 and c 3 is two, we choose f 0 = (2, 2, . . ., 2) T as the initial guess for the minimization of the functions corresponding to the TV and WTV functionals, which is not too far from the exact solution.The exact solution and the right-hand side are shown in Figure 1a,b, two noisy right-hand sides are shown in Figure 1c,d, respectively.Numerical solutions obtained by different methods, as well as the point-wise error for all numerical solutions are shown in Figure 2. The relevant parameters, the numbers of iterations and the relative errors are given in Tables 1 and 2.      We have tried a one-stage method for solving Example 1.The best numerical solution we obtained by minimizing the function 1  2 J(f) + αTV 1,β (f) + γ 2 M(f) for Example 1 is shown in Figure 3.We observe from all figures in Figures 2 and 3 that the numerical solution obtained via the WTVMM method is the best.We make the following remarks based on the above numerical results.
(1) For the WTVMM method, the main cost lies in obtaining the numerical solution f * TV , i.e., minimization of Φ TV (f) defined by ( 8).
(2) The numerical solutions obtained by using the WTV method are better than those obtained by the TV method.The numerical solutions obtained by using the WTVMM method can be quite accurate, which can be clearly seen from the relative errors shown in Tables 1 and 2, the numerical solutions shown in Figure 2c and the point-wise errors shown in Figure 2d.(3) If the position of the break points in the solution can be identified by the previous step, the approximate solution obtained by the WTVMM method can be very accurate.In the other hand, since the local minimizer of Φ M obtained by the DN method is sensitive to the initial guess, we may not be able to obtain a good numerical solution if the approximate solution obtained in the previous step cannot identify the break points well.
To demonstrate the influence of the values of parameters on numerical solutions, we present the relative errors of the numerical solutions of the three methods for a considerably different set of the parameter combinations (β, θ) in Tables 3 and 4. We can see from the tables that if the values of the parameters are not too far from the optimal ones, we can obtain satisfactory numerical solutions.Moreover, both values of β and θ affect the quality of the numerical solutions, and a correct choice of β seems more important (see the column corresponding to β = 2.5 × 10 −3 in Table 3 and the one for β = 2.5 × 10 −2 in Table 4).Moreover, we can see that satisfying numerical solutions can be obtained by using different parameter combinations.
In our numerical tests, we choose σ = 0.01.The intensity of the printed barcode and the exact right-hand side are shown in Figure 4a,b, and two noisy right-hand sides are shown in Figure 4c,d, respectively.
Since the average of c 1 and c 2 is 0.5, we choose f 0 = (0.5, 0.5, . . ., 0.5) T as initial guess for minimization of the functions corresponding to the TV and WTV functionals.Numerical solutions obtained by different methods, as well as the the point-wise error for all numerical solutions are shown in Figure 5.The relevant parameters, the number of iterations and the relative errors are given in Tables 5 and 6.Again, from Figure 5 and Tables 5 and 6, one can observe that the numerical solutions obtained by using the WTVMM method are much better than those obtained by using the TV and WTV methods.In fact, we can obtain an approximate barcode of very high quality by using the WTVMM method.

Concluding Remarks
We presented a two-stage numerical method for estimating piecewise-constant solutions for Fredholm integral equations of the first kind.The main work of the two-stage method is the minimization of Φ TV (f), which may require many iterations.Numerical results showed that if the relevant parameters are chosen suitably, the proposed method can obtain satisfying numerical solutions.We will study regularization parameter choosing methods, as well as efficient algorithms for relevant minimization problems for weighted total variation (WTV) methods, including iterative scheme for parameter determination.

Figure 1 .
Figure 1.The exact solution, the right-hand side and noisy right-hand sides for Example 1.(a) The exact solution; (b) the right-hand side; (c) a noisy right-hand side containing 1% of white noise; and (d) a noisy right-hand side containing 10% of white noise.

Figure 4 .
Figure 4.The exact solution, the right-hand side, and noisy right-hand sides for Example 2. (a) The exact solution; (b) the right-hand side; (c) a noisy right-hand side containing 1% of white noise; and (d) a noisy right-hand side containing 10% of white noise.

Figure 5 .
Figure 5.The exact solution (dot), numerical solutions (x) and point-wise errors for the three methods for Example 2 with noisy right-hand sides containing 1% (left) and 10% (right) of white noise.(a) Numerical solutions of TV; (b) numerical solutions of WTV; (c) numerical solutions of WTVMM; and (d) point-wise errors of the numerical solutions.

Table 1 .
Parameters used, number of iterations and relative errors for different methods for Example 1 where the right-hand side contains 1% of white noise.TV: total variation; WTV: weighted total variation; WTVMM: weighted total variation Modica-Mortola.

Table 2 .
Parameters used, number of iterations and relative errors for different methods for Example 1 where the noisy right-hand side contains 10% of white noise.

Table 3 .
Relative errors for TV (1st element), WTV (2nd element) and WTVMM (3rd element) with different β and θ for Example 1 with the right-hand side containing 1% of white noise.

Table 5 .
Parameters used, number of iterations and relative errors for different methods for Example 2 where the right-hand side contains 1% of white noise.

Table 6 .
Parameters used, number of iterations and relative errors for different methods for Example 2 where the right-hand side contains 10% of white noise.