## 3 Methodology
The minimization of 3 with respect to the density model will result in the system  
\[\left[\mathbf{A}^{\text{T}}\mathbf{C}_{\mathbf{d}}^{\text{T}}\mathbf{C}_{ \mathbf{d}}\mathbf{A}+\alpha^{2}\mathbf{D}^{\text{T}}\mathbf{D}+\mathbf{C}_{ \mathbf{m}}^{\text{T}}\mathbf{C}_{\mathbf{m}}\right]\mathbf{m}^{\rho}=\mathbf{ A}^{\text{T}}\mathbf{C}_{\mathbf{d}}^{\text{T}}\mathbf{C}_{\mathbf{d}} \mathbf{d}^{\text{obs}}+\mathbf{C}_{\mathbf{m}}^{\text{T}}\mathbf{C}_{\mathbf{ m}}\mathbf{m}_{\text{apr}}^{\text{ge}}, \tag{4}\]  
where \(\mathbf{D}\) is the discrete operator of \(\
abla\) and \(\mathbf{A}\) is the density kernel equal to the derivative of Equation (2) with respect to \(\rho\) (linear). For depth estimation of the model, the model parameter \(\mathbf{m}\equiv\mathbf{m}^{z}\), where \(z\) is the depth of the prism (Figure 2b) and \(z\) does not denote exponent. By inspecting Equation (2), it is evident that the inversion by depth is a non-linear problem. This issue is solved by performing a linearization around a starting model \(\mathbf{m}_{0}^{\text{z}}\), resulting in the following system:  
\[\left[\mathbf{J}^{\text{T}}\mathbf{C}_{\mathbf{d}}^{\text{T}}\mathbf{C}_{ \mathbf{d}}\mathbf{J}+\alpha^{2}\mathbf{D}^{\text{T}}\mathbf{D}+\mathbf{C}_{ \mathbf{m}}^{\text{T}}\mathbf{C}_{\mathbf{m}}\right]\mathbf{m}^{z}=\mathbf{J} ^{\text{T}}\mathbf{C}_{\mathbf{d}}^{\text{T}}\mathbf{C}_{\mathbf{d}}(\mathbf{ d}^{\text{obs}}+\mathbf{J}\mathbf{m}_{0}^{z})+\mathbf{C}_{\mathbf{m}}^{\text{T}} \mathbf{C}_{\mathbf{m}}\mathbf{m}_{\text{apr}}^{\text{z}}. \tag{5}\]  
In this case, the matrix \(\mathbf{J}\) is the Jacobian corresponding to the derivative of Equation (2) with respect to \(z\) evaluated on the starting model \(\mathbf{m}_{0}^{z}\), with this, Frechet derivatives are avoided, e.g, methods such as perturbations and sensitivity approaches, which can be computationally expensive. The solution of Equations (4) and (5) can be achieved by inverting the symmetric square matrix of the left size, commonly called \(\mathbf{A}^{\text{T}}\mathbf{A}\) inversion. However, this implementation requires the storage and computation of such a matrix, which results in a large computational time, and sometimes the memory is surpassed for large scale problems. One solution for this problem is to perform conjugate gradient least squares to solve a system in the form of \(\mathbf{G}\mathbf{m}=\mathbf{d}\). This method internally computes the  
Figure 2: (a) Single prism of constant density in Cartesian coordinates for gravity measurements and (b) an example of \(n_{x}\times n_{y}\times 2=64\) juxtapositioned prisms for depth estimation. square matrix of Equation (4), allocating only the matrix \(\mathbf{A}\) and the sparse matrix \(\mathbf{D}\), \(\mathbf{C}_{\mathbf{m}}\), and \(\mathbf{C}_{d}\). For example, for depth inversion, the system \(\mathbf{G}\mathbf{m}=\mathbf{d}\) is solved as  
\[\left[ \begin{array}{c}\mathbf{C}_{\mathbf{d}}\mathbf{J}\\ \mathbf{\alpha}_{\mathbf{m}\mathbf{g}}\mathbf{D}\\ \mathbf{C}_{\mathbf{m}}\end{array} \right]\mathbf{m}^{z}=\left[ \begin{array}[] {c}\mathbf{C}_{\mathbf{d}}(\mathbf{d}^{\text{obs}}+\mathbf{J}\mathbf{m}_{0}^ {z})\\ \mathbf{0}\\ \mathbf{C}_{\mathbf{m}}\mathbf{m}_{\text{apr}}^{z}\end{array} \right], \tag{6}\]  
In this case, the augmented matrix \(\mathbf{G}\) were large and sparse due to the discrete operations for the Tikhonov regularization, and the model vector \(\mathbf{m}^{\rho}\) was not modified. Notice that multiplication of Equation (6) by \(\mathbf{G}^{T}\) yields to Equation (5). The algorithm to solve the system \(\mathbf{G}\mathbf{m}=\mathbf{d}^{\text{CG}}\) is presented in Algorithm 1, modified from Nocedal & Wright [34]. Notice that an efficient implementation of this algorithm requires that all the matrices be stored in a sparse representation. We used coordinate list (COO) for storing such sparse matrices, given that the matrix and its transpose have to be accessed simultaneously [35].