27 research outputs found
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
Application of the principle and unbiased predictive risk estimator for determining the regularization parameter in 3D focusing gravity inversion
The 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 inversion of gravity data
Sparse inversion of gravity data based on -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 -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 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 . 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 , for data
sets of size . 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 . Moreover,
it is sufficient to use projected spaces of size when is large, , to reconstruct volumes with . 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
Regularization Parameter Estimation for Underdetermined problems by the principle with application to focusing gravity inversion
The -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 -distribution with
degrees of freedom, model matrix and regularizer
. It is proved that the result holds also for when .
A Newton root-finding algorithm is used to find the regularization parameter
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 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 -principle, moreover the L-curve and MDP are
generally outperformed. Furthermore, the MS is of general use for the
-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
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
-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
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 -D gravity inverse problem using a randomized generalized singular value decomposition
We present a fast algorithm for the total variation regularization of the
-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
IGUG: A MATLAB package for D inversion of gravity data using graph theory
The open source MATLAB package IGUG for 3D inversion of gravity data is
presented. It is based on methodology that was introduced by Bijani et al
(2015), in which a homogeneous subsurface body is modeled by an ensemble of
simple point masses. The model parameters are the Cartesian coordinates of the
point masses and their total mass. Associating the point masses to the vertices
of a weighted full graph with weights computed by the Euclidean pairwise
distances separating vertices, Kruskal's algorithm is used to find the minimum
spanning tree of the graph. Stabilization is achieved using an equidistance
function that restricts the spatial distribution of the masses and favors a
homogeneous subsurface structure. A regularization parameter is
introduced to balance the two terms of the objective function, and reasonable
physically-relevant bound constraints are imposed on the model parameters. Then
the bound constrained objective function is solved using a genetic algorithm. A
new diagnostic approach is presented for determining a suitable choice for
, requiring a limited number of solutions for a small set of
without using the L-curve as suggested by Bijani et al. Simulations for
synthetic examples demonstrate the efficiency and effectiveness of the
implementation of the algorithm. It is verified that the constraints on the
model parameters are not restrictive. Included in the package is the script
GMD.m which is used for generating synthetic data and for putting measurement
data in the format required for the inversion implemented by IGUG.m The script
Diagnostic_Results.m is included for analyzing and visualizing the results. The
software can be used to verify the simulations and the analysis of real data
that is presented here, The real data set uses gravity data from the Mobrun ore
body, north east of Noranda, Quebec, Canada