The fatiando package has been deprecated. Please check out the new tools in the Fatiando a Terra website: www.fatiando.org

Data misfit functions (fatiando.inversion.misfit)

Defines base classes to represent a data-misfit functions (l2-norm, etc)

These classes can be used to implement parameter estimation problems (inversions). They automate most of the boiler plate required and provide direct access to ready-made optimization routines and regularization.

For now, only implements an l2-norm data misfit:

  • Misfit: an l2-norm data-misfit function

See the documentation for fatiando.inversion for examples of using Misfit.


class fatiando.inversion.misfit.Misfit(data, nparams, islinear, cache=True)[source]

Bases: fatiando.inversion.base.OptimizerMixin, fatiando.inversion.base.OperatorMixin

An l2-norm data-misfit function.

This is a kind of objective function that measures the misfit between observed data \(\bar{d}^o\) and data predicted by a set of model parameters \(\bar{d} = \bar{f}(\bar{p})\).

The l2-norm data-misfit is defined as:

\[\phi (\bar{p}) = \bar{r}^T \bar{r}\]

where \(\bar{r} = \bar{d}^o - \bar{d}\) is the residual vector and \(N\) is the number of data.

When subclassing this class, you must implement the method:

  • predicted(self, p): calculates the predicted data \(\bar{d}\) for a given parameter vector p

If you want to use any gradient-based solver (you probably do), you’ll need to implement the method:

  • jacobian(self, p): calculates the Jacobian matrix of \(\bar{f}(\bar{p})\) evaluated at p

If \(\bar{f}\) is linear, then the Jacobian will be cached in memory so that it is only calculated once when using the class multiple times. So solving the same problem with different methods or using an iterative method doesn’t have the penalty of recalculating the Jacobian.

Warning

When subclassing, be careful not to set the following attributes: data, nparams, islinear, nparams, ndata, and (most importantly) regul_param and _regularizing_parameter. This could mess with internal behavior and break things in unexpected ways.

Parameters:

  • data
    : 1d-array

    The observed data vector \(\bar{d}^o\)

  • nparams
    : int

    The number of parameters in parameter vector \(\bar{p}\)

  • islinear
    : True or False

    Whether \(\bar{f}\) is linear or not.

  • cache
    : True

    Whether or not to cache the output of some methods to avoid recomputing matrices and vectors when passed the same input parameter vector.

config(method, **kwargs)

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

  • method
    : string

    The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

  • newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)
  • acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

copy(deep=False)[source]

Make a copy of me together with all the cached methods.

estimate_

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a fatiando.mesher object.

fit()

Solve for the parameter vector that minimizes this objective function.

Uses the optimization method and parameters defined using the config method.

The estimated parameter vector can be accessed through the p_ attribute. A (possibly) formatted version (converted to a more manageable type) of the estimate can be accessed through the property estimate_.

fmt_estimate(p)

Called when accessing the property estimate_.

Use this to convert the parameter vector (p) to a more useful form, like a geometric object, etc.

Parameters:

  • p
    : 1d-array

    The parameter vector.

Returns:

  • formatted

    Pretty much anything you want.

gradient(p)[source]

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

  • p
    : 1d-array

    The parameter vector where the gradient is evaluated

Returns:

  • gradient
    : 1d-array

    The gradient vector.

hessian(p)[source]

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

  • p
    : 1d-array

    The parameter vector where the Hessian is evaluated

Returns:

  • hessian
    : 2d-array

    The Hessian matrix

predicted(p=None)[source]

Calculate the predicted data for a given parameter vector.

Parameters:

  • p
    : 1d-array or None

    The parameter vector used to calculate the predicted data. If None, will use the current estimate stored in estimate_.

Returns:

  • predicted
    : 1d-array or list of 1d-arrays

    The predicted data. If this is the sum of 1 or more Misfit instances, will return the predicted data from each of the summed misfits in the order of the sum.

regul_param

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

residuals(p=None)[source]

Calculate the residuals vector (observed - predicted data).

Parameters:

  • p
    : 1d-array or None

    The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

  • residuals
    : 1d-array or list of 1d-arrays

    The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)[source]

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

  • weights
    : 1d-array or 2d-array or None

    Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)[source]

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

  • p
    : 1d-array or None

    The parameter vector.

Returns:

  • value
    : float

    The value of the misfit function.