Simulations#

Generally speaking, any inversion require at least one forward modelling function, i.e. a function that can predict data based on a given model. It is usually represented by a function \(\mathbf{f}: \mathbb{R}^M \rightarrow \mathbb{R}^N\), and it’s usually incorporated to the inversion through the data misfit term \(\phi_\text{d}(\mathbf{m})\):

\[\phi_\text{d}(\mathbf{m}) = \sum\limits_{i=1}^N \frac{ \left\lvert f_i(\mathbf{m}) - d_i^\text{obs} \right\rvert^2 }{ \epsilon_i^2 }\]

where \(\mathbf{m}\) is the model vector, \(d_i^\text{obs}\) is the \(i\)-th observed datum, \(f_i(\mathbf{m})\) is the \(i\)-th element of the \(\mathbf{f}(\mathbf{m})\) vector, and \(\epsilon_i\) is the uncertainty of the \(i\)-th datum.

The inversion will require the ability to evaluate that function \(\mathbf{f}\) on any model vector \(\mathbf{m}\), and also its Jacobian matrix \(\mathbf{J}_\mathbf{f}(\mathbf{m})\), i.e. the matrix of first derivatives:

\[\begin{split}\mathbf{J}_\mathbf{f}(\mathbf{m}) = \begin{bmatrix} \frac{\partial f_1(\mathbf{m})}{\partial m_1} & \dots & \frac{\partial f_1(\mathbf{m})}{\partial m_M} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_N(\mathbf{m})}{\partial m_1} & \dots & \frac{\partial f_N(\mathbf{m})}{\partial m_M} \end{bmatrix}\end{split}\]

where \(\frac{\partial f_i(\mathbf{m})}{\partial m_j}\) is the derivative of the \(i\)-th element of the forward modelling function with respect to the \(j\)-th element of the model vector evaluated on a given model \(\mathbf{m}\).

In the inversion framework the forward modelling functions are represented through Simulation classes. These classes always provide two attributes:

n_params: the number of elements in the expected model array.
n_data: the number of elements of the array with predicted data.

and two methods:

__call__(): evaluates the function \(\mathbf{f}\) on the provided model array (with n_params elements). It will always return an array with n_data elements.
jacobian(): returns the Jacobian matrix of the function \(\mathbf{f}\) evaluated on the provided model array (with n_params elements). It can either return a dense or sparse 2D array (with shape (n_data, n_params)), or a LinearOperator object.

For example, let’s define a inversion_ideas.LinearRegressor simulation with a random matrix:

import inversion_ideas as ii

n_data = 5
n_params = 10
simulation = ii.LinearRegressor.create_random(n_data, n_params, seed=42)

We can access the n_params and n_data attributes of the simulation:

simulation.n_data, simulation.n_params

(5, 10)

We can call the simulation with a given model vector to evaluate the forward modelling function:

import numpy as np

rng = np.random.default_rng(seed=43)
model = rng.uniform(low=-1.0, high=1.0, size=n_params)

simulation(model)

array([-0.68606714, -0.07853333,  0.95821454, -0.50970221, -1.40535098])

We can also get its Jacobian matrix for the same model:

jacobian = simulation.jacobian(model)
jacobian

array([[ 0.5479121 , -0.12224312,  0.71719584,  0.39473606, -0.8116453 ,
         0.9512447 ,  0.5222794 ,  0.57212861, -0.74377273, -0.09922812],
       [-0.25840395,  0.85352998,  0.28773024,  0.64552323, -0.1131716 ,
        -0.54552256,  0.10916957, -0.87236549,  0.65526234,  0.2633288 ],
       [ 0.51617548, -0.29094806,  0.94139605,  0.78624224,  0.55676699,
        -0.61072258, -0.06655799, -0.91239247, -0.69142102,  0.36609791],
       [ 0.48952431,  0.93501946, -0.34834928, -0.25908059, -0.06088838,
        -0.62105728, -0.74015699, -0.04859015, -0.5461813 ,  0.33962799],
       [-0.12569616,  0.66535639,  0.4005302 , -0.37526672,  0.6645196 ,
         0.60952871, -0.22504324, -0.42334379,  0.36499101, -0.72049503]])

jacobian.shape

(5, 10)

Alternatively, we could configure the LinearRegressor to return a LinearOperator as a Jacobian. In this example, we need to pass build_jacobian=False as argument while creating the simulation:

simulation = ii.LinearRegressor.create_random(n_data, n_params, seed=42, build_jacobian=False)

Hint

The build_jacobian option in LinearRegressor forces its jacobian() method to return a LinearOperator for learning and testing purposes, but it will still store the whole matrix in memory under the hood.

Let’s now generate its Jacobian matrix for a given model:

jacobian = simulation.jacobian(model)
jacobian

<5x10 MatrixLinearOperator with dtype=float64>

As we can see it no longer returns a dense matrix, but a MatrixLinearOperator, which is a child class of LinearOperator. We can use this object to compute dot products with any other array, as if it were a regular matrix:

vector = rng.random(size=n_params)
jacobian @ vector

array([ 0.41488308,  0.61800282, -0.32520017,  0.55317812,  0.21128147])

We can also compute the dot product of its transpose against a vector:

vector = rng.random(size=n_data)
jacobian.T @ vector

array([ 0.66408941,  2.12031029,  1.24977814,  0.64936908,  0.43495609,
       -0.3933718 , -0.60517756, -1.5562028 , -0.30267928,  0.01912028])

Important

The Jacobian matrix of any simulation will contain n_data * n_params elements. For large problems with high n_params and/or high n_data, those matrices can grow to big that they don’t fit in memory. In those cases it’s preferrable to work with LinearOperator objects that are capable of computing dot products on the fly, without the need to store all the matrix elements in memory.

The capability of each simulation my depend on its implementation, and users are encouraged to read the specific documentation of each one of them.