Linear Least-Squares Inversion¶

Here we demonstrate the basics of inverting data with SimPEG by considering a linear inverse problem. We formulate the inverse problem as a least-squares optimization problem. For this tutorial, we focus on the following:

Defining the forward problem

Defining the inverse problem (data misfit, regularization, optimization)

Specifying directives for the inversion

Recovering a set of model parameters which explains the observations

Import Modules¶

import numpy as np
import matplotlib.pyplot as plt

from discretize import TensorMesh

from SimPEG import (
    simulation,
    maps,
    data_misfit,
    directives,
    optimization,
    regularization,
    inverse_problem,
    inversion,
)

# sphinx_gallery_thumbnail_number = 3

Defining the Model and Mapping¶

Here we generate a synthetic model and a mappig which goes from the model space to the row space of our linear operator.

nParam = 100  # Number of model paramters

# A 1D mesh is used to define the row-space of the linear operator.
mesh = TensorMesh([nParam])

# Creating the true model
true_model = np.zeros(mesh.nC)
true_model[mesh.vectorCCx > 0.3] = 1.0
true_model[mesh.vectorCCx > 0.45] = -0.5
true_model[mesh.vectorCCx > 0.6] = 0

# Mapping from the model space to the row space of the linear operator
model_map = maps.IdentityMap(mesh)

# Plotting the true model
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
ax.plot(mesh.vectorCCx, true_model, "b-")
ax.set_ylim([-2, 2])

Out:

(-2, 2)

Defining the Linear Operator¶

Here we define the linear operator with dimensions (nData, nParam). In practive, you may have a problem-specific linear operator which you would like to construct or load here.

# Number of data observations (rows)
nData = 20

# Create the linear operator for the tutorial. The columns of the linear operator
# represents a set of decaying and oscillating functions.
jk = np.linspace(1.0, 60.0, nData)
p = -0.25
q = 0.25


def g(k):
    return np.exp(p * jk[k] * mesh.vectorCCx) * np.cos(
        np.pi * q * jk[k] * mesh.vectorCCx
    )


G = np.empty((nData, nParam))

for i in range(nData):
    G[i, :] = g(i)

# Plot the columns of G
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
for i in range(G.shape[0]):
    ax.plot(G[i, :])

ax.set_title("Columns of matrix G")

Out:

Text(0.5, 1.0, 'Columns of matrix G')

Defining the Simulation¶

The simulation defines the relationship between the model parameters and predicted data.

sim = simulation.LinearSimulation(mesh, G=G, model_map=model_map)

Out:

/Users/josephcapriotti/opt/anaconda3/envs/py36/lib/python3.6/site-packages/SimPEG/simulation.py:531: UserWarning: G has not been implemented for the simulation
  warnings.warn("G has not been implemented for the simulation")

Predict Synthetic Data¶

Here, we use the true model to create synthetic data which we will subsequently invert.

# Standard deviation of Gaussian noise being added
std = 0.01
np.random.seed(1)

# Create a SimPEG data object
data_obj = sim.make_synthetic_data(true_model, relative_error=std, add_noise=True)

Define the Inverse Problem¶

The inverse problem is defined by 3 things:

Data Misfit: a measure of how well our recovered model explains the field data

Regularization: constraints placed on the recovered model and a priori information

Optimization: the numerical approach used to solve the inverse problem

# Define the data misfit. Here the data misfit is the L2 norm of the weighted
# residual between the observed data and the data predicted for a given model.
# Within the data misfit, the residual between predicted and observed data are
# normalized by the data's standard deviation.
dmis = data_misfit.L2DataMisfit(simulation=sim, data=data_obj)

# Define the regularization (model objective function).
reg = regularization.Tikhonov(mesh, alpha_s=1.0, alpha_x=1.0)

# Define how the optimization problem is solved.
opt = optimization.InexactGaussNewton(maxIter=50)

# Here we define the inverse problem that is to be solved
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)

Define Inversion Directives¶

Here we define any directiveas that are carried out during the inversion. This includes the cooling schedule for the trade-off parameter (beta), stopping criteria for the inversion and saving inversion results at each iteration.

# Defining a starting value for the trade-off parameter (beta) between the data
# misfit and the regularization.
starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-2)

# Setting a stopping criteria for the inversion.
target_misfit = directives.TargetMisfit()

# The directives are defined as a list.
directives_list = [starting_beta, target_misfit]

Setting a Starting Model and Running the Inversion¶

To define the inversion object, we need to define the inversion problem and the set of directives. We can then run the inversion.

# Here we combine the inverse problem and the set of directives
inv = inversion.BaseInversion(inv_prob, directives_list)

# Starting model
starting_model = np.zeros(nParam)

# Run inversion
recovered_model = inv.run(starting_model)

Out:

SimPEG.InvProblem will set Regularization.mref to m0.

        SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
        ***Done using same Solver and solverOpts as the problem***
model has any nan: 0
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
6.04e+02  1.00e+05  0.00e+00  1.00e+05    1.26e+06      0
6.04e+02  4.68e+04  3.48e-01  4.70e+04    7.96e+04      0
6.04e+02  3.19e+04  1.30e+00  3.27e+04    5.73e+04      0
6.04e+02  1.77e+04  4.62e+00  2.05e+04    5.14e+04      0   Skip BFGS
6.04e+02  1.27e+04  4.60e+00  1.55e+04    7.73e+04      0
6.04e+02  8.97e+03  6.89e+00  1.31e+04    4.36e+04      0
6.04e+02  5.35e+03  9.47e+00  1.11e+04    5.29e+04      0
6.04e+02  4.67e+03  9.62e+00  1.05e+04    2.35e+04      0
6.04e+02  4.17e+03  1.01e+01  1.03e+04    6.71e+04      0
6.04e+02  3.72e+03  1.05e+01  1.01e+04    4.59e+04      0
6.04e+02  3.29e+03  1.10e+01  9.94e+03    6.23e+04      0
6.04e+02  2.77e+03  1.16e+01  9.80e+03    4.95e+04      0
6.04e+02  2.82e+03  1.15e+01  9.76e+03    4.80e+04      0
6.04e+02  2.66e+03  1.17e+01  9.74e+03    4.28e+04      0
6.04e+02  2.60e+03  1.17e+01  9.64e+03    1.70e+04      0   Skip BFGS
6.04e+02  2.55e+03  1.17e+01  9.63e+03    4.77e+03      0
6.04e+02  2.49e+03  1.18e+01  9.62e+03    1.52e+04      0   Skip BFGS
6.04e+02  2.51e+03  1.18e+01  9.62e+03    2.04e+03      0
6.04e+02  2.64e+03  1.15e+01  9.61e+03    3.80e+03      0   Skip BFGS
6.04e+02  2.54e+03  1.17e+01  9.60e+03    1.07e+04      0
6.04e+02  2.48e+03  1.18e+01  9.60e+03    1.34e+04      0   Skip BFGS
6.04e+02  2.53e+03  1.17e+01  9.60e+03    9.47e+03      0
6.04e+02  2.51e+03  1.17e+01  9.58e+03    1.02e+04      0   Skip BFGS
6.04e+02  2.55e+03  1.16e+01  9.58e+03    7.48e+03      0
6.04e+02  2.55e+03  1.16e+01  9.58e+03    6.37e+03      0   Skip BFGS
6.04e+02  2.55e+03  1.16e+01  9.58e+03    1.04e+04      0
6.04e+02  2.54e+03  1.17e+01  9.58e+03    6.78e+03      0   Skip BFGS
6.04e+02  2.54e+03  1.16e+01  9.58e+03    6.88e+03      0
6.04e+02  2.53e+03  1.17e+01  9.57e+03    6.89e+03      0
6.04e+02  2.51e+03  1.17e+01  9.56e+03    5.25e+03      0   Skip BFGS
6.04e+02  2.51e+03  1.17e+01  9.56e+03    5.39e+03      0
6.04e+02  2.51e+03  1.17e+01  9.56e+03    3.31e+03      0
6.04e+02  2.52e+03  1.17e+01  9.56e+03    1.64e+03      0
6.04e+02  2.52e+03  1.17e+01  9.56e+03    2.35e+03      0
6.04e+02  2.52e+03  1.17e+01  9.56e+03    6.80e+02      0   Skip BFGS
6.04e+02  2.52e+03  1.17e+01  9.56e+03    2.17e+03      0
6.04e+02  2.52e+03  1.17e+01  9.56e+03    2.08e+03      0   Skip BFGS
6.04e+02  2.52e+03  1.17e+01  9.56e+03    2.19e+03      0
6.04e+02  2.53e+03  1.16e+01  9.56e+03    3.74e+03      0   Skip BFGS
6.04e+02  2.53e+03  1.16e+01  9.56e+03    2.91e+03      0
6.04e+02  2.54e+03  1.16e+01  9.56e+03    2.78e+03      0   Skip BFGS
6.04e+02  2.53e+03  1.16e+01  9.56e+03    2.96e+03      0
6.04e+02  2.53e+03  1.16e+01  9.56e+03    2.25e+03      0   Skip BFGS
6.04e+02  2.53e+03  1.17e+01  9.56e+03    3.03e+03      0
6.04e+02  2.53e+03  1.16e+01  9.56e+03    2.09e+03      0
6.04e+02  2.53e+03  1.16e+01  9.56e+03    1.55e+03      0   Skip BFGS
6.04e+02  2.53e+03  1.16e+01  9.56e+03    2.07e+03      0
6.04e+02  2.54e+03  1.16e+01  9.56e+03    2.64e+03      0   Skip BFGS
6.04e+02  2.54e+03  1.16e+01  9.56e+03    2.25e+03      0
6.04e+02  2.54e+03  1.16e+01  9.56e+03    2.17e+03      0   Skip BFGS
6.04e+02  2.54e+03  1.16e+01  9.56e+03    2.20e+03      0
------------------------- STOP! -------------------------
: |fc-fOld| = 2.8706e-04 <= tolF*(1+|f0|) = 1.0000e+04
: |xc-x_last| = 3.8037e-04 <= tolX*(1+|x0|) = 1.0000e-01
: |proj(x-g)-x|    = 2.2012e+03 <= tolG          = 1.0000e-01
: |proj(x-g)-x|    = 2.2012e+03 <= 1e3*eps       = 1.0000e-02
: maxIter   =      50    <= iter          =     50
------------------------- DONE! -------------------------

Plotting Results¶

# Observed versus predicted data
fig, ax = plt.subplots(1, 2, figsize=(12 * 1.2, 4 * 1.2))
ax[0].plot(data_obj.dobs, "b-")
ax[0].plot(inv_prob.dpred, "r-")
ax[0].legend(("Observed Data", "Predicted Data"))

# True versus recovered model
ax[1].plot(mesh.vectorCCx, true_model, "b-")
ax[1].plot(mesh.vectorCCx, recovered_model, "r-")
ax[1].legend(("True Model", "Recovered Model"))
ax[1].set_ylim([-2, 2])