Objective function#
Defining objective functions#
Any objective function is defined as a scalar field \(\phi: \mathbb{R}^M \rightarrow \mathbb{R}\), i.e. it takes a model vector \(\mathbf{m}\) with \(M\) elements and returns a single real value.
In the new SimPEG’s inversion framework, objective functions are represented by a child of the inversion_ideas.base.Objective, like the inversion_ideas.DataMisfit for example.
import numpy as np
import inversion_ideas as ii
# Define a linear regressor as a simulation
n_data, n_params = 3, 5
simulation = ii.simulations.LinearRegressor.create_random(n_data, n_params, seed=42)
# Create some synthetic data and uncertainties
true_model = np.array([-5., 2., -3., 4., 1.])
std = 1e-2
data = simulation(true_model) + np.random.default_rng(seed=42).normal(size=n_data, scale=std)
uncertainties = std * np.ones_like(data)
# Define a data misfit function
data_misfit = ii.DataMisfit(data, uncertainties, simulation)
data_misfit
Once defined, the objective functions can only be evaluated on model vectors with n_params elements. We can check that property on every objective function:
print("Expected size of model vector:", data_misfit.n_params)
Expected size of model vector: 5
Objective functions can be evaluated by calling them:
# Define some model
model = np.array([-2., -1., 0., 1., 2.])
# Evaluate the data_misfit on that model
data_misfit(model)
np.float64(516777.97897255165)
All objective functions have two methods to compute its derivatives:
gradient()computes the gradient of the objective function: \(\nabla \phi(\mathbf{m})\).hessian()computes the Hessian of the objective function, i.e. the matrix of the second derivatives: \(\nabla^2 \phi(\mathbf{m})\).
# Get gradient
gradient = data_misfit.gradient(model)
gradient
array([ 145012.76749299, -28865.37179219, 63841.4074464 ,
-117950.76338047, -35073.33393289])
# Get Hessian
hessian = data_misfit.hessian(model)
hessian
<5x5 _ProductLinearOperator with dtype=float64>
Important
By default, the DataMisfit and many other
objective functions in the framework won’t build the whole Hessian matrix, but
return a LinearOperator instead.
Hessian matrices are n_params \(\times\) n_params square matrices
that can grow very large when n_params is big.
A LinearOperator allows us to operate with such
matrix as if it were an array, but without storing it entirely in memory.
For example, we can compute the dot product of the hessian with the
model vector:
hessian @ model
array([-88614.46048782, -28378.74681096, -69639.8526546 , 33764.94681921,
42242.44031587])
Renaming a function#
An objective function can be renamed using the set_name() method. For example:
data_misfit.set_name("data_misfit")
data_misfit
Let’s rename it to its previous name:
data_misfit.set_name("d")
data_misfit
Tip
The set_name() method returns the same objective function, so we can use it directly when defining a new one. For example:
data_misfit = ii.DataMisfit(data, uncertainties, simulation).set_name("data_misfit")
Scaling a function#
We can easily multiply an Objective function by a scalar, obtaining a inversion_ideas.base.Scaled object.
The inversion_ideas.base.Scaled is another type of objective functions that represents an objective function (encapsulated in the function attribute) multiplied by a scalar (encapsulated in the multiplier attribute).
For example, we can build a new scaled function by multiplying the data_misfit by a float:
scaled = 5.0 * data_misfit
scaled
We can still access the underlying objects:
scaled.multiplier
5.0
scaled.function
Important
The function variable stores a reference to the original objective function. The
Scaled doesn’t perform any copy of it. Therefore scaled.function is the data_misfit object:
scaled.function is data_misfit
True
We can evaluate the scaled function and get its gradient and Hessian as with any other objective function:
scaled(model)
np.float64(2583889.894862758)
scaled.gradient(model)
array([ 725063.83746496, -144326.85896097, 319207.03723201,
-589753.81690234, -175366.66966446])
scaled.hessian(model)
<5x5 _ScaledLinearOperator with dtype=float64>
Combining functions#
Add two functions#
It’s possible to define any linear combination of objective functions.
Consider we have a SimpleSmallness regularization function:
smallness = ii.SimpleSmallness(n_params)
smallness
We can add it to the data_misfit one by summing them together to obtain a new objective function:
phi = data_misfit + smallness
phi
This inversion_ideas.base.Combo is an objective function that works as its mathematical counterpart \(\phi(\mathbf{m}) = \phi_\text{d}(\mathbf{m}) + \phi_\text{0}(\mathbf{m})\) does: we can evaluate it on a given model, get its gradient and Hessian:
phi(model)
np.float64(516787.97897255165)
phi.gradient(model)
array([ 145008.76749299, -28867.37179219, 63841.4074464 ,
-117948.76338047, -35069.33393289])
phi.hessian(model)
<5x5 _SumLinearOperator with dtype=float64>
Important
The inversion_ideas.base.Combo class is not meant to be instantiated.
Add two or more objective functions together to obtain one.
The inversion_ideas.base.Combo works as a collection of objective functions. This means that we can get the number of objective functions it contains through the len() function, and we can also access each one of its terms by index:
len(phi)
2
print(phi[0])
print(phi[1])
φd(m)
φ0(m)
for i, term in enumerate(phi):
print(f"Term {i}: {term}")
Term 0: φd(m)
Term 1: φ0(m)
print(data_misfit in phi)
print(smallness in phi)
True
True
Add multiple functions#
We can also add more than two objective functions together. For example, consider that we want to add the data_misfit and the smallness regularization with a SimpleFlatness regularization:
flatness = ii.SimpleFlatness(n_params).set_name("1")
flatness
phi = data_misfit + smallness + flatness
phi
By default, when adding more than two objective functions together, we’ll get a nested inversion_ideas.base.Combo object: the first element will be a inversion_ideas.base.Combo containing the data_misfit and the smallness regularization, and the second element will be the flatness regularization.
print(len(phi))
print(phi[0])
print(phi[1])
2
φd(m) + φ0(m)
φ1(m)
And since the first element is a Combo, we can also treat it as a collection:
print(len(phi[0]))
print(phi[0][0])
print(phi[0][1])
2
φd(m)
φ0(m)
If we don’t want to have this kind of structure, we can
flatten any Combo through the
flatten() method:
flat_phi = phi.flatten()
flat_phi
len(flat_phi)
3
print(flat_phi[0])
print(flat_phi[1])
print(flat_phi[2])
φd(m)
φ0(m)
φ1(m)
data_misfit in flat_phi
True
smallness in flat_phi
True
flatness in flat_phi
True
Is in or contains?#
When using the in statement to check whether an objective function is part of a Combo, we are only checking if that function is one of its elements, but not a recursive search through its nested Combo.
For example, in we know that flatness is the second element of phi:
phi
phi[1]
So, we know that flatness in phi should be True:
flatness in phi
True
But, the flatness element of phi is a Combo containing data_misfit and smallness, so data_misfit in phi should be False:
data_misfit in phi
False
And the same goes for smallness:
smallness in phi
False
Alternatively, we can use the contains() method to check recursively for a particular term:
phi.contains(data_misfit)
True
phi.contains(smallness)
True
phi.contains(flatness)
True
Linear combination of functions#
Since it’s possible to scale objective functions and also to add them up, we can define any linear combination of objective functions. For example, we could build an classic objective function like:
as follows:
beta = 1e2
alpha = 3.0
phi = data_misfit + beta * (smallness + alpha * flatness)
phi
We can evaluate the phi function and get its gradient and Hessian as with any other objective function:
phi(model)
np.float64(518977.97897255165)
phi.gradient(model)
array([ 144012.76749299, -29065.37179219, 63841.4074464 ,
-117750.76338047, -34073.33393289])
phi.hessian(model)
<5x5 _SumLinearOperator with dtype=float64>
Internal structure#
The phi object is a Combo with a nested structure. It contains two elements: the data_misfit and a Scaled object that represents the \(\beta \left[ phi_\text{0}(\mathbf{m}) + \alpha \phi_\text{1}(\mathbf{m}) \right]\):
len(phi)
2
phi[0]
phi[1]
Because the phi[1] is a Scaled, it has a multiplier (the value of \(\beta\)) and a function:
phi[1].multiplier
100.0
phi[1].function
The phi[1].function is another Combo that contains the smallness and a scaled version of flatness.
phi[1].function[0]
phi[1].function[1]
We can use the contains() to perform a recursive check for a given objective function within the phi function. This recursive check also explores
Scaled objects:
phi.contains(data_misfit)
True
phi.contains(smallness)
True
phi.contains(flatness)
True