Regularization Parameter Estimation for Underdetermined problems by the χ2\chi^2 principle with application to 2D2D focusing gravity inversion

The $\chi^2$-principle generalizes the Morozov discrepancy principle (MDP) to the augmented residual of the Tikhonov regularized least squares problem. Weighting of the data fidelity by a known Gaussian noise distribution on the measured data, when the regularization term is weighted by unknown inverse covariance information on the model parameters, the minimum of the Tikhonov functional is a random variable following a $\chi^2$-distribution with $m+p-n$ degrees of freedom, model matrix $G:$ $m \times n$ and regularizer $L:$ $p\times n$. It is proved that the result holds also for $m<n$ when $m+p\ge n$. A Newton root-finding algorithm is used to find the regularization parameter $\alpha$ which yields the optimal inverse covariance weighting in the case of a white noise assumption on the mapped model data. It is implemented for small-scale problems using the generalized singular value decomposition. Numerical results verify the algorithm for the case of regularizers approximating zero to second order derivative approximations, contrasted with the methods of generalized cross validation and unbiased predictive risk estimation. The inversion of underdetermined $2D$ focusing gravity data produces models with non-smooth properties, for which typical implementations in this field use the iterative minimum support (MS) stabilizer and both regularizer and regularizing parameter are updated each iteration. For a simulated data set with noise, the regularization parameter estimation methods for underdetermined data sets are used in this iterative framework, also contrasted with the L-curve and MDP. Experiments demonstrate efficiency and robustness of the $\chi^2$-principle, moreover the L-curve and MDP are generally outperformed. Furthermore, the MS is of general use for the $\chi^2$-principle when implemented without the knowledge of a mean value of the model

Convergence of Regularization Parameters for Solutions Using the Filtered Truncated Singular Value Decomposition

The truncated singular value decomposition may be used to find the solution of linear discrete ill-posed problems in conjunction with Tikhonov regularization and requires the estimation of a regularization parameter that balances between the sizes of the fit to data function and the regularization term. The unbiased predictive risk estimator is one suggested method for finding the regularization parameter when the noise in the measurements is normally distributed with known variance. In this paper we provide an algorithm using the unbiased predictive risk estimator that automatically finds both the regularization parameter and the number of terms to use from the singular value decomposition. Underlying the algorithm is a new result that proves that the regularization parameter converges with the number of terms from the singular value decomposition. For the analysis it is sufficient to assume that the discrete Picard condition is satisfied for exact data and that noise completely contaminates the measured data coefficients for a sufficiently large number of terms, dependent on both the noise level and the degree of ill-posedness of the system. A lower bound for the regularization parameter is provided leading to a computationally efficient algorithm. Supporting results are compared with those obtained using the method of generalized cross validation. Simulations for two-dimensional examples verify the theoretical analysis and the effectiveness of the algorithm for increasing noise levels, and demonstrate that the relative reconstruction errors obtained using the truncated singular value decomposition are less than those obtained using the singular value decomposition. This is a pre-print of an article published in BIT Numerical Mathematics. The final authenticated version is available online at: https://doi.org/10.1007%2Fs10543-019-00762-7

Application of the χ2\chi^2 principle and unbiased predictive risk estimator for determining the regularization parameter in 3D focusing gravity inversion

The $\chi^2$ principle and the unbiased predictive risk estimator are used to determine optimal regularization parameters in the context of 3D focusing gravity inversion with the minimum support stabilizer. At each iteration of the focusing inversion the minimum support stabilizer is determined and then the fidelity term is updated using the standard form transformation. Solution of the resulting Tikhonov functional is found efficiently using the singular value decomposition of the transformed model matrix, which also provides for efficient determination of the updated regularization parameter each step. Experimental 3D simulations using synthetic data of a dipping dike and a cube anomaly demonstrate that both parameter estimation techniques outperform the Morozov discrepancy principle for determining the regularization parameter. Smaller relative errors of the reconstructed models are obtained with fewer iterations. Data acquired over the Gotvand dam site in the south-west of Iran are used to validate use of the methods for inversion of practical data and provide good estimates of anomalous structures within the subsurface

3-D Projected L1 L_{1} inversion of gravity data

Sparse inversion of gravity data based on $L_1$-norm regularization is discussed. An iteratively reweighted least squares algorithm is used to solve the problem. At each iteration the solution of a linear system of equations and the determination of a suitable regularization parameter are considered. The LSQR iteration is used to project the system of equations onto a smaller subspace that inherits the ill-conditioning of the full space problem. We show that the gravity kernel is only mildly to moderately ill-conditioned. Thus, while the dominant spectrum of the projected problem accurately approximates the dominant spectrum of the full space problem, the entire spectrum of the projected problem inherits the ill-conditioning of the full problem. Consequently, determining the regularization parameter based on the entire spectrum of the projected problem necessarily over compensates for the non-dominant portion of the spectrum and leads to inaccurate approximations for the full-space solution. In contrast, finding the regularization parameter using a truncated singular space of the projected operator is efficient and effective. Simulations for synthetic examples with noise demonstrate the approach using the method of unbiased predictive risk estimation for the truncated projected spectrum. The method is used on gravity data from the Mobrun ore body, northeast of Noranda, Quebec, Canada. The $3$-D reconstructed model is in agreement with known drill-hole information

A fast methodology for large-scale focusing inversion of gravity and magnetic data using the structured model matrix and the 2D2D fast Fourier transform

Focusing inversion of potential field data for the recovery of sparse subsurface structures from surface measurement data on a uniform grid is discussed. For the uniform grid the model sensitivity matrices exhibit block Toeplitz Toeplitz block structure, by blocks for each depth layer of the subsurface. Then, through embedding in circulant matrices, all forward operations with the sensitivity matrix, or its transpose, are realized using the fast two dimensional Fourier transform. Simulations demonstrate that this fast inversion algorithm can be implemented on standard desktop computers with sufficient memory for storage of volumes up to size $n \approx 1M$. The linear systems of equations arising in the focusing inversion algorithm are solved using either Golub Kahan bidiagonalization or randomized singular value decomposition algorithms in which all matrix operations with the sensitivity matrix are implemented using the fast Fourier transform. These two algorithms are contrasted for efficiency for large-scale problems with respect to the sizes of the projected subspaces adopted for the solutions of the linear systems. The presented results confirm earlier studies that the randomized algorithms are to be preferred for the inversion of gravity data, and that it is sufficient to use projected spaces of size approximately $m/8$, for data sets of size $m$. In contrast, the Golub Kahan bidiagonalization leads to more efficient implementations for the inversion of magnetic data sets, and it is again sufficient to use projected spaces of size approximately $m/8$. Moreover, it is sufficient to use projected spaces of size $m/20$ when $m$ is large, $m \approx 50000$, to reconstruct volumes with $n \approx 1M$. Simulations support the presented conclusions and are verified on the inversion of a practical magnetic data set that is obtained over the Wuskwatim Lake region in Manitoba, Canada.Comment: 39 pages. 16 figure

Hybrid and iteratively reweighted regularization by unbiased predictive risk and weighted GCV for projected systems

Tikhonov regularization for projected solutions of large-scale ill-posed problems is considered. The Golub-Kahan iterative bidiagonalization is used to project the problem onto a subspace and regularization then applied to find a subspace approximation to the full problem. Determination of the regularization parameter for the projected problem by unbiased predictive risk estimation, generalized cross validation and discrepancy principle techniques is investigated. It is shown that the regularized parameter obtained by the unbiased predictive risk estimator can provide a good estimate for that to be used for a full problem which is moderately to severely ill-posed. A similar analysis provides the weight parameter for the weighted generalized cross validation such that the approach is also useful in these cases, and also explains why the generalized cross validation without weighting is not always useful. All results are independent of whether systems are over or underdetermined. Numerical simulations for standard one dimensional test problems and two dimensional data, for both image restoration and tomographic image reconstruction, support the analysis and validate the techniques. The size of the projected problem is found using an extension of a noise revealing function for the projected problem Hn\u etynkov\'a, Ple\u singer, and Strako\u s, [\textit{BIT Numerical Mathematics} {\bf 49} (2009), 4 pp. 669-696.]. Furthermore, an iteratively reweighted regularization approach for edge preserving regularization is extended for projected systems, providing stabilization of the solutions of the projected systems and reducing dependence on the determination of the size of the projected subspace

Automatic estimation of the regularization parameter in 2-D focusing gravity inversion: an application to the Safo manganese mine in northwest of Iran

We investigate the use of Tikhonov regularization with the minimum support stabilizer for underdetermined 2-D inversion of gravity data. This stabilizer produces models with non-smooth properties which is useful for identifying geologic structures with sharp boundaries. A very important aspect of using Tikhonov regularization is the choice of the regularization parameter that controls the trade off between the data fidelity and the stabilizing functional. The L-curve and generalized cross validation techniques, which only require the relative sizes of the uncertainties in the observations are considered. Both criteria are applied in an iterative process for which at each iteration a value for regularization parameter is estimated. Suitable values for the regularization parameter are successfully determined in both cases for synthetic but practically relevant examples. Whenever the geologic situation permits, it is easier and more efficient to model the subsurface with a 2-D algorithm, rather than to apply a full 3-D approach. Then, because the problem is not large it is appropriate to use the generalized singular value decomposition for solving the problem efficiently. The method is applied on a profile of gravity data acquired over the Safo mining camp in Maku-Iran, which is well known for manganese ores. The presented results demonstrate success in reconstructing the geometry and density distribution of the subsurface source

Total variation regularization of the 33-D gravity inverse problem using a randomized generalized singular value decomposition

We present a fast algorithm for the total variation regularization of the $3$-D gravity inverse problem. Through imposition of the total variation regularization, subsurface structures presenting with sharp discontinuities are preserved better than when using a conventional minimum-structure inversion. The associated problem formulation for the regularization is non linear but can be solved using an iteratively reweighted least squares algorithm. For small scale problems the regularized least squares problem at each iteration can be solved using the generalized singular value decomposition. This is not feasible for large scale problems. Instead we introduce the use of a randomized generalized singular value decomposition in order to reduce the dimensions of the problem and provide an effective and efficient solution technique. For further efficiency an alternating direction algorithm is used to implement the total variation weighting operator within the iteratively reweighted least squares algorithm. Presented results for synthetic examples demonstrate that the novel randomized decomposition provides good accuracy for reduced computational and memory demands as compared to use of classical approaches

A fast algorithm for regularized focused 3-D inversion of gravity data using the randomized SVD

A fast algorithm for solving the under-determined 3-D linear gravity inverse problem based on the randomized singular value decomposition (RSVD) is developed. The algorithm combines an iteratively reweighted approach for $L_1$-norm regularization with the RSVD methodology in which the large scale linear system at each iteration is replaced with a much smaller linear system. Although the optimal choice for the low rank approximation of the system matrix with m rows is q=m, acceptable results are achievable with q<<m. In contrast to the use of the LSQR algorithm for the solution of the linear systems at each iteration, the singular values generated using the RSVD yield a good approximation of the dominant singular values of the large scale system matrix. The regularization parameter found for the small system at each iteration is thus dependent on the dominant singular values of the large scale system matrix and appropriately regularizes the dominant singular space of the large scale problem. The results achieved are comparable with those obtained using the LSQR algorithm for solving each linear system, but are obtained at reduced computational cost. The method has been tested on synthetic models along with the real gravity data from the Morro do Engenho complex from central Brazil

Model Error Correction for Linear Methods of Reversible Radioligand Binding Measurements in PET Studies

Graphical analysis methods are widely used in positron emission tomography quantification because of their simplicity and model independence. But they may, particularly for reversible kinetics, lead to bias in the estimated parameters. The source of the bias is commonly attributed to noise in the data. Assuming a two-tissue compartmental model, we investigate the bias that originates from model error. This bias is an intrinsic property of the simplified linear models used for limited scan durations, and it is exaggerated by random noise and numerical quadrature error. Conditions are derived under which Logan's graphical method either over- or under-estimates the distribution volume in the noise-free case. The bias caused by model error is quantified analytically. The presented analysis shows that the bias of graphical methods is inversely proportional to the dissociation rate. Furthermore, visual examination of the linearity of the Logan plot is not sufficient for guaranteeing that equilibrium has been reached. A new model which retains the elegant properties of graphical analysis methods is presented, along with a numerical algorithm for its solution. We perform simulations with the fibrillar amyloid-beta radioligand [11C] benzothiazole-aniline using published data from the University of Pittsburgh and Rotterdam groups. The results show that the proposed method significantly reduces the bias due to model error. Moreover, the results for data acquired over a 70 minutes scan duration are at least as good as those obtained using existing methods for data acquired over a 90 minutes scan duration