# 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:

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:

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:

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}\)):

```
# 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
```