In this study, a GPU-enabled implementation of an algorithm for the solution of a real-life optimization problem arising in the geoacoustic inversion is proposed. In the inversion algorithm, a single-hydrophone recording of a pulse acoustic signal is used for the estimation of the waveguide parameters from the dispersion data. We devised and implemented a hierarchical optimization scheme where and efficient GPUenabled implementation of the objective function calculation is combined with the CPU implementation of the search algorithm. The GPU-enabled implementation turned out to be several dozen times faster than the pure CPU one. Exploiting this advantage we were able to successfully solve the considered problem in reasonable time.
Black-box optimization methods can be employed to solve real-life problems from different areas
Let us give a brief outline of the paper. In Section 2 we describe an black-box optimization algorithm, used to solve the problem. In Section 3 we describe a problem of geoacoustic inversion. In Section 4 we describe two implementations of the proposed algorithm, and show the results of computational experiments on the test problem. In the rest of the paper we draw conclusions. 2
In black-box optimization problems derivatives are either unavailable or very hard to compute. Thus, solution methods should not rely on derivatives. There are a lot of local search techniques which can be used for this purpose: Hooke-Jeeves method [Hooke & Jeeves, 1961], pseudo-gradient approaches [Evtushenko et al., 2016] and a variety of coordinate descent techniques. We developed a simple modification of the coordinate descent method presented in Fig. 1. 1 while maxi=1;:::;n i > " do 2 for i = 1 to n do 3 x′ = x + iei 4 if f (x′) < f (x) then 5 x := x′ 6 i := i
The algorithm suggested in Section 2 is now applied to the solution of the following (synthetic) problem of geoacoustic inversion.
Consider a homogeneous 2D geoacoustic waveguide f(r; z)jz > 0g [Jensen et al., 2011], where the sound speed
and density depend on the depth z but do not depend on r
Assume that a pulse acoustic signal is emitted by a point source located at R = 0; z = zs and
received by a hydrophone at the distance R from the source. The exact value of the receiver range
R is also considered unknown. Performing some time-frequency analysis of the received signal
Using the mode theory of sound propagation in shallow-water waveguides we can also compute the theoretical arrival times th(f; m) for any given set A = (cb; b; R) of the values of the unknown parameters. Introducing a fitness that quantifies the agreement between the theoretical and experimental arrival times E(A) = ∑M m=1 ∑Nm n=1
mth(fn; A)
∑M
m=1 Nm
mexp(fn) 2
;
(1)
we turn the geoacoustic inversion problem into a global minimization problem
For our test case we simulated the signal at the receiver using the mode theory, then we applied a warping-based algorithm to the resulting time series, and subsequently used the resulting syntetic dispersion data exp(f; m) as the input for our inversion algorithm. The sound speed profile in the water column and the pulse signal emitted by the source was the same as in [Zaikin & Petrov, 2017]. The true values of the unknown parameters and the search space (a cuboid) are presented in Table 1. 4
For many classes of computational problems, graphics processing units (GPUs) show significant speedup over
central processing units (CPUs). In particular, GPUs can help to solve global optimization problems
In our computational experiments we considered the geoacoustic inversion problem described in Section 3. So, the ASN search was applied for minimizing the objective function (1).
The GPU programming model is based on data parallelism. To achieve peak efficiency, modern GPUs should simultaneously run about 10 000 computational threads [Corporation, 2017, Bulavintsev, 2015]. In the GPUbased implementation of our search procedure, we exploit data parallelism found in the parallel calculation of
modal group velocities for a single point. Level 1 sub-algorithm executes on the CPU, while level 2-4 subalgorithms execute on the GPU (Table 2). A typical residue calculation for a single search space point requires calculation of modal group velocities for thousands of frequencies. Thus, data parallelism exposed in this way should be enough to efficiently exploit modern GPUs.
Traditionally, geoacoustics computations are performed in double precision floating point arithmetics (FP64). However, modern consumer-grade GPUs suffer a great (for some devices, up to 32 times [Corporation, 2017]) drop in performance when using FP64 instructions instead of their single-precision (FP32) counterparts. This fact led us to investigate the consequences of changing the algorithm to FP32. We considered the model problem of geoacoustic inversion with 3 parameters (R; b; cb) (see Section 3). The results are presented in Table 3. All experiments were conducted on the Intel Core i7 930 CPU and the Nvidia GTX 1050 GPU, with the use of CUDA 8.0 SDK and GCC 5.4.0 compiler. The NVCC compiler was forced to produce native Compute Capability 6.1 code. Both NVCC and GCC compilers were run with -O3 optimization flag.
The most computationally intensive part of our algorithm is the bisection sub-algorithm at level no. 4 in Table 2. We implemented it as described in [Demmel, 1997]) and [Corporation, 2017], producing single- and double-precision versions, algorithmically identical on both GPU and CPU platforms. For reference, we provide the double-precision version based on the AlgLib library [Bochkanov & Bystritsky, 2016], that we used in our previous work [Zaikin et al., press]. The AlgLib library includes a state-of-the-art implementation of the bisection algorithm, thoroughly tuned to produce the most accurate results possible with the modern CPUs floating point units (FPUs). Our CPU-based implementation of bisection algorithm can not boast such accuracy. However, it is notably faster than AlgLib due to its simplicity. That is why, our and AlgLib-based algorithm’s results are different. The discrepancy between the outputs (residue) of the same algorithm on the CPU and GPU is the result of the different implementations of floating-point units on these platforms.
Table 3 clearly indicates the superiority of the GPU performance, even in FP64. The FP32 mode of the GPU is 150 times faster than the FP64 mode, though the former’s final residue (0.0215154) is much worse than the latter one’s (0.0131284). However, recomputation of this final point in FP64 gives the result (0.0134315, see Table 4) that is on par with that of the FP64 CPU implementation.
This leads us to suggest the mixed-mode, two-stage algorithm, that uses the ”coarse” FP32 computation mode
during Stage I, and then refines its result in Stage II by resuming the search in the ”fine” FP64 computation
mode. The performance of this algorithm is presented in Table 4. For reference, the ”Initial residue” column of
Table 4 shows residue of the best point found by Stage I recomputed by various Stage II algorithms. It is easy to
see that the GPU-based FP32 algorithm for Stage I is best complemented for Stage II by the GPU-based FP64
algorithm. The total runtime of this mixed-mode algorithm is 127 seconds. This is faster than runtime of the
pure GPU-based FP64 algorithm (777 seconds), and its final residue (0.0130635) is almost equal to that of the
pure CPU-based FP64
The FP64 performance of our GPU-based algorithm is only 2,5 times lower than the performance of our FP32 algorithm on the same GPU (see Table 3). However, according to [Corporation, 2017] this performance drop should be at least 4x times for the GTX1050 GPU used in our experiments. This means that the GPU performance is likely to be bottlenecked by suboptimal multiprocessor occupancy [Corporation, 2017] and/or excessive memory transfers, rather than the FP64 FPUs of the device. The higher GPU performance we perceived in our previous experiments with the parallelization scheme based on the simultaneous computing of residues for many search space points supports this theory. This means that a thorough optimization of the GPU code could increase the GPU performance of our geoacoustic inversion algorithm at least 2-3 times.
The source code of our application is available online1. 5
In the present paper, we suggested a modification of a coordinate descent algorithm. This algorithm was implemented in two versions. The first one is designed for launching on a CPU. In the second one, the algorithm itself works on a CPU, while the calculation of the objective function is performed on a GPU. Both versions were applied for solving one synthetic problem of geoacoustic inversion. Being little less accurate (in the sense of found solutions), the GPU-based version turned out to be much faster, than the CPU one. Acknowledgements This study was partially supported by the Council for Grants of the President of the Russian Federation (grants No. MK-2262.2017.5, No. NSh-8081.2016.9, No. NSh-8860.2016.1), the Russian Foundation for Basic research (grants No. 16-05-01074 a, No. 16-07-00155 a and No. 17-07-00510 a), Presidium of RAS programs I.33, I.5, and the POI FEB RAS Program “Nonlinear dynamical processes in the ocean and atmosphere”. 1https://github.com/mposypkin/acouwater/tree/f gpu [Demmel, 1997] Demmel, J. W. (1997). Applied numerical linear algebra. SIAM.