Regularization

If there is one model that has a misfit that equals the desired tolerance, then there are infinitely many other models which can fit to the same degree. The challenge is to find that model which has the desired characteristics and is compatible with a priori information. A single model can be selected from an infinite ensemble by measuring the length, or norm, of each model. Then a smallest, or sometimes largest, member can be isolated. Our goal is to design a norm that embodies our prior knowledge and, when minimized, yields a realistic candidate for the solution of our problem. The norm can penalize variation from a reference model, spatial derivatives of the model, or some combination of these.

Tikhonov Regularization

Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:

\[R(m) = \int_\Omega \frac{\alpha_x}{2}\left(\frac{\partial m}{\partial x}\right)^2 + \frac{\alpha_y}{2}\left(\frac{\partial m}{\partial y}\right)^2 \partial v\]

Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:

\[\begin{split}R(m) \approx \sum_{ij} \left[\frac{\alpha_x}{2}\left[\left(\frac{m_{i+1,j} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i-1,j}}{h}\right)^2\right] \\ + \frac{\alpha_y}{2}\left[\left(\frac{m_{i,j+1} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i,j-1}}{h}\right)^2\right] \right]h^2\end{split}\]

If we let D_1 be the derivative matrix in the x direction

\[\mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1\]
\[\mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1\]

Where d_1 is the one dimensional derivative:

\[\begin{split}\mathbf{d}_1 = \frac{1}{h} \left[ \begin{array}{cccc} -1 & 1 & & \\ & \ddots & \ddots&\\ & & -1 & 1\end{array} \right]\end{split}\]
\[R(m) \approx \mathbf{v}^\top \left[\frac{\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \frac{\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \right]\]

Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)

\[\mathbf{a\odot b} = \text{diag}(\mathbf{a})\mathbf{b} = \text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}\]

and the transpose also is true (but the sizes have to make sense…):

\[\mathbf{a}^\top\text{diag}(\mathbf{b}) = \mathbf{b}^\top\text{diag}(\mathbf{a})\]

So R(m) can simplify to:

\[R(m) \approx \mathbf{m}^\top \left[\frac{\alpha_x}{2}\mathbf{D}_1^\top \text{diag}(\mathbf{A}_1^\top\mathbf{v}) \mathbf{D}_1 + \frac{\alpha_y}{2}\mathbf{D}_2^\top \text{diag}(\mathbf{A}_2^\top \mathbf{v}) \mathbf{D}_2 \right] \mathbf{m}\]

We will define W_x as:

\[\mathbf{W}_x = \sqrt{\alpha_x}\text{diag}\left(\sqrt{\mathbf{A}_1^\top\mathbf{v}}\right) \mathbf{D}_1\]

And then W as a tall matrix of all of the different regularization terms:

\[\begin{split}\mathbf{W} = \left[ \begin{array}{c} \mathbf{W}_s\\ \mathbf{W}_x\\ \mathbf{W}_y\end{array} \right]\end{split}\]

Then we can write

\[R(m) \approx \frac{1}{2}\mathbf{m^\top W^\top W m}\]

The API

class SimPEG.Regularization.RegularizationMesh(mesh, **kwargs)

Regularization Mesh

This contains the operators used in the regularization. Note that these are not necessarily true differential operators, but are constructed from a SimPEG Mesh.

param discretize.base.BaseMesh mesh

problem mesh

param numpy.ndarray indActive

bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells

Required Properties:

  • indActive (Array): active indices in mesh, a list or numpy array of <class ‘bool’>, <class ‘int’> with shape (*)

regularization_type = None
property indActive

indActive (Array): active indices in mesh, a list or numpy array of <class ‘bool’>, <class ‘int’> with shape (*)

property vol

reduced volume vector

Return type

numpy.ndarray

Returns

reduced cell volume

property nC

reduced number of cells

Return type

int

Returns

number of cells being regularized

property dim

dimension of regularization mesh (1D, 2D, 3D)

Return type

int

Returns

dimension

property Pac

projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)

Return type

scipy.sparse.csr_matrix

Returns

active cell projection matrix

property Pafx

projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )

Return type

scipy.sparse.csr_matrix

Returns

active face-x projection matrix

property Pafy

projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )

Return type

scipy.sparse.csr_matrix

Returns

active face-y projection matrix

property Pafz

projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )

Return type

scipy.sparse.csr_matrix

Returns

active face-z projection matrix

property aveFx2CC

averaging from active cell centers to active x-faces

Return type

scipy.sparse.csr_matrix

Returns

averaging from active cell centers to active x-faces

property aveCC2Fx

averaging from active x-faces to active cell centers

Return type

scipy.sparse.csr_matrix

Returns

averaging matrix from active x-faces to active cell centers

property aveFy2CC

averaging from active cell centers to active y-faces

Return type

scipy.sparse.csr_matrix

Returns

averaging from active cell centers to active y-faces

property aveCC2Fy

averaging from active y-faces to active cell centers

Return type

scipy.sparse.csr_matrix

Returns

averaging matrix from active y-faces to active cell centers

property aveFz2CC

averaging from active cell centers to active z-faces

Return type

scipy.sparse.csr_matrix

Returns

averaging from active cell centers to active z-faces

property aveCC2Fz

averaging from active z-faces to active cell centers

Return type

scipy.sparse.csr_matrix

Returns

averaging matrix from active z-faces to active cell centers

property cellDiffx

cell centered difference in the x-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the x-direction

property cellDiffy

cell centered difference in the y-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the y-direction

property cellDiffz

cell centered difference in the z-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the z-direction

property faceDiffx

x-face differences

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active faces in the x-direction

property faceDiffy

y-face differences

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active faces in the y-direction

property faceDiffz

z-face differences

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active faces in the z-direction

property cellDiffxStencil

cell centered difference stencil (no cell lengths include) in the x-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the x-direction

property cellDiffyStencil

cell centered difference stencil (no cell lengths include) in the y-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the y-direction

property cellDiffzStencil

cell centered difference stencil (no cell lengths include) in the y-direction

Return type

scipy.sparse.csr_matrix

Returns

differencing matrix for active cells in the y-direction