Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples

TL;DR

This paper introduces a new adaptive sampling criterion for Bayesian regression that efficiently identifies the output statistics of unknown functions with minimal samples, especially in high-dimensional spaces.

Contribution

The paper proposes a novel output-weighted sampling criterion that considers existing output values, improving sample efficiency in Bayesian regression for high-dimensional problems.

Findings

The new method outperforms existing criteria in sample efficiency.

It effectively handles high-dimensional input spaces.

The approach adapts to known and unknown output variance scenarios.

Abstract

For many important problems the quantity of interest is an unknown function of the parameters, which is a random vector with known statistics. Since the dependence of the output on this random vector is unknown, the challenge is to identify its statistics, using the minimum number of function evaluations. This problem can been seen in the context of active learning or optimal experimental design. We employ Bayesian regression to represent the derived model uncertainty due to finite and small number of input-output pairs. In this context we evaluate existing methods for optimal sample selection, such as model error minimization and mutual information maximization. We show that for the case of known output variance, the commonly employed criteria in the literature do not take into account the output values of the existing input-output pairs, while for the case of unknown output variance this dependence can be very weak. We introduce a criterion that takes into account the values of the output for the existing samples and adaptively selects inputs from regions of the parameter space which have important contribution to the output. The new method allows for application to high-dimensional inputs, paving the way for optimal experimental design in high-dimensions.

Equations

$$\begin{array}{rl}\mathbf{y}& =\mathbf{Ax}+\mathbf{e},\\ \mathbf{e}& \sim \mathcal{N}(0,\mathbf{V}),\\ p(y|\mathbf{x},\mathbf{A},\mathbf{V})& =\mathcal{N}(\mathbf{Ax},\mathbf{V}).\end{array}$$

$$p(\mathbf{A})$$

$$p(\mathbf{A}|D,\mathbf{V})$$

$$\begin{array}{rl}{\mathbf{S}}_{xx}& =\mathbf{X}{\mathbf{X}}^T+\mathbf{K},\\ {\mathbf{S}}_{yx}& =\mathbf{Y}{\mathbf{X}}^T+MK.\end{array}$$

$${\mathbf{S}}_{xx}$$

$${\mathbf{S}}_{yx}$$

$$\begin{array}{rl}p(\mathbf{y}|\mathbf{x},D,\mathbf{V})& =\mathcal{N}({\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}\mathbf{x},\mathbf{V}(1+c)),\\ c& ={\mathbf{x}}^T{\mathbf{S}}_{xx}^{-1}\mathbf{x}.\end{array}$$

$$p(y|D,\mathbf{V})=\int p(\mathbf{y}|\mathbf{x},D,\mathbf{V})p(\mathbf{x})d\mathbf{x}.$$

$$\hat{\mathbf{R}}=[{\hat{\mathbf{r}}}_1|...|{\hat{\mathbf{r}}}_m]\in {\mathbb{R}}^{m\times m}\text{ and }{\sigma }_i^2,i=1,...,m.$$

$${\mathbf{x}}^T{\mathbf{S}}_{xx}^{-1}\mathbf{x}=\sum _{i=1}^m\frac{{\chi }_i^2}{{\sigma }_i^2}.$$

$${\mathbf{S}}_{xx}^{\prime }=\sum _{i=1}^N{\mathbf{x}}_i{\mathbf{x}}_i^T={\mathbf{S}}_{xx}+\mathbf{h}{\mathbf{h}}^T.$$

$${\mathbf{S}}_{xx}^{\prime }={\mathbf{S}}_{xx}+{\kappa }^2{\mathbf{r}}_j{\mathbf{r}}_j^T.$$

$$\begin{array}{rl}{\mathbf{S}}_{yx}^{\prime }{\mathbf{S}}_{xx}^{\prime -1}& =[{\mathbf{S}}_{yx}+{\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}{\mathbf{hh}}^T][{\mathbf{S}}_{xx}+{\mathbf{hh}}^T{]}^{-1}\\ & ={\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}[{\mathbf{S}}_{xx}+{\mathbf{hh}}^T][{\mathbf{S}}_{xx}+{\mathbf{hh}}^T{]}^{-1}\\ & ={\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}.\end{array}$$

$${\mu }_c$$

$${\sigma }_c^2$$

$${\mu }_c(\mathbf{h})=\text{tr}[{\mathbf{S}}_{xx}^{\prime -1}{\mathbf{R}}_{xx}]=\sum _{i=1}^m\frac{{\sigma }_i^2+{\mu }_{x_i}^2}{n_i+{\delta }_{ik}},h_i={\delta }_{ik}\in {\mathbb{S}}^{m-1},$$

$$h_i={\delta }_{ik},k=\arg \underset{i}{\min }\left(-\frac{{\sigma }_i^2+{\mu }_{x_i}^2}{n_i^2}\right).$$

$$\mathcal{I}(\mathbf{x},\mathbf{y}|D^{\prime },\mathbf{V})={\mathcal{E}}_x+{\mathcal{E}}_{y|D^{\prime }}-{\mathcal{E}}_{x,y|D^{\prime }}.$$

$${\mathcal{E}}_x=\int p(\mathbf{x})\log p(\mathbf{x})d\mathbf{x},{\mathcal{E}}_{y|D^{\prime }}=\iint p(\mathbf{y}|D^{\prime },\mathbf{V})\log p(\mathbf{y}|D^{\prime },\mathbf{V})d\mathbf{y},$$

$${\mathcal{E}}_{x,y|D^{\prime }}=\iint p(\mathbf{y},\mathbf{x}|D^{\prime },\mathbf{V})\log p(\mathbf{y},\mathbf{x}|D^{\prime },\mathbf{V})d\mathbf{x}dy.$$

$${\mathbb{E}}^y[D_{KL}[p(\mathbf{x}|\mathbf{y},D^{\prime })||p(\mathbf{x})]]$$

$$={\int }_y{\int }_{x}p(\mathbf{x},\mathbf{y}|D^{\prime },\mathbf{V})\log \frac{p(\mathbf{x},\mathbf{y}|D^{\prime },\mathbf{V})}{p(\mathbf{x})p(\mathbf{y}|D^{\prime },\mathbf{V})}d\mathbf{x}d\mathbf{y}$$

$$=\mathcal{I}(\mathbf{x},\mathbf{y}|D^{\prime },\mathbf{V}).$$

$${\mathcal{E}}_{x,y}(\mathbf{h})$$

$$=\iint p(\mathbf{y}|\mathbf{x},D^{\prime },\mathbf{V})p(\mathbf{x})\log p(\mathbf{y}|\mathbf{x},D^{\prime },\mathbf{V})d\mathbf{x}d\mathbf{y}+\iint p(\mathbf{y}|\mathbf{x},D^{\prime },\mathbf{V})p(\mathbf{x})\log p(\mathbf{x})d\mathbf{x}d\mathbf{y}$$

$$=\int {\mathcal{E}}_{y|x}(\mathbf{x};\mathbf{h})p(\mathbf{x})d\mathbf{x}+\int p(\mathbf{x})\log p(\mathbf{x})d\mathbf{x}$$

$$={\mathbb{E}}^x[{\mathcal{E}}_{y|x}(\mathbf{x};\mathbf{h})]+{\mathcal{E}}_x.$$

$${\mathcal{E}}_{y|x}(\mathbf{x};\mathbf{h})=\frac{1}{2}\log (1+c{)}^d|2\pi e\mathbf{V}|=\frac{d}{2}\log (1+c(\mathbf{x};\mathbf{h}))+\frac{1}{2}\log |2\pi e\mathbf{V}|.$$

$${\mathbb{E}}^x[{\mathcal{E}}_{y|x}(\mathbf{x};\mathbf{h})]=\frac{d}{2}{\mathbb{E}}^x[\log (1+c(\mathbf{x};\mathbf{h}))]+\frac{1}{2}\log |2\pi e\mathbf{V}|.$$

$$\mathcal{I}(\mathbf{x},\mathbf{y}|D^{\prime },\mathbf{V})={\mathcal{E}}_y(\mathbf{h})-\frac{d}{2}{\mathbb{E}}^x[\log (1+c(\mathbf{x};\mathbf{h}))]-\frac{1}{2}\log |2\pi e\mathbf{V}|.$$

$${\mathbf{C}}_{yy}(D^{\prime },\mathbf{V})$$

$${\mathbf{C}}_{yy}(D^{\prime },\mathbf{V})=\mathbf{V}(1+{\mu }_c(\mathbf{h}))+{\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}{\mathbf{C}}_{xx}{\mathbf{S}}_{xx}^{-1}{\mathbf{S}}_{yx}^T.$$

$${\mathcal{E}}_y(\mathbf{h})$$

Full text

Keywords: Optimal experimental design; Rare extreme events; Bayesian regression; Optimal sampling; Active learning

1 Introduction

For a wide range of problems in engineering and science it is essential to quantify the statistics of specific quantities of interest (or output) that depend on uncertain parameters (or input) with known statistical characteristics. The main obstacle towards this goal is that this dependence is not known a priori and numerical or physical experiments need to be performed in order to specify it. If the problem at hand allows for the generation of many input-output pairs then one can employ standard regression methods to machine learn the input-output map over the support of the input parameters and subsequently compute the statistics of the output.

However, for several problems of interest it is not possible to simulate even a moderate size of input parameters. In this case it is critical to choose the input samples carefully so that they provide the best possible information for the output of interest [4, 9, 1]. A class of problems that belong in this family is the probabilistic quantification of extreme or rare events rising from high dimensional complex systems such as turbulence [16, 7, 19, 2, 3], networks [20], waves [5, 13, 11], and materials or structures [21, 6]. Of course the considered setup is not limited to extreme or rare events but it is also relevant for any problem where the aim is to quantify the input-output relationship with very few but carefully selected data points.

The described setup is a typical example of an optimal experimental design or active learning problem [4]. Specifically, we will assume that we have already a sequence of input-output data and our goal will be to sequentially identify the next most informative input or experiment that one should perform in order to achieve fastest possible convergence for the output statistics. The problem has been studied extensively using criteria relying on mutual information theory or the Kullback–Leibler divergence (see e.g. [15]). More recently another criterion was introduced focusing on the rapid convergence of the output statistics [14]. A common characteristic of these methods is the large computational cost associated with the resulted optimization problem that constrains applicability to low-dimensional input or parameter spaces.

The first objective of this work is to understand some fundamental limitations of popular selection criteria widely used for optimal experimental design (beyond the large computational cost). Specifically, we will examine how well these criteria distinguish and promote the input parameters that have the most important influence to the quantities of interest. The second objective is the formulation of a new, output-weighted selection approach that explicitly and in a controllable manner takes into account, beyond the uncertainty of each input parameter, its effect on the output variables, i.e. the quantities of interest. This is an important characteristic as it is often the case that a small number of input parameters controls a specific quantity of interest. The philosophy of the developed criterion is to exploit the existing samples in order to estimate which input parameters are the most influential for the input and then bias the sampling process using this information.

Beyond its intuitive and controllable character on selecting input samples according to their effect to the output statistics, the new criterion has a numerically tractable form which allows for easy computation of each value and gradient. The latter property allows for the employment of gradient optimization methods and therefore the applicability of the approach even in high-dimensional input spaces. We demonstrate ideas through several examples ranging from linear to nonlinear maps with low and high dimensional input spaces. In particular, we show that the important dependencies of given quantities of interest can be identified and quantified using a very small number of input-output pairs, allowing also for quantification of rare event statistics with minimal computational cost.

2 Setup

Let the input vector $\mathbf{x}\in\mathbb{R}^{m}$ denote the set of parameters or system variables and $\mathbf{y}\in\mathbb{R}^{d}$ be the output vector describing the quantities of interest. The input vector can be thought as high-dimensional with known statistics described by the probability density function (pdf) $p(\mathbf{x})$ that corresponds to mean value $\mu_{x}$ and covariance $\mathbf{C}_{xx}$ (or correlation $\mathbf{R}_{xx}$). In what follows we will use $p$ to denote pdf and an index will be used only if the random variable is not automatically implied by the argument.

A map from the input to the output variables, $\mathbf{y}=\mathbf{T(x)}$, exists and our aim is to approximate the statistics of the output, $p(\mathbf{y})$, using the smallest possible number of evaluations of the map $\mathbf{T}$. We will assume that we have already obtained some input-output pairs which we employ in order to optimize the selection of the next input that one should evaluate. This problem can be seen as an optimal experimental design problem where the experimental parameters that one is optimizing coincide with the random input parameters. All the results/methods presented in this work can be formulated in the standard setup of optimal experimental design in a straightforward way.

The first step of the approach is to employ a Bayesian regression model to represent the map $\mathbf{T}$. Our choice of the Bayesian framework is dictated by our need to have a priori estimates for the model error, as those will be employed for the sample selection criteria. For simplicity we will present our results for linear regression models, although the extension for regression schemes with nonlinear basis functions or Gaussian process regression schemes is straightforward. We formulate a linear regression model with an input vector $\mathbf{x}$ that multiplies a coefficient matrix $\mathbf{A}\in\mathbb{R}^{d\times m}$ to produce an output vector $\mathbf{y}$, with Gaussian noise added:

[TABLE]

We emphasize that for what follows we consider the case of a known noise variance $\mathbf{V}\in\mathbb{R}^{d\times d}$. The case of unknown covariance matrix $\mathbf{V}$ is discussed in Appendix D. The basic setup involves a given data set of pairs $D=\left\{(\mathbf{y}_{1},\mathbf{x}_{1}),(\mathbf{y}_{2},\mathbf{x}_{2}),...,(\mathbf{y}_{N},\mathbf{x}_{N})\right\}$. We set $\mathbf{Y}=[\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N}]$ and $\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N}]$.

For the matrix $\mathbf{A}$ we assume a Gaussian prior with mean $\mathbf{M}\in\mathbb{R}^{d\times m}$ and covariance $\mathbf{K}\in\mathbb{R}^{m\times m}$ for the columns, and $\mathbf{V}$ for the rows. This has the form:

[TABLE]

Then one can obtain the posterior for the matrix $\mathbf{A}$ [17, 12]

[TABLE]

where,

[TABLE]

Essentially, $\mathbf{X}\mathbf{X}^{T}$ is the data correlation of the sample input points $\mathbf{x}_{i},\ i=1,...,N$. We choose $\mathbf{K}=\alpha\mathbf{I}$ ($\mathbf{I}$ is the identity matrix) and $\mathbf{M}=0$, where $\alpha$ is an empirical parameter that will be optimized. Therefore, the above relations take the form

[TABLE]

Based on the above we obtain the probability distribution function for new inputs $\mathbf{x}$:

[TABLE]

Then we can obtain an estimate for the probability density function of the output as

[TABLE]

It is important to emphasize that the output $\mathbf{y}$ is random due to two sources: i) the uncertainty of the input vector $\mathbf{x}$, and ii) the uncertainty due to the model error expressed by the term $c$. The latter is directly related to the choice of data-points $\mathbf{x}_{i},\ i=1,...,N$ and the goal is to choose these points in such a way so that the statistics of $\mathbf{y}$ converge most rapidly.

The most notable property for this model is the fact that the model error is independent of the expected output value of the system. This fact holds true also for Gaussian Process regression (GPR) schemes or regression models that use nonlinear basis functions. This property will have very important consequences when it comes to the optimal input sample selection for the modeling of the input-output relation.

2.1 Properties of the data correlation $\mathbf{S}_{xx}$

One can compute the eigenvectors, $\mathbf{\hat{r}}_{i}$, and eigenvalues, $\sigma_{i}^{2}$, of the data correlation matrix, $\mathbf{S}_{xx}$,

[TABLE]

By applying a linear transformation to the $\mathbf{S}_{xx}$ eigendirections, $\mathcal{X}=\mathbf{\hat{R}}^{T}\mathbf{x}$ we have

[TABLE]

Thus, the eigendirections of $\mathbf{S}_{xx}$ indicate the principal directions of maximum confidence for the linear model. The eigenvalues quantify this confidence: the larger the eigenvalue the slower the uncertainty increases (quadratically) as $\chi_{i}^{2}$ increases. For a new, arbitrary point, $\mathbf{x}_{N+1}=\mathbf{h}$ , added to the family of $\mathbf{x}$ points we will have $\mathbf{X}^{\prime}=[\mathbf{X}|\mathbf{x}_{N+1}]$. By direct computation we obtain

[TABLE]

If the new point belongs to the $j$ eigendirection, $\mathbf{x}_{N+1}=\kappa\mathbf{r}_{j}$, where $\kappa\in\mathbb{R}$ then the new data correlation will be,

[TABLE]

It can be easily checked that under this assumption the new matrix $\mathbf{S}^{\prime}_{xx}$ will have the same eigenvectors. Moreover the $j$ eigenvalue will be ${\sigma^{\prime}}_{j}^{2}=\sigma_{j}^{2}+\kappa^{2}$, while all other eigenvalues will remain invariant. Therefore, by adding one more data-point along a principal direction will increase the confidence along this direction by the magnitude of this new point.

The larger the magnitude of any point we add, the larger its impact on the covariance. One can trivially increase the magnitude of the new points but this does not offer any real benefit. Moreover, in a typical realistic scenario there will be magnitude constraints. To avoid this ambiguity, typical of linear regression problems, we will fix the magnitude of the input points, i.e. $\mathbf{x}\in\mathbb{S}^{m-1}=\{\mathbf{x}\in\mathbb{R}^{m},\left\|\mathbf{x}\right\|=1\}$, so that we can assess the direction of the new points, without being influenced by the magnitude. For nonlinear problems the input points should be chosen from a compact set, typically defined by the mechanics of the specific problem.

3 Fundamental limitations of standard optimal experimental design criteria

Here we consider two popular criteria that can be employed for the selection of the next most informative input sample $\mathbf{x}_{N+1}$. The first one is based on the minimization of the model error expressed by the parameter $c$ (eq. (5)), while the second one is the Kullback–Leibler divergence or equivalently the maximization of the mutual information between input and output variables, which is the standard approach in the optimal experimental design literature [4].

We hypothesize a new input point $\mathbf{x}_{N+1}=\mathbf{h}$. As the corresponding output is not a priori known we will assume that it is given by the mean regression model, $\mathbf{y}_{N+1}=\mathbf{S}_{yx}\mathbf{S}_{xx}^{-1}\mathbf{x}_{N+1}$. The new pairs of data points will be $D^{\prime}=\{D,(\mathbf{x}_{N+1},\mathbf{y}_{N+1})\}$. Under this setup the new model error will be given by $c(\mathbf{x};\mathbf{h})=\mathbf{x}^{T}\mathbf{S}_{xx}^{\prime-1}\mathbf{x}$, where the new data correlation matrix is given by (7). In addition, the mean estimate of the new model will remain invariant since,

[TABLE]

3.1 Minimization of the mean model error

The first approach we will employ is to select $\mathbf{h}$ by minimizing the mean value of the uncertainty parameter $c$ (eq. (5)). Using standard expressions for quadratic forms of a random variable [18] we obtain a closed expression, valid for any input distribution. More specifically, we will have:

[TABLE]

Moreover for the case of Gaussian input we also obtain [18]

[TABLE]

We note that the model uncertainty depends only on the statistics of the input $\mathbf{x}$ (expressed through the covariance $\mathbf{C}_{xx}$) and the samples $\mathbf{X}$ (expressed through the constant (i.e. non-dependent on $\mathbf{x}$) matrix $\mathbf{S}^{\prime}_{xx}$). In other words, the matrix $\mathbf{Y}$ and the output distribution play no role on the mean model uncertainty.

To understand the mechanics of selecting input samples using the mean model error we assume that $\mathbf{R}_{xx}$ is diagonal with eigenvalues $\sigma_{i}^{2}+\mu_{x_{i}}^{2}\ \ i=1,...,d,$ arranged with increasing order. We also assume that samples are collected only along the principal directions of the input covariance. In this case the quantity that is minimized takes the form

[TABLE]

where $n_{i}$ denotes the number of samples in the $i^{th}$ direction. One should choose $\mathbf{h,}$ or equivalently $k,$ according to value of the derivative of $\mu_{c}(\mathbf{h})$. In particular,

[TABLE]

If all directions have been sampled with equal number (e.g. each of the directions have $n_{i}=1$), sampling will continue with the most uncertain direction. After sufficient sampling in this direction, the addition of a new sample will cause smaller effect than sampling the next most important direction and this is when the scheme will change sampling direction. This behavior guarantees that the scheme will never get ‘trapped’ in one direction. It will continuously evolve, as more samples in one direction lead to very small eigenvalue of $\mathbf{S}_{xx}^{\prime-1}\mathbf{R}_{xx}$ along this direction, and therefore sampling along another input direction will cause bigger contribution to the trace.

It is clear that sampling based on the uncertainty parameter $c$ searches only in $\mathbf{x}-$directions with important uncertainty, while the impact of each input variable is completely neglected. Therefore, even directions that have zero effect on the output variable will still be sampled as long as they are uncertain.

3.2 Maximization of the mutual information

An alternative approach for the selection of a new sample, $\mathbf{x}_{N+1}=\mathbf{h}\in\mathbb{S}^{m-1}$, is maximizing the entropy transfer or mutual information between the input and output variables, when a new sample is added [4]:

[TABLE]

where each of the entropies above are defined as,

[TABLE]

[TABLE]

This is also equivalent with maximizing the mean value of the Kullback–Leibler (KL) divergence [9]

[TABLE]

We first compute the entropy of $p(\mathbf{x,y}|D^{\prime},\mathbf{V})$:

[TABLE]

We focus on computing the first term on the right hand side. For the linear regression model, the conditional error follows a Gaussian distribution. From standard expressions about the entropy of a multivariate Gaussian we have

[TABLE]

Therefore,

[TABLE]

In the general case, we cannot compute the entropy of the output, conditional on $D^{\prime}$. To this end, the mutual information of the input and output, conditioned on $D^{\prime},$ takes the form

[TABLE]

This expression is valid for any input distribution and relies only on the assumption of Bayesian linear regression. To compute the involved terms, one has to perform a Monte-Carlo or importance sampling approach, even for linear regression models and Gaussian inputs. This, of course, limits the applicability of the approach to very low-dimensional input spaces. We note that the above expression is valid for the case of known noise variance, $\mathbf{V}$. The case of unknown variance $\mathbf{V}$ is considered in Appendix D.

3.2.1 Gaussian approximation of the output

To overcome this computational obstacle one can consider an analytical approximation of the mutual entropy, assuming Gaussian statistics for the output. This assumption is not true in general, even for Gaussian input, because of the multiplication of the (Gaussian) uncertain model parameters (matrix $\mathbf{A}$) with the the Gaussian input (vector $\mathbf{x}$).

We focus on the computation of the entropy of the output $\mathbf{y}$, so that we can derive an expression for the mutual information. We will approximate the pdf for $\mathbf{y}$ through its second order statistics. Given that the input variable is Gaussian and the exact model is linear the Gaussian approximation for the output is asymptotically accurate. Still, it will help us to obtain an understanding of how the criterion works to select new samples.

We express the covariance of the output variable using the law of total variance

[TABLE]

The first term is the average of the updated conditional covariance of the output variables and it is capturing the regression error. The second term expresses the covariance due to the uncertainty of the input variable $\mathbf{x}$, as measured by the estimated regression model using the input data in $D^{\prime}$.

As we pointed out earlier the mean model using either $D$ or $D^{\prime}$ remains invariant. Therefore, we have

[TABLE]

In this way we have the approximated entropy of the output variable using a Gaussian approximation, which is also an upper bound for any other non-Gaussian distribution with the same second order statistics

[TABLE]

Therefore, we have the second-order statistics approximation of the mutual information in terms of the new sample $\mathbf{h}\in\mathbb{S}^{m-1}$, denoted as $\mathcal{I}_{G}$:

[TABLE]

We observe that the second-order approximation of the mutual information criterion has minimial dependence on the output samples $\mathbf{Y}$. Specifically, (14) depends on the uncertainty parameter $c$ and its statistical moments, as well as the term $\mathbf{S}_{yx}\mathbf{S}_{xx}^{-1}\mathbf{C}_{xx}\mathbf{S}_{xx}^{-1}\mathbf{S}^{T}_{yx}$. However, the latter, is not coupled with the new hypothetical point $\mathbf{h}$ and to this end the minimization of this criterion does not guarantee that the output values will be taken into account in a meaningful way. Instead, the selection of the new sample, depends primarily on minimizing $\mu_{c}=\text{tr}[\mathbf{S}_{xx}^{\prime-1}\mathbf{R}_{xx}]$, always under the constraint $\left\|\mathbf{h}\right\|=1$, a process that depends exclusively on the current samples $\mathbf{X}$ and the statistics of the input $\mathbf{x}$.

Therefore, regions of $\mathbf{x}$ associated with large or important values of the output $\mathbf{y}$ are not emphasized by this sampling approach and the emphasis is given in regions that minimize the mean model error $\mu_{c}$. We note that these conclusions are valid for the second-order approximation of the mutual information criterion, with known output variance, $\sigma_{V}^{2}$. When one considers the mutual information criterion with unknown output variance which has to be inferred using conjugate priors, there is dependence of the criterion on the output vector $\mathbf{Y}$. However, this dependence may be very weak or even zero depending on inference parameters that are optimized based on the data. See Appendix D for details.

3.3 Nonlinear basis regression

Similar conclusions can be made for the case where one utilizes nonlinear basis functions. In this case we assume that the input points ‘live’ within a compact set. Specifically, let the input $\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{m},$ be expressed as a function of another input $\mathbf{z}\in\mathbb{\mathcal{Z}\subset R}^{s}$ where the input value has distribution $p(\mathbf{z})$ and $\mathcal{Z}$ be a compact set. One can choose a set of basis functions

[TABLE]

In this case the distribution of the output values will be given by :

[TABLE]

The mean of the model uncertainty parameter $c=\phi(\mathbf{z})\mathbf{S}^{-1}_{\phi\phi}\phi(\mathbf{z})^{T}$ will become

[TABLE]

where

[TABLE]

and

[TABLE]

Following the same steps as we did for the linear model, we will have, first for the conditional entropy (assuming that the model noise in the nonlinear case is Gaussian)

[TABLE]

Therefore,

[TABLE]

The exact expression for the mutual information for the nonlinear case will be:

[TABLE]

To perform, the second-order statistical approximation for the entropy $\mathcal{E}_{y}$, we follow the same steps as for the linear model case to obtain

[TABLE]

The sampling strategy is more complicated in this case due to the nonlinearity of the basis elements. However, even in the present setup the sampling depends exclusively on the statistics of the input variable $\mathbf{z}$ and the form of the basis elements $\phi$. The measured output values of the modeled process do not enter explicitly into the optimization procedure for the next sample, in the same fashion with the linear model.

4 Optimal sample selection considering the output values

We saw that selecting input samples based on either the mean model error or the mutual information does not effectively take into account the output values of the existing samples. Our goal is to develop an approach that i) will give emphasis on the output values of the existing samples, and ii) will be computationally tractable. In [14] a similar problem was considered where the goal was to design a sampling method that will accelerate the convergence of the pdf in regions associated with rare events. In particular the following steps were followed in [14]:

Using the existing samples the authors obtain an estimate of the map (denote it as $\mathbf{y}_{0}(\mathbf{x})$), as well as, the map-estimation-error, $\sigma_{y_{0}}^{2}(\mathbf{x})$, at every point $\mathbf{x}$. 2. 2.

This map-estimate and its error are then used to estimate the output pdf, denoted as $p_{y_{0}}(\mathbf{y})$, as well as, the pdf of the perturbed map along the direction of the map-estimation-error, denoted as $p_{y_{0}^{+}}(\mathbf{y}).$ 3. 3.

A new hypothetical sample point, $\mathbf{h}$, was assumed and its impact first on the map-estimation-error, $\sigma_{y_{0}}^{2}(\mathbf{x;\mathbf{h}}),$ and then on the pdf of the perturbed map, $p_{y_{0}^{+}}(\mathbf{y;h})$, was quantified. 4. 4.

Then the goal was to select the new sample that minimizes the distance of the two pdfs, i.e. between $p_{y_{0}^{+}}(\mathbf{y;h})$ and $p_{y_{0}}(\mathbf{y})$. As a distance the authors considered the $L_{1}$ difference of the logarithms of the two pdfs, instead of the KL divergence. The reason was that in the KL divergence the difference of the logarithms is multiplied with the pdf itself and therefore rare events play a less important role on the value of the criterion. By considering only the difference of the logarithms gave more emphasis in the regions associated with rare events.

The approach was very effective on computing the rare event properties (tails of the pdf) for arbitrary quantities of interest with a very few samples. However, it was limited by the large cost related to the computation of the two pdfs mentioned above, which was performed with direct Monte-Carlo methods. For this reason the method could be applied in problems with relatively low-input dimensionality.

In the present work we are going to build on [14] to derive a new criterion that follows the same principles as the one just described but it is also computationally tractable and can be used beyond the context of rare events, i.e. for general optimal experimental design problems with a large number of parameters.

Specifically, we are going to apply the following steps:

We will employ an asymptotic form of the criterion in [14] that will provide, analytically, the distance of the pdf logarithms, i.e. the pdf of the map-estimate and the pdf of the asymptotically-perturbed map-estimate. 2. 2.

Using standard inequalities for norms of derivatives we will bound the asymptotic form of the criterion by a more intuitive and tractable form. This new form has an interesting interpretation as it naturally weights the importance of the estimated map-error by the pdf of the input but also the inverse of the pdf of the output. In this way more emphasis is given to inputs associated with large values of the output. 3. 3.

The final step of our analysis is to demonstrate how the derived bound can be analytically approximated in terms of second-order properties of the input pdf, as well as, second-order properties of the estimated output. This last step is the most cumbersome but also crucial in order to apply the method in high-dimensional problems, as the analytical approximation of the criterion allows for the application of gradient optimization methods.

In Figure 1 we provide a sketch of the main idea. A two dimensional input space is shown where the output is a function that primarily depends on one input variable. The input variable has important variance in both dimensions. While methods relying on mutual information give emphasis primarily on the statistics of the input variable, resulting an-equally-good approximation of the map over all input dimensions through a uniform coverage of all input dimensions, an output-weighted approach assess the important of each candidate input sample by its effect on the statistics of the output, i.e. the quantity of interest. In this way the output values for the observable of interest are taken into account explicitly and in a controllable manner.

Derivation of the output-weighted criterion

Our goal is to compute samples that accelerate the convergence of the output statistics, expressed by the probability density function, $p(\mathbf{y})$. To measure how well this convergence has occurred, we are going to rely on the distance between the probability density function of the mean model

[TABLE]

and the perturbed model along the most important direction of the model uncertainty (dominant eigenvector of $\mathbf{V}$), denoted as $\mathbf{r_{V}}$:

[TABLE]

where $\beta\ll 1$ is a small scaling factor. The corresponding probability density functions, $p_{{y}_{0}}(\mathbf{y})$ and $p_{{y}_{+}}(\mathbf{y})$ will differ only due to errors of the Bayesian regression, which vary as $\mathbf{h}$ changes. It is therefore meaningful to select the next sample that will minimize their distance. Moreover, as we are interested on capturing the probability density function equally well in regions of low and large probability we will consider the difference between the logarithms. Specifically, we define, [14],

[TABLE]

where $S_{y}$ is a finite domain over which we are interested to define the criterion. Note that the latter has to be finite in order to have a bounded value for this distance. It can be chosen so that it contains several standard deviations of the output process. The defined criterion focuses exactly on our goal which is the convergence of the output statistics, while the logarithm guarantees that even low probability regions have converged as well. This criterion for selecting samples was first defined in [14] and was shown that it results in a very effective strategy for sampling processes related to low-probability extreme events. However, it is also associated with a very expensive optimization problem that has to be solved in order to minimize this distance. Apart from the cost, its complicated form does not allow for the application of gradient methods for optimization and therefore it is practical only for low-dimensional input spaces where non-gradient methods can be applied.

Here one of our goals is to study its relationship with existing criteria. We are also aiming to bound it by a more trackable form that is applicable for gradient optimization methods.

To study the relationship of the criterion (22) with the KL divergence, we note that for bounded probability density functions the following inequality holds

[TABLE]

where $\kappa$ is a constant. To this end, the criterion based on the difference of the logarithms is more conservative (i.e. harder to minimize) compared with the KL divergence (defined over the same domain).

Our next goal is to bound the $D_{Log^{1}}$ criterion by one that is more tractable to optimize. We consider the criterion (22) for an asymptotically small value of $\beta$. This form of the criterion essentially expresses the infinitesimal difference between the mean model and the infinitesimally perturbed model by $\beta{\sigma_{{y}}^{2}}$. To compute analytically the value of the criterion for $\beta\rightarrow 0$ we employ an asymptotic result originally obtained in [14] for the study of the criterion for a large number of input samples, i.e. very small ${\sigma_{{y}}^{2}}$. For this case, or equivalently the case where $\beta$ is very small that we are interested here, we have the asymptotic form (Theorem 1 in [14])

[TABLE]

where,

[TABLE]

is the conditional variance (on $\mathbf{x}$) if the output is scalar or the trace of the output conditional covariance matrix in the general case, while $\mathbf{y}_{0}(\mathbf{x})$ is the mean model from the input-output data collected so far.

Using standard inequalities for the derivatives of differential functions one can bound the derivative in (24). Specifically, if the function $\mathbb{E}[{\sigma^{2}_{{y}}}(\mathbf{x})\cdot\boldsymbol{1}_{{{y}_{0}}({x})=\text{s}}]$ has uniformly bounded second derivative (with respect to a hypothetical new point $\mathbf{h}$), and $p_{{y}_{0}}(\mathbf{s})$ has not zeros or singular points, there exists a constant $\kappa_{0}$ such that ([10], Theorem 3.13, p. 109)

[TABLE]

Moreover,

[TABLE]

where $S_{x}$ is the inverse image of the domain $S_{y}$ through the map, $\mathbf{y}_{0}(\mathbf{x})$. Based on this we obtain the output-weighted model-error criterion, which bounds (i.e. it is more conservative) the original criterion (22) as well as the information-based criterion:

[TABLE]

In practice, $S_{x}$ is chosen as $\mathbb{R}^{m}$ or the support of the input pdf, $p_{x}$. Because of the inequality (25) we can conclude that converge of $Q[\sigma_{{y}}^{2}]$ also implies convergence of the metric $D_{Log^{1}}(\mathbf{y}_{+}\|\mathbf{y}_{0})$. However, the $Q$ criterion is much easier to compute compared with $D_{Log^{1}}(\mathbf{y}_{+}\|\mathbf{y}_{0})$ and it can be employed even in high-dimensional input spaces. With the modified criterion the output data and their pdf is taken into account explicitly. In particular, the conditional variance (or uncertainty) of the model at each input point $\mathbf{x}$ is weighted by the probability of the input at this point, $p_{x}(\mathbf{x})$, as well as the inverse of the estimated probability of the output at the same input point, $p_{{y}_{0}}({\mathbf{y}_{0}}(\mathbf{x}))$.

The term in the denominator comes as a result of considering the distance between the logarithms in (22). If we had started with $D_{KL}(\mathbf{y}_{0}\|\mathbf{y}_{+};{S_{y}})$, we would have cancellation of this important term. Note that a relevant approach, based on the heuristic superposition of the outcome and the mutual information criterion was presented in [22]. However, there is no clear way how the two terms should be weighted or in what sense the outcome can be superimposed to the information content. We emphasize that the presented framework is not restricted to linear regression problems and it can also be applied to Bayesian deep learning problems (a task that will not be considered in this work). In addition, we have not made any assumption for the distribution of the input $\mathbf{x}$.

A simple demonstration

To illustrate the properties of the new criterion we consider the map

[TABLE]

Note that the $x_{2}$ variable is more important than $x_{1}$ in determining the value of the output, given that the two input variables have the same variance. It is therefore intuitive to require more accuracy for the second direction. However, the information distance or entropy based criteria take into account only the input variable statistics to select the next sample, in which case, both directions will have equal importance. This is illustrated in Figure 2, where we present contours of the exact map, $T(\mathbf{x})$, as well as, of the input pdf, $p(\mathbf{x})$. We also present the contours of the output pdf conditional on the input, $p_{y}(T(\mathbf{x}))$ (bottom left), and the weight that is used in the criterion $Q$. Clearly, relying on the sampling criterion that uses only information about the input will not be able to approximate the map in the most important directions. On the other hand, we observe that the weight used in the $Q$ criterion takes into account explicitly the importance of both the input variable statistics but also the information that one has estimated so far from the input-output samples. Here we used the exact map $T(\mathbf{x})$ to demonstrate the weight function but in a realistic scenario the estimated mean model $y_{0}(\mathbf{x})$ will be used to approximate the output pdf.

4.1 Approximation of the criterion for symmetric output pdf

Our efforts will now focus on the efficient approximation of the criterion $Q$. To simplify the presentation we will focus on the scalar output case, $d=1$. The first step of the approximation focuses on the denominator. This term introduces the dependence to the output data and acts as a weight to put more emphasis on regions associated with large deviations of $y(\mathbf{x})$ from its mean. We will approximate the weight, $p^{-1}_{{y}_{0}}({y})$, by a quadratic function that optimally represents it over the region of interest, $S_{y}$. Therefore, for the scalar case we will have

[TABLE]

where $p_{1},p_{2}$ are constants chosen so that the above expression approximates the inverse of the output pdf optimally over the region of interest. We use this expression into the definition of $Q$ (eq. (26)) and obtain the approximation

[TABLE]

Note that the first term does not depend on the output values but only on the input process. It is essentially the same term that appears in the entropy based criteria. The second term however depends explicitly on the deviation of the output process from its mean and therefore on the output data. Specifically, it has large values in regions of $\mathbf{x}$ where the output process has important deviations from its mean, essentially promoting the sampling of these regions. The two constants $p_{1},p_{2}$ provide the relative weight between the two contributions. They are computed for a Gaussian approximation of the output pdf in Appendix B. For the case where the pdf $p_{y}$ is expected to have important skewness, i.e. asymmetry around its mean, a linear term can be included in the expansion of $p^{-1}_{{y}_{0}}({y})$, so that this asymmetry is reflected in the sampling process.

4.2 Linear regression with Gaussian input

For the case of liner regression the first term in the criterion (28) will take the form

[TABLE]

where we have considered the case of a scalar output with $\mathbf{V}=\sigma_{V}^{2}$. The second term of the criterion (28) will take the form

[TABLE]

where $c_{0}$ is a constant that does not depend on $\mathbf{h}$

[TABLE]

We observe that the second term depends on fourth order moments of the input process $\mathbf{x}$ but also on the output values of the samples $Y$. This term can be computed in a closed form for the case of Gaussian input. Specifically,

[TABLE]

where $\mathbf{x}^{\prime}=\mathbf{x}-\mu_{x}$ and $p(\mathbf{x}^{\prime})$ is the zero-mean translation pdf of the original one. The second and third term on the right hand side vanish as they consist of third order central moments of a Gaussian random variable. For the first term we employ a theorem for the covariance of quadratic forms which gives for two symmetric matrices, $A$ and $B$ [18]:

[TABLE]

Therefore,

[TABLE]

From this equation, it follows,

[TABLE]

In addition, the last term becomes

[TABLE]

We collect all the computed terms and obtain

[TABLE]

This is the form of the $Q$ criterion under the assumption of Gaussian input for the case of linear regression. For the case of zero mean input it becomes

[TABLE]

The coefficients $p_{1},p_{2}$ are determined using the output pdf of the estimated model through the samples $D$ (eq. (27)), i.e. the pdf of $y_{0}(\mathbf{x})$. Note that the exact form of the output pdf, used in the criterion, is not important at this stage as it only defines the weights of the criterion $Q(\mathbf{h})$. For a Gaussian approximation of the output process the coefficients are given in Appendix B.

4.3 Nonlinear regression with Gaussian input

For the case of regression with non-linear basis the first term in the criterion (28) will take the form

[TABLE]

where we have considered the case of a scalar output with $\mathbf{V}=\sigma_{V}^{2}$. The second term of the criterion (28) will take the form

[TABLE]

where we expressed the integral using the pdf for the basis elements $\phi$. In this way the integral is now expressed exactly as in the linear regression case. To obtain a closed approximation we approximate the pdf for $\phi$ through its second-order statistics (i.e. approximate $p(\phi)$ with a Gaussian pdf). The analysis shown for the linear case with Gaussian input is then valid leading to the following expression for the $Q$ criterion:

[TABLE]

So, for a given basis $\phi(\mathbf{z})$ one needs first to obtain the mean vector $\mu_{\phi}$ and covariance $C_{\phi\phi}$ using the expressions (3.3) and then follow the same steps as in the linear case. The expression for the gradient of the $Q$ criterion, under general choice of $\phi(\mathbf{z})$ is given in Appendix A.

5 Examples

5.1 Linear map with a 2d input space

To demonstrate the properties of the new criterion we first consider the two-dimensional problem

[TABLE]

We consider two cases of parameters

• •

Case I : $\hat{a}_{1}=0.8,\hat{a}_{2}=1.3$, and $\sigma_{1}^{2}=1.4,\sigma_{2}^{2}=0.6.$

• •

Case II: $\hat{a}_{1}=0.01,\hat{a}_{2}=2.0$, and $\sigma_{1}^{2}=2.0,\sigma_{2}^{2}=0.2.$

The two cases are presented in Figure 3 in polar coordinates. The black arrows indicate the principal directions of the input covariance, scaled according to the eigenvalues of the covariance matrix, while the green arrow indicates the direction of the gradient of the map $T(\mathbf{x})$. While for the first case the contributions of both input variables to the output are comparable and thus sampling is important for both of them, for Case II the contribution of the first input variable is negligible. However, this input variable is the one with the highest uncertainty.

For each case we assess four adaptive sampling strategies according to the criteria:

The directly computed mutual information, $\mathcal{I}(\mathbf{x,y}|D^{\prime})$ is maximized in $\mathbb{S}^{1}$, 2. 2.

The second-order statistical approximation of the mutual information, $\mathcal{I}_{G}(\mathbf{x,y}|D^{\prime})$ is maximized in $\mathbb{S}^{1}$, 3. 3.

The uncertainty parameter, $\mu_{c}(\mathbf{h})$ is minimized in $\mathbb{S}^{1}$, 4. 4.

The output-weighted model error criterion $Q(\mathbf{h})$ is minimized in $\mathbb{S}^{1}$.

For the $Q$ criterion we choose $p_{1}=0$ and $p_{2}=1$ to emphasize the role of the second, new term that takes into account the output samples. This case of parameters corresponds to the case where we optimally approximate $p_{y_{0}}^{-1}$ over the full real axis, i.e. $\beta\rightarrow\infty$ using the notation of Appendix B. We denote this criterion as $Q_{\infty}$. We also compare with a Monte-Carlo approach where samples are randomly generated from the input distribution of $\mathbf{x}$ and then normalized so they belong in $\mathbb{S}^{1}$.

For the adaptive strategies based on $\mu_{c}$ and $Q_{\infty}$ we use the analytical expressions (9) and (31), respectively together with their gradient computed in Appendix A. This allowed us to utilize gradient based optimization methods. For the adaptive strategies based on the mutual information, $\mathcal{I}(\mathbf{x,y}|D,\mathbf{h})$ and its second-order approximation, $\mathcal{I}_{G}(\mathbf{x,y}|D,\mathbf{h})$, we used a random sampling approach and equations (52) and (14), respectively. Specifically, we generated $10^{5}$ samples from the input distribution $\mathbf{x}$ and utilized the exact expression:

[TABLE]

which was numerically computed as an ensemble average. For the computation of $\mathcal{I}$ we also generated $10^{5}$ realizations of the vector $\mathbf{a}=(a_{1},a_{2})$, which according to the linear regression method follows the normal distribution (3):

[TABLE]

Next, we compute a pdf approximation for $y$ using the generated samples and the kernel smoothing functions method [8], and we approximated the entropy of the resulted distribution by direct numerical integration. Note that this additional step, required for $\mathcal{I}$, has a vast computational cost. Most importantly, because of the absence of an analytical expression for the gradient of $\mathcal{I}$, its application to high dimensional inputs is impossible. For this example the next sample vector was parametrized as $\mathbf{h}=[\cos(\theta),\sin(\theta)]$ and the criterion was optimized by direct selection of the maximum value over a one dimensional grid for $\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$.

All four adaptive strategies are initiated with $4$ random samples drawn from the $\mathbf{x}$ distribution. For each case we present the average error curve over $400$ experiments, i.e. experiments with different sets of initial samples (Figure 4: left panels). The standard deviation for each error curve is also presented in the shaded region. For the four adaptive sampling strategies we also present the pdf of the orientation of the samples $\mathbf{h}=(h_{1},h_{2})$, i.e. $\theta=\arctan(h_{2}/h_{1})$.

In both cases of parameters we observe that the strategy based on the $Q_{\infty}$ criterion outperforms the other three adaptive strategies. This difference in performance is even more pronounced for Case II, where one of the input variables has negligible contribution but large uncertainty. An interesting observation is that the Monte-Carlo strategy, performs as good as the $\mu_{c}$ criterion and the mutual information criteria, $\mathcal{I}$ and $\mathcal{I}_{G}$. This is not a surprise given that $\mathcal{I}_{G}$ depends primarily on $\mu_{c}$ and the latter is designed to give more emphasis on input directions with large uncertainty, without taking into account their expected contributions to the output, similarly with Monte-Carlo. The same conclusions hold even for the full mutual information criterion, an indication that although $\mathcal{I}$ partially incorporates the output samples, it does not do it in a useful way.

Similar observations can be made if we examine the pdf for the input samples obtained from the four adaptive strategies. We can see that for each case of parameters the strategies based on $\mu_{c}$, $\mathcal{I}$, and $\mathcal{I}_{G}$ behave very similarly and tend to place input samples in the direction of larger uncertainty. On the other hand, the $Q_{\infty}$ based strategy is placing more samples in directions that compromise between large expected impact to the output variable but also with important uncertainty. We note that all the cases presented here correspond to a known output variance, $\sigma_{V}^{2}$. Results for this example corresponding to unknown variance $\sigma_{V}^{2}$ are presented in Appendix D.

5.2 A high-dimensional linear problem

The next problem to demonstrate the optimal sampling approach is a 20 dimensional linear function. Note that for this case optimization of the mutual information is an impossible task given the fact that the full expression for the mutual information is hard to optimize in the absence of expressions for its gradient.

Specifically, we consider the system

[TABLE]

where the coefficients and input variances are chosen as

[TABLE]

This system represents a typical high dimensional case, where we have some very influential degrees of freedom and some that have negligible impact to the output variable. The energy of these coefficients is typically not related to their influence to the output variable. In Figure 5 we present the coefficients and input variances.

For the observation noise we consider two cases:

• •

Case I: $\sigma_{\epsilon}^{2}=0.05$ (accurate observations)

• •

Case II: $\sigma_{\epsilon}^{2}=0.5$ (noisy observations)

Given that $\sum_{m=1}^{20}\hat{a}_{m}^{2}\sigma_{m}^{2}=0.0272$ the first case corresponds to relatively accurate observations while the second is a highly noisy observations case. We expect the adaptive sampling approach to be more valuable for the second case, given that for the first case we need very few samples anyway.

We apply the adaptive criteria after we have obtained one sample per input direction, to guarantee that the matrix $\mathbf{S}_{xx}$ is invertible. Then we run each numerical experiment $L=400$ times to make sure that the randomness due to the observation noise does not favor any method.

In Figure 6 we present the performance of the sampling approach based on $\mu_{c}$ and $Q_{\infty}$, as well as a direct Monte-Carlo approach. For the first case (accurate observations), shown in the left plot, we note a clear advantage of the $Q_{\infty}$ sampling approach that takes into account the output samples. This advantage is more pronounced in the second case of noisy observations, where the approach using the $Q_{\infty}$ criterion obtains an order of magnitude higher accuracy from the very first samples. Note that the $\mu_{c}$ sampling strategy is comparable with the Monte-Carlo approach, since it does not take into account the output samples.

The same conclusions can be obtained if we observe the variance of the $\mathbf{h}_{N}$ components, over different runs of the numerical experiments, $l=1,...,L$ (here $L=400$), i.e. over different realizations of the observation noise:

[TABLE]

Results are shown in Figure 7. Sampling according to $\mu_{c}$ results in a distribution that is following the shape of the variance $\sigma_{m}$. Specifically, the scheme iteratively changes input directions based on their variance, starting from the most energetic ($m=1$ and $m=20$) and moving towards the less energetic ones ($m=10$). Then the loop begins again, until all the input directions are equally well sampled, after which point the sampling is random.

Sampling according to $Q_{\infty}$, on the other hand, is performed in one loop starting from the most energetic directions, but giving more emphasis in the input directions close to $m=20$ that have both high energy and large contribution to the output $y$. This ‘asymmetry’ in the sampling results in significantly faster convergence compared with the Monte-Carlo method or the $\mu_{c}$ criterion.

The effect of the domain selection $S_{y}$ is discussed in detail in Appendix C, where a parametric study is also shown. Specifically, if we utilize a finite number of standard deviations to optimally approximate $p_{y_{0}}^{-1}$ (eq. (27)), by employing $Q_{\beta}$ (with $\beta$ finite), the term $\mu_{c}$ in the criterion improves the behavior of sampling for large $N$ (Figure 11).

5.3 A 2d nonlinear problem with nonlinear basis functions

The next application involves a nonlinear map with a 2D input space. Specifically, we consider the two-dimensional nonlinear problem

[TABLE]

We consider two cases of parameters

• •

Case I : $\hat{a}_{1}=10^{-2},\hat{a}_{2}=5,\hat{a}_{3}=0,\hat{a}_{4}=10^{2}$, and $\sigma_{1}^{2}=2.10^{-1},\sigma_{2}^{2}=5.10^{-3}.$

• •

Case II: $\hat{a}_{1}=10,\hat{a}_{2}=5,\hat{a}_{3}=0,\hat{a}_{4}=10^{2}$, and $\sigma_{1}^{2}=2.10^{-3},\sigma_{2}^{2}=5.10^{-3}.$

In the first case the output has very weak dependence on the first variable although the latter has very large variance. Moreover, the second variable has significantly smaller variance but plays the dominant role for the output. On the other hand, for the second case, both input variables play an important role and their variance is also comparable. The exact pdf computed with an expensive Monte-Carlo simulation is shown in Figure 8. Both distributions are characterized by heavy tails.

We setup a non-linear Bayesian regression scheme with the following odd basis functions:

[TABLE]

This set of basis functions contains all the odd monomials with order less or equal to 3. We observe that for the nonlinear case although the input space is two-dimensional, the regression is performed in a five-dimensional space. To avoid an ill-conditioned matrix $\mathbf{S}_{\phi\phi}$ we assume a prior with covariance $\mathbf{K}=\alpha\mathbf{I}$ (see eq. (4)). For the cases considered we set $\alpha=10^{-1}$. For each case we first choose randomly (using the distribution of $\mathbf{z}$) two samples. Then we use the criteria based on $\mu_{c}$ (equation (17)) and $Q_{0.01}$ (equation (33)) to optimize 100 samples. For each step we employed a gradient-based optimization using the expressions presented in Appendix A and we restricted the samples to in the disk: $|\mathbf{z}|\leq 2$. For each criterion we performed 200 optimization cycles, i.e. we computed for each criterion the full sequence of 100 samples, 200 times and we computed the statistics of the error (mean and standard deviation), so that the results are not sensitive on the randomness due to observational noise or the initial samples.

The convergence analysis for each criterion is presented in Figure 9. The left plot shows the convergence of the two methods for the parameters of the first numerical experiment (Case I). Each curve is the mean error computed from the 200 optimization cycles, while the shaded area indicates the spread across different runs. Note that for Case I there is only input variable ($z_{2}$) that plays a dominant role on the output, while the other input variable has negligible effect but important variance. In agreement with the results of the linear problems, the samples based on the $Q_{0.01}$ criterion, achieve better performance as they rapidly align with the direction that has the most important influence on the output.

This is not the case for the samples based on the $\mu_{c}$ criterion that align primarily with the directions of importance variance, resulting in a slower convergence. For the Case II parameters both input directions have comparable variance and comparable effect to the output. In this case, as expected, the two criteria have comparable performance. This is clearly demonstrated by the right plot in Figure 9. Finally, in Figure 10 we demonstrate the convergence of the pdfs for the first case of parameters.

6 Conclusions

We have analyzed fundamental limitations of popular criteria for samples selection, employed in the optimal experimental design community. These criteria are based on maximization of entropy-based quantities, typically having the form of mutual information between input and output variables. Specifically, we have shown that beyond the large computational cost associated with these criteria that restricts their applicability to very low-dimensional problems, there is weak dependence of the induced sampling process to the output values of the existing samples when the variance of the output noise is assumed to be known. Even for the case of unknown variance though the dependence on the output values is not controllable and can become very weak. In this way, directions of the parameter space that contribute the most to the output may not be emphasized.

Motivated by these limitations, we have presented a new criterion for optimally selecting training samples that significantly accelerates the convergence of Bayesian regression schemes with respect to the state of the art. The criterion explicitly takes into account the fact that different input parameters have different impact to the output of interest, with some of them being much more influential than others. In this way, it places more samples towards the influential parameters, which are also characterized by important uncertainty. In addition, the introduced criterion is more practical to compute, compared with mutual information criteria, as its gradient can be analytically derived, allowing for the employment of gradient optimization methods. Therefore, the new method allows for the optimization of samples, even for a large number of input parameters, paving the way for optimal experimental design and active learning in high-dimensions.

Future work will focus on the formulation of the presented framework on the training of deep neural networks. The presented approach is expected to have important impact in application areas such as optimal experimental design for systems where very few experiments are available (e.g. biology), adaptive sampling in complex environments with multiple objectives, uncertainty quantification and extreme event statistics in challenging problems such as fatigue-crack, coastal flooding, critical network events, and others.

Acknowledgments

The author would like to thank Prof. Munther Dahleh, Prof. Sanjoy Mitter and Dr Mustafa Mohamad for several stimulating discussions. This work was initiated during a sabbatical visit of the author at ETH, hosted by Prof. George Haller, which is greatfully acknowledged. It was later supported by a Doherty Career Development Chair, a Mathworks Faculty Research Innovation Fellowship and the ARO MURI Grant W911NF-17-1-0306. The detailed comments of the referees led to several improvements and are also highly appreciated.

Appendix A: Gradient of trace criteria

Several criteria in this work take the form

[TABLE]

where $\mathbf{C}$ is a symmetric matrix and $\mathbf{S}_{xx}^{\prime}=\mathbf{S}_{xx}+\mathbf{hh}^{T}.$ The gradient of this expression can be explicitly computed. We first note that

[TABLE]

where,

[TABLE]

In this way we will have

[TABLE]

For the case of nonlinear regression $\mathbf{h}=\phi(\mathbf{z})$. Then

[TABLE]

In this way we will have

[TABLE]

Appendix B: Optimal approximation of $p_{y}^{-1}$

To approximate the inverse of the output pdf, $\frac{1}{p_{y}}$ over $S_{y}$ we are going to use a least square approach. Specifically, we will assume a symmetric output density and we will employ the approximation

[TABLE]

The constant $p_{2}$ is chosen so that we have optimal least square approximation over the interval $[\mu_{y},\mu_{y}+\beta\sigma_{y}]$ where $\beta$ is a fixed parameter that defines the output region of interest, $S_{y}$. By direct minimization we obtain

[TABLE]

For the case of Gaussian output the above expression takes the form

[TABLE]

In this way the least square approximation over the interval $[\mu_{y},\mu_{y}+\beta\sigma_{y}]$ will be

[TABLE]

We denote the criterion with the coefficients obtained from this approximation as $Q_{\beta}$. Specifically,

[TABLE]

For large values of $\beta$ we have $\lim_{\beta\rightarrow\infty}\kappa=\infty$ and the $Q$ criterion is essentially dominated by the output-dependent term. For the case of very small $\beta$ we have $\lim_{\beta\rightarrow 0}\kappa=\frac{\sqrt{2\pi}}{\sigma_{y}}$.

Appendix C: Effect of the weights in the $Q$ criterion

Here we present additional results for Case I of the high dimensional system. Specifically, in Figure 11 we present the performance of the sampling algorithm according to various choices of the $\beta$ parameter. The corresponding sampling patterns are shown in (Figure 12). All cases presented are averaged over $L=500$ numerical experiments to remove the effect of observational noise.

We observe that for $\beta=2$ or $\beta=3$ the performance for small $N$ is very close to the one obtained with $Q_{\infty}$. For larger $N$, however, the performance with finite $\beta$ is improved. In addition to the finite $\beta$ case, we also present two cases with fixed $p_{1},p_{2}$. Specifically, we have $Q_{0.01}$ representing the case $p_{1}=0.01$ and $p_{2}=1$, while $Q_{0.001}$ represents the case $p_{1}=0.001$ and $p_{2}=1$. It is interesting to observe that $Q_{0.001}$ has better performance for small $N$ compared with all other criteria.

Appendix D: The case of unknown $\sigma_{V}^{2}$

Here we consider the case of a priori unknown covariance $\sigma_{V}^{2}$. To simplify the presentation we will restrict our analysis to scalar output and vector input. We formulate a linear regression model with an input vector $\mathbf{x}$ that multiplies a coefficient vector $\mathbf{a}$ to produce an output scalar ${y}$, with Gaussian noise added that has unknown variance:

[TABLE]

For the vector $\mathbf{a}$ we assume a Gaussian prior with mean $\mathbf{m=0}$ and covariance $\mathbf{K=I}\alpha$, where $\alpha$ is a parameter that will be optimized with respect to the evidence. This has the form:

[TABLE]

A conjugate prior for $\sigma_{V}^{2}$ is the inverse Gamma (or inverse Wishart in the multi-dimensional case):

[TABLE]

where $q$ is a normalization constant: $q=\frac{\sigma_{0}^{\nu}}{2^{\frac{\nu}{2}}\Gamma(\nu/2)}$, $\sigma_{0}^{2}$ is a prior value for $\sigma_{V}^{2}$ and $\nu$ is a parameter that is optimized via empirical Bayes.

The posterior for the unknown coefficients will be given by eq. (3), while the posterior for $\sigma_{V}^{2}$ takes the form

[TABLE]

where, $\sigma_{\mathbf{Y}|\mathbf{X}}^{2}=\mathbf{Y}\mathbf{Y}^{T}-(\alpha+1)^{{-1}}\mathbf{S}_{yx}\mathbf{S}_{{xx}}^{-1}\mathbf{S}^{T}_{yx}$. Multiplying the predictive distribution (5) with the posterior for $\sigma^{2}$ and integrating over this argument we have the predictive pdf ($t-$distributed):

[TABLE]

where the parameters $(\alpha,\nu)$ are chosen by maximizing the evidence (set the gradient of $p(\mathbf{Y}|\mathbf{X})$ equal to zero), which leads to the following fixed-point problem:

[TABLE]

where $\Psi(x)=\frac{d\log{\Gamma{(x)}}}{dx}$ is the digamma function.

Selecting inputs by maximization of the mutual information

We follow the same steps as in the known variance case. We hypothesize a new sample, $\mathbf{x}_{N+1}=\mathbf{h}\in\mathbb{S}^{m-1}$ and the goal is maximizing the entropy transfer or mutual information between the input and output variables, when this new sample is added. For this case of a priori unknown variance $\sigma_{V}^{2}$ we denote the mutual information as $\mathcal{\hat{I}}$. We will have

[TABLE]

Following the same steps with section 3.2 we have for the entropy of $p(\mathbf{x,y}|D^{\prime})$:

[TABLE]

We focus on computing the first term on the right hand side. In this case of a priori unknown variance the conditional output follows the $t-$student distribution (eq. (48)). Using standard expressions for its entropy we have (setting $N^{\prime}=\nu+N+1$)

[TABLE]

where,

[TABLE]

Note that in the second equality we used eq. (8). In general, we cannot compute analytically the entropy of the output, conditional on $D^{\prime}$. To this end, the mutual information of the input and output, conditioned on $D^{\prime},$ takes the form

[TABLE]

where $R$ are terms that do not depend on the new point $\mathbf{h}$. The second and third terms are computed for each $\mathbf{h}$ with direct Monte-Carlo using $10^{5}$ samples. It is important to emphasize that the mutual information criterion with unknown variance, $(\ref{info_full})$, depends on the output values $\mathbf{Y}$ (through the third term on the right hand side), in contrary to the known variance case (section 3.2). However, as it can be seen from eq. (51) this dependence can be very weak or even zero depending on the value of the parameter $\alpha$ which is chosen based on maximization of the evidence.

In the last expression the term with the highest computational cost is the entropy of the output (first term) as one needs to estimate the histogram of $y$. An analytical approximation can be obtained based on a Gaussian assumption for the output, $y$. In this case the mutual information takes the form:

[TABLE]

where $\sigma_{y}^{2}(\mathbf{h})$ is the estimated variance of the output after a new candidate input $\mathbf{h}$. This is computed with Monte-Carlo simulation of $y$, i.e. generate $10^{5}$ random realizations using (47), (3) and (44).

$Q-$criterion with unknown output variance

We emphasize that the selection approach using the $Q$ criterion is not modified at all for the case of unknown output variance. This is because the unknown variance $\sigma_{V}^{2}$ appears as a multiplication factor in the $Q$ criterion (eq. (31)), i.e. $\sigma_{V}^{2}$ is re-estimated each time a new data point is added using equation (31) but its value does not modify the optimal sample $\mathbf{h}$.

Numerical comparison for the 2d linear problem

Results and direct comparison with the case of known $\mathbf{\sigma}_{V}^{2}$ are shown in Figure 13 for the linear problem of section 5.1 with a two-dimensional input. For all methods shown we have run $400$ experiments (as we did for Fig. 4). In the higher dimensional problem the approach based on direct computation of mutual information is not applicable due to the vast computational cost. As we can observe the selection process based on mutual information with unknown output variance ($\mathcal{\hat{I}}_{G}$ or $\mathcal{\hat{I}}$), has slightly improved performance compared with the case of known $\sigma_{V}^{2}$ ($\mathcal{I}_{G}$ or $\mathcal{I}$) but this improvement is observed only for very small number of samples. As $N$ increases the criterion with unknown output variance is comparable with mutual information with given output variance ($\mathcal{I}_{G}$ or $\mathcal{I}$). This is because the optimal value of $\alpha$ tends to zero as $N$ increases and therefore the criteria with known and unknown output variance become practically identical.

It is important to emphasize that beyond the faster convergence of the $Q-$criterion, its main advantage is the low computational cost (many orders of magnitude smaller compared with the methods based on mutual information either with direct computation or through a Gaussian approximation) which allows applicability to higher dimensional problems.