Specification of additional information for solving stochastic inverse problems

TL;DR

This paper explores how to identify the probability distribution of a random vector in stochastic inverse problems by incorporating additional information like moments or distribution families, using methods such as Bayes' theorem and maximum entropy.

Contribution

It introduces a framework for specifying additional information needed to recover the true law of a random vector in stochastic inverse problems, extending existing methods.

Findings

Methods can incorporate moments or distribution families as additional information.

Proposed approaches ensure solutions are consistent with known information.

Framework improves the ability to predict unobserved quantities of interest.

Abstract

Methods have been developed to identify the probability distribution of a random vector $Z$ from information consisting of its bounded range and the probability density function or moments of a quantity of interest, $Q(Z)$. The mapping from $Z$ to $Q(Z)$ may arise from a stochastic differential equation whose coefficients depend on $Z$. This problem differs from Bayesian inverse problems as the latter is primarily driven by observation noise. We motivate this work by demonstrating that additional information on $Z$ is required to recover its true law. Our objective is to identify what additional information on $Z$ is needed and propose methods to recover the law of $Z$ under such information. These methods employ tools such as Bayes' theorem, principle of maximum entropy, and forward uncertainty quantification to obtain solutions to the inverse problem that are consistent with information on $Z$ and $Q(Z)$. The additional information on $Z$ may include its moments or its family of distributions. We justify our objective by considering the capabilities of solutions to this inverse problem to predict the probability law of unobserved quantities of interest.

Equations

$$P_Z(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}{\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}{\rho }_{X_{\mathcal{L}}}(x_{\mathcal{L}})d{\mu }_{X_{\mathcal{L}}}$$

$${\mu }_Z(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}d{\mu }_{X_{\mathcal{L}}}$$

$${\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}}):={\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})={\left({\int }_{{\pi }^{-1}(x_{\mathcal{L}})}d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}\right)}^{-1}$$

$$P_Z(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}{f}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})d{x}_{\mathcal{C}}{f}_{X_{\mathcal{L}}}(x_{\mathcal{L}})d{x}_{\mathcal{L}}$$

$$f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}})={\left({\int }_{{\pi }^{-1}(x_{\mathcal{L}})}d{x}_{\mathcal{C}}\right)}^{-1}.$$

$$P(Q\le q)={\int }_0^{1}P\left(Z_2\le \frac{q}{z_1}\right)f_{Z_1}(z_1)d{z}_1={\int }_0^{q}d{z}_1+{\int }_q^1\frac{q}{z_1}d{z}_1=q-q\log (q)$$

$$f_{X_{\mathcal{L}}}(x_{\mathcal{L}})=-x_{\mathcal{L}}\cdot \log (x_{\mathcal{L}}^2/2).$$

$$f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}})={\left({\int }_{x_{\mathcal{L}}^2/2}^1\sqrt{1+\frac{x_{\mathcal{L}}^4}{4z_1^4}}d{z}_1\right)}^{-1}.$$

$$P_Z^{Unif}(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}{f}_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Unif}(x_{\mathcal{C}}|x_{\mathcal{L}})f_{X_{\mathcal{L}}}^{Unif}(x_{\mathcal{L}})d{x}_{\mathcal{C}}d{x}_{\mathcal{L}}$$

$$d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}=f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Unif}(x_{\mathcal{C}}|x_{\mathcal{L}})d{x}_{\mathcal{C}}$$

$$P_Z^{Beta}(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}{f}_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Beta}(x_{\mathcal{C}}|x_{\mathcal{L}})f_{X_{\mathcal{L}}}^{Beta}(x_{\mathcal{L}})d{x}_{\mathcal{C}}d{x}_{\mathcal{L}}$$

$$P_Z^{ansatz}(A)={\int }_{\pi (A)}{\int }_{{\pi }^{-1}(x_{\mathcal{L}})\cap A}{\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}}){\rho }_{X_{\mathcal{L}}}(x_{\mathcal{L}})d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}d{\mu }_{X_{\mathcal{L}}}$$

$${\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}})d{\mu }_{X_{\mathcal{C}}|X_{\mathcal{L}}}=f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Beta}(x_{\mathcal{C}}|x_{\mathcal{L}})d{x}_{\mathcal{C}}$$

$${\rho }_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{ansatz}(x_{\mathcal{C}}|x_{\mathcal{L}})=\frac{f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Beta}(x_{\mathcal{C}}|x_{\mathcal{L}})}{f_{X_{\mathcal{C}}|X_{\mathcal{L}}}^{Unif}(x_{\mathcal{C}}|x_{\mathcal{L}})}$$

$${\rho }_Z^{post}(Z)={\rho }_Z^{prior}(Z)\cdot \frac{{\rho }_Q(Q(Z))}{{\rho }_Q^{Q(prior)}(Q(Z))}$$

$${\tilde{A}}^N(x,y)=\sum _{j=1}^{N}v(x,{\mathbf{y}}_j){\psi }_j(y),x\in {\mathbb{R}}^d,y\in \Gamma \subset {\mathbb{R}}^k$$

$$-\nabla \cdot (A(x,\omega )\nabla U(x,\omega ))=f(x),U(x,\omega )=0\text{on}\partial D$$

$$\text{min}\frac{1}{2}\Vert U-\hat{U}{\Vert }^2+\frac{\beta }{2}\Vert A{\Vert }^2$$

$$\hat{U}(x,\omega )={\hat{u}}_0(x)+\sum _{k=1}^{\infty }\sqrt{{\hat{\lambda }}_k}{\hat{\varphi }}_k(x)Y_k(\omega )$$

$$A^{\ast }(x,\omega )=a_0(x)+\sum _{k=1}^{\infty }{a}_k(x)Y_k(\omega ).$$

$$-\frac{d}{dx}(A(x,\omega )\cdot \frac{d}{dx}U(x,\omega ))=0,U(0,\omega )=0,U(1,\omega )={\int }_0^1\frac{1}{A(y,\omega )}dy.$$

$${\text{min}}_{\tilde{A}}{\int }_0^{1}E\left[{\left|U(x,\cdot )-{\int }_0^x\frac{1}{\tilde{A}(y,\cdot )}dy\right|}^2\right]dx={\text{min}}_{\tilde{B}}{\int }_0^{1}E\left[{\left|U(x,\cdot )-{\int }_0^x\tilde{B}(y,\cdot )dy\right|}^2\right]dx.$$

$$\tilde{B}(x,\omega )=b_0(x)+\sum _{k=1}^M{b}_k(x)Y_k(\omega )$$

$$\sqrt{{\lambda }_k}{\varphi }_k(x)={\int }_0^x{b}_k(y)dy\text{for}k=0,\dots ,M$$

$$\tilde{A}(x,\omega )=a_0(x)+\sum _{k=1}^M{a}_k(x)Y_k(\omega )$$

$${\text{min}}_{\tilde{A}}{\int }_0^{1}E[|A(x,\cdot )-\tilde{A}(x,\cdot ){|}^2]dx$$

$$\frac{d}{dx}U(x,\omega )=A(x,\omega ),U(0,\omega )=0$$

$$A(x,\omega )=\sum _{k=1}^{\infty }\sqrt{{\lambda }_k^G}{\varphi }_k^G(x)Y_k(\omega ),U(x,\omega )=\sum _{k=1}^{\infty }\sqrt{{\lambda }_k^U}{\varphi }_k^U(x){\tilde{Y}}_k(\omega )$$

$$\underset{f_Z}{\mathrm{minimize}}$$

$${\int }_{\Gamma }{g}_k(z)f_Z(z)dz={\mu }_k,k=1,\dots ,N,$$

$${\int }_{\Gamma }{f}_Z(z)dz=1.$$

Full text

1 Introduction

Inverse problems emerge from applications in science and engineering when information about inputs to a system is sought given measurements of observable quantities. Suppose that the physical system is modeled by the mapping $A(x,Z)\mapsto U(x,Z)$ where $x\in\mathbb{R}^{d}$ is the spatial variable, $Z\in\mathbb{R}^{n}$ is a parameter, $A(x,Z)$ is a deterministic function, and $U(x,Z)$ is the response. One of the most commonly investigated aspects of this mapping involves deterministic inverse problems. They deal with the estimation of the unknown parameter $Z$ provided measurements $\{U_{x_{i}}\}_{i=1}^{N_{obs}}$ of $U$ at spatial points $x_{i}\in\mathbb{R}^{d}$ which are referred to as quantities of interest. The ill-posedness of this problem is due to the non-uniqueness of the solution for $Z$ and is commonly addressed via two well-established methods. Optimization approaches [4] solve for $Z$ by minimizing the objective function $\sum_{i=1}^{N_{obs}}|U_{x_{i}}-U(x_{i},Z)|^{2}+\lambda^{2}\|Z\|^{2}$ where the regularization term $\lambda^{2}\|Z\|^{2}$ suppresses noisy solutions. In contrast, Bayesian approaches [5] construct the solution as a probability density function (pdf) for $Z$ instead of a point estimate, i.e. we acquire a probabilistic solution to a deterministic problem, by specifying a prior pdf on $Z$ based on available information on the unknown parameter. The observations $\{U_{x_{i}}\}_{i=1}^{N_{obs}}$ are typically assumed to be polluted by random noise whose law coupled with the prior on $Z$ yield the posterior pdf on $Z$. Despite the connections between both approaches, they emphasize that additional information on $Z$ is required to address the ill-posedness that arises from solving the inverse problem.

In this work, we focus instead on a different class of inverse problems where the unknown quantity is inherently stochastic. Using the same mapping as above, let $Z$ be a random vector defined on the probability space $(\Omega,\mathcal{F},P)$ such that $A(x,Z)$ is a random field while $U(x,Z)$ is the stochastic response. The stochastic inverse problem we address is the following [9, 10, 11, 22]:

Determine the probability law of $Z$ given probabilistic information (such as pdf) of the quantity of interest $Q(Z)\in\mathbb{R}^{m}$ represented by functionals of the stochastic response $U(x,Z)$.

By definition, the quantity of interest $Q$ is a function of the response $U$. Examples include $Q(Z)=\text{max}_{x}|U(x,Z)|$ and $Q(Z)=(U(x_{1},Z),\dots,U(x_{m},Z))$ where $x_{1},\dots,x_{m}\in\mathbb{R}^{d}$ are fixed. Also, the output $Q(Z)$ may or may not contain observation noise; regardless, it is still a random quantity as it depends on $Z$ which distinguishes it from what is encountered in the Bayesian formulation above in which $Q(Z)$, without the additive observation noise, is deterministic. This inverse problem is the direct reverse of forward uncertainty propagation which is concerned with obtaining probabilistic information of $Q(Z)$ given the probability law of $Z$.

Without additional information on $Z$, solving this stochastic inverse problem becomes challenging since distinct probability laws on $Z$ can produce the same law for $Q(Z)$ cf. [9, pp. 1839-1840]. This ill-posedness is then compensated by imposing assumptions on the law of $Z$ [9, 10, 11, 15, 18, 19, 22]. We mainly scrutinize two methods proposed in literature in which the only information known about $Z$ comprises, at most, its bounded range. The first approach [9, 10, 11] uses the disintegration theorem for probability measures to obtain a pdf for $Z$ given the pdf of $Q(Z)$. The second approach [15] aims to approximate the unknown random field $A(x,\omega)$, $\omega\in\Omega$, by solving an optimization problem; it is shown that this is conceptually identical to the stochastic inverse problem we consider above.

These existing methodologies are examined and the inadequacy of their basic implementations in recovering the true law of $Z$ in the absence of further information on $Z$ is used to motivate the following contribution of this work:

Identification of additional information to recover the true law of $Z$.

This work offers supplemental analysis to [9, 10, 11, 23] in that while the same inverse problem is tackled, the desired properties that we seek in designing the solution are different. Our focus here is on the reconstruction of the true law which is made possible by incorporating various types of additional information on $Z$ not considered in [9, 10, 11, 23]. We investigate what type of information $Z$ needs to be equipped with and the corresponding tools that can be employed to solve the inverse problem consistently provided such additional information. Furthermore, we underscore the importance of recovering the true law of $Z$ so that the resulting distribution on $Z$ can be used to characterize unobserved quantities of interest other than those to which it was calibrated. As such, the contributions of this work aid in improving predictions based on solutions to stochastic inverse problems. The examples analyzed here are not meant to be critiques of existing literature. Instead, they offer further insight to better understand how these methods operate and suggest clues on issues that need to be accommodated to enhance their applicability.

2 Absence of information on the unknown random quantity

We survey existing methods based on: the disintegration theorem (Section 2.1) and an optimization approach (Section 2.2), to tackle the inverse problem. The information available on the unknown random quantity is limited to its bounded domain. It is argued that this is insufficient to recover the true law of the unknown quantity.

2.1 Disintegration of probability measures on generalized contours

The disintegration theorem is summarized in Section 2.1.1 which serves as the basis of the works [9, 10, 11]. Two design approaches for the disintegration theorem are then elaborated in Sections 2.1.2 and 2.1.3 with each containing an example applying the method. The examples underscore that the ansatz imposed by this methodology – i.e., the pdf on the generalized contour is uniform, may be restrictive as the pdf on the contours can possess complex behavior. They also show that the true pdf of the unknown quantity is under/overestimated if this ansatz is accepted.

2.1.1 Review of methodology

For the random vector $Z\in\mathbb{R}^{n}$ defined on the probability space $(\Omega,\mathcal{F},P)$, let $Q$ be the mapping $Q:\Gamma\rightarrow\mathcal{D}$ where $\Gamma\coloneqq Z(\Omega)$ is a compact subset, $\mathcal{D}\coloneqq Q(\Gamma)\subset\mathbb{R}^{m}$ with $m<n$, and $Q(z)=q(U(x,z))$ for $z\in\Gamma$ and for some function $q(\cdot)$ that is locally differentiable. In the rest of this work, distinction is made between pdfs constructed with respect to Lebesgue and non-Lebesgue measures. For measurable $A\subset\Gamma$ and $B\subset\mathcal{D}$, denote by $P_{Z}(A)=\displaystyle\int_{A}\rho_{Z}(z)\,d\mu_{Z}$ and $P_{Q}(B)=\displaystyle\int_{B}\rho_{Q}(q)\,d\mu_{Q}$ the probability measures on $Z$ and $Q$, respectively, where $\rho_{Z}(z)$ and $\rho_{Q}(q)$ are the corresponding pdfs with respect to the specified measures $\mu_{Z}$ and $\mu_{Q}$. Unlike in calculus-based probability, $\mu_{Z}$ and $\mu_{Q}$ are not restricted to be Lebesgue. For notation purposes, if the pdfs are constructed with respect to the Lebesgue measure, they are denoted by $f$ instead of $\rho$.

The method proposed in [9, 10, 11] then addresses the following:

Given the pdf of $Q(Z)$ and the bounded range $\Gamma$ of $Z$, estimate the pdf of $Z$.

The cited approaches compute probabilities in $\Gamma$ by viewing the domain from a coordinate system of contours instead of the traditional Cartesian system. A summary of [10] is outlined below.

For $d\in\mathcal{D}$, define the set of points $\{z\in\Gamma|Q(z)=d\}$ to be a generalized contour. This is a generalization of contour curves when $n=2,m=1$. Generalized contours are equivalence classes with the relation $a\sim b$ if and only if $Q(a)=Q(b)$. Consequently, as $\mathcal{D}$ is the range of $Q$, the domain $\Gamma$ is a union of generalized contours. A representative element from each generalized contour can then be selected which serves as an indexing mechanism across all contours. This indexing set is a $m$-dimensional manifold that intersects each contour once and is called a transverse parameterization. In other words, there is a bijection between points in the transverse parameterization and the generalized contours. To clarify these concepts, an example is shown in Figure 1(a) below for $n=2,m=1$ and $Q(z_{1},z_{2})=z_{1}\cdot z_{2}$. The red thin curves are the contours of $Q$ while the solid blue curve marked $\mathcal{L}$ ($z_{2}=z_{1}$) and the dashed green curve marked $\mathcal{L}^{\prime}$ ($z_{2}=z_{1}^{2}$) are different transverse parameterizations. Because the transverse parameterization is not unique in general, we will only consider one parameterization in what follows and denote it by $\mathcal{L}$.

Every $z\in\Gamma$ can then be expressed under this new coordinate system. Let $\pi:\Gamma\rightarrow\mathcal{L}$ be the onto mapping such that for $x_{\mathcal{L}}\in\mathcal{L}$, $\pi^{-1}(x_{\mathcal{L}})$ is the corresponding generalized contour and that for $A\subset\Gamma$, $\pi(A)$ is the portion of the transverse parameterization $\mathcal{L}$ which intersects the generalized contours contained in $A$. We associate every $z\in\Gamma$ with $(x_{\mathcal{L}},x_{\mathcal{C}})$ where $x_{\mathcal{L}}\in\mathcal{L}$ represents the generalized contour in which $z$ resides and $x_{\mathcal{C}}$ is the coordinate along the generalized contour $\pi^{-1}(x_{\mathcal{L}})$. Figure 1(b) exhibits this change of coordinate system where for this example, $x_{\mathcal{L}}$ parametrizes the arc length of the transverse from the origin while $x_{\mathcal{C}}$ parametrizes the arc length of the contour from $z_{2}=1$. The transverse $\mathcal{L}$ is parameterized by $z_{2}=z_{1}$. The green thick curve is the contour $\pi^{-1}(x_{\mathcal{L}})$ that is indexed by $x_{\mathcal{L}}=0.4$, that is, $x_{\mathcal{L}}\in\mathcal{L}$ is 0.4 units from the origin. Meanwhile, the $x_{\mathcal{C}}$-coordinate of each of the 3 magenta circles on this contour is obtained by measuring the arc length of said contour from $z_{2}=1$ up to each of the 3 marked circles.

With this new coordinate system, the randomness in $Z$ implies that there are random vectors $X_{\mathcal{L}}$ and $X_{\mathcal{C}}$ associated with $x_{\mathcal{L}}$ and $x_{\mathcal{C}}$, respectively. In particular, for measurable $K\subset\mathcal{L}$, denote by $P_{X_{\mathcal{L}}}(K)=\displaystyle\int_{K}\rho_{X_{\mathcal{L}}}(x_{\mathcal{L}})\,d\mu_{X_{\mathcal{L}}}$ the probability measure on $X_{\mathcal{L}}$ where $\rho_{X_{\mathcal{L}}}$ is the pdf on $\mathcal{L}$ with respect to $\mu_{X_{\mathcal{L}}}$ that is specified. To solve the inverse problem, the probability of a measurable set $A\subset\Gamma$ can therefore be computed using the disintegration theorem for probability measures as follows [10, Corollary 4.1, Theorems 4.4, 4.5].

Theorem 1** (Disintegration theorem)**

Let $P_{Z}$ be the probability measure defined by the law on $Z$ and $Q$ the measurable mapping between $\Gamma$ and $\mathcal{D}$. For measurable $A\subset\Gamma$, $P_{Z}$ admits the following disintegration

[TABLE]

where $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ is the conditional pdf on the generalized contour corresponding to $X_{\mathcal{L}}$ while the measure $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ satisfies

[TABLE]

*with $\mu_{Z}$ and $\mu_{X_{\mathcal{L}}}$ specified. *

Figure 2 illustrates the region of integration in (2.1) for $A\subset\Gamma$.

In Theorem 1, we have assumed that the disintegration of $P_{Z}$ into a marginal and conditional family of measures is such that the latter two measures are absolutely continuous with respect to $\mu_{X_{\mathcal{L}}}$ and $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ to admit the pdfs given by $\rho_{X_{\mathcal{L}}}(x_{\mathcal{L}})$ and $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$. If the probability density functions $\rho_{X_{\mathcal{L}}}$ and $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ are specified, the pdf of $Z$ can then be estimated. One such computational approach is to partition $\Gamma$ into Voronoi cells and estimate (2.1) on each partition as in [11]. Due to the bijection between $\mathcal{L}$ and $\mathcal{D}$, the disintegration theorem guarantees that when the pdf of $Z$ constructed in this manner is propagated (pushed forward) through the model, the resulting pdf on $Q$ matches the prescribed $\rho_{Q}$. This is because $\rho_{X_{\mathcal{L}}}$ is fully specified given knowledge of $\rho_{Q}$ from the relation $\displaystyle\int_{K}\rho_{X_{\mathcal{L}}}(x_{\mathcal{L}})\,d\mu_{X_{\mathcal{L}}}=\int_{Q(K)}\rho_{Q}(q)\,d\mu_{Q}$ for measurable $K\subset\mathcal{L}$. On the other hand, $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ cannot be identified solely relying on information from $\rho_{Q}(q)$, rendering the inverse problem ill-posed. Two distinct choices for $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ yield distinct pdfs on $Z$ whose resulting pdf on $Q$ when propagated through the model both match $\rho_{Q}$. It has been argued that it is reasonable to assume that $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ is uniform over the generalized contour, i.e.

[TABLE]

(see [10, equation 4.5]). The ansatz (2.3) has been employed with apparently satisfactory results in examples related to recovering the true probability law of $Z$ in PDE models [10] as well as in inverse problems that arise in storm-surge applications [10], hydrodynamic models [11], and groundwater contamination [12].

To summarize, the method discussed constructs a consistent probability measure on $Z$, i.e. the push-forward of this measure matches the observed probability measure on $Q(Z)$. In other words, it constructs a particular pullback measure of an observed measure. It is therefore reasonable to use this method to attempt to “recover” the true probability measure/pdf on $Z$ assuming information on $Q(Z)$ only. Of course, this can be successfully accomplished if $m=n$ and if $Q$ is bijective through a standard change of variables. When this is not the case, the ansatz (2.3) is likely insufficient if the objective is beyond consistency.

In order to apply the method above, the measures $\mu_{Z},\mu_{Q},\mu_{X_{\mathcal{L}}}$ have to be specified. In the next 2 sections, we examine design approaches for these measures as presented in literature [10, 11]. We also investigate the issues that may arise from application of the ansatz (2.3).

2.1.2 Pdfs with respect to Lebesgue measures

Following Figures 3, 4 and Equations 4.3, 4.4 of [11], we interpret $\mu_{Z},\mu_{Q},\mu_{X_{\mathcal{L}}}$ to be chosen as the Lebesgue measure so that $d\mu_{Z}(z)=dz$, $\mu_{Q}(q)=dq$, $d\mu_{X_{\mathcal{L}}}(x_{\mathcal{L}})=dx_{\mathcal{L}}$ and the respective densities are $f_{Z},f_{Q},f_{X_{\mathcal{L}}}$. Using geometric arguments, it follows from (2.2) that $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ is also Lebesgue so that $d\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}=dx_{\mathcal{C}}$. As a consequence, the disintegration theorem (2.2) can be rewritten as

[TABLE]

while the ansatz (2.3) becomes

[TABLE]

Note the change in notation to emphasize that the pdfs are now with respect to the Lebesgue measure.

Demonstration of method on an example. Since all measures in the disintegration theorem have been specified, the method can now be applied to solve stochastic inverse problems. In the following, we consider a simple example to emphasize two points. First, assuming that the pdf on the generalized contours is uniform may not always enable the recovery of the true pdf of $Z$ for any mapping. The calculations carried out that indicate that the true pdf is under/overestimated is supported by Figure 4. Second, the pdf along the generalized contours can have substantial variation for different types of laws on $Z$. Figures 5 and 6 show how complicated the behavior can be of the pdf on the contour.

Example 1

Suppose that the inverse problem is independent of the spatial discretization $x$. Consider $Q(Z)=Z_{1}\cdot Z_{2},\,Z=(Z_{1},Z_{2})$ where $Z_{1},Z_{2}\sim U(0,1)$ and are independent. It is shown that given the pdf $f_{Q}$ of $Q$, (1) the methodology above is unable to recover the probability law of $Z$ and (2) the pdf along the generalized contours of $Q$ are not uniform.

Consider the contours of $Q$ (red thin curves) with the transverse parameterization $\mathcal{L}$ (thick blue line) where $0\leq x_{\mathcal{L}}\leq\sqrt{2}$ in Figure 1(a). In order to apply the proposed methodology, suppose that the support $\Gamma\coloneqq Z(\Omega)=[0,1]\times[0,1]$ is known and that $f_{Q}$ is given. The pdf $f_{Q}$ can be computed as follows: for $q\in(0,1]$,

[TABLE]

where $f_{Z_{1}}$ is the pdf of $Z_{1}$. Thus, $f_{Q}(q)=-\log(q)$ for $0<q\leq 1$. Given $f_{Q}$, it remains to determine $f_{X_{\mathcal{L}}}$ and $f_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ in order to compute (2.4).

As mentioned above, $f_{X_{\mathcal{L}}}$ can be uniquely obtained from $f_{Q}$. For $x_{\mathcal{L}}\in\mathcal{L}$, the generalized contour $\pi^{-1}(x_{\mathcal{L}})$ corresponding to $x_{\mathcal{L}}$ passes through the point $(z_{1},z_{2})=\left(\frac{x_{\mathcal{L}}}{\sqrt{2}},\frac{x_{\mathcal{L}}}{\sqrt{2}}\right)$. Hence, the contour $\pi^{-1}(x_{\mathcal{L}})$ can be parameterized as $z_{2}=\frac{x_{\mathcal{L}}^{2}}{2z_{1}}$, $x_{\mathcal{L}}^{2}/2\leq z_{1}\leq 1$, and we deduce the relationship $Q=\frac{X_{\mathcal{L}}^{2}}{2}$. We therefore have that for $x_{\mathcal{L}}\in(0,\sqrt{2}]$, $P(X_{\mathcal{L}}\leq x_{\mathcal{L}})=P(2Q\leq x_{\mathcal{L}}^{2})=F_{Q}(x_{\mathcal{L}}^{2}/2)$ which yields

[TABLE]

Figure 3 shows a plot of $f_{X_{\mathcal{L}}}$ (thick magenta line) over $\mathcal{L}$. Using (2.5) and the arc length formula applied to the parameterization of the contour $\pi^{-1}(x_{\mathcal{L}})$ for $x_{\mathcal{L}}\in\mathcal{L}$, we obtain for any $x_{\mathcal{C}}\in\pi^{-1}(x_{\mathcal{L}})$:

[TABLE]

With the components of (2.4) specified, $P_{Z}(A)$ can now be approximated for any measurable $A\subset\Gamma=[0,1]^{2}$. We partition $\Gamma$ into 4 measurable regions of equal area as shown in Figure 1(a), namely the northwest $(A_{1})$, northeast $(A_{2})$, southwest $(A_{3})$, and southeast $(A_{4})$ regions. If the proposed methodology is able to recover the true probability law of $Z$ given $f_{Q}$, we expect the approximations of $P_{Z}(A_{i})$ to be close to 0.25. Numerical calculations reveal that $P_{Z}(A_{i})\approx 0.2886,0.2459,0.1770,0.2886$ for $i=1,\dots,4$, respectively which shows that $P_{Z}(A_{3})$ in particular underestimates the true probability by a significant amount.

To understand the values obtained for $P_{Z}(A_{i})$, we plot each term in the integrand in (2.4), namely $\displaystyle\int_{\pi^{-1}(x_{\mathcal{L}})\cap A_{i}}f^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})\,dx_{\mathcal{C}}$ for $i=1,\dots,4$ and $f_{X_{\mathcal{L}}}(x_{\mathcal{L}})$, both as a function of $x_{\mathcal{L}}$ in Figure 4. Each subplot corresponds to each quadrant. The first term of the integrand amounts to the proportion of the generalized contour inside the quadrant and is displayed with red solid curves. The blue dashed curves meanwhile refer to $f_{x_{\mathcal{L}}}$ for values of $x_{\mathcal{L}}$ contained within the quadrant.

It is worth noting that $P_{Z}(A_{1})$ and $P_{Z}(A_{4})$ are the largest because of the following reasons. Firstly, the regions $A_{1}$ and $A_{4}$ are spanned by $0\leq x_{\mathcal{L}}\leq 1$ whereas $f_{X_{\mathcal{L}}}(x_{\mathcal{L}})$ attains its maximum in the neighborhood $0.2\leq x_{\mathcal{L}}\leq 0.8$. Secondly, these two quadrants possess more contours than regions $A_{2}$ and $A_{3}$ because every contour passing through the 2 former regions always passes through one of the latter regions. This implies that $A_{1}$ and $A_{4}$ comprise a wider range of values of $x_{\mathcal{L}}$ as evidenced by the domains of the corresponding subplots in Figure 4.

Finally, we numerically demonstrate that for this example, the actual conditional pdf along the generalized contours is not necessarily uniform (see Appendix A for the methodology). The actual conditional pdf results from disintegrating the measure $P_{Z}$ associated with the uniform distribution on $Z$. The left panel of Figure 5 shows three contours of $Q$, the middle panel shows the corresponding actual pdfs along each contour, whereas for comparison, the right panel shows the pdf according to the ansatz (2.5). As the concavity of the contours increases, the actual pdf along the contour becomes less uniform. We see that the pdf on each contour resulting from the disintegration is not necessarily uniform even though $Z$ is uniformly distributed on $\Gamma$. In addition, Figure 5 corroborates the results in Figure 4: the majority of the contours residing in $A_{3}$ exhibit high concavity which lead to the underestimation of $P_{Z}(A_{3})$. If instead the $f_{X_{\mathcal{C}}\,|\,X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ as obtained in the middle panel of Figure 5 were used in computing $P_{Z}(A)$ in (2.4), the resulting $Z$ would be uniformly distributed. Example 1 therefore underscores the point made at the beginning of this section that if the objective is to recover the true pdf on $Z$, then the ansatz (2.5), which only ensures that the constructed measure/density is a pullback measure, may be insufficient to recover the structure of the true pdf in directions not informed by the data.

Furthermore, we evaluate how (2.5) fares as a means of regularizing across other pdfs on $X_{\mathcal{C}}|X_{\mathcal{L}}$ in the absence of information on $Z$. Suppose instead that $Z_{1},Z_{2}$ are independent with $Z_{1}\sim Beta(\nu_{1},\nu_{2})$ and $Z_{2}\sim Beta(\tau_{1},\tau_{2})$. Notice that the pdf of $Z$ in Example 1 is a special case with $\nu_{1},\nu_{2},\tau_{1},\tau_{2}$ all being equal to 1. The 3 subplots of Figure 6 show the conditional pdf along the solid blue contour in Figure 5 (contour 2) for different combinations of the parameters for $Z_{1},Z_{2}$. The plots reveal that the pdfs on the contours can be very complex and suggest that using the pdf along the contour as the only means of regularizing against other solutions to this inverse problem may be insufficient.

2.1.3 Pdfs with respect to non-Lebesgue measures

In contrast to Section 2.1.2, [10] chose the measures in (2.1) as follows: $\mu_{Z}$ is taken to be Lebesgue, $\mu_{Q}$ is defined to be the pushforward measure of $\mu_{Z}$ through $Q$ (i.e. for measurable $B\subset\mathcal{D}$, $\mu_{Q}(B)=\mu_{Z}(Q^{-1}(B))$), while $\mu_{X_{\mathcal{L}}}$ is computed as $\mu_{X_{\mathcal{L}}}(K)=\mu_{Q}(Q(A))$ where $K=\pi(A)$. Note that $\mu_{Q}$ may not be Lebesgue especially if $Q$ is nonlinear which implies that the same holds for $\mu_{X_{\mathcal{L}}}$. As before, $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ results from (2.2) since $\mu_{Z},\mu_{X_{\mathcal{L}}}$ are now specified; however, if $\mu_{X_{\mathcal{L}}}$ is not be the Lebesgue measure, then neither is $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$.

To further distinguish this design approach from that of Section 2.1.2, we emphasize that for fixed $x_{\mathcal{L}}$, $\rho^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ is constant with respect to $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$. If we were to express this pdf with respect to the Lebesgue measure, i.e. we seek $f^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ such that $\rho^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})\,d\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}=f^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})\,dx_{\mathcal{C}}$, then $f^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ may not be constant unlike in (2.5). In addition, this approach offers a computationally efficient approximation for the case when $Z,Q(Z)$ are high-dimensional that does not require explicitly constructing $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ from (2.2) [10, Algorithm 1]. The publicly available code111https://github.com/UT-CHG/BET is based on this design.

Demonstration of method on an example. Even with this selection of the measures in (2.1), it is shown that the ansatz (2.3) may still be unable to recover the true law of $Z$. We show this in the next example where the true distribution of $Z$ in Example 1 is modified. The established results in Figures 5 and 6 of Section 2.1.2 will be invoked to aid in the discussion.

Example 2

We revisit the setup in Example 1 but with the true distribution of $Z$ altered to $Z_{1}\sim Beta(\nu_{1},\nu_{2}),Z_{2}\sim Beta(\tau_{1},\tau_{2})$, independent. Given the pdf $\rho_{Q}$ of $Q$, it is shown that the probability law of $Z$ cannot be recovered with the choice of ansatz (2.3).

In order to apply the disintegration theorem (2.1), $\mu_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ needs to be first identified from (2.2). Let $P^{Unif}_{Z}$ be the probability measure of the uniform distribution on $Z$ where $Z_{1},Z_{2}$ are independent. For measurable $A\subset\Gamma$, (2.4) and the results in Section 2.1.2 yield the disintegration

[TABLE]

where $f^{Unif}_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ are the pdfs displayed in the middle panel of Figure 5. Since $\mu_{Z}$ is Lebesgue, $P_{Z}^{Unif}(A)=\mu_{Z}(A)$ which implies that

[TABLE]

by comparing (2.7) with (2.2).

We now return to the setup in Example 2 where $Z_{1},Z_{2}$ are actually independent beta random variables. Denote by $P_{Z}^{Beta}$ the probability measure on $Z$ corresponding to its true law. From the results in Section 2.1.2, its actual disintegration is

[TABLE]

where $f^{Beta}_{X_{\mathcal{C}}|X_{\mathcal{L}}}$ is obtained similarly as the pdfs plotted in the panels of Figure 6. On the other hand, applying the ansatz (2.3) to (2.1), the inverse problem solution is

[TABLE]

where $\rho^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ is a constant for fixed $x_{\mathcal{L}}$. If (2.10) is able to recover the true law of $Z$, it must be that $P^{Beta}_{Z}(A)=P^{ansatz}_{Z}(A)$ for any measurable $A\subset\Gamma$. Since the probability measure on $X_{\mathcal{L}}$ is the same regardless of the choice of $\mu_{X_{\mathcal{L}}}$, it follows that $\rho_{X_{\mathcal{L}}}(x_{\mathcal{L}})\,\,d\mu_{X_{\mathcal{L}}}=f^{Beta}_{X_{\mathcal{L}}}(x_{\mathcal{L}})\,\,dx_{\mathcal{L}}.$ Hence, by comparing (2.10) with (2.9), the ansatz is able to recover the true law if and only if

[TABLE]

which is equivalent to

[TABLE]

using (2.8).

We now proceed by contradiction. Consider contour 2 (solid blue) in the left panel of Figure 5 and let $\nu_{1},\nu_{2},\tau_{1},\tau_{2}$ be any of the values utilized in the plots of Figure 6. The pdf $f^{Beta}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ on this contour could then be any of these plotted pdfs. It then follows from the middle panel of Figure 5 and each of the panels of Figure 6 that $\frac{f^{Beta}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})}{f^{Unif}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})}$ is not constant for $x_{\mathcal{L}}$ fixed which contradicts the assumption on $\rho^{ansatz}_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$. The ansatz is therefore unable to recover the true law of $Z$.

Despite the examples presented in Sections 2.1.2 and 2.1.3, the methodology in [9, 10, 11] can still be useful in physical applications in that it can serve as a first model for the unknown pdf of $Z$. This can be later tuned or scrutinized for plausibility depending on available information on $Z$.

2.1.4 Incorporating prior information

An alternative to the method in Section 2.1.1 under the same specifications on the inverse problem has been developed in [23]. It assumes that some information on $Z$ is available in the form of a prior pdf $\rho_{Z}^{prior}(Z)$. Analogous to Bayes’ theorem, the solution to the inverse problem is a posterior pdf $\rho_{Z}^{post}$ on $Z$ that is constructed as

[TABLE]

for $Z\in\Gamma$ where $\rho_{Q}(\cdot)$ is the given pdf of $Q$ while $\rho_{Q}^{Q(prior)}(\cdot)$ is the pdf of $Q$ that is obtained by propagating $\rho_{Z}^{prior}$ through the model. This solution was derived using the disintegration theorem based on conditional densities which is more general than Theorem 1 based on generalized contours. Although this method does not explicitly deal with contours, it is related to the method in Section 2.1 in that $\rho_{Z}^{prior}(Z)$ implies a pdf $\rho_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ on the generalized contours that is not necessarily uniform. We note that a sufficient condition for this method to recover the true pdf $\rho_{Z}$ on $Z$ would be if $\rho_{Z}^{prior}=\rho_{Z}$, signifying no gain in information. A more general sufficient condition only requires that the conditional pdf along the contours arising from the disintegration of $\rho_{Z}^{prior}$ and the true pdf $\rho_{Z}$ need to be equal [23, Sections 3, 7.3]. The methodology proposed in [23] only stresses the need for additional information to be specified on $Z$ in order to solve the inverse problem.

2.2 Parametric representations of the unknown random field

This section elaborates on the second approach [15] which estimates the coefficient field of a differential equation given observations of the solution field. Section 2.2.1 reviews the methodology and establishes its similarity with the above inverse problem. Section 2.2.2 meanwhile explores two examples utilizing this method. It is shown that the probability law and the truncation level of the random variables arising from the Karhunen-Loève expansion of the solution field may be inadequate to characterize the coefficient field.

2.2.1 Review of methodology

Consider the stochastic equation $\mathcal{L}(U(x,\omega))=0$ defined on the probability space $(\Omega,\mathcal{F},P)$ for $\omega\in\Omega,x\in D$ where the operator $\mathcal{L}$ characterizes a stochastic differential equation that depends on the random field $A(x,\omega)$. The inverse problem tackled by [19, 18, 15, 17] arising from this set-up is:

Given observations $\hat{U}$ of the solution $U(x,\omega)$, estimate the unknown field $A(x,\omega)$.

Although this problem appears different from the one posed in Section 1, they are in fact conceptually identical. In practice, the spatial domain is discretized so that the random fields $A$ and $U$ are represented as random vectors characterized by the finite-dimensional distributions of the random fields. The inverse problem is now reminiscent of the problem above. The use of the finite-dimensional distribution of $A(x,\omega)$ to characterize the field itself can be justified under the mild assumption that $A(x,\omega)$ has almost surely continuous sample paths or that it satisfies the Hölder continuity condition [1, Theorem 3.1]; see [1, 6] for more details. In what follows, it will be assumed that either assumption on $A$ holds.

Despite this connection, proposed methodologies [19, 18, 15, 17] estimate the unknown random field by finding a finite-dimensional noise approximation, that is, $A(x,\omega)\approx A(x,Z(\omega))$ where $Z$ is a random vector. Essentially, these methodologies construct a parametric representation of the unknown field $A$ in which the law, the dimension of $Z$, and the functional form of $A(x,Z)$ are to be determined. In [17], each observed sample of $U$ is used to acquire samples of $A$ via optimization; these samples of $A$ then serve to calibrate a polynomial chaos expansion (PCE) [1] for $A$. This optimization procedure may result in non-unique global minima yet its implications on the fitted PCE were not addressed. In [19], observed samples of $U$ are utilized to obtain a truncated PCE of $U$. The unknown $A$ is then expressed in terms of this stochastic basis. The limitations of PCE are well known [13, 14] and furthermore, as the stochastic basis for $U$ is used as the stochastic basis for the unknown $A$, the truncation level of the PCE for $U$ can be insufficient for the PCE of $A$. The following illustrations tackles some of these issues.

In order to address the limitations of a PCE model for $A$, [18, 15] proposed to express the unknown random field as a sparse grid representation $\widetilde{A}^{N}(x,Y)$ of level $N$. This is represented as

[TABLE]

for $\Gamma$ bounded. Here, $Y\in\Gamma$ is a random vector whose dimension is not necessarily identical to that of $Z$ and whose law has to be specified. In addition, $\{\mathbf{y}_{j}\}_{j=1}^{N}\subset\Gamma$ are the $N$ sparse grid nodes, $\psi_{j}(y)$ are specified interpolating functions on $\Gamma$, while $v(x,\mathbf{y}_{j})$ is a deterministic function of $x$. It only remains to address the following: How must the dimension and the law of the random vector $Y$ be specified for the nodal values $v(x,\mathbf{y}_{j})$ to be approximated through optimization?

Both [18, 15] are similar in that they employ the finite-dimensional model (2.13) for $A$ yet they differ in how the dimension and the law of $Y$ are prescribed. In [18], the dimension of $Y$ is postulated to be $k=1,2$ in the numerical examples while the choice of law for $Y$ is downplayed. This was demonstrated with a numerical example in which various distributions, with various support, were chosen for $Y\in\mathbb{R}^{1}$ to compute moments of $\widetilde{A}^{N}(x,Y)$. In contrast, [15] pursued an approach in selecting the dimension and law of $Y$ based on the Karhunen-Loève (KL) expansion [1] and is detailed as follows. Consider the elliptic system

[TABLE]

where $x\in D\subset\mathbb{R}^{d},\,\,\omega\in\Omega$. Given observations $\hat{U}$ (possibly noisy) of $U$, the objective is to approximate the unknown field $A$ by solving the optimization problem

[TABLE]

for $U$ and $A$ under the constraint that they satisfy (2.14) and that $A$ satisfies ellipticity constraints. The second term in (2.15) regularizes the solution $A$ through the parameter $\beta>0$ while the norms in (2.15) are formed using a tensor product of norms for Sobolev spaces and the $L^{2}(\Omega)$ norm for the probability space. A KL expansion of $\hat{U}$ is then performed to obtain for $\omega\in\Omega$

[TABLE]

where $E[Y_{k}]=0,\,E[Y_{k}^{2}]=1\,\,\forall k$, $E[Y_{k}Y_{j}]=0$ for $k\neq j$ and $\hat{\lambda}_{k},\,\hat{\phi}_{k}(x)$ are the eigenvalues and eigenfunctions of the covariance function of $\hat{U}(x,\omega)$. It is assumed [15, p. 10] that $\{Y_{i}\}_{i=1}^{\infty}$ form a basis for $L^{2}(\Omega)$ in the sense that every random variable with finite variance can be expressed as a linear combination of $\{Y_{i}\}_{i=1}^{\infty}$. As a consequence, since the optimal random field $A^{*}$ for $A$ in (2.15) satisfies $A^{*}(x,\cdot)\in L^{2}(\Omega)\,\,\forall x\in D$, $A^{*}$ can be written as

[TABLE]

The deterministic functions $\{a_{k}(x)\}_{k=0}^{\infty}$ are determined by solving (2.15) while $\{Y_{k}(\omega)\}_{k=1}^{\infty}$ are obtained from (2.16). Due to the dependence of $A^{*}$ on an infinite number of random variables, (2.17) is usually referred to as the solution to the infinite-dimensional problem.

For practical numerical implementation, the KL expansion in (2.16) is truncated to only consider $\{Y_{k}\}_{k=1}^{M},\,M\ll\infty$ based on the decay of $\lambda_{k}$. The optimal solution then takes on the form $A^{\dagger}(x,\omega)=A^{\dagger}(x,Y_{1}(\omega),\dots,Y_{M}(\omega))$, i.e. the optimal $A^{\dagger}$ that is sought from (2.15) is a function of $Y_{1},\dots,Y_{M}$ only. $A^{\dagger}$ is typically referred to as the solution to the finite noise problem. It was shown in [15] that the sequence of minimizers of the finite noise problem has a subsequence that converges weakly to a minimizer of the infinite-dimensional problem under the assumption (2.17). Due to this truncation, $A^{\dagger}(x,Y_{1}(\omega),\dots,Y_{M}(\omega))$ is not expected to be linear in $\{Y_{k}\}_{k=1}^{M}$ as in (2.17); as such, a sparse grid representation (2.13) for $A^{\dagger}(x,Y_{1}(\omega),\dots,Y_{M}(\omega))$ using linear hat functions for $\psi_{j}$ was considered in [15] to accommodate smoothness conditions on $A^{\dagger}$ as a function of $Y_{1},\dots,Y_{M}$.

In summary, [15] parameterized the unknown field $A(x,\omega)$ by a random vector $Y$ whose dimension and distribution are based on the random variables arising from the truncated KL expansion of the observed samples $\hat{U}$. The approach pursued by [15] in specifying the random vector $Y$ is reasonable because it is shown in [21] that if $A$ in (2.14) depends on $Y$, the response $U$ is analytical in $Y$. But challenges may surface from this approach which we demonstrate through examples in the next section. We clarify that the issue does not lie with the use of a sparse grid approximation for $A$ but with the choice of the random vector $Y$ used to construct $A^{*}(x,\omega)$ in (2.17).

2.2.2 Demonstration of method on examples

This section investigates the implications of the approach in [15]. First, unlike the PCE, the infinite set of random variables $\{Y_{k}(\omega)\}_{k=1}^{\infty}$ in (2.16) do not form a stochastic basis for $L^{2}(\Omega)$. For this reason, if $A\in L^{2}(\Omega)$, the truncated set $\{Y_{k}(\omega)\}_{k=1}^{M}$ might not even be adequate to characterize the true field $A$ because it is unlikely that $A\in L^{2}(\Omega)\cap\text{span}(\{Y_{k}\}_{k=1}^{M})$, as Example 3 will clarify. Second, even if it can be analytically shown that the random variables $\{Y_{k}\}_{k=1}^{\infty}$ in the KL expansion of $U$ also characterize $A$, the truncation level employed for the practical implementation of $U$ might not be sufficient for the implementation of $A$. Example 4 lends support to this claim.

Example 3

Consider the stochastic ODE for $x\in[0,1],\,\omega\in\Omega$:

[TABLE]

The true coefficient random field is modeled as a translation process [1] $A(x,\omega)=\alpha+(\beta-\alpha)\cdot F^{-1}_{beta}(\Phi(G(x,\omega)))$ where $\alpha=4,\beta=20$, $F^{-1}_{beta}$ is the inverse cumulative distribution function (cdf) of $Beta(1,3)$, $\Phi$ is the cdf of $N(0,1)$, and $G(x,\omega)$ is a zero mean, unit variance stationary Gaussian process with Matérn covariance, i.e. $E[G(s,\cdot)G(t,\cdot)]=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}|s-t|}{\ell}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|s-t|}{\ell}\right),\,\,s,t\in[0,1]$ where $K_{\nu}$ is the modified Bessel function of the second kind and $\nu=\frac{5}{2},\ell=0.03$. To formulate the inverse problem, we approximate $A(x,\omega)$ using observed samples of $U(x,\omega)$. It is shown that the random variables arising from the KL expansion of $U(x,\omega)$ are inadequate to characterize $A(x,\omega)$.

We generate 10000 noiseless samples of $U(x,\omega)$ by solving (2.18) for each sample of $A(x,\omega)$. The analytical solution is given by $U(x,\omega)=\displaystyle\int_{0}^{x}\frac{1}{A(y,\omega)}\,dy=\displaystyle\int_{0}^{x}B(y,\omega)\,dy$ where $B(y,\omega)\coloneqq\frac{1}{A(y,\omega)}$. A KL expansion of $U(x,\omega)$ is then performed to obtain (2.16). Figure 7 shows samples of $A(x,\omega)$ together with their respective samples $U(x,\omega)$ while Figure 8 displays histograms of the samples $Y_{k}$ in the KL expansion of $U(x,\omega)$ corresponding to the four largest eigenvalues. Note that the integral equations that need to be solved to obtain the eigenvalues and eigenfunctions of the covariance function of $U(x,\omega)$ have to be discretized solely for numerical implementation. Hence, we are essentially performing an eigendecomposition of the covariance matrix $\mathbf{K}$ of $\mathbf{U}=(U(x_{1},\omega),\dots,U(x_{M},\omega))$ where $\mathbf{K}_{ij}=Cov(U(x_{i}),U(x_{j}))$ for $\{x_{i}\}_{i=1}^{M}\subset[0,1]$. In the following, we used an extremely fine mesh with $x_{i+1}-x_{i}=0.005$ such that $M=201$. Further decreasing the mesh size does not alter the conclusions that follow.

In order to approximate the unknown field $A(x,\omega)$, we solve the optimization problem by minimizing the $L^{2}([0,1])\otimes L^{2}(\Omega)$ norm, i.e.

[TABLE]

We tackle the infinite-dimensional problem to obtain solutions of the form (2.17) instead of (2.13) because we do not perform truncation. As a first attempt, we characterize the random field $B$ (instead of $A$) in the spirit of (2.17) to yield the expression

[TABLE]

where $\{b_{k}(x)\}_{k=0}^{M}$ are to be determined. The quantity $M$ here is not the truncation level according to the decay of the eigenvalues of the covariance function of $U(x,\omega)$ but results from the discrete implementation of the KL expansion as discussed earlier.

It is observed from the optimization problem on the right-side of (2.19) that the minimum can be obtained if the unknown deterministic functions $b_{k}(x)$ satisfy

[TABLE]

where $\lambda_{k},\,\phi_{k}(x)$ are eigenvalues and eigenfunctions of $Cov(U(x),U(y))$. In other words, the minimizer is achieved if $\widetilde{B}(x,\omega)=\frac{dU(x,\omega)}{dx}$ in the mean square sense [1] which results in $Cov(U(x),U(y))=\displaystyle\int_{0}^{x}\int_{0}^{y}Cov(\widetilde{B}(s),\widetilde{B}(t))\,ds\,dt$ using second moment calculus. This equality is what we would have obtained from the analytic solution for $U$. This implies that by setting $b_{k}(x)$ as in (2.21), the approach pursued through (2.20) recovers the second-order statistics of $B$ and hence that of $A$.

We verify this numerically by noticing that solving for $\{b_{k}(x)\}_{k=1}^{M}$ in (2.19) results in a matrix least squares problem given the 10000 realizations of $U(x,\omega)$. Figure 9 shows plots of the second-order statistics of the true field $A(x,\omega)$ together with the second-order statistics of the numerical solution to the inverse problem $\widetilde{A}(x,\omega)=(\widetilde{B}(x,\omega))^{-1}$ with condition (2.21) on $b_{k}(x)$. In particular, in the left panel of Figure 9, the red dotted and blue dashed curves represent $E[\widetilde{A}(x)]$ and $E[A(x)]$ respectively, while the magenta solid and green dashed curves represent $Var[\widetilde{A}(x)]$ and $Var[A(x)]$ respectively. Only 2 curves are visible because the mean and variance of $A(x,\omega)$ and $\widetilde{A}(x,\omega)$ are almost indistinguishable. On the other hand, the right panel shows the absolute value of the difference between $c(s,t)\coloneqq Cov(A(s),A(t))$ and $\widetilde{c}(s,t)\coloneqq Cov(\widetilde{A}(s),\widetilde{A}(t))$ over $s,t\in[0,1]$ which are in good agreement.

As a second attempt, we pursue a more typical approach in which we characterize the unknown field $A$ directly instead of its inverse, consistent with [15]. To parameterize $A$, we consider

[TABLE]

where $\widetilde{A}$ solves the optimization problem (2.19) and $\{a_{k}(x)\}_{k=0}^{\infty}$ are to be determined. The unknown field is approximated by solving

[TABLE]

for $\widetilde{A}(x)$ under the norm $L^{2}([0,1])\otimes L^{2}(\Omega)$ where $A(x,\cdot)$ refers to the true field specified in Example 3. By using the mean value theorem on $U(x,\cdot)-\displaystyle\int_{0}^{x}\frac{1}{\widetilde{A}(y,\cdot)}\,dy$, the boundedness of $A(x,\cdot)$, and assumptions on the boundedness of $\widetilde{A}(x,\cdot)$, minimizing (2.23) under the constraint that $\widetilde{A}(y,\cdot)>0$ provides an upper bound for the objective function in (2.19).

Similar to the first attempt, a matrix least squares problem utilizing the 10000 samples of the true $A(x,\omega)$ was solved to optimize (2.23). Figure 10 shows, with identical legend to Figure 9, the second order statistics of $\widetilde{A}(x,\omega)$ under the second approach in comparison with that of $A(x,\omega)$. While $E[A(x)]$ and $E[\widetilde{A}(x)]$ are almost similar, there is a considerable difference in the variance and covariance functions between the 2 random fields. Since we minimized (2.23) instead of the optimization problem on the left-side of (2.19), whatever optimal solution we acquire from (2.19) of the form (2.22) cannot do better in matching the statistics of $A(x,\omega)$ than what we have already achieved. The derivations we have presented in the first attempt above do not hold anymore for the current attempt due to the nonlinearity of $U$ with respect to $A$. In addition, increasing the number of samples of $A(x,\omega)$ to 100000 does not alter the results. This underscores that the random variables $\{Y_{k}\}_{k=1}^{\infty}$ arising from the KL expansion of $U(x,\omega)$ do not necessarily form a stochastic basis for $L^{2}(\Omega)$. This implies that $\{Y_{k}\}_{k=1}^{\infty}$ may be inadequate to characterize the unknown $A(x,\omega)$, even when no truncation on the basis of decay of eigenvalues is performed. However, even when both $U(x,\omega)$ and $A(x,\omega)$ are analytically characterized by the same random variables, the KL expansion of $U$ has to be truncated for numerical implementation. In [15], the finite set $\{Y_{k}\}_{k=1}^{M}$ resulting from this truncation is then used to characterize $A$. The following example shows that the truncation level for $U$ may not be sufficient for $A$ such that the solution to the inverse problem under this approach underestimates statistics of $A$.

Example 4

Consider the stochastic ODE for $x\in[0,1],\,\,\omega\in\Omega$:

[TABLE]

whose analytical solution is given by $U(x,\omega)=\displaystyle\int_{0}^{x}A(y,\omega)\,dy$. Set $A(x,\omega)=G(x,\omega)$ where $G(x,\omega)$ is a zero mean, unit variance stationary Gaussian process with spectral density [2, p. 196] $s_{G}(\nu)=\frac{F}{(\nu^{2}-\nu_{0}^{2})^{2}+(2\zeta\nu\nu_{0}^{2})^{2}}$ in which $\nu_{0}=20$, $\zeta=0.1$, and $F$ is a scaling factor such that the correlation function of $A$, $r_{G}(\tau)$, satisfies $r_{G}(\tau)=\displaystyle\int_{-\infty}^{\infty}e^{i\nu\tau}s_{G}(\nu)\,d\nu=1.$ To formulate the inverse problem, we approximate $A(x,\omega)$ using observed samples of $U(x,\omega)$.

As $A(x,\omega)$ is a Gaussian process, the response $U(x,\omega)$ is also a Gaussian process with mean $m_{U}=\displaystyle\int_{0}^{x}E[A(y)]\,dy=0$ and correlation function $r_{U}(x,y)=\displaystyle\int_{0}^{x}\int_{0}^{y}r_{G}(s-t)\,ds\,dt$ [3]. Consequently, the KL expansion of both $U(x,\omega)$ and $A(x,\omega)$ can be expressed as

[TABLE]

where $Y_{k}(\omega),\widetilde{Y}_{k}(\omega)\sim N(0,1)\,\,\forall k$ and $\lambda^{G}_{k},\,\phi^{G}_{k}(x)$ and $\lambda^{U}_{k},\,\phi^{U}_{k}(x)$ are eigenvalues and eigenfunctions of $r_{G}$ and $r_{U}$, respectively. Thus, $A$ and $U$ are characterized by random variables with the same law unlike in Example 3. In the solution approach in [15], the truncation level of the KL expansion of $U$ is used to determine the number of random variables to characterize $A$. The truncation level $M$ of the KL expansion is usually deduced from the total variance of the random field which is given by $\displaystyle\int_{0}^{1}Var(U(x))\,dx=\sum_{k=1}^{\infty}\lambda_{k}^{U}$. The value of $M$ is then chosen to be the smallest integer such that $\frac{\sum_{k=1}^{M}\lambda_{k}^{U}}{\sum_{k=1}^{\infty}\lambda_{k}^{U}}\geq\alpha$ with $\alpha$ being close to 1.

The above principle is applied to choose $M$ using $\alpha=0.95$. To simulate solving the inverse problem, we do not generate samples of $U$ through samples of $A$ in (2.24) and estimate $r_{U}(x,y)$ from samples of $U$; rather, $r_{U}$ is obtained from the relationship between $r_{U}$ and $r_{G}$ above. Figure 11 exhibits the behavior of the eigenvalues of $r_{U}$. The left panel shows the first 30 eigenvalues $\lambda_{k}^{U}$ while the right panel displays $\frac{\sum_{k=1}^{M}\lambda_{k}^{U}}{\sum_{k=1}^{\infty}\lambda_{k}^{U}}$ as a function of $M$. It is evident that the truncation level for the KL expansion of $U$ according to the above procedure is $M=6$.

Figure 12 displays the behavior of the eigenvalues $\lambda_{k}^{G}$ of $r_{G}$ with the same legend as in Figure 11. In this case, the truncation level for the KL expansion of $A$ is $M=9$. Note that the eigenvalues of $r_{G}$ decay slower than that of $r_{U}$; intuitively, this is because large values of the frequency $\nu$ are required to capture the total energy of the spectral density $\displaystyle 2\int_{0}^{\infty}s_{G}(\nu)\,d\nu$. As $U(x,\omega)$ is obtained by integrating $A(x,\omega)$, the variation in $A(x,\omega)$ is diminished which yields a faster decay of eigenvalues for $r_{U}$.

Hence, if the truncation level of the KL expansion of $U$ is used to characterize $A$, the inverse problem solution would be $\widetilde{A}(x,\omega)=\sum_{k=1}^{6}\sqrt{\lambda_{k}^{G}}\phi_{k}^{G}(x)Y_{k}(\omega)$. Based on the discussion above, $\widetilde{A}$ may not be sufficient to capture the statistics of $A$ and we confirm this in Figure 13. The four subplots in this figure represent the first 4 moments of $\sum_{k=1}^{M}\sqrt{\lambda_{k}^{G}}\phi_{k}^{G}(x)Y_{k}(\omega)$ for $M=6,9,101$ with $M=101$ being a sufficient approximation for $M=\infty$. We notice that for $M=6$, $\widetilde{A}(x,\omega)$ underestimates the statistics of $A(x,\omega)$, especially for moments of even order. The underestimation can be avoided provided appropriate selection of the truncation level. Perhaps deducing the truncation level for $A$ using other information than the observations $U$ can resolve this issue.

To summarize, Examples 3 and 4 demonstrate that in parameterizing an unknown random field $A(x,\omega)$ by random variables $\{Y_{k}\}_{k=1}^{M}$ obtained from the response $U(x,\omega)$, the value chosen for $M$ and the probability law chosen for $Y_{k}$ can significantly affect the accuracy of the approximation for $A(x,\omega)$. A strategy to tackle the ill-posedness could instead require additional information being specified on the dimension of the random variables parameterizing $A$ and that they belong to a family of distributions subject to unknown parameters. This will be elaborated next.

3 Required additional information on the unknown random quantity

Section 2 showed that the existing methods in general cannot recover the true probability law of a random vector $Z$ if the information is limited to its bounded domain. We therefore devote this section to addressing the objective of this work: identify realistic additional information on $Z$ that is required to characterize its law. The mapping in Example 1, $Q(Z)=Z_{1}\cdot Z_{2}$, is revisited in which the true law is set as $Z_{1},Z_{2}\sim U(0,1)$, independent. All pdfs from this section onwards are constructed with respect to the Lebesgue measure and are denoted by $f$.

It is impossible to devise a general method to determine the minimum amount of additional information that is required on $Z$. Different applications possess forward models with unique properties and specific information on the unknown. Hence, we outline a few scenarios with their corresponding solution methodologies. These are categorized based on what is known about $Z$: moment information (Section 3.1) or the family of distribution to which it belongs (Section 3.2). In each category, further subcategories are considered depending on the given information on $Q(Z)$: pdf (Sections 3.1.1 and 3.2.1) or samples (Sections 3.1.2 and 3.2.2).

3.1 Information on moments of $Z$

If information about moments of the random vector $Z$ is available, the principle of maximum entropy [7] can be employed to determine the pdf of $Z$, assuming that this philosophy is accepted and that it is believed that the solution to the inverse problem resides in the subspace of maximum entropy pdfs. The principle of maximum entropy constructs the pdf $f_{Z}(z)$ of $Z$ by solving

[TABLE]

for some functions $g_{k}$, whose solution is derived as

[TABLE]

where the Lagrange multipliers $\lambda_{k}$ satisfy the relationship $\displaystyle\frac{\partial}{\partial\lambda_{k}}\int_{\Gamma}\exp[\lambda_{1}g_{1}(z)+\dots+\lambda_{N}g_{N}(z)]\,dz=\mu_{k}$ for $k=1,\dots,N$. If these multipliers exist, it can be shown that (3.2) is the unique minimizer satisfying the above constraints [7]. The succeeding sections detail how (3.2) can be used to solve the inverse problem given the pdf of $Q(Z)$ (Section 3.1.1) or samples of $Q(Z)$ (Section 3.1.2).

3.1.1 Pdf of $Q(Z)$

Assume that the pdf $f_{Q}(q)=-\log(q),\,q\in(0,1]$ of $Q$ were known. The example below illustrates the above construction.

Example 5

Suppose that the only information known about $Z$ aside from $Z\in\Gamma$ are that 1) $Z_{1}$ and $Z_{2}$ are independent and that 2) the first-order moments of $Z_{1}$ and $Z_{2}$ are within a certain range, i.e. $E[Z_{1}]=\mu_{1}\in[0,0.75],E[Z_{2}]=\mu_{2}\in[0.4,1].$ It is demonstrated how the inverse problem can be solved using an entropy-based pdf.

For a fixed value of $(\mu_{1},\mu_{2})$, the principle of maximum entropy yields the conditional pdf of $Z$ as

[TABLE]

where there is a bijective relationship between $\mu_{i}$ and $\lambda_{i}$ for $i=1,2$ via $\mu_{i}=\frac{1}{1-e^{-\lambda_{i}}}-\frac{1}{\lambda_{i}}$, $\mu_{i}\in(0,1)$. To estimate $(\mu_{1},\mu_{2})$, an optimization problem can be solved which measures the discrepancy between the given pdf $f_{Q}$ of $Q$ and the one obtained by propagating the pdf in (3.3) through the forward model $Q$ which we denote by $\widetilde{f}_{Q}(\cdot|\mu_{1},\mu_{2})$. Mathematically, this is expressed as

[TABLE]

for some distance function $d$ such as the $L^{p}$ error, Kullback-Leibler divergence, etc. Let $F_{Z_{i}}(\cdot|\mu_{1},\mu_{2})$ and $f_{Z_{i}}(\cdot|\mu_{1},\mu_{2})$ represent the marginal cdf and pdf, respectively, of $Z_{i}$ for $i=1,2$ based on (3.3). Elementary calculations similar to (2.6) yield

[TABLE]

for $q\in[0,1]$. Figure 14(a) displays the logarithm of the $L^{1}$ error $\|f_{Q}(q)-\widetilde{f}_{Q}(q|\mu_{1},\mu_{2})\|_{L^{1}}$ for $(\mu_{1},\mu_{2})$ in the specified ranges above. This discrepancy is minimized at $(\mu_{1},\mu_{2})=(0.5,0.5)$ which corresponds to $(\lambda_{1},\lambda_{2})=(0,0)$, thereby recovering the true pdf of $Z$: $f_{Z}(z_{1},z_{2}|\mu_{1},\mu_{2})=\mathbbm{1}_{(z_{1},z_{2})\in\Gamma}$ with $\mathbbm{1}$ being the indicator function.

For this simple example, the pdf of $Q$ given the pdf of $Z$ can be computed analytically. In general, however, if $f_{Z}(\cdot|\theta)$ for $\theta\in\Theta$ represents the pdf of $Z$, propagating this through more complicated forward models $Q$ implies evaluation of $Q$ multiple times. An approach to ameliorate this computational burden is to use a surrogate model [1] for $Q$ as a function of $Z$. This can be supplemented by the following procedure if $\Theta$ is low-dimensional, as done in [26]:

• •

Select $M$ points $\{\theta_{i}\}_{i=1}^{M}\subset\Theta$.

• •

For each $\theta_{i}$, $i=1,\dots,M$, propagate $f_{Z}(\cdot|\theta_{i})$ through $Q$ to approximate the pdf $\widetilde{f}_{Q}(\cdot|\theta_{i})$ of $Q$.

• •

Using $\{\widetilde{f}_{Q}(\cdot|\theta_{i})\}_{i=1}^{M}$, construct an interpolant for $\widetilde{f}_{Q}(\cdot|\theta)$ over $\Theta$.

3.1.2 Samples of $Q(Z)$

In contrast to the previous section, consider that the available information on the quantity of interest $Q$ is its $N_{s}$ samples represented by $\{q^{i}\}_{i=1}^{N_{s}}$ instead of the pdf of $Q$. Such information is what is typically encountered in practical applications of stochastic inverse problems [26, 19]. This naturally leads to employing the Bayesian framework [5] to solve the inverse problem. If information is available on moments of $Z$ in the form of a prior pdf, the principle of maximum entropy can be utilized to construct the likelihood function in Bayes’ theorem as the next example elaborates.

Example 6

We postulate that only the following information is known about $Z$: 1) $Z\in\Gamma$, 2) $Z_{1},Z_{2}$ are independent, and that 3) $(\mu_{1},\mu_{2})\coloneqq(E[Z_{1}],E[Z_{2}])$ is equally likely to take any value in the range $[0.25,0.75]^{2}$. Bayes’ theorem coupled with the principle of maximum entropy provide an approach to address the inverse problem.

Knowledge that $(\mu_{1},\mu_{2})$ is equally likely in $[0.25,0.75]^{2}$ translates to a prior pdf on $(\mu_{1},\mu_{2})$ denoted by $f_{\mu}^{prior}(\mu_{1},\mu_{2})=\mathbbm{1}_{(\mu_{1},\mu_{2})\in[0.25,0.75]^{2}}\frac{1}{0.5^{2}}$. By Bayes’ theorem, the posterior pdf on $(\mu_{1},\mu_{2})$ is

[TABLE]

in which $\ell(\mu_{1},\mu_{2}|\{q^{i}\}_{i=1}^{N_{s}})$ symbolizes the likelihood function. Given $(\mu_{1},\mu_{2})$, let $f_{Z}(\cdot|\mu_{1},\mu_{2})$ as in (3.3) be the entropy-based pdf of $Z$ while $\widetilde{f}_{Q}(\cdot|\mu_{1},\mu_{2})$ as in (3.5) be the resulting pdf when the entropy-based pdf is propagated through $Q$. The likelihood function is then established as

[TABLE]

by the independence of the samples $\{q^{i}\}_{i=1}^{N_{s}}$.

Figure 15(a) displays the posterior pdf $f_{\mu}^{post}(\cdot|\{q^{i}\}_{i=1}^{N_{s}})$ on $(\mu_{1},\mu_{2})$ using $N_{s}=100$ samples of $Q$. The maximum a posteriori (MAP) estimate is around $(\mu_{1},\mu_{2})=(0.5,0.5)$ which recovers the true pdf of $Z$. We remark that the accuracy of this approach is in part influenced by the observed number of samples $N_{s}$ of $Q$.

Does the entropy-based pdf for $Z$ guarantee that the true pdf of $Z$ can be recovered, or that the optimization problem formulated or the posterior density have a unique global minimum? The answer depends on the specific application. If not, additional information of the moments of the unknown $Z$ need to be specified to obtain a solution to the inverse problem that can be used for prediction as elaborated in Section 4.2. Sufficient conditions exist which impose criteria that the moments of $Z$ have to satisfy to uniquely determine its distribution; see [8, Theorem 3.3.11] for an example.

3.2 Parametric family of distributions of $Z$

Another strategy to combat ill-posedness is to require that the practitioner has information about the family of distributions in which the law of $Z$ resides, subject to unknown parameters $\theta$. This information is represented as $f_{Z}(\cdot|\theta)$ in which the functional form of the pdf of $Z$ is known. We illustrate how this can be employed to solve the inverse problem if the pdf of $Q(Z)$ (Section 3.2.1) or samples of $Q(Z)$ (Section 3.2.2) is given.

3.2.1 Pdf of $Q(Z)$

The next example expounds on the above idea assuming that the pdf $f_{Q}(q)=-\log(q),q\in(0,1]$ of $Q$ were known.

Example 7

Suppose that it is known that $Z_{1},Z_{2}$ are independent with $Z_{1}\sim Beta(\nu_{1},\nu_{1})$ and $Z_{2}\sim Beta(\nu_{2},\nu_{2})$ whose joint pdf is denoted by $f_{Z}(\cdot|\nu_{1},\nu_{2})$. The objective is to estimate $(\nu_{1},\nu_{2})$ such that the pdf $f_{Q}$ of $Q$ matches the pdf obtained by propagating $f_{Z}(\cdot|\nu_{1},\nu_{2})$ through $Q$.

The solution methodology for this approach is identical to the optimization procedure in (3.4) in which we write the propagated pdf of $f_{Z}(\cdot|\nu_{1},\nu_{2})$ through $Q$ as $\widetilde{f}_{Q}(\cdot|\nu_{1},\nu_{2})$. By nature of the beta distribution, $(\nu_{1},\nu_{2})\in(0,\infty)^{2}$. Figure 14(b) displays the $L^{1}$ error $\|f_{Q}(q)-\widetilde{f}_{Q}(q|\nu_{1},\nu_{2})\|_{L^{1}}$ for $(\nu_{1},\nu_{2})\in(0,10]^{2}$ where the global minimum in this domain is attained at $(\nu_{1},\nu_{2})=(1,1)$, thereby recovering the true pdf of $Z$. Heuristic arguments can be made to deduce that no other values for $(\nu_{1},\nu_{2})$ outside $(0,10]^{2}$ yield the global minimum. If $\nu_{1},\nu_{2}$ are simultaneously large, $Z_{1},Z_{2}\rightarrow\frac{1}{2}$ a.s. which implies that $Q\rightarrow\frac{1}{4}$ a.s. Likewise, if only $\nu_{1}$ is large then $Q\rightarrow\frac{1}{2}Z_{2}$ a.s. whose pdf does not match $f_{Q}(q)=-\log(q)$, and vice versa.

3.2.2 Samples of $Q(Z)$

In contrast to the previous example, suppose instead that the available information on $Q$ pertains to its $N_{s}$ samples $\{q^{i}\}_{i=1}^{N_{s}}$. Information on the family of distributions enables the construction of the likelihood function that is required to find the posterior pdf on the unknown parameters using Bayes’ theorem. The next example highlights this idea of standard Bayesian inversion.

Example 8

Consider a model in which the following information is at the practitioner’s disposal: 1) $Z_{1},Z_{2}$ are independent, 2) $Z_{1}\sim Beta(1,\nu_{1}),Z_{2}\sim Beta(1,\nu_{2})$, and 3) the prior pdf on $(\nu_{1},\nu_{2})$ is characterized by $f_{\nu}^{prior}(\nu_{1},\nu_{2})=\mathbbm{1}_{(\nu_{1},\nu_{2})\in[\frac{1}{3},3]^{2}}\frac{1}{0.5^{2}}\frac{1}{(\nu_{1}+1)^{2}}\frac{1}{(\nu_{2}+1)^{2}}$. The posterior pdf on $(\nu_{1},\nu_{2})$ results directly from Bayes’ theorem.

Since $E[Z_{i}]=\frac{1}{1+\nu_{i}}$ for $i=1,2$, the prior pdf $f_{\nu}^{prior}$ translates to a uniform prior pdf on $(E[Z_{1}],E[Z_{2}])$ with values in the range $[0.25,0.75]^{2}$. The construction of the likelihood function and the posterior pdf $f_{\nu}^{post}(\nu_{1},\nu_{2}|\{q^{i}\}_{i=1}^{N_{s}})$ is identical to that in (3.7) and (3.6), respectively, wherein the specified pdf on $Z$ conditioned on $(\nu_{1},\nu_{2})$ is propagated through $Q$ to obtain $\widetilde{f}_{Q}(\cdot|\nu_{1},\nu_{2})$. Figure 15(b) exhibits $f_{\nu}^{post}(\nu_{1},\nu_{2}|\{q^{i}\}_{i=1}^{N_{s}})$ using the same $N_{s}=100$ samples of $Q$ generated in Section 3.1.2. The MAP estimate hovers close to $(\nu_{1},\nu_{2})=(1,1)$; with more samples of $Q$, the contours of $f_{\nu}^{post}(\nu_{1},\nu_{2}|\{q^{i}\}_{i=1}^{N_{s}})$ center more at this point.

The above approach can be extended if practitioner believes that the true pdf of $Z$ belongs to multiple families of distributions, each with its own set of parameters, i.e. $f^{1}_{Z}(z|\theta^{1}),\dots,f^{N_{m}}_{Z}(z|\theta^{N_{m}})$. Bayesian model selection [28] offers a strategy to solve the inverse problem.

Since the practitioner possesses information about the family of distributions in which $Z$ resides, this approach guarantees that the true pdf of $Z$ can be recovered. Does this approach ensure that the optimization problem or the posterior density have a unique global minimum? As before, the answer is case dependent. For example, [25] considers a specific model in stochastic homogenization where it was analytically shown that their specified distribution on $Z$ and the optimization problem they formulated accommodate a unique solution on the parameters of the pdf of $Z$.

Otherwise, additional information such as moments of $Z$ need to be supplied to regularize against other plausible parameter values in the pdf of $Z$. To clarify this, assume instead that the ranges of $Z_{1},Z_{2}$ are unknown yet $f_{Q}(q)=-\log(q)$, $Z_{1}=\lambda Z_{1}^{\prime}$, $Z_{2}=\frac{1}{\lambda}Z_{2}^{\prime}$ where $\lambda>0,Z_{1}^{\prime}\sim Beta(\nu_{1},\nu_{1}),Z_{2}^{\prime}\sim Beta(\nu_{2},\nu_{2})$ are specified information. In this new model, $\lambda,\nu_{1},\nu_{2}$ are parameters to be approximated. Without additional information, the inverse problem possesses infinitely many solutions of the form $(\nu_{1},\nu_{2},\lambda)=(1,1,\lambda)$ for any $\lambda>0$. Specifying additional moment information on $Z$ such as $E[Z_{1}]$ and $E[Z_{2}]$ would resolve such ill-posedness.

4 Remarks

Section 3 dealt with the objective of this work motivated by the discussion in Section 2. As we have seen, existing methods may not succeed in recovering the true pdf on the random quantity in the absence of further information. As such, we argued on how the inverse problem should be formulated and suggested appropriate solution methods. Here, we support and clarify aspects on the model construction (Section 4.1) and on our objective (Section 4.2). In particular, the latter section highlights that in the absence of information on the unknown random quantity, the solution to the inverse problem may not be suitable to predict the law of quantities of interest other than the one it was calibrated to.

4.1 Posing the stochastic inverse problem

In the following, three remarks are made about the various formulations of the inverse problem considered in Section 3. The same notation above is used.

• •

Different models on the pdf of $Z$ result in different solutions to the inverse problem. To clarify this point, we revisit the examples presented in Sections 3.1.2 and 3.2.2. The parameters between both models are related through $\mu_{i}=\frac{1}{1+\nu_{i}}$ and the same prior pdf characterized both sets of parameters. With the same set of $N_{s}=100$ samples, we plot the posterior distributions under each model in the $(\mu_{1},\mu_{2})$ space (Figures 15(a) and 15(d)) and in the $(\nu_{1},\nu_{2})$ space (Figures 15(b) and 15(c)). It is clear that the contours manifest distinct behavior.

• •

Difference in information content between knowing the pdf of $Q(Z)$ vs only having samples of $Q(Z)$. Compared to having the pdf of $Q$ as given information, only possessing samples of $Q$ provides less information about the quantity of interest. A possible consequence of this lower information content includes posterior densities (3.6) that are not sharp about the MAP estimate. This is remedied by requiring a large number of samples of $Q(Z)$.

• •

Distinction in applying the principle of maximum entropy when moment information on $Z$ or $Q(Z)$ is supplied. Finally, we remark that the principle of maximum entropy has also been invoked to solve a stochastic inverse problem formulated differently from that of Section 3. Consider the mapping $Q:D\times\Gamma\rightarrow\mathcal{D}$ with $\Gamma\coloneqq Z(\Omega)$, $D$ being the physical domain, and $Q(x,z)=q(U(x,z))$ for $x\in D,z\in\Gamma$ and some function $q$. The methodology developed in [22] seeks to address the inverse problem described as follows:

Determine the pdf $f_{Z}$ of $Z$ given the bounded range $\Gamma$ of $Z$ and the observed moments of $Q$ up to order $N$ for $x\in D$ denoted by $\hat{\mu}^{p}(x)$, $p=1,\dots,N$.

Since the $p$-th order moments of $Q$ can be expressed as integrals on $Z$, i.e. $\displaystyle\int_{\Gamma}Q(x,z)^{p}f_{Z}(z)\,dz$ for $x\in D$, this naturally leads to a solution based on the principle of maximum entropy in which $f_{Z}$ is estimated via the optimization problem (3.1) subject to the constraint

$\displaystyle\int_{\Gamma}Q(x,z)^{p}f_{Z}(z)\,dz=\hat{\mu}^{p}(x)\,\,\,\forall x\in D,\,\,1\leq p\leq N$.

Although the principle of maximum entropy has been used as a regularizer to infer the pdf of a random vector provided information about its moments, for problems involving forward models, the absence of information on $Z$ raises issues mentioned earlier. There is no guarantee that the proposed method is able to recover the true pdf of $Z$, in contrast to the authors’ comment on p. B761. To see this, let $f^{1}_{Z},f^{2}_{Z}$ be two pdfs of $Z$ such that when propagated through the model, the pdf of $Q$ is $f_{Q}$. The $p$-th order moments of $Q$ under both pdfs of $Z$ are identical even though $f^{1}_{Z}$ has larger entropy while $f^{2}_{Z}$ is the true pdf of $Z$ or vice versa. Specifying all moments of $Q$ neither resolves the issue.

4.2 Validation

In Section 3, we considered types of information required on the unknown random quantity in order to solve the stochastic inverse problem. While some of this required information may be exigent, we argue that they are necessary to obtain solutions such that the resulting law of $Z$ can be used to characterize other quantities of interest $\widetilde{Q}$. We remark that in relation to the methods tackled in Section 2, such additional information may not be necessary if the structure of the contours of $Q$ and $\widetilde{Q}$ are similar. If this is not the case, without additional information, methods such as in Section 2.1.4 may result in a posterior pdf for $Z$ whose predicted probability measure on the new quantity of interest $\widetilde{Q}$ is similar to the predicted measure on $\widetilde{Q}$ produced by the prior. The field of optimal experimental design for prediction addresses these concerns.

Here, we revisit the solutions obtained from the methods described in Sections 2.1.2, 3.1.2, and 3.2.2. It is demonstrated that in the absence of information on $Z$, the resulting solution may be inadequate to characterize the law of quantities of interest to which it was not calibrated, thereby limiting its use in practical applications.

Example 9

We revisit the forward mapping $Q(Z_{1},Z_{2})=Z_{1}\cdot Z_{2}$ where $Z_{1},Z_{2}\sim U(0,1)$, independent, characterizes the true law on $Z$. The following methodologies are employed to solve the inverse problem on approximating the pdf of $Z$ depending on available information on $Z$ and $Q$. The resulting pdf on $Z$ is then used to predict the pdf on an unobserved quantity of interest $\widetilde{Q}(Z_{1},Z_{2})=Z_{1}+Z_{2}$. The domain $Z\in[0,1]^{2}$ is assumed for all methods.

• •

Method based on the disintegration theorem using an ansatz as in Section 2.1.2 in which $f_{Q}(q)=-\log(q)$ is given and no other information on $Z$ is required.

• •

Bayes’ theorem with entropy-based pdf for the likelihood as in Section 3.1.2 in which $N_{s}=100$ samples $\{q^{i}\}_{i=1}^{N_{s}}$ of $Q$ are available and the following is known about $Z$: $Z_{1},Z_{2}$ are independent and $f_{\mu}^{prior}(\mu_{1},\mu_{2})=\mathbbm{1}_{(\mu_{1},\mu_{2})\in[0.25,0.75]^{2}}\frac{1}{0.5^{2}}$ is the prior pdf on $(\mu_{1},\mu_{2})=(E[Z_{1}],E[Z_{2}])$.

• •

Bayes’ theorem with known family of distributions for the likelihood as in Section 3.2.2 in which the same $N_{s}=100$ samples $\{q^{i}\}_{i=1}^{N_{s}}$ of $Q$ are available as above and the following is known about $Z$: $Z_{1},Z_{2}$ are independent, $Z_{1}\sim Beta(1,\nu_{1}),Z_{2}\sim Beta(1,\nu_{2})$, and $f_{\nu}^{prior}(\nu_{1},\nu_{2})=\mathbbm{1}_{(\nu_{1},\nu_{2})\in[0.75,1.25]^{2}}\frac{1}{0.5^{2}}$ is the prior pdf on $(\nu_{1},\nu_{2})$.

The solution approach for the latter 2 methods has already been discussed. The pdf on the unobserved quantity of interest $\widetilde{Q}$ then results by computing the posterior predictive distribution. Denote by $f_{\Theta}^{post}(\theta|\{q^{i}\}_{i=1}^{N_{s}})$ the obtained posterior pdf on the corresponding parameter space which qualifies as the solution to the inverse problem upon application of either of the latter 2 methods. The pdf on $\widetilde{Q}(Z_{1},Z_{2})=Z_{1}+Z_{2}$ is obtained through

[TABLE]

where $f_{\widetilde{Q}}(\widetilde{q}|\theta)$ is the pdf on $\widetilde{Q}$ obtained by propagating the conditional pdf $f_{Z}(\cdot|\theta)$ on $Z$ through $\widetilde{Q}$. Elementary calculations show that

[TABLE]

with $f_{Z_{i}}(\cdot|\theta)$ being the marginal pdf of $Z_{i}$, $i=1,2$. In particular, under the true law on $Z$, $\widetilde{Q}$ has a triangular distribution whose pdf is $f_{\widetilde{Q}}(\widetilde{q})=\widetilde{q}$ for $\widetilde{q}\in[0,1]$ and $f_{\widetilde{Q}}(\widetilde{q})=2-\widetilde{q}$ for $\widetilde{q}\in[1,2]$.

On the other hand, the pdf $f_{Z}^{ansatz}$ arising from the method based on the disintegration theorem results through this procedure:

Partition $\Gamma=[0,1]^{2}$ into $N_{sq}=10000$ squares $\{A_{i}\}_{i=1}^{N_{sq}}$ with area $(0.01)^{2}$. 2. 2.

Compute $P_{Z}(A_{i})$ using (2.4) and (2.5) for $i=1,\dots,N_{sq}$. 3. 3.

For each $A_{i}$, let $(z_{1,i}^{*},z_{2,i}^{*})$ be its center. The approximate pdf of $Z$ is then calculated as $f_{Z}^{ansatz}(z_{1,i}^{*},z_{2,i}^{*})\simeq\frac{P_{Z}(A_{i})}{(0.01)^{2}}$.

If $(z_{1,i}^{SW},z_{2,i}^{SW})$ and $(z_{1,i}^{NE},z_{2,i}^{NE})$ represent the lower left and the upper right vertices, respectively, of $A_{i}$, the parameter along the transverse curve $x_{\mathcal{L}}$ is bounded by $\sqrt{2z_{1,i}^{SW}z_{2,i}^{SW}}\leq x_{\mathcal{L}}\leq\sqrt{2z_{1,i}^{NE}z_{2,i}^{NE}}$.

The left panel of Figure 16 displays the approximate pdf of $Z$ produced by the method in Section 2.1.2 whereas the right panel shows 25000 samples of $(Z_{1},Z_{2})$ drawn from this pdf through rejection sampling and interpolation. We confirm that when $f_{Z}^{ansatz}$ is propagated through the forward model $Q$, we recover the specified pdf $f_{Q}(q)=-\log(q)$ as guaranteed by Theorem 1. Figure 17 exhibits the histogram of $Z_{1}\cdot Z_{2}$ using 52154 samples of $Z$ obtained from $f_{Z}^{ansatz}$ together with $f_{Q}$.

With the pdf $f_{Z}^{ansatz}$ on $Z$ at hand, one approach to propagate this pdf through the unobserved quantity of interest $\widetilde{Q}$ is to construct a discrete pdf approximation to $\widetilde{Q}$ using the centers of each square $(z_{1,i}^{*},z_{2,i}^{*})$ and their corresponding probabilities $P_{Z}(A_{i})$. This yields distinct outcomes $\widetilde{q}_{j}$ of $\widetilde{Q}$ with weights $P(\widetilde{Q}=\widetilde{q}_{j})=\sum_{i:z_{1,i}^{*}+z_{2,i}^{*}=\widetilde{q}_{j}}P_{Z}(A_{i})$ that need to normalized to obtain the approximate pdf $f_{\widetilde{Q}}^{ansatz}$ of $\widetilde{Q}$.

We now evaluate the performance of the solution from each method to predict the probability law of the unobserved $\widetilde{Q}$. Figure 18 contains 3 subplots, one for each method, plotting $f_{\widetilde{Q}}^{ansatz}$ or $f_{\widetilde{Q}}(\cdot|\{q^{i}\}_{i=1}^{N_{s}})$ together with the true pdf of $\widetilde{Q}$. A stem plot was used for the leftmost subplot to emphasize that the pdf of a discrete random variable was used to accurately approximate $f_{\widetilde{Q}}^{ansatz}$. Quantitatively, the discrepancy between the true pdf and the simulated pdf’s is: $\text{max}|f_{\widetilde{Q}}^{ansatz}-f_{\widetilde{Q}}|\approx 0.2141$ whereas $\text{max}|f_{\widetilde{Q}}(\cdot|\{q^{i}\}_{i=1}^{N_{s}})-f_{\widetilde{Q}}|\approx 0.0786$ for Method 2 (Bayes’+entropy) while $\text{max}|f_{\widetilde{Q}}(\cdot|\{q^{i}\}_{i=1}^{N_{s}})-f_{\widetilde{Q}}|\approx 0.0415$ for Method 3 (Bayes’+known family). It was also observed that increasing the number of available samples $\{q^{i}\}_{i=1}^{N_{s}}$ of the observed quantity of interest $Q$ decreased the discrepancy for the latter 2 methods. The objective of this example was not to conclude which method is better since the available information for each was not identical. Rather, this example justifies the need to specify additional information on $Z$ as was carried out in Sections 3.

We conclude this section by showing the repercussions that may arise if the pdf $f_{Z}^{ansatz}$ on $Z$ is used to predict the pdf of more complicated quantities of interest $\widetilde{Q}$. If we consider $\widetilde{Q}(Z_{1},Z_{2})=Z_{1}^{2}+Z_{2}^{2}$ where $Z_{1},Z_{2}\sim U(0,1)$ and independent, it was proven in [27] that

[TABLE]

Additionally, we also consider $\widetilde{Q}(Z_{1},Z_{2})=\exp(-(Z_{1}^{2}+Z_{2}^{2}))$. Figure 19 contains subplots comparing the true pdf $f_{\widetilde{Q}}$ and histogram of $f_{\widetilde{Q}}^{ansatz}$ computed via samples of $Z$ drawn from $f_{Z}^{ansatz}$ for the two aforementioned unobserved $\widetilde{Q}$’s. We see that $f_{\widetilde{Q}}^{ansatz}$ under/overestimates probabilistic properties of $\widetilde{Q}$ such as $P(\widetilde{Q}\leq 0.5)$ for the left subplot and $P(\widetilde{Q}\geq 0.8)$ for the right subplot, among other properties such as moments of $\widetilde{Q}$.

The above examples are not intended to discredit existing methods for solving stochastic inverse problems. Instead, they underscore why additional information should be specified on $Z$ to yield solutions that can be useful in applications. In fact, existing works such as [24, 25, 26] impose probability distributions on the unknown random quantities after which standard mathematical tools were applied to infer the parameters that characterize these distributions.

5 Conclusion

This work dealt with the stochastic inverse problem of identifying the distribution of $Z$ given probabilistic information of the quantity of interest $Q(Z)$. We surveyed general methods that have been developed to tackle the problem in which no information other than the bounded range of $Z$ is assumed. These methods coped with the ill-posedness of this inverse problem in the following ways: [9, 10, 11] suggested an ansatz that the pdf on the generalized contours of $Q$ is uniform while [15] obtained the probability law of $Z$ from the Karhunen-Loève expansion of the solution to the stochastic PDE. We motivated this work by showing that this lack of additional information entails that the true pdf of $Z$ may not be recovered.

Consequently, we argued that it is necessary for this inverse problem to be posed such that further information on $Z$ is specified to attain solutions that are of practical use. We demonstrated that this specified information can take the form of moments of $Z$ or the family of distributions in which $Z$ resides subject to unknown parameters, among others. Using these information, a conjunction of tools such as Bayes’ theorem, the principle of maximum entropy, and forward uncertainty propagation were utilized to solve the inverse problem in a manner that is consistent with the present information on the quantity of interest and on $Z$. Issues arising from this framework were also highlighted. Finally, we emphasized the need for this specified information by assessing how well the resulting solutions from the discussed methods can predict the probability law of unobserved quantities of interest. It is observed that solving the inverse problem without additional information on $Z$ can lead to solutions that may be unreliable for prediction.

The intention of this work was not to discredit existing contributions in this area but to stress on how we believe the inverse problem must be posed for use in practical applications.

Appendix A Computing the pdf on 1-dimensional contours for Example 1

Fix $x_{\mathcal{L}}$ and consider the contour $\pi^{-1}(x_{\mathcal{L}})$. Let $x_{\mathcal{C}}$ parameterize the arc length of $\pi^{-1}(x_{\mathcal{L}})$ and denote by $\mu(\pi^{-1}(x_{\mathcal{L}}))$ the arc length of $\pi^{-1}(x_{\mathcal{L}})$ so that $0\leq x_{\mathcal{C}}\leq\mu(\pi^{-1}(x_{\mathcal{L}}))$. The conditional pdf $f_{X_{\mathcal{C}}|X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})$ along the contour can be approximated as follows:

• •

Select equidistant points $\{x_{\mathcal{C}}^{(i)}\}_{i=0}^{N}\subset[0,\mu(\pi^{-1}(x_{\mathcal{L}}))]$ such that $x_{\mathcal{C}}^{(0)}=0,x_{\mathcal{C}}^{(N)}=\mu(\pi^{-1}(x_{\mathcal{L}}))$ and $x_{\mathcal{C}}^{(i+1)}-x_{\mathcal{C}}^{(i)}=\frac{\mu(\pi^{-1}(x_{\mathcal{L}}))}{N}$ for $i=0,\dots,N-1$.

• •

From $P(X_{\mathcal{C}}\in(x_{\mathcal{C}}^{(i)},x_{\mathcal{C}}^{(i+1)})\,|\,X_{\mathcal{L}}=x_{\mathcal{L}})=\displaystyle\int_{x_{\mathcal{C}}^{(i)}}^{x_{\mathcal{C}}^{(i+1)}}f_{X_{\mathcal{C}}\,|\,X_{\mathcal{L}}}(x_{\mathcal{C}}|x_{\mathcal{L}})\,dx_{\mathcal{C}}$, we deduce the approximation

[TABLE]

for $x_{\mathcal{C}}\in(x_{\mathcal{C}}^{(i)},x_{\mathcal{C}}^{(i+1)})$ assuming that $N$ is sufficiently large.

• •

To approximate the numerator in (A.1), we construct an infinitesimal region $R$ in $\Gamma=[0,1]^{2}$ bounded by the contours $\pi^{-1}(x_{\mathcal{L}})$ and $\pi^{-1}(x_{\mathcal{L}}+\epsilon)$ for $\epsilon$ sufficiently small and partition $R$ into regions $\{R_{i}\}_{i=1}^{N}\subset\Gamma$ with $R_{i}$ corresponding to the arc lengths between $x_{\mathcal{C}}^{(i-1)}$ and $x_{\mathcal{C}}^{(i)}$ for $i=1,\dots,N$. It then follows that

[TABLE]

where $f_{Z}$ is the joint pdf of $Z_{1},Z_{2}$.