Regularization - concepts explained#

In geophysical inversion, we minimize the data objective functional as the L2 norm of the misfit between data $d$ and the forward response $f$ of the model $m$, weighted by the data error $epsilon$:

$\Phi_d = \sum\limits_i^N \left(\frac{d_i-f_i(m)}{\epsilon_i}\right)^2=\|W_d(d-f(m))\|^2$

As this minimization problem is non-unique and ill-posed, we introduce a regularization term $Phi$, weighted by a regularization parameter $lambda$:

$\Phi = \Phi_d + \lambda \Phi_m$

The regularization strength $lambda$ should be chosen so that the data are fitted within noise, i.e. $chi^2=Phi_d/N=1$.

In the term $Phi-m$ we put our expectations to the model, e.g. to be close to any prior model. In many cases we do not have much information and aim for the smoothest model that is able to fit our data. We decribe it by the operator $W_m$:

$\Phi_m=\|W_m (m-m_{ref})\|^2$

The regularization operator is defined by some constraint operator $C$ weighted by some weighting function $w$ so that $W_m=mbox{diag}(w) C$. The operator $C$ can be a discrete smoothness operator, or the identity to keep the model close to the reference model $m_{ref}$.




import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.math.matrix import GeostatisticConstraintsMatrix
from pygimli.core.math import symlog


Regularization drives the model where the data are too weak to constrain the model. In order to explain different kinds of regularization (also called constraints), we use a very simple mapping forward operator: The values at certain positions are picked.

from pygimli.frameworks import PriorModelling


Implementation 1. determine the indices where the cells are

ind = [mesh.findCell(po).id() for po in pos]

1. forward response: take the model at indices

response = model[ind]

1. Jacobian matrix

J = pg.SparseMapMatrix()
J.resize(len(ind), mesh.cellCount())
for i, n in enumerate(self.ind):
self.J.setVal(i, n, 1.0)


We exemplify this on behalf of a simple triangular mesh in a rectangular domain.

rect = mt.createRectangle(start=[0, -10], end=[10, 0])
mesh = mt.createMesh(rect, quality=34.5, area=0.3)
print(mesh)

Mesh: Nodes: 377 Cells: 684 Boundaries: 1060


We define two positions where we associate two arbitrary values.

pos = [[3, -3], [7, -7]]
vals = np.array([20., 15.])
fop = PriorModelling(mesh, pos)


We set up an inversion instance with the forward operator and prepare the keywords for running the inversion always the same way: - the data vector - the error vector (as relative error) - a starting model value (could also be vector)

inv = pg.Inversion(fop=fop, verbose=False)
invkw = dict(dataVals=vals, errorVals=np.ones_like(vals)*0.03, startModel=10)


Classical smoothness constraints#

inv.setRegularization(cType=1)  # the default
result = inv.run(**invkw)
pg.show(mesh, result);

(<matplotlib.axes._subplots.AxesSubplot object at 0x7fe8f5515af0>, <matplotlib.colorbar.Colorbar object at 0x7fe8997a9a90>)


We will have a closer look at the regularization matrix $C$.

C = fop.constraints()
print(C.rows(), C.cols(), mesh)
ax, _ = pg.show(fop.constraints(), markersize=1)

row = C.row(111)
nz = np.nonzero(row)[0]
print(nz, row[nz])

992 684 Mesh: Nodes: 377 Cells: 684 Boundaries: 1060
[ 48 116] 2 [1.0, -1.0]


How does that change the regularization matrix $C$?

inv.setRegularization(cType=1, zWeight=0.2)  # the default
result = inv.run(**invkw)
pg.show(mesh, result)

RM = fop.regionManager()
cw = RM.constraintWeights()
print(min(cw), max(cw))

0.19999999999999996 1.0


Now we try some other regularization options.

inv.setRegularization(cType=0)  # damping difference to starting model
result = inv.run(**invkw)
ax, _ = pg.show(mesh, result)


Obviously, the damping keeps the model small ($log 1=0$) as the starting model is NOT a reference model by default. We will enable this by specifying the isReference switch.

invkw["isReference"] = True
result = inv.run(**invkw)
ax, cb = pg.show(mesh, result)


cType=10 means a mix between 1st order smoothness (1) and damping (0)

inv.setRegularization(cType=10)  # mix of 1st order smoothing and damping
result = inv.run(**invkw)
ax, _ = pg.show(mesh, result)


In the matrix both contributions are under each other

C = fop.constraints()
print(C.rows(), C.cols())
print(mesh)
ax, _ = pg.show(fop.constraints(), markersize=1)

1676 684
Mesh: Nodes: 377 Cells: 684 Boundaries: 1060


We see that we have the first order smoothness and the identity matrix below each other. We can also use a second-order (-1 2 -1) smoothness operator by cType=2.

inv.setRegularization(cType=2)  # 2nd order smoothing
result = inv.run(**invkw)
ax, _ = pg.show(mesh, result)


We have a closer look at the constraints matrix

C = fop.constraints()
print(C.rows(), C.cols(), mesh)
ax, _ = pg.show(C, markersize=1)

684 684 Mesh: Nodes: 377 Cells: 684 Boundaries: 1060


It looks like a Laplace operator and seems to have a wider range compared to first-order smoothness.

Geostatistical regularization#

The idea is that not only neighbors are correlated to each other but to have a wider correlation by using an operator

$\textbf{C}_{\text{M},ij}=\sigma^{2}\exp{\left( -\sqrt{ \left(\frac{\textbf{H}^x_{ij}}{I_{x}}\right)^{2}+ \left(\frac{\textbf{H}^y_{ij}}{I_{y}}\right)^{2} }\right)}.$

More details can be found in https://www.pygimli.org/_tutorials_auto/3_inversion/plot_6-geostatConstraints.html

We generate such a matrix and multiply it with a zero vector of just one 1. For displaying the wide range of magnitudes we use the symlog function

C = GeostatisticConstraintsMatrix(mesh=mesh, I=[8, 4], dip=-20)
print(C)

vec = pg.Vector(mesh.cellCount())
vec[mesh.findCell([5, -5]).id()] = 1.0
ax, _ = pg.show(mesh, symlog(C*vec, 1e-2), cMin=-2, cMax=2, cMap="bwr")

<pygimli.math.matrix.GeostatisticConstraintsMatrix object at 0x7fe8da959540>


For comparison, we use a much finer mesh and compute the same matrix

fineMesh = mt.createMesh(rect, area=0.03)
Cfine = GeostatisticConstraintsMatrix(mesh=fineMesh, I=[8, 4], dip=-20)
vec = pg.Vector(fineMesh.cellCount())
vec[fineMesh.findCell([5, -5]).id()] = 1.0
ax, _ = pg.show(fineMesh, symlog(Cfine*vec, 1e-2), cMin=-1, cMax=1, cMap="bwr")


Application#

We can pass the correlation length directly to the inversion instance

inv.setRegularization(correlationLengths=[2, 2, 2])
result = inv.run(**invkw)
ax, cb = pg.show(mesh, result)


This look structurally similar to the second-order smoothness, but can drive values outside the expected range in regions of no data coverage. We change the correlation lengths and the dip to be inclining

inv.setRegularization(correlationLengths=[2, 0.5, 2], dip=-20)
result = inv.run(**invkw)
ax, cb = pg.show(mesh, result)


We now add many more points.

N = 30
x = np.random.rand(N) * 10
y = -np.random.rand(N) * 10
v = np.random.rand(N) * 10 + 10
plt.plot(x, y, "*")

[<matplotlib.lines.Line2D object at 0x7fe89a0de0a0>]


and repeat the above computations

fop = PriorModelling(mesh, zip(x, y))
inv = pg.Inversion(fop=fop, verbose=True)
inv.setRegularization(correlationLengths=[4, 4])
result = inv.run(v, np.ones_like(v)*0.03, startModel=10)
ax, cb = pg.show(mesh, result)

fop: <pygimli.frameworks.modelling.PriorModelling object at 0x7fe8fa2bae50>
Data transformation: <pygimli.core._pygimli_.RTrans object at 0x7fe8da954520>
Model transformation (cumulative):
0 <pygimli.core._pygimli_.RTransLogLU object at 0x7fe8da954520>
min/max (data): 10.28/19.93
min/max (error): 3%/3%
min/max (start model): 10/10
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 6.16 (dPhi = 17.15%) lam: 20
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 6.15 (dPhi = 0.01%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 6.15 (dPhi = 0.0%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 0.0 (< 2.0%)                  #
################################################################################


Comparing the data with the model response is always a good idea.

plt.plot(v, inv.response, "*")

[<matplotlib.lines.Line2D object at 0x7fe899f1f9a0>]


Individual regularization operators#

Say you want to combine geostatistic operators with a damping, you can create a block matrix pasting the matric vertically.

C = pg.matrix.BlockMatrix()
G = pg.matrix.GeostatisticConstraintsMatrix(mesh=mesh, I=[2, 0.5], dip=-20)
I = pg.matrix.IdentityMatrix(mesh.cellCount(), val=0.1)
ax, _ = pg.show(C)


Note that in pg.matrix you find a lot of matrices and matrix generators.

We set this matrix directly and do the inversion.

fop.setConstraints(C)
result = inv.run(v, np.ones_like(v)*0.03, startModel=10, isReference=True)
ax, cb = pg.show(mesh, result)

fop: <pygimli.frameworks.modelling.PriorModelling object at 0x7fe8fa2bae50>
Data transformation: <pygimli.core._pygimli_.RTrans object at 0x7fe8ab142be0>
Model transformation (cumulative):
0 <pygimli.core._pygimli_.RTransLogLU object at 0x7fe8997ae940>
min/max (data): 10.28/19.93
min/max (error): 3%/3%
min/max (start model): 10/10
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 5.95 (dPhi = 23.43%) lam: 20
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 5.95 (dPhi = 0.02%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 5.95 (dPhi = 0.0%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 0.0 (< 2.0%)                  #
################################################################################


If you are using a method manager, you access the inversion instance by mgr.inv and the forward operator by mgr.fop.

Note

Take-away messages

• regularization drives the model where data are weak

• think and play with your assumptions to the model

• there are several predefined options

• geostatistical regularization can be superior, because: - it is mesh-independent - it better fills the data gaps (e.g. 3D inversion of 2D profiles)

Total running time of the script: ( 0 minutes 31.201 seconds)

Gallery generated by Sphinx-Gallery