Rosemary A. Renaut

Rosemary A. Renaut, Saeed Vatankhah, Vahid E. Ardestani

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á, 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.

Saeed Vatankhah, Vahid E. Ardestani, Rosemary A. Renaut

The $χ^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.

Saeed Vatankhah, Rosemary A. Renaut, Vahid E. Ardestani

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.

Rosemary A. Renaut, Michael Horst, Yang Wang et al.

The solution, $x$, of the linear system of equations $A x\approx b$ arising from the discretization of an ill-posed integral equation with a square integrable kernel $H(s,t)$ is considered. The Tikhonov regularized solution $ x(λ)$ is found as the minimizer of $J( x)=\{ \|A x - b\|_2^2 + λ^2 \|L x\|_2^2\}$. $ x(λ)$ depends on regularization parameter $λ$ that trades off the data fidelity, and on the smoothing norm determined by $L$. Here we consider the case where $L$ is diagonal and invertible, and employ the Galerkin method to provide the relationship between the singular value expansion and the singular value decomposition for square integrable kernels. The resulting approximation of the integral equation permits examination of the properties of the regularized solution $ x(λ)$ independent of the sample size of the data. We prove that estimation of the regularization parameter can be obtained by consistently down sampling the data and the system matrix, leading to solutions of coarse to fine grained resolution. Hence, the estimate of $λ$ for a large problem may be found by downsampling to a smaller problem, or to a set of smaller problems, effectively moving the costly estimate of the regularization parameter to the coarse representation of the problem. Moreover, the full singular value decomposition for the fine scale system is replaced by a number of dominant terms which is determined from the coarse resolution system, again reducing the computational cost. Numerical results illustrate the theory and demonstrate the practicality of the approach for regularization parameter estimation using generalized cross validation, unbiased predictive risk estimation and the discrepancy principle applied for both the system of equations, and the augmented system of equations.

Saeed Vatankhah, Rosemary A. Renaut, Vahid E. Ardestani

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.