Fork me on GitHub

geoscience / scientific computing / visualization / web design / thoughts

simpeg / SimPEG Tutorial - Linear Problem

by Rowan Cockett

Objective:

In this tutorial we will create a simple linear problem from scratch using the SimPEG framework.

The linear operator that we will be concerned with in this case will be a matrix with rows that are made up of the function:

$$\text{row}_i = \exp(-p i x)\cos(2 \pi q i x)$$

Where $p$ and $q$ are constants (in this case they are both set to $\frac{1}{4}$) and $x \in [0,1]$.

In [1]:
from SimPEG import *
Warning: mumps solver not available.

In [2]:
M = Mesh.TensorMesh([1000]) # Make a 1D mesh from 0-1 with 1000 cells

nk = 20       #: The number of kernals
p = q = 0.25  #: Some constants

# Create a function for each row of the matrix
g = lambda k, x: np.exp(-p*k*x)*np.cos(2*np.pi*q*k*x)

# Concatenate each row, and reshape
G = np.concatenate([g(i,M.vectorCCx) for i in range(nk)])
G.shape = (nk, M.nC)

figure(figsize=(12,4))
plot(G.T)
title('Rows of G')
Out[2]:
<matplotlib.text.Text at 0x10c764710>

Now that we have a matrix, we can create a few models that live at the cell-centers of our mesh. We will multiply these models by the matrix that we just created to create some predicted data:

$$ \mathbf{d}_{\text{pred}} = \mathbf{Pu} = \mathbf{PGm}$$

Here we will create data by projecting the field, $\mathbf{u}$, onto the data-space using a linear projection $\mathbf{P}$. We will create the field by multiplying a model, $\mathbf{m}$, by the matrix we just created $\mathbf{G}$. In this case, we will just assume that our linear projection is the identity (we sample the data everywhere!). This reduces the equation to:

$$ \mathbf{d}_{\text{pred}} = \mathbf{Gm}$$

An interesting question is to see if different models give different data (very important if we ever want to predict the model!).

Note: in numpy, sometimes you get weird results when you just multiply, so if you intend to do a matrix multiply, use np.dot()

In [3]:
mtrue = np.zeros(M.nC)  # Create a zeros vector, here nC means: number of cells
mtrue[M.vectorCCx > 0.3] = 1.
mtrue[M.vectorCCx > 0.45] = -0.5
mtrue[M.vectorCCx > 0.6] = 0
figure(figsize=(12,4))
subplot(121)
# Create a few models
models = np.c_[mtrue, np.r_[mtrue[200:],mtrue[:200]]*0.5 - 0.2,
               M.vectorCCx - 0.5, -M.vectorCCx + 0.75,
               Utils.ModelBuilder.randomModel(1000,seed=751,its=20000)-0.3]
plot(M.vectorCCx, models)
title('A few synthetic model.')
xlabel('x');ylabel('Amplitude')
data = G.dot(models)  #: this is matrix multiplication!!
subplot(122)
plot(data)
title('The data created by each model.')
xlim([0,19]);xlabel('Kernal Number, i');ylabel('g(i) * x')
Out[3]:
<matplotlib.text.Text at 0x10db27ad0>

What is important to see here is that each model gives us different data, which means that if we just have the data, we might be able to recover the model!

To do this efficiently (i.e. optimize) we will need the derivatives of this problem, luckily, the matrix that creates our data does not depend on our model (i.e. the inverse problem is linear). Our problem is as follows:

$$ \mathbf{d}_{\text{pred}} = \mathbf{Gm}$$

Where $\mathbf{m} \in \mathcal{R}^n$ is our model vector, $\mathbf{G} \in \mathcal{R}^{k,n}$ is our matrix of $k$ kernal functions that samples our model to produce our predicted data ($\mathbf{d}_\text{pred}$).

The question that we want to know is: "How sensitive is our data to a slight change in our model?" To answer this we shall use calculus! *gasp*

$$ \mathbf{J} = \frac{\partial \mathbf{d}_{\text{pred}}}{\partial \mathbf{m}} = \frac{\partial \mathbf{Gm} }{\partial \mathbf{m}} = \mathbf{G}\frac{\partial \mathbf{m} }{\partial \mathbf{m}} = \mathbf{G}$$

The derivative describes how the predicted data changes with small changes in our model, and in this case it is fairly straight forward.

All that SimPEG needs to know when you create a problem is:

  1. how to create data ($\mathbf{Gm}$),
  2. the sensitivities effect on a vector, $\mathbf{v}$, ($\mathbf{Jv = Gv}$), and
  3. the transpose sensitivities effect on a vector, $\mathbf{v}$, ($\mathbf{J^\top v = G^\top v}$).

To learn more about the problem class in SimPEG and how you can use it, check out the documentation at: http://simpeg.rtfd.org

In [4]:
class LinearSurvey(Survey.BaseSurvey):
    def projectFields(self, u):
        return u

class LinearProblem(Problem.BaseProblem):
    
    surveyPair = LinearSurvey
    
    def __init__(self, model, G, **kwargs):
        Problem.BaseProblem.__init__(self, model, **kwargs)
        self.G = G
        
    def fields(self, m):
        return self.G.dot(m)
        
    def Jvec(self, m, v, u=None):
        return self.G.dot(v)

    def Jtvec(self, m, v, u=None):
        return self.G.T.dot(v)
    

Once we have our problem, we can use the inversion tools in SimPEG to run our inversion:

In [5]:
# Create a model, by default this is a model on the mesh
model = Model.BaseModel(M)
# Create the problem
prob = LinearProblem(model, G)
# create a simpeg data object, using a synthetic model and 5% noise
survey = prob.createSyntheticSurvey(models[:,4],std=0.05)
# Create an optimization program
opt = Optimization.InexactGaussNewton(maxIter=20)
# Create a regularization program
reg = Regularization.Tikhonov(model)
# Create an objective function
beta = Parameters.BetaSchedule(beta0=1e-1)
obj = ObjFunction.BaseObjFunction(survey, reg, beta=beta)
# Create an inversion object
inv = Inversion.BaseInversion(obj, opt)
# Start the inversion with a model of zeros, and run the inversion
m0 = survey.mtrue*0
mopt = inv.run(m0)
Regularization has not set mref. SimPEG will set it to m0.
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  1.00e-01  3.94e+03  0.00e+00  3.94e+03    3.80e+03      0              
Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.
   1  1.00e-01  1.14e+03  3.65e-01  1.14e+03    5.76e+02      0              
   2  1.00e-01  6.38e+02  1.68e+00  6.38e+02    6.92e+02      0              
   3  5.00e-02  3.74e+02  4.55e+00  3.74e+02    6.00e+02      0   Skip BFGS  
   4  5.00e-02  2.65e+02  4.18e+00  2.65e+02    7.62e+02      0              
   5  5.00e-02  2.09e+02  5.70e+00  2.10e+02    4.60e+02      0              
   6  2.50e-02  1.37e+02  8.22e+00  1.37e+02    8.01e+02      0              
   7  2.50e-02  1.21e+02  9.29e+00  1.21e+02    5.88e+02      0              
   8  2.50e-02  1.17e+02  9.44e+00  1.17e+02    6.65e+02      0              
   9  1.25e-02  1.14e+02  9.81e+00  1.14e+02    6.18e+02      0              
  10  1.25e-02  1.10e+02  1.06e+01  1.10e+02    6.70e+02      0              
  11  1.25e-02  1.09e+02  1.05e+01  1.09e+02    6.29e+02      0              
  12  6.25e-03  1.08e+02  1.05e+01  1.08e+02    6.80e+02      0              
  13  6.25e-03  1.07e+02  1.02e+01  1.07e+02    6.86e+02      0   Skip BFGS  
  14  6.25e-03  1.06e+02  1.04e+01  1.06e+02    6.65e+02      0              
  15  3.13e-03  1.99e+01  1.19e+01  1.99e+01    1.46e+02      0   Skip BFGS  
  16  3.13e-03  1.56e+01  1.12e+01  1.56e+01    1.14e+02      0              
  17  3.13e-03  1.50e+01  1.13e+01  1.50e+01    1.42e+02      0              
  18  1.56e-03  1.22e+01  1.06e+01  1.22e+01    1.03e+02      0   Skip BFGS  
  19  1.56e-03  1.16e+01  1.09e+01  1.16e+01    8.47e+01      0              
  20  1.56e-03  1.15e+01  1.10e+01  1.15e+01    9.08e+01      0   Skip BFGS  
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.4209e-01 <= tolF*(1+|f0|) = 3.9393e+02
0 : |xc-x_last| = 1.0186e-01 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 9.0818e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 9.0818e+01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------

Explore the documentation to see what other parameters you can tweak in the different elements of the inversion, but let's check how well we recovered the model just by using the default parameters:

In [6]:
plot(M.vectorCCx, np.c_[survey.mtrue, mopt])
xlabel('x')
legend(('true model','recovered model'))
Out[6]:
<matplotlib.legend.Legend at 0x10db80290>

Hopefully you now have an idea of how to create a Problem class in SimPEG, and how this can be used with the other tools available.