The idea of the RFMP is the iterative construction of a sequence of approximations (_Fn)n. This means that we add a basis function d_k from the dictionary to the approximation in each step. This basis function is furthermore equipped with a coefficient α_k.

Since the considered inverse problem is ill-posed, we use the Tikhonov-Phillips regularization, that is, our task is to find a function F which minimizes

That means, if we have the approximation F_n up to step n, our greedy algorithm chooses α_n+1 ∈ ℝ and dn+1∈D such that

is minimized. Here, λ > 0 is the regularization parameter.

We can state the following algorithm for the RFMP.

The maximization, which is necessary to get d_n+1, is implemented by evaluating the fraction for all d∈D in each iteration and picking a maximizer. Since many involved terms can be calculated in a preprocessing, the numerical expenses can be kept low (see [22]). For a convergence proof of the RFMP, we refer to [25]. Briefly, under certain conditions, one can show that the sequence (_Fn)n converges to the solution F_∞ of the Tikhonov-regularized normal equation

where I is the identity operator and F* is the adjoint operator to F.

The ROFMP is an advancement of the RFMP from the previous section.

The basic idea is to project the residual onto the span of the chosen vectors, i.e.,

and then adjust the previously chosen coefficients in such a way that the residual is afterwards contained in the orthogonal complement of the span. Since this so-called backfitting (cf. [14, 15]) might not be optimal, we implement the so-called prefitting (cf. [16]), where the next function and all coefficients are chosen simultaneously to guarantee optimality at every single stage of the algorithm. Moreover, let Wn:=Vn⊥ and let the orthogonal projections on Vn and Wn be denoted by PVn and PWn, respectively. All in all, our aim is to find

Here,

and, thereby, we set

The updated coefficients for the expansion at step n + 1 are given by

The ROFMP is summarized in Algorithm 2.2. For practical details of the implementation, see [13].

Remark 2.1. If we choose di∈D and α_i as in Algorithm 2.2 and update the coefficients as in Equations (13) and (14), we obtain for the regularized case (λ > 0) that R_n is, in general, not orthogonal to Vn for all n ∈ ℕ₀, that means

In [13], it was shown that, with the assumptions from Remark 2.1, there exists a number N: = N(λ) such that

That means we get a stagnation of the residual. This is a problem since the ROFMP is not able to reconstruct a certain contribution of the signal which lies in Vn. Therefore, we modify the algorithm based on iterated Tikhonov-Phillips regularization. That means we run the algorithm for a given number of iterations (in our case K > 0), then interrupt the process and start the algorithm with the latest residual R_K. This action is called restart (or repetition) and it is recurringly executed after every K iterations. In order not to lose information we want to maintain the entire expansion F_n during the process. For this purpose, we need an additional notation: we add a further subscript j to the expansion F_n. Note that we have two levels of iterations here. The upper level is associated to the restart procedure and is enumerated by the second subscript j. The lower iteration level is the previously described ROFMP iteration with the first subscript n. We denote the current expansion by

where F_0,1:= 0 and F_0,j:= F_K,j−1. In analogy to the previous definitions, the residual can be defined in the following way:

That means, after K iterations, we keep the previously chosen coefficients fixed and restart the ROFMP with the residual of the step before. All in all, we have to solve

and update the coefficients in the following way

We summarize for the expansion F_K,m, which is the approximation produced by the ROFMP after m restarts:

In analogy to the RFMP, we obtain a similar convergence result for the ROFMP. That is, under certain technical conditions, the sequence (_Tm)m converges in the Hilbert space H. For further details, we refer to [13, 24].

The Earth Gravitational Model 2008 (EGM2008, see [12]) is a spherical harmonics model of the gravitational potential of the Earth up to degree 2190 and order 2159. We use this model up to degree 100 for the solution F in our numerical tests. For checking the parameter choice methods, we generate different test cases that means test scenarios which vary in the satellite height, the noise-to-signal ratio and the data grid. Based on the chosen function F, our dictionary contains all spherical harmonics up to degree 100 that means our approximation F from the algorithm has the following representation

This is a strong limitation, but higher degrees would essentially enlarge the computational expenses.

Moreover, for the stabilization of the solution, we use the norm of the Sobolev space H:=H((an)n;Ω) which is constructed with

see [26]. This Sobolev space contains all functions F on Ω which fulfill

The inner product of functions F,G∈H is given by

The particular sequence (an)n=((n+12)2)n was chosen, because preliminary numerical experiments showed that the associated regularization term yielded results with an appropriate smoothness.

In our test scenarios, we use a finite set {_{λk}k=1,…,100} of 100 regularization parameters (for details, see Section 4.4). The approximate solution of the inverse problem as an output of the RFMP/ROFMP corresponding to the regularization parameter λ_k and the data vector y is denoted by x_k. This notation is introduced to avoid confusions with the functions F_n which occur at the n-th step of the iteration within the RFMP.

In practice, we deal with noisy data y^ε where the noise level ε is defined by

where l is the length of the data vector y and N2S is called the noise-to-signal ratio. The corresponding result of the RFMP/ROFMP for the regularization parameter λ_k and the noisy data vector y^ε is called xkε. Due to the convergence results for the RFMP/ROFMP (see Equation 8), we introduce the linear regularization operators Rk:ℝl→H,

and assume x_k to be Rky and xkε to be Rkyε, though this could certainly only be guaranteed for an infinite number of iterations.

Due to the importance of the regularization parameter, we summarize in the next section some methods for the choice of this parameter λ. For the comparison of the methods, we have to define the optimal regularization parameter λ_kopt. We do this by minimizing the difference between the exact solution x and the regularized solution xkε corresponding to the parameter λ_k and noisy data.

Then we evaluate the results by computing the so-called inefficiency by

where λ_k* is the regularization parameter selected by the considered parameter choice method. For the computation of the inefficiency, we use the L²(Ω)-norm, since our numerical results led to a better distinction of the different inefficiencies than by using the H-norm. However, the tendency regarding “good” and “bad” parameters were the same in both cases.

In practical cases, the exact solution is usually unknown such that k_opt is a theoretical object only. Parameter choice methods, such as those which are used for the experiments here, provide us with a way to estimate (in some sense) a regularization parameter λ_k* which yields a solution which is expected to be close to the unknown exact solution. The inefficiency quantifies the accuracy of this estimate by comparing xk*ε with the theoretically optimal regularized solution xkoptε. The closer the obtained inefficiency is to 1, the better the parameter choice method performs.

Note that the norms which occur in the several parameter choice methods can be computed with the help of the singular value decomposition. However, we use the singular values of Ψ (see Equations 1 and 3) for this purpose, because the singular value decomposition of F is unavailable. This certainly causes an inaccuracy in our calculations, but appears to be unavoidable for the sake of practicability.

Table 1 shows the different parameter choice methods we tested and introduces their abbreviations. The tuning parameters are chosen in accordance to [18, 19]. For the choice of the maximal index K^, see Section 4.5.

The DP of [27–29] tries to balance the norm of the residual and the noise level ε for a good regularization. It requires knowledge of ε which we provide in our tests. The TDP of [30, 31] tries to deal with the instability of the DP if the noise level is not known. It is designed to work with an estimate ε^ of ε.

The QOC goes back to [32, 33] and its continuous version tries to minimize the norm ‖λddλxλε‖. By using the difference quotient instead of the derivative for our discrete parameters λk=λ0qλk (see Section 4.4), we obtain the method as in Table 1.

The LC (cf. [34–36]) is based on the fact that a log-log plot of (‖Fxkε-yε‖,‖xkε‖) often is L-shaped. The points on the vertical part correspond to under-smoothed solutions, where the horizontal part consists of points due to over-smoothed solutions. This suggests that the “corner” defines a good value for the regularization parameter. Note that usually the method is applied manually and its automating can be difficult.

The EEM of [37, 38] chooses the regularization parameter by minimizing an estimate of the error ‖x-xkε‖ found by an extrapolation procedure. The RM of [39] is based on minimizing a certain weighted form of the norms of the residuals, where the weighting penalizes under-smoothing parameter values.

It is known that the Tikhonov regularized solution for discrete data with independent Gaussian errors can be interpreted as a Bayes estimate of x if x is endowed with the prior of a certain zero mean Gaussian stochastic process (cf. [47]). With this interpretation, the GML estimate was derived by [40] and formulated for stochastic white noise where knowledge of the noise level ε is not needed.

The GCV of [41] uses the following idea: consider all the “leave-one-out” regularized solutions and choose the parameter which minimizes the average of the prediction errors using each solution to predict the data value that was left out. Note that the computation of all these regularized solutions is not necessary for this. A detailed derivation can be found in [47]. Several modifications and extensions of the method such as RGCV (cf. [42, 43]), SRGCV (cf. [44]) and MGCV (cf. [45, 46]) have been developed to deal with instabilities of the GCV. Note that for certain choices of the robustness/stabilization parameter these variants correspond to GCV.

In Sections 2.1 and 2.2, we mentioned that we need to define stopping criteria for our algorithm. We state the following stopping criteria for the RFMP and ROFMP (see also Algorithms 2.1 and 2.2).

In the case of the ROFMP, we choose K = 200 for the restart.

Figure 1 shows two data grids which we use for our experiments. First of all, the Reuter grid (see [48]) is an example of a regular data grid on the sphere. Second, we have a set of irregularly distributed data points on a grid which we refer to as the scattered grid in the following and which was first used in [13]. The latter grid tries to imitate the distribution of measurements along the tracks of a satellite. It possesses additional shorter tracks and, thus, a higher accumulation of data points at the poles and only fewer tracks in a belt around the equator.

For our various scenarios, we get our noisy data if we add white noise to our data values or we add colored noise that is obtained by an autoregression process. Additionally, we test some local noise.

We constructed the admissible values λ_k for the parameter choice as a monotonically decreasing sequence with 100 values from λ₁ = 1 to λ100=10-14 and a logarithmically equal spacing in the following way

Here, λ₀ = 1.3849 and q_λ = 0.7221. The test scenarios are chosen such that the parameter range lies between 1 and 10⁻¹⁴ and includes the optimal parameter away from the boundaries λ₁ and λ₁₀₀. For the choice of the parameters λ_k, we refer to [18, 19].

Most parameter choice methods either increase the index k until a certain condition is satisfied or minimize a certain function for all regularization parameters λ, i.e., after our discretization (see Equation 33) they minimize for all k (see Table 1). For some methods like the quasi-optimality criterion, the values of k have to be constrained by a suitable maximal index K^ which must be chosen such that kopt<K^. To increase computational efficiency, such a maximal index can be used for other methods as well without changing their performance. As in [18, 19], we define this maximal index by

where E‖xk-xkε||2=ε2ρ2(k) is the variance of the regularized solution corresponding to noisy data and ε²ρ²(∞) is its largest value. It is well-known that, in the case of white noise, ρ(k) for the Tikhonov-Phillips regularization is generally given by

Since our singular values occur with a multiplicity of 2n + 1 and we restrict our tests to n = 0, …, 100, the sum above in our tests is given by

For any colored noise, we use the estimate (cf. [19])

with two independent data sets y1ε, y2ε for the same regularization parameter λ_k. Note that xk,1ε, xk,2ε are the regularized solutions corresponding to the parameter λ_k and the noisy data sets y1ε, y2ε.

We can see from Figure 4 that the DP leads to results which are in the range from good to acceptable in all test cases. It yields better results with a more uniformly distributed grid.

The results for the TDP (see Figure 4) for all test cases are rather poor. We can remark that the results get better with a more uniformly distributed data grid. Furthermore, the colored noise leads to slightly bigger boxes than the white noise.

In Figure 5, the inefficiencies of the QOC show that the performance of this method is rather poor. In the case of the Reuter grid, the results reach from good to mediocre in contrast to the scattered grid.

The LC (see Figure 5) yields good results in all test cases. We can remark that there are a few outliers and bigger boxes for the test cases with colored noise and the scattered grid.

The EEM yields acceptable to rather poor results (see Figure 6). We cannot observe any dependency on the grid or the kind of noise related to the acceptable results. Moreover, in the test case with a height of 300 km and an N2S of 5% with colored noise we have some outliers for the RFMP and a large box for the ROFMP.

The results for the RM (see Figure 6) are good to acceptable in all test cases. We only have a few minor outliers.

In Figure 7, we can see that the GML leads to acceptable results only in the case of the Reuter grid. In all cases of the scattered grid, its performance is rather bad.

From Figure 7, we can observe that the GCV yields good results in all test cases. We only have, in the case of the ROFMP, some minor outliers. It yields the best results with a more regularly distributed data grid.

The RGCV yields good to acceptable results (see Figure 8) which get slightly worse and show a larger variance for a higher N2S or colored noise scenarios.

The SRGCV (see Figure 8) has good to acceptable results in all the test cases which are a little bit worse than for the RGCV. The Reuter grid leads to good results whereas the scattered grid seems to be more difficult to handle by the method.

The inefficiencies for the MGCV (see Figure 9) for the test cases with white noise and the Reuter grid are good. In particular, in several of the cases with colored noise the boxes are so big that they partially do not fit in the figure. Obviously, we get here a very large distribution of the inefficiencies. These cases seem to be very hard to handle for this method.

In this section, we show briefly the approximations of the gravitational potential which we obtain by the RFMP and the ROFMP for one typical noisy data set considering a good or a rather poor parameter choice.

For the test case (500 km, 5%, colored noise, scattered grid) with α = 0.54 for the AR(1)-process, Figure 10 shows the approximation which we obtain by the RFMP for the optimal regularization parameter λ₂₉ and the difference to the EGM2008 up to degree 100. Figure 11 shows the approximation belonging to the regularization parameter λ₂₂ which is chosen by the GCV. In Figure 12, we can see the approximation belonging to the parameter λ₄₃ which is chosen by the MGCV. We can see that the MGCV chooses the regularization parameter too small and with this choice we obtain a solution which is underregularized. North-South oriented anomalies occur in the reconstruction which appear to be artifacts due to the noise along the simulated satellite tracks. In contrast, the approximation of the potential for the GCV-based parameter is only slightly worse than the result for the optimal parameter.

Furthermore, we show the same test case as above but with the approximation from the ROFMP with α = 0.56 in the AR(1)-process. Figure 13 shows the approximation for the optimal parameter λ₂₉ and the difference to EGM2008 up to degree 100. In Figure 14, we see the approximation which belongs to the parameter λ₂₂ which is chosen by the GCV. Figure 15 shows the approximation with the regularization parameter λ₉ which is chosen by the GML.

Here, the GML chooses a regularization parameter which is too large, that means our approximation is overregularized. We get less information and details about the gravitational potential. Essential details such as signals due to the Andes or the region around Indonesia occur in the difference plot—much stronglier than for the other examples. Again the parameter choice of the GCV yields a good approximation for the gravitational potential.

Finally, Figures 16, 17 show the difference ‖x-xkε‖L2(Ω) between the original solution (i.e., EGM2008 up to degree 100) and the approximation xkε obtained for the different regularization parameters which were chosen by the considered strategies. The horizontal axis states the index k of the regularization parameter λ_k. The plots refer to the same scenario as Figures 10–15. The arrows show the parameters which are chosen by the methods. The diagrams confirm our observations that the GCV and the LC yield parameters which are closest to the (theoretical) optimal parameter. We obtain almost equally good results for the DP, the RM and the RGCV.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.