# Stochastic static fault slip inversion from geodetic data with non-negativity and bound constraints

J-M Nocquet

► To cite this version:

J-M Nocquet. Stochastic static fault slip inversion from geodetic data with non-negativity and bound constraints. Geophysical Journal International, 2018, 214 (1), pp.366 - 385. 10.1093/gji/ggy146. hal-01856117

HAL Id: hal-01856117

https://hal.science/hal-01856117

Submitted on 7 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stochastic static fault slip inversion from geodetic data with non-negativity and bound constraints

J.-M. Nocquet<fnref n="1" />,<fnref n="2" />

Accepted 2018 April 11. Received 2018 April 5; in original form 2017 June 22

## SUMMARY

Despite surface displacements observed by geodesy are linear combinations of slip at faults in an elastic medium, determining the spatial distribution of fault slip remains a ill-posed inverse problem. A widely used approach to circumvent the illness of the inversion is to add regularization constraints in terms of smoothing and/or damping so that the linear system becomes invertible. However, the choice of regularization parameters is often arbitrary, and sometimes leads to significantly different results. Furthermore, the resolution analysis is usually empirical and cannot be made independently of the regularization. The stochastic approach of inverse problems provides a rigorous framework where the a priori information about the searched parameters is combined with the observations in order to derive posterior probabilities of the unknown parameters. Here, I investigate an approach where the prior probability density function (pdf) is a multivariate Gaussian function, with single truncation to impose positivity of slip or double truncation to impose positivity and upper bounds on slip for interseismic modelling. I show that the joint posterior pdf is similar to the linear untruncated Gaussian case and can be expressed as a truncated multivariate normal (TMVN) distribution. The TMVN form can then be used to obtain semi-analytical formulae for the single, 2-D or n-D marginal pdf. The semi-analytical formula involves the product of a Gaussian by an integral term that can be evaluated using recent developments in TMVN probabilities calculations. Posterior mean and covariance can also be efficiently derived. I show that the maximum posterior (MAP) can be obtained using a non-negative least-squares algorithm for the single truncated case or using the bounded-variable least-squares algorithm for the double truncated case. I show that the case of independent uniform priors can be approximated using TMVN. The numerical equivalence to Bayesian inversions using Monte Carlo Markov chain (MCMC) sampling is shown for a synthetic example and a real case for interseismic modelling in Central Peru. The TMVN method overcomes several limitations of the Bayesian approach using MCMC sampling. First, the need of computer power is largely reduced. Second, unlike Bayesian MCMC-based approach, marginal pdf, mean, variance or covariance are obtained independently one from each other. Third, the probability and cumulative density functions can be obtained with any density of points. Finally, determining the MAP is extremely fast.

**Key words:** Inverse theory; Earthquake source observation; Satellite Geodesy.

## 1 INTRODUCTION

Surface deformation measured by geodetic data is the most direct observation of slip along faults at depth. Once given a discretized fault geometry, displacements, tilts and strain are all quantities linearly related to slip components at the subfault in an elastic medium. Despite its linearity, determining the finite fault slip distribution at depth from surface measurement remains an ill-posed problem. The illness of the inverse problem comes from a combination of the partial spatial sampling of the deformation field and a rapidly decreasing sensitivity of displacement, strain and tilt to slip when the distance from the sources increases.

A widely used approach to solve the inverse problem is to add regularization constraints, making the linear system invertible so that a single solution is obtained. Regularization constraints are also usually tuned so that they exclude too oscillatory unphysical solutions and avoid overfitting the data. Nonetheless, several difficulties arise from this approach. First, it is uncertain whether the obtained solutions provide the complete range of models allowed by the data, even when different regularization parameters are explored. Second, the resolution analysis

<footnote n="1">¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Geoazur, F06560 Valbonne, France. E-mail: nocquet@geoazur.unice.fr</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, CNRS, 75238 Paris, France</footnote>
usually remains empirical (for instance, using checkerboard tests) and cannot be made independently of the regularization parameters and inversion schemes. Since models should help to answer questions about physical processes, the regularization approach difficultly provides a definitive answer to questions like ‘how much slip occurred over a given area? or do co- and post-seismic slips overlap?’

The Bayesian approach (e.g. Tarantola & Valette 1982; Yabuki & Matsu'ura 1992; Tarantola 2005; Minson et al. 2013) of inverse problems provides a rigorous framework where prior information on the model parameters can be specified through probability density functions (pdfs). The prior pdfs are then combined to the pdf in the data space in order to derive posterior pdfs on model parameters.

Applied to the finite slip inverse problem, several methods have been developed to derive the posterior pdf and useful associated quantities. Minson et al. (2013) use a Monte Carlo Markov Chain (MCMC) sampler based on a cascading Metropolis algorithm to approximate the joint posterior pdf. Within this approach, any type of prior can be handled, including the less informative prior as uniform pdf, enabling for instance to impose the constraint of slip non-negativity. Although this approach overcomes the issue of the subjectivity of regularization, it also has several drawbacks. First, the posterior pdfs are directly related to the discretization size of the fault, which needs to be chosen *a priori*. This problem can be overcome by adding the number of subfaults as unknown in the inversion and letting the data determine an optimal parametrization (Dettmer et al. 2014), at the price however of additional computer burden. Second, posterior pdfs are approximated using a necessary large number of samples, which also requires large computing power.

Although less generic than the MCMC approach, very fast solutions exist for specific choice of the prior pdf, enabling a fine discretization of the fault. Under the assumption of a linear inverse problem with a Gaussian prior on model parameters, Tarantola & Valette (1982) showed that the joint posterior pdf is also Gaussian. Tarantola & Valette (1982) further provide algebraic expressions for the model expectation, the posterior model covariance and the resolution operator. This implementation has been used for static slip inversion from Global Positioning System (GPS) data by Radiguet et al. (2011) and Villegas-Lanza et al. (2015). Yabuki & Matsu'ura (1992) extended the previous approach for the case where a scaling factor applying to the data covariance matrix and a scaling factor controlling the prior pdf of model parameters are unknown in addition to the slip components. They used Akaike Bayesian information criterion (ABIC) (Akaike 1980) to determine the scaling factors. Fukahata & Wright (2008) also provide algebraic formulations for the posterior covariance and resolution matrices for this case.

For both the Tarantola & Valette (1982) and Yabuki & Matsu'ura (1992) approaches, there is a potential difficulty to express the priors so that the results are independent from the discretization size. In order to solve this problem, Radiguet et al. (2011) use a multivariate Gaussian prior defined through a covariance matrix scaled by a factor ensuring a constant weight independent from the size of the subfaults. Yabuki & Matsu'ura (1992) show that for different discretization or basis functions, the ABIC approach provides different optimal scaling factors for the prior, but the resulting slip model remains almost the same and is independent from the discretization. A more problematic issue is that a Gaussian prior allows negative slip. Negative slip can be significant in some configurations, requiring to artificially modify the priors so that too much unphysical negative slip is avoided. One view is that negative slip, if significantly greater than estimation errors, reveals problems in the physical model and/or an improper error model (Yagi & Fukahata 2011). For instance, Yagi & Fukahata (2011) empirically find that a more realistic error model accounting for uncertainty in the Green's function results in no negative slip. The opposite view adopted by, for example, Minson et al. (2014) or Nocquet et al. (2014) is that non-negative slip is an information relying on reasonable physical assumption that should be accounted for as a prior in the inversion rather than used as a validation check during the post-analysis.

Defining priors as reflecting the true state of knowledge about slip variation is not straightforward. Because stress must remain finite close to the fault, slip must somehow be smooth (Yabuki & Matsu'ura 1992). So, adding smoothing constraints as a prior seems reasonable from a physical point of view. However, translating a physical smoothness constraint into a prior pdf is difficult. Minson et al. (2013) adopt the radical position that no information about smoothness is available and does not use any smoothing, except the one imposed by the discretization, implicitly imposing constant slip at each subfault, usually chosen to be relatively large to ensure a reasonable resolution. Oppositely, Radiguet et al. (2011) and Yabuki & Matsu'ura (1992) let the data defining the level of smoothing either empirically through an L-curve or using the ABIC criterion. In both latter cases, the level of smoothing corresponds to an average over the fault. A consequence of this approach is that the posterior covariance does not quantify any more the true state of knowledge about the estimated parameters because it also includes a contribution coming from the regularization constraints required to obtain a physically meaningful model and avoiding overfitting the observations.

Here, I adopt the view that non-negativity or bound constraints are reasonable from a physical point of view: earthquakes or slow slip events release stress previously accumulated and must induce slip in the direction opposite to the relative plate motion. Afterslip is driven by the response of areas nearby the rupture to the stress increment induced by the coseismic slip (e.g. Perfettini et al. 2010; Avouac 2015) and therefore should be in the same sense as the earthquake itself. Bounds on the interseismic steady slip defined by the relative motion of plates is also reasonable under the assumptions of the backslip model (Savage 1983; Kanda & Simons 2010).

I therefore describe a Bayesian approach that allows the non-negativity or bound constraints to be imposed as a prior as in Minson et al. (2013) and Minson et al. (2014), but where posterior pdfs can be calculated numerically without the need of performing an MCMC sampling. I study both cases of imposing or discarding smoothness of slip as a prior. The central idea in this paper is that if we use the Tarantola & Valette (1982) approach but now restrict the model space to allow positive or bounded slip only, that is if the prior pdf is a multivariate truncated Gaussian, then the joint posterior pdf is also a multivariate truncated Gaussian. Hence, I can use recent advances in truncated multivariate normal (TMVN) probability calculations (Genz & Bretz 2009, see Appendix A) to derive posterior marginal pdf, variance, covariance and other statistics efficiently. I show that the case of uniform independent priors (i.e. no regularization) can be approximated using this approach.
The main advantage of this method is that various quantities can be obtained independently one from each other without the need of running a full MCMC sampling before getting the results.

The paper is organized as follows. First, I focus on the case where the prior information pdf on model parameters is specified as a TMVN. I show that the form of the posterior pdf is similar to the unbounded untruncated Gaussian case. Second, I show that every marginal pdf can be calculated semi-analytically independently of each other. Third, I show that the maximum posterior (MAP) model can be obtained numerically using non-negative or bounded least-squares algorithm. Fourth, I investigate the case where the prior is chosen as being uniform over a set of intervals as classically done in Bayesian inversion when no better prior information is available. I show that for the overdetermined problem, the posterior pdf can be exactly determined and exact marginal pdf can be obtained. For the underdetermined case, approximate solutions can still be derived. The method is illustrated through a simple synthetic 2-D case. Finally, a real case of interseismic modelling is used to demonstrate the equivalence to the Bayesian approach using MCMC sampling.

## 2 LINEAR INVERSE PROBLEM WITH TRUNCATED GAUSSIAN PRIORS

### 2.1 Problem setup and notations

I consider the case of a linear space model $\mathcal{M}$ with dim($\mathcal{M}$) = $p$, corresponding to the unbounded inverse problem ($\mathcal{M} = \mathbb{R}^p$), the case where $\mathcal{M}_+ = (\mathbb{R}^+)^p$ for the inverse problem with non-negativity constraints corresponding to the static inversion of coseismic, post-seismic or transient slow slip and $\mathcal{M}_b$ the space model associated with the bounded problem for interseismic modelling, where $\mathcal{M}_b = [0; b_i]^p$ with $b^i$ being the plate convergence rate at subfault $i$. For the unbounded problem, I consider the case where the a priori information about the unknown model vector $m \in \mathcal{M}$ is a sample of a known Gaussian probability density $\rho_{\mathcal{M}}(m)$ with mean $m_0$ and associated covariance $C_m$, that is, $\rho_{\mathcal{M}}(m) = K_m \exp(-\frac{1}{2}(m - m_0)^T C_m^{-1} (m - m_0))$ with $K_m = \frac{1}{\sqrt{2\pi}^p \sqrt{|C_m|}}$, $|C_m|$ being the determinant of $C_m$. Noting $2P(m) = (m - m_0)^T C_m^{-1} (m - m_0)$, for the non-negativity inverse problem, we have

$$ \rho_{\mathcal{M}_+}(m) = K_m^+ \exp(-P(m)) $$

and for the bounded case,

$$ \rho_{\mathcal{M}_b}(m) = K_m^b \exp(-P(m)), $$

where $K_m^+$ and $K_m^b$ are constant so that the integrals of $\rho_{\mathcal{M}_+}(m)$ over $\mathcal{M}_+$ and $\rho_{\mathcal{M}_b}(m)$ over $\mathcal{M}_b$ are equal to 1.

I consider a linear data space $\mathbb{D}$ with dim($\mathbb{D}$) = $n$. The observation vector $d_{\text{obs}} \in \mathbb{D}$ is a sample of a known Gaussian probability density $\rho_{\mathbb{D}}(d)$ with mean $d_{\text{obs}}$ and associated covariance $C_d$, that is, $\rho_{\mathbb{D}}(d) = K_d \exp(-\frac{1}{2}(d - d_{\text{obs}})^T C_d^{-1} (d - d_{\text{obs}}))$ with $K_d = \frac{1}{\sqrt{2\pi}^p \sqrt{|C_d|}}$, $|C_d|$ being the determinant of $C_d$. Finally, I consider the linear case where observations are related to parameters by a model matrix $G$, $d = G(m) = G_m$ and note $2L(m) = (G_m - d_{\text{obs}})^T C_m^{-1} (G_m - d_{\text{obs}})$. In the following, I will assume that $C_d$ correctly represents the pdf over the data space. In other words, $C_d$ not only represents errors on the data, but also accounts for the physical model prediction errors, coming from our imperfect knowledge of the fault geometry and of the Earth's structure convolved with the slip distribution. The reader is referred to Duputel et al. (2014) and Yagi & Fukahata (2011) for a discussion and methodology to evaluate the covariance matrix associated with model prediction errors.

### 2.2 Joint posterior probability density function

The general form of the joint posterior pdf is the product of the prior pdf with the likelihood function $\rho_{\mathbb{D}}(G(m))$ (Tarantola 2005, p. 35)

$$ \sigma_{\mathcal{M}}(m) = k \rho_{\mathcal{M}}(m) \rho_{\mathbb{D}}(G(m)), \qquad (1) $$

which reduces for the linear case considered here to

$$ \sigma_{\mathcal{M}}(m) = k \rho_{\mathcal{M}}(m) \rho_{\mathbb{D}}(Gm). \qquad (2) $$

Replacing $\rho_{\mathcal{M},\mathcal{M}_+,\mathcal{M}_b}(m)$ and $\rho_{\mathbb{D}}(Gm)$ by their explicit expressions and defining

$$ \begin{aligned} 2S(m) &= 2L(m) + 2P(m) \\ &= (G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}) + (m - m_0)^T C_m^{-1} (m - m_0) \end{aligned} \qquad (3) $$

we have for the non-negativity case

$$ \sigma_{\mathcal{M}_+}(m) = k_+ \exp(-S(m)) \qquad (4) $$

with $k_+ = 1/\int_{\mathcal{M}_+} \exp(-S(m))dm$, and for the bounded problem

$$ \sigma_{\mathcal{M}_b}(m) = k_b \exp(-S(m)) \qquad (5) $$

with $k_b = 1/\int_{\mathcal{M}_b} \exp(-S(m))dm$.
Eqs (4) and (5) show that the joint posterior pdf for the non-negative and bounded inverse problems specified by truncated Gaussian priors has the same expression as for the unbounded fully Gaussian case, except that the normalization constant and the definition domain for $m$ are changed.

I show in the Appendix B1 that $2S(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})^T + K_o$, where $K_o = (G_{m0} - d_{obs})^T(GC_mG^T + C_d)^{-1}(G_{m0} - d_{obs})$ (eq. B1). $\tilde{m}$ and $C_{\tilde{m}}$ are the expectation and posterior covariance matrix for the linear untruncated Gaussian inverse problem, respectively. $\tilde{m}$ and $C_{\tilde{m}}$ are given by any of the formulae from Tarantola & Valette (1982):

$$ \tilde{m} = (G^T C_d^{-1} G + C_m^{-1})^{-1} (G^T C_d^{-1} d_{obs} + C_m^{-1} m_0) \\ = m_0 + C_m G^T (G C_m G^T + C_d)^{-1} (d_{obs} - G m_0) \qquad (6) $$

$$ C_{\tilde{m}} = (G^T C_d^{-1} G + C_m^{-1})^{-1} \\ = C_m - C_m G^T (G C_m G^T + C_d)^{-1} G C_m. \qquad (7) $$

We then have for the non-negativity case

$$ \sigma_{M_+}(m) = k_+ \exp\left(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})\right) \exp\left(-\frac{K_0}{2}\right) \\ = K_+ \exp\left(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})\right) \qquad (8) $$

with $K_+ = k_+ \exp(-K_0/2)$, and for the bounded case

$$ \sigma_{M_b}(m) = K_b \exp\left(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})\right) \qquad (9) $$

with $K_b = k_b \exp(-K_0/2)$.

Eqs (8) and (9) show that the joint posterior pdf for the non-negative and bounded linear inverse problems specified by truncated Gaussian priors is a truncated Gaussian. The integrals $\int_{M_+} \exp(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})) dm$ and $\int_{M_b} \exp(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})) dm$ can be obtained by the approximation from Genz & Bretz (2009) so that eqs (8) and (9) allow us to calculate the value of the posterior pdf for any vector $m \in M_+$ or $m \in M_b$.

## 2.3 Marginal probability density functions

1-D or multidimensional marginal pdfs for the untruncated Gaussian case are trivial since they are Gaussian. In order to derive the marginal pdf for the truncated Gaussian case, I partition the vectors $m$ and $\tilde{m}$ into two subvectors $m = \begin{bmatrix} m_1 \\ m_2 \end{bmatrix}$, $\tilde{m} = \begin{bmatrix} \tilde{m}_{11} \\ \tilde{m}_{12} \end{bmatrix}$. The associated partitioning for the posterior covariance matrix is noted $C_{\tilde{m}} = \begin{bmatrix} C_{\tilde{m}_{11}} & C_{\tilde{m}_{21}} \\ C_{\tilde{m}_{12}} & C_{\tilde{m}_{22}} \end{bmatrix}$.

Using the decomposition described in Appendix C, we have

$$ (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) = Q_1(m_1) + Q_2(m_1, m_2) \qquad (10) $$

with

$$ Q_1(x) = (x - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (x - \tilde{m}_1) $$

$$ Q_2(x_1, x_2) = (x_2 - b(x_1))^T A^{-1} (x_2 - b(x_1)) $$

$$ A = C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^T C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}} $$

$$ b(x) = \tilde{m}_2 + C_{\tilde{m}_{12}}^T C_{\tilde{m}_{11}}^{-1} (x - \tilde{m}_1). $$

Taking the bounded inverse problem case, the joint pdf can now be written as

$$ \sigma_{M_b}(m) = K_b \exp(-\frac{1}{2} Q_1(m_1)) \exp(-\frac{1}{2} Q_2(m_1, m_2)). $$

Noting $M_{b1}$ and $M_{b2}$ the subspaces of $M_b$ corresponding to the domain of definition for $m_1$ and $m_2$ respectively, the marginal pdf of $m_1$ is

$$ \sigma_{M_{b1}}(m_1) = K_b \exp(-\frac{1}{2} Q_1(m_1)) \int_{M_{b2}} \exp(-\frac{1}{2} Q_2(m_1, m_2)) dm_2 \\ = K_b \exp\left(-\frac{1}{2}(m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1}(m_1 - \tilde{m}_1)\right) \\ \qquad \int_{M_{b2}} \exp\left(-\frac{1}{2}(m_2 - b(m_1))^T A^{-1}(m_2 - b(m_1))\right) dm_2. \qquad (11) $$
For the non-negative inverse problem, similarly we have

$$ \sigma_{M_{+1}}(m_1) = K_b \exp\left(-\frac{1}{2}(m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)\right) \int_{M_{+2}} \exp\left(-\frac{1}{2}(m_2 - b(m_1))^T A^{-1} (m_2 - b(m_1))\right) dm_2. \quad (12) $$

The integral term $\int_{M_{p2}, M_{+2}} \exp(-\frac{1}{2}(m_2 - b(m_1))^T A^{-1} (m_2 - b(m_1)) dm_2$ can again be evaluated using the approximation from Genz & Bretz (2009) for any $m_1$ in order to draw the marginal pdf for $m_1$. As noted by Horrace (2005) and Cartinhour (1990), eqs (11) and (12) show that the marginal pdfs of a TMVN are not TMVN anymore. This is because the term $\int \exp(-\frac{1}{2}(m_2 - b(m_1))^T A^{-1} (m_2 - b(m_1)) dm_2$ is not constant, but depends on $m_1$ and was named a 'skewness' function by Cartinhour (1990). At this stage, note that $m_1$ is any subset of the components of $m$. Eqs (11) and (12) can therefore be used to derive single (1-D) posterior pdf but also posterior marginal joint pdf for any subset of components of $m$. From eqs (11) and (12), we can obtain the conditional distribution of $m_2$ given $m_1$ which is

$$ \sigma_{m_2|m_1}(m_2|m_1) = \frac{\sigma_{M_b}(m_1, m_2)}{\sigma_{M_{b1}}(m_1)}. $$

Horrace (2005) demonstrates that the conditional pdf is a TMVN.

Although eqs (11) and (12) can be used to numerically derive posterior individual variance or posterior bivariate marginal pdf and covariances, more efficient algorithms have been developed for these specific cases. Leppard & Tallis (1989) proposed an algorithm for a more direct evaluation of the mean and covariance of single TMVN, extended by Manjunath & Wilhelm (2009) for the double-truncated case. A summary of TMVN probability calculation methods and information for practical implementation are given in Appendix A.

## 2.4 Numerical solution for the maximum a posteriori (MAP) probability estimates

From the posterior joint pdf given in eqs (8) and (9), the value $\tilde{m}_{MAP}$ which maximizes the posterior joint pdf verifies

$$ \min \|S(m)\|_2 \text{ with } 0 \le m_i, i \in \{1, p\} \quad (13) $$

$$ \min \|S(m)\|_2 \text{ with } 0 \le m_i \le b_i, i \in \{1, p\}, \quad (14) $$

for the non-negative and bounded problems, respectively. $\tilde{m}_{MAP}$ can be numerically obtained using the non-negative least-squares algorithm from Lawson & Hanson (1974) for the non-negative case or the Bounded-Variable Least-Squares (BVLS) algorithm from Stark & Parker (1995) for the bounded problem. Both algorithms do not directly include a covariance matrix, requiring a change of variable (Snieder & Trampert 1999). For instance, the BVLS algorithm solves the following problem

$$ \min \|A_m - B\|_2 \text{ with } lb_i \le m_i \le ub_i, i \in \{1, n\}. \quad (15) $$

I define the square roots $C_d^{-1/2}$ and $C_m^{-1/2}$ so that

$$ (C_d^{-1/2})^T C_d^{-1/2} = C_d^{-1} \quad (16) $$

$$ (C_m^{-1/2})^T C_m^{-1/2} = C_m^{-1}. \quad (17) $$

Inserting $C_d^{-1/2}$ and $C_m^{-1/2}$ into eq. (3), we have

$$ \begin{aligned} 2S(m) &= (G_m - d_{obs})^T C_d^{-1} (G_m - d_{obs}) + (m - m_0)^T C_m^{-1} (m - m_0) \\ &= (C_d^{-1/2} G_m - C_d^{-1/2} d_{obs})^T (C_d^{-1/2} G_m - C_d^{-1/2} d_{obs}) + (C_m^{-1/2} m - C_m^{-1/2} m_0)^T (C_m^{-1/2} m - C_m^{-1/2} m_0). \end{aligned} $$

Taking $A$ and $B$ as

$$ A = \begin{bmatrix} C_d^{-1/2} & 0 \\ 0 & C_m^{-1/2} \end{bmatrix} \begin{bmatrix} G \\ I_p \end{bmatrix} = \begin{bmatrix} C_d^{-1/2} G \\ C_m^{-1/2} \end{bmatrix} \quad (18) $$

$$ B = \begin{bmatrix} C_d^{-1/2} & 0 \\ 0 & C_m^{-1/2} \end{bmatrix} \begin{bmatrix} d_{obs} \\ m_0 \end{bmatrix} = \begin{bmatrix} C_d^{-1/2} d_{obs} \\ C_m^{-1/2} m_0 \end{bmatrix}, \quad (19) $$

where $I_p$ is the identity matrix of order $p$, $2S(m)$ becomes

$$ 2S(m) = (A_m - B)^T (A_m - B). \quad (20) $$

Taking $lb_i = 0$ and $ub_i = b_i$, $i \in \{1, p\}$, eq. (20) shows that solving for eq. (14) or (15) is equivalent using $A$ and $B$ defined by eqs (18) and (19).

For the non-negative case, Lawson & Hanson (1974) show that (1) there is always a solution to eq. (20), (2) the solution is unique and (3) the solution is obtained by finding the subset of $m_i$ which are set to 0, the other $m_i$ being the results of a classical least-squares solution. Stark & Parker (1995) indicate that these three properties also hold for the bounded case. Because of the equivalence of eqs (14) and (15), for
linear inverse problems with truncated Gaussian priors, the results of Lawson & Hanson (1974) and Stark & Parker (1995) demonstrate that the MAP always exists and is unique.

In pratice, $C_d^{-1/2}$ and $C_m^{-1/2}$ can be obtained using the Cholesky decomposition and taking the inverse. $C_m$ may sometimes be a badly conditioned matrix. Alternative approaches can still be used in that case to obtain approximate square root of the inverse (Benzi et al. 1996, 2000; Kharchenko et al. 2001).

# 3 EQUIVALENCE TO BAYESIAN INVERSION WITH INDEPENDENT BOUNDED UNIFORM PRIORS

So far, all results have been derived for a prior pdf being a TMVN. I now describe how this approach can handle the case where prior pdfs are specified using the less informative pdf being a constant over a set of intervals.

## 3.1 Exact solution for the overdetermined inverse problem

The joint pdf for bounded uniform priors is

$$ \rho_{\mathcal{M}_b} = \mu_{\mathcal{M}_b}(m) = \prod_{i=1}^{p} \mu_{[0,b_i]}(m_i) \qquad (21) $$

with $\mathcal{M}_b = [0; b_i]_{i \in \{1, p\}}$ as before and $\mu_{[0,b_i]}(m_i) = \frac{1}{b_i}$.

Then, $\sigma_{\mathcal{M}_b}(m) = k \rho_{\mathcal{M}_b}(m) \rho_{\mathcal{D}}(G_m) = k \prod_{i=1}^{p} \mu_{[0,b_i]}(m_i) K_d \exp(-\frac{1}{2}(G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}})) = K_d \exp(-\frac{1}{2}(G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}))$ because $k = \prod_{i=1}^{p} b_i$.

From Appendix B2, for the overdetermined case, we can write

$$ (G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) + K_o $$

with $C_{\tilde{m}} = (G^T C_d^{-1} G)^{-1}$, $\tilde{m} = C_{\tilde{m}} G^T C_d^{-1} d_{\text{obs}}$. Thus

$$ \sigma_{\mathcal{M}_b}(m) = K_{uo} \exp \left( -\frac{1}{2} (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) \right) \qquad (22) $$

with $K_{uo} = \exp(-\frac{1}{2} K_o K_d)$.

Eq. (22) shows that for a linear bounded inverse problem with uniform independent priors, the posterior joint pdf is exactly a TMVN. As a consequence, 1- and 2-D marginal pdf can be obtained using the results described in Section 2.3. The individual means and variances and covariance functions for pair of parameters can be obtained using the algorithms implemented by Wilhelm (2015). The MAP can also be obtained using the results from Section 2.4.

## 3.2 Approximate solution for the underdetermined problem

For the underdetermined inverse problem, $G^T C_d^{-1} G$ is no longer invertible and the results from the previous paragraph cannot be used. Nonetheless, the results obtained for truncated Gaussian priors can be used to approximate posterior pdf of the bounded uniform independent priors inverse problem. The basic idea is that (1) from eq. (9), we have an exact solution for the bounded inverse problem with TMVN prior and (2) individual marginal uniform prior pdfs are rectangle functions that can be approximated by a truncated Gaussian as shown in Fig. 1.

From

$$ \rho_{\mathcal{M}_b}(m) = K_b \exp \left( -\frac{1}{2} (m - m_0)^T C_m^{-1} (m - m_0) \right) \qquad (23) $$

taking $C_m = \sigma^2 I_p$, $C_m^{-1} = 1/\sigma^2 I_p$ and $m_0 = (m_{0i})$, $i \in \{1, p\}$, with $m_{0i} = \frac{b_i}{2}$, the prior pdf becomes

$$ \rho_{\mathcal{M}_b}(m) = \prod_{i=1}^{p} \frac{\exp(-\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2})}{\int_0^{b_i} \exp(-\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2}) dm_i} \qquad (24) $$

By a shift of $b_i/2$,

$$ \int_0^{b_i} \exp \left( -\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2} \right) dm_i = \int_{-b_i/2}^{b_i/2} \exp \left( -\frac{1}{2} \frac{x^2}{\sigma^2} \right) dx = \operatorname{erf} \left( \frac{\sqrt{2}b_i}{4\sigma} \right) \sqrt{2\pi}\sigma. $$

If $\sigma \gg \max_{i \in \{1, p\}}(b_i/2)$, then $\sigma \gg \max_{i \in \{1, p\}}(|m_i - m_{0i}|)$ and $\exp(-\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2}) \approx 1 + \frac{m_i - m_{0i}}{\sigma^2} \exp(-\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2})$. Since $\operatorname{erf}(x) \to \frac{2}{\sqrt{\pi}}x$ when $x \to 0$, $\int_0^{b_i} \exp(-\frac{1}{2} \frac{(m_i - m_{0i})^2}{\sigma^2}) dm_i \approx b_i$, for $\sigma \gg \max(\frac{b_i}{2})$. So $\rho_{[0;ub_i]}(m_i) \approx \mu_{[0;ub_i]}(m_i)$ and $\rho_{\mathcal{M}_b}(m) \approx \mu_{\mathcal{M}_b}(m)$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Approximation of a rectangular function (grey curve) by truncated Gaussian functions. As the parameter $\alpha = \frac{\sigma}{b}$ increases, better approximation of the rectangle function is obtained.</figcaption></figure>

This result shows that for the bounded inverse problem with uniform priors often used in Bayesian inversions, an approximate solution of the posterior pdf can be obtained using $C_m = \sigma^2 I_p$ for large $\sigma(\sigma \gg \max(\frac{b_i}{2}))$ and an a priori model vector $m_0 = (m_{0i})$, $i \in \{1, p\}$, with $m_{0i} = \frac{b_i}{2}$. This result directly takes benefit that the derivative of a Gaussian at its centre is null. It offers an alternative to Bayesian approaches for fast computation of 1- and 2-D marginal pdf, means, variances, covariances and MAP determination using the results described in Sections 2.3 and 2.4.

## 3.3 Convergence of the TMVN approximation to uniform priors and numerical stability

The results from the previous section show that the parameter $\alpha = \frac{\sigma}{\max(\frac{b_i}{2})}$ controls how the TMVN $\rho_{M_b}(m)$ in eq. (24) approximates the uniform prior $\mu_{M_b}(m)$ defined in eq. (21). From the mathematical point of view, using a large $\alpha$ provides a better approximation to the results from an inversion with true uniform priors. However, from a numerical point of view, we see that the use of eq. (9) requires $C_m$ to be evaluated. As $\alpha$ increases, $C_m^{-1}$ tends to be smaller and smaller and $(G^T C_d^{-1} G + C_m^{-1})$ gets more and more badly conditioned, leading to enhanced numerical errors and oscillatory behavior of the obtained posterior pdf. Therefore, there is a trade-off between the numerical accuracy and the mathematical accuracy.

The convergence of the TMVN approximation towards the results from a true inversion using uniform priors can be empirically observed by the fact that for $\alpha$ greater than a threshold value, there is no change in the shape of the individual posterior marginal pdf, as illustrated in the next Section and Fig. 2. Starting with very small $\alpha$ results in narrow Gaussian-like shapes, which progressively spread out as $\alpha$ increases, before reaching a stable shape. This evolution in change of the marginal pdf is a consequence of the relative weight of the prior and the likelihood in eq. (2). For small $\alpha$, the prior TMVN dominates and controls the shape of the posterior pdf. For large $\alpha$, the likelihood dominates. Therefore, running the inversion with increasing values of $\alpha$ provides an empirical way to obtain the best trade-off between mathematical and numerical accuracies.

Another quantity useful to evaluate the convergence is the individual variance reduction, defined for each model parameter $m_i$ by

$$vr_i(\alpha_i) = \frac{\rho_i^2(\alpha_i) - \gamma_i^2(\alpha_i)}{\rho_i^2(\alpha_i)}, \qquad (25)$$

where $\rho_i^2(\alpha_i)$ and $\gamma_i^2(\alpha_i)$ are respectively the prior and posterior variance for parameter $m_i$ with $\alpha_i = \frac{\sigma_i}{b_i}$, $\sigma_i^2$ being the prior variance now taken different for each $m_i$. The variance reduction quantifies the information gain for an individual parameter. When $\alpha_i$ increases, $vr_i(\alpha_i)$ increases and tends towards the variance reduction that would be obtained for an inversion using true uniform priors. There is also a threshold value $\alpha_i^{\text{threshold}}$ for each parameter $i$ where the variance reduction only increases marginally. This dependence of the variance reduction to $\alpha_i$ can be used to design a simple algorithm providing an optimal trade-off between the mathematical and numerical accuracy. Taking benefit from the speed of calculation of marginal pdf using eq. (12), we can first run a first set of simulations with all $\alpha_i$ increasing equally. From this simulation, for ‘good’ parameters, we see a threshold value $\alpha_i^{\text{threshold}}$ of $\alpha$ for each parameter $i$ where the variance reduction only increases marginally. For ‘bad’ parameters, numerical instabilities occur before convergence is achieved. I then fix $\alpha_{i0} = \alpha_{i0}^{\text{threshold}}$ where $i0$ is the parameter having the smallest $\alpha_i^{\text{threshold}}$, that is, $i0 = \text{argmin}(\alpha_i^{\text{threshold}})$. Then a new series of equally increasing $\alpha_i$ is run, leading to fix a second $\alpha_{i1}$. The procedure is iterated to fix all $\alpha_i$. The validation of this approach is illustrated in the case of interseismic modelling in Section 5.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the equivalence of the TMVN method to the Bayesian approach for a simple synthetic case $dim(M_b) = 2$ (see text for details). (a) Likelihood $L(m)$ over the full space as a function of parameters $m_1$ and $m_2$. The white box indicates $M_b$, that is, the bounded domain considered for the inverse problem. (b) Joint posterior pdf given by eq. (5). (c–e) Approximate joint posterior pdf using the TMVN method with $\alpha = 0.5, 1.0$ and $8.0$. As $\alpha$ increases, a better approximation of the true posterior joint pdf is obtained. (f) Approximated posterior joint pdf using Monte Carlo sampling. Small dots show the sampling whose densities are shown by colour codes by $0.05 \times 0.05$ wide bins. (g and h) 1-D marginal pdf for $m_1$ and $m_2$ parameters. Yellow histograms are the marginal pdf derived from the MCMC sampling. For $\alpha$ larger than 4, the recovered pdfs are visually indistinct from the true pdf and are not represented.</figcaption></figure>

4 SYNTHETIC CASE

In order to visualize the methodology previously described, I use a simple synthetic case where $G = \begin{bmatrix} -7 & -4 \\ 1 & 10 \\ 2 & -11 \end{bmatrix}$, $d = \begin{bmatrix} 10 \\ 3 \\ -5 \end{bmatrix}$, $C_d = 5^2 I_3$,
$M_b = [0, 1] \times [0, 1]$, with uniform priors. These values are arbitrary and have been chosen for visualization purposes. Fig. 2(a) showing the likelihood $L(m) = \frac{1}{2}(G_m - d_{obs})^T C_d^{-1} (G_m - d_{obs})$ over a domain wider than $M_b$ indicates that $L(m)$ has a minimum located outside $M_b$. For this 2-D inverse problem, we can visualize the true posterior pdf over $M_b$ using eq. (5), shown in Fig. 2(b). We then compare the true posterior pdf with its approximation provided by eq. (8) obtained for $\sigma = \alpha \max(\frac{b_i}{2})$, with $\alpha = 0.5, 1.0$ and $8.0$. Figs 2(c–e) show that the approximated posterior pdf rapidly converges towards the true pdf. We further compare the results with the posterior pdf derived from a Monte Carlo sampling. The Monte Carlo results have been obtained using the PyMC package (https://github.com/pymc-devs/pymc), using
<table><caption>**Table 1.** Mean, standard deviation and maximum posterior of parameters ($m_1/m_2$) for the different methods described in the text. std, standard deviation; MAP, $m_1, m_2$ values corresponding to the maximum posterior joint pdf.</caption><thead><tr><th>Method</th><th>Mean</th><th>Std</th><th>MAP</th></tr></thead><tbody><tr><td>True pdf</td><td>0.229/0.328</td><td>0.200/0.219</td><td>0.000/0.190</td></tr><tr><td>App. pdf α = 0.5</td><td>0.251/0.346</td><td>0.200/0.211</td><td>0.000/0.282</td></tr><tr><td>App. pdf α = 1.0</td><td>0.233/0.332</td><td>0.201/0.217</td><td>0.000/0.219</td></tr><tr><td>App. pdf α = 8.0</td><td>0.229/0.328</td><td>0.200/0.219</td><td>0.000/0.190</td></tr><tr><td>Monte Carlo</td><td>0.228/0.326</td><td>0.199/0.218</td><td>0.000/0.186</td></tr></tbody></table>

$5 \times 10^5$ samples, with burn and thinning values of 1000 and 10 (Fig. 2f). Figs 2(g) and (h) show the posterior marginal pdf obtained for each of the two parameters using eq. (11) to be compared with the Monte Carlo approach and the true value. Again, we see the rapid convergence of the approximate pdf towards the true value. For $\alpha = 3$, the marginal pdf obtained using the TMVN approximation is undistinguishable from the true pdf.

Table 1 shows the mean and the standard deviation numerically derived from the marginal pdf, together with the values of $m_1$ and $m_2$ corresponding to the maximum value of the posterior joint pdf (MAP). For the approximate method, the MAP is obtained using the method described in Section 2.4 with the BVLS algorithm from Stark & Parker (1995). Again, we see a very rapid convergence towards the true value. For $\alpha = 8.0$, we have convergence at a level better than $10^{-3}$, while the Monte Carlo approach still suffers from small (a few $10^{-3}$) errors. In terms of speed, the Monte Carlo inversion took 63 s for $5 \times 10^5$ samples against 0.15 s to calculate all the marginal pdfs with 20 points and 0.0004 s to obtain the MAP.

## 5 AN UNDERDETERMINED BOUNDED PROBLEM WITH INDEPENDENT UNIFORM PRIORS: INTERSEISMIC MODELLING IN CENTRAL PERU

In order to further validate the methodology, I apply the TMVN method to a real case of interseismic modelling in Central Peru, where the Nazca plate converges towards the peruvian Inca sliver (Nocquet et al. 2014) at a velocity of $\sim 57$ mm yr$^{-1}$ (Villegas-Lanza et al. 2016). Central Peru delimits a $\sim 500$ km long segment of high interseismic coupling from a depth of $\sim 50$ km possibly extending updip to the trench as found from seafloor geodesy measurements (Gagnon et al. 2005) with lower interseismic coupling north and south of it (Villegas-Lanza et al. 2016). Rather than providing a new model for this area, the aim here is to evaluate the TMVN methodology on a real case. I selected 13 horizontal GPS velocities providing 26 observations and discretized a 700 km long fault following the trench (strike 35°W) with a dip of 15° into 30 70 × 70 km square subfaults (Fig. 3). I use a fixed rake of 113° in agreement with the relative motion between the Nazca plate and the Inca sliver as defined in Villegas-Lanza et al. (2016). The inverse problem is underdetermined with $dim(M_b) = 30$ and $dim(D) = 26$ and $M_b = [0, v_i]_{i \in (1,30)}$, where $v_i$ is the relative plate motion calculated a subfault $i$. $C_d$ is simply taken as the diagonal matrix containing the individual variance of velocity components coming the GPS analysis.

The backslip approach (Savage 1983) is used and the Green's function relating the GPS velocity components to the unit slip at a subfault in a direction opposite to the plate convergence is calculated for a semi-infinite homogeneous elastic half-space (Okada 1992). The Bayesian approach was implemented using the PyMC package (https://github.com/pymc-devs/pymc), using $10^7$ samples, with burn and thinning values of 5000 and 10 respectively, leaving $9.95 \times 10^5$ posterior samples. Uniform prior bounded between 0 and the plate velocity were specified for the 30 parameters. For the TMVN approach, the algorithm described in Section 3.3 has been applied. A convergence tolerance of 2 per cent has been used, leading to values for $\alpha_i$ between 1.4 and 8.7.

In the following, we compare the TVMN results to the MCMC results for various quantities usually presented for Bayesian inversions. The aim is to assess the ability of the TMVN method to reproduce Bayesian results obtained using MCMC sampling for independent uniform priors for an underdetermined problem. I chose this case because it is the most difficult: the underdetermined case prevents an exact solution and uniform priors require the priors to be also approximated. When comparing the results between the TMVN and Bayesian MCMC inversions, we need to remember that the MCMC is also an approximation to the true quantities. We will see that part of the few differences noted can partly be attributed to errors in the MCMC with respect to the truth.

### 5.1 1- and 2-D marginal pdf

Fig. 4 shows the comparison between the 1-D marginal pdfs obtained from the two methods. The 1-D pdfs were obtained using the dtmvnorm.marginal function from the tmvtnorm R package (Wilhelm 2015) which implements the method proposed by Cartinhour (1990). Fig. 4 shows that the chosen inverse problem includes a great diversity in terms of parameter resolution. Some parameters are well defined and centred around a value between the bounds and show nearly Gaussian shape (parameters 12, 13, 14, 17, 22 and 23), while for other parameters, the maximum of the pdf has its maximum close to either 0 per cent of coupling (e.g. 21) or 100 per cent (e.g. 1–6). Parameters 7, 8 and 9 are badly determined, as indicated by a widely spread pdf. We can note that the resolution indicated by the marginal pdf is not a direct function of the distance between the subfault and the GPS sites as indicated by the very good resolution for subfaults 1–5. In general, the TMVN and MCMC approaches show a very good agreement in terms of the shape and the values of the marginal pdf. The agreement is
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Inverse problem setup for interseismic modelling in central Peru. White squares indicate individual subfaults together with their number in the inversion. Yellow arrows are interseismic velocities with respect to the overriding plate, with their 95 per cent confidence level ellipses.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** 1-D pdf for the coupling coefficient in per cent. Yellow histograms show the MCMC results. Red curves show the marginal pdf from the TMVN approximation. Numbers indicate the subfault index as in Fig. 3</figcaption></figure>

even better for well resolved parameters. Differences increase as parameter are less resolved (e.g. 7, 9 and 29), but are always less than 20 per cent of the pdf value.

Fig. 5 shows the posterior pdf for pair of parameters, that is, the 2-D covariance functions. The posterior 2-D covariance functions were obtained using the dtmvnorm.marginal2 function from the tmvtnorm R package (Wilhelm 2015) which implements the method proposed by Leppard & Tallis (1989), modified by Wilhelm (2015). Fig. 5 shows that the TMVN method is also able to provide the covariance information with a similar accuracy as the MCMC approach. Slight differences can be noted for pair of less resolved parameters (for instance 00-07, 07-09, 07-29 and 09-29), but the shape of the 2-D covariance function always remains preserved.

## 5.2 Mean, maximum posterior and median models

Fig. 6 shows the mean, MAP and median models obtained using the MCMC and TMVN approaches in a map view. Fig. 7 shows the comparison in an MCMC value versus TMVN value plot. Both figures show that the mean and median models obtained for both approaches are the same within 5 per cent. However, for the MAP larger differences are found, the largest being for subfault #7 where the value is 28 per
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 1- and 2-D marginal pdf obtained using the MCMC and TMVN approaches. Subplots on the diagonal show the 1-D pdf as in Fig. 4. Subplots in the upper triangle part are the marginal pdf from the MCMC approach for pairs of parameters. Plots in the lower triangle part of the plot matrix are their equivalent derived using the TMVN approach. Both representations use evaluation of the covariance function over a 30 × 30 points grid. For the non-diagonal subplots, the *x*-axis corresponds to the coupling value for the subfault with the higher index, while the *y*-axis corresponds to coupling value for the subfault with the lowest index. The colour code is from white to black through yellow and red, white being the lowest probability value and black the highest. The colour scale is optimized for each subplot and therefore only reflects relative values within each subplot.</figcaption></figure>

cent for the MCMC approach and 97 per cent for the TMVN approach. Looking at the marginal pdf for subfault #7 (Fig. 4), this value for the MCMC approach looks surprising, because the marginal pdf would rather suggest a maximum posterior value closer to 100 per cent. While the MAP for the TMVN method has been obtained using the method described in Section 2.4, the MAP for the MCMC approach has simply been taken as the sample from the MCMC results having the maximum posterior probability. The MAP value from the MCMC therefore likely represents an error probably due to the partial exploration of parameters by the MCMC sampler for a badly determined parameter, for which different values induce very little change in the likelihood. For the MCMC approach, more reliable MAP estimates are obtained from the maximum of the posterior 1-D marginal pdf that would in this case agree within 2–4 per cent with the TMVN results. The difference between the TMVN and MCMC approaches observed for a badly determined parameter also suggests that the TMVN method can be used as an independent verification for MCMC results.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Results for interseismic coupling inversion. Bayesian MCMC results are shown on the left column and TMVN results on the right column. Top: mean model. Middle: maximum posterior (MAP) model. Bottom: median model.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison between MCMC and TMVN results. For the median, standard deviation, median, MAP and information gain, results from the TMVN approach are along the x-axis and results from the MCMC approach along the y-axis. The grey line indicates deviation of ±5 per cent to the equality between the MCMC and TMVN results. The bottom right subplot shows the 95 per cent confidence level for each parameter. X-axis indicates the subfault number as shown in Fig. 3. Confidence intervals obtained from the TMVN approach are in light pink. Confidence intervals obtained from the MCMC approach are shown by the thin black lines.</figcaption></figure>

## 5.3 Standard deviation, confidence interval and information gain

Fig. 8 shows that the standard deviation derived from both methods are equivalent within 5 per cent. Fig. 8 also shows that the 95 per cent confidence intervals are similar. It is interesting to note that for instance parameter #7, which showed the largest MAP difference between the MCMC and TMVN estimates, shows an identical 95 per cent confidence interval for both method. This suggests that although MCMC and the TMVN methods can sometimes provide differences for the MAP for badly resolved parameters, the inferences made using statistical criterion will not be different.

I finally compare the information gain, here simply defined as the normalized variance reduction using eq. (25). An information gain of 1 (or 100 per cent) indicates that the data enable the parameter to be perfectly resolved while an information gain close to 0 indicates that the data do not help to improve the prior knowledge we had on the parameter. In the chosen example, the information gain is lower close to the trench and at the northern edge of the chosen fault. For instance, the parameter #7 has an extremely low information gain (10 per cent found for both approaches), while parameter #23, close to the GPS data has an information gain better than 90 per cent. We see that the TMVN method usually provides information gain estimates equivalent to the MCMC estimates at the 10 per cent level (maximum difference for parameter #3 is 25 per cent). In general, we also note that the information gain from the TMVN is lower than the MCMC value. This can be understood in the light of Fig. 1. Indeed, the prior variance for a truncated Gaussian function is always smaller than for the true uniform rectangle function, leading to a smaller variance reduction. Therefore, the TMVN approach provides slightly biased (more pessimistic) estimates of the true information gain.

## 5.4 Comparison of computation time

The MCMC approach took 3.5 hr on a 2.9 GHz Intel Core i7 CPU with a RAM of 8 Gb. On the same computer, the TMVN approach took less than 10 s to calculate an individual marginal pdf with 30 points (i.e. 300 s to generate the pdfs shown in Fig. 4). It, however, took in average 6 mn to get the 2-D pdf over 900 points (30 × 30), making the TMVN approach not efficient to derive 2-D pdf compared to the Bayesian approach if a single CPU is used. The main advantage still remains that every parameter it can be calculated ‘on demand’ very quickly. Calculating the mean independently from the pdf took 17 s and getting the MAP took 0.1 s.
<figure><img src="image_10.png" /><figcaption>information gain from MCMC (%)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>coupling model standard deviation from TMVN</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>information gain from TMVN (%)</figcaption></figure>

<figure><img src="image_13.png" /></figure>

**Figure 8.** Standard deviation and information gain from the Bayesian MCMC (left column) and TMVN (right column) approaches. Top: standard deviation on coupling estimates at each subfault. Bottom: information gain (in per cent) for each subfault. Values are colour-coded according to the scale shown on the left of each subplot.

## 6 CONCLUSIONS

This study builds upon the framework of the stochastic approach of inverse problem proposed by Tarantola & Valette (1982). It extends their results for linear inverse problems for which non-negativity or bounds are added as prior information to the model parameters. The most useful result is that 1- or 2-D marginal pdf can be obtained, as well as different quantities characterizing the inversion results. The non-negativity or bounds are reasonable physical assumptions and define by themselves strong regularization constraints to the inverse problem. Incorporating these constraints further changes the resolution analysis and potentially the inferences made from the inversion results. This study also shows that the case of low-informative uniform independent priors can be approximated. The Bayesian approach had been so far addressed by samplings, requiring intensive computer power. The approach proposed in this study reduces by far the required computer power. Among its advantage, it allows the marginal pdf to be obtained independently one from each other, allowing new strategies to be developed in future. For instance, the resolution dependence to size of the subfault in a given subspace of the model can be assessed without the need of running the full inversion. New development in the mathematical field of multivariate normal probabilities might make this approach even faster in future and useful for a range of geophysical and physical inverse problems.

## ACKNOWLEDGEMENTS

I wish to pay tribute to A. Tarantola's contribution to the inverse problem theory. This work benefits from seeds he sown about two decades ago at IPG Paris. I am indebted to B. Valette for a course given at the Mathematics Department, Escuela Politecnica Nacional, Quito and for discussions that helped me to better formalize the ideas presented in this paper. I also thank M. Simons for a fruitful and stimulating exchange during a short visit at Caltech. Finally, I thank two anonymous reviewers for their comments that contributed to improve the paper. Most of the work presented here was done during my staying in Quito and I thank the Instituto Geofísico, Escuela Politécnica Nacional, Quito, Ecuador.
This work has been supported by the Agence National de la Recherche (ANR) under grants REMAKE ANR-15-CE04-0004 and E-POST ANR-14-CE03-0002.

## REFERENCES

Akaike, H., 1980. Likelihood and the Bayes procedure, in *Bayesian Statistics*, pp. 27–41, eds Bernardo, J.M., De Groot, M.H., Lindley, D.V. & Smith, A.F.M., University Press.

Avouac, J.-P., 2015. From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle, *Ann. Rev. Earth Planet. Sci.*, **43**, 233–271.

Benzi, M., Meyer, C.D. & Tùma, M., 1996. A sparse approximate inverse preconditioner for the conjugate gradient method, *SIAM J. Sci. Comput.*, **17**(5), 1135–1149.

Benzi, M., Cullum, J.K. & Tuma, M., 2000. Robust approximate inverse preconditioning for the conjugate gradient method, *SIAM J. Sci. Comput.*, **22**(4), 1318–1332.

Cartinhour, J., 1990. One-dimensional marginal density functions of a truncated multivariate normal density function, *Comm. Stat. - Theory Methods*, **19**(1), 197–203.

Dettmer, J., Benavente, R., Cummins, P.R. & Sambridge, M., 2014. Trans-dimensional finite-fault inversion, *J. Geophys. Int.*, **199**(2), 735–751.

Duputel, Z., Agram, P.S., Simons, M., Minson, S.E. & Beck, J.L., 2014. Accounting for prediction uncertainty when inferring subsurface fault slip, *J. Geophys. Int.*, **197**(1), 464–482.

Fukahata, Y. & Wright, T.J., 2008. A non-linear geodetic data inversion using ABIC for slip distribution on a fault with an unknown dip angle, *J. Geophys. Int.*, **173**(2), 353–364.

Gagnon, K., Chadwell, C.D. & Norabuena, E., 2005. Measuring the onset of locking in the Peru–Chile trench with GPS and acoustic measurements, *Nature*, **520**(2004), 205–208.

Genz, A., 1992. Numerical computation of multivariate normal probabilities, *J. Comput. Graph. Stat.*, **1**(2), 141–149.

Genz, A. & Bretz, F., 1999. Numerical computation of multivariate t-probabilities with application to power calculation of multiple contrasts, *J. Stat. Comput. Simul.*, **63**(4), 103–117.

Genz, A. & Bretz, F., 2009. *Computation of Multivariate Normal and t Probabilities*, Vol. **195**, Springer Science & Business Media.

Genz, A., Bretz, F. & Hothorn, T., 2006, MVTNORM: an R package for computing multivariate normal and t probabilities, quantiles and densities. Available from cran.r-project.org/web/pac.kages/mvtnorm, last accessed May 2017.

Horrace, W.C., 2005. Some results on the multivariate truncated normal distribution, *J. Multivariate Anal.*, **94**(1), 209–221.

Hothorn, T., Bretz, F. & Genz, A., 2001. On multivariate t and Gauss probabilities in R, *SIGMA*, **1**(2), 1–6.

Kanda, R.V.S. & Simons, M., 2010. An elastic plate model for interseismic deformation in subduction zones, *J. geophys. Res.*, **115**(B3), B03405, doi:10.1029/2009JB006611.

Kharchenko, S.A., Kolotilina, L.Y., Nikishin, A.A. & Yeremin, A.Y., 2001. A robust AINV-type method for constructing sparse approximate inverse preconditioners in factored form, *Numer. Linear Algebr. Appl.*, **8**(3), 165–179.

Lawson, C. L. & Hanson, R.J., 1974. *Solving Least Squares Problems*, Prentice-Hall.

Leppard, P. & Tallis, G., 1989. Algorithm as 249: evaluation of the mean and covariance of the truncated multinormal distribution, *J. R. Stat. Soc. Series C (Appl. Stat.)*, **38**(3), 543–553.

Manjunath, B.G. & Wilhelm, S., 2009. Moments calculation for the double truncated multivariate normal density, *SSRN eLibrary*, (1963), 1–11.

Minson, S.E., Simons, M. & Beck, J.L., 2013. Bayesian inversion for finite fault earthquake source models I-theory and algorithm, *J. Geophys. Int.*, **194**(3), 1701–1726.

Minson, S.E. et al., 2014. Bayesian inversion for finite fault earthquake source models - II: The 2011 great Tohoku-oki, Japan earthquake, *J. Geophys. Int.*, **198**(2), 922–940.

Nocquet, J.-M. et al., 2014. Motion of continental slivers and creeping subduction in the northern Andes, *Nature Geoscience*, **7**(8), 287–292.

Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **82**(2), 1018–1040.

Perfettini, H. et al., 2010. Seismic and aseismic slip on the central peru megathrust, *Nature*, **465**(7294), 78–81.

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostoglodov, V. & Cotte, N., 2011. Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero Slow Slip Event, *J. Geophys. Int.*, **184**(2), 816–828.

Savage, J. C., 1983. A dislocation model of strain accumulation and release at a subduction zone, *J. geophys. Res.*, **88**(B6), 4984–4996.

Snieder, R. & Trampert, J., 1999. "Inverse problems in geophysics." Wavefield inversion, *Springer Vienna*, **398**, 119–190.

Stark, P. & Parker, R., 1995. Bounded-variable least-squares: an algorithm and applications, *Comput. Stat.*, **10**, 129–141.

Tallis, G.M., 1961. The moment generating function of the truncated multivariate distribution, *J. R. Stat. Soc. Series B (Methodological)*, **23**, 1, 223–229.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, SIAM.

Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, *Rev. Geophys.*, **20**(2), 219, doi:10.1029/RG020i002p00219.

Villegas-Lanza, J.C. et al., 2015. A mixed seismic-aseismic stress release episode in the andean subduction zone, *Nature Geoscience*, **9**(2), 150–154.

Villegas-Lanza, J.C., Chlieh, M., Cavalié, O., Tavera, H., Baby, P., Chire-Chira, J. & Nocquet, J.M., 2016. Active tectonics of Peru: Heterogeneous interseismic coupling along the Nazca megathrust, rigid motion of the Peruvian Sliver, and Subandean shortening accommodation, *J. geophys. Res.*, **121**(10), 7371–7394.

Wilhelm, S. & Manjunath, B. G., 2015. Package 'tmvtnorm': Truncated Multivariate Normal and Student t Distribution, TMVTNORM: Truncated Multivariate Normal and Student t Distribution. Available from CRAN.R-project.org/package=tmvtnorm, last accessed May 2017.

Yabuki, T. & Matsu'ura, M., 1992. Geodetic data inversion using a bayesian information criterion for spatial distribution of fault slip, *J. Geophys. Int.*, **109**(2), 363–375.

Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of Green's function into waveform inversion for seismic source processes, *J. Geophys. Int.*, **186**(2), 711–720.

## APPENDIX A: TMVN PROBABILITIES CALCULATION AND PRACTICAL IMPLEMENTATION

The core of this paper relies on the possibility to evaluate multivariate normal integral over an *m*-dimensional axis-aligned hyper-rectangle, that is, an integral of the form

$$ \int_{\mathbb{A}} \exp \left( -\frac{1}{2}(x - \mu)^T \Sigma^{-1} (x - \mu) \right) dx, \qquad (A1) $$
where $\mu$ is the untruncated multivariate normal mean vector of length $p$ and $\Sigma$ is a symmetric positive definite covariance matrix. For the bounded case (eq. 12) $\mathbb{A} = [0; b_i]^p$ and for the non-negative case (eq. 11), the hyper-rectangle is unbounded on one side, that is, $b_i = +\infty$.

Most available routines to evaluate eq. (A1) implement methods and codes developed by Alan Genz, involving two steps. In a first step, a series of three transformations is applied to the original integral making it more suitable for subsequent numerical integration. These transformations were introduced in Genz (1992) and can be summarized as follows: (1) a Cholesky decomposition of the covariance matrix $\Sigma$ results in a new problem involving the integration of independent standard univariate normal distribution functions and their inverse, but over a new integration region where integration bounds are not constant anymore; (2) the new integration bounds are expressed as functions of the cumulative distribution function (cdf) and inverse cdf of a univariate Gaussian and (3) the problem is finally transformed to an equivalent integral over the unit hypercube $[0; 1]^p$. The obtained equivalent problem is more suitable for numerical integration because it has ordered the importance of variables for the integration. In a second step, the numerical multidimensional integration is performed over the unit hypercube. For this step, different schemes have been proposed. Both methods tested in this paper to produce the results shown in Sections 4 and 5 use a so-called randomized lattice rule, where projections of integration points onto the axis produce an equidistant division of the axis, before a Monte Carlo integration (Genz & Bretz 1999; Hothorn et al. 2001).

Several implementations are available to users for scientific programming languages. Matlab and FORTRAN codes are available at Alan Genz web page (last accessed on 2017 December 14) http://www.math.wsu.edu/faculty/genz/software/software.html. Python wrappers to Genz FORTRAN routine MVNDST.f are available either within the numerical Scipy library (scipy.stats.mvn) or statsmodels module. The most elaborated library is available for the R language with the original mvtnorm package (Genz et al. 2006). Aside a wrapper to the mvtnorm package, the tmvtnorm package (Wilhelm 2015) bundles routines for 1- and 2-D marginal pdfs calculation. Furthermore, Manjunath & Wilhelm (2009) extending previous results from Tallis (1961) show that the mean for a TMVN over a hyper-rectangle can be expressed as a linear sum of the values of the individual marginal pdf at their bounds (Manjunath & Wilhelm 2009, eq. 11). Although more complicated, a similar relationship holds for the covariance (Manjunath & Wilhelm 2009, eq. 16). This approach for calculating the first two moments (mean and covariance) of a TMVN over a axis-aligned hyper-rectangle has been implemented in the mtvtnorm function, as part of the tmvtnorm package (Wilhelm 2015).

# APPENDIX B: CANONICAL FORMS OF THE COST FUNCTION

## B1 Underdetermined problem with Gaussian priors

We want to demonstrate the following relationship:

$$2S(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) + K_u, \qquad (B1)$$

where $S(m)$, $\tilde{m}$ and $C_{\tilde{m}}$ are defined in eqs (3), (6) and (7) respectively and $K_u = (Gm_0 - d_{obs})^T(GC_m G^T + C_d)^{-1}(Gm_0 - d_{obs})$.

* Proof

$$\begin{aligned} 2S(m) &= (G_m - d_{obs})^T C_d^{-1} (G_m - d_{obs}) + (m - m_0)^T C_m^{-1} (m - m_0) \\ &= m^T G^T C_d^{-1} G_m - m^T G^T C_d^{-1} d_{obs} - d_{obs}^T C_d^{-1} G_m + d_{obs}^T C_d^{-1} d_{obs} \\ &\quad + m^T C_m^{-1} m - m^T C_m^{-1} m_0 - m_0^T C_m^{-1} m + m_0^T C_m^{-1} m_0 \\ &= m^T (G^T C_d^{-1} G + C_m^{-1}) m - m^T (G^T C_d^{-1} d_{obs} + C_m^{-1} m_0) \\ &\quad - (G^T C_d^{-1} d_{obs} + C_m^{-1} m_0)^T m + d_{obs}^T C_d^{-1} d_{obs} + m_0^T C_m^{-1} m_0 \\ &= m^T C_{\tilde{m}}^{-1} m - m^T a - a^T m + d_{obs}^T C_d^{-1} d_{obs} + m_0^T C_m^{-1} m_0 \end{aligned}$$

with $a = G^T C_d^{-1} d_{obs} + C_m^{-1} m_0$. Noting that $m - \tilde{m} = C_{\tilde{m}}(C_{\tilde{m}}^{-1} m - a)$, we have

$$\begin{aligned} (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) &= (m^T C_{\tilde{m}}^{-1} - a^T) C_{\tilde{m}} (C_{\tilde{m}}^{-1} m - a) \\ &= m^T C_{\tilde{m}}^{-1} m - m^T a - a^T m + a^T C_{\tilde{m}} a. \end{aligned}$$

Thus,

$$2S(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m})^T + K$$

(Tarantola 2005, p. 66), with $K = -a^T C_{\tilde{m}} a + d_{obs}^T C_d^{-1} d_{obs} + m_0^T C_m^{-1} m_0$. K is a constant independent from $m$. We will show that $K = (Gm_0 - d_{obs})^T(GC_m G^T + C_d)^{-1}(Gm_0 - d_{obs})$

$$\begin{aligned} K &= -a^T C_{\tilde{m}} a + d_{obs}^T C_d^{-1} d_{obs} + m_0^T C_m^{-1} m_0 \\ &= -d_{obs}^T C_d^{-1} G C_{\tilde{m}} G^T C_d^{-1} d_{obs} - d_{obs}^T C_d^{-1} G C_{\tilde{m}} C_m^{-1} m_0 - m_0^T C_m^{-1} C_{\tilde{m}} G^T C_d^{-1} d_{obs} - m_0^T C_m^{-1} C_{\tilde{m}} C_m^{-1} m_0 \\ &\quad + d_{obs}^T C_d^{-1} d_{obs} + m_0^T C_m^{-1} m_0 \\ &= d_{obs}^T (C_d^{-1} - C_d^{-1} G C_{\tilde{m}} G^T C_d^{-1}) d_{obs} + m_0^T (C_m^{-1} - C_m^{-1} C_{\tilde{m}} C_m^{-1}) m_0 \\ &\quad - d_{obs}^T C_d^{-1} G C_{\tilde{m}} C_m^{-1} m_0 - (d_{obs}^T C_d^{-1} G C_{\tilde{m}} C_m^{-1} m_0)^T. \end{aligned}$$
We first show that $(C_d^{-1} - C_d^{-1}GC_{\tilde{m}}G^T C_d^{-1}) = (GC_m G^T + C_d)^{-1}$.

$$ 
\begin{align*}
(C_d^{-1} - C_d^{-1}GC_{\tilde{m}}G^T C_d^{-1})(GC_m G^T + C_d) &= C_d^{-1}GC_m G^T + I - C_d^{-1}GC_{\tilde{m}}G^T C_d^{-1}GC_m G^T - C_d^{-1}GC_{\tilde{m}}G^T \\
&= I + C_d^{-1}GC_{\tilde{m}}(C_{\tilde{m}}^{-1}C_m - G^T C_d^{-1}GC_m - I)G^T \\
&= I + C_d^{-1}GC_{\tilde{m}}((G^T C_d^{-1}G + C_m^{-1})C_m - G^T C_d^{-1}GC_m - I)G^T \\
&= I
\end{align*}
$$

We now show that $C_m^{-1} - C_m^{-1}C_{\tilde{m}}C_m^{-1} = G^T(GC_m G^T + C_d)^{-1}G$

$$ 
\begin{align*}
C_m^{-1} - C_m^{-1}C_{\tilde{m}}C_m^{-1} &= C_m^{-1} - C_m^{-1} + G^T(GC_m G^T + C_d)^{-1}GC_m C_m^{-1} \\
&= G^T(GC_m G^T + C_d)^{-1}G
\end{align*}
$$

Finally, I show that $d_{\text{obs}}^T C_d^{-1}GC_{\tilde{m}}C_m^{-1}m_0 = d_{\text{obs}}^T(GC_m G^T + C_d)^{-1}m_0$

$$ 
\begin{align*}
d_{\text{obs}}^T C_d^{-1}GC_{\tilde{m}}C_m^{-1}m_0 &= d_{\text{obs}}^T(C_d^{-1}G - C_d^{-1}GC_m G^T(GC_m G^T + C_d)^{-1}G)m_0 \\
&= d_{\text{obs}}^T(C_d^{-1} - C_d^{-1}GC_m G^T(GC_m G^T + C_d)^{-1})Gm_0 \\
&= d_{\text{obs}}^T(C_d^{-1}(GC_m G^T + C_d) - C_d^{-1}GC_m G^T)(GC_m G^T + C_d)^{-1}Gm_0 \\
&= d_{\text{obs}}^T(GC_m G^T + C_d)^{-1}Gm_0,
\end{align*}
$$

so we have

$$ 
\begin{align*}
K &= d_{\text{obs}}^T(GC_m G^T + C_d)^{-1}d_{\text{obs}} + (Gm_0)^T(GC_m G^T + C_d)^{-1}(Gm_0) \\
&\quad - d_{\text{obs}}^T(GC_m G^T + C_d)^{-1}Gm_0 - (d_{\text{obs}}^T(GC_m G^T + C_d)^{-1}Gm_0)^T \\
&= (Gm_0 - d_{\text{obs}})^T(GC_m G^T + C_d)^{-1}(Gm_0 - d_{\text{obs}}) \\
&= K_u.
\end{align*}
$$

## B2 Overdetermined problem without priors

In the case of an overdetermined problem, without prior, the cost function simplifies to

$$ 2S(m) = (G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}). $$

We want to show that in this specific case

$$ 2S(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) + K_o \quad (B2) $$

with $K_o = d_{\text{obs}}^T C_d^{-1}(d_{\text{obs}} - \tilde{d})$ and where $C_{\tilde{m}} = (G^T C_d^{-1} G)^{-1}$, $\tilde{m} = (G^T C_d^{-1} G)^{-1} G^T C_d^{-1} d_{\text{obs}} = C_{\tilde{m}} G^T C_d^{-1} d_{\text{obs}}$ and $\tilde{d} = G \tilde{m}$.

Using the same derivation as in Appendix B1,

$$ 
\begin{align*}
(G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}) &= m^T G^T C_d^{-1} G_m - m^T G^T C_d^{-1} d_{\text{obs}} - d_{\text{obs}}^T C_d^{-1} G_m + d_{\text{obs}}^T C_d^{-1} d_{\text{obs}} \\
&= m^T C_{\tilde{m}}^{-1} m - m^T a - a^T m + d_{\text{obs}}^T C_d^{-1} d_{\text{obs}}
\end{align*}
$$

with $a$ now being $a = G^T C_d^{-1} d_{\text{obs}}$.

We still have $m - \tilde{m} = C_{\tilde{m}}(C_{\tilde{m}}^{-1}m - a)$ and

$$ 
\begin{align*}
(m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) &= (m^T C_{\tilde{m}}^{-1} - a^T) C_{\tilde{m}} (C_{\tilde{m}}^{-1} m - a) \\
&= m^T C_{\tilde{m}}^{-1} m - m^T a - a^T m + a^T C_{\tilde{m}} a,
\end{align*}
$$

so

$$ 2S(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) + d_{\text{obs}}^T C_d^{-1} d_{\text{obs}} - a^T C_{\tilde{m}} a = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m}) + d_{\text{obs}}^T C_d^{-1} (d_{\text{obs}} - \tilde{d}). $$

## APPENDIX C: PARTITIONING OF THE QUADRATIC FUNCTION

For the calculation of the marginal pdfs, I seek a partitioning of the quadratic function $Q(m) = (m - \tilde{m})^T C_{\tilde{m}}^{-1} (m - \tilde{m})$ so that I can evaluate the integral for a subset of components of $m$.

I define the following partitions for vectors and matrices:

$$ m = \begin{bmatrix} m_1 \\ m_2 \end{bmatrix}, \quad \tilde{m} = \begin{bmatrix} \tilde{m}_1 \\ \tilde{m}_2 \end{bmatrix}, \quad C_{\tilde{m}} = \begin{bmatrix} C_{\tilde{m}11} & C_{\tilde{m}12} \\ C_{\tilde{m}21} & C_{\tilde{m}22} \end{bmatrix} \quad \text{and} \quad C_{\tilde{m}}^{-1} = \begin{bmatrix} \Sigma_{\tilde{m}11} & \Sigma_{\tilde{m}12} \\ \Sigma_{\tilde{m}21} & \Sigma_{\tilde{m}22} \end{bmatrix}. $$

Because $C_{\tilde{m}}$ is symmetric, $C_{\tilde{m}21} = C_{\tilde{m}12}^T$ and $\Sigma_{\tilde{m}21} = \Sigma_{\tilde{m}12}^T$.
With these partitions, we have $Q(m) = Q(m_1, m_2)$ and

$$Q(m_1, m_2) = [(m_1 - \tilde{m}_1)^T, (m_2 - \tilde{m}_2)^T] \begin{bmatrix} \Sigma_{\tilde{m}_{11}} & \Sigma_{\tilde{m}_{21}} \\ \Sigma_{\tilde{m}_{21}} & \Sigma_{\tilde{m}_{22}} \end{bmatrix} \begin{bmatrix} m_1 - \tilde{m}_1 \\ m_2 - \tilde{m}_2 \end{bmatrix} \\ = (m_1 - \tilde{m}_1)^T \Sigma_{\tilde{m}_{11}} (m_1 - \tilde{m}_1) + 2(m_1 - \tilde{m}_1)^T \Sigma_{\tilde{m}_{12}} (m_2 - \tilde{m}_2) + (m_2 - \tilde{m}_2)^T \Sigma_{\tilde{m}_{22}} (m_2 - \tilde{m}_2).$$

For the inverse of a partitioned symmetric matrix, we have the following identities:

$$\Sigma_{\tilde{m}_{11}} = (C_{\tilde{m}_{11}} - C_{\tilde{m}_{12}} C_{\tilde{m}_{22}} C_{\tilde{m}_{12}}^{-1})^{-1} = C_{\tilde{m}_{11}}^{-1} + C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}} (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1} C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1}$$

$$\Sigma_{\tilde{m}_{22}} = (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1} = C_{\tilde{m}_{22}}^{-1} + C_{\tilde{m}_{22}}^{-1} C_{\tilde{m}_{12}}^{T} (C_{\tilde{m}_{11}} - C_{\tilde{m}_{12}} C_{\tilde{m}_{22}}^{-1} C_{\tilde{m}_{12}}^{T})^{-1} C_{\tilde{m}_{12}} C_{\tilde{m}_{22}}^{-1}$$

$$\Sigma_{\tilde{m}_{12}} = (\Sigma_{\tilde{m}_{21}})^T = -C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}} (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1}$$

$$\Sigma_{\tilde{m}_{21}} = (\Sigma_{\tilde{m}_{12}})^T = -C_{\tilde{m}_{22}}^{-1} C_{\tilde{m}_{12}}^{T} (C_{\tilde{m}_{11}} - C_{\tilde{m}_{12}} C_{\tilde{m}_{22}}^{-1} C_{\tilde{m}_{12}}^{T})^{-1}$$

Substituting the second expression for $\Sigma_{\tilde{m}_{11}}$, the first expression for $\Sigma_{\tilde{m}_{22}}$ and the expression of $\Sigma_{\tilde{m}_{12}}$ into $Q(m_1, m_2)$, we get

$$Q(m_1, m_2) = (m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1) \\ + (m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}} (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1} C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1) \\ - 2(m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}} (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1} (m_2 - \tilde{m}_2) \\ + (m_2 - \tilde{m}_2)^T (C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}})^{-1} (m_2 - \tilde{m}_2).$$

Defining $A = C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}}$ and rewriting $(m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}}$ as $[C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)]^T$, we recognize the developed form of a quadratic function

$$Q(m_1, m_2) = (m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1) \\ + [C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)]^T A^{-1} [C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)] \\ - 2[C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)]^T A^{-1} (m_2 - \tilde{m}_2) \\ + (m_2 - \tilde{m}_2)^T A^{-1} (m_2 - \tilde{m}_2) \\ = (m_1 - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1) \\ + [(m_2 - \tilde{m}_2) - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)]^T A^{-1} [(m_2 - \tilde{m}_2) - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (m_1 - \tilde{m}_1)] \\ = Q_1(m_1) + Q_2(m_1, m_2),$$

where

$$Q_1(x) = (x - \tilde{m}_1)^T C_{\tilde{m}_{11}}^{-1} (x - \tilde{m}_1) \\ Q_2(x_1, x_2) = (x_2 - b(x_1))^T A^{-1} (x_2 - b(x_1)) \\ b(x) = \tilde{m}_2 + C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} (x - \tilde{m}_1) \\ A = C_{\tilde{m}_{22}} - C_{\tilde{m}_{12}}^{T} C_{\tilde{m}_{11}}^{-1} C_{\tilde{m}_{12}}.$$

## APPENDIX D: MIXED GAUSSIAN AND TRUNCATED GAUSSIAN PROBLEMS

So far, the results presented in this paper were for the bounded or non-negative problem assuming that slip occurs only along an *a priori* fixed rake (although this rake may be different for each subfault). Earthquake coseismic slips sometimes show departures from the main rake direction. In this case, for each of the *p* subfaults, the slip can be decomposed into its component along the main rake and its perpendicular components. Noting $m_\mu$ and $m_\perp$ the vector of slip along the main rake and perpendicular to it respectively, the vector of unknowns is now $m = [m_\mu^T, m_\perp^T]^T$ with 2*p* components.

For $m_\perp$, we can use a prior either as a TMVN, as a uniform or as an unbounded multivariate Gaussian, all being centred on $m_{0\perp} = 0^p$. For the case of TMVN priors, all previous results for the bounded case remain valid, with the only change that the bounds for the *ith* component of $m_\perp$ will now be $[-b_i/2, +b_i/2]$. The case of uniform priors will also be approximated using the TMVN prior with large values of $\sigma_i/b_i$ as in Section 3. The case of $m_\perp$ with unbounded multivariate Gaussian priors can also be approximated by choosing large bound values $b_i$ for the components of $m_\perp$, that is, $\sigma_i/b_i$ should be small (e.g. bounds of $[-100, 100]$ m and $\sigma_i = 1$ m). For marginal pdf, the remaining integral term in eqs (11) and (12) can be evaluated with infinite bounds that are compatible with the method from Genz & Bretz (2009) method.

Although these approximations are valid, it is interesting to scrutinize how some of the equations previously derived behave for the mixed Gaussian/truncated Gaussian case. In the following, I will assume that the prior pdf for $m_\mu$ is a TMVN with covariance $C_{m_\mu}$ and a priori $m_{0\mu}$, while the prior pdf for $m_\perp$ is an untruncated centred Gaussian, with covariance $C_{m_\perp}$ and mean $m_{0\perp} = 0^p$. I will also make the reasonable assumption that the priors for $m_\mu$ and $m_\perp$ are independent (Minson et al. 2013; Minson et al. 2014). I note $M_*$ the associated space model with dim($M_*) = 2p$. $M_\parallel$ and $M_\perp$ are the subspaces of $M_*$ corresponding to the main rake and main rake perpendicular slip components so that $M_* = M_\parallel \times M_\perp$. $M_\parallel$ can be $[0, b_i]^p$ for the bounded case or $(\mathbb{R}^+)^p$ for the non-negative case as before and $M_\perp = \mathbb{R}^p$.
## D1 Marginal probability density functions for the mixed Gaussian/truncated Gaussian case

Under these assumptions, we can write the prior pdf of $m$ as

$$ \rho_{M_*}(m) = K_{m_*} \exp\left(-\frac{1}{2}(m_{\parallel} - m_{0_{\parallel}})^T C_{m_{\parallel}}^{-1}(m_{\parallel} - m_{0_{\parallel}})\right) \exp\left(-\frac{1}{2}(m_{\perp} - m_{0_{\perp}})^T C_{m_{\perp}}^{-1}(m_{\perp} - m_{0_{\perp}})\right). \quad (D1) $$

The posterior joint pdf from eq. (8) is

$$ \sigma_{M_*}(m) = K_* \exp\left(-\frac{1}{2}(m - \tilde{m})^T C_{\tilde{m}}^{-1}(m - \tilde{m})\right). \quad (D2) $$

Using the partitioning $\tilde{m} = \begin{bmatrix} \tilde{m}_{\parallel} \\ \tilde{m}_{\perp} \end{bmatrix}$ and $C_{\tilde{m}}^{-1} = \begin{bmatrix} C_{\tilde{m}_{\parallel}} & C_{\tilde{m}_{\parallel\perp}} \\ C_{\tilde{m}_{\perp}}^T & C_{\tilde{m}_{\perp\perp}} \end{bmatrix}$. From eq. (11), we can write

$$ \begin{aligned} \sigma_{M_{\parallel}}(m_{\parallel}) &= K_* \exp\left(-\frac{1}{2} Q_1(m_{\parallel}) \int_{M_{\parallel b_2}} \exp\left(-\frac{1}{2} Q_2(m_{\parallel}, m_{\perp})\right) dm_{\perp} \right) \\ &= K_* \exp\left(-\frac{1}{2} (m_{\parallel} - \tilde{m}_{\parallel})^T C_{\tilde{m}_{\parallel 11}}^{-1} (m_{\parallel} - \tilde{m}_{\parallel}) \right) \end{aligned} \quad (D3) $$

$$ \int_{M_{\perp}} \exp\left(-\frac{1}{2}(m_{\perp} - b(m_{\parallel}))^T A^{-1}(m_{\perp} - b(m_{\parallel}))\right) dm_{\perp}. \quad (D4) $$

Now $M_{\perp} = R^p$ so that

$$ \int_{M_{\perp}} \exp\left(-\frac{1}{2}(m_{\perp} - b(m_{\parallel}))^T A_{\perp}^{-1}(m_{\perp} - b(m_{\parallel}))\right) dm_{\perp} = \sqrt{2\pi}^n \sqrt{|A_{\perp}|} \quad (D5) $$

with $A_{\parallel} = C_{\tilde{m}_{\perp}} - C_{\tilde{m}_{\parallel\perp}}^T C_{\tilde{m}_{\parallel}} C_{\tilde{m}_{\parallel\perp}}$ and $|A_{\parallel}|$ is the determinant of $A_{\perp}$. This result shows that the calculation of posterior marginal pdf for a subset of $\tilde{m}_{\parallel}$, therefore also including the 1-D posterior pdf, can be reduced because the contribution of $\tilde{m}_{\perp}$ is simply a constant. Thus, we are now again in the case of Section 2.3.

If we now want to calculate the marginal pdf for $\tilde{m}_{\perp}$, eq. (8) gives

$$ \sigma_{M_{\perp}}(m_{\perp}) = K_* \exp\left(-\frac{1}{2} Q_1(m_{\perp}) \int_{M_{\perp}} \exp\left(-\frac{1}{2} Q_2(m_{\perp}, m_{\parallel})\right) dm_{\parallel} \right) \quad (D6) $$

$$ = K_* \exp\left(-\frac{1}{2}(m_{\perp} - \tilde{m}_{\perp})^T C_{\tilde{m}_{\perp}}^{-1} (m_{\perp} - \tilde{m}_{\perp}) \right) \cdot \int_{M_{\parallel}} \exp\left(-\frac{1}{2}(m_{\parallel} - b(m_{\perp}))^T A_{\parallel}^{-1} (m_{\parallel} - b(m_{\perp}))\right) dm_{\parallel} \quad (D7) $$

with $A_{\parallel} = C_{\tilde{m}_{\perp}} - C_{\tilde{m}_{\parallel\perp}} C_{\tilde{m}_{\perp}} C_{\tilde{m}_{\parallel\perp}}^T$.

So the posterior of parameters with prior Gaussian are not Gaussian anymore but also are the product of a Gaussian with a skewness function.

## D2 MAP for the mixed Gaussian/truncated Gaussian case

From eq. (2), we have

$$ \sigma_{M_*}(m) = k \rho_{M_*}(m) \rho_D(G_m). \quad (D8) $$

Inserting the expression of the prior from eq. (D1) gives

$$ \sigma_{M_*}(m) = K_{m_*} \exp\left(-\frac{1}{2}(m_{\parallel} - m_{0_{\parallel}})^T C_{m_{\parallel}}^{-1}(m_{\parallel} - m_{0_{\parallel}})\right) \exp\left(-\frac{1}{2}(m_{\perp} - m_{0_{\perp}})^T C_{m_{\perp}}^{-1}(m_{\perp} - m_{0_{\perp}})\right) \exp\left(-\frac{1}{2}(G_m - d_{\text{obs}})^T C_d^{-1}(G_m - d_{\text{obs}})\right). $$

So, the cost function becomes

$$ 2S(m) = (G_m - d_{\text{obs}})^T C_d^{-1} (G_m - d_{\text{obs}}) + (m_{\parallel} - m_{0_{\parallel}})^T C_{m_{\parallel}}^{-1} (m_{\parallel} - m_{0_{\parallel}}) + (m_{\perp} - m_{0_{\perp}})^T C_{m_{\perp}}^{-1} (m_{\perp} - m_{0_{\perp}}). \quad (D9) $$

Defining $G = [G_{\perp}, G_{\parallel}]$ with $G_{\perp}$ and $G_{\parallel}$ being the Green's function relating respectively the unit slip in the main rake and its perpendicular to the observation vector $d_{\text{obs}}$ and

$$ C_m^{-1/2} = \begin{bmatrix} C_{m_{\parallel}}^{-1/2} & 0 \\ 0 & C_{m_{\perp}}^{-1/2} \end{bmatrix}. $$
I write

$$A = \begin{bmatrix} C_d^{-1/2} G_{\parallel} & C_d^{-1/2} G_{\perp} \\ C_{m_{\parallel}}^{-1/2} & 0 \\ 0 & C_{m_{\perp}}^{-1/2} \end{bmatrix} = \begin{bmatrix} C_d^{-1/2} G \\ C_m^{-1/2} \end{bmatrix}$$

and

$$B = \begin{bmatrix} C_d^{-1/2} d_{\text{obs}} \\ C_{m_{\parallel}}^{-1/2} m_{0\parallel} \\ C_{m_{\perp}}^{-1/2} m_{0\perp} \end{bmatrix} = \begin{bmatrix} C_d^{-1/2} d_{\text{obs}} \\ C_m^{-1/2} m_0 \end{bmatrix}.$$

We have $2S(m) = (A_m - B)^T(A_m - B)$ which can be used as input for non-negative or bounded least-squares minimization. The algorithms from Lawson & Hanson (1974) or Stark & Parker (1995) do not explicitly include an option for infinite bounds. This is not a problem since the solution is achieved by finding the subset of $m_i$ which are set to their bounds, the other $m_i$ being the results of a classical least-squares solution (Lawson & Hanson 1974). So, for parameters with untruncated Gaussian priors, setting sufficiently large bounds that will never be reached provides the correct solution for the MAP.