Massively Parallel Lagrangian Relaxation Algorithm for Solving Large-Scale Spatial Optimization Problems Using GPGPU

Lagrangian Relaxation (LR) is an effective method for solving spatial optimization problems in geospatial analysis and GIS. Among others, it has been used to solve the classic p-median problem that served as a unified local model in GIS since the 1990s. Despite its efficiency, the LR algorithm has seen limited usage in practice and is not as widely used as off-the-shelf solvers such as OPL/CPLEX or GPLK. This is primarily because of the high cost of development, which includes (i) the cost of developing a full gradient descent algorithm for each optimization model with various tricks and modifications to improve the speed, (ii) the computational cost can be high for large problem instances, (iii) the need to test and choose from different relaxation schemes, and (iv) the need to derive and compute the gradients in a programming language. In this study, we aim to solve the first three issues by utilizing the computational power of GPGPU and existing facilities of modern deep learning (DL) frameworks such as PyTorch. Based on an analysis of the commonalities and differences between DL and general optimization, we adapt DL libraries for solving LR problems. As a result, we can choose from the many gradient descent strategies (known as “optimizers”) in DL libraries rather than reinventing them from scratch. Experiments show that implementing LR in DL libraries is not only feasible but also convenient. Gradient vectors are automatically tracked and computed. Furthermore, the computational power of GPGPU is automatically used to parallelize the optimization algorithm (a long-term difficulty in operations research). Experiments with the classic p-median problem show that we can solve much larger problem instances (of more than 15,000 nodes) optimally or nearly optimally using the GPU-based LR algorithm. Such capabilities allow for a more fine-grained analysis in GIS. Comparisons with the OPL solver and CPU version of the algorithm show that the GPU version achieves speedups of 104 and 12.5, respectively. The GPU utilization rate on an RTX 4090 GPU reaches 90%. We then conclude with a summary of the findings and remarks regarding future work.