We consider in this paper a 1D Coefficient Inverse Problem (CIP) for a hyperbolic PDE and construct a globally convergent numerical method for it. Unlike all previous works of this group on the convexification method for CIPs, which are cited below in this section, we construct in this paper a sequence of linearized boundary value problems with overdetermined boundary data. Each of these problems is solved using a weighted version of the Quasi-Reversibility Method (QRM). The weight is the Carleman Weight Function (CWF), i.e. this is the function, which is involved as the weight in the Carleman estimate for the corresponding PDE operator, see, e.g. [7, 4, 5, 22, 18, 31, 38, 53] for Carleman estimates. Thus, we call this method “Carleman Quasi-Reversibility Method” (CQRM). The key convergence estimate for this sequence is similar with the one of the conventional contraction mapping principle, see the first item of Remarks 8.1 in section 8. This explains the title of our paper. That estimate implies the global convergence of that sequence to the true solution of our CIP. Furthermore, that estimate implies a rapid convergence. As a result, computations are done here in real time, also, see Remark 11.1 in section 11 as well as section 12. The real time computation is obviously important for our target application, which is described below. The numerical method of this paper can be considered as the second generation of the convexification method for CIPs.
The convexification method ends up with the minimization of a weighted globally strictly convex Tikhonov functional, in which the weight is the CWF. The convexification has been developed by our research group since the first two publications in 1995  and 1997 . Only the theory was presented in those two originating papers. The start for numerical studies of the convexification was given in the paper . Thus, in follow up publications [24, 25, 30, 16, 14, 15, 24, 25, 27, 28, 29, 32, 33, 30, 40, 50, 51] the theoritical results for the convexification are combined with the numerical ones. In particular, some of these papers work with experimentally collected backscattering data [14, 15, 26, 32, 33]. We also refer to the recently published book of Klibanov and Li , in which many of these results are described.
The CIP of our paper has a direct application in the problem of the standoff detection and identification of antipersonnel land mines and improvised explosive devices (IEDs). Thus, in the numerical part of this publication, we present results of the numerical performance of our method for both computationally simulated and experimentally collected data for targets mimicking antipersonnel land mines and IEDs. The experimental data were collected in the field, which is more challenging than the case of the data collected in a laboratory [14, 15, 26, 33]. The data used in this paper were collected by the forward looking radar of the US Army Research Laboratory . Since these data were described in the previous publications of our research group [13, 23, 24, 25, 30, 35, 36, 51], then we do not describe them here. Those previous publications were treating these experimental data by several different numerical methods for 1D CIPs. We also note that the CIP of our paper was applied in [32, 33] to the nonlinear synthetic aperture (SAR) imaging, including the cases of experimentally collected data.
Our goal is to compute approximate values of dielectric constants of targets. We point out that our experimental data are severely underdetermined
. Indeed, any target is a 3D object. On the other hand, we have only one experimentally measured time resolved curve for each target. Therefore, we can compute only a sort of an average value of the dielectric constant of each target. This explains the reason of our mathematical modeling here by a 1D hyperbolic PDE rather than by its 3D analog. This mathematical model, along with its frequency domain analog, was considered in all above cited previous publications of our research group, which were working with the experimental data of this paper by a variety of globally convergent inversion techniques.
An estimate of the dielectric constant of a target cannot differentiate with a good certainty between an explosive and a non explosive. However, we hope that these estimates, combined with current targets classification procedures, might potentially result in a decrease of the false alarm rate for those targets. We believe therefore that our results for experimental data might help to address an important application to the standoff detection and identification of antipersonnel land mines and IEDs.
Given a numerical method for a CIP, we call this method globally convergent if:
There is a theorem claiming that this method delivers at least one point in a sufficiently small neighborhood of the true solution of that CIP without relying on a good initial guess about this solution.
This theorem is confirmed numerically.
CIPs are both highly nonlinear and ill-posed. As a result, conventional least squares cost functionals for them are non convex and typically have many local minima and ravines, see, e.g.  for a numerical example of this phenomenon. Conventional numerical methods for CIPs are based on minimizations of those functionals, see, e.g. [8, 11, 12]. The goal of the convexification is to avoid local minima and ravines and to provide accurate and reliable solutions of CIPs.
In the convexification, a CWF is a part of a certain least squares cost functional , where is the parameter of the CWF. The presence of the CWF ensures that, with a proper choice of the parameter the functional is strictly convex on a certain bounded convex set of an arbitrary diameter in a Hilbert space. It was proven later in  that the existence and uniqueness of the minimizer of on are also guaranteed. Furthermore, the gradient projection method of the minimization of on converges globally to that minimizer if starting at an arbitrary point of . In addition, if the level of noise in the data tends to zero, then the gradient projection method delivers a sequence, which converges to the true solution. Thus, since is an arbitrary number, then this is the global convergence, as long as numerical studies confirm this theory. Furthermore, it was recently established in [33, 40] that the gradient descent method of the minimization of on also converges globally to the true solution.
First, we come up with the same boundary value problem (BVP) for a nonlinear and nonlocal hyperbolic PDE as we did in the convexification method for the same 1D CIP in . The solution of this BVP immediately provides the target unknown coefficient. However, when solving this BVP, rather than constructing a globally strictly convex cost functional for this BVP, as it was done in  and other works on the convexification, we construct a sequence of linearized BVPs. Just as the original BVP, each BVP of this sequence involves a hyperbolic operator with non local terms and overdetermined boundary conditions.
Because of non local terms and the overdetermined boundary conditions, we solve each BVP of this sequence by the new version of the QRM, the Carleman Quasi-Reversibility Method (CQRM). Indeed, the classical QRM solves ill-posed/overdetermined BVPs for linear PDEs, see  for the originating work as well as, e.g. [19, 31] for updates. Convergence of the QRM-solution to the true solution is proven via an appropriate Carleman estimates [19, 31]. However, a CWF is not used in the optimizing quadratic functional of the conventional QRM of [19, 31]. Unlike this, a CWF is involved in CQRM. It is exactly the presence of the CWF, which enables us to derive the above mentioned convergence estimate, which is an analog of the one of the contraction principle. This allows us, in turn to establish convergence rate of the sequence of CQRM-solutions to the true solution of that BVP. The latter result ensures that our method is a globally convergent one.
Various versions of the second generation of the convexification method involving CQRM were previously published in [2, 3, 45, 44, 39]. In these publications, globally convergent numerical methods for some nonlinear inverse problems are presented. In [2, 3] CIPs for hyperbolic PDEs in are considered. The main difference between [2, 3] and our work is that it is assumed in [2, 3] that one of initial conditions is not vanishing everywhere in the closed domain of interest. In other words, papers [2, 3] work exactly in the framework of the Bukhgeim-Klibanov method, see  for the originating work on this method and, e.g. [4, 5, 17, 22, 18, 31, 53] for some follow up publications. On the other hand, in all above cited works on the convexification, including this one, either the initial condition is the function, or a similar requirement holds for the Helmholtz equation. The only exception is the paper , which also works within the framework of the method of . We refer to papers [10, 34] for some numerical methods for CIPs with the Dirichlet-to-Neumann map data. In these works the number of free variables in the data exceeds the number of free variables in the unknown coefficient, In our paper Also, in all other above cited works on the convexification.
There is a classical Gelfand-Levitan method  for solutions of 1D CIPs for some hyperbolic PDEs. This method does not rely on the optimization, and, therefore, avoids the phenomenon of local minima. It reduces the original CIP to a linear integral equation of the second kind. This equation is called “Gelfand-Levitan equation” (GL). However, the questions about uniqueness and stability of the solution of GL for the case of noisy data is open, see, e.g. Lemma 2.4 in the book [48, Chapter 2]. This lemma is valid only in the case of noiseless data. However, the realistic data are always noisy. In addition, it was demonstrated numerically in  that GL cannot work well for exactly the same experimental data as the ones we use in the current paper. On the other hand, it was shown in [24, 25, 51] that the convexification works well with these data. The same is true for the method of the current paper.
Uniqueness and Lipschitz stability theorems of the CIP considered here are well known. Indeed, it was shown in, e.g.  that, using the same change of variables as the one considered in section 2 below, one can reduce our CIP to a similar CIP for the equation with the unknown coefficient We refer to [48, Theorem 2.6 of Chapter 2] for the Lipschitz stability estimate for the latter CIP. In addition, both uniqueness and Lipschitz stability results for our CIP actually follow from Theorem 8.1 below as well as from the convergence analysis in two other works of our research group on the convexification method for this CIP [30, 51].
This paper is arranged as follows. In section 2 we state both forward and inverse problems. In section 3 we derive a nonlinear BVP with non local terms. In section 4 we describe iterations of our method to solve that BVP. In section 5 we formulate the Carleman estimate for the principal part of the PDE operator of that BVP. In section 6 we prove the strong convexity of a functional of section 5 on an appropriate bounded set. In section 7 we formulate two methods for finding the unique minimizer of that functional: gradient descent method and gradient projection method and prove the global convergence to the minimizer for each of them. In section 8 we establish the contraction mapping property and prove the global convergence of the method of section 4. In section 9 we formulate two more global convergence theorems, which follow from results of sections 7 and 8. Numerical results with simulated and experimental data are presented in Section 10 and 11 respectively. Concluding remarks are given in section 12.
2 Statements of Forward and Inverse Problems
Below all functions are real valued ones. Let be a known number, be the spatial variable and the function , represents the spatially distributed dielectric constant. We assume that
where is a small number. Let be a positive number. We consider the following Cauchy problem for a 1D hyperbolic PDE with a variable coefficient in the principal part of the operator:
Let be the travel time needed for the wave to travel from the point source at to the point ,
By (2.5) the following 1D analog of the eikonal equation is valid:
Let be the Heaviside function,
Lemma 2.1 [30, Lemma 2.1]. For , the function has the form
where the function Furthermore,
Lemma 2.2 [Absorbing boundary condition] [30, 50]. Let be the number in (2.1). Let and be two arbitrary numbers. Then the solution of Forward Problem (2.3), (2.4) satisfies the absorbing boundary condition at and , i.e.
We are interested in the following inverse problem:
Coefficient Inverse Problem (CIP). Suppose that the following two functions , are known:
Determine the function for , assuming that the number in (2.1) is known.
3 A Boundary Value Problem for a Nonlinear PDE With Non-Local Terms
We now introduce a change of variable
Furthermore, it follows from (2.8)
Equation, (3.5) is a nonlinear PDE with respect to the function with nonlocal terms and . We now need to obtain boundary conditions for the function By (2.2) and (2.5) for Hence, (2.11), (2.12) and (3.1) lead to
We will solve equation (3.5) in the rectangle
In addition, we need and we also need to bound the norm from the above. Let be an arbitrary number. We define the set as
We assume below that
By embedding theorem and
where the constant depends only on the domain Using (3.10), we define the function as:
Thus, we focus below on the numerical solution of the following problem:
4 Numerical Method for Problem 3.1
4.1 The function
BVP (4.2), (4.3) has overdetermined boundary conditions (4.3). Typically, QRM works well for BVPs for PDEs with overdetermined boundary conditions [19, 31]. Therefore, we solve BVP (4.2), (4.3) via CQRM. This means that we consider the following minimization problem:
Minimization Problem Number . Assuming (3.12), minimize the functional on the set
where is the parameter, and is the regularization parameter. Both parameters will be chosen later.
Theorem 6.1 guarantees that, for appropriate values of parameters there exists unique minimizer of the functional
4.2 The function for
Assume that functionals are defined for and their minimizers functions are constructed already, all for the same values of parameters . Replace in (3.17) with with with and with Then problem (3.17), (3.18) becomes a linear one with respect to the function
Minimization Problem Number . Assuming (3.12), minimize the functional on the set
5 The Carleman Estimate
In this section we formulate the Carleman estimate, which is the main tool of our construction. This estimate follows from Theorem 3.1 of  as well as from (3.15). Let be an arbitrary function and let the function be constructed from the function as in (3.14). Consider the operator
Theorem 5.1 (Carleman estimate ). There exists a number depending only on such that for any there exists a sufficiently large number depending only on , such that for all and for all functions the following Carleman estimate holds:
Remark 5.1. Here and everywhere below denotes different constants depending only on listed parameters.
6 Strict Convexity of Functional (4.8) on the Set , Existence and Uniqueness of Its Minimizer
Functional (4.8) is quadratic. We prove in this section that it is strictly convex on the set In addition, we prove existence and uniqueness of its minimizer on this set. Although similar results were proven in many of the above cited publications on the convexification, see, e.g.  for the closest one, there are some peculiarities here, which are important for our convergence analysis, see Remarks 6.1, 6.2 below.
Introduce the subspace as:
Denote the scalar product in the space Also, denote
Theorem 6.1. 1. For any set of parameters and for any the functional has the Fréchet derivative The formula for is:
2. This derivative is Lipschitz continuous in , i.e. there exists a constant such that
3. Let and be the numbers of Theorem 5.1 Then there exists a sufficiently large constant
depending only on listed parameters such that for all and for all the functional is strictly convex on the set More precisely, let be an arbitrary function and also let the function Then the following inequality holds:
4. For any there exists unique minimizer
of the functional on the set and the following inequality holds:
where denotes the scalar product in
Remark 6.1. Even though the expression in the right hand side of (6.4) is linear with respect to the function we cannot use Riesz theorem here to prove existence and uniqueness of the minimizer at which Rather, all what we can prove is (6.9). This is because we need to ensure that the function .
Proof of Theorem 6.1. Since both function satisfy the same boundary conditions, then
The expression in the first line of (6.11) coincides with expression in the right hand side of (6.4). In fact, this is a bounded linear functional mapping in . Therefore, by Riesz theorem, there exists unique function such that
In addition, it is clear from (6.11) and (LABEL:5.9) that
By trace theorem there exists a constant depending only on the domain such that
Since the regularization parameter then we can choose so large that
Hence, for these values of the expression in the last line of (6.13) can be estimated from the below as:
This proves (6.7). The existence and uniqueness of the minimizer and inequality (6.9) follow from (6.7) as well as from a combination of Lemma 2.1 and Theorem 2.1 of , also see [42, Chapter 10, section 3].
Remark 6.2. Since the functional is quadratic, then its strict convexity on the whole space follows immediately from the presence of the regularization term in it. However, in addition to the claim of its strict convexity, we actually need in our convergence analysis below those terms in the right hand side of the strict convexity estimate (6.7), which are different from the term These terms are provided by Carleman estimate (6.7). The condition of Theorem 6.1 is imposed to dominate the negative term in the last line of (6.7).
7 How to Find the Minimizer
Since we search for the minimizer of functional (4.8) on the bounded set rather than on the whole space then we cannot just use Riesz theorem to find this minimizer, also see  and [42, Chapter 10, section 3] for the case of finding a minimizer of a strictly convex functional on a bounded set. Two ways of finding the minimizer are described in this section.
7.1 Gradient descent method
Theorem 7.1 claims the global convergence of the gradient descent method (7.1) in the case when
Theorem 7.1. Let the number be the one defined in (6.6). Let . For this value of let be the unique minimizer of the functional on the set with the regularization parameter (Theorem 6.1) Assume that the function For each choose the starting point of the gradient descent method (7.1) as Then, there exists a number such that for any all functions Furthermore, there exists a number such that the following convergence estimates are valid:
7.2 Gradient projection method
Suppose now that there is no information on whether or not the function In this case we construct the gradient projection method. We introduce the function since it is easy to construct the projection operator on a ball with the center at
Suppose that there exists a function such that
This function exists, for example, in the case when the function Indeed, as one of examples of this function (among many), one can take, e.g.
where the function is such that
The existence of such functions is well known from the Analysis course. Denote