Improving the use of the randomized singular value decomposition for the inversion of gravity and magnetic data

TL;DR

This paper enhances large-scale gravity and magnetic data inversion by employing randomized singular value decomposition, improving computational efficiency and approximation quality, especially for magnetic data, through power iteration techniques.

Contribution

It introduces a comprehensive methodology using randomized SVD for large-scale gravity and magnetic data inversion, highlighting differences and improvements for magnetic problems.

Findings

Randomized SVD effectively accelerates large-scale inversion computations.

Magnetic inversion requires larger rank approximations than gravity inversion.

Power iteration improves approximation quality for magnetic data inversion.

Abstract

The large-scale focusing inversion of gravity and magnetic potential field data using $L_1$-norm regularization is considered. The use of the randomized singular value decomposition methodology facilitates tackling the computational challenge that arises in the solution of these large-scale inverse problems. As such the powerful randomized singular value decomposition is used for the numerical solution of all linear systems required in the algorithm. A comprehensive comparison of the developed methodology for the inversion of magnetic and gravity data is presented. These results indicate that there is generally an important difference between the gravity and magnetic inversion problems. Specifically, the randomized singular value decomposition is dependent on the generation of a rank $q$ approximation to the underlying model matrix, and the results demonstrate that $q$ needs to be larger, for equivalent problem sizes, for the magnetic problem as compared to the gravity problem. Without a relatively large $q$ the dominant singular values of the magnetic model matrix are not well-approximated. The comparison also shows how the use of the power iteration embedded within the randomized algorithm is used to improve the quality of the resulting dominant subspace approximation, especially in magnetic inversion, yielding acceptable approximations for smaller choices of $q$. The price to pay is the trade-off between approximation accuracy and computational cost. The algorithm is applied for the inversion of magnetic data obtained over a portion of the Wuskwatim Lake region in Manitoba, Canada

Tables

Table 1. Table 1: Results of the inversion algorithms applied on gravity data of Fig.  2 .

Method	$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$R{G}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
Hybrid-FSVD	$-$	$61712$	$32.2$	$0.7232$	$-$	$12$	$811.8$	$31.5$
Hybrid-RSVD	$100$	$24074$	$44.5$	$0.7550$	$0.0521$	$14$	$935.8$	$37.2$
	$150$	$29524$	$37.4$	$0.7441$	$0.0497$	$12$	$928.7$	$32.7$
	$200$	$33726$	$34.5$	$0.7221$	$0.0448$	$12$	$905.6$	$32.3$
	$300$	$40360$	$33.4$	$0.7153$	$0.0257$	$12$	$901.4$	$33.0$
	$500$	$49504$	$29.5$	$0.7077$	$0.0143$	$12$	$801.9$	$35.0$
	$700$	$56105$	$31.1$	$0.7061$	$0.0082$	$12$	$908.1$	$35.9$
	$900$	$61712$	$32.2$	$0.7232$	$2.8510e^{-14}$	$12$	$811.8$	$37.9$

Table 2. Table 2: Results of the inversion algorithms applied on magnetic data of Fig.  2 .

Method	$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$R{G}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
Hybrid-FSVD	$-$	$21086$	$4106.4$	$0.8454$	$-$	$8$	$898.6$	$20.7$
Hybrid-RSVD	$100$	$9866$	$4372.3$	$1.0651$	$0.3735$	$50$	$2880.6$	$125.0$
	$150$	$11247$	$3089.6$	$0.9904$	$0.4143$	$50$	$2592.5$	$127.9$
	$200$	$12415$	$3018.0$	$0.9050$	$0.3701$	$50$	$1673.7$	$134.9$
	$300$	$14371$	$3030.9$	$0.8606$	$0.2948$	$50$	$1120.5$	$135.0$
	$500$	$17011$	$3861.2$	$0.8426$	$0.0453$	$9$	$917.6$	$26.9$
	$700$	$19071$	$3854.1$	$0.8470$	$0.0268$	$8$	$918.2$	$24.8$
	$900$	$21086$	$4106.4$	$0.8454$	$2.8165e^{-14}$	$8$	$898.6$	$26.1$

Table 3. Table 3: Results of the Hybrid-RSVD algorithm via power iterations applied on gravity data of Fig.  2 .

$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$R{G}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
$100$	$20629$	$54.3$	$0.7018$	$0.0177$	$11$	$913.4$	$30.4$
$150$	$25862$	$45.5$	$0.7288$	$0.0152$	$11$	$909.6$	$29.5$
$200$	$30218$	$40.8$	$0.7185$	$0.0129$	$11$	$904.7$	$30.9$
$300$	$37279$	$34.8$	$0.7188$	$0.0101$	$11$	$898.1$	$31.8$
$500$	$47447$	$29.6$	$0.7099$	$0.0063$	$12$	$851.0$	$38.0$
$700$	$54892$	$29.0$	$0.7093$	$0.0047$	$12$	$893.6$	$39.6$
$900$	$61712$	$32.2$	$0.7232$	$1.5534e^{-15}$	$12$	$811.8$	$41.3$

Table 4. Table 4: Results of the Hybrid-RSVD algorithm via power iterations applied on magnetic data of Fig.  2 .

$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$R{G}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
$100$	$8375$	$3127.0$	$0.9187$	$0.2880$	$50$	$1310.6$	$128.3$
$150$	$9920$	$3206.5$	$0.8464$	$0.2294$	$50$	$955.5$	$128.7$
$200$	$11198$	$5098.5$	$0.8235$	$0.0542$	$12$	$915.4$	$33.2$
$300$	$13235$	$4447.1$	$0.8310$	$0.0323$	$10$	$848.9$	$28.7$
$500$	$16152$	$3719.6$	$0.8421$	$0.0269$	$9$	$814.5$	$29.0$
$700$	$18528$	$4036.8$	$0.8405$	$0.0165$	$8$	$887.5$	$26.7$
$900$	$21086$	$4106.4$	$0.8454$	$2.0710e^{-15}$	$8$	$898.6$	$28.7$

Table 5. Table 5: Results of the inversion algorithms applied on the gravity data of Fig.  14 .

Method	$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
Hybrid-RSVD	$1000$	$40276$	$38.4$	$0.7389$	$21$	$8108.5$	$376.8$
	$1500$	$48847$	$44.6$	$0.7168$	$20$	$7706.2$	$524.3$
	$2500$	$61439$	$39.1$	$0.7111$	$19$	$7893.3$	$946.6$
	$4200$	$76008$	$37.9$	$0.7013$	$20$	$6943.7$	$2114.0$
Hybrid-RSVD	$1000$	$34324$	$38.1$	$0.6927$	$21$	$7915.9$	$617.0$
with power iterations	$1500$	$43000$	$36.7$	$0.6976$	$20$	$7985.6$	$886.3$
	$2500$	$56392$	$40.0$	$0.6942$	$20$	$7561.5$	$1775.4$
	$4200$	$72451$	$38.2$	$0.6986$	$20$	$7989.2$	$4178.8$

Table 6. Table 6: Results of the inversion algorithms applied on the magnetic data of Fig.  14 .

Method	$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$R{E}^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
Hybrid-RSVD	$1000$	$12600$	$10774.9$	$1.2110$	$50$	$26292.4$	$887.9$
	$1500$	$14564$	$8947.8$	$1.1720$	$50$	$20789.4$	$1312.4$
	$2500$	$17376$	$9253.8$	$1.1063$	$50$	$11891.4$	$2520.0$
	$4200$	$20695$	$8778.1$	$1.0975$	$13$	$7750.3$	$1435.9$
Hybrid-RSVD	$1000$	$10705$	$11128.1$	$1.1315$	$50$	$13480.7$	$1459.1$
with power iterations	$1500$	$12744$	$8992.1$	$1.0910$	$50$	$9570.6$	$2263.9$
	$2500$	$15837$	$8425.2$	$1.0947$	$14$	$7863.0$	$1260.4$
	$4200$	$19518$	$8665.2$	$1.1080$	$12$	$7299.8$	$2685.4$

Table 7. Table 7: Results of the inversion algorithms applied on the magnetic data of Fig.  19 .

Method	$q$	${\alpha }^{\left(1\right)}$	${\alpha }^{\left(K\right)}$	$K$	${\chi }^2$	Time (s)
Hybrid-RSVD	$2100$	$262575$	$3630.2$	$15$	$4007.1$	$559.7$
Hybrid-RSVD with power iterations	$1100$	$215644$	$2906.7$	$16$	$4089.7$	$464.0$

Equations

$${\mathbf{d}}_{\mathrm{obs}}=G\mathbf{m},$$

$$\mathbf{m}=\text{arg}\underset{\mathbf{m}}{\min }\{P^{\alpha }(\mathbf{m})\}=\text{arg}\underset{\mathbf{m}}{\min }\{\Vert W_{\mathbf{d}}(G\mathbf{m}-{\mathbf{d}}_{\mathrm{obs}}){\Vert }_2^2+{\alpha }^2\Vert W(\mathbf{m}-{\mathbf{m}}_{\mathrm{apr}}){\Vert }_2^2\}.$$

$$P^{\alpha }(\mathbf{h})=\Vert \widetilde{\stackrel{~}{G}}\mathbf{h}-\tilde{\mathbf{r}}{\Vert }_2^2+{\alpha }^2\Vert \mathbf{h}{\Vert }_2^2,$$

$$\mathbf{h}(\alpha )=({\widetilde{\stackrel{~}{G}}}^T\widetilde{\stackrel{~}{G}}+{\alpha }^2{I}_n{)}^{-1}{\widetilde{\stackrel{~}{G}}}^T\tilde{\mathbf{r}},$$

$$\mathbf{m}(\alpha )={\mathbf{m}}_{\mathrm{apr}}+W^{-1}\mathbf{h}(\alpha ).$$

$$U(\alpha )=\sum _{i=1}^{m^{\ast }}{\left(\frac{1}{{\sigma }_i^2{\alpha }^{-2}+1}\right)}^2{\left({\mathbf{u}}_i^T\tilde{\mathbf{r}}\right)}^2+2\left(\sum _{i=1}^{m^{\ast }}\frac{{\sigma }_i^2}{{\sigma }_i^2+{\alpha }^2}\right)-m^{\ast }.$$

$${\widetilde{\stackrel{~}{G}}}^{(s)}=\widetilde{\stackrel{~}{G}}({\widetilde{\stackrel{~}{G}}}^T\widetilde{\stackrel{~}{G}}{)}^{(s)},$$

$$E\Vert \widetilde{\stackrel{~}{G}}-{\widetilde{\stackrel{~}{G}}}_q\Vert \le {\left[1+\sqrt{\frac{q}{p-1}}+\frac{e\sqrt{q+p}}{p}.\sqrt{min(m,n)-q}\right]}^{\frac{1}{2s+1}}{\sigma }_{q+1},$$

$$SNR=20{\log }_{10}\frac{\Vert {\mathbf{d}}_{\mathrm{exact}}{\Vert }_2}{\Vert {\mathbf{d}}_{\mathrm{obs}}-{\mathbf{d}}_{\mathrm{exact}}{\Vert }_2},$$

$$({\chi }^2{)}^{(k)}=\Vert W_{\mathbf{d}}({\mathbf{d}}_{\mathrm{obs}}-{\mathbf{d}}_{\mathrm{pre}}^{(k)}){\Vert }_2^2,$$

$${\alpha }^{(1)}={\left(\frac{n}{m}\right)}^{3.5}\frac{{\sigma }_1}{\mathrm{mean}({\sigma }_i)},$$

$$R{E}^{(k)}=\frac{\Vert {\mathbf{m}}_{\mathrm{exact}}-{\mathbf{m}}^{(k)}{\Vert }_2}{\Vert {\mathbf{m}}_{\mathrm{exact}}{\Vert }_2}$$

$$R{G}^{(k)}=\frac{\Vert \widetilde{\stackrel{~}{G}}-{\widetilde{\stackrel{~}{G}}}_q^{(k)}{\Vert }_2}{\Vert \widetilde{\stackrel{~}{G}}{\Vert }_2}$$

Full text

{summary}

The large-scale focusing inversion of gravity and magnetic potential field data using $L_{1}$-norm regularization is considered. The use of the randomized singular value decomposition methodology facilitates tackling the computational challenge that arises in the solution of these large-scale inverse problems. As such the powerful randomized singular value decomposition is used for the numerical solution of all linear systems required in the algorithm. A comprehensive comparison of the developed methodology for the inversion of magnetic and gravity data is presented. These results indicate that there is generally an important difference between the gravity and magnetic inversion problems. Specifically, the randomized singular value decomposition is dependent on the generation of a rank $q$ approximation to the underlying model matrix, and the results demonstrate that $q$ needs to be larger, for equivalent problem sizes, for the magnetic problem as compared to the gravity problem. Without a relatively large $q$ the dominant singular values of the magnetic model matrix are not well-approximated. The comparison also shows how the use of the power iteration embedded within the randomized algorithm is used to improve the quality of the resulting dominant subspace approximation, especially in magnetic inversion, yielding acceptable approximations for smaller choices of $q$. The price to pay is the trade-off between approximation accuracy and computational cost. The algorithm is applied for the inversion of magnetic data obtained over a portion of the Wuskwatim Lake region in Manitoba, Canada

1 Introduction

Potential field surveys, gravity and magnetic, have been used for many years for a wide range of studies including oil and gas exploration, mining applications, and mapping basement topography [\citenameNabighian et al. 2005, \citenameBlakely 1996]. The inversion of acquired data is one of the important steps in the interpretation process [\citenameLi & Oldenburg 1998, \citenamePortniaguine & Zhdanov 1999, \citenameBoulanger & Chouteau 2001, \citenameSilva & Barbosa 2006, \citenameFarquharson 2008, \citenameLiu et al. 2013, \citenameLiu et al. 2018]. The problem is ill-posed and then, generally, the solution is obtained via the minimization of a global objective function that consists of two terms, the data misfit term and the stabilizing, or regularizing, term. These two terms are balanced by a scalar regularization parameter that weights the contribution of the stabilizing term to the solution. Extensive background on the modeling and the solution of the regularized objective function is provided in the literature, e.g. [\citenameLi & Oldenburg 1998, \citenamePortniaguine & Zhdanov 1999, \citenameVatankhah et al. 2015]. While the concept of regularization is well-known, the numerical solution of a large-scale problem, in conjunction with effective and efficient estimation of a suitable regularization parameter, continues to be computationally challenging. For large-scale problems, the most effective and widely-used strategy is to transform the problem from the original large space to a much smaller subspace. The resulting subspace solution can then be projected back to the original full space, under reasonable assumptions that the subspace problem sufficiently captures the characteristics of the full space problem. For example the LSQR algorithm, which is based on the Golub-Kahan bidiagonalization of the model matrix, is frequently used for the inversion of geophysical data. Using the LSQR algorithm the large-scale problem is projected onto a Krylov subspace of smaller dimension and the subspace solution is then obtained relatively efficiently using a standard factorization such as the singular value decomposition (SVD), [\citenamePaige & Saunders 1982a, \citenamePaige & Saunders 1982b, \citenameKilmer & O’Leary 2001, \citenameChung et al. 2008, \citenameRenaut et al. 2017, \citenameVatankhah et al. 2017]. Still the need to solve ever larger problems so as to provide greater resolution of the subsurface structures, while also automatically estimating a suitable regularization parameter, presents a challenge computationally. Even with the continually-increasing computational power and memory that is available, it is sometimes impossible, or computationally prohibitive to obtain effective solutions with these traditional computational algorithms.

Recently, the powerful concept of randomization has been introduced as an alternative strategy for dealing with large-scale inverse problems [\citenameHalko et al. 2011, \citenameXiang & Zou 2013, \citenameVoronin et al. 2015, \citenameXiang & Zou 2015, \citenameWei et al. 2016, \citenameVatankhah et al. 2018a, \citenameVatankhah et al. 2018b]. The fundamental idea is that some amount of randomness can be employed to obtain a matrix of smaller dimension that effectively captures the essential spectral properties of the original model matrix, and thus provides an approximately optimal rank-$q$ randomized singular value decomposition (RSVD) of the original system. Effectively, the approach consists of the first stochastic step which finds a random but orthonormal matrix that samples rows and columns of a given matrix, and a second completely deterministic step which then finds the eigen-decomposition of the subsampled matrix and hence yields a rank-$q$ approximation of the original matrix.

For the $3$D inversion of gravity data, Vatankhah et al. \shortciteVRA:2018b developed a fast inversion methodology that combined an $L_{1}$-norm regularization strategy with the RSVD, for the generation of a focused image of the subsurface. For the under-determined problem in which there are $m$ data measurements to be used to find the subsurface structures on a volume with $n$ model parameters, $m\ll n$, their results indicate that acceptable results, nearly equivalent to those using the full SVD (FSVD) for the SVD calculations in the algorithm, are achievable using a target rank $q$ with $q\gtrsim(m/6)$.111The solutions are the same when $q=m$. Here we will denote the algorithms that use the RSVD and FSVD, respectively, for all steps in the $L_{1}$-norm regularization strategy the Hybrid-RSVD and Hybrid-FSVD algorithms. For the Hybrid-RSVD algorithm it is important to estimate an appropriate lower bound on $q$ in order that acceptable solutions are provided using an algorithm that is computationally efficient with respect to both time and memory.

The purpose of this discussion is to assess the application of the Hybrid-RSVD algorithm, as presented in Vatankhah et al. \shortciteVRA:2018b, for the inversion of magnetic potential data. Our results demonstrate that direct application of the Hybrid-RSVD algorithm does not lead to acceptable solutions. Rather, for a problem of equivalent size as given by the pair $(m,n)$, we find that $q$ must be larger; the lower bound $q\gtrsim m/6$ does not yield acceptable solutions. On the other hand, by introducing power iterations into the determination of the rank-$q$ approximation, improved results are achievable for smaller choices of $q$, for both magnetic and gravity inversion problems. Both cases are analyzed carefully and useful estimates for the appropriate choices of $q$ for both magnetic and gravity problems, when solved using power iterations to improve the approximate RSVD, are provided.

The remainder of the paper is organized as follows. In Section 2 we briefly describe the standard $L_{1}$-norm stabilized inversion algorithm. A short explanation of the methodology used for estimating the regularization parameter is presented in Section 2.0.1. Then, the concept of the RSVD is reviewed in Section 2.1, with the extension applying power iterations presented in Section 2.1.1. In Section 3 we show the results of using the Hybrid-RSVD algorithm on two different synthetic models. A small model in Section 3.1 is used for further analysis of the algorithm and the impact of using power iterations to obtain the RSVD is also demonstrated. To further examine the approach we consider a model with multiple bodies in Section 3.2. The inversion of magnetic data from a meta-sedimentary gneiss belt in Canada are used to illustrate the application of the presented methodology for real data in Section 4. Conclusions and discussions on future work are given in Section 5.

2 Inversion methodology

We discuss the general case for the linear inversion of potential field data in which the subsurface is discretized into a large number of cells of fixed size but with unknown physical properties. The unknown parameters of the cells are stacked in a vector $\mathbf{m}\in\mathcal{R}^{n}$, and measured potential field data are also stacked in a vector $\mathbf{d}_{\mathrm{obs}}\in\mathcal{R}^{m}$. The measurements are connected to the model parameters via the model matrix $G\in\mathcal{R}^{m\times n}$, $m\ll n$, which is the forward model operator which maps from model to data spaces, yielding the under-determined linear system

[TABLE]

Element $G_{ij}$ represents the effect of the unit model parameter at cell $j$ on the data at location $i$.

In the case of the inversion of magnetic data, vector $\mathbf{m}$ collects the values of the unknown susceptibilities of the cells and $\mathbf{d}_{\mathrm{obs}}$ collects the total magnetic field. For the gravity problem these are the cell densities and vertical component of gravity field, respectively. Model matrix $G$ depends on the problem. Rao $\&$ Babu \shortciteRaBa:91 provided a fast approach for computing the total magnetic field anomaly of a cube that is then used here to form the elements of model matrix $G$ for the magnetic problem. For gravity inversion, the elements of $G$ are computed using the formula developed by Haáz \shortciteHaaz, see for example Boulanger & Chouteau \shortciteBoCh:2001 for more details. The spectral properties of these two matrices impact the condition of the under-determined systems, and hence the performance of any algorithm that is used for data inversion.

Stabilization, or regularization, is required to find an acceptable solution of the ill-posed system given by (1); its solution is neither unique nor stable. Here we consider the $L_{1}$-norm stabilized solution of (1) as presented by Vatankhah et al. \shortciteVRA:2017. The approach also includes weighting of the data based on the knowledge, or estimate, of the standard deviation of the independent noise in the data, weighting for the depths of the cells, and the inclusion of prior knowledge on the solution using provided information which may come from geology, logging or previous geophysical surveys, or may be taken to be the vector of [math] values when no prior information is available. Overall, the formulation finds the minimum of the global objective function $P^{\alpha}(\mathbf{m})$ as given by

[TABLE]

The diagonal entries of the data weighting matrix $W_{\mathbf{d}}$ are estimates for the inverse of the standard deviations of the independent noise in the data, $\mathbf{m}_{\mathrm{apr}}$ encodes the prior knowledge of the solution, and stabilizing matrix $W$ is the product of three diagonal matrices, $W_{\mathrm{depth}}$ which is a depth-weighting matrix with diagonal entries $z^{-\beta}$ at depth $z$, $W_{\mathrm{h}}$ which is a hard constraint matrix, and $W_{{{L}}_{1}}$ which is a matrix that arises from the approximation of the $L_{1}$-norm stabilizer via an $L_{2}$-norm term. Details of all aforementioned matrices are provided in Vatankhah et al. \shortciteVRA:2017,VRA:2018b, but we note that it is $W_{{{L}}_{1}}$ which enables the algorithm to produce non-smooth and focused images of the subsurface that are more consistent with real geological structures. The scalar regularization parameter $\alpha$ balances the two terms in the objective function and is discussed further in Section 2.0.1

Noting that the inverse of a diagonal matrix is obtained at effectively zero computational cost, the objective function (2) is easily transformed to the standard Tikhonov form

[TABLE]

see e.g. Vatankhah et al. \shortciteVAR:2015, where $\tilde{\tilde{G}}={W_{\mathbf{d}}}GW^{-1}$, $\tilde{\mathbf{r}}={W_{\mathbf{d}}}(\mathbf{d}_{\mathrm{obs}}-G\mathbf{m}_{\mathrm{apr}})$, and $\mathbf{h}=W(\mathbf{m}-\mathbf{m}_{\mathrm{apr}})$. The solution of (3), dependent on the choice of $\alpha$, is then

[TABLE]

and the model update is given by

[TABLE]

An iteratively reweighted approach is used to find $\mathbf{m}(\alpha)$, dependent not only on $\alpha$, changing with each iteration, but also on matrix $W$ which changes at each step through the update of matrix $W_{{{L}}_{1}}$. The complete approach is described in Vatankhah et al. \shortcite[Algorithm $1$]VRA:2017 and Vatankhah et al. \shortcite[Algorithm $2$]VRA:2018b. At each stage of the algorithm upper and lower bounds on the physical parameters are imposed in order that the recovered model is reliable within known acceptable ranges, and regularization parameter $\alpha$ is adjusted automatically. The iteration is terminated when either the data predicted by the reconstructed model satisfies the observed data to within a $\chi^{2}$ value relating to the noise level, or a maximum number of iterations $K_{\mathrm{max}}$ is reached without the predicated data satisfying the $\chi^{2}$ estimate.

When both $m$ and $n$ are small, $\mathbf{h}(\alpha)$ can be found using the FSVD $\tilde{\tilde{G}}=U\Sigma V^{T}$, where matrices $U\in\mathcal{R}^{m\times m}$ and $V\in\mathcal{R}^{n\times n}$ are orthogonal with columns $\mathbf{u}_{i}$ and $\mathbf{v}_{i}$, and $\Sigma\in\mathcal{R}^{m\times n}$ is the matrix of singular values $\sigma_{i}$, ordered from large to small. The solution obtained using the FSVD is described in Vatankhah et al \shortcite[Algorithm $1$]VRA:2017, and see also Paoletti et al. \shortcitePHHF:2014 and Vatankhah et al. \shortciteVAR:2015. The availability of the FSVD makes it possible to estimate regularization parameter $\alpha$ cheaply using standard parameter-choice techniques [\citenameXiang & Zou 2013, \citenameChung & Palmer 2015]. Unfortunately the calculation of the SVD for large, or even moderate, under-determined systems is not practical; the cost is approximately $6nm^{2}+20m^{3}$, Golub $\&$ Van Loan \shortciteGoLo:2013. The traditional alternative is the use of a hybrid method such as the iterative LSQR algorithm that can be used to project (3) onto a Krylov subspace of smaller dimension and for which an SVD of the projected problem is then efficiently used to yield the subspace solution. This solution is projected back to the original full space at minimal additional computational cost. The SVD for the projected problem also facilitates the use of parameter-choice algorithms to find an optimal $\alpha$, denoted by $\alpha_{\mathrm{opt}}$, also at minimal additional computational cost.

2.0.1 Regularization parameter-choice method

Here we use the method of unbiased predictive risk estimation (UPRE) to find $\alpha_{\mathrm{opt}}$. This a-posteriori rule for choosing the Tikhonov regularization parameter method is well-described in Vogel \shortciteVogel:2002, and has been extensively applied for the inversion of data when an estimate of the noise covariance is available, including for the inversion of geophysical data, [\citenameVatankhah et al. 2015, \citenameVatankhah et al. 2017] and for more general inversion problems [\citenameRenaut et al. 2017, \citenameChung & Palmer 2015]. Using the SVD of matrix $\tilde{\tilde{G}}$, the UPRE function to be minimized is given by

[TABLE]

Here, ${m^{*}}$ indicates the number of non-zero singular values. This is the numerical rank of the matrix, namely $m^{*}=m$ when UPRE is applied using the FSVD for the underdetermined problem. The optimum regularization parameter, $\alpha_{\mathrm{opt}}$, is found by evaluating (6) on a range of $\alpha$, between minimum and maximum $\sigma_{i}$; equivalently $\alpha_{\mathrm{opt}}=\textnormal{arg}\min_{\sigma_{\mathrm{min}}\leq\alpha\leq\sigma_{\mathrm{max}}}\{U(\alpha)\}$.

2.1 Randomized Singular Value Decomposition

While the hybrid-LSQR methodology is effective and practical; acceptable results are obtained as compared with the Hybrid-FSVD solution, see for example Renaut et al. \shortciteRVA:2017 and Vatankhah et al. \shortciteVRA:2017, the approach is still cost-limited due to the need to build a relatively large Krylov space for the solution when both $m$ and $n$ are large. Recent approaches based on randomization have presented an interesting alternative for tackling problems requiring high resolution of the subsurface, [\citenameXiang & Zou 2013, \citenameVoronin et al. 2015, \citenameVatankhah et al. 2018a, \citenameVatankhah et al. 2018b]. In this case, random sampling is used to construct a low-dimensional subspace that approximates the column space of the model matrix and maintains the most dominant spectrum of the original matrix, [\citenameHalko et al. 2011]. Standard deterministic matrix decomposition methods such as the SVD, or eigen-decomposition, can then be used to compute a low-rank approximation of the original matrix. Specifically, in the context of potential field inversion, it is desirable to find a $q$-rank matrix $\tilde{\tilde{G}}_{q}$, which is as close as possible to $\tilde{\tilde{G}}$ in the least-squares sense, while at the same time the target rank $q$ is as small as possible in order that the inversion process is fast. We note that of course the best rank $q$ approximation in the least squares sense is given by the exact truncated SVD of $\tilde{\tilde{G}}$ with $q$ terms [\citenameGolub & van Loan 2013]. It is generally, however, not practical to calculate the truncated SVD for large scale problems. Thus here we focus on the RSVD and carefully describe the approach that was presented in Vatankhah et al \shortciteVRA:2018b for the solution of the gravity inversion problem.

The fundamental aspects of an RSVD algorithm consist of two stages. Here we present this for the under-determined matrix $\tilde{\tilde{G}}$. (i) A low-dimensional subspace is constructed that approximates the column space of $\tilde{\tilde{G}}^{T}$. The aim is to find a matrix $Q\in\mathcal{R}^{n\times q}$ with orthonormal columns such that $\tilde{\tilde{G}}\approx\tilde{\tilde{G}}QQ^{T}$; (ii) Given the near-optimal basis spanned by the columns of $Q$, a smaller matrix $B=\tilde{\tilde{G}}Q\in\mathcal{R}^{m\times q}$ is formed. This means that $\tilde{\tilde{G}}$ is restricted to the smaller subspace spanned by the basis from the columns of $Q$. Moreover, $B$ is a projection from high-dimensional space into a low-dimensional space which preserves the geometric structure of the matrix in a Euclidean sense [\citenameErichson et al. 2016]. Subsequently, $B$ can then be used to compute an approximate matrix decomposition for $\tilde{\tilde{G}}$ using a traditional algorithm. Step (i) is completely random and depends on the selection of a specific approach to find $Q$, while (ii) is deterministic. The fundamental approach, as presented for the gravity problem by Vatankhah et al. \shortcite[Algorithm $1$]VRA:2018b, as summarized in Algorithm 1, but with the inclusion of power iterations, is now discussed.

In step 1 a random test matrix $\Omega\in\mathcal{R}^{\ell\times m}$ is generated from a standard normal distribution. Probability theory guarantees that the columns of $\Omega$ are linearly independent. Here, $q+p=\ell\ll m$ and $p$ is a small oversampling parameter that provides the flexibility and effectiveness of the algorithm [\citenameHalko et al. 2011, \citenameErichson et al. 2016]. At step 2 a set of independent randomly-weighted linear combinations of the rows of $\tilde{\tilde{G}}$, or equivalently columns of $\tilde{\tilde{G}}^{T}$, are formed. The sketch matrix $Y$ has a much smaller number of rows than $\tilde{\tilde{G}}$. Ignoring for the moment power iterations at step 3, step 4 constructs $Q\in\mathcal{R}^{n\times\ell}$. The columns of $Q$ form an orthonormal basis for the range of $Y^{T}$; equivalently for the range of $\tilde{\tilde{G}}^{T}$. Given near-optimal $Q$, a smaller matrix $B$ is reconstructed in step 5. Therefore, large matrix $\tilde{\tilde{G}}$ is projected onto a low-dimensional space that captures most of the action of $\tilde{\tilde{G}}$. Choosing an optimal $q$ is highly dependent on the task. Generally, $q$ should be as small as possible so that the algorithm is fast and efficient, but simultaneously $q$ should be large enough that the dominant spectral properties of $\tilde{\tilde{G}}$ are accurately captured. We discuss the effect of the choice of $q$ for the inversion of gravity and magnetic data in Section 3.

Having obtained matrix $B$ at step 5, traditional SVD algorithms could then be used to compute the approximations to the first $q$ left singular vectors as well as the corresponding singular values for matrix $\tilde{\tilde{G}}$. The approximate right singular vectors could also be recovered, see e.g. Xiang and Zou \shortciteXiZo:2013. Alternatively, as discussed, with proof, in Vatankhah et al. \shortciteVRA:2018b, the much smaller matrix $B^{T}B\in\mathcal{R}^{\ell\times\ell}$ can be used to find the required SVD components for $B$ using the eigen-decomposition of $B^{T}B$. It was also demonstrated that the computational cost of this algorithm, excluding power iterations, is $O(\ell mn)$. Thus, it is feasible to compute the large singular values of a given matrix efficiently.

2.1.1 Randomized Singular Value Decomposition with Power Iterations

The quality of the RSVD which is obtained using Algorithm 1 without step 3 depends on the quality of the basis matrix $Q$ as providing a basis for the column space of $\tilde{\tilde{G}}^{T}$. Halko et al \shortcite[p224]Halko:2011 suggested an improvement of their Proto Algorithm for forming the basis matrix $Q$ and the associated RSVD that should lead to an improved approximation of the dominant spectrum. In the Prototype Algorithm for power iterations [\citenameHalko et al. 2011, p227], the sketch matrix $Y$ is obtained after first preprocessing matrix $\tilde{\tilde{G}}$ to give

[TABLE]

where integer $s$ specifies the number of power iterations. While this is specifically implemented as [\citenameHalko et al. 2011, Algorithm 4.3], Halko et al \shortciteHalko:2011 note that the approach is sensitive to floating point rounding errors, which reduces the quality of the $Q$ basis. Thus, Halko et al \shortciteHalko:2011 suggested their Algorithm 4.4 which orthonormalizes the columns of $Y$ between each application of $\tilde{\tilde{G}}$ and $\tilde{\tilde{G}}^{T}$. It is this latter approach, described in Algorithm 2, that we implement for step 3 within Algorithm 1, but here with our extension of the approach for the under-determined case.

To analyze the effect of power iterations on improving the accuracy of the computed matrix, a simplified upper bound on the expected error between the original and the $q$-rank computed matrices is given by

[TABLE]

[\citenameMartinsson 2016, \citenameErichson et al. 2016]. Here $e$ is Euler’s number, $\sigma_{q+1}$ is the $q+1$ largest singular value of the matrix $\tilde{\tilde{G}}$, $E$ denotes the expectation operator, and it is assumed that $p\geq 2$. Upper bound (8) indicates how parameters $p$, $q$, and $s$ can be used to control the approximation error. Note immediately that for $q=min(m,n)$ then $\sigma_{q+1}=0$ and $E\|\tilde{\tilde{G}}-\tilde{\tilde{G}}_{q}\|=0$. With increasing oversampling $p$, the second and third terms in the bracket tend toward zero which means that the bound approaches the theoretically optimum value $\sigma_{q+1}$ [\citenameErichson et al. 2016]. For larger values of the subspace iteration parameter $s$, $1/(2s+1)$ goes to zero and the error bound is reduced.

In Section 3 we will illustrate the application of Algorithm 1 for the inversion of synthetic gravity and magnetic data sets without power iterations at step 3. These results, especially for the magnetic data, suggest that power iterations as indicated in Algorithm 2 are needed to improve the approximation of the dominant spectral space. Our tests have shown that setting $s=1$ yields a good performance that trades-off between accuracy and computational time. We will see that the magnitude of the singular values of the magnetic kernel are larger for a given $q$ than their counterparts for the gravity kernel, and thus the upper bound estimate in (8) for the approximation error is consistently higher for the magnetic problem. Equivalently this leads to the need to implement a power iteration scheme to reduce the factor multiplying $\sigma_{q+1}$ in (8).

3 Synthetic examples

In Section 3.1 we first evaluate the performance of the RSVD algorithm without power iterations, comparing its performance for the solution of relatively small-scale gravity and magnetic problems under the same configuration, but the appropriate choice of model matrix $\tilde{\tilde{G}}$. The results are compared to those obtained using the Hybrid-FSVD in each case. To understand the performance of the algorithm we examine the spectrum of the approximate operator, as compared to that of $\tilde{\tilde{G}}$, for each problem in Section 3.1.1, and then examine the improvement obtained using the power iterations for $s=1$ in Section 3.1.2. A more complicated configuration for a structure with multiple bodies is then examined in Section 3.2 and confirms the conclusions obtained for the dipping dike example in Section 3.1.

In the simulations for the generation of the total field anomaly, the intensity of the geomagnetic field, the inclination, and the declination are selected as $47000$ nT, $50^{\circ}$ and $2^{\circ}$, respectively. The density contrast and the susceptibilities of the model structures embedded in a homogeneous non-susceptible background are $\rho=1$ g cm*-3* and $\kappa=0.1$ (SI unit), respectively. The bound constraints $0=\rho_{\mathrm{min}}\leq\rho\leq\rho_{\mathrm{max}}=1$ in units  g cm*-3* and $0=\kappa_{\mathrm{min}}\leq\kappa\leq\kappa_{\mathrm{max}}=0.1$ in SI units are imposed at each iteration of the gravity and magnetic inversions, respectively. Further, for all simulations we add Gaussian noise with zero mean and standard deviation $(\tau_{1}~{}|\mathbf{d}_{\mathrm{exact}}|_{i}+\tau_{2}~{}\mathrm{max}|\mathbf{d}_{\mathrm{exact}}|)$ to datum $i$, for chosen pairs $(\tau_{1},\tau_{2})$, where $\mathbf{d}_{\mathrm{exact}}$ is the exact data set, yielding $\mathbf{d}_{\mathrm{obs}}$ with a known distribution of noise for the error. This standard deviation is used to generate the matrix $W_{\mathbf{d}}$ in the data fit term. The values of $(\tau_{1},\tau_{2})$ are selected such that the signal to noise ratios, as given by

[TABLE]

are close for both gravity and magnetic noise-contaminated data. The values for $(\tau_{1},\tau_{2})$ and the resulting SNRs are specified in the captions of figures associated with the results for each data set. Then to test convergence of the update $\mathbf{m}^{(k)}$ at iteration $k$ we calculate the $\chi^{2}$ estimate,

[TABLE]

where $\mathbf{d}_{\mathrm{pre}}^{(k)}=G\mathbf{m}^{(k)}$, and which assesses the predictive capability of the current solution. When $(\chi^{2})^{(k)}\leq m+\sqrt{2m}$ the iteration terminates. Otherwise, the iteration is allowed to proceed in all cases to a maximum number of iterations $K_{\mathrm{max}}=50$. The regularization parameter $\alpha$ is adjusted with iteration $k$ and for $k>1$ is obtained in all cases using the UPRE function (6) with the appropriate approximate SVD terms. Based on the suggestion of Farquharson & Oldenburg \shortciteFar:2004, a large $\alpha$ is always used at the first iteration,

[TABLE]

as used in Vatankhah et al. \shortcite[eq.19]VRA:2017. In all simulations we use a fixed oversampling parameter, $p=10$, assume $\mathbf{m}_{\mathrm{apr}}=\mathbf{0}$, impose physically reasonable constraints on the model parameters, and for $W_{\mathrm{depth}}$ take $\beta=0.8$ and $1.4$, for gravity and magnetic inversions, respectively, with values close to those suggested by Li $\&$ Oldenburg \shortciteLiOl:98. In the results we examine the dependence of the solution on the choice of $q$ for the rank $q$ approximation and record (i) the number of iterations $K$ required, (ii) the relative error progression with increasing $k$ as given by

[TABLE]

(iii) the relative error in the rank $q$ approximation to $\tilde{\tilde{G}}$,

[TABLE]

(iv) $\alpha^{(k)}$ with increasing $k$, and (v) all values at the final iteration $K$ as well as the time for the iterations. Moreover, in recording the computational time to convergence, in all cases we do not include the calculation of the original model matrix $G$. Computations are performed on a desktop computer with Intel(R) Xeon(R) W-2133 CPU 3.6 GHz processor and 32 GB RAM.

3.1 Small-scale model consisting of two dipping dikes

The small-scale but complicated structure of two dipping dikes, illustrated in Fig. 1, makes it computationally feasible to compare the solutions for the gravity and magnetic inverse problems using Algorithm 1 without power iteration, with the solutions obtained using the Hybrid-FSVD. The data for the problem, the vertical component of the gravity and the total magnetic field, are generated on the surface on $30\times 30=900$ grid with grid spacing $50$ m. The noisy gravity and magnetic data in each case are illustrated in Figs. 2 and 2.

The subsurface volume is discretized into $9000$ cubes of size $50$ m in each dimension, corresponding to $10$ slices in depth for a cross section of $30$ by $30$ cells. The resulting matrix $\tilde{\tilde{G}}$ is of size $900\times 9000$ and thus facilitates solutions using the FSVD for comparison with the RSVD in the inversion algorithm. The results obtained using the Hybrid-FSVD and Hybrid-RSVD, with the choices $q=100$, $150$, $200$, $300$, $500$, $700$, and $900$, are detailed in Table 1 and Table 2, respectively. The results for the inversions using the Hybrid-FSVD for the gravity and magnetic data inversions are given in Fig. 3 and Fig. 5, respectively. For comparison the results using $q=200$ are illustrated in Fig. 4 and Fig. 6 for the gravity inversion and magnetic data inversion, respectively.

The results presented for the inversion of the gravity data using the Hybrid-FSVD show that the iteration terminates after just $12$ iterations, and as indicated in Fig. 3 the reconstructed model is in good agreement with the original model. A sharp and focused image of the subsurface is obtained, and while the depths to the top of the structures are consistent with those of the original model, the extensions of the dikes are overestimated for the left dike and underestimated for the right dike. Figs. 3 and 3 illustrate the progression of the relative error and regularization parameter at each iteration, respectively, and Fig. 3 shows the UPRE function at the final iteration. These figures are presented for comparison with the results obtained using the Hybrid-RSVD algorithm.

For the same gravity problem the results using the Hybrid-RSVD algorithm for very small values of $q$, $q<100$, are not acceptable. With increasing $q$ the solution improves until at $q=m$ the solution matches the Hybrid-FSVD solution. For the reported choices of $q$ all the solutions terminate prior to $K_{\mathrm{max}}$, with $K=12$ for $q\geq 150$, and demonstrating that the $\chi^{2}$ estimate, (10) is satisfied. We can see that for a suitable value of $q$, the RSVD leads to a solution that is close to that achieved using the Hybrid-FSVD. These conclusions confirm the results in Vatankhah et al. \shortciteVRA:2018b; the Hybrid-RSVD algorithm can be used with $q\gtrsim(m/6)$ for the inversion of gravity data. We illustrate the results of the inversion using $q=200$ in Fig. 4.

The results presented for the inversion of the magnetic data using the Hybrid-FSVD show that the iteration terminates after $8$ iterations, and as indicated in Fig. 5 the reconstructed model is in reasonable agreement with the original model. Again a focused image of the subsurface is obtained, but the overestimation of the depth of the left dike is greater as compared with the results in Fig. 3. Correspondingly, the relative error is larger than that achieved in the inversion of the gravity data.

For the same magnetic problem the results using the Hybrid-RSVD algorithm are not acceptable with reasonable choices for $q$, indeed the algorithm does not terminate prior to $K=K_{\mathrm{max}}$ for $q<500$, and the cost therefore increases significantly. As compared to the inversion of gravity data, larger values of $q$, are required to yield acceptable results. Observe, for example, that $q=200$ is not a suitable choice because the relative error of the reconstructed model is unacceptable, and the predicted data does not satisfy the observed data at the given noise level; (10) is not satisfied for $k\leq 50$. Using $q=500$, the reconstructed model relative error is reduced and the inversion terminates at $9$ iterations with an acceptable $\chi^{2}$ value. We deduce that for the inversion of the magnetic data it is necessary to take $q$ larger than we would use for the inversion of the gravity data, as is indicated by the relatively larger estimates for the rank-$q$ error, $RG^{(K)}$ in Tables 1 and 2, respectively. To demonstrate the impact of the choice of $q$ we show the results of the inversion for $q=200$ and $q=500$, in Figs. 6 and 7, respectively.

3.1.1 The Spectral Properties

We now compare the singular values $\sigma_{i}^{(k)}$ and $(\sigma_{i}^{(k)})_{q}$ for the matrices $\tilde{\tilde{G}}^{(k)}$ and $\tilde{\tilde{G}}^{(k)}_{q}$, respectively, for the gravity and the magnetic problems in Figs. 8 and 9, respectively. In both figures, we show the values for $q=200$ and $q=500$ at iteration $k=8$. It is immediate from these plots that the Hybrid-RSVD algorithm does not capture the dominant spectrum of the original matrix for the magnetic problem as closely as is the case for the gravity problem. This clarifies why the magnetic inversion requires larger values of $q$ in order for the inversion to converge. It is also evident that the singular values of both gravity and magnetic problems decay rather slowly, after an initial fast decay. Generally, the RSVD algorithm is more efficient when applied to matrices with rapidly decaying singular values, as can be seen from the error estimate (8).

3.1.2 Power Iterations

Now as discussed in Section 2.1 the error estimate (8) will decrease with increasing $s$. Having seen that the lack of power iterations leads to lack of convergence for the inversion of the magnetic data unless $q$ is taken relatively large, $q\gtrsim m/2$, as compared to just $q\gtrsim m/6$ for the gravity problem, we investigate the power iteration step given by Algorithm 2 to improve the column space approximation of $\tilde{\tilde{G}}^{T}$. We therefore repeat the simulations for the data of the two dike problem illustrated in Figs. 2-2 but employing a power iteration with $s=1$. The results are presented in Tables 3 and 4 for gravity and magnetic data, respectively, and indicate improvements as compared to the results presented without power iterations. For both problems the number of iterations $K$ is generally reduced, once $q$ is large enough, and the error is generally decreased for a result that used the same number of iterations $K$ with and without power iterations. Moreover, the results are achieved without a large increase in computational cost, where the algorithm converged both with and without power iterations. But the major impact is that it is possible to take a much smaller $q$ to obtain convergence of the inversion of the magnetic problem with reasonable computational cost. In both cases it is sufficient to now take $q=200$ to obtain converged solutions for a relatively small $K$, $11$ and $12$, respectively.

The results for magnetic data inversion with $q=200$ and $s=1$ are illustrated in Fig. 10 for comparison with Fig. 6 obtained without the power iterations. A noticeable improvement in the reconstructed model is obtained. To further show the impact of the power iterations we also show the singular values for the power iterations with $s=1$ and $q=200$ for both gravity and magnetic problems in Figs. 11-11, comparing with Fig. 8 and Fig. 9, respectively. These plots demonstrate that the power iterations have indeed improved the accuracy of the estimated singular values.

Naturally these results raise the question as to whether it is better to apply the Hybrid-RSVD algorithm without power iterations and a large choice for $q$ or to use power iterations and take a smaller $q$. But the main purpose of using Algorithm 1 is to make it feasible, in terms of both memory and computational cost, to find accurate solutions of large scale problems. Indeed the aim is to solve problems which are either too expensive to solve using the Hybrid-FSVD or cannot be solved at all using the Hybrid-FSVD. We discuss this further for a larger problem in Section 3.2

3.2 Model of multiple bodies

We now study the application of the Hybrid-RSVD algorithm for the solution of a larger and more complex model consisting of six different bodies with different shapes, dimensions and depths, as shown in the perspective view in Fig. 12 and the six plane-sections in Fig. 13. The data for the problem, the vertical component of the gravity and the total magnetic field, are generated on the surface on a $100\times 80$ grid with $100$ m spacing. The noisy gravity and magnetic data in each case are illustrated in Figs. 14 and 14.

To perform the inversion, the subsurface volume is discretized into $100\times 80\times 10$ cubes of size $100$ m in each dimension. The resulting matrix $\tilde{\tilde{G}}$ is of size $8000\times 80000$ and is too large for efficient use of the Hybrid-FSVD. The results obtained using the Hybrid-RSVD algorithm with the choices $q=1000$, $1500$, $2500$, and $4200$ are detailed in Tables 5 and 6, respectively, both with and without power iterations.

As with the inversion of the two-dike problem, it is immediate that the inversions for the gravity problem are acceptable, without power iterations, for far smaller $q$ than for the magnetic case. Applying the power iterations reduces the size of $q$ that is required to obtain convergence and excellent results are obtained with just $q=1000$ and a time that is little more than is required for $q=1500$ and no power iterations. This suggests that we may use $q\gtrsim m/8$ with $s=1$ for the power iterations. As for the two-dike problem the magnetic inversion iteration does not converge to the required $\chi^{2}$ level within $50$ iterations, except when we take $q=4200>m/2$. Further, the relative errors are large and the computational cost is high. Including the power iterations yields convergence at $K=14$ when $q=2500$ and a smaller relative error for an acceptable computational cost. These results suggest that it may be acceptable to take $q\gtrsim m/4$ in the inversion of magnetic data with the Hybrid-RSVD algorithm combined with $s=1$ power iterations. This choice yields run times for the magnetic inversion that are comparable to those for the gravity inversion.

The cross-sections for the inversions using the power iterations and $q=2500$ for the gravity and magnetic data are given in Fig. 15 and Fig. 16, respectively. The perspective view of these solutions is also given in Figs. 17 and 18, respectively. We observe that the horizontal borders of the reconstructed models, in both inversions, are in good agreement with those of the original model, but that additional structures appear at depth. Here, the reconstructed model of the magnetic susceptibility exhibits more artifacts at depth and has a higher relative error than achieved for the gravity reconstruction. On the other hand, the magnetic structure better illustrates the dip of both dikes which is significant for accurate geophysical interpretation of the structures. Moreover, these results indicate that joint interpretation of the individual magnetic and gravity inversions may improve the quality of the final subsurface model.

4 Real data

The Hybrid-RSVD algorithm with and without power iterations, is now applied for the inversion of total magnetic field data that have been obtained over a portion of the Wuskwatim Lake region in Manitoba, Canada, Fig. 19111Results showing the application of the methodology for real gravity data were already presented in Vatankhah et al. \shortciteVRA:2018b.. The given area lies within a poorly exposed meta-sedimentary gneiss belt consisting of paragneiss, amphibolite, and migmatite derived from Proterozoic volcanic and sedimentary rocks [\citenamePilkington 2009]. The total field anomaly shows magnetic targets elongated in the NE-SW direction. A data-space inversion algorithm with a Cauchy norm sparsity constraint on model parameters was applied by Pilkington \shortcitePi:09 for the inversion of this data set. Furthermore, the results of the inversion algorithm of Li $\&$ Oldenburg \shortciteLiOl:98 are also presented in Pilkington \shortcitePi:09. Therefore, the results presented here can be compared with the results of presented in Pilkington \shortcitePi:09 and for consistency we thus use a grid of $64\times 64$ data points with $100$ m spacing and a uniform subsurface discretization of $64\times 64\times 20=81920$ blocks. The intensity of the geomagnetic fields, the inclination and the declination are $60000$ nT, $78.5^{\circ}$, $5.3^{\circ}$, respectively. As for the simulations we set $K_{\mathrm{max}}=50$, $\mathbf{m}_{\mathrm{apr}}=\mathbf{0}$, and impose bound constraints on the model parameters. In this case, these are $0=\kappa_{\mathrm{min}}\leq\kappa\leq\kappa_{\mathrm{max}}=0.2$ SI unit [\citenamePilkington 2009]. The values of parameter $q$ are selected based on the presented rules in Section 3. We select $q=1100>m/4$ and $q=2100>m/2$ with and without power iterations, respectively.

The results of the inversions, as given in Table 7, demonstrate that the methodology generates converged solutions using a limited number of iterations and at computational cost on the order of a few minutes only. Overall these results demonstrate the feasibility of using the Hybrid-RSVD algorithm for the inversion of large-scale geophysical data sets. Three plane-sections of the reconstructed model obtained using $s=1$ for the power iterations are illustrated in Fig. 20. The anomaly produced by this model is shown in Fig. 19. The progression of the regularization parameter at each iteration and the UPRE functional at the final iteration are also presented in Fig. 21. Furthermore, Fig. 22 illustrates a $3$-D view of the model for cells with $\kappa>0.05$. Our results indicate that, generally, there are three main subsurface targets. The target in the South-East of the area starts from $300$ m and extends to $400$ m; it is not as deep as the other two targets. The target in the central part of the area, is elongated in the SW-NE direction, and starts from about $300$ m and extends to about $800$ m in depth. This target in its northeastern part is divided into two sub-parallel targets. The third main target, located in the central north part of the area, is the deepest target which starts at about $400$ m and extends to about $900$ m. The results in the shallow and intermediate layers are in agreement with the results of Pilkington \shortcitePi:09, but at depth the results presented here are more focused.

5 Conclusions

We present an algorithm for fast implementation of large-scale focusing inversion of both magnetic and gravity data. The algorithm is based on a combination of the $L_{1}$-norm regularization strategy with the RSVD. For large-scale problems, the powerful concept of the RSVD provides an attractive and indeed fast alternative to methods such as the LSQR algorithm. Here we have presented a comprehensive comparison of the Hybrid-RSVD methodology with power iterations for the inversion of gravity and magnetic data. Generally, we have shown that there is an important difference between gravity and magnetic inverse problems when approximating a $q$-rank matrix from the original matrix. For the inversion of magnetic data it is necessary to take larger values of $q$, as compared with the inversion of gravity data, in order that a suitable approximation of the system matrix is obtained. Furthermore, including power iterations within the algorithm improves the approximation quality. Indeed the RSVD obtained using the power iterations step yields approximation of the dominant singular space even for small choices of $q$. Thus, the RSVD with power iterations yields an efficient strategy when the singular values of input matrix decay slowly. In particular, the presented methodology can be used for other geophysical data sets and that the choice of the rank $q$ approximation will depend on the spectral properties of the relevant kernel matrices. If the RSVD without power iteration does not approximate the dominant singular values then power iterations should be included to improve the quality of estimated singular values. Our results also demonstrate that it is possible to obtain a reasonable lower estimate for $q$, for both gravity and magnetic data inversions, based on the number of data measurements, $m$. In conclusion, we have demonstrated that it is feasible to use an efficient RSVD methodology for problems that are too large to be handled using the full SVD. Application of the RSVD for the joint inversion of gravity and magnetic data is a topic for future work.

Acknowledgements.

The authors would like to thank Dr. Mark Pilkington for providing data from Wuskwatim Lake area. This study is supported by the National Key R&D Program (NO. 2016YFC0600109) and NSF of China (NO. 41604087 and 41874122).