Improving the use of the randomized singular value decomposition for the inversion of gravity and magnetic data
Saeed Vatankhah Shuang Liu Rosemary A. Renaut Xiangyun Hu Jamaledin Baniamerian

Nonnegative sum-symmetric matrices, optimal-score partitions, and optimal resource allocation 



    Iosif Pinelis       



    Received: date / Accepted: date
=================================================================================================




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






Inverse theory; Numerical approximation and analysis;  Gravity anomalies and Earth structure;  Magnetic anomalies: modeling and interpretation.





§ 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 <cit.>. The inversion of acquired data is one of the important steps in the interpretation process <cit.>. 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. <cit.>. 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), <cit.>. 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 <cit.>. 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 3D inversion of gravity data, Vatankhah et al. VRA: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 ≪ 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 ≳ (m/6).[The 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  and  algorithms.  For the   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  algorithm, as presented in  Vatankhah et al. VRA:2018b, for the  inversion of magnetic potential data.  Our results demonstrate that direct application of the  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 ≳ 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 <ref> 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 <ref>. Then, the concept of the RSVD is reviewed in Section <ref>, with the  extension applying power iterations  presented in Section <ref>. In Section <ref> we show the results of using the  algorithm on two different synthetic models. A small model in Section <ref> 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 <ref>. 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 <ref>. Conclusions and discussions on future work are given in Section <ref>.






§  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 ∈, and  measured potential field data are also  stacked in a vector ∈. The measurements are connected to the model parameters via the model matrix G ∈, m≪ n, which is the forward model operator which maps from model to data spaces, yielding the under-determined linear  system

    = G,

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  collects the values of the unknown susceptibilities of the cells and  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 RaBa: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 Haaz, see for example Boulanger & Chouteau BoCh: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  (<ref>); its solution is neither unique nor stable. Here we consider the L_1-norm stabilized  solution of (<ref>) as  presented by Vatankhah et al. VRA: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 0 values when no prior information is available. Overall, the formulation finds the minimum of the global objective function P^α() as given by

    = {P^α()} ={(G-)  _2^2 + α^2 W(-) _2^2}.
 
The diagonal entries of the data weighting matrix  are estimates for the  inverse of the standard deviations of the  independent noise in the data,  encodes the prior knowledge of the solution, and stabilizing matrix  W is the product of three diagonal matrices,  which is a depth-weighting matrix with diagonal entries z^-β at depth z,  which is a hard constraint matrix, and  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. VRA:2017,VRA:2018b, but we note that it is  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 α balances the two terms in the objective function and is discussed further in Section <ref>

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

    P^α()= - _2^2 + α^2 _2^2,

see e.g. Vatankhah et al. VAR:2015, where =GW^-1, =(-G), and =W(-). The solution of (<ref>), dependent on the choice of α, is then

    (α)=(^T+α^2 I_n)^-1^T,

and the model update is given by

    (α)=+W^-1 (α).

An iteratively reweighted approach is used to find  (α), dependent not only on α, changing with each iteration, but also on matrix W which changes at each step through the update of matrix .  The complete approach is described in  Vatankhah et al. [Algorithm 1]VRA:2017 and Vatankhah et al. [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 α is adjusted automatically. The iteration is terminated when either the  data predicted by the reconstructed model satisfies the observed data to within a χ^2 value relating to the noise level, or a maximum number of iterations   is reached without the predicated data satisfying  the χ^2 estimate. 

When both m and n are small,  (α) can be found  using the FSVD = U Σ V^T, where matrices U∈ and V∈ are orthogonal with columns _i and _i, and Σ∈ is the matrix of singular values σ_i, ordered from large to small. The solution obtained using the  FSVD is  described in Vatankhah et al [Algorithm 1]VRA:2017, and see also Paoletti et al. PHHF:2014 and  Vatankhah et al. VAR:2015. The availability of the FSVD makes it  possible to estimate regularization parameter α cheaply using standard parameter-choice techniques <cit.>. Unfortunately the calculation of the  SVD for large, or even moderate, under-determined systems is not practical;  the cost is approximately 6nm^2+20 m^3, Golub & Van Loan GoLo:2013. The traditional alternative is the use of a hybrid method such as the iterative LSQR algorithm that can be used to project (<ref>) 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 α, denoted by ,  also at minimal additional computational cost. 



  §.§.§ Regularization parameter-choice method

Here  we use the method of unbiased predictive risk estimation (UPRE) to find . This a-posteriori rule for choosing the Tikhonov regularization parameter method is well-described in Vogel Vogel: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, <cit.> and for more general inversion problems <cit.>. Using the SVD of matrix , the UPRE function to be minimized is given by

    U(α)=∑_i=1^m^*( 1/σ_i^2 α^-2 + 1)^2 (_i^T)^2 + 2 ( ∑_i=1^m^*σ_i^2/σ_i^2+α^2) - m^*.
  
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, , is found by evaluating (<ref>) on a range of α, between minimum and maximum σ_i; equivalently =≤α≤{U(α)}. 





 §.§ Randomized Singular Value Decomposition

While the hybrid-LSQR methodology is effective and practical; acceptable results are obtained as compared with the  solution, see for example Renaut et al. RVA:2017 and Vatankhah et al. VRA: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,  <cit.>. 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, <cit.>. 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 _q, which is as close as possible to   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  with q terms <cit.>. 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 VRA: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 .  (i) A low-dimensional subspace is constructed that approximates the column space of ^T. The aim is  to find a matrix Q ∈ℛ^n × q with orthonormal columns such that ≈ Q Q^T; (ii) Given the near-optimal basis spanned by the columns of Q, a smaller matrix B= Q ∈ℛ^m × q is formed. This means that  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 <cit.>. Subsequently, B can then be used to compute an approximate  matrix decomposition for  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.  [Algorithm 1]VRA:2018b, as summarized in Algorithm <ref>, but with the inclusion of power iterations,  is now discussed.

In step <ref> a random test matrix Ω∈ℛ^ℓ× m is generated from a standard normal distribution. Probability theory guarantees that the columns of Ω are linearly independent. Here, q+p=ℓ≪ m and p is a small oversampling parameter  that provides the flexibility and effectiveness of the algorithm <cit.>.  At step <ref> a set of independent randomly-weighted linear combinations of the rows of , or equivalently columns of ^T, are formed. The sketch matrix Y has a much smaller number of rows than . Ignoring for the moment power iterations at step <ref>,  step <ref> constructs  Q ∈ℛ^n ×ℓ. The columns of Q form an orthonormal basis for the range of Y^T; equivalently for the  range of ^T. Given near-optimal Q, a smaller matrix B is reconstructed in step <ref>. Therefore, large matrix  is projected onto a low-dimensional space that captures most of the action of . 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  are accurately captured. We discuss the effect of the choice of q for the inversion of  gravity and magnetic data in Section <ref>.

Having obtained matrix B at step <ref>, 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 . The approximate right singular vectors could also be recovered,  see e.g. Xiang and Zou XiZo:2013. Alternatively, as discussed, with proof, in Vatankhah et al. VRA:2018b, the much smaller matrix B^TB∈ℛ^ℓ×ℓ can be used to find the required SVD components for B using the eigen-decomposition of B^TB. It was also demonstrated that the computational cost of this algorithm, excluding power iterations,  is  O(ℓ mn). Thus, it is feasible to compute the large singular values of  a given matrix efficiently. 



 




  §.§.§ Randomized Singular Value Decomposition with Power Iterations

The quality of the RSVD which is obtained using Algorithm <ref> without step <ref> depends on the quality of the basis matrix Q as providing a basis for the column space of ^T. Halko et al [p224]Halko:2011  suggested an improvement of their  Algorithm for forming the basis matrix Q and the associated RSVD that should lead to an improved approximation of the dominant spectrum. In the  Algorithm for power iterations <cit.>, the sketch matrix Y is obtained after  first preprocessing   matrix  to give 

    ^(s)=(^T)^(s),

where integer s specifies the number of power iterations. While this is specifically implemented as <cit.>, Halko et al Halko:2011 note that the approach is sensitive to floating point rounding errors, which reduces the quality of the Q basis. Thus, Halko et al Halko:2011 suggested their  which orthonormalizes the columns of Y between each application of  and ^T. It is this latter approach, described in Algorithm <ref>,  that we implement for step <ref> within Algorithm <ref>,  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

    E - _q ≤[1+√(q/p-1)+e√(q+p)/p.√(min(m,n)-q)]^1/2s+1σ_q+1,

 <cit.>. Here e is Euler's number, σ_q+1 is the q+1 largest singular value of the matrix 
, E denotes the expectation operator, and it is assumed that p≥ 2. Upper bound (<ref>) indicates  how parameters p, q, and s can be used to control the approximation error.  Note immediately that for q=min(m,n)  then σ_q+1=0 and  E - _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 σ_q+1 <cit.>. For larger values of the subspace iteration parameter  s, 1/(2s+1) goes to zero and      the error bound is reduced. 

 In Section <ref> we will illustrate the application of   Algorithm <ref> for the inversion of synthetic  gravity and magnetic data sets without power iterations at step <ref>. These results, especially for the magnetic data, suggest that power iterations as indicated in Algorithm <ref> 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 (<ref>) 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 σ_q+1 in (<ref>).
 
 






§ SYNTHETIC EXAMPLES

In Section <ref>  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 . The results are compared to those obtained using the  in each case. To understand the performance of the algorithm we examine the spectrum of the approximate operator, as compared to that of , for each problem in Section <ref>, and then examine the improvement obtained using  the power iterations for s=1 in Section <ref>. A more complicated configuration for a structure with multiple bodies is then examined in  Section <ref> and confirms the conclusions obtained for the dipping dike example in Section <ref>.

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^∘ and 2^∘, respectively. The density contrast  and the susceptibilities of the model structures embedded in a homogeneous non-susceptible background are  ρ=1 g cm^-3 and κ=0.1 (SI unit), respectively.  The bound constraints  0=ρ_min≤ρ≤ρ_max=1 in units  g cm^-3 and 0=κ_min≤κ≤κ_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 (τ_1 ||_i + τ_2 max|  |) to datum i, for chosen pairs (τ_1, τ_2), where  is the exact data set, yielding  with a known distribution of noise for the error. This standard deviation is used to generate the  matrix  in the data fit term. The values of 
(τ_1,τ_2) are selected such that the signal to noise ratios, as given by 

    SNR=20 log_10_2/-_2,

are close for both gravity and magnetic noise-contaminated data. The values for (τ_1, τ_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 ^(k) at iteration k we calculate the χ^2 estimate, 

    (χ^2)^(k)=(-^(k))_2^2,
 
where ^(k) = G^(k), and which assesses the predictive capability of the current solution. When (χ^2)^(k)≤ m+√(2m) the iteration terminates. Otherwise, the iteration is allowed to proceed in all cases to a maximum number of iterations =50.   The regularization parameter α is adjusted with iteration k and for k>1 is obtained in all cases using the UPRE function (<ref>) with the appropriate approximate SVD terms. Based on the suggestion of Farquharson & Oldenburg Far:2004, a large α is always used at the first iteration, 

    α^(1) = (n/m)^3.5σ_1/mean(σ_i),
 
as used in Vatankhah et al. [eq.19]VRA:2017. In all simulations we use a fixed oversampling parameter, p=10, assume =0, impose physically reasonable constraints on the model parameters, and for   take β=0.8 and 1.4,  for gravity and magnetic inversions, respectively,  with values close to those suggested by Li & Oldenburg LiOl: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 

    RE^(k)=_exact- ^(k)_2/_exact_2

(iii) the relative error in the rank q approximation to , 

    RG^(k)=-_q^(k)_2/_2

(iv) α^(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.

 
 

 §.§ Small-scale model consisting of two dipping dikes

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




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   is of size 900 × 9000 and thus facilitates solutions using the  FSVD  for comparison with the RSVD in the inversion algorithm.  The results obtained using the  and , with the choices q=100, 150, 200, 300, 500, 700, and 900, are detailed in Table <ref> and Table <ref>, respectively. The results for the  inversions using the  for the gravity and magnetic data inversions are given in Fig. <ref> and Fig. <ref>, respectively. For comparison the results using q=200 are illustrated in Fig. <ref> and Fig. <ref> for the gravity inversion and  magnetic data inversion, respectively. 


The results presented for the inversion of the gravity data using the  show that the iteration terminates after just 12 iterations, and as indicated in Fig. <ref>  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. <ref> and  <ref> illustrate the progression of the relative error and regularization parameter at each iteration, respectively, and Fig. <ref> shows the UPRE function at the final iteration. These figures are presented for comparison with the results obtained using the  algorithm. 

For the same gravity problem the results using the  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  solution. For the reported choices of q all the solutions terminate prior to , with K=12 for q≥ 150, and demonstrating that the χ^2 estimate, (<ref>) 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 . These conclusions confirm the results in Vatankhah et al. VRA:2018b; the  algorithm can be used with q ≳  (m/6) for the inversion of gravity data.  We illustrate the results of the inversion using q=200 in Fig. <ref>. 









The results presented for the inversion of the magnetic  data using the  show that the iteration terminates after 8 iterations, and as indicated in Fig. <ref>  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. <ref>. 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  algorithm  are not acceptable with reasonable choices for q, indeed the algorithm does not terminate prior to K= 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; (<ref>) is not satisfied for k≤ 50. Using q=500, the reconstructed model  relative error is reduced and the  inversion terminates at 9 iterations with an acceptable χ^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 <ref> and <ref>, 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. <ref> and <ref>, respectively. 
















  §.§.§ The Spectral Properties

We now compare the singular values σ_i^(k) and (σ_i^(k))_q for the matrices ^(k)  and ^(k)_q, respectively, for the gravity and the magnetic problems in Figs. <ref> and <ref>, 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  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 (<ref>). 


 




  §.§.§ Power Iterations

Now as discussed in Section <ref> the error estimate (<ref>) 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 ≳ m/2, as compared to just q≳ m/6 for the gravity problem, we investigate the power iteration step given by Algorithm <ref> to improve the column space approximation of  ^T. We therefore repeat the simulations for the data of the two dike problem illustrated in Figs. <ref>-<ref> but employing a power iteration with s=1. The results are presented in Tables <ref> and <ref> 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. <ref> for comparison with Fig. <ref> 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. <ref>-<ref>, comparing with Fig. <ref> and Fig. <ref>, 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  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 <ref>  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  or cannot be solved at all using the . We discuss this further for a larger problem in Section <ref>




 §.§ Model of multiple bodies

We now study the application of the  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. <ref> and the six plane-sections in Fig. <ref>. The data for the problem, the vertical component of the gravity and the total magnetic field, are generated on the surface on a   100 × 80 grid with 100 m spacing.  The noisy gravity and magnetic data in each case are illustrated in Figs. <ref> and <ref>. 






To perform the inversion, the subsurface volume is discretized into 100 × 80 × 10  cubes of size 100 m in each dimension. The resulting matrix   is of size 8000 × 80000 and is too large for efficient use of the .  The results obtained using the  algorithm with the choices q=1000, 1500, 2500, and 4200 are detailed in Tables <ref> and <ref>, 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≳ 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 χ^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≳ m/4 in the inversion of magnetic data with the  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. <ref> and Fig. <ref>, respectively. The perspective view of these solutions is also given in Figs. <ref> and <ref>, 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. 






















§ REAL DATA

The  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. <ref>[Results showing the application of the methodology for real gravity data were already presented in Vatankhah et al. VRA: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 <cit.>. 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 Pi:09 for the inversion of this data set. Furthermore, the results of the inversion algorithm of  Li & Oldenburg LiOl:98 are also presented in Pilkington Pi:09. Therefore, the results presented here can be compared with the results of presented in Pilkington Pi:09 and for  consistency we thus use a grid of 64 × 64 data points with 100 m spacing and a uniform subsurface discretization of 64 × 64 × 20 = 81920 blocks. The intensity of the geomagnetic fields, the inclination and the declination are 60000 nT, 78.5^∘, 5.3^∘, respectively. As for the simulations we  set =50, =0, and impose bound constraints on the model parameters. In this case, these are 0=κ_min≤κ≤κ_max=0.2 SI unit <cit.>. The values of parameter q are selected based on the presented rules in Section <ref>. 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 <ref>, 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  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. <ref>. The anomaly produced by this model is  shown in Fig. <ref>. The progression of the regularization parameter at each iteration and the UPRE functional at the final iteration are also presented in Fig. <ref>. Furthermore,  Fig. <ref> illustrates a 3-D view of the model for cells with κ>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 Pi:09, but at depth the  results presented here are more focused.
























§ 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  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.





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).








[Blakely 1996]Blakely Blakely, R. J. 1996, Potential Theory in Gravity & Magnetic Applications, Cambridge University Press, Cambridge, United Kingdom. 
[Boulanger & Chouteau 2001]BoCh:2001 Boulanger, O. & Chouteau, M., 2001. Constraint in 3D gravity inversion, Geophysical prospecting, 49, 265-280.

[Chung et al. 2008]ChNaOl:2008 Chung, J., Nagy, J. & O'Leary, D. P., 2008. A weighted GCV method for Lanczos hybrid regularization, ETNA, 28, 149-167.
[Chung & Palmer 2015]ChPa:2015 Chung, J. & Palmer, K., 2015. A hybrid
LSMR algorithm for large-scale Tikhonov regularization, SIAM J. Sci. Comput., 37 (5), S562-S580. 


[Erichson et al. 2016]Eri:2016 Erichson, N. B., Voronin, S., Brunton, S. L.  & Kutz, J. N., 2016. Randomized matrix decompositions using R, arXiv preprint arXiv:1608.02148.
[Farquharson 2008]Far:2008 Farquharson, C. G., 2008. Constructing piecwise-constant models in  multidimensional minimum-structure inversions, Geophysics, 73(1), K1-K9.

[Farquharson & Oldenburg 2004]Far:2004 Farquharson, C. G. &  Oldenburg, D. W., 2004. A comparison of Automatic techniques for estimating the regularization parameter in non-linear inverse problems, Geophys. J. Int., 156, 411-425.



[Golub & van Loan 2013]GoLo:2013 Golub, G. H. &  van Loan, C., 2013. Matrix Computations, John Hopkins Press Baltimore,  4th ed.

[Haáz 1953]Haaz Haáz, I. B., 1953. Relations between the potential of the attraction of the mass contained in a finite rectangular prism and its first and second derivatives, Geofizikai Kozlemenyek, 2(7), (in Hungarian). 
[Halko et al. 2011]Halko:2011 Halko, N., Martinsson, P. G., & Tropp, J. A., 2011. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Rev, 53(2), 217-288.







[Kilmer & O'Leary 2001]KiOl:2001 Kilmer, M. E. & O'Leary, D. P., 2001. Choosing regularization parameters in iterative methods for ill-posed problems, SIAM journal on Matrix Analysis and Application, 22, 1204-1221.


[Li & Oldenburg 1998]LiOl:98 Li, Y. &  Oldenburg, D. W., 1998. 3-D inversion of gravity data, Geophysics, 63, 109-119.



[Liu et al. 2013]Liu:13 Liu, S., Hu, X., Liu, T., Feng, J., Gao, W., & Qiu, L., 2013. Magnetization vector imaging for borehole magnetic data based on magnitude magnetic anomaly, Geophysics, 78(6), D429-D444. 
[Liu et al. 2018]Liu:18 Liu, S., Fedi, M., Hu, X., Ou, Y., Baniamerian, J., Zuo, B., Liu, Y., & Zhu, R., 2018. Three-dimensional inversion of magnetic data in the simultaneous presence of significant remanent magnetization and self-demagnetization: example from Daye iron-ore deposit, Hubei province, China. Geophys. J. Int., 215(1), 614-634.









[Martinsson  2016]Martinsson:2016 Martinsson, P. G., 2016. Randomized Methods for Matrix Computations and Analysis of
High Dimensional Data.  arXiv preprint arXiv:1607.01649. <https://arxiv.org/abs/1607.01649>.

[Nabighian et al.  2005]nabighian:2005 Nabighian, M. N., Ander, M. E., Grauch, V. J. S., Hansen, R. O., Lafehr, T. R., Li, Y., Pearson, W. C., Peirce, J. W., Philips, J. D. &  Ruder, M. E., 2005. Historical development of the gravity method in exploration, Geophysics, 70, 63-89.

[Paoletti et al. 2014]PHHF:2014 Paoletti, V., Hansen, P. C.,  Hansen, M. F. & Fedi, M., 2014. A computationally efficient tool for assessing the depth resolution in potential-field inversion,  Geophysics, 79(4), A33-A38.

[Paige & Saunders 1982a]PaSa:1982a Paige, C. C. & Saunders, M. A., 1982a. LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Software, 8, 43-71.
[Paige & Saunders 1982b]PaSa:1982b  Paige, C. C. & Saunders, M. A., 1982b.
 ALGORITHM 583 LSQR: Sparse linear equations and least squares problems, ACM Trans. Math. Software, 8, 195-209.

[Pilkington 2009]Pi:09 Pilkington, M., 2009. 3D magnetic data-space inversion with sparseness constraints, Geophysics, 74, L7-L15.
[Portniaguine & Zhdanov 1999]PoZh:99 Portniaguine, O. & Zhdanov, M. S., 1999. Focusing geophysical inversion images, Geophysics, 64, 874-887

[Rao & Babu 1991]RaBa:91 Rao, D. B. & Babu, N. R., 1991. A rapid
methods for three-dimensional modeling of magnetic anomalies, Geophysics, 56, 1729-1737.
[Renaut et al. 2017]RVA:2017Renaut, R. A., Vatankhah, S. & Ardestani, V. E., 2017. Hybrid and iteratively reweighted regularization  by unbiased predictive risk and weighted GCV for projected systems, SIAM J. Sci. Comput., 39(2), B221-B243. 
[Silva & Barbosa  2006]SiBa:2006 Silva, J. B. C. & Barbosa, V. C. F., 2006. Interactive gravity inversion, Geophysics,  71(1), J1-J9




[Vatankhah et al. 2015]VAR:2015 Vatankhah, S., Ardestani, V. E. &  Renaut, R. A., 2015. Application of the χ^2 principle and unbiased predictive risk estimator for determining the regularization parameter in 3D focusing gravity inversion, Geophys. J. Int., 200, 265-277. 

[Vatankhah et al. 2017]VRA:2017 Vatankhah, S.,  Renaut, R. A.   &  Ardestani, V. E., 2017. 3-D Projected L_1 inversion of gravity data using truncated unbiased predictive risk estimator for regularization parameter estimation, Geophys. J. Int.,  210(3), 1872-1887.
[Vatankhah et al. 2018a]VRA:2018a Vatankhah, S.,  Renaut, R. A.   &  Ardestani, V. E., 2018a. Total variation regularization of the 3-D gravity inverse problem
using a randomized generalized singular value decomposition, Geophys. J. Int., 213, 695-705.  
[Vatankhah et al. 2018b]VRA:2018b Vatankhah, S.,  Renaut, R. A.   &  Ardestani, V. E., 2018b. A fast algorithm for regularized focused 3D inversion of gravity data using the randomized singular-value decomposition, Geophysics, 83(4), G25-G34.  
[Vogel 2002]Vogel:2002 Vogel, C. R., 2002. Computational Methods for Inverse Problems, SIAM Frontiers in Applied Mathematics, SIAM Philadelphia U.S.A.

[Voronin et al. 2015]VMN:2015 Voronin, S.,  Mikesell, D.   &  Nolet, G., 2015. Compression approaches for the regularized solutions of linear systems from large-scale inverse problems, Int. J. Geomath, 6, 251-294.
[Wei et al. 2016]WXZ:2016 Wei, Y.,  Xie, P.   &  Zhang, L., 2016. Tikhonov regularization and randomized GSVD, SIAM J. Matrix Anal. Appl., 37(2), 649-675.

[Xiang & Zou 2013]XiZo:2013 Xiang, H. & Zou, J. 2013. Regularization with randomized SVD for large-scale discrete inverse problems,   Inverse Problems,  29, 085008.
[Xiang & Zou 2015]XiZo:2015 Xiang, H. & Zou, J. 2015. Randomized algorithms for large-scale inverse problems with general Tikhonov regularization, Inverse Problems,  31, 085008.

University of Geosciences (Beijing), China.






 
