Building a Jacobian matrix#

In several applications, like 3D inversions, we need to build the Jacobian matrix (a.k.a sensitivity matrix) of our forward model for every observation point, i.e. how much does our field changes when I apply an infinitesimal change to the physical property of each source in the model.

Let’s take for example the upward component of the gravity acceleration of prisms. Given \(N\) observation points and \(M\) prisms, the gravity field on the \(i\)-th observation point \(\mathbf{p}_i\) can be computed as:

\[g_u(\mathbf{p}_i) = \sum\limits_{j=1}^M G \, \rho_j \, u_j(\mathbf{p}_i)\]

where \(\rho_j\) is the density of the \(j\)-th prism, \(G\) is the Universal Gravitational Constant and \(u_j(\mathbf{p}_i)\) comprehends the analytical solution for the forward modelling of the \(j\)-th rectangular prisms on the \(i\)-th observation point:

\[u_j(\mathbf{p}_i) = \Bigg\lvert \Bigg\lvert \Bigg\lvert k_z(x, y, z) \Bigg\rvert_{X_1}^{X_2} \Bigg\rvert_{Y_1}^{Y_2} \Bigg\rvert_{Z_1}^{Z_2}\]

See also

See notes in choclo.prism.gravity_u for more details on the definition of the kernel function \(k_z\) and the \(x\), \(y\), \(z\) coordinates.

The Jacobian matrix \(\mathbf{J}\) is an \(N \times M\) matrix whose elements are the partial derivatives of \(g_u(\mathbf{p_i})\) with respect to the density of the \(j\)-th prism:

\[J_{ij} = \frac{\partial g_u(\mathbf{p}_i)}{\partial \rho_j} = G \, u_j(\mathbf{p}_i)\]

In most potential field cases, the forward model is linear on the physical property (density in this case), so the Jacobian elements are constant.

In order to build the sensitivity matrix we must need to evaluate \(u_j(\mathbf{p}_i)\) on every observation point and for every prism. We can easily do so with Choclo forward modelling functions, considering that the prisms have unit density.

Let’s write a function that can build a sensitivity matrix for a set of observation points and a set of prisms. Since this operation is as demanding as forward modelling our entire set of prisms on every observation point, we want it to run fast and optionally in parallel. Therefore, we are going to write a Numba function with parallelization enabled:

import numba
import numpy as np
from choclo.prism import gravity_u

@numba.jit(nopython=True, parallel=True)
def build_jacobian(coordinates, prisms):
    """
    Build a sensitivity matrix for gravity_u of a prism
    """
    # Unpack coordinates of the observation points
    easting, northing, upward = coordinates[:]
    # Initialize an empty 2d array for the sensitivity matrix
    n_coords = easting.size
    n_prisms = prisms.shape[0]
    jacobian = np.empty((n_coords, n_prisms), dtype=np.float64)
    # Compute the gravity_u field that each prism generate on every observation
    # point, considering that they have a unit density
    for i in numba.prange(len(easting)):
        for j in range(prisms.shape[0]):
            jacobian[i, j] = gravity_u(
                easting[i],
                northing[i],
                upward[i],
                prisms[j, 0],
                prisms[j, 1],
                prisms[j, 2],
                prisms[j, 3],
                prisms[j, 4],
                prisms[j, 5],
                1.0,
            )
    return jacobian

Note

The numpy.empty function creates an empty array. This means it allocates the memory for this array, but it doesn’t write any values in it. Instead, the jacobian array is filled with garbage values after being initialized.

By using numpy.empty we are saving some time, avoiding to fill the jacobian array with values that we will soon overwrite in the following for loops.

Let’s try this function by defining some prisms and observation points:

easting = np.linspace(-5.0, 5.0, 21)
northing = np.linspace(-4.0, 4.0, 21)
easting, northing = np.meshgrid(easting, northing)
upward = 10 * np.ones_like(easting)

coordinates = (easting.ravel(), northing.ravel(), upward.ravel())

prisms = np.array(
    [
        [-10.0, 0.0, -7.0, 0.0, -15.0, -10.0],
        [-10.0, 0.0, 0.0, 7.0, -25.0, -15.0],
        [0.0, 10.0, -7.0, 0.0, -20.0, -13.0],
        [0.0, 10.0, 0.0, 7.0, -12.0, -8.0],
    ]
)

And run it:

jacobian = build_jacobian(coordinates, prisms)
jacobian

array([[-4.49966911e-11, -4.76014375e-11, -3.80054986e-11,
        -2.83231448e-11],
       [-4.49659418e-11, -4.75828152e-11, -3.86898648e-11,
        -2.90433576e-11],
       [-4.48738920e-11, -4.75270202e-11, -3.93578743e-11,
        -2.97540033e-11],
       ...,
       [-3.19387365e-11, -4.58884222e-11, -4.10526880e-11,
        -4.48751701e-11],
       [-3.12924999e-11, -4.52460654e-11, -4.11114162e-11,
        -4.49880072e-11],
       [-3.06338007e-11, -4.45849449e-11, -4.11310225e-11,
        -4.50257165e-11]])

Warning

Jacobian matrices can be very big. Large number of observation points and sources can lead to Jacobian matrices that cannot fit in the available memory of your system.

Now that we have defined our Jacobian matrix, we can use it to forward model the gravity field of our prisms on every observation point by just computing a dot product between it and the density vector of the prisms (\(\mathbf{m}\)):

\[\begin{split}\mathbf{g_u} = \begin{bmatrix} g_u({\mathbf{p}_1}) \\ \vdots \\ g_u({\mathbf{p}_N}) \\ \end{bmatrix} = \begin{bmatrix} J_{11} & \cdots & J_{1M} \\ \vdots & \ddots & \vdots \\ J_{N1} & \cdots & J_{NM} \end{bmatrix} \cdot \begin{bmatrix} \rho_1 \\ \vdots \\ \rho_M \\ \end{bmatrix} = \mathbf{J} \cdot \mathbf{m}\end{split}\]

# Define densities for the prisms
densities = np.array([200.0, 300.0, -100.0, 400.0])

# Compute result
g_u = jacobian @ densities

Note

The @ operator performs a matrix product. It’s a shorthand of the numpy.matmul function.

We can check that this result is right by comparing it with the output of the gravity_u_parallel function we defined in the How to use Choclo:

expected = gravity_upward_parallel(coordinates, prisms, densities)
np.allclose(g_u, expected)

True

Download Python script

Download Jupyter notebook