Optimization of gradient algorithm for solving the nonlinear inverse potential problem

Abstract

The paper is devoted to developing the new and memory-efficient parallel algorithms based on the nonlinear conjugate gradient method. The algorithms are aimed to solving the nonlinear inverse potential problem of finding a boundary surface between two layers with constant densities. The algorithms utilize the structure of the Jacobian matrix of integral operator by excluding the less significant elements. The efficient implementation for multicore CPU was developed. Investigation of efficiency and speedup of the parallel algorithm was performed. A model problem with synthetic gravitational data was solved. It was shown that the new algorithms reduce the computation time in comparison with the unmodified one.

Keywords

Potential theory inverse problem nonlinear conjugate gradient method band matrix parallel algorithmPACS: 02.60.Nm 93.85.Hj 07.05.Tp 89.20.Ff

1. Introduction

The problem of finding the interface between two layers with different densities using observed gravitational data arises in studying the structure of the Earth crust. This problem was studied in many works [1, 2, 3].

The problem is described by the integral equation of the first kind; so, it is the ill-posed problem. After the discretization and approximation, this problem is reduced to solving a large ill-conditioned system of nonlinear algebraic equations.

The real observations are performed on the large areas. To increase the accuracy and degree of details, it is essential to use larger grids producing high dimensional data sets. The way to reduce the computation time is to use parallel algorithms and high-performance computing systems.

In works [4, 5], the effective parallel algorithm for solving such inverse problem was constructed and implemented for the multicore CPU using the OpenMP technology. The algorithm is based on the iterative nonlinear conjugate gradient method.

In this paper, we construct and implement modifications of this algorithm based on the structure of Jacobian matrix. The model problem with a large grid is solved. We also investigate efficiency and speedup of the algorithm.

2. Statement of the problem and solving method

The problem is to find the contact surface $\zeta=\zeta\left({x,y}\right)$ , which divides two layers with densities $\sigma_{1}$ and $\sigma_{2}$ (as shown in Fig. 1) using the gravitational field $\Delta g\left({x,y,0}\right)$ measured on the Earth surface, the density contrast $\Delta\sigma=\sigma_{1}-\sigma_{2}$ , and the depth of an asymptotic plane $h$ , i.e., $\mathop{\lim}\limits_{x,y\to\pm\infty}\zeta\left({x,y}\right)=h$ .

This surface should satisfy the following equation:

$\displaystyle\Delta g(x,\;y,\;0)=f\Delta\sigma\int\limits_{-\infty}^{\infty}% \int\limits_{-\infty}^{\infty}\left[\frac{1}{((x^{\prime}-x)^{2}+(y^{\prime}-y% )^{2}+\zeta^{2}(x^{\prime},y^{\prime}))^{1/2}}\right.$

(1) $\displaystyle\left.-\frac{1}{((x^{\prime}-x)^{2}+(y^{\prime}-y)^{2}+h^{2})^{1/% 2}}\right]dx^{\prime}dy^{\prime},$

where $f$ is the gravitational constant.

As a preliminary processing of gravity observations, we use the heightwise transformations technique described in [6]. The technique will allow us to suppress the effects of shallow and deep objects and extract approximately the signal $\Delta g\left({x,y,0}\right)$ of the considered layer $\Pi=\left\{{x,y,z:z_{\min}\leqslant z\leqslant z_{\max}}\right\}$ .

We can rewrite Eq. (2) as

$\displaystyle A\left(\zeta\right)=\Delta g,$ (1-c)

where $A:L_{2}\to L_{2}$ is a Frechét differentiable operator.

Figure 1.
Model of a two-layer medium for the gravimetry problem.

Problem (1-1a) is a nonlinear integral equation of the first kind. This problem belongs to the class of ill-posed problems. After the discretization of the area into $n=M\times N$ grid and approximation of the integral operator using quadrature rules, we obtain the system of nonlinear algebraic equations

$\displaystyle A\left(z\right)=F,$ (2)

where $z$ and $F$ are the $n$ -dimensional vectors.

To solve system Eq. (2), we use the following iterative linearized conjugate gradient method [5]:

$\displaystyle z^{k+1}=z^{k}-\psi\frac{\left\langle{p^{k},S_{\alpha}\left({z^{k% }}\right)}\right\rangle}{\left\|{{A}^{\prime}\left({z^{k}}\right)p^{k}}\right% \|^{2}+\alpha\left\|{p^{k}}\right\|^{2}}p^{k},\;\;p^{k}=S_{\alpha}\left({z^{k}% }\right)+\beta^{k}p^{k-1},\;\;p^{0}=S_{\alpha}\left({z^{0}}\right),$ $\displaystyle\beta^{k}=\max\left\{{0,\frac{\left\langle{S_{\alpha}\left({z^{k}% }\right),S_{\alpha}\left({z^{k}}\right)-S_{\alpha}\left({z^{k-1}}\right)}% \right\rangle}{\left\|{S_{\alpha}\left({z^{k-1}}\right)}\right\|^{2}}}\right\}% ,S_{\alpha}\left(z\right)={A}^{\prime}\left(z\right)^{\text{T}}\left({A\left(z% \right)-F}\right)+\alpha(z-z^{0}),$ (3)

where $z^{k}$ is a solution estimate on the $k$ -th iteration, ${A}^{\prime}\left(z\right)$ is the Jacobian matrix of the operator $A\left(z\right)$ , $0<\psi\leqslant 1$ is the damping factor, and $\alpha>0$ is the regularization parameter.

The condition

$\frac{\left\|{A\left({z^{k}}\right)-F}\right\|}{\left\|F\right\|}\leqslant\varepsilon$

for some sufficiently small $\varepsilon$ is taken as a stopping rule.

As the initial estimation, $z^{0}\equiv h$ is used.
3. Matrix structure and storage algorithm

The variant of optimization of this algorithm is based on the idea used in work [7]. Storing the value of the Jacobian matrix ${A}^{\prime}\left(z\right)$ can be very memory consuming for large grids; hence, it is worthwhile investigating the matrix structure to optimize the storage method.

Let us assume that ${A}^{\prime}\left(z\right)$ is a block matrix. Then, its elements can be defined as

$a_{k,p,l,q}=a_{(k-1)M+p,(l-1)M+q}=f\Delta\sigma\Delta x\Delta y\frac{-z_{(k-1)% M+p}}{\left({\left({x_{k}-x_{l}}\right)^{2}+\left({y_{p}-y_{q}}\right)^{2}+z_{% (k-1)M+p}^{2}}\right)^{3/2}},$

where $k,l=1..M$ are the block indices and $p,q=1..N$ are the indices of the elements inside each block.

The elements $a_{k,p,l,q}$ of the Jacobian matrix ${A}^{\prime}\left(z\right)$ depend on $\left({x_{k}-x_{l}}\right)^{2}+\left({y_{p}-y_{q}}\right)^{2}$ and become smaller as the distance to the main diagonal increases. The structure of the matrix is shown in Fig. 2.

Figure 2.

Structure of the Jacobian matrix ${A}^{\prime}\left(z\right)$ .

The first way to reduce the number of stored elements is to drop off small elements and to approximate the matrix ${A}^{\prime}\left(z\right)$ with a band matrix

$\displaystyle\bar{A}^{\prime}(z)=\left\{{{\begin{array}[]{cc}{a_{ij}}\hfill&{% \left|{i-j}\right|\leqslant\mu,}\\ {0,}\hfill&{\left|{i-j}\right|>\mu,}\\ \end{array}}}\right.$ (4)

where $0<\mu\leqslant n$ is the bandwidth parameter. This approach significantly reduces memory requirements and computation time.

Figure 3 shows the band matrix scheme with respect to the parameter $\mu$ .

To adjust the parameter $\mu,$ we use the error of approximation of the Jacobian matrix:

$\displaystyle\tau=\frac{\left\|{{A}^{\prime}\left(z\right)}\right\|-\left\|{{A% }^{\prime}_{\text{approx}}\left(z\right)}\right\|}{\left\|{{A}^{\prime}\left(z% \right)}\right\|}.$ (5)

We set the parameter $\mu$ in such a way that the resulting error $\tau$ is lower than some threshold.

The second way is to retain the elements of ${A}^{\prime}\left(z\right)$ that correspond some neighborhood of the observation point. For the program implementation, the most adequate criterion for the neighborhood is the chessboard distance; i.e., for some radius $0<\lambda\leqslant n$ , we calculate the elements with indices $(k,l,p,q)$ : $\max(\left|{k-p}\right|,\left|{l-q}\right|)<\lambda$ . We set the parameter $\lambda$ in such a way that the resulting error of approximation $\tau$ is lower than some threshold. For the observation point $(k,l)$ , we calculate the coefficients with indices $(p,q)$ : $\left\{{\max(1,k-\lambda)\leqslant p\leqslant\min(M,k+\lambda),\max(1,l-% \lambda)\leqslant q\leqslant\min(N,l+\lambda)}\right\}$ . The approximate matrix has the following form:

$\displaystyle\overline{\overline{A}}^{\prime}(z)=\left\{{{\begin{array}[]{ll}{% a_{ij},}&{\left({k-\lambda\leqslant p\leqslant k+\lambda}\right)\wedge\left({l% -\lambda\leqslant q\leqslant l+\lambda}\right),}\\ {0,}&{\left({p<k-\lambda}\right)\vee\left({p>k+\lambda}\right)\vee\left({q<l-% \lambda}\right)\vee\left({q>l+\lambda}\right).}\\ \end{array}}}\right.$ (6)

Figure 3.

The numbers show the bandwidth parameter $\mu$ for coefficients in 24 $\times$ 24 matrix ${A}^{\prime}\left(z\right)$ for 6 $\times$ 4 grid.

The scheme of the matrix coefficients with respect to the parameter $\lambda$ is shown in Fig. 4. Each row of the Jacobian matrix corresponds to one observation point in the computational grid. Parameter $\lambda$ is adjusted with respect to the error of approximation $\tau$ of the Jacobian matrix Eq. (5).

Figure 4.

The numbers show the distance $\lambda$ for coefficients in the 24 $\times$ 24 matrix ${A}^{\prime}\left(z\right)$ for 6 $\times$ 4 grid.

Figure 5.

The model (synthetic) gravitational field $\Delta g$ .

This algorithm reduces the number of calculated elements in comparison with the band matrix Eq. (4).

4. Parallel implementation on gpus and numerical experiments

The test problem of finding the interface between two layers for the area 40 $\times$ 45 km ${}^{2}$ is considered.

Figure 5 shows the model gravitational field $\Delta g\left(x,y,\,0\right)$ . This field was obtained by solving the direct problem using the original synthetic surface

$\hat{\zeta}=5-3e^{-(x/15-2.6)^{2}-(y/16.6-2.7)^{2}}-2.2e^{-(x/10-4)-(y/8.6-7)^% {2}}-3e^{-(x/8-6)^{2}-(y/9-4.4)^{2}}.$

This surface is shown in Fig. 6. The noise with 10% peak value and zero mean value was added to the field. Density contrast was $\Delta g=$ 0.1 g/cm ${}^{3}$ . The asymptotic plane was $h=$ 5 km. The grid size was 320 $\times$ 360.

The stopping rule was $\varepsilon=$ 0.08, the regularization parameter was $\alpha=$ 0.1 and the damping factor was $\psi=$ 1.

The problem was solved by three variants of the algorithm:

•
calculating the entire matrix ${A}^{\prime}\left(z\right)$ ;
•
approximation $\bar{A}^{\prime}(z)$ using the band matrix Eq. (4);
•
approximation $\overline{\overline{A}}^{\prime}(z)$ using the neighborhood Eq. (6).

For the band matrix with the bandwidth parameter $\mu=$ 23 and neighborhood with radius $\lambda=$ 4, the error of approximation was $\tau<10^{-15}$ , and the double precision format gives 15 significant digits.

Table 1
Results of the numerical experiments

Method $T_{1}$ , minutes $T_{8}$ , minutes Speedup Efficiency

Full matrix 115 15 7.6 0.96

Band matrix (4) 70 9 7.7 0.97

Neighbourhood (6) 50 6.5 7.7 0.97

Figure 6.
The original surface $\hat{\zeta}\left({x,y}\right)$ .

Figure 7.
The reconstructed surface $\zeta\left({x,y}\right)$ .

Figure 8.
The difference $\hat{\zeta}-\zeta$ between the original and reconstructed surfaces.

Figure 7 shows the reconstructed interface. The solving process took 10 iterations. Resulting relative error is $\|{\hat{\zeta}-\zeta}\|_{\infty}/\|\hat{\zeta}\|_{\infty}<$ 0.02, where the maximum norm is $\left\|\zeta\right\|_{\infty}=\mathop{\max}\limits_{i}\left|{\zeta_{i}}\right|$ . The difference $\hat{\zeta}-\zeta$ between the original and reconstructed surfaces is shown in Fig. 8.

Table 1 shows the execution times $T_{1}$ on one and $T_{8}$ on eight cores of the Intel E5-2650 processor. The third column is the speedup $S_{m}=\frac{T_{1}}{T_{m}}$ , where $T_{1}$ is the execution time for 1 core and $T_{m}$ is the execution time for $m$ cores. The fourth column is the efficiency $E_{m}=\frac{S_{m}}{m}$ .
5. Conclusion

Method	$T_{1}$ , minutes	$T_{8}$ , minutes	Speedup	Efficiency
Full matrix	115	15	7.6	0.96
Band matrix (4)	70	9	7.7	0.97
Neighbourhood (6)	50	6.5	7.7	0.97

The optimized parallel algorithms for solving the nonlinear inverse potential problem of finding a boundary surface between two layers were constructed and investigated. The algorithms are based on the linearized conjugate gradient method.

The main idea of optimizations consists in approximating the Jacobian matrix of the integral operator by excluding the less significant elements and using either the band matrix or neighborhood of the observation point.

The parallel algorithms were developed for solving the inverse gravimetry problem of finding a boundary surface and were implemented for multicore CPU using the OpenMP technology.

A model problem with synthetic gravitational data was solved for a large grid. The experiment shows that the proposed optimizations reduce the computational time by 40% using the band matrix and to 55% using the neighborhood.

The speedup and efficiency of the parallel algorithms were studied. The algorithms show good scaling with 97% efficiency.

Footnotes

Acknowledgments

This work was partially supported by the Center of Excellence “Geoinformation technologies and geophysical data complex interpretation” of the Ural Federal University Program.

References

Numerov

B.V.

, Interpretation of gravitational observations in the case of one contact surface, Doklady Akad. Nauk SSSR. pp. 1930, 569–574,.

Martyshko

P.S.

Pyankov

V.A.

Akimova

E.N.

Vasin

V.V.

and Misilov

V.E.

, On solving a structural gravimetry problem on supercomputer “Uran” for the Bashkir Predural’s area, 12th International Conference on Geoinformatics: Theoretical and Applied Aspects 2013. http://dx.doi.org/10.3997/2214-4609.20142423.

Akimova

E.N.

and Misilov

V.E.

, A fast componentwise gradient method for solving structural inverse gravity problem, 15th International Multidisciplinary Scientific Geoconference and EXPO, SGEM 2015, 3(4), 775–781.

Akimova

E.N.

Martyshko

P.S.

and Misilov

V.E.

, A fast parallel gradient algorithm for solving structural inverse gravity problem, ICNAAM-2014, Rhodos, AIP Conference Proceedings, 1648, 850063. http://dx.doi.org/10.1063/1.4913118

Martyshko

P.S.

Akimova

E.N.

and Misilov

V.E.

, Solving the structural inverse gravity problem by the modified gradient methods, Izvestiya, Physics of the Solid Earth 52(5) (2016), 704–708. http:/dx.doi.org/10.1134/S1069351316050098

Martyshko

P.S.

Fedorova

N.V.

Akimova

E.N.

and Gemaidinov

D.V.

, Studying the structural features of the lithospheric magnetic and gravity fields with the use of parallel algorithms, Izvestiya, Physics of the Solid Earth 50(4) (2014), 508–513. http://dx.doi.org/10.1134/S1069351314040090

Akimova

Miniakhmetova

and Martyshko

, Optimization and parallelization of Newton type methods for solving structural gravimetry and magnetometry inverse problems, 13th International Conference on Geoinformatics – Theoretical and Applied Aspects 2014, 178–182. http://dx.doi.org/10.3997/2214-4609.20140408