Table of contents

In this paper we review 30 years of developments and applications of the variable projection method for solving separable nonlinear least-squares problems. These are problems for which the model function is a linear combination of nonlinear functions. Taking advantage of this special structure, the method of variable projections eliminates the linear variables obtaining a somewhat more complicated function that involves only the nonlinear parameters. This procedure not only reduces the dimension of the parameter space but also results in a better-conditioned problem. The same optimization method applied to the original and reduced problems will always converge faster for the latter. We present first a historical account of the basic theoretical work and its various computer implementations, and then report on a variety of applications from electrical engineering, medical and biological imaging, chemistry, robotics, vision, and environmental sciences. An extensive bibliography is included. The method is particularly well suited for solving real and complex exponential model fitting problems, which are pervasive in their applications and are notoriously hard to solve.

In this paper we consider the inverse boundary problem for the heat equation Δu(x, t) = ρ(x)∂_tu(x, t) in a bounded domain Ω ⊂ R². We develop and test numerically an algorithm of an approximate reconstruction of the unknown ρ(x). This algorithm is based on the moments method for the heat equation developed by Kawashita, Kurylev and Soga.

We propose an extension to the functional modelling methods described by Dawid and Stone (1982 Ann. Stat.10 1119–38) that leads naturally to a method for selecting vague parameter priors for Bayesian analyses involving stochastic population models. Motivated by applications from quantum optics and epidemiology, we focus on analysing observed sequences of event times obeying a non-homogeneous Poisson process, although the techniques are more widely applicable. The extended functional modelling approach is illustrated for the particular case of Bayesian estimation of the death rate in the immigration–death model from observation of the death times only. It is shown that the prior selected naturally leads to a well defined posterior density for parameters and avoids some undesirable pathologies reported by Gibson and Renshaw (2001a Inverse Problems17 455–66, 2001b Stat. Comput.11 347–58) for the case of exponential priors. Some limitations of the approach are also discussed.

We consider the inverse scattering problem of determining the shape of an infinite cylinder having an open arc as cross section from a knowledge of the TM-polarized scattered electromagnetic field corresponding to time-harmonic incident plane waves propagating from arbitrary directions. We assume that the arc is a (possibly) partially coated perfect conductor and develop the linear sampling method, which was originally developed for solving the inverse scattering problem for obstacles with nonempty interior, to include the above case of obstacles with empty interior.

We consider the problem of determining the shear modulus of a linear-elastic, incompressible medium given boundary data and one component of the displacement field in the entire domain. The problem is derived from applications in quantitative elasticity imaging. We pose the problem as one of minimizing a functional and consider the use of gradient-based algorithms to solve it. In order to calculate the gradient efficiently we develop an algorithm based on the adjoint elasticity operator. The main cost associated with this algorithm is equivalent to solving two forward problems, independent of the number of optimization variables. We present numerical examples that demonstrate the effectiveness of the proposed approach.

We consider the inverse problem of recovering a 2D periodic structure from near-field measurements above the structure. First, following Bruckner et al (2001 Preprint No 682 Weierstrass Institute, Berlin; 2003 Proc. 3rd ISAAC Congress (Berlin, 2001) at press), the inverse problem is reformulated as an optimization problem which consists of two parts: a linear severely ill-posed problem and a nonlinear well-posed one. Then, unlike the work of Bruckner et al, the two problems are solved separately to diminish the computational effort by exploiting their special properties. Numerical results for exact and noisy data demonstrate the practicability of the inversion algorithm.

A sequential minimization algorithm for the numerical solution of inverse problems of frequency sounding is presented. This algorithm is based on the concept of convexification of a multiextremal objective function proposed recently by the present authors. A key point in the sequential minimization algorithm is that, unlike conventional layer-stripping algorithms, it provides the stable approximate solution via minimization of a finite sequence of strictly convex objective functions resulting from applying the nonlinear weighted least squares method with Carleman's weight functions (CWFs). Another advantage of the proposed algorithm is that the starting vectors for the descent methods of minimization are directly determined from the data eliminating the uncertainty inherent to the local methods, such as the gradient or Newton-like methods. The one-dimensional inverse model of magnetotelluric frequency sounding is selected to demonstrate its computational feasibility. Based on the localizing property of CWFs, it is proven that the distance between the approximate and 'exact' solutions is small if the approximation error is small. The computational experiments with several realistic and synthetic marine shallow water configurations are presented to demonstrate the computational feasibility of the proposed algorithm.

This paper proposes a method for reconstructing the positions, strengths, and number of point sources in a three-dimensional (3D) Poisson field from boundary measurements. Algebraic relations are obtained, based on multipole moments determined by the sources and data on the boundary of a domain. To solve for the source parameters with efficient use of data, we select the necessary number of equations from them in the following two ways: (1) the use of those starting from lower-degree multipole moments; and (2) the use of combined ones involving infinitely higher-degree multipole moments. We show that both methods are based on the projection of 3D sources onto a two-dimensional space: the xy-plane for the first one and the Riemann sphere which is set to contain the domain for the second one. We also show that they share the same fundamental equations which can be solved by a procedure proposed by El-Badia and Ha-Duong (2000 Inverse Problems16 651–63). Numerical simulations show that projection onto the xy-plane is more appropriate for sources scattered in the middle of the domain, whereas projection onto the Riemann sphere is more appropriate for sources concentrated close to the boundary of the domain. We also give an appropriate method of measurement for the Riemann sphere projection.

We present the asymptotic expansion of the solution to a diffusion equation with a finite number of absorbing inclusions of small volume. We use the first few terms in this expansion measured at the domain boundary to reconstruct the absorption parameters of the inclusions and certain geometrical characteristics. We demonstrate theoretically and numerically that the number of inclusions, their location and their capacity can be reconstructed in a stable way even from moderately noisy data. The reconstruction of the absorption parameter, which is important in optical tomography to discriminate between healthy and unhealthy tissues, requires us however to have far less noisy data. Since the reconstruction of absorption maps from boundary measurements is an extremely ill posed problem, the method of asymptotic expansions of small volume inclusions provides a useful framework to decide which information can be reconstructed from boundary measurements with a given noise level.

Optical tomography (OT) recovers the cross-sectional distribution of optical parameters inside a highly scattering medium from information contained in measurements that are performed on the boundaries of the medium. The image reconstruction problem in OT can be considered as a large-scale optimization problem, in which an appropriately defined objective function needs to be minimized. In the simplest case, the objective function is the least-square error norm between the measured and the predicted data. In biomedical applications that apply near-infrared light as the probing tool the predictions are obtained from a model of light propagation in tissue. Gradient techniques are commonly used as optimization methods, which employ the gradient of the objective function with respect to the optical parameters to find the minimum. Conjugate gradient (CG) techniques that use information about the first derivative of the objective function have shown some good results in the past. However, this approach is frequently characterized by low convergence rates. To alleviate this problem we have implemented and studied so-called quasi-Newton (QN) methods, which use approximations to the second derivative. The performance of the QN and CG methods are compared by utilizing both synthetic and experimental data.

The solution of an Urison-type integral equation is the fixed point of an associated integral operator. The direct problem involves determining this fixed point when the operator is specified. We consider instead the following inverse problem: determine such an operator, in a restricted class, with fixed point close to a given target function, which may be an approximation or an interpolation of data points. We establish that these operators are contractive on an appropriate complete metric space and use Banach's fixed point theorem and related results, including the collage theorem, to develop a method of solution for the inverse problem. The technique involves the minimization of the distance between the target solution and its image under the integral operator, called the 'collage distance'. We discuss some ways to further reduce this distance.

Localization of an acoustic source in the ocean is often limited by lack of knowledge of the physical properties of the environment, such as seabed geoacoustic parameters and water-column sound-speed profile (SSP). Quantifying environmental uncertainties and how they transfer to uncertainties for source localization represent important problems that are addressed in this paper using Bayesian inference theory. Metropolis Gibbs' sampling is applied to estimate the uncertainties for environmental inversion in the form of marginal probability distributions, covariances and credibility intervals. Heat-bath sampling is applied to source localization with environmental uncertainties, with the resulting localization uncertainty quantified in terms of the joint marginal probability distribution for source range and depth (i.e. a probability ambiguity surface). Localization uncertainties are examined as a function of uncertainties in geoacoustic parameters and SSP.

We report on a new iterative method for regularizing a nonlinear operator equation in Hilbert spaces. The proposed TIGRA algorithm is a combination of Tikhonov regularization and a gradient method for minimizing the Tikhonov functional. Under the assumptions that the operator F is twice continuous Fréchet differentiable with a Lipschitz-continuous first derivative and that the solution of the equation F (x) = y fulfils a smoothness condition, we will give a convergence rate result. Finally we present some applications and a numerical result for the reconstruction of the activity function in single-photon-emission computed tomography.

In this paper, we prove a uniqueness theorem for the potential V (x) of the following Schrödinger operator H = −Δ + q(|x|) + V (x) in ², where q(|x|) is a known increasing radial potential satisfying lim _|x|→+∞q(|x|) = +∞ and V (x) is a bounded potential.

The mathematical formulation and numerical aspects of the inverse problem of a three-dimensional (3D) x-ray imaging method are introduced. The imaging method relies on measuring radiation scattered at two directions orthogonal to an incident beam that scans one side of the object, in addition to the traditionally recorded transmitted radiation. From this set of measurements, the inverse problem reconstructs three images: the attenuation coefficient at the incident energy, the attenuation coefficient at a scattering energy and the scattering coefficients (with the latter directly related to the electron, and hence mass density). The algorithms for this triple imaging process are presented and their numerical characteristics are discussed. With the aid of regularization and an iterative learning process, stable and bounded solutions are obtained. This 3D image-reconstruction process was able to recover the domain of interest and distinguish between different materials.