Search

Understanding the orientation of geological structures is crucial for analyzing the complexity of the Earths' subsurface. For instance, information about geological structure orientation can be incorporated into local anisotropic regularization methods as a valuable tool to stabilize the solution of inverse problems and produce geologically plausible solutions. We introduce a new variational method that employs the alternating direction method of multipliers within an alternating minimization scheme to jointly estimate orientation and model parameters in both 2D and 3D inverse problems. Specifically, the proposed approach adaptively integrates recovered information about structural orientation, enhancing the effectiveness of anisotropic Tikhonov regularization in recovering geophysical parameters. The paper also discusses the automatic tuning of algorithmic parameters to maximize the new method's performance. The proposed algorithm is tested across diverse 2D and 3D examples, including structure-oriented denoising and trace interpolation. The results show that the algorithm is robust in solving the considered large and challenging problems, alongside efficiently estimating the associated tilt field in 2D cases and the dip, strike, and tilt fields in 3D cases. Synthetic and field examples show that the proposed anisotropic regularization method produces a model with enhanced resolution and provides a more accurate representation of the true structures.

Regularization techniques are necessary to compute meaningful solutions to discrete ill-posed inverse problems. The well-known 2-norm Tikhonov regularization method equipped with a discretization of the gradient operator as regularization operator penalizes large gradient components of the solution to overcome instabilities. However, this method is homogeneous, i.e., it does not take into account the orientation of the regularized solution and therefore tends to smooth the desired structures, textures and discontinuities, which often contain important information. If the local orientation field of the solution is known, a possible way to overcome this issue is to implement local anisotropic regularization by penalizing weighted directional derivatives. In this paper, considering problems that are inherently two-dimensional, we propose to automatically and simultaneously recover the regularized solution and the local orientation parameters (used to define the anisotropic regularization term) by solving a bilevel optimization problem. Specifically, the lower level problem is Tikhonov regularization equipped with local anisotropic regularization, while the objective function of the upper level problem encodes some natural assumptions about the local orientation parameters and the Tikhonov regularization parameter. Application of the proposed algorithm to a variety of inverse problems in imaging (such as denoising, deblurring, tomography and Dix inversion), with both real and synthetic data, shows its effectiveness and robustness.

In this paper, we describe TRIPs-Py, a new Python package of linear discrete inverse problems solvers and test problems. The goal of the package is two-fold: 1) to provide tools for solving small and large-scale inverse problems, and 2) to introduce test problems arising from a wide range of applications. The solvers available in TRIPs-Py include direct regularization methods (such as truncated singular value decomposition and Tikhonov) and iterative regularization techniques (such as Krylov subspace methods and recent solvers for $\ell_p$-$\ell_q$ formulations, which enforce sparse or edge-preserving solutions and handle different noise types). All our solvers have built-in strategies to define the regularization parameter(s). Some of the test problems in TRIPs-Py arise from simulated image deblurring and computerized tomography, while other test problems model realistic problems in dynamic computerized tomography. Numerical examples are included to illustrate the usage as well as the performance of the described methods on the provided test problems. To the best of our knowledge, TRIPs-Py is the first Python software package of this kind, which may serve both research and didactical purposes., Comment: 27 pages, 10 figures, 3 tables

Variational regularization is commonly used to solve linear inverse problems, and involves augmenting a data fidelity by a regularizer. The regularizer is used to promote a priori information and is weighted by a regularization parameter. Selection of an appropriate regularization parameter is critical, with various choices leading to very different reconstructions. Classical strategies used to determine a suitable parameter value include the discrepancy principle and the L-curve criterion, and in recent years a supervised machine learning approach called bilevel learning has been employed. Bilevel learning is a powerful framework to determine optimal parameters and involves solving a nested optimization problem. While previous strategies enjoy various theoretical results, the well-posedness of bilevel learning in this setting is still an open question. In particular, a necessary property is positivity of the determined regularization parameter. In this work, we provide a new condition that better characterizes positivity of optimal regularization parameters than the existing theory. Numerical results verify and explore this new condition for both small and high-dimensional problems., Comment: 34 pages, 11 figures. Version for publication

Full Waveform Inversion (FWI) is a standard algorithm in seismic imaging. Its implementation requires the a priori choice of a number of "design parameters", such as the positions of sensors for the actual measurements and one (or more) regularisation weights. In this paper we describe a novel algorithm for determining these design parameters automatically from a set of training images, using a (supervised) bilevel learning approach. In our algorithm, the upper level objective function measures the quality of the reconstructions of the training images, where the reconstructions are obtained by solving the lower level optimisation problem -- in this case FWI. Our algorithm employs (variants of) the BFGS quasi-Newton method to perform the optimisation at each level, and thus requires the repeated solution of the forward problem -- here taken to be the Helmholtz equation. This paper focuses on the implementation of the algorithm. The novel contributions are: (i) an adjoint-state method for the efficient computation of the upper-level gradient; (ii) a complexity analysis for the bilevel algorithm, which counts the number of Helmholtz solves needed and shows this number is independent of the number of design parameters optimised; (iii) an effective preconditioning strategy for iteratively solving the linear systems required at each step of the bilevel algorithm; (iv) a smoothed extraction process for point values of the discretised wavefield, necessary for ensuring a smooth upper level objective function. The algorithm also uses an extension to the bilevel setting of classical frequency-continuation strategies, helping avoid convergence to spurious stationary points. The advantage of our algorithm is demonstrated on a problem derived from the standard Marmousi test problem.

This paper presents a couple of preconditioning techniques that can be used to enhance the performance of iterative regularization methods applied to image deblurring problems with a variety of point spread functions (PSFs) and boundary conditions. More precisely, we first consider the anti-identity preconditioner, which symmetrizes the coefficient matrix associated to problems with zero boundary conditions, allowing the use of MINRES as a regularization method. When considering more sophisticated boundary conditions and strongly nonsymmetric PSFs, the anti-identity preconditioner improves the performance of GMRES. We then consider both stationary and iteration-dependent regularizing circulant preconditioners that, applied in connection with the anti-identity matrix and both standard and flexible Krylov subspaces, speed up the iterations. A theoretical result about the clustering of the eigenvalues of the preconditioned matrices is proved in a special case. The results of many numerical experiments are reported to show the effectiveness of the new preconditioning techniques, including when considering the deblurring of sparse images.

Combinations of spectroscopic analysis and microscopic techniques are used across many disciplines of scientific research, including material science, chemistry and biology. X-ray spectromicroscopy, in particular, is a powerful tool used for studying chemical state distributions at the micro and nano scales. With the beam fixed, a specimen is typically rastered through the probe with continuous motion and a range of multimodal data is collected at fixed time intervals. The application of this technique is limited in some areas due to: long scanning times to collect the data, either because of the area/volume under study or the compositional properties of the specimen; and material degradation due to the dose absorbed during the measurement. In this work, we propose a novel approach for reducing the dose and scanning times by undersampling the raster data. This is achieved by skipping rows within scans and reconstructing the x-ray spectromicroscopic measurements using low-rank matrix completion. The new method is robust and allows for x 5-6 reduction in sampling. Experimental results obtained on real data are illustrated., Comment: 22 pages,11 figures

This paper considers large-scale linear ill-posed inverse problems whose solutions can be represented as sums of smooth and piecewise constant components. To solve such problems we consider regularizers consisting of two terms that must be balanced. Namely, a Tikhonov term guarantees the smoothness of the smooth solution component, while a total-variation (TV) regularizer promotes blockiness of the non-smooth solution component. A scalar parameter allows to balance between these two terms and, hence, to appropriately separate and regularize the smooth and non-smooth components of the solution. This paper proposes an efficient algorithm to solve this regularization problem by the alternating direction method of multipliers (ADMM). Furthermore, a novel algorithm for automatic choice of the balancing parameter is introduced, using robust statistics. The proposed approach is supported by some theoretical analysis, and numerical experiments concerned with different inverse problems are presented to validate the choice of the balancing parameter.

In this work, we investigate various approaches that use learning from training data to solve inverse problems, following a bi-level learning approach. We consider a general framework for optimal inversion design, where training data can be used to learn optimal regularization parameters, data fidelity terms, and regularizers, thereby resulting in superior variational regularization methods. In particular, we describe methods to learn optimal $p$ and $q$ norms for ${\rm L}^p-{\rm L}^q$ regularization and methods to learn optimal parameters for regularization matrices defined by covariance kernels. We exploit efficient algorithms based on Krylov projection methods for solving the regularized problems, both at training and validation stages, making these methods well-suited for large-scale problems. Our experiments show that the learned regularization methods perform well even when there is some inexactness in the forward operator, resulting in a mixture of model and measurement error.

We consider efficient methods for computing solutions to dynamic inverse problems, where both the quantities of interest and the forward operator (measurement process) may change at different time instances but we want to solve for all the images simultaneously. We are interested in large-scale ill-posed problems that are made more challenging by their dynamic nature and, possibly, by the limited amount of available data per measurement step. To remedy these difficulties, we apply regularization methods that enforce simultaneous regularization in space and time (such as edge enhancement at each time instant and proximity at consecutive time instants) and achieve this with low computational cost and enhanced accuracy. More precisely, we develop iterative methods based on a majorization-minimization (MM) strategy with quadratic tangent majorant, which allows the resulting least squares problem to be solved with a generalized Krylov subspace (GKS) method; the regularization parameter can be defined automatically and efficiently at each iteration. Numerical examples from a wide range of applications, such as limited-angle computerized tomography (CT), space-time image deblurring, and photoacoustic tomography (PAT), illustrate the effectiveness of the described approaches., Comment: 30 pages, 10 figure, 3 tables

This paper is concerned with the regularization of large-scale discrete inverse problems by means of inexact Krylov methods. Specifically, we derive two new inexact Krylov methods that can be efficiently applied to unregularized or Tikhonov-regularized least squares problems, and we study their theoretical properties, including links with their exact counterparts and strategies to monitor the amount of inexactness. We then apply the new methods to separable nonlinear inverse problems arising in blind deblurring. In this setting inexactness stems from the uncertainty in the parameters defining the blur, which may be recovered using a variable projection method leading to an inner-outer iteration scheme (i.e., one cycle of inner iterations is performed to solve one linear deblurring subproblem for any intermediate values of the blurring parameters computed by a nonlinear least squares solver). The new inexact solvers can naturally handle varying inexact blurring parameters while solving the linear deblurring subproblems, allowing for a much reduced number of total iterations and substantial computational savings with respect to their exact counterparts.

This paper surveys an important class of methods that combine iterative projection methods and variational regularization methods for large-scale inverse problems. Iterative methods such as Krylov subspace methods are invaluable in the numerical linear algebra community and have proved important in solving inverse problems due to their inherent regularizing properties and their ability to handle large-scale problems. Variational regularization describes a broad and important class of methods that are used to obtain reliable solutions to inverse problems, whereby one solves a modified problem that incorporates prior knowledge. Hybrid projection methods combine iterative projection methods with variational regularization techniques in a synergistic way, providing researchers with a powerful computational framework for solving very large inverse problems. Although the idea of a hybrid Krylov method for linear inverse problems goes back to the 1980s, several recent advances on new regularization frameworks and methodologies have made this field ripe for extensions, further analyses, and new applications. In this paper, we provide a practical and accessible introduction to hybrid projection methods in the context of solving large (linear) inverse problems.

We present a new inner-outer iterative algorithm for edge enhancement in imaging problems. At each outer iteration, we formulate a Tikhonov-regularized problem where the penalization is expressed in the 2-norm and involves a regularization operator designed to improve edge resolution as the outer iterations progress, through an adaptive process. An efficient hybrid regularization method is used to project the Tikhonov-regularized problem onto approximation subspaces of increasing dimensions (inner iterations), while conveniently choosing the regularization parameter (by applying well-known techniques, such as the discrepancy principle or the ${\mathcal L}$-curve criterion, to the projected problem). This procedure results in an automated algorithm for edge recovery that does not involve regularization parameter tuning by the user, nor repeated calls to sophisticated optimization algorithms, and is therefore particularly attractive from a computational point of view. A key to the success of the new algorithm is the design of the regularization operator through the use of an adaptive diagonal weighting matrix that effectively enforces smoothness only where needed. We demonstrate the value of our approach on applications in X-ray CT image reconstruction and in image deblurring, and show that it can be computationally much more attractive than other well-known strategies for edge preservation, while providing solutions of greater or equal quality.

This paper introduces new solvers for the computation of low-rank approximate solutions to large-scale linear problems, with a particular focus on the regularization of linear inverse problems. Although Krylov methods incorporating explicit projections onto low-rank subspaces are already used for well-posed systems that arise from discretizing stochastic or time-dependent PDEs, we are mainly concerned with algorithms that solve the so-called nuclear norm regularized problem, where a suitable nuclear norm penalization on the solution is imposed alongside a fit-to-data term expressed in the 2-norm: this has the effect of implicitly enforcing low-rank solutions. By adopting an iteratively reweighted norm approach, the nuclear norm regularized problem is reformulated as a sequence of quadratic problems, which can then be efficiently solved using Krylov methods, giving rise to an inner-outer iteration scheme. Our approach differs from the other solvers available in the literature in that: (a) Kronecker product properties are exploited to define the reweighted 2-norm penalization terms; (b) efficient preconditioned Krylov methods replace gradient (projection) methods; (c) the regularization parameter can be efficiently and adaptively set along the iterations. Furthermore, we reformulate within the framework of flexible Krylov methods both the new inner-outer methods for nuclear norm regularization and some of the existing Krylov methods incorporating low-rank projections. This results in an even more computationally efficient (but heuristic) strategy, that does not rely on an inner-outer iteration scheme. Numerical experiments show that our new solvers are competitive with other state-of-the-art solvers for low-rank problems, and deliver reconstructions of increased quality with respect to other classical Krylov methods.

This paper derives a new class of adaptive regularization parameter choice strategies that can be effectively and efficiently applied when regularizing large-scale linear inverse problems by combining standard Tikhonov regularization and projection onto Krylov subspaces of increasing dimension (computed by the Golub-Kahan bidiagonalization algorithm). The success of this regularization approach heavily depends on the accurate tuning of two parameters (namely, the Tikhonov parameter and the dimension of the projection subspace): these are simultaneously set using new strategies that can be regarded as special instances of bilevel optimization methods, which are solved by using a new paradigm that interlaces the iterations performed to project the Tikhonov problem (lower-level problem) with those performed to apply a given parameter choice rule (higher-level problem). The discrepancy principle, the GCV, the quasi-optimality criterion, and Regi\'{n}ska criterion can all be adapted to work in this framework. The links between Gauss quadrature and Golub-Kahan bidiagonalization are exploited to prove convergence results for the discrepancy principle, and to give insight into the behavior of the other considered regularization parameter choice rules. Several numerical tests modeling inverse problems in imaging show that the new parameter choice strategies lead to regularization methods that are reliable, and intrinsically simpler and cheaper than other strategies already available in the literature.

GMRES is one of the most popular iterative methods for the solution of large linear systems of equations that arise from the discretization of linear well-posed problems, such as Dirichlet boundary value problems for elliptic partial differential equations. The method is also applied to iteratively solve linear systems of equations that are obtained by discretizing linear ill-posed problems, such as many inverse problems. However, GMRES does not always perform well when applied to the latter kind of problems. This paper seeks to shed some light on reasons for the poor performance of GMRES in certain situations, and discusses some remedies based on specific kinds of preconditioning. The standard implementation of GMRES is based on the Arnoldi process, which also can be used to define a solution subspace for Tikhonov or TSVD regularization, giving rise to the Arnoldi-Tikhonov and Arnoldi-TSVD methods, respectively. The performance of the GMRES, the Arnoldi-Tikhonov, and the Arnoldi-TSVD methods is discussed. Numerical examples illustrate properties of these methods.

In this paper we develop flexible Krylov methods for efficiently computing regularized solutions to large-scale linear inverse problems with an $\ell_2$ fit-to-data term and an $\ell_p$ penalization term, for $p\geq 1$. First we approximate the $p$-norm penalization term as a sequence of $2$-norm penalization terms using adaptive regularization matrices in an iterative reweighted norm fashion, and then we exploit flexible preconditioning techniques to efficiently incorporate the weight updates. To handle general (non-square) $\ell_p$-regularized least-squares problems, we introduce a flexible Golub-Kahan approach and exploit it within a Krylov-Tikhonov hybrid framework. The key benefits of our approach compared to existing optimization methods for $\ell_p$ regularization are that efficient projection methods replace inner-outer schemes and that expensive regularization parameter selection techniques can be avoided. Theoretical insights are provided, and numerical results from image deblurring and tomographic reconstruction illustrate the benefits of this approach, compared to well-established methods. Furthermore, we show that our approach for $p=1$ can be used to efficiently compute solutions that are sparse with respect to some transformations.

This paper describes a new MATLAB software package of iterative regularization methods and test problems for large-scale linear inverse problems. The software package, called IR Tools, serves two related purposes: we provide implementations of a range of iterative solvers, including several recently proposed methods that are not available elsewhere, and we provide a set of large-scale test problems in the form of discretizations of 2D linear inverse problems. The solvers include iterative regularization methods where the regularization is due to the semi-convergence of the iterations, Tikhonov-type formulations where the regularization is explicitly formulated in the form of a regularization term, and methods that can impose bound constraints on the computed solutions. All the iterative methods are implemented in a very flexible fashion that allows the problem's coefficient matrix to be available as a (sparse) matrix, a function handle, or an object. The most basic call to all of the various iterative methods requires only this matrix and the right hand side vector; if the method uses any special stopping criteria, regularization parameters, etc., then default values are set automatically by the code. Moreover, through the use of an optional input structure, the user can also have full control of any of the algorithm parameters. The test problems represent realistic large-scale problems found in image reconstruction and several other applications. Numerical examples illustrate the various algorithms and test problems available in this package.

This paper introduces and analyzes an original class of Krylov subspace methods that provide an efficient alternative to many well-known conjugate-gradient-like (CG-like) Krylov solvers for square nonsymmetric linear systems arising from discretizations of inverse ill-posed problems. The main idea underlying the new methods is to consider some rank-deficient approximations of the transpose of the system matrix, obtained by running the (transpose-free) Arnoldi algorithm, and then apply some Krylov solvers to a formally right-preconditioned system of equations. Theoretical insight is given, and many numerical tests show that the new solvers outperform classical Arnoldi-based or CG-like methods in a variety of situations.

Constrained least squares problems arise in a variety of applications, and many iterative methods are already available to compute their solutions. This paper proposes a new efficient approach to solve nonnegative linear least squares problems. The associated KKT conditions are leveraged to form an adaptively preconditioned linear system, which is then solved by a flexible Krylov subspace method. The new method can be easily applied to image reconstruction problems affected by both Gaussian and Poisson noise, where the components of the solution represent nonnegative intensities. {Theoretical insight is given, and} numerical experiments and comparisons are displayed in order to validate the new method, which delivers results of equal or better quality than many state-of-the-art methods for nonnegative least squares solvers, with a significant speedup.

This paper introduces a new strategy for setting the regularization parameter when solving large-scale discrete ill-posed linear problems by means of the Arnoldi-Tikhonov method. This new rule is essentially based on the discrepancy principle, although no initial knowledge of the norm of the error that affects the right-hand side is assumed; an increasingly more accurate approximation of this quantity is recovered during the Arnoldi algorithm. Some theoretical estimates are derived in order to motivate our approach. Many numerical experiments, performed on classical test problems as well as image deblurring are presented.