# Robust 3D Joint Inversion of Gravity and Magnetic Data: A High-Performance Computing Approach
## 1 Introduction
High-performance computing (HPC) techniques to obtain solutions to inversion problems have been used since the late 1980s, as mentioned before [54]. However, starting in the second decade of the 2000s, a new era of HPC has developed thanks to the arrival of the new graphic processing units (GPU), its high parallelization capacity, and easy access outside the business or academic field due to its low cost. The most commonly used techniques in geophysics have focused mainly on electromagnetic and seismic methods and can be generally classified as parallelization of numerical methods [55; 56; 57; 58; 59; 60; 61], based on neural networks and deep learning [62; 63; 64; 65], and more recently, based on quantum computing [66]. In this case, the variational modeling approach has been employed to construct an objective function or parametric functional consisting of (1) a mismatch term measuring the distance between \(F(x)\) and \(y\), usually via the L2 norm, and (2) a regularization functional favoring appropriate minimizers and penalizing potential solutions with undesirable structures [67] using the form:  
\[P^{a}(x)=\|F(x)-y\|^{2}+\alpha R(x) \tag{1}\]  
The regularization functional \(R(x)\) represents a family of linear operators that should stabilize the problem and bound desirable solutions within the infinite space of possible solutions, while \(\alpha\) is a scalar called the regularization parameter that gives weight to the regularization functional. The higher this value, the more influence it has on the solution; conversely, the closer to zero this value is, the closer it is to a solution without regularization. The purpose of regularization tools is to provide the user with easy-to-use routines based on numerically robust and efficient algorithms for doing experiments with regularization of discrete, ill-posed problems [68]. Two widely used regularizers are the Laplacian, which adds smoothness to the models because it is known that (general) changes in properties occur gradually; therefore, smoothed distributions are expected. On the other hand, sometimes there is subsurface information obtained directly or indirectly in previous studies, and the solution should incorporate such information, so an _a priori_ model is used as a second regularizer. In addition, it is possible to know the predominant orientations of structures in the study area and some directional preference regularizer may be useful; therefore, the implementation of such regularizers allows robust solutions capable of fitting a wide variety of needs according to geological-structural context. Joint inversion as a regularizer to mitigate the non-uniqueness of solutions has proven to be very effective [69; 70; 71; 72; 73; 74; 75; 76]. This assumes that information obtained from different methods, whether directly or indirectly, comes from the same source interacting with different physical fields. Therefore, it is reasonable to expect inversion results where the distribution of different physical properties in the subsurface of the same study area have some relationship, as well as to propose the generation of models that feedback through each other through correspondence relationships that can be (1) intrinsic to rock physics or (2) structural. [PERSON] et al. [73] shows that in a Gramian space, the norm is a function that provides a measurement of the correlation between that function and the different parameters of the model within a Gram matrix. Due to its quadratic form similar to the L2 norm, it can be used as a regularizer of structural similarity, with the advantage that this function is differentiable, thus avoiding resorting to some linearization as in the case of crossed gradients [41]. To perform a 3D inversion of multiple types of geophysical data together is a task with an extremely high computational cost that grows exponentially as the resolution of the solutions increases.