Repository: simpeg/simpeg
Pull request #943 (2020-11-05 22:55:13)

Cross-gradient joint inversion
This PR provides the code to implement cross-gradient / normalized cross-gradient joint inversion. 

A "coupling.py" was added in the regularization folder, where we can compute cross-gradient, normalized cross-gradient, Jacobian and Hessian of cross-gradient.

Also, add a few classes in the "directives.py". These classes are borrowed from SimPEG existing classes and modified to output the joint inversion results, say,  two different data misfit terms, two regularization terms and coupling terms. I think this part of the code can be highly optimized and combined with SimPEG original classes.

In addition, just a few minor changes, e.g., timing in the forward modeling, output iteration number of CG, etc.

Future works, a). I am going to do several tests for the "coupling.py". b). write the tutorial or documentation with more mathematical details about cross-gradient. c). when tests are done, I will write one or two examples.


Patch:
diff --git a/.gitignore b/.gitignore
index 9a9beab62..40a19f316 100644
--- a/.gitignore
+++ b/.gitignore
@@ -64,7 +64,8 @@ tutorials/06-ip/dcip2d/*
 tutorials/06-ip/dcip3d/*
 tutorials/06-ip/outputs/*
 tutorials/12-seismic/seismic/*
-tutorials/13-pgi/pgi/*
+tutorials/13-joint_inversion/pgi/*
+tutorials/13-joint_inversion/cross_gradient_data/*
 
 # jupyter
 *.ipynb
diff --git a/SimPEG/directives/__init__.py b/SimPEG/directives/__init__.py
index bc23d4e9b..aed37d2e1 100644
--- a/SimPEG/directives/__init__.py
+++ b/SimPEG/directives/__init__.py
@@ -24,3 +24,11 @@ from .pgi_directives import (
     PGI_BetaAlphaSchedule,
     PGI_AddMrefInSmooth,
 )
+
+from .sim_directives import (
+    SimilarityMeasureInversionDirective,
+    SimilarityMeasureSaveOutputEveryIteration,
+    PairedBetaEstimate_ByEig,
+    PairedBetaSchedule,
+    MovingAndMultiTargetStopping,
+)
diff --git a/SimPEG/directives/directives.py b/SimPEG/directives/directives.py
index 0561afb64..290193e8e 100644
--- a/SimPEG/directives/directives.py
+++ b/SimPEG/directives/directives.py
@@ -23,6 +23,7 @@ from ..regularization import (
     SimpleSmoothDeriv,
     SparseDeriv,
     PGIwithRelationships,
+    BaseSimilarityMeasure,
 )
 from ..utils import (
     mkvc,
@@ -35,6 +36,7 @@ from ..utils import (
     eigenvalue_by_power_iteration,
 )
 from ..utils.code_utils import deprecate_property
+from .. import optimization
 
 
 class InversionDirective(properties.HasProperties):
@@ -275,7 +277,7 @@ class BetaSchedule(InversionDirective):
 class AlphasSmoothEstimate_ByEig(InversionDirective):
     """
     Estimate the alphas multipliers for the smoothness terms of the regularization
-     as a multiple of the ratio between the highest eigenvalue of the
+    as a multiple of the ratio between the highest eigenvalue of the
     smallness term and the highest eigenvalue of each smoothness term of the regularization.
     The highest eigenvalue are estimated through power iterations and Rayleigh quotient.
     """
@@ -1637,14 +1639,16 @@ class UpdateSensitivityWeights(InversionDirective):
         # Normalize and threshold weights
         wr = np.zeros_like(self.invProb.model)
         for reg in self.reg.objfcts:
-            wr += reg.mapping.deriv(self.invProb.model).T * (
-                (reg.mapping * jtj_diag) / reg.objfcts[0].regmesh.vol ** 2.0
-            )
+            if not isinstance(reg, BaseSimilarityMeasure):
+                wr += reg.mapping.deriv(self.invProb.model).T * (
+                    (reg.mapping * jtj_diag) / reg.objfcts[0].regmesh.vol ** 2.0
+                )
         wr /= wr.max()
         wr += self.threshold
         wr **= 0.5
         for reg in self.reg.objfcts:
-            reg.cell_weights = reg.mapping * wr
+            if not isinstance(reg, BaseSimilarityMeasure):
+                reg.cell_weights = reg.mapping * wr
 
     def validate(self, directiveList):
         # check if a beta estimator is in the list after setting the weights
diff --git a/SimPEG/directives/sim_directives.py b/SimPEG/directives/sim_directives.py
new file mode 100644
index 000000000..3ff8a512e
--- /dev/null
+++ b/SimPEG/directives/sim_directives.py
@@ -0,0 +1,381 @@
+import numpy as np
+from ..regularization import BaseSimilarityMeasure
+from ..utils import eigenvalue_by_power_iteration
+from ..optimization import IterationPrinters, StoppingCriteria
+from .directives import InversionDirective, SaveEveryIteration
+
+
+###############################################################################
+#                                                                             #
+#              Directives of joint inversion                                  #
+#                                                                             #
+###############################################################################
+class SimilarityMeasureInversionPrinters:
+    betas = {
+        "title": "betas",
+        "value": lambda M: ["{:.2e}".format(elem) for elem in M.parent.betas],
+        "width": 26,
+        "format": "%s",
+    }
+    lambd = {
+        "title": "lambda",
+        "value": lambda M: M.parent.lambd,
+        "width": 10,
+        "format": "%1.2e",
+    }
+    phi_d_list = {
+        "title": "phi_d",
+        "value": lambda M: ["{:.2e}".format(elem) for elem in M.parent.phi_d_list],
+        "width": 26,
+        "format": "%s",
+    }
+    phi_m_list = {
+        "title": "phi_m",
+        "value": lambda M: ["{:.2e}".format(elem) for elem in M.parent.phi_m_list],
+        "width": 26,
+        "format": "%s",
+    }
+    phi_sim = {
+        "title": "phi_sim",
+        "value": lambda M: M.parent.phi_sim,
+        "width": 10,
+        "format": "%1.2e",
+    }
+    iterationCG = {
+        "title": "iterCG",
+        "value": lambda M: M.cg_count,
+        "width": 10,
+        "format": "%3d",
+    }
+
+
+class SimilarityMeasureInversionDirective(InversionDirective):
+    """
+    Directive for two model similiraty measure joint inversions. Sets Printers and
+    StoppingCriteria.
+
+    Notes
+    -----
+    Methods assume we are working with two models, and a single similarity measure.
+    Also, the SimilarityMeasure objective function must be the last regularization.
+    """
+
+    printers = [
+        IterationPrinters.iteration,
+        SimilarityMeasureInversionPrinters.betas,
+        SimilarityMeasureInversionPrinters.lambd,
+        IterationPrinters.f,
+        SimilarityMeasureInversionPrinters.phi_d_list,
+        SimilarityMeasureInversionPrinters.phi_m_list,
+        SimilarityMeasureInversionPrinters.phi_sim,
+        SimilarityMeasureInversionPrinters.iterationCG,
+    ]
+
+    def initialize(self):
+        if not isinstance(self.reg.objfcts[-1], BaseSimilarityMeasure):
+            raise TypeError(
+                f"The last regularization function must be an instance of "
+                f"BaseSimilarityMeasure, got {type(self.reg.objfcts[-1])}."
+            )
+
+        # define relevant attributes
+        self.betas = self.reg.multipliers[:-1]
+        self.lambd = self.reg.multipliers[-1]
+        self.phi_d_list = []
+        self.phi_m_list = []
+        self.phi_sim = 0.0
+
+        # pass attributes to invProb
+        self.invProb.betas = self.betas
+        self.invProb.num_models = len(self.betas)
+        self.invProb.lambd = self.lambd
+        self.invProb.phi_d_list = self.phi_d_list
+        self.invProb.phi_m_list = self.phi_m_list
+        self.invProb.phi_sim = self.phi_sim
+
+        self.opt.printers = self.printers
+        self.opt.stoppers = [StoppingCriteria.iteration]
+
+    def validate(self, directiveList):
+        # check that this directive is first in the DirectiveList
+        dList = directiveList.dList
+        self_ind = dList.index(self)
+        if self_ind != 0:
+            raise IndexError(
+                "The CrossGradientInversionDirective must be first in directive list."
+            )
+        return True
+
+    def endIter(self):
+        # compute attribute values
+        phi_d = []
+        for dmis in self.dmisfit.objfcts:
+            phi_d.append(dmis(self.opt.xc))
+
+        phi_m = []
+        for reg in self.reg.objfcts:
+            phi_m.append(reg(self.opt.xc))
+
+        # pass attributes values to invProb
+        self.invProb.phi_d_list = phi_d
+        self.invProb.phi_m_list = phi_m[:-1]
+        self.invProb.phi_sim = phi_m[-1]
+        self.invProb.betas = self.reg.multipliers[:-1]
+        # Assume last reg.objfct is the coupling
+        self.invProb.lambd = self.reg.multipliers[-1]
+
+
+class SimilarityMeasureSaveOutputEveryIteration(SaveEveryIteration):
+    """
+    SaveOutputEveryIteration for Joint Inversions.
+    Saves information on the tradeoff parameters, data misfits, regularizations,
+    coupling term, number of CG iterations, and value of cost function.
+    """
+
+    header = None
+    save_txt = True
+    betas = None
+    phi_d = None
+    phi_m = None
+    phi_sim = None
+    phi = None
+
+    def initialize(self):
+        if self.save_txt is True:
+            print(
+                "CrossGradientSaveOutputEveryIteration will save your inversion "
+                "progress as: '###-{0!s}.txt'".format(self.fileName)
+            )
+            f = open(self.fileName + ".txt", "w")
+            self.header = "  #          betas            lambda         joint_phi_d                joint_phi_m            phi_sim       iterCG     phi    \n"
+            f.write(self.header)
+            f.close()
+
+        # Create a list of each
+        self.betas = []
+        self.lambd = []
+        self.phi_d = []
+        self.phi_m = []
+        self.phi = []
+        self.phi_sim = []
+
+    def endIter(self):
+
+        self.betas.append(["{:.2e}".format(elem) for elem in self.invProb.betas])
+        self.phi_d.append(["{:.3e}".format(elem) for elem in self.invProb.phi_d_list])
+        self.phi_m.append(["{:.3e}".format(elem) for elem in self.invProb.phi_m_list])
+        self.lambd.append("{:.2e}".format(self.invProb.lambd))
+        self.phi_sim.append(self.invProb.phi_sim)
+        self.phi.append(self.opt.f)
+
+        if self.save_txt:
+            f = open(self.fileName + ".txt", "a")
+            i = self.opt.iter
+            f.write(
+                " {0:2d}  {1}  {2}  {3}  {4}  {5:1.4e}  {6:d}  {7:1.4e}\n".format(
+                    i,
+                    self.betas[i - 1],
+                    self.lambd[i - 1],
+                    self.phi_d[i - 1],
+                    self.phi_m[i - 1],
+                    self.phi_sim[i - 1],
+                    self.opt.cg_count,
+                    self.phi[i - 1],
+                )
+            )
+            f.close()
+
+    def load_results(self):
+        results = np.loadtxt(self.fileName + str(".txt"), comments="#")
+        self.betas = results[:, 1]
+        self.lambd = results[:, 2]
+        self.phi_d = results[:, 3]
+        self.phi_m = results[:, 4]
+        self.phi_sim = results[:, 5]
+        self.f = results[:, 7]
+
+
+class PairedBetaEstimate_ByEig(InversionDirective):
+    """
+    Estimate the trade-off parameter, beta, between pairs of data misfit(s) and the
+    regularization(s) as a multiple of the ratio between the highest eigenvalue of the
+    data misfit term and the highest eigenvalue of the regularization.
+    The highest eigenvalues are estimated through power iterations and Rayleigh
+    quotient.
+
+    Notes
+    -----
+    This class assumes the order of the data misfits for each model parameter match
+    the order for the respective regularizations, i.e.
+
+    >>> data_misfits = [phi_d_m1, phi_d_m2, phi_d_m3]
+    >>> regs = [phi_m_m1, phi_m_m2, phi_m_m3]
+
+    In which case it will estimate regularization parameters for each respective pair.
+    """
+
+    beta0_ratio = 1.0  #: the estimated ratio is multiplied by this to obtain beta
+    n_pw_iter = 4  #: number of power iterations for estimation.
+    seed = None  #: Random seed for the directive
+
+    def initialize(self):
+        """
+        The initial beta is calculated by comparing the estimated
+        eigenvalues of JtJ and WtW.
+        To estimate the eigenvector of **A**, we will use one iteration
+        of the *Power Method*:
+
+        .. math::
+            \\mathbf{x_1 = A x_0}
+        Given this (very course) approximation of the eigenvector, we can
+        use the *Rayleigh quotient* to approximate the largest eigenvalue.
+
+        .. math::
+            \\lambda_0 = \\frac{\\mathbf{x^\\top A x}}{\\mathbf{x^\\top x}}
+        We will approximate the largest eigenvalue for both JtJ and WtW,
+        and use some ratio of the quotient to estimate beta0.
+
+        .. math::
+            \\beta_0 = \\gamma \\frac{\\mathbf{x^\\top J^\\top J x}}{\\mathbf{x^\\top W^\\top W x}}
+        :rtype: float
+        :return: beta0
+        """
+        if self.seed is not None:
+            np.random.seed(self.seed)
+
+        if self.debug:
+            print("Calculating the beta0 parameter.")
+
+        m = self.invProb.model
+        dmis_eigenvalues = []
+        reg_eigenvalues = []
+        dmis_objs = self.dmisfit.objfcts
+        reg_objs = [
+            obj
+            for obj in self.reg.objfcts
+            if not isinstance(obj, BaseSimilarityMeasure)
+        ]
+        if len(dmis_objs) != len(reg_objs):
+            raise ValueError(
+                f"There must be the same number of data misfit and regularizations."
+                f"Got {len(dmis_objs)} and {len(reg_objs)} respectively."
+            )
+        for dmis, reg in zip(dmis_objs, reg_objs):
+            dmis_eigenvalues.append(
+                eigenvalue_by_power_iteration(dmis, m, n_pw_iter=self.n_pw_iter,)
+            )
+
+            reg_eigenvalues.append(
+                eigenvalue_by_power_iteration(reg, m, n_pw_iter=self.n_pw_iter,)
+            )
+
+        self.ratios = np.array(dmis_eigenvalues) / np.array(reg_eigenvalues)
+        self.invProb.betas = self.beta0_ratio * self.ratios
+        self.reg.multipliers[:-1] = self.invProb.betas
+
+
+class PairedBetaSchedule(InversionDirective):
+    """
+    Directive for beta cooling schedule to determine the tradeoff
+    parameters when using paired data misfits and regularizations for a joint inversion.
+    """
+
+    chifact_target = 1.0
+    beta_tol = 1e-1
+    update_beta = True
+    cooling_rate = 1
+    cooling_factor = 2
+    dmis_met = False
+
+    @property
+    def target(self):
+        if getattr(self, "_target", None) is None:
+            nD = np.array([survey.nD for survey in self.survey])
+
+            self._target = nD * 0.5 * self.chifact_target
+
+        return self._target
+
+    @target.setter
+    def target(self, val):
+        self._target = val
+
+    def initialize(self):
+        self.dmis_met = np.zeros_like(self.invProb.betas, dtype=bool)
+
+    def endIter(self):
+
+        # Check if target misfit has been reached, if so, set dmis_met to True
+        for i, phi_d in enumerate(self.invProb.phi_d_list):
+            self.dmis_met[i] = phi_d < self.target[i]
+
+        # check separately if misfits are within the tolerance,
+        # otherwise, scale beta individually
+        for i, phi_d in enumerate(self.invProb.phi_d_list):
+            if self.opt.iter > 0 and self.opt.iter % self.cooling_rate == 0:
+                target = self.target[i]
+                ratio = phi_d / target
+                if (
+                    self.update_beta
+                    and self.dmis_met[i]
+                    and np.abs(1.0 - ratio) > self.beta_tol
+                ):
+                    if ratio <= 1:
+                        ratio = np.minimum(1.5, ratio)
+                    else:
+                        ratio = np.maximum(0.75, ratio)
+
+                    self.invProb.betas[i] *= ratio
+                elif not self.dmis_met[i]:
+                    self.invProb.betas[i] /= self.cooling_factor
+
+        self.reg.multipliers[:-1] = self.invProb.betas
+
+
+class MovingAndMultiTargetStopping(InversionDirective):
+    r"""
+        Directive for setting stopping criteria for a joint inversion.
+        Ensures both that all target misfits are met and there is a small change in the
+        model. Computes the percentage change of the current model from the previous model.
+
+        ..math::
+            \frac {\| \mathbf{m_i} - \mathbf{m_{i-1}} \|} {\| \mathbf{m_{i-1}} \|}
+    """
+
+    tol = 1e-5
+    beta_tol = 1e-1
+    chifact_target = 1.0
+
+    @property
+    def target(self):
+        if getattr(self, "_target", None) is None:
+            nD = []
+            for survey in self.survey:
+                nD += [survey.nD]
+            nD = np.array(nD)
+
+            self._target = nD * 0.5 * self.chifact_target
+
+        return self._target
+
+    @target.setter
+    def target(self, val):
+        self._target = val
+
+    def endIter(self):
+        for phi_d, target in zip(self.invProb.phi_d_list, self.target):
+            if np.abs(1.0 - phi_d / target) >= self.beta_tol:
+                return
+        if (
+            np.linalg.norm(self.opt.xc - self.opt.x_last)
+            / np.linalg.norm(self.opt.x_last)
+            > self.tol
+        ):
+            return
+
+        print(
+            "stopping criteria met: ",
+            np.linalg.norm(self.opt.xc - self.opt.x_last)
+            / np.linalg.norm(self.opt.x_last),
+        )
+        self.opt.stopNextIteration = True
diff --git a/SimPEG/objective_function.py b/SimPEG/objective_function.py
index 7292f716f..b80b60b3e 100644
--- a/SimPEG/objective_function.py
+++ b/SimPEG/objective_function.py
@@ -140,11 +140,12 @@ class BaseObjectiveFunction(BaseSimPEG):
                 x = np.random.randn(self.nP)
 
         v = x + 0.1 * np.random.rand(len(x))
+        expectedOrder = kwargs.pop("expectedOrder", 1)
         return checkDerivative(
             lambda m: [self.deriv(m).dot(v), self.deriv2(m, v=v)],
             x,
             num=num,
-            expectedOrder=1,
+            expectedOrder=expectedOrder,
             plotIt=plotIt,
             **kwargs
         )
diff --git a/SimPEG/optimization.py b/SimPEG/optimization.py
index d313d7d5f..d8537791c 100644
--- a/SimPEG/optimization.py
+++ b/SimPEG/optimization.py
@@ -228,6 +228,13 @@ class IterationPrinters(object):
         "format": "%1.2e",
     }
 
+    iterationCG = {
+        "title": "iterCG",
+        "value": lambda M: M.cg_count,
+        "width": 10,
+        "format": "%3d",
+    }
+
 
 class Minimize(object):
     """
@@ -1237,7 +1244,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
         findSearchDirection()
         Finds the search direction based on projected CG
         """
-
+        self.cg_count = 0
         Active = self.activeSet(self.xc)
         temp = sum((np.ones_like(self.xc.size) - Active))
 
diff --git a/SimPEG/potential_fields/gravity/simulation.py b/SimPEG/potential_fields/gravity/simulation.py
index 17557dd0b..4f719d77f 100644
--- a/SimPEG/potential_fields/gravity/simulation.py
+++ b/SimPEG/potential_fields/gravity/simulation.py
@@ -81,7 +81,6 @@ class Simulation3DIntegral(BasePFSimulation):
         Gravity forward operator
         """
         if getattr(self, "_G", None) is None:
-
             self._G = self.linear_operator()
 
         return self._G
diff --git a/SimPEG/regularization/__init__.py b/SimPEG/regularization/__init__.py
index 950de7c68..1b39d6261 100644
--- a/SimPEG/regularization/__init__.py
+++ b/SimPEG/regularization/__init__.py
@@ -1,4 +1,4 @@
-from .base import BaseRegularization, BaseComboRegularization
+from .base import BaseRegularization, BaseComboRegularization, BaseSimilarityMeasure
 from .regularization_mesh import RegularizationMesh
 from .tikhonov import (
     SimpleSmall,
@@ -16,3 +16,5 @@ from .pgi import (
     PGIwithNonlinearRelationshipsSmallness,
     PGIwithRelationships,
 )
+from .cross_gradient import CrossGradient
+from .correspondence import LinearCorrespondence
diff --git a/SimPEG/regularization/base.py b/SimPEG/regularization/base.py
index 652a0f2ac..d60b6400c 100644
--- a/SimPEG/regularization/base.py
+++ b/SimPEG/regularization/base.py
@@ -141,14 +141,14 @@ class BaseRegularization(BaseObjectiveFunction):
 
         .. math::
 
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top
+            R(m) = \\frac{1}{2}\\mathbf{(m-m_\\text{ref})^\\top W^\\top
                    W(m-m_\\text{ref})}
 
         So the derivative is straight forward:
 
         .. math::
 
-            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+            R(m) = \\mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
 
@@ -383,3 +383,90 @@ class BaseComboRegularization(SimpleComboRegularization):
                 ), "cell_weights must be length {} not {}".format(
                     self._nC_residual, len(change["value"])
                 )
+
+
+###############################################################################
+#                                                                             #
+#                        Base Coupling Regularization                         #
+#                                                                             #
+###############################################################################
+class BaseSimilarityMeasure(BaseRegularization):
+
+    """
+    Base class for the similarity term in joint inversions. Inherit this for building
+    your own similarity term.  The BaseSimilarityMeasure assumes two different
+    geophysical models through one similarity term. However, if you wish
+    to combine more than two models, e.g., 3 models,
+    you may want to add a total of three coupling terms:
+
+    e.g., lambda1*(m1, m2) + lambda2*(m1, m3) + lambda3*(m2, m3)
+
+    where, lambdas are weights for coupling terms. m1, m2 and m3 indicate
+    three different models.
+    """
+
+    wire_map = properties.Instance(
+        "Wire Map for the two coupled parameters", maps.Wires,
+    )
+
+    def __init__(self, mesh, wire_map, **kwargs):
+        super().__init__(mesh, wire_map=wire_map, **kwargs)
+        # do this as a hack to make TreeMesh work.
+        self.regmesh.regularization_type = "Tikhonov"
+
+    @properties.validator("wire_map")
+    def _wire_map_validator(self, change):
+        map = change["value"]
+        try:
+            m1, m2 = map.maps  # Assume a map has been passed for each model.
+        except ValueError:
+            ValueError("Wire map must have two model mappings")
+
+        if m1[1].shape[0] != m2[1].shape[0]:
+            raise ValueError(
+                f"All models must be the same size! Got {m1[1].shape[0]} and {m2[1].shape[0]}"
+            )
+
+    @property
+    def nP(self):
+        """
+        number of model parameters
+        """
+        return self.wire_map.nP
+
+    def deriv(self):
+        """
+        First derivative of the coupling term with respect to individual models.
+        Returns an array of dimensions [k*M,1],
+        k: number of models we are inverting for.
+        M: number of cells in each model.
+
+        """
+        raise NotImplementedError(
+            "The method deriv has not been implemented for {}".format(
+                self.__class__.__name__
+            )
+        )
+
+    def deriv2(self):
+        """
+        Second derivative of the coupling term with respect to individual models.
+        Returns either an array of dimensions [k*M,1] (v is not None), or
+        sparse matrix of dimensions [k*M, k*M] (v is None).
+        k: number of models we are inverting for.
+        M: number of cells in each model.
+
+        """
+        raise NotImplementedError(
+            "The method _deriv2 has not been implemented for {}".format(
+                self.__class__.__name__
+            )
+        )
+
+    def __call__(self):
+        """ Returns the computed value of the coupling term. """
+        raise NotImplementedError(
+            "The method __call__ has not been implemented for {}".format(
+                self.__class__.__name__
+            )
+        )
diff --git a/SimPEG/regularization/correspondence.py b/SimPEG/regularization/correspondence.py
new file mode 100644
index 000000000..34d8364b2
--- /dev/null
+++ b/SimPEG/regularization/correspondence.py
@@ -0,0 +1,105 @@
+import numpy as np
+import scipy.sparse as sp
+import properties
+
+from .. import utils
+from .base import BaseSimilarityMeasure
+
+
+class LinearCorrespondence(BaseSimilarityMeasure):
+    """
+    The petrophysical linear constraint for joint inversions.
+
+    ..math::
+        \\phi_c({\\mathbf m}_{\\mathbf1},{\\mathbf m}_{\\mathbf2})=\\lambda\\sum_{i=1}^M
+        (k_1*m_1 + k_2*m_2 + k_3)
+
+    Assuming that we are working with two models only.
+
+    """
+
+    coefficients = properties.Array(
+        "coefficients for the linear relationship between parameters",
+        shape=(3,),
+        default=np.array([1.0, -1.0, 0.0]),
+    )
+
+    def relation(self, model):
+        """
+        Computes the values of petrophysical linear relationship between two different
+        geophysical models.
+
+        The linear relationship is defined as:
+
+        f(m1, m2)  = k1*m1 + k2*m2 + k3
+
+        :param numpy.ndarray model: stacked array of individual models
+                                    np.c_[model1, model2,...]
+
+        :rtype: float
+        :return: linearly related petrophysical values of two different models,
+                  dimension: M by 1, :M number of model parameters.
+
+        """
+        m1, m2 = self.wire_map * model
+        k1, k2, k3 = self.coefficients
+
+        return k1 * m1 + k2 * m2 + k3
+
+    def __call__(self, model):
+        """
+        Computes the sum of values of petrophysical linear relationship
+        between two different geophysical models.
+
+        :param numpy.ndarray model: stacked array of individual models
+                                    np.c_[model1, model2,...]
+
+        :rtype: float
+        :return: a scalar value.
+        """
+
+        result = self.relation(model)
+        return 0.5 * result.T @ result
+
+    def deriv(self, model):
+        """Computes the Jacobian of the coupling term.
+
+        :param list of numpy.ndarray ind_models: [model1, model2,...]
+
+        :rtype: numpy.ndarray
+        :return: result: gradient of the coupling term with respect to model1, model2,
+                 :dimension 2M by 1, :M number of model parameters.
+        """
+        k1, k2, k3 = self.coefficients
+        r = self.relation(model)
+        dc_dm1 = k1 * r
+        dc_dm2 = k2 * r
+
+        result = np.r_[dc_dm1, dc_dm2]
+
+        return result
+
+    def deriv2(self, model, v=None):
+        """Computes the Hessian of the linear coupling term.
+
+        :param list of numpy.ndarray ind_models: [model1, model2, ...]
+        :param numpy.ndarray v: vector to be multiplied by Hessian
+        :rtype: scipy.sparse.csr_matrix if v is None
+                numpy.ndarray if v is not None
+        :return Hessian matrix: | h11, h21 | :dimension 2M*2M.
+                                |          |
+                                | h12, h22 |
+        """
+
+        k1, k2, k3 = self.coefficients
+        if v is not None:
+            v1, v2 = self.wire_map * v
+            p1 = k1 ** 2 * v1 + k2 * k1 * v2
+            p2 = k2 * k1 * v1 + k2 ** 2 * v2
+            return np.r_[p1, p2]
+        else:
+            n = self.regmesh.nC
+            A = utils.sdiag(np.ones(n) * (k1 ** 2))
+            B = utils.sdiag(np.ones(n) * (k2 ** 2))
+            C = utils.sdiag(np.ones(n) * (k1 * k2))
+            return sp.bmat([[A, C], [C, B]], format="csr")
diff --git a/SimPEG/regularization/cross_gradient.py b/SimPEG/regularization/cross_gradient.py
new file mode 100644
index 000000000..0b4179afb
--- /dev/null
+++ b/SimPEG/regularization/cross_gradient.py
@@ -0,0 +1,256 @@
+import numpy as np
+import scipy.sparse as sp
+import properties
+
+from .base import BaseSimilarityMeasure
+
+
+###############################################################################
+#                                                                             #
+#                                Cross-Gradient                               #
+#                                                                             #
+###############################################################################
+
+
+class CrossGradient(BaseSimilarityMeasure):
+    """
+    The cross-gradient constraint for joint inversions.
+
+    ..math::
+        \\phi_c(\\mathbf{m_1},\\mathbf{m_2}) = \\lambda \\sum_{i=1}^{M} \\|
+        \\nabla \\mathbf{m_1}_i \\times \\nabla \\mathbf{m_2}_i \\|^2
+
+    All methods assume that we are working with two models only.
+
+    """
+
+    # reset this here to clear out the properties attribute
+    cell_weights = None
+
+    # These are not fully implemented yet
+    # grad_tol = properties.Float(
+    #     "tolerance for avoiding the exteremly small gradient amplitude", default=1e-10
+    # )
+    # normalized = properties.Bool(
+    #     "whether to implement normalized cross-gradient", default=False
+    # )
+
+    approx_hessian = properties.Bool(
+        "whether to use the semi-positive definate approximation for the hessian",
+        default=True,
+    )
+
+    def __init__(self, mesh, wire_map, **kwargs):
+        super().__init__(mesh, wire_map=wire_map, **kwargs)
+
+        regmesh = self.regmesh
+
+        if regmesh.mesh.dim not in (2, 3):
+            raise ValueError("Cross-Gradient is only defined for 2D or 3D")
+        self._G = regmesh.cell_gradient
+        self._Av = sp.diags(np.sqrt(regmesh.vol)) * regmesh.average_face_to_cell
+
+    def _calculate_gradient(self, model, normalized=False, rtol=1e-6):
+        """
+        Calculate the spatial gradients of the model using central difference.
+
+        Concatenates gradient components into a single array.
+        [[x_grad1, y_grad1, z_grad1],
+         [x_grad2, y_grad2, z_grad2],
+         [x_grad3, y_grad3, z_grad3],...]
+
+        :param numpy.ndarray model: model
+
+        :rtype: numpy.ndarray
+        :return: gradient_vector: array where each row represents a model cell,
+                 and each column represents a component of the gradient.
+
+        """
+        regmesh = self.regmesh
+        Avs = [regmesh.aveFx2CC, regmesh.aveFy2CC]
+        if regmesh.dim == 3:
+            Avs.append(regmesh.aveFz2CC)
+        Av = sp.block_diag(Avs)
+        gradient = (Av @ (self._G @ model)).reshape((-1, regmesh.dim), order="F")
+
+        if normalized:
+            norms = np.linalg.norm(gradient, axis=-1)
+            ind = norms <= norms.max() * rtol
+            norms[ind] = 1.0
+            gradient /= norms[:, None]
+            gradient[ind] = 0.0
+            # set gradient to 0 if amplitude of gradient is extremely small
+
+        return gradient
+
+    def calculate_cross_gradient(self, model, normalized=False, rtol=1e-6):
+        """
+        Calculates the cross-gradients of the models at each cell center.
+
+        Parameters
+        ----------
+        model : numpy.ndarray
+            The input model, which will be automatically separated into the two
+            parameters internally
+        normalized : bool, optional
+            Whether to normalize the gradient
+        rtol : float, optional
+            relative cuttoff for small gradients in the normalization
+
+        Returns
+        -------
+        cross_grad : numpy.ndarray
+            The norm of the cross gradient vector in each active cell.
+        """
+        m1, m2 = self.wire_map * model
+        # Compute the gradients and concatenate components.
+        grad_m1 = self._calculate_gradient(m1, normalized=normalized, rtol=rtol)
+        grad_m2 = self._calculate_gradient(m2, normalized=normalized, rtol=rtol)
+
+        # for each model cell, compute the cross product of the gradient vectors.
+        cross_prod = np.cross(grad_m1, grad_m2)
+        if self.regmesh.dim == 3:
+            cross_prod = np.linalg.norm(cross_prod, axis=-1)
+
+        return cross_prod
+
+    def __call__(self, model):
+        """
+        Computes the sum of all cross-gradient values at all cell centers.
+
+        :param numpy.ndarray model: stacked array of individual models
+                                    np.c_[model1, model2,...]
+        :param bool normalized: returns value of normalized cross-gradient if True
+
+        :rtype: float
+        :returns: the computed value of the cross-gradient term.
+
+
+        ..math::
+
+            \\phi_c(\\mathbf{m_1},\\mathbf{m_2})
+
+            = \\lambda \\sum_{i=1}^{M} \\|\\nabla \\mathbf{m_1}_i \\times \\nabla \\mathbf{m_2}_i \\|^2
+
+            = \\sum_{i=1}^{M} \\|\\nabla \\mathbf{m_1}_i\\|^2 \\ast \\|\\nabla \\mathbf{m_2}_i\\|^2
+                - (\\nabla \\mathbf{m_1}_i \\cdot \\nabla \\mathbf{m_2}_i )^2
+
+            = \\|\\phi_{cx}\\|^2 + \\|\\phi_{cy}\\|^2 + \\|\\phi_{cz}\\|^2 (optional strategy, not used in this script)
+
+
+        """
+        m1, m2 = self.wire_map * model
+        Av = self._Av
+        G = self._G
+        g_m1 = G @ m1
+        g_m2 = G @ m2
+        return 0.5 * np.sum(
+            (Av @ g_m1 ** 2) * (Av @ g_m2 ** 2) - (Av @ (g_m1 * g_m2)) ** 2
+        )
+
+    def deriv(self, model):
+        """
+        Computes the Jacobian of the cross-gradient.
+
+        :param list of numpy.ndarray ind_models: [model1, model2,...]
+
+        :rtype: numpy.ndarray
+        :return: result: gradient of the cross-gradient with respect to model1, model2
+
+        """
+        m1, m2 = self.wire_map * model
+
+        Av = self._Av
+        G = self._G
+        g_m1 = G @ m1
+        g_m2 = G @ m2
+
+        return np.r_[
+            (((Av @ g_m2 ** 2) @ Av) * g_m1) @ G
+            - (((Av @ (g_m1 * g_m2)) @ Av) * g_m2) @ G,
+            (((Av @ g_m1 ** 2) @ Av) * g_m2) @ G
+            - (((Av @ (g_m1 * g_m2)) @ Av) * g_m1) @ G,
+        ]
+
+    def deriv2(self, model, v=None):
+        """
+        Computes the Hessian of the cross-gradient.
+
+        :param list of numpy.ndarray ind_models: [model1, model2, ...]
+        :param numpy.ndarray v: vector to be multiplied by Hessian
+
+        :rtype: scipy.sparse.csr_matrix if v is None
+                numpy.ndarray if v is not None
+        :return Hessian matrix if v is None
+                Hessian multiplied by vector if v is not No
+
+        """
+        m1, m2 = self.wire_map * model
+
+        Av = self._Av
+        G = self._G
+
+        g_m1 = G @ m1
+        g_m2 = G @ m2
+
+        if v is None:
+            A = (
+                G.T
+                @ (
+                    sp.diags(Av.T @ (Av @ g_m2 ** 2))
+                    - sp.diags(g_m2) @ Av.T @ Av @ sp.diags(g_m2)
+                )
+                @ G
+            )
+
+            C = (
+                G.T
+                @ (
+                    sp.diags(Av.T @ (Av @ g_m1 ** 2))
+                    - sp.diags(g_m1) @ Av.T @ Av @ sp.diags(g_m1)
+                )
+                @ G
+            )
+
+            B = None
+            BT = None
+            if not self.approx_hessian:
+                # d_m1_d_m2
+                B = (
+                    G.T
+                    @ (
+                        2 * sp.diags(g_m1) @ Av.T @ Av @ sp.diags(g_m2)
+                        - sp.diags(g_m2) @ Av.T @ Av @ sp.diags(g_m1)
+                        - sp.diags(Av.T @ Av @ (g_m1 * g_m2))
+                    )
+                    @ G
+                )
+                BT = B.T
+
+            return sp.bmat([[A, B], [BT, C]], format="csr")
+        else:
+            v1, v2 = self.wire_map * v
+
+            Gv1 = G @ v1
+            Gv2 = G @ v2
+
+            p1 = G.T @ (
+                (Av.T @ (Av @ g_m2 ** 2)) * Gv1 - g_m2 * (Av.T @ (Av @ (g_m2 * Gv1)))
+            )
+            p2 = G.T @ (
+                (Av.T @ (Av @ g_m1 ** 2)) * Gv2 - g_m1 * (Av.T @ (Av @ (g_m1 * Gv2)))
+            )
+
+            if not self.approx_hessian:
+                p1 += G.T @ (
+                    2 * g_m1 * (Av.T @ (Av @ (g_m2 * Gv2)))
+                    - g_m2 * (Av.T @ (Av @ (g_m1 * Gv2)))
+                    - (Av.T @ (Av @ (g_m1 * g_m2))) * Gv2
+                )
+
+                p2 += G.T @ (
+                    2 * g_m2 * (Av.T @ (Av @ (g_m1 * Gv1)))
+                    - g_m1 * (Av.T @ (Av @ (g_m2 * Gv1)))
+                    - (Av.T @ (Av @ (g_m2 * g_m1))) * Gv1
+                )
+            return np.r_[p1, p2]
diff --git a/SimPEG/regularization/regularization_mesh.py b/SimPEG/regularization/regularization_mesh.py
index 0e876b974..68b74731b 100755
--- a/SimPEG/regularization/regularization_mesh.py
+++ b/SimPEG/regularization/regularization_mesh.py
@@ -19,7 +19,7 @@ class RegularizationMesh(props.BaseSimPEG):
 
     This contains the operators used in the regularization. Note that these
     are not necessarily true differential operators, but are constructed from
-    a SimPEG Mesh.
+    a `discretize` 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
@@ -175,6 +175,15 @@ class RegularizationMesh(props.BaseSimPEG):
                     self._Pafz = utils.speye(self.mesh.nFz)[:, indActive_Fz]
         return self._Pafz
 
+    @property
+    def average_face_to_cell(self):
+        if self.dim == 1:
+            return self.aveFx2CC
+        elif self.dim == 2:
+            return sp.hstack([self.aveFx2CC, self.aveFy2CC])
+        else:
+            return sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC])
+
     @property
     def aveFx2CC(self):
         """
@@ -288,6 +297,15 @@ class RegularizationMesh(props.BaseSimPEG):
                 )
         return self._aveCC2Fz
 
+    @property
+    def cell_gradient(self):
+        if self.dim == 1:
+            return self.cellDiffx
+        elif self.dim == 2:
+            return sp.vstack([self.cellDiffx, self.cellDiffy])
+        else:
+            return sp.vstack([self.cellDiffx, self.cellDiffy, self.cellDiffz])
+
     @property
     def cellDiffx(self):
         """
diff --git a/docs/content/release/0.15.0-notes.rst b/docs/content/release/0.15.0-notes.rst
index 21780504c..57cf94905 100644
--- a/docs/content/release/0.15.0-notes.rst
+++ b/docs/content/release/0.15.0-notes.rst
@@ -19,7 +19,7 @@ PGI
 ===
 We have merged @thast's petrophysically guided inversion (PGI) methodology into SimPEG!
 There any many new directives and objective functions related to this methodology of
-joint inversion. There are also now new :ref:`tutorials <sphx_glr_content_tutorials_13-pgi>` showing how you can implement this
+joint inversion. There are also now new :ref:`tutorials <sphx_glr_content_tutorials_13-joint_inversion>` showing how you can implement this
 in both joint inversions and single-domain inversions to incorporate petrophysical
 information into the inversion! Check out the multiple PR's related to this branch
 addition for more details.
diff --git a/tutorials/13-joint_inversion/README.txt b/tutorials/13-joint_inversion/README.txt
new file mode 100644
index 000000000..eebd9357d
--- /dev/null
+++ b/tutorials/13-joint_inversion/README.txt
@@ -0,0 +1,2 @@
+Joint Inversion
+===============
diff --git a/tutorials/13-pgi/plot_inv_1_joint_pf_pgi_full_info_tutorial.py b/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
similarity index 97%
rename from tutorials/13-pgi/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
rename to tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
index fd1eba942..f77e8051b 100644
--- a/tutorials/13-pgi/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
+++ b/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
@@ -230,18 +230,12 @@ magmap = actvMap * wires.sus
 idenMap = maps.IdentityMap(nP=nactv)
 # Grav problem
 simulation_grav = pf.gravity.simulation.Simulation3DIntegral(
-    survey=data_grav.survey,
-    mesh=mesh,
-    rhoMap=wires.den,
-    actInd=actv,
+    survey=data_grav.survey, mesh=mesh, rhoMap=wires.den, actInd=actv,
 )
 dmis_grav = data_misfit.L2DataMisfit(data=data_grav, simulation=simulation_grav)
 # Mag problem
 simulation_mag = pf.magnetics.simulation.Simulation3DIntegral(
-    survey=data_mag.survey,
-    mesh=mesh,
-    chiMap=wires.sus,
-    actInd=actv,
+    survey=data_mag.survey, mesh=mesh, chiMap=wires.sus, actInd=actv,
 )
 dmis_mag = data_misfit.L2DataMisfit(data=data_mag, simulation=simulation_mag)
 
@@ -355,20 +349,12 @@ Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=alpha0_ratio, verbos
 # initialize beta and beta/alpha_s schedule
 beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-2)
 betaIt = directives.PGI_BetaAlphaSchedule(
-    verbose=True,
-    coolingFactor=2.0,
-    tolerance=0.2,
-    progress=0.2,
+    verbose=True, coolingFactor=2.0, tolerance=0.2, progress=0.2,
 )
 # geophy. and petro. target misfits
-targets = directives.MultiTargetMisfits(
-    verbose=True,
-)
+targets = directives.MultiTargetMisfits(verbose=True,)
 # add learned mref in smooth once stable
-MrefInSmooth = directives.PGI_AddMrefInSmooth(
-    wait_till_stable=True,
-    verbose=True,
-)
+MrefInSmooth = directives.PGI_AddMrefInSmooth(wait_till_stable=True, verbose=True,)
 # update the parameters in smallness (L2-approx of PGI)
 update_smallness = directives.PGI_UpdateParameters(
     update_gmm=False  # keep GMM model fixed
diff --git a/tutorials/13-pgi/plot_inv_2_joint_pf_pgi_no_info_tutorial.py b/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
similarity index 97%
rename from tutorials/13-pgi/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
rename to tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
index b4cf7c293..6bc26d832 100644
--- a/tutorials/13-pgi/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
+++ b/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
@@ -231,18 +231,12 @@ magmap = actvMap * wires.sus
 idenMap = maps.IdentityMap(nP=nactv)
 # Grav problem
 simulation_grav = pf.gravity.simulation.Simulation3DIntegral(
-    survey=data_grav.survey,
-    mesh=mesh,
-    rhoMap=wires.den,
-    actInd=actv,
+    survey=data_grav.survey, mesh=mesh, rhoMap=wires.den, actInd=actv,
 )
 dmis_grav = data_misfit.L2DataMisfit(data=data_grav, simulation=simulation_grav)
 # Mag problem
 simulation_mag = pf.magnetics.simulation.Simulation3DIntegral(
-    survey=data_mag.survey,
-    mesh=mesh,
-    chiMap=wires.sus,
-    actInd=actv,
+    survey=data_mag.survey, mesh=mesh, chiMap=wires.sus, actInd=actv,
 )
 dmis_mag = data_misfit.L2DataMisfit(data=data_mag, simulation=simulation_mag)
 
@@ -370,21 +364,14 @@ Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=alpha0_ratio, verbos
 # initialize beta and beta/alpha_s schedule
 beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-4)
 betaIt = directives.PGI_BetaAlphaSchedule(
-    verbose=True,
-    coolingFactor=2.0,
-    tolerance=0.2,
-    progress=0.2,
+    verbose=True, coolingFactor=2.0, tolerance=0.2, progress=0.2,
 )
 # geophy. and petro. target misfits
 targets = directives.MultiTargetMisfits(
-    verbose=True,
-    chiSmall=0.5,  # ask for twice as much clustering (target value is /2)
+    verbose=True, chiSmall=0.5,  # ask for twice as much clustering (target value is /2)
 )
 # add learned mref in smooth once stable
-MrefInSmooth = directives.PGI_AddMrefInSmooth(
-    wait_till_stable=True,
-    verbose=True,
-)
+MrefInSmooth = directives.PGI_AddMrefInSmooth(wait_till_stable=True, verbose=True,)
 # update the parameters in smallness (L2-approx of PGI)
 update_smallness = directives.PGI_UpdateParameters(
     update_gmm=True,  # update the GMM each iteration
diff --git a/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py b/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py
new file mode 100755
index 000000000..b4238806f
--- /dev/null
+++ b/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py
@@ -0,0 +1,508 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Cross-gradient Joint Inversion of Gravity and Magnetic Anomaly Data
+===================================================================
+
+Here we simultaneously invert gravity and magentic data using cross-gradient
+constraint. The recovered density and susceptibility models are supposed to have
+structural similarity. For this tutorial, we focus on the following:
+
+    - Defining the survey from xyz formatted data
+    - Generating a mesh based on survey geometry
+    - Including surface topography
+    - Defining the inverse problem via combmaps (2 data misfit terms,
+        2 regularization terms, a coupling term and optimization)
+    - Specifying directives for the inversion
+    - Plotting the recovered model and data misfit
+
+
+Although we consider gravity and magnetic anomaly data in this tutorial,
+the same approach can be used to invert gradiometry and other types of geophysical data.
+
+"""
+
+#########################################################################
+# Import modules
+# --------------
+#
+
+import os
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import tarfile
+
+from discretize import TensorMesh
+
+from SimPEG.utils import plot2Ddata, surface2ind_topo
+from SimPEG.potential_fields import gravity, magnetics
+from SimPEG import (
+    maps,
+    data,
+    data_misfit,
+    inverse_problem,
+    regularization,
+    optimization,
+    directives,
+    inversion,
+    utils,
+)
+
+np.random.seed(0)
+
+#############################################
+# Define File Names
+# -----------------
+#
+# File paths for assets we are loading. To set up the inversion, we require
+# topography and field observations. The true model defined on the whole mesh
+# is loaded to compare with the inversion result. These files are stored as a
+# tar-file on our google cloud bucket:
+# "https://storage.googleapis.com/simpeg/doc-assets/gravity.tar.gz"
+
+# # storage bucket where we have the data
+data_source = (
+    "https://storage.googleapis.com/simpeg/doc-assets/cross_gradient_data.tar.gz"
+)
+
+# # download the data
+downloaded_data = utils.download(data_source, overwrite=True)
+
+# unzip the tarfile
+tar = tarfile.open(downloaded_data, "r")
+tar.extractall()
+tar.close()
+
+# path to the directory containing our data
+dir_path = downloaded_data.split(".")[0] + os.path.sep
+
+# files to work with
+topo_filename = dir_path + "topo.txt"
+model_filename = dir_path + "true_model.txt"
+
+
+#############################################
+# Load Data and Plot
+# ------------------
+#
+# Here we load and plot synthetic gravity anomaly data. Topography is generally
+# defined as an (N, 3) array. Gravity data is generally defined with 4 columns:
+# x, y, z and data.
+#
+
+# Load topography
+xyz_topo = np.loadtxt(topo_filename)
+
+# Load field data
+dobs_grav = np.loadtxt(dir_path + "gravity_data.obs")
+dobs_mag = np.loadtxt(dir_path + "magnetic_data.obs")
+
+# Define receiver locations and observed data
+receiver_locations = dobs_grav[:, 0:3]
+
+dobs_grav = dobs_grav[:, -1]
+dobs_mag = dobs_mag[:, -1]
+
+# Plot
+mpl.rcParams.update({"font.size": 12})
+
+# gravity data
+fig = plt.figure(figsize=(7, 5))
+
+ax1 = fig.add_axes([0.1, 0.1, 0.73, 0.85])
+plot2Ddata(receiver_locations, dobs_grav, ax=ax1, contourOpts={"cmap": "bwr"})
+ax1.set_title("Gravity Anomaly")
+ax1.set_xlabel("x (m)")
+ax1.set_ylabel("y (m)")
+
+ax2 = fig.add_axes([0.8, 0.1, 0.03, 0.85])
+norm = mpl.colors.Normalize(
+    vmin=-np.max(np.abs(dobs_grav)), vmax=np.max(np.abs(dobs_grav))
+)
+cbar = mpl.colorbar.ColorbarBase(
+    ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr, format="%.1e"
+)
+cbar.set_label("$mgal$", rotation=270, labelpad=15, size=12)
+
+# magnetic data
+fig = plt.figure(figsize=(7, 5))
+ax1 = fig.add_axes([0.1, 0.1, 0.73, 0.85])
+plot2Ddata(receiver_locations, dobs_mag, ax=ax1, contourOpts={"cmap": "bwr"})
+ax1.set_title("Magnetic Anomaly")
+ax1.set_xlabel("x (m)")
+ax1.set_ylabel("y (m)")
+
+ax2 = fig.add_axes([0.8, 0.1, 0.03, 0.85])
+norm = mpl.colors.Normalize(
+    vmin=-np.max(np.abs(dobs_mag)), vmax=np.max(np.abs(dobs_mag))
+)
+cbar = mpl.colorbar.ColorbarBase(
+    ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr, format="%.1e"
+)
+cbar.set_label("$nT$", rotation=270, labelpad=15, size=12)
+
+plt.show()
+
+
+#############################################
+# Assign Uncertainties
+# --------------------
+#
+# Inversion with SimPEG requires that we define standard deviation on our data.
+# This represents our estimate of the noise in our data. For gravity inversion,
+# a constant floor value is generally applied to all data. For this tutorial,
+# the standard deviation on each datum will be 1% of the maximum observed
+# gravity anomaly value. For magnetic inversion, the same strategy is performed.
+#
+
+maximum_anomaly_grav = np.max(np.abs(dobs_grav))
+uncertainties_grav = 0.01 * maximum_anomaly_grav * np.ones(np.shape(dobs_grav))
+
+maximum_anomaly_mag = np.max(np.abs(dobs_mag))
+uncertainties_mag = 0.01 * maximum_anomaly_mag * np.ones(np.shape(dobs_mag))
+
+
+#############################################
+# Defining the Survey
+# -------------------
+#
+# Here, we define survey that will be used for this tutorial. Gravity
+# surveys are simple to create. The user only needs an (N, 3) array to define
+# the xyz locations of the observation locations. From this, the user can
+# define the receivers and the source field.
+#
+
+
+# Define the receivers. The data consist of vertical gravity anomaly measurements.
+# The set of receivers must be defined as a list.
+receiver_grav = gravity.receivers.Point(receiver_locations, components="gz")
+
+# Define the source field and survey for gravity data
+source_field_grav = gravity.sources.SourceField(receiver_list=[receiver_grav])
+survey_grav = gravity.survey.Survey(source_field_grav)
+
+
+# Define the component(s) of the field we want to simulate as a list of strings.
+# Here we simulation total magnetic intensity data.
+components = ["tmi"]
+
+# Use the observation locations and components to define the receivers. To
+# simulate data, the receivers must be defined as a list.
+receiver_mag = magnetics.receivers.Point(receiver_locations, components=components)
+
+# Define the inducing field H0 = (intensity [nT], inclination [deg], declination [deg])
+inclination = 90
+declination = 0
+strength = 50000
+inducing_field = (strength, inclination, declination)
+
+# Define the source field and survey for gravity data
+source_field_mag = magnetics.sources.SourceField(
+    receiver_list=[receiver_mag], parameters=inducing_field
+)
+survey_mag = magnetics.survey.Survey(source_field_mag)
+
+
+#############################################
+# Defining the Data
+# -----------------
+#
+# Here is where we define the data that are inverted. The data are defined by
+# the survey, the observation values and the standard deviation.
+#
+
+data_object_grav = data.Data(
+    survey_grav, dobs=dobs_grav, standard_deviation=uncertainties_grav
+)
+data_object_mag = data.Data(
+    survey_mag, dobs=dobs_mag, standard_deviation=uncertainties_mag
+)
+
+
+#############################################
+# Defining a Tensor Mesh
+# ----------------------
+#
+# Here, we create the tensor mesh that will be used to invert gravity anomaly
+# data. If desired, we could define an OcTree mesh.
+#
+
+dh = 5.0
+hx = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
+hy = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
+hz = [(dh, 5, -1.3), (dh, 15)]
+mesh = TensorMesh([hx, hy, hz], "CCN")
+
+
+########################################################
+# Starting/Reference Model and Mapping on Tensor Mesh
+# ---------------------------------------------------
+#
+# Here, we create starting and/or reference models for the inversion as
+# well as the mapping from the model space to the active cells. Starting and
+# reference models can be a constant background value or contain a-priori
+# structures. Here, the backgrounds are 1e-6 g/cc and 1e-6 SI for density and
+# susceptibility models, respectively. Note that the background values could
+# be different for density and susceptibility models.
+#
+
+# Define density contrast values for each unit in g/cc.
+background_dens, background_susc = 1e-6, 1e-6
+
+# Find the indicies of the active cells in forward model (ones below surface)
+ind_active = surface2ind_topo(mesh, xyz_topo)
+
+# Define mapping from model to active cells
+nC = int(ind_active.sum())
+model_map = maps.IdentityMap(nP=nC)  # model consists of a value for each active cell
+
+# Create Wires Map that maps from stacked models to individual model components
+# m1 refers to density model, m2 refers to susceptibility
+wires = maps.Wires(("density", nC), ("susceptibility", nC))
+
+# Define and plot starting model
+starting_model = np.r_[background_dens * np.ones(nC), background_susc * np.ones(nC)]
+
+
+##############################################
+# Define the Physics
+# ------------------
+#
+# Here, we define the physics of the gravity and magnetic problems by using the simulation
+# class.
+#
+
+simulation_grav = gravity.simulation.Simulation3DIntegral(
+    survey=survey_grav, mesh=mesh, rhoMap=wires.density, actInd=ind_active
+)
+
+simulation_mag = magnetics.simulation.Simulation3DIntegral(
+    survey=survey_mag,
+    mesh=mesh,
+    model_type="scalar",
+    chiMap=wires.susceptibility,
+    actInd=ind_active,
+)
+
+
+#######################################################################
+# Define the Inverse Problem
+# --------------------------
+#
+# The inverse problem is defined by 4 things:
+#
+#     1) Data Misfit: a measure of how well our recovered model explains the field data
+#     2) Regularization: constraints placed on the recovered model and a priori information
+#     3) Coupling: a connection of two different physical property models
+#     4) 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_grav = data_misfit.L2DataMisfit(data=data_object_grav, simulation=simulation_grav)
+dmis_mag = data_misfit.L2DataMisfit(data=data_object_mag, simulation=simulation_mag)
+
+# Define the regularization (model objective function).
+reg_grav = regularization.Simple(mesh, indActive=ind_active, mapping=wires.density)
+reg_mag = regularization.Simple(
+    mesh, indActive=ind_active, mapping=wires.susceptibility
+)
+
+# Define the coupling term to connect two different physical property models
+lamda = 2e12  # weight for coupling term
+cross_grad = regularization.CrossGradient(mesh, wires, indActive=ind_active)
+
+# combo
+dmis = dmis_grav + dmis_mag
+reg = reg_grav + reg_mag + lamda * cross_grad
+
+# Define how the optimization problem is solved. Here we will use a projected
+# Gauss-Newton approach that employs the conjugate gradient solver.
+opt = optimization.ProjectedGNCG(
+    maxIter=10,
+    lower=-2.0,
+    upper=2.0,
+    maxIterLS=20,
+    maxIterCG=100,
+    tolCG=1e-3,
+    tolX=1e-3,
+)
+
+# 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.PairedBetaEstimate_ByEig(beta0_ratio=1e0)
+# starting_beta.n_pw_iter = 10
+
+# Defining the fractional decrease in beta and the number of Gauss-Newton solves
+# for each beta value.
+beta_schedule = directives.PairedBetaSchedule(cooling_factor=5, cooling_rate=1)
+
+# Options for outputting recovered models and predicted data for each beta.
+save_iteration = directives.SimilarityMeasureSaveOutputEveryIteration(save_txt=False)
+
+joint_inv_dir = directives.SimilarityMeasureInversionDirective()
+
+stopping = directives.MovingAndMultiTargetStopping(tol=1e-6)
+
+sensitivity_weights = directives.UpdateSensitivityWeights(everyIter=False)
+
+# Updating the preconditionner if it is model dependent.
+update_jacobi = directives.UpdatePreconditioner()
+
+
+# The directives are defined as a list.
+directives_list = [
+    joint_inv_dir,
+    sensitivity_weights,
+    stopping,
+    starting_beta,
+    beta_schedule,
+    save_iteration,
+    update_jacobi,
+]
+
+#####################################################################
+# 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)
+
+# Run inversion
+recovered_model = inv.run(starting_model)
+
+
+############################################################
+# Plotting True Model and Recovered Model
+# ---------------------------------------
+#
+
+# Load the true model (was defined on the whole mesh) and extract only the
+# values on active cells.
+
+true_model_dens = np.loadtxt(dir_path + "true_model_dens.txt")
+true_model_dens[~ind_active] = np.NaN
+
+true_model_susc = np.loadtxt(dir_path + "true_model_susc.txt")
+true_model_susc[~ind_active] = np.NaN
+
+# Plot True Model
+fig = plt.figure(figsize=(9, 8))
+ax1 = plt.subplot(211)
+
+(im,) = mesh.plot_slice(true_model_dens, normal="Y", ax=ax1, grid=True)
+ax1.set_title("True density model slice at y = 0 m")
+cbar = plt.colorbar(im, format="%.1e")
+cbar.set_label("g/cc", rotation=270, labelpad=15, size=12)
+
+ax2 = plt.subplot(212)
+(im,) = mesh.plot_slice(
+    true_model_susc, normal="Y", ax=ax2, grid=True, pcolor_opts={"cmap": "inferno"}
+)
+
+ax2.set_title("True susceptibility model slice at y = 0 m")
+cbar = plt.colorbar(im, format="%.1e")
+cbar.set_label("SI", rotation=270, labelpad=15, size=12)
+plt.tight_layout()
+plt.show()
+
+# Plot Recovered Model
+m_dens_joint, m_susc_joint = wires * recovered_model
+plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)
+
+fig = plt.figure(figsize=(9, 8))
+ax1 = plt.subplot(211)
+(im,) = mesh.plotSlice(
+    plotting_map * m_dens_joint, normal="Y", ax=ax1, clim=(-0.04, 0.03),
+)
+ax1.set_title("Density model slice at y = 0 m")
+cbar = plt.colorbar(im)
+cbar.set_label("g/cc", rotation=270, labelpad=15, size=12)
+
+ax2 = plt.subplot(212)
+(im,) = mesh.plotSlice(
+    plotting_map * m_susc_joint, normal="Y", ax=ax2, pcolor_opts={"cmap": "inferno"}
+)
+ax2.set_title("Susceptibility model slice at y = 0 m")
+cbar = plt.colorbar(im)
+cbar.set_label("SI", rotation=270, labelpad=15, size=12)
+
+plt.tight_layout()
+plt.show()
+
+
+############################################################
+# Comparing jointly and separatly recovered models
+# ---------------------------------------
+#
+
+# Normalized Cross Gradient of Jointly Recovered Susceptibility and Density Models
+ncg = cross_grad.calculate_cross_gradient(recovered_model, normalized=True)
+
+fig = plt.figure(figsize=(9, 4))
+ax = plt.subplot(111)
+(im,) = mesh.plot_slice(plotting_map * ncg, normal="Y", ax=ax, grid=True,)
+ax.set_title("Normalized cross gradient for joint inversion slice at y = 0 m")
+cbar = plt.colorbar(im, format="%.1e")
+cbar.set_label("|cross grad|", rotation=270, labelpad=15, size=12)
+plt.show()
+
+# Normalized Cross Gradient of Separately Recovered Susceptibility and Density Models
+m_dens_single = np.loadtxt(dir_path + "single_model_dens.txt")
+m_susc_single = np.loadtxt(dir_path + "single_model_susc.txt")
+m_separate = np.r_[m_dens_single[ind_active], m_susc_single[ind_active]]
+
+ncg_single = cross_grad.calculate_cross_gradient(m_separate, normalized=True)
+
+fig = plt.figure(figsize=(9, 4))
+ax = plt.subplot(111)
+(im,) = mesh.plot_slice(plotting_map * ncg_single, normal="Y", ax=ax, grid=True,)
+ax.set_title("Normalized cross gradient for separate inversion slice at y = 0 m")
+cbar = plt.colorbar(im, format="%.1e")
+cbar.set_label("|cross grad|", rotation=270, labelpad=15, size=12)
+
+plt.show()
+
+
+# Cross Plots Recovered Susceptibility and Density Models
+fig = plt.figure(figsize=(14, 5))
+ax0 = plt.subplot(121)
+ax0.scatter(
+    plotting_map * m_dens_joint, plotting_map * m_susc_joint, s=4, c="black", alpha=0.1
+)
+
+ax0.set_xlabel("Density", size=12)
+ax0.set_ylabel("Susceptibility", size=12)
+ax0.tick_params(labelsize=12)
+ax0.set_title("Joint inversion")
+
+ax1 = plt.subplot(122)
+ax1.scatter(m_dens_single, m_susc_single, s=4, c="black", alpha=0.1)
+
+ax1.set_xlabel("Density", size=12)
+ax1.set_ylabel("Susceptibility", size=12)
+ax1.tick_params(labelsize=12)
+ax1.set_title("Separate inversion")
+
+plt.show()
diff --git a/tutorials/13-pgi/README.txt b/tutorials/13-pgi/README.txt
deleted file mode 100644
index f4f4f13ff..000000000
--- a/tutorials/13-pgi/README.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-PGI: Petrophysically and Geologically Guided Geophysical Inversion
-==================================================================


Test patch:
diff --git a/tests/base/test_correspondance.py b/tests/base/test_correspondance.py
new file mode 100644
index 000000000..b183ef3b2
--- /dev/null
+++ b/tests/base/test_correspondance.py
@@ -0,0 +1,57 @@
+import unittest
+
+import numpy as np
+
+from discretize import TensorMesh
+from SimPEG import (
+    maps,
+    regularization,
+)
+
+np.random.seed(10)
+
+
+class LinearCorrespondenceTest(unittest.TestCase):
+    def setUp(self):
+
+        dh = 1.0
+        nx = 12
+        ny = 12
+
+        hx = [(dh, nx)]
+        hy = [(dh, ny)]
+        mesh = TensorMesh([hx, hy], "CN")
+
+        # reg
+        actv = np.ones(len(mesh), dtype=bool)
+
+        # maps
+        wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
+
+        corr = regularization.LinearCorrespondence(
+            mesh, wire_map=wires, indActive=actv,
+        )
+
+        self.mesh = mesh
+        self.corr = corr
+
+    def test_order_full_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of linear correspondance with approx_hessian=True
+
+        """
+        corr = self.corr
+        self.assertTrue(corr._test_deriv())
+        self.assertTrue(corr._test_deriv2(expectedOrder=2))
+
+    def test_deriv2_no_arg(self):
+        m = np.random.randn(2 * len(self.mesh))
+
+        corr = self.corr
+
+        v = np.random.rand(len(m))
+
+        W = corr.deriv2(m)
+        Wv = corr.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
diff --git a/tests/base/test_cross_gradient.py b/tests/base/test_cross_gradient.py
new file mode 100644
index 000000000..f9ed0aa37
--- /dev/null
+++ b/tests/base/test_cross_gradient.py
@@ -0,0 +1,303 @@
+import unittest
+
+import numpy as np
+
+from discretize import TensorMesh, TreeMesh
+from SimPEG import (
+    maps,
+    regularization,
+)
+
+np.random.seed(10)
+
+
+class CrossGradientTensor2D(unittest.TestCase):
+    def setUp(self):
+
+        dh = 1.0
+        nx = 12
+        ny = 12
+
+        hx = [(dh, nx)]
+        hy = [(dh, ny)]
+        mesh = TensorMesh([hx, hy], "CN")
+
+        # reg
+        actv = np.ones(len(mesh), dtype=bool)
+
+        # maps
+        wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
+
+        cros_grad = regularization.CrossGradient(mesh, wire_map=wires, indActive=actv,)
+
+        self.mesh = mesh
+        self.cross_grad = cros_grad
+
+    def test_order_approximate_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = True
+        self.assertTrue(cross_grad.test())
+
+    def test_order_full_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = False
+        self.assertTrue(cross_grad._test_deriv())
+        self.assertTrue(cross_grad._test_deriv2(expectedOrder=2))
+
+    def test_deriv2_no_arg(self):
+        m = np.random.randn(2 * len(self.mesh))
+
+        cross_grad = self.cross_grad
+
+        v = np.random.rand(len(m))
+
+        cross_grad.approx_hessian = False
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+        cross_grad.approx_hessian = True
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+    def test_cross_grad_calc(self):
+        mesh = self.mesh
+        # get index of center point
+        midx = int(mesh.nCx / 2)
+        midy = int(mesh.nCy / 2)
+
+        # create a model1
+        m1 = np.zeros((mesh.nCx, mesh.nCy))
+        m1[(midx - 3) : (midx + 3), (midy - 3) : (midy + 3)] = 1
+
+        # create a model2
+        m2 = np.zeros((mesh.nCx, mesh.nCy))
+        m2[(midx - 5) : (midx + 1), (midy - 5) : (midy + 1)] = 1
+
+        m1 = m1.reshape(-1, order="F")
+        m2 = m2.reshape(-1, order="F")
+
+        # stack the true models
+        m = np.r_[m1, m2]
+
+        cross_grad = self.cross_grad
+
+        v1 = 0.5 * np.sum(np.abs(cross_grad.calculate_cross_gradient(m)))
+        v2 = cross_grad(m)
+        self.assertEqual(v1, v2)
+
+
+class CrossGradientTensor3D(unittest.TestCase):
+    def setUp(self):
+
+        dh = 1.0
+        nx = 12
+        ny = 12
+        nz = 12
+
+        hx = [(dh, nx)]
+        hy = [(dh, ny)]
+        hz = [(dh, nz)]
+        mesh = TensorMesh([hx, hy, hz], "CNN")
+
+        # reg
+        actv = np.ones(len(mesh), dtype=bool)
+
+        # maps
+        wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
+
+        cros_grad = regularization.CrossGradient(mesh, wire_map=wires, indActive=actv,)
+
+        self.mesh = mesh
+        self.cross_grad = cros_grad
+
+    def test_order_approximate_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = True
+        self.assertTrue(cross_grad.test())
+
+    def test_order_full_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = False
+        self.assertTrue(cross_grad._test_deriv())
+        self.assertTrue(cross_grad._test_deriv2(expectedOrder=2))
+
+    def test_deriv2_no_arg(self):
+        m = np.random.randn(2 * len(self.mesh))
+
+        cross_grad = self.cross_grad
+
+        v = np.random.rand(len(m))
+
+        cross_grad.approx_hessian = False
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+        cross_grad.approx_hessian = True
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+    def test_cross_grad_calc(self):
+        m = np.random.randn(2 * len(self.mesh))
+        cross_grad = self.cross_grad
+
+        m1, m2 = cross_grad.wire_map * m
+        v1 = cross_grad._calculate_gradient(m1)
+        v2 = cross_grad._calculate_gradient(m1, normalized=True)
+        np.testing.assert_allclose(v1 / np.linalg.norm(v1, axis=-1)[:, None], v2)
+
+        # test that calling the calculate_cross_gradient function works in 3D
+        cross_grad.calculate_cross_gradient(m)
+
+
+class CrossGradientTree2D(unittest.TestCase):
+    def setUp(self):
+
+        dh = 1.0
+        nx = 16
+        ny = 16
+
+        hx = [(dh, nx)]
+        hy = [(dh, ny)]
+        mesh = TreeMesh([hx, hy], "CN")
+        mesh.insert_cells([5, 5], [4])
+
+        # reg
+        actv = np.ones(len(mesh), dtype=bool)
+
+        # maps
+        wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
+
+        cross_grad = regularization.CrossGradient(mesh, wire_map=wires, indActive=actv)
+
+        self.mesh = mesh
+        self.cross_grad = cross_grad
+
+    def test_order_approximate_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = True
+        self.assertTrue(cross_grad.test())
+
+    def test_order_full_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = False
+        self.assertTrue(cross_grad._test_deriv())
+        self.assertTrue(cross_grad._test_deriv2(expectedOrder=2))
+
+    def test_deriv2_no_arg(self):
+        m = np.random.randn(2 * len(self.mesh))
+
+        cross_grad = self.cross_grad
+
+        v = np.random.rand(len(m))
+
+        cross_grad.approx_hessian = False
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+        cross_grad.approx_hessian = True
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+
+class CrossGradientTree3D(unittest.TestCase):
+    def setUp(self):
+
+        dh = 1.0
+        nx = 16
+        ny = 16
+        nz = 16
+
+        hx = [(dh, nx)]
+        hy = [(dh, ny)]
+        hz = [(dh, nz)]
+        mesh = TreeMesh([hx, hy, hz], "CNN")
+        mesh.insert_cells([5, 5, 5], [4])
+
+        # reg
+        actv = np.ones(len(mesh), dtype=bool)
+
+        # maps
+        wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
+
+        cross_grad = regularization.CrossGradient(mesh, wire_map=wires, indActive=actv)
+
+        self.mesh = mesh
+        self.cross_grad = cross_grad
+
+    def test_order_approximate_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = True
+        self.assertTrue(cross_grad.test())
+
+    def test_order_full_hessian(self):
+        """
+
+        Test deriv and deriv2 matrix of cross-gradient with approx_hessian=True
+
+        """
+        cross_grad = self.cross_grad
+        cross_grad.approx_hessian = False
+        self.assertTrue(cross_grad._test_deriv())
+        self.assertTrue(cross_grad._test_deriv2(expectedOrder=2))
+
+    def test_deriv2_no_arg(self):
+        m = np.random.randn(2 * len(self.mesh))
+
+        cross_grad = self.cross_grad
+
+        v = np.random.rand(len(m))
+
+        cross_grad.approx_hessian = False
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+        cross_grad.approx_hessian = True
+        W = cross_grad.deriv2(m)
+        Wv = cross_grad.deriv2(m, v)
+        np.testing.assert_allclose(Wv, W @ v)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index e6225c406..9a1b892e6 100644
--- a/tests/base/test_regularization.py
+++ b/tests/base/test_regularization.py
@@ -22,12 +22,15 @@ np.random.seed(639)
 IGNORE_ME = [
     "BaseRegularization",
     "BaseComboRegularization",
+    "BaseSimilarityMeasure",
     "SimpleComboRegularization",
     "BaseSparse",
     "PGI",
     "PGIwithRelationships",
     "PGIwithNonlinearRelationshipsSmallness",
     "PGIsmallness",
+    "CrossGradient",
+    "LinearCorrespondence",
 ]
 
 
@@ -284,7 +287,8 @@ class RegularizationTests(unittest.TestCase):
         )
 
         objfct = objective_function.L2ObjectiveFunction(
-            W=utils.sdiag(np.sqrt(cell_weights * mesh.cell_volumes)), mapping=wires.sigma
+            W=utils.sdiag(np.sqrt(cell_weights * mesh.cell_volumes)),
+            mapping=wires.sigma,
         )
 
         self.assertTrue(reg(m) == objfct(m))
diff --git a/tests/examples/test_tutorials_3.py b/tests/examples/test_tutorials_3.py
index 8fdd70ea3..02c66318d 100644
--- a/tests/examples/test_tutorials_3.py
+++ b/tests/examples/test_tutorials_3.py
@@ -10,7 +10,14 @@ matplotlib.use("Agg")
 
 dirname, filename = os.path.split(os.path.abspath(__file__))
 example_dir = dirname.split(os.path.sep)[:-2] + ["tutorials"]
-dirs_to_test = ["08-tdem", "09-nsem", "10-vrm", "11-flow", "12-seismic", "13-pgi"]
+dirs_to_test = [
+    "08-tdem",
+    "09-nsem",
+    "10-vrm",
+    "11-flow",
+    "12-seismic",
+    "13-joint_inversion",
+]
 
 
 class ExampleTest(unittest.TestCase):
diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py
index c7b95c20e..15163f814 100644
--- a/tests/flow/test_Richards.py
+++ b/tests/flow/test_Richards.py
@@ -100,7 +100,7 @@ class BaseRichardsTest(unittest.TestCase):
         passed = checkDerivative(
             lambda m: [self.prob.dpred(m), lambda v: self.prob.Jvec(m, v)],
             self.mtrue,
-            num=4,
+            num=3,
             plotIt=False,
         )
         self.assertTrue(passed, True)