Observable dictionary learning for high-dimensional statistical inference

TL;DR

This paper presents a novel observable dictionary learning method for high-dimensional inference, enabling efficient and accurate estimation of quantities from limited observations using a Bayesian framework and tailored dictionaries.

Contribution

It introduces a new observable dictionary learning approach that aligns the basis with the observation chain for improved inference accuracy.

Findings

Outperforms PCA and K-SVD in velocity field estimation

Provides closed-form posterior distribution for Bayesian inference

Demonstrates effectiveness on cavity flow velocity estimation

Abstract

This paper introduces a method for efficiently inferring a high-dimensional distributed quantity from a few observations. The quantity of interest (QoI) is approximated in a basis (dictionary) learned from a training set. The coefficients associated with the approximation of the QoI in the basis are determined by minimizing the misfit with the observations. To obtain a probabilistic estimate of the quantity of interest, a Bayesian approach is employed. The QoI is treated as a random field endowed with a hierarchical prior distribution so that closed-form expressions can be obtained for the posterior distribution. The main contribution of the present work lies in the derivation of \emph{a representation basis consistent with the observation chain} used to infer the associated coefficients. The resulting dictionary is then tailored to be both observable by the sensors and accurate in approximating the posterior mean. An algorithm for deriving such an observable dictionary is presented. The method is illustrated with the estimation of the velocity field of an open cavity flow from a handful of wall-mounted point sensors. Comparison with standard estimation approaches relying on Principal Component Analysis and K-SVD dictionaries is provided and illustrates the superior performance of the present approach.

Equations

$$y\left(\tau \right)\approx y^{\left(n_D\right)}\left(\tau \right)=\sum _{l=1}^{n_D}{x}^l{d}_l\left(\tau \right).$$

$$\mathbf{s}=F\left(y^{\left(n_D\right)}\right)+\mathbf{\eta },$$

$${\Vert u\Vert }_{\mathcal{W}}^2:={⟨u,u⟩}_{\mathcal{W}}<+\infty ,{⟨u,v⟩}_{\mathcal{W}}:={\int }_{\Theta }u\left(\theta \right)v\left(\theta \right)\mathrm{d}\mu \left(\theta \right),\forall \left(u,v\right)\in L^2\left(\Theta ,\mu \right).$$

$$y\left(\tau ,\theta \right)=\overline{y}\left(\tau \right)+\sum _{l=1}^{\infty }{\lambda }_l^{1/2}{\phi }_l\left(\tau \right){\psi }_l\left(\theta \right),\overline{y}\left(\tau \right)\equiv {\mathcal{E}}_{\mu }\left[y\left(\tau ,\theta \right)\right]:={\int }_{\Theta }y\left(\tau ,\theta \right)\mathrm{d}\mu \left(\theta \right).$$

$${\kappa }_y\left(\tau ,{\tau }^{\prime }\right):={\mathcal{E}}_{\mu }\left[\left(y\left(\tau ,\theta \right)-\overline{y}\left(\tau \right)\right)\left(y\left({\tau }^{\prime },\theta \right)-\overline{y}\left({\tau }^{\prime }\right)\right)\right].$$

$${\Vert y-y^{\left(n_D\right)}\Vert }_{\mathcal{H}\otimes \mathcal{W}}=\underset{{\phi }_l\in \mathcal{H},{\psi }_l\in \mathcal{W}}{\min }{\Vert y-\overline{y}-\sum _{l=1}^{n_D}{\lambda }_l^{1/2}{\phi }_l{\psi }_l\Vert }_{\mathcal{H}\otimes \mathcal{W}}.$$

$$C_y:={\left(n_{\mathrm{train}}-1\right)}^{-1}Y{Y}^{\mathsf{T}}⪰0.$$

$${\Vert Y-Y^{\left(n_D\right)}\Vert }_F=\underset{{\mathbf{\phi }}_l\in {\mathcal{H}}_h,{\mathbf{\psi }}_l\in {\mathcal{W}}_h}{\min }{\Vert Y-\sum _{l=1}^{n_D}{\lambda }_l^{1/2}{\mathbf{\phi }}_l\otimes {\mathbf{\psi }}_l\Vert }_F.$$

$$\left(Y^{\mathsf{T}}Y\right)\Psi \stackrel{\mathrm{eigen}}{=}\Psi \Lambda ,D=Y\Psi {\Lambda }^{-1/2}=\left({\mathbf{d}}_1{\mathbf{d}}_2\dots {\mathbf{d}}_{n_D}\right).$$

$$\mathbf{y}\approx \hat{\mathbf{y}}=D\mathbf{x},\mathrm{with}\mathbf{x}\in \underset{\tilde{\mathbf{x}}\in {\mathbb{R}}^{n_D}}{\mathrm{argmin}}{\Vert \mathbf{s}-F\left(D\tilde{\mathbf{x}}\right)\Vert }_2^2.$$

$$\mathbf{x}\in \underset{\tilde{\mathbf{x}}\in {\mathbb{R}}^{n_D}}{\mathrm{argmin}}{\Vert \mathbf{s}-FD\tilde{\mathbf{x}}\Vert }_2^2.$$

$$\mathbf{y}\approx \hat{\mathbf{y}}=D\mathbf{x},\mathrm{with}\mathbf{x}={\left(FD\right)}^{+}\mathbf{s},$$

$$\mathbf{y}\approx \hat{\mathbf{y}}=D_{\mathcal{I}}{\mathbf{x}}^{\mathcal{I}},\mathrm{with}{\mathbf{x}}^{\mathcal{I}}={\left({\left(F{D}_{\mathcal{I}}\right)}^{\mathsf{T}}F{D}_{\mathcal{I}}\right)}^{-1}{\left(F{D}_{\mathcal{I}}\right)}^{\mathsf{T}}\mathbf{s}.$$

$$\mathbf{x}\in \underset{\tilde{\mathbf{x}}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert A\tilde{\mathbf{x}}-\mathbf{b}\Vert }_2^2+\rho {\Vert \tilde{\mathbf{x}}\Vert }_1,$$

$$\left\{D,X\right\}\in \underset{\tilde{D},\tilde{X}}{\mathrm{argmin}}{\Vert Y-\tilde{D}\tilde{X}\Vert }_F^{2}s.t.{\Vert {\mathbf{x}}_i\Vert }_0\le K,\forall 1\le i\le n_{\mathrm{train}},$$

$$Y\approx D{D}^{+}Y,$$

$$Y\approx D{\left(FD\right)}^{+}FY.$$

$$\left\{\Phi ,\psi \right\}$$

$$=\underset{\tilde{\Phi },\tilde{\psi }}{\mathrm{argmin}}{\mathcal{E}}_{\mu }\left[{\Vert y\left(\tau ,\theta \right)-\sum _{l=1}^{n_D}{\tilde{\phi }}_l\left(\tau \right){\tilde{\psi }}_l\left(\theta ;\mathbf{s}\left(\theta \right)\right)\Vert }_{\mathcal{H}}^2\right],$$

$${\mathcal{E}}_{\mu ,{\mu }_{\eta }}\left[{\Vert y\left(\tau ,\theta \right)-y^{\left(n_D\right)}\left(\tau ,\theta ;\mathbf{s}\left(\theta ,{\theta }_{\eta }\right)\right)\Vert }_{\mathcal{H}}^2\right].$$

$$\left\{D,X\right\}$$

$$X$$

$$D$$

$$\left\{P,X\right\}$$

$$D$$

$${\mathbf{x}}_i\in \underset{\tilde{\mathbf{x}}\in {\mathbb{R}}^{n_D}}{\mathrm{argmin}}{\Vert {\mathbf{s}}_i-P\tilde{\mathbf{x}}\Vert }_2^{2}s.t.{\Vert {\mathbf{x}}_i\Vert }_0\le n_s,\forall 1\le i\le n_{\mathrm{train}}.[\mathbf{Sparse}\mathbf{Coding}]$$

$$D=Y{X}^{+}.[\mathbf{Estimation}\mathbf{Codebook}\mathbf{Update}]$$

$$E^{\left(Y,l\right)}\stackrel{\mathrm{SVD}}{=}\Phi \Sigma {\Psi }^{\mathsf{H}},{\tilde{\mathbf{x}}}_{{\mathcal{I}}^{(l)}}^l\propto {\Psi }^{\mathsf{H},1}.$$

$${\mathbf{p}}_l\propto E^{\left(S,l\right)}{\left({\tilde{\mathbf{x}}}_{{\mathcal{I}}^{(l)}}^l\right)}^{\mathsf{T}},{\Vert {\mathbf{p}}_l\Vert }_2=1.[\mathbf{Features}\mathbf{Codebook}\mathbf{Update}]$$

$$B^{\left(D\right)}\in \underset{{\tilde{B}}^{\left(D\right)}\in {\mathbb{R}}^{n_{\mathrm{train}}\times n_D}}{\mathrm{argmin}}{\Vert Q^{\left(Y\right)}{R}^{\left(Y\right)}-Q^{\left(Y\right)}{R}^{\left(Y\right)}{\tilde{B}}^{\left(D\right)}X\Vert }_F^2,D=Y{B}^{\left(D\right)},$$

$$C^{\left(D\right)}\in \underset{{\tilde{C}}^{\left(D\right)}\in {\mathbb{R}}^{n_{\mathrm{train}}\times n_D}}{\mathrm{argmin}}{\Vert R^{\left(Y\right)}-{\tilde{C}}^{\left(D\right)}X\Vert }_F^2.$$

$$Y^{\mathsf{T}}Y\stackrel{\mathrm{eigen}}{=}{V}^{(Y)}{\left({\Sigma }^{(Y)}\right)}^2{\left(V^{(Y)}\right)}^{\mathsf{T}},$$

Full text

\savesymbol

AND

1 Introduction

State estimation and parameter inference are very common problems occurring in many situations. As an example, when controlling a physical system, state feedback control achieves the best results but the state vector is rarely fully mesurable and it is hence necessary to estimate it from the scarce observations (sensors) available. In another context, numerical simulations in general have now reached a degree of maturity and reliability such that the bottleneck of the accuracy of the output of the simulations typically is the quality of the input parameters. These parameters can be initial conditions, physical properties, geometry, etc. Accurately estimating these input parameters is key to reliable simulations. In both the situations discussed above, the estimation of the hidden variables typically has to be done from scarce, indirect, observations of the system at hand while the information of interest is often a high-, or even infinite-, dimensional quantity. This leads to mathematically severely ill-posed problems.

A uniform grid-based description of the quantity of interest (QoI) is often employed and yields a representation under the form of a (very) high-dimensional vector, expressing a high degree of redundancy and suffering from over-parameterization. Identification of such a grid-based description through inverse techniques is a challenge. Several approaches can be followed to address this issue. A popular one is to reduce the number of degrees of freedom to describe the quantity of interest, at the price of an approximation. If expertise knowledge is available, it can be exploited to guide the choice of the approximation space. For instance, physical quantities have finite-energy, motivating an approximation in, say, a $L^{2}$-space. Other typical examples include smoothness, known class of probability distribution, etc. This naive approach can advantageously be improved by an approach exploiting correlations of the QoI among the entries of its grid-based description. For instance, in the case the QoI is a spatial field of physical properties, as encountered, for example, in subsurface reconstruction problems, this amounts to exploiting spatial correlations and the fact that the field to be inferred typically enjoys some degree of smoothness.111Note that for physical quantities exhibiting sharp discontinuities, dedicated approaches have been proposed relying on zonation or multi-resolution, see for instance Ellam et al. (2016), or multi-point statistics, Strebelle (2002). A representation basis is chosen and characterizing the QoI reduces to estimating the coefficients associated with the elements of the basis. This certainly is the most popular approach to inverse problems where the estimation typically reformulates as finding the small set of scalar-valued coefficients associated with the dominant modes of the approximation basis. Among many efforts from the literature, this strategy has been adopted in Buffoni et al. (2008); Yang et al. (2010) or Leroux et al. (2015), where a Principal Component Analysis (PCA), also termed Proper Orthogonal Decomposition (POD), Holmes et al. (2012); Sirovich (1987), was used as an approximation basis. Once the approximation basis is determined, estimating the associated coefficients is done by minimizing the data misfit with the available observations. With very few observations, this data misfit minimization can itself be ill-posed and the solution method must involve some kind of regularization. A typical choice, often justified from physics-related arguments, is to look for the least energy solution, i.e., that with the lowest $L^{2}$-norm: the ubiquitous Tikhonov regularization. Another approach is to determine an approximation basis in which the QoI is likely to be sparse. The data misfit minimization then includes a penalization of the cardinality of the solution, in the form of a $L^{0}$-(quasi-)norm, and is efficiently solved using sparse coding algorithms such as Orthogonal Matching Pursuit (OMP), Mallat & Zhang (1993), CoSaMP, Needell & Tropp (2008), or Iterative Hard Thresholding (IHT), Blumensath & Davies (2009), to cite only a few. For instance, this approach is followed by Khaninezhad et al. (2012a, b) who have used a K-SVD algorithm, Aharon et al. (2006); Rubinstein et al. (2008), to determine such a sparsity-promoting approximation basis in a subsurface reconstruction context.

Another strategy to address the ill-posedness of inverse problems is to explicitly introduce prior knowledge. For instance, a deterministic QoI can be treated as a random field and endowed with a prior distribution. A Bayesian approach can then be employed to “assimilate” the QoI from the prior to the posterior distribution using the available observations. This is a rigorous approach and results in a complete probabilistic estimation of the parameters at hand, Evans & Stark (2002); Kaipio & Somersalo (2005); Marzouk & Najm (2009); Stuart (2010). The QoI is then described in a probabilistic sense, with closed-form expressions in some specific situations (linear forward model from the input to the observations, Gaussian distribution of the likelihood, etc.) or, more generally, by sampling the input space according to the posterior measure density with a Markov-chain Monte-Carlo (MCMC) technique, Brooks et al. (2013).

Note that the two strategies can be combined and most efficient methods use both an approximation basis for the input space and a Bayesian approach. As an example, a Karhunen-Loève-based approximation of the input space, subsequently sampled with a MCMC scheme, can be used, Efendiev et al. (2006); Marzouk & Najm (2009).

Besides the quality of the approximation basis or the relevance of the expert knowledge introduced as a prior in the estimation method, the performance of the recovery procedure critically depends on observations, both quantity and quality of which are key factors for an accurate inference. Among other things, including the type of acquired information, the sensors location is a critical point. In an operational situation, sensors are limited in number and are subjected to a series of constraints. For instance, in a fluid mechanics context, their are usually mounted on a solid wall, as opposed to measuring in the bulk flow. Optimal location of the sensors requires an exhaustive search and constitutes an NP-hard problem. Many approaches have been proposed to determine a reasonably good solution in a computationally acceptable way. These include Effective Independence, Kammer (1991), Optimal Driving Point, Sensor Set Expansion, Kinetic Energy Method and Variance Method, Bakir (2011); Meo & Zumpano (2005), the maximization of the smallest singular value of the Fisher Information Matrix (FIM), Alonso et al. (2004), compressed-sensing-inspired approaches, Bright et al. (2013); Sargsyan et al. (2015); Brunton et al. (2016), framesense, Ranieri et al. (2014) or the sensorspace method introduced by the authors of the present work, Kasper et al. (2015).

The quality of the reconstruction from the sensor information largely depends on the observability of the QoI from the sensors. Clearly, if the trace of the QoI on the support of the point sensors vanishes, the observations do not bear information from the QoI and the estimation is not possible. In the case of linear dynamical systems, observability of the state vector from the sensors is quantified by studying the spectrum of the observability Gramian matrix. A similar situation occurs when the QoI is estimated within an approximation basis. In this case, the quality of the estimation of the approximation coefficients largely depends on the observability of the elements of the basis from the sensors. It may happen that the sensors are well informed by some approximation modes while being blind to others. As will be seen below in our examples, this situation is quite common and strongly impacts the recovery performance.

The motivation for the present work is the following: we advocate that an efficient estimation procedure should make use of prior expertise allowing an educated choice of a suitable approximation space whose basis functions should be consistent with the way the associated coefficients will be estimated from the sensor measurements. As will be shown below, deriving an approximation basis disregarding the way the coefficients are actually determined, as is common practice with the use of “generic” bases such as PCA/POD or K-SVD, may significantly affect the performance of the recovery method.

The paper is structured as follows. The general framework of inverse problems is introduced in Sec. 2, together with the main notations of the paper. Approximation of both an infinite- and finite-dimensional quantity is introduced and solution methods to estimate the coefficients associated with the approximation basis from observations are discussed. The need for an approximation basis aware of the sensors subsequently used to inform the coefficients is advocated and our observability-constrained basis learning approach is presented in Sec. 3. Both a deterministic and a stochastic estimation methods are presented, leading to a deterministic estimate or to a full probabilistic description of the QoI when considered as a random field. The proposed basis learning method allows both cases. Computational aspects are briefly discussed and an algorithm to learn the approximation basis is presented, together with the full solution method. Our approach is illustrated with the flow field estimation from a few wall-mounted point sensors in a configuration of numerically simulated two-dimensional flow over an open cavity, Sec. 4. Robustness of the recovery is studied with respect to noise in the measurements. Comparisons with standard recovery techniques relying on both PCA/POD and K-SVD show the superior performance of our approach. The paper concludes with closing remarks in Sec. 5.

2 Inverse problem

2.1 Problem statement

Let us introduce the problem and associated notations. In inverse problems, one typically wants to estimate a real-valued quantity of interest $y\in\mathcal{V}_{y}$ from a set of information $\mathcal{S}_{s}$. Few pieces of information are typically available so that the $n_{s}$ real-valued elements of $\mathcal{S}_{s}$ are collected in a vector $\boldsymbol{s}\in\mathcal{V}_{\boldsymbol{s}}\subset\mathbb{R}^{n_{s}}$ of observations. The QoI can be a set of variables or a distributed quantity such as a field. The dimension of $\mathcal{V}_{y}$ can be finite or infinite and we note ${n_{y}}:=\operatorname*{{\mathrm{dim}}}\left(\mathcal{V}_{y}\right)$ where it is understood that ${n_{y}}$ can $\to+\infty$ in which case $\mathcal{V}_{y}$ is assumed to be a Hilbert space $\mathcal{H}$ of real-valued functions defined over a domain $\Omega$, $\mathcal{H}\equiv L^{2}\left(\Omega;\mathbb{R}\right)$, and endowed with an inner product $\left<{\cdot},{\cdot}\right>_{\mathcal{H}}$.

Observations are typically scarce while the representation of the field of interest involves a large number of degrees of freedom (DOFs), so that $n_{s}\ll{n_{y}}$.

In many situations, a good approximation of $y$ is sufficient for the application of interest. Further, the structure of the QoI along the coordinates $\tau\in\Omega$, the field $y\left(\tau\right)$ is indexed upon, can often be exploited to lower the number of DOFs necessary for a good description. The field of interest can be approximated in a low dimensional linear subspace $\mathcal{V}_{{D}}$ on the Grassmann manifold $\mathrm{Gr}\left(n_{{D}},\mathcal{V}_{y}\right)$: $\displaystyle y\approx{{y}^{\left(n_{{D}}\right)}}\in\mathcal{V}_{{D}}\subset\mathbb{R}^{n_{y}}$. Let $\left\{d_{l}\right\}_{l}$ be a set of modes defining a basis of $\mathcal{V}_{{D}}$, the field of interest is approximated as

[TABLE]

The set ${D}:=\left\{d_{1},d_{2},\ldots d_{n_{{D}}}\right\}$ gathers the approximation modes and is termed the dictionary. Under this approximation, the field of interest is completely characterized in $\mathcal{V}_{{D}}$ by the set of coefficients $\boldsymbol{x}\in\mathbb{R}^{n_{{D}}}$. The inverse problem then reformulates as, given a dictionary ${D}$, inferring $\boldsymbol{x}$ from $\boldsymbol{s}$.

For the inference problem to admit a solution, observations $\boldsymbol{s}$ must be informed by the field of interest. Introducing the observation form $F^{\infty}:y\in\mathcal{V}_{y}\mapsto\boldsymbol{s}\in\mathcal{V}_{\boldsymbol{s}}$, observations are related to $y$ via $\boldsymbol{s}=F^{\infty}\left(y\right)$. In practice, as discussed above, the QoI is approximated in a low dimensional space $\mathcal{V}_{{D}}$. Let $F$ be a finite-dimensional continuous operator, $F:{{y}^{\left(n_{{D}}\right)}}\in\mathcal{V}_{{D}}\mapsto\boldsymbol{s}\in\mathcal{V}_{\boldsymbol{s}}$, the observations obey

[TABLE]

where $\boldsymbol{\eta}\in\mathcal{V}_{\boldsymbol{s}}$ denotes the contribution of $y^{\backslash\mathcal{V}_{{D}}}\in\mathcal{V}_{y}\backslash\mathcal{V}_{{D}}$ through $F^{\infty}$ and measurement noise, if any.

2.2 Approximation basis

To allow for a reduction of the dimensionality of the quantity to infer, an approximation basis must be chosen, see Eq. (1). It should allow for an accurate approximation of $y$ while being low dimensional. Prior or expertise knowledge on the QoI is typically used to guide the derivation of a suitable basis. In the context of an unknown spatially-distributed quantity over a spatial domain $\Omega$, the QoI is typically considered as a statistical realization of a random field. Consider the probability space $\left(\Theta,\sigma,\mu\right)$, where $\Theta$ is the set of random elementary events $\theta$, $\sigma$ the $\sigma$-algebra of the events and $\mu$ a probability measure. The space $\mathcal{W}\equiv L^{2}\left(\Theta,\mu\right)$ of real-valued, second-order, random variables with inner product $\left<{\cdot},{\cdot}\right>_{\mathcal{W}}$ is such that

[TABLE]

If the QoI is a realization $y\left(\tau;\theta\right)$ of a second-order random field, it can be efficiently approximated in $\mathcal{H}\otimes\mathcal{W}$ using the Hilbert-Karhunen-Loève decomposition, Levy & Rubinstein (1999):

[TABLE]

The random variables $\left\{\psi_{l}\right\}_{l}$ are uncorrelated in the sense of $\left<{\cdot},{\cdot}\right>_{\mathcal{W}}$, $\left<{\psi_{l}},{\psi_{l^{\prime}}}\right>_{\mathcal{W}}=0$ if $l\neq l^{\prime}$, and $\left(\lambda_{l},\varphi_{l}\right)\in\mathbb{R}^{+}\times\mathcal{H}$ is the ordered (decreasing $\lambda_{l}$) sequence of eigenvalues and eigenfunctions of the linear operator $\displaystyle\mathscr{T}_{y}\left(u\right):\mathcal{H}\to\mathcal{H}:=\left<{\kappa_{y}},{u}\right>_{\mathcal{H}}=\mathscr{E}_{\mu}\left[{\left(y-\overline{y}\right)\,\left<{y-\overline{y}},{u}\right>_{\mathcal{H}}}\right]$. The correlation kernel $\kappa_{y}\in\mathcal{H}\otimes\mathcal{H}$ is defined as

[TABLE]

For a known kernel $\kappa$, a popular choice of approximation basis of a realization in $\mathcal{H}$ is the set of eigenfunctions $\left\{\varphi_{l}\right\}_{l}$. This is justified by the classic optimality result: letting $\displaystyle{{y}^{\left(n_{{D}}\right)}}\left(\tau,\theta\right):=\overline{y}\left(\tau\right)+\sum_{l=1}^{n_{{D}}}{\lambda^{1/2}_{l}\,\varphi_{l}\left(\tau\right)\,\psi_{l}\left(\theta\right)}$ be the $n_{{D}}$-term truncation of the series in Eq. (4) and $\displaystyle\left\|\cdot\right\|_{\mathcal{H}\otimes\mathcal{W}}$ the norm induced by the scalar product $\displaystyle\left<{\cdot},{\cdot}\right>_{\mathcal{H}\otimes\mathcal{W}}:=\mathscr{E}_{\mu}\left[{\left<{\cdot},{\cdot}\right>_{\mathcal{H}}}\right]$, one infers

[TABLE]

For a given number of terms retained, the approximation error is then of minimal norm, in the sense of Eq. (6). Approximating a given instance $y\left(\tau;\theta\right)$ of the random field in the linear $\operatorname*{{\mathrm{span}}}\left\{\varphi_{l}\right\}_{l}$ hence minimizes the average spatial error norm $\left\|{y-{{y}^{\left(n_{{D}}\right)}}}\right\|_{\mathcal{H}}$ in the sense of $\mu$.

In practice, the correlation kernel is rarely known. However, if realizations of the random field are available, the correlation kernel can be approximated by the empirical correlation. These available realizations are usually finite-dimensional, as produced by numerical simulations or from past inference problems. A realization $y\left(\tau;\theta\right)\in\mathcal{V}_{y}$ can then be identified with $\boldsymbol{y}\in\mathbb{R}^{{n_{y}}}$ and the empirical covariance is determined from a collection of $n_{{\mathrm{train}}}$ known instances of the field $\left\{\boldsymbol{y}_{i}\right\}_{i=1}^{n_{{\mathrm{train}}}}$. Letting $\displaystyle{Y}:=\left(\boldsymbol{y}^{\circ}_{1}\>\boldsymbol{y}^{\circ}_{2}\ldots\boldsymbol{y}^{\circ}_{n_{{\mathrm{train}}}}\right)\in\mathbb{R}^{{n_{y}}\times n_{{\mathrm{train}}}}$ be the training matrix of centered known realizations $\boldsymbol{y}^{\circ}_{i}:=\boldsymbol{y}_{i}-\overline{\boldsymbol{y}}$, $\displaystyle\overline{\boldsymbol{y}}:=\left(n_{{\mathrm{train}}}\right)^{-1}\Sigma_{i=1}^{n_{{\mathrm{train}}}}{\boldsymbol{y}_{i}}$, the empirical covariance $C_{y}\in\mathbb{R}^{{n_{y}}\times{n_{y}}}$ is given by

[TABLE]

From the eigenvalues decomposition of the symmetric positive definite matrix $C_{y}\stackrel{{\scriptstyle\rm eigen}}{{=}}\varPhi\>\varLambda\>\varPhi^{\mathsf{T}}$, the approximation basis can be defined by the dictionary of the eigenvectors of $C_{y}$: ${D}=\varPhi=\left(\boldsymbol{\varphi}_{1}\,\boldsymbol{\varphi}_{2}\ldots\boldsymbol{\varphi}_{{n_{y}}}\right)$ and the training matrix can then express as ${Y}=\varPhi\,\varLambda^{1/2}\,\varPsi^{\mathsf{T}}$ with a suitable matrix $\varPsi\in\mathbb{R}^{n_{{\mathrm{train}}}\times{n_{y}}}$. However, one typically has a limited training set while the realizations are finely discretized so that $n_{{\mathrm{train}}}<{n_{y}}$. The matrix $C_{y}$ is hence not full-rank and leads to a poor approximation of the true covariance. Further, it involves the dimension ${n_{y}}$ of the discretized field, which typically is large, leading to a significant computational burden.

From now on, and without loss of generality, the data are assumed centered. With a collection of finite-dimensional (centered) training instances $\left\{\boldsymbol{y}_{i}\right\}_{i}$, and in the context of a second-order QoI, $y\in L^{2}\left(\mathbb{R}^{{n_{y}}}\right)\otimes L^{2}\left(\mathbb{R}^{n_{{\mathrm{train}}}}\right)\equiv\mathcal{H}_{h}\otimes\mathcal{W}_{h}$, $\mathcal{H}_{h}\subset\mathcal{V}_{y}$, the expression in Eq. (6) is reformulated in

[TABLE]

The Eckart-Young theorem, Eckart & Young (1936), allows to directly identify $\left\{\boldsymbol{\varphi}_{l}\right\}_{l}$ with the $n_{{D}}$ dominant left singular vectors of ${Y}$ and $\left\{\boldsymbol{\psi}_{l}\right\}_{l}$ and $\left\{\lambda_{l}\right\}_{l}$ with the corresponding right singular vectors and singular values. Letting $\varSigma=\varLambda^{1/2}$ be the diagonal $n_{{D}}\times n_{{D}}$ matrix of the positive square root of (positive) eigenvalues $\left\{\lambda_{l}\right\}_{l}$, the training matrix is best approximated with a rank-$n_{{D}}$ matrix as ${Y}\stackrel{{\scriptstyle{\rm rank}\mathchar 45\relax{n_{{D}}}}}{{\approx}}\varPhi\,\varSigma\,\varPsi^{\mathsf{H}}$. The dictionary ${D}$ is then given by the dominant left singular vectors and is efficiently determined from the solution $\varPsi$ of the following $n_{{\mathrm{train}}}$-dimensional eigen problem:

[TABLE]

The resulting approximation modes $\left\{\boldsymbol{d}_{l}\right\}_{l}$ are sometimes termed principal components (PCA) or empirical modes in the context of the Proper Orthogonal Decomposition (POD), Holmes et al. (2012).

2.3 Direct approach

Once an approximation basis for the QoI is determined, the problem becomes to identify the coefficients associated with the approximation, $\displaystyle\boldsymbol{y}\approx{D}\,\boldsymbol{x}$, from the only pieces of available information $\boldsymbol{s}$. The link between measurements $\boldsymbol{s}$ and the QoI $\boldsymbol{y}$ is similar as in Eq. (2) and the estimation $\widehat{\boldsymbol{y}}$ can formulate as minimizing the $L^{2}$-norm of the data misfit:

[TABLE]

In the present study, $F:\mathcal{H}_{h}\to L^{2}\left(\mathbb{R}^{n_{s}}\right)$ is a bounded finite-dimensional linear operator and the measurements have finite energy, $\boldsymbol{s}\in L^{2}\left(\mathbb{R}^{n_{s}}\right)\equiv\mathcal{S}_{h}\subset\mathcal{V}_{\boldsymbol{s}}$. Minimizing the data misfit then reduces to

[TABLE]

Since one has limited measurements compared to the size of the dictionary, $n_{s}\leq n_{{D}}$, the linear system of normal equations associated with Eq. (11) is underdetermined. Additional knowledge or constraints must hence be considered to select among the infinite number of solutions. A popular choice is to retain the minimum $L^{2}$-norm solution. This is partly justified by the heuristic that one looks for the QoI of minimal energy compatible with the observations. The inferred field is then given by

[TABLE]

where the superscript + denotes the Moore-Penrose pseudo-inverse. The lift-up operator $L:\mathcal{S}_{h}\to\mathcal{H}_{h}$ from measurements to the estimated QoI is then $L={D}\,\left(F\,{D}\right)^{+}$ and $\widehat{\boldsymbol{y}}=L\,\boldsymbol{s}$.

2.4 Parsimonious approach

An alternative choice among the infinite number of solutions to problem (11) is to select a solution involving at most $n_{s}$ distinct non-zero coefficients. Upon knowledge of the corresponding set of indices $\mathcal{I}\subseteq\left\{1,2,\ldots,n_{{D}}\right\}\in\mathbb{N}^{\ast}$, $|{\mathcal{I}}|\leq n_{s}$, one can define the coefficient vector restricted to these non-zero coefficients: $\boldsymbol{x}^{\mathcal{I}}\in\mathbb{R}^{|{\mathcal{I}}|}$. Let ${D}_{\mathcal{I}}\in\mathbb{R}^{{n_{y}}\times|{\mathcal{I}}|}$ be the restriction of the dictionary to its columns indexed by $\mathcal{I}$, the inferred field is given by the least squares solution

[TABLE]

This approach results in a well-posed problem for $\boldsymbol{x}$ once $\mathcal{I}$ is determined.

Determining a sparse solution to a linear system of equations has attracted significant efforts in the recent past. A decisive step was made with the Compressed Sensing theory, Candès & Tao (2004); Candès et al. (2005); Donoho (2006); Candès et al. (2006), which has rigorously justified the formulation of an underdetermined linear problem $A\,\boldsymbol{x}\approx\boldsymbol{b}$, $A$ being a tall $n\times p$ matrix, $n<p$, with a sparse solution as a constrained convex optimization problem (basis pursuit denoising problem, BPDN):

[TABLE]

where $\rho>0$ is a regularization parameter.

In case the cardinality of $\boldsymbol{x}$ is known or imposed, dedicated algorithms relying on a greedy approach, such as Orthogonal Matching Pursuit (OMP), Mallat & Zhang (1993), CoSaMP, Needell & Tropp (2008) or Iterative Hard Thresholding (IHT), Blumensath & Davies (2009), can be employed.

2.5 Dictionary learning

As discussed in the previous section, if the solution to a linear system is sparse, efficient solution methods exist. However, the approximation of the field of interest $\boldsymbol{y}$ with a sparse coefficient vector $\boldsymbol{x}$ associated with the dictionary ${D}$ is poor, in the general case. As an example, the field of interest $\boldsymbol{y}$ is unlikely to admit a sparse representation in the PCA basis. The benefit of a well-posed problem for estimating the set of (sparse) coefficients is then lost by the poor resulting approximation. This situation asks for a dictionary specifically tailored so that the associated coefficients are more likely to be sparse. That is, instead of using a generic approximation basis such as PCA, one should determine the basis222Strictly speaking, this is then not a basis since the elements $\left\{\boldsymbol{d}_{l}\right\}_{l}$ are not linearly independent and constitute an overcomplete set. with the specific aim of favoring sparse solutions. Many algorithms exist to this aim and we here briefly describe the K-SVD method, Aharon et al. (2006); Rubinstein et al. (2008), as it will be useful hereafter in the paper. For a given training set ${Y}\in\mathbb{R}^{{n_{y}}\times n_{{\mathrm{train}}}}$, the dictionary learning problem consists in finding ${D}\in\mathbb{R}^{{n_{y}}\times n_{{D}}}$ and ${X}\in\mathbb{R}^{n_{{D}}\times n_{{\mathrm{train}}}}$ such that

[TABLE]

where $K\in\mathbb{N}^{\ast}$ is the maximum sparsity of the approximation for any $\boldsymbol{y}_{i}$ from the training set.

The K-SVD algorithm essentially relies on an alternating minimization strategy and consists of two steps. In a Sparse Coding (SC) step, the matrix of coefficients ${X}$ is updated solving Eq. (15) given ${D}$. While different methods can be used for this step, the original implementation of K-SVD relies on the OMP algorithm. By construction, each column vector of ${X}$ is, at most, $K$-sparse. Once the coefficients are estimated, the dictionary ${D}$ is updated in a Codebook Update (CU) step, exploiting the sparse structure of ${X}$ and one iterates between these two steps until a convergence criterion is satisfied, see Aharon et al. (2006) for details.

2.6 Discussion

Two approximation bases have been briefly discussed above. If the correlation kernel of the QoI is known, or the empirical covariance of a representative training set can be estimated, the (Hilbert-) Karhunen-Loève decomposition can be used and leads to an approximation which is optimal in the sense of Eq. (8). Minimization of the data misfit allows to estimate the associated vector of coefficients $\boldsymbol{x}$ from the measurements. This estimation method is hereafter referred to as the Principal Component Analysis-based (PCA) method. A major limitation of this approach is that, without additional hypotheses on the solution, the estimation problem is only well-posed for a coefficient vector of size $n_{{D}}\leq n_{s}$. Since measurements are typically scarce, this can severely limit the versatility of the approximation basis and the accuracy of this approach.

The additional sparsity constraint introduced in the K-SVD approach, and its variants, allows to use a redundant dictionary in the sense that $\left\{\boldsymbol{d}_{l}\right\}_{l}$ constitutes a frame for $\mathcal{H}_{h}$. This redundancy property is key as it allows to rely on dictionary with more modes than measurements while still dealing with a solvable problem for $\boldsymbol{x}$ thanks to the suitable sparsity constraint that one should have $K\leq n_{s}\leq n_{{D}}$. While a K-SVD-based dictionary is not optimal at approximating elements of the training set $\left\{\boldsymbol{y}_{i}\right\}_{i}$ for a given size $n_{{D}}$, it however potentially allows to represent a QoI $\boldsymbol{y}_{i}$ with elements from a larger dictionary, $n_{{D}}\geq n_{s}$ than with the PCA approach where, in practice, one should have $n_{{D}}\leq n_{s}$. A potentially better approximation can hence be achieved.

3 Observability-constrained dictionary learning

3.1 Observability issue

While enjoying optimality in a different sense, dictionaries derived with the PCA or the K-SVD approaches discussed in the previous section, and more generally most dictionary learning techniques, rely on a minimization of some norm of the approximation residual of the (centered) training set, $\displaystyle\left\|{Y}-{D}\,{X}\right\|$. Once the dictionary is determined, the unknown QoI $\boldsymbol{y}$ is estimated as $\widehat{\boldsymbol{y}}={D}\,\boldsymbol{x}$, with $\boldsymbol{x}$ inferred from the measurements $\boldsymbol{s}$ by minimization of the data misfit.

This is the standard approach but it nonetheless suffers an important drawback. Indeed, the dictionary is determined a priori from the training fields $\left\{\boldsymbol{y}_{i}\right\}_{i}$ solely, while the actual QoI is, in fine, estimated from the measurements $\boldsymbol{s}$: $\boldsymbol{y}\approx\widehat{\boldsymbol{y}}={D}\,\boldsymbol{x}\left(\boldsymbol{s}\right)$ where the dependance of $\boldsymbol{x}$ upon $\boldsymbol{s}$ has been explicitly written. Instead, standard dictionary learning approaches rely on a training step where a given field is implicitly approximated as $\boldsymbol{y}\approx\widehat{\boldsymbol{y}}={D}\,\boldsymbol{x}\left(\boldsymbol{y}\right)$. For instance, in the CU step of K-SVD, modes $\left\{\boldsymbol{d}_{l}\right\}_{l}$ are updated given exact coefficient vectors $\left\{\boldsymbol{x}_{i}\right\}_{i}$ as evaluated from the SC step. This is significantly different from the actual situation where the coefficient vector can only be evaluated from the data misfit.

From this situation, it results that the modes of the dictionary may be poorly informed by the available data in practice. The modes of the dictionary are determined without accounting for the observation operator $F$. In particular, the standard approach essentially consists in looking for ${D}$ such that ${Y}\approx{D}\,{X}$ (CU-like step) with ${X}={D}^{+}{Y}$ (SC-like step) and it can be interpreted as looking for a $n_{{D}}$-dimensional autoencoder ${D}$ such that

[TABLE]

where the low-rank approximation is understood in the sense of (6) or (8).

Instead, a more consistent approach is to look for a $n_{{D}}$-dimensional autoencoder ${D}$ such that ${Y}\approx{D}\,{X}$ (CU-like step) with ${X}=\left(F\,{D}\right)^{+}{S}$ (SC-like step), i.e.,

[TABLE]

As an illustration of the potential issues of standard approaches of Sec. 2.3 and 2.4, consider the situation where the observation operator $F$ is simply a restriction of the field of interest to a small set of points (point measurements). The vector of observations $\boldsymbol{s}$ is then a subset of $\boldsymbol{y}$ at indices $\mathcal{J}\subseteq\left\{1,2,\ldots,{n_{y}}\right\}\in\mathbb{N}^{\ast}$, $|{\mathcal{J}}|=n_{s}$: $\boldsymbol{s}=\boldsymbol{y}^{\mathcal{J}}$. Since the dictionary is learned from ${Y}$ only, it may happen that the modes vanish over the support of the sensors, so that $F\,{D}={D}^{\mathcal{J}}\approx\left[0\right]_{n_{s}\times n_{{D}}}$. In this situation, the coefficients $\boldsymbol{x}$ cannot be estimated from minimization of the data misfit in (11): the dictionary is not observable from the sensors.

3.2 Goal-oriented basis learning

To address the limitations of standard estimation techniques discussed above, we now introduce a Goal-Oriented BAsis Learning (GOBAL) solution method. In a general, possibly infinite-dimensional, setting, we look for a redundant dictionary $\varPhi\left(\tau\right)=\left(\varphi_{1}\,\varphi_{2}\ldots\varphi_{n_{{D}}}\right)$, and an estimation technique providing the function $\psi$ from the measurements, minimizing the associated Bayes risk:

[TABLE]

with $y$ the centered data, $\mathscr{E}_{\mu}\left[{y\left(\cdot,\theta\right)}\right]=0$.

Remark 1

This definition of the Bayes risk implicitly assumes perfect measurements $\boldsymbol{s}=F\,y$. In practice, noise affects the measurements so that this equation does not hold anymore and Eq. (2) should be accounted for instead. The definition of the Bayes risk above can be readily extended to account for a random noise with a given measure. Assuming the noise is statistically independent from the realizations $\theta$, one can introduce the probability space it belongs to, with an associated set of random elementary events $\theta_{\eta}$ and probability measure $\mu_{\eta}$. Accordingly extending the definition of the expectation operator, the Bayes risk is then defined as

[TABLE]

Restricting to the finite-dimensional framework, we let ${D}\equiv\varPhi=\left(\boldsymbol{\varphi}_{1}\,\boldsymbol{\varphi}_{2}\ldots\boldsymbol{\varphi}_{n_{{D}}}\right)$. From a collection of a finite number $n_{{\mathrm{train}}}$ of independent identically distributed (iid) realizations $\left\{\boldsymbol{y}_{i}\right\}_{i}$ and associated (noisy) measurements $\left\{\boldsymbol{s}_{i}\right\}_{i}$, problem (18) can reformulate as looking for the dictionary that minimizes the empirical Bayes risk, with the Bayes risk definition (19):

[TABLE]

where ${Y}=\left(\boldsymbol{y}_{1}\,\boldsymbol{y}_{2}\ldots\boldsymbol{y}_{n_{{\mathrm{train}}}}\right)$ and ${S}=\left(\boldsymbol{s}_{1}\,\boldsymbol{s}_{2}\ldots\boldsymbol{s}_{n_{{\mathrm{train}}}}\right)$. The dictionary learning problem then admits a formulation similar to those discussed in Sec. 2. The coefficient vector $\boldsymbol{x}$ associated with the estimation $\widehat{\boldsymbol{y}}={D}\,\boldsymbol{x}$ is determined by minimization of the data misfit. Since the dictionary is redundant, $n_{{D}}\geq n_{s}$, constraints have to be considered for the estimation problem to be well-posed. Similarly to the parsimonious approach of Sec. 2.4, a sparsity constraint on $\boldsymbol{x}$ is enforced and Eq. (20) finally reformulates as:

[TABLE]

The decomposition (22) does not impose the basis elements to be orthogonal, in contrast with other methods such as PCA. This brings more flexibility in deriving an efficient approximation basis and results in superior potential performance. The dictionary ${D}$ minimizes the training field error given the coefficients ${X}$ (CU step) while the coefficients are such that they minimize the data misfit with sparsity constraints for a given dictionary (SC step).

As such, this formulation would however lead to a poor recovery procedure. Indeed, in the general case, it is unlikely that $\boldsymbol{s}_{i}$ admits a good $n_{s}$-sparse approximation in the frame defined by $F\,{D}$. The naive approach of Eqs. (21) and (22) would hence lead to a large residual norm in the SC step and a poor overall performance. Instead, one could explicitly introduce a sparsifying predictor operator ${P}=\left(\boldsymbol{p}_{1}\,\boldsymbol{p}_{2}\ldots\boldsymbol{p}_{n_{{D}}}\right):\mathbb{R}^{n_{{D}}}\to\mathcal{S}_{h}\equiv L^{2}\left(\mathbb{R}^{n_{s}}\right)$ learned from the training set of measurements such that $\boldsymbol{s}_{i}$ admits a good $n_{s}$-sparse approximation in the linear span of $\left\{\boldsymbol{p}_{l}\right\}_{l}$, $\forall\,1\leq i\leq n_{{\mathrm{train}}}$. The recovery procedure then relies on a pair of dictionaries $\left\{{P},{D}\right\}$ learned from training data with

[TABLE]

The general procedure discussed above is hereafter termed the GOBAL method and is now discussed in more details.

3.3 Implementation

As seen above, a given ${n_{y}}$-dimensional field $\boldsymbol{y}_{i}$ is approximated in a $n_{{D}}$-dimensional space, $\boldsymbol{y}_{i}\approx{D}\,\boldsymbol{x}_{i}$, with $n_{{D}}\leq{n_{y}}$. The $n_{{D}}$-dimensional coefficient vector $\boldsymbol{x}_{i}$ is informed by the measurements $\boldsymbol{s}_{i}\in\mathbb{R}^{n_{s}}$, with $n_{{D}}\geq n_{s}$. For a well-posed problem, $\boldsymbol{x}_{i}$ is then imposed to be $n_{s}$-sparse, at most. For a given predictor ${P}$, it is then the solution of a sparse coding problem

[TABLE]

Since ${X}$ is entirely determined from ${S}$ for a given ${P}$, the solution to (24) immediately follows:

[TABLE]

While ${D}$ as given above minimizes the Frobenius norm of the residual ${Y}-{D}\,{X}$ among all matrices of rank $n_{{D}}\leq\min\left(n_{{\mathrm{train}}},{n_{y}}\right)$, the Bayes risk can still be large. One should then exploit the remaining DOFs offered by ${P}$ to improve upon this situation.

From (23), ${P}$ is seen to minimize $\left\|{{S}-{P}\,\widetilde{{X}}}\right\|_{F}^{2}$ for a given $\widetilde{{X}}$. We here adopt an approach similar to the Codebook Update step in K-SVD. For a given atom index $l\in\left[1,n_{{D}}\right]\in\mathbb{N}^{\ast}$, let $\mathcal{I}^{(l)}$ be the set of column indices of non-vanishing coefficients, $\mathcal{I}^{(l)}=\left\{i,~{}1\leq i\leq n_{{\mathrm{train}}};~{}\widetilde{{X}}^{l}_{i}\neq 0\right\}$. Let $E^{\left({Y},l\right)}:={Y}_{\mathcal{I}^{(l)}}-{D}_{{\backslash l}}\,\widetilde{{X}}_{\mathcal{I}^{(l)}}^{{\backslash l}}\in\mathbb{R}^{{n_{y}}\times|{\mathcal{I}^{(l)}}|}$. The vector $\widetilde{\boldsymbol{x}}^{l}_{\mathcal{I}^{(l)}}$ can then be obtained from the rank-1 approximation of $E^{\left({Y},l\right)}$:

[TABLE]

Iterating at random over the atoms $1\leq l\leq n_{{D}}$, one can estimate the coefficients $\widetilde{{X}}$ informed from the QoI recovery error $E^{\left({Y},l\right)}$.

As for the coefficients ${X}$ above, the predictor ${P}$, solution to (23), is then determined in the spirit of the K-SVD algorithm. For any atom index $l$, let the measurement error matrix be $E^{\left({S},l\right)}:={S}_{\mathcal{I}^{(l)}}-{{P}}_{{\backslash l}}\,{X}_{\mathcal{I}^{(l)}}^{{\backslash l}}\in\mathbb{R}^{n_{s}\times|{\mathcal{I}^{(l)}}|}$. The predictor $\boldsymbol{p}_{l}$ is then determined as the minimizer of $\left\|{E^{\left({S},l\right)}-\boldsymbol{p}_{l}\,\widetilde{\boldsymbol{x}}^{l}_{\mathcal{I}^{(l)}}}\right\|_{F}$, with the coefficient vector $\widetilde{\boldsymbol{x}}^{l}_{\mathcal{I}^{(l)}}$ achieving a low estimation error $\left\|{E^{\left({Y},l\right)}}\right\|_{F}$ as given by (27). It yields

[TABLE]

This procedure leads to a predictor ${P}$ minimizing the data misfit while accounting for associated coefficients ${X}$ achieving a low QoI estimation error.

3.4 Computational aspects

While the procedure described above solely relies on simple linear algebra operations, it involves potentially large scale vectors, typically the training elements $\boldsymbol{y}_{i}\in\mathbb{R}^{n_{y}}$. In order to improve its computational efficiency, a number of modifications to the plain procedure is now discussed.

The estimation of the dictionary ${D}$ is given by the product of the training field matrix ${Y}$ and the Moore-Penrose pseudo-inverse of the real-valued $n_{{D}}\times n_{{\mathrm{train}}}$ matrix of coefficients ${X}$, Eq. (26). Efficient algorithms are available to compute ${X}^{+}$. However, ${X}$, as given by the Sparse Coding step in (25), is sparse, i.e., each of its columns only has $n_{s}$ non-zero entries, at most. This structure can be exploited when computing its pseudo-inverse. In this work, we rely on the QRginv procedure introduced in Katsikis et al. (2011) which explicitly accounts for the sparse structure of ${X}$ to speed up computations.

In the same objective of handling large matrices, it is important to notice that ${D}$, as given by Eq. (26), lies in the linear span of the columns of ${Y}$. Each dictionary atom $\boldsymbol{d}_{l}$ then lies in the linear span of $\left\{\boldsymbol{y}_{i}\right\}_{i=1}^{n_{{\mathrm{train}}}}$. Let consider the QR decomposition of the training matrix, ${Y}\stackrel{{\scriptstyle\rm QR}}{{=}}Q^{\left({Y}\right)}\,R^{\left({Y}\right)}$. Since $\boldsymbol{d}_{l}\in\operatorname*{{\mathrm{span}}}\left\{\boldsymbol{y}_{i}\right\}_{i}$, there exists a matrix ${B^{\left({D}\right)}}$ such that ${D}={Y}\,{B^{\left({D}\right)}}$. The dictionary is the solution to Eq. (24), which reformulates as

[TABLE]

or, with $C^{\left({D}\right)}$ such that ${D}=Q^{\left({Y}\right)}\,C^{\left({D}\right)}$,

[TABLE]

Since $n_{{\mathrm{train}}}\leq{n_{y}}$, the residual in Eq. (31) involves smaller matrices than in Eq. (24) and is hence faster to evaluate. The dimension ${n_{y}}$ of $\boldsymbol{y}$ then entirely vanishes from the algorithm, allowing to handle QoIs of arbitrary size, provided their QR decomposition can be computed in a prior (offline) step.

Remark 2

The idea above is to decompose the, potentially very large, training matrix ${Y}$ in a unitary matrix $Q^{\left({Y}\right)}$ and a smaller $n_{{\mathrm{train}}}\times n_{{\mathrm{train}}}$ matrix $R^{\left({Y}\right)}$, so that the estimation procedure involves $R^{\left({Y}\right)}$ only. A potentially faster alternative to the QR decomposition discussed above relies on an eigenvalue decomposition of a small $n_{{\mathrm{train}}}\times n_{{\mathrm{train}}}$ symmetric matrix:

[TABLE]

so that there exists $R^{\left({Y}\right)}$ such that ${Y}=Q^{\left({Y}\right)}\,R^{\left({Y}\right)}$, with $Q^{\left({Y}\right)}$ a unitary matrix:

[TABLE]

with $\varSigma^{({Y})}$ a diagonal matrix with non-negative entries.

Disregarding the specific decomposition used, the generic terms $Q^{\left({Y}\right)}$ and $R^{\left({Y}\right)}$ will be employed in the remainder of the paper.

The last implementation aspect concerns the SVD in Eq. (27) as part of the Features Codebook Update step. The whole GOBAL procedure relies on a block-coordinate descent approach, with alternate minimizations between the three steps. As such, high accuracy in not necessarily required at early steps of the iterative procedure and the SVD can safely be approximated. Further, one should notice that only the dominant (right) singular vector $\varPsi^{\mathsf{H},1}$ is required. An inexpensive approximation of the solution to Eq. (27) is then given by a 1-iteration alternate minimization:

[TABLE]

Using the ${Y}=Q^{\left({Y}\right)}\,R^{\left({Y}\right)}$ decomposition discussed above, the approximation of the solution to (27) can be derived even more efficiently:

[TABLE]

From this coefficient vector $\widetilde{\boldsymbol{x}}^{l}_{\mathcal{I}^{(l)}}$ as estimated from the residual $E^{\left(R,l\right)}$, the predictor $\boldsymbol{p}_{l}$ follows from (28).

3.5 Sparse Bayesian Learning

3.5.1 A probabilistic estimate

The GOBAL method as discussed so far allows to learn a dictionary leading to a good sparse approximation of a set of data. From given observations $\boldsymbol{s}$, one gets a deterministic estimation $\widehat{\boldsymbol{y}}={D}\,\boldsymbol{x}\left(\boldsymbol{s}\right)$. Instead, it is often useful to know the uncertainty associated with an outcome, in a broad sense. It allows, for instance, to define confidence intervals for the estimated QoI and identify spatial regions where the estimation is likely to be poor, hence possibly asking for more local information, such as additional sensors. To extend the general framework of the GOBAL approach introduced in the previous section, we now adopt a Bayesian viewpoint.

In the general case where the sparsity of the solution is not known beforehand, the sparse coding step given the predictor operator ${P}$ can formulate under the form of the BPDN problem, see also Eq. (14), where the compressible coefficient vector is given by:

[TABLE]

Observations are assumed to be affected with additive errors $\boldsymbol{\eta}$ independent from the prediction:

[TABLE]

The components of the error are modeled as zero-mean iid Gaussian random variables. The measurement noise is assumed to affect similarly and independently each component of $\boldsymbol{s}$ and the error $\boldsymbol{\eta}$ is modeled as a multivariate, zero-mean, Gaussian random variable with diagonal covariance matrix ${\beta^{-1}\,I_{n_{s}}}$: $\boldsymbol{\eta}~{}\sim~{}\mathcal{N}\left(\boldsymbol{0},{\beta^{-1}\,I_{n_{s}}}\right)$.

The coefficient vector is modeled as a random vector endowed with a Laplace prior density:

[TABLE]

with $\lambda>0$. In a Bayesian perspective, the coefficient vector as estimated from the observations is described by its posterior density

[TABLE]

This Bayesian formulation provides a complete description of the coefficient vector. In particular, since $\displaystyle\operatorname*{{\mathrm{arg\,min}}}_{\widetilde{\boldsymbol{x}}}\,\left\|\boldsymbol{s}-{P}\,\widetilde{\boldsymbol{x}}\right\|_{2}^{2}+\rho\,\left\|\widetilde{\boldsymbol{x}}\right\|_{1}=\operatorname*{{\mathrm{arg\,max}}}_{\widetilde{\boldsymbol{x}}}\,\mathrm{p}_{\boldsymbol{s}}\left({\boldsymbol{s}\mid\widetilde{\boldsymbol{x}},\beta}\right)\,\mathrm{p}_{\boldsymbol{x}}\left({\widetilde{\boldsymbol{x}}\mid\lambda}\right)$, with $\lambda=\rho\,\beta$, the solution to the BPDN problem (36) is equivalently obtained as the maximum a posteriori (MAP).

3.5.2 Hierarchical prior

A priori values for the parameters $\beta$ and $\lambda$ are difficult to set and they should be learned from the data instead. We then adopt a hyperprior formulation and rely on an approach inspired from the Sparse Bayesian Learning (SBL) described in Tipping & Faul (2006) and Babacan et al. (2010). In particular, while the Laplace prior (38) is not conjugate to the noise model (37), a suitable three-stage hierarchical prior formulation allows to derive a closed-form expression of the MAP via a constructive sequential procedure.

To this end, the prior (38) is substituted with a product of univariate Gaussian random variables:

[TABLE]

with $\boldsymbol{\gamma}:=\left(\gamma_{1}\,\gamma_{2}\ldots\gamma_{n_{{D}}}\right)$ and the variances $\gamma_{l}\geq 0$, $1\leq l\leq n_{{D}}$, modeled as

[TABLE]

where $\lambda\geq 0$ has a Gamma prior with parameter $\nu>0$

[TABLE]

The hyperparameter $\lambda$ has the same role as in Eq. (38) since

[TABLE]

is of similar form as in Eq. (38). The hierarchical prior formulation then results in a Laplace prior on $\left(\boldsymbol{x}\mid\lambda\right)$. It is important to note that $\boldsymbol{\gamma}$ is directly responsible for the sparsity of the solution vector $\boldsymbol{x}$ since, for $\gamma_{l}\to 0$, the prior variance of $x^{l}$ tends to zero and the posterior $\mathrm{p}_{x^{l}}\left({x^{l}\mid\boldsymbol{s}}\right)$ will hence be zero, almost surely.

The posterior then follows a multivariate Gaussian distribution which mean and variance admit a closed-form expression, Babacan et al. (2010),

[TABLE]

with $\varLambda:=\text{diag}\left(1/\gamma_{l}\right)$, $1\leq l\leq n_{{D}}$.

Hyperparameters $\left\{\beta,\boldsymbol{\gamma},\lambda,\nu\right\}$ are determined from the data $\boldsymbol{s}$ via a type-II evidence procedure by maximizing the (log)-marginal likelihood, i.e., by integrating out $\boldsymbol{x}$:

[TABLE]

A careful analysis of this optimization problem allows an analytical treatment and a sequential solution procedure where the sparse quantities $\boldsymbol{\mu}$ and $\varSigma$ are iteratively updated, in a memory and computationally efficient way. In particular, manipulations only act on the active set, i.e., predictors $\boldsymbol{p}_{l}$ associated with a non-vanishing entry $\mu^{l}$ through principled sequential additions or deletions of basis element candidates. Our implementation heavily borrows from the constructive Sparse Bayesian Learning (SBL) algorithm, Tipping & Faul (2006) and Babacan et al. (2010).

Once the coefficient random vector $\boldsymbol{x}$ is determined, the full probability distribution $\mathrm{p}\left({\widehat{\boldsymbol{y}}\mid\boldsymbol{s}}\right)$ of the inferred QoI can then be estimated and leads to the probabilistic estimate

[TABLE]

3.6 Algorithm

The dictionary learning procedure is summarized in Algorithm 1.

In the online stage, i.e., once in situ with only some measurements $\boldsymbol{s}$ are available to estimate the QoI $\boldsymbol{y}$, dictionaries ${P}$ and ${D}$ learned in the offline (training) stage are used. The Sparse Bayesian Learning (SBL) procedure described in Sec. 3.5 can be employed to estimate the posterior mean $\boldsymbol{\mu}$ and associated covariance $\varSigma$ of the coefficients $\boldsymbol{x}$ associated with the atoms of the dictionary ${D}$. The estimate $\widehat{\boldsymbol{y}}\left(\boldsymbol{s};{P},{D}\right)$ is finally given as in Eq. (46) in the form of a random field.

Remark 3

If one is not interested in the statistics associated with the estimation, the (deterministic) coefficients $\boldsymbol{x}\left(\boldsymbol{s};{P}\right)$ can be more directly estimated as the solution to the Sparse Coding problem, Eq. (25), typically obtained with an Orthogonal Matching Pursuit (OMP) algorithm. The then deterministic estimate is hence $\widehat{\boldsymbol{y}}\left(\boldsymbol{s};{P},{D}\right)={D}\,\boldsymbol{x}$.

3.7 Numerical complexity

The complexity of the GOBAL method, in its deterministic version, is now discussed.

We assume that the action of a $m\times n$ matrix on a vector is of complexity $2\,m\,n$ and denote the number of iterations of the algorithm by $r$. We let $n_{\mathrm{up}}\leftarrow n_{{D}}$.

First, the training set matrix ${Y}$ is decomposed into $Q^{\left({Y}\right)}$ and $R^{\left({Y}\right)}$ via the eigenvalue decomposition of ${Y}^{\mathsf{T}}\,{Y}$, see Eq. (32). The cost associated with forming the symmetric matrix is ${n_{y}}\,n_{{\mathrm{train}}}^{2}$, while the eigenvalue decomposition takes $\mathcal{O}\left(n_{{D}}\,n_{{\mathrm{train}}}^{2}\right)$ operations. The cost of the initial decomposition of the training set is then

[TABLE]

In the Estimation Codebook Update step, the numerical complexity associated with computing the pseudo-inverse of the full row-rank $n_{{D}}\times n_{{\mathrm{train}}}$-matrix ${X}$ is:

[TABLE]

at leading order, Möller (1999). It is cubic in the small number of atoms and only linear in the number of training elements, hence allowing a fast evaluation. Note that this is the complexity of the plain vanilla pseudo-inverse and that the used QRginv algorithm has a lower complexity and enjoys superior computational performance.

Relying on a Cholesky-based variant of the OMP algorithm for the Sparse Coding step, Cotter et al. (1999); Blumensath & Davies (2008), the number of operations is found to be, Rubinstein et al. (2008),

[TABLE]

Matrix operations involved in the Features Codebook Update step are of similar nature to those of the approximate K-SVD discussed in Rubinstein et al. (2008), so that one can borrow from their analysis. The number of operations per signal is here:

[TABLE]

The framework of this study is limited observations compared to the size of training set. Further, the size of the dictionary is a small multiple $m\in\mathbb{N}^{\ast}$ of the number of sensors. One then has $n_{s}\leq n_{{D}}\leq n_{{\mathrm{train}}}\leq{n_{y}}$ and the total cost $\mathscr{J}$ approximates as

[TABLE]

Initialization of the training set aside (line 2 of the algorithm), the complexity of our core dictionary learning procedure is then seen to mainly scale with the square of the number of sensors and the square of the size of the training set. The bottleneck of the numerical complexity is the Sparse Coding step, which scales with $n_{{\mathrm{train}}}^{2}$. Note however that this step is embarrassingly parallelizable among the $n_{{\mathrm{train}}}$ different elements of the training set, so that the walltime cost of this step, carried-out over $n_{{\mathrm{train}}}$ cores, only scales with $n_{{\mathrm{train}}}$.

For a very large training set, alternatives to the present approach can be derived, relying on an online, or mini-batch, treatment of the training set, e.g., inspired from the work of Mairal et al. (2010).

4 Numerical experiments

4.1 Configuration

The methodology introduced above is now demonstrated on a numerical simulation of the two-dimensional incompressible flow over an open cavity using our in-house code olorin. The cavity aspect ratio is 2 (deep cavity). The numerical grid comprises $296\times 128$ nodes on a structured, non-uniform, staggered mesh. The flow is from left to right on the plots and a laminar boundary layer flow profile is prescribed at the entrance of the numerical domain. The flow Reynolds number $Re$, based on the cavity length and the uniform flow rate velocity, is $Re=7500$. The Navier-Stokes equations are solved using a prediction-projection method and a finite-volume second-order centered scheme in space. Viscous terms are evaluated implicitly while convective fluxes are estimated with a linear Adams-Bashford scheme. Details about the numerical settings may be found in Rizi et al. (2014); Podvin et al. (2006).

In the aim of being as realistic as possible, the retained sensors mimic actual devices routinely used in the flow control community. In particular, the sensors should be physically mounted on a solid wall and hence cannot be located in the bulk flow. They are here modeled as noisy point devices measuring the local wall shear stress, i.e.:

[TABLE]

with $\boldsymbol{z}_{s}^{(j)}$ the location of the $j$-th sensor $\boldsymbol{u}\left(\boldsymbol{z}\right)$ the fluid velocity vector and $n$ the distance to the wall along the local normal direction.

Sensor placement is of critical importance for the efficiency of any recovery procedure and is a well known problem. Many methods have been introduced to address this issue, as briefly discussed in Sec. 1. Here, no sensor placement method is employed. In an attempt to convincingly demonstrate the performance of our GOBAL method against alternatives approaches such as PCA and K-SVD, the sensors all these methods rely on are selected such that the recovery performance of K-SVD is maximized. Sets of $n_{s}=10$ sensors locations $\left\{\boldsymbol{z}_{s}^{(j)}\right\}_{j=1}^{n_{s}}$ are sampled uniformly at random from the geometry of the cavity wall. For each realization of the sensor set, the recovery performance of the K-SVD-based approach is evaluated. The sensor location finally retained is the one leading to the best recovery performance. The performance of the K-SVD-based approach is thus significantly favored by this choice.

The objective of the estimation procedure is to recover the field of velocity magnitude. The training set comprises $n_{{\mathrm{train}}}=1000$ instances, $\left\{\boldsymbol{s}_{i},\boldsymbol{y}_{i}\right\}_{i=1}^{n_{{\mathrm{train}}}}$, and the recovery performance is quantified in terms of relative error norm:

[TABLE]

with $\left\{\boldsymbol{s}_{i}^{\mathrm{CV}},\boldsymbol{y}_{i}^{\mathrm{CV}}\right\}_{i=1}^{n_{\mathrm{CV}}}$ a set of $n_{\mathrm{CV}}=246$ instances not contained in the training set.

Unless explicitly specified, the GOBAL-based estimate is the posterior mean.

4.2 Dictionaries

In addition to our method GOBAL, the performance of a recovery using both a PCA- and a K-SVD-based dictionary is studied. The PCA dictionary is given by the set of $n_{{D}}$ dominant eigenvectors of the empirical correlation matrix ${Y}^{\mathsf{T}}{Y}$, see Eq. (9), while the K-SVD dictionary is determined as the solution of Eq. (15), with $K=n_{s}$ since this is the maximum number of atoms which can be informed by the available measurements. With either of these dictionaries, determined solely from ${Y}$, and not from $\left\{{S},{Y}\right\}$, the estimation is given by $\widehat{\boldsymbol{y}}_{i}={D}\,\widehat{\boldsymbol{x}}_{i}$, $1\leq i\leq n_{\mathrm{CV}}$, with $\widehat{\boldsymbol{x}}_{i}$ minimizing the data misfit in the sense of Eq. (11) (PCA approach) or via a Sparse Coding step, Eq. (21) (K-SVD approach).

Since K-SVD and GOBAL are iterative procedures, the PCA-based dictionary is used as initial condition in both these methods. Note that the mean field $\overline{\boldsymbol{y}}$ has been subtracted from the training set so that dictionaries are learned from centered data.

In our cavity flow simulation, the flow is horizontal, from left to right. Spatial oscillations develop in the shear layer above the cavity, while the left-most part of the domain essentially remains still and the flow varies very little about its mean $\overline{\boldsymbol{y}}$. A large-scale structure within the cavity as produced by the flow recirculation. One can refer to Rizi et al. (2014) for a detailed discussion of the physics of this flow configuration.

The first and 10-th modes of the dictionaries from the different learning strategies are plotted in Fig. 1 for illustration. The main structures of the modes tend to reproduce the salient features of the flow. While the PCA and K-SVD dictionary look alike, at least for these two modes, the GOBAL-based dictionary significantly departs from these two, while still exhibiting similar trends.

4.3 Dependence on $n_{{D}}$

A distinct feature of both the K-SVD and the present GOBAL approaches is that they rely on redundant approximations and allow to consistently inform a dictionary with more atoms than the number of sensors. At most $n_{s}$ atoms are inferred from the data while the other atoms are associated with a vanishing coefficient. In contrast, the PCA-based approach informs the modes via Eq. (11) and, among the infinite number of solutions when $n_{{D}}>n_{s}$, selects the one with the least $L^{2}$-norm.

The recovery error $\varepsilon$ achieved by these different methods is plotted in Fig. 2 as a function of the size $n_{{D}}$ of the dictionary for $n_{s}=10$ sensors. When $n_{{D}}=10=n_{s}$, all dictionary modes can be informed from the available data. Yet, both the PCA- and the K-SVD-based approaches achieve a poor performance, with $\varepsilon\approx 0.8$. The conditioning of the observation of the dictionary atoms via $F\,{D}$ is poor and prevents to accurately estimate ${X}^{\mathrm{CV}}$ from ${S}^{\mathrm{CV}}$. In this particular configuration, the sparsity constraint $\left\|{\boldsymbol{x}_{i}}\right\|_{0}\leq n_{s}$ is always inactive and the sparse coding step, Eq. (21), essentially reduces to the minimization of the data misfit (11). However, this minimization is achieved via the greedy OMP algorithm and potentially deteriorates the recovery performance compared to (11). It results that, for $n_{{D}}=n_{s}$, the K-SVD approach performs worst. In contrast, the present GOBAL approach enjoys a well conditioned predictor operator ${P}$ and accurately estimates ${X}^{\mathrm{CV}}$ and ${Y}^{\mathrm{CV}}$, with $\varepsilon\approx 0.2$.

When the size of the dictionary grows beyond $n_{{D}}=n_{s}$, the performance of the PCA-based approach essentially gradually deteriorates as estimating $\boldsymbol{x}_{i}$ from $\boldsymbol{s}_{i}^{\mathrm{CV}}$ via (11) becomes increasingly ill-posed. In contrast, the performance of the K-SVD-based approach remains roughly constant, $\varepsilon\approx 0.7$. However, since ${D}$ stems from ${Y}$ solely, observations $\boldsymbol{s}_{i}^{\mathrm{CV}}$ do not necessarily admit a sparse representation in the linear span of the columns of $F\,{D}$ and the estimation of the coefficients via the sparse coding step (21) is likely to lead to a poor set of coefficients, hence a poor $\varepsilon$. The GOBAL approach, however, precisely addresses this issue and explicitly determines the predictor operator such that ${S}^{\mathrm{CV}}$ is likely to admit a compressible representation in $\operatorname*{{\mathrm{span}}}\left\{\boldsymbol{p}_{l}\right\}_{l}$. The performance is seen to be good and to (slightly) improve when the size of the dictionary increases, allowing for better tailored atoms to be used for reconstructing a given field $\boldsymbol{y}^{\mathrm{CV}}$ from observations $\boldsymbol{s}^{\mathrm{CV}}$.

The recovery performance is also illustrated in Fig. 3 where the different methods are employed to infer a field chosen at random from $\left\{\boldsymbol{y}_{i}^{\mathrm{CV}}\right\}_{i}$ from associated observations, as given by the shear stress measurement at the locations of the sensors. The sensors are denoted with red dots in the figure. The recovered field is plotted for the PCA- (top right), K-SVD- (bottom right) and GOBAL-based (bottom left) methods, along with the actual field (top left). The superior performance of the GOBAL approach is clearly visible.

4.4 Dependence on $n_{s}$

The quality and quantity of information is pivotal to a reliable estimation. The impact of the quantity is here studied in terms of number of sensors. From the 10 sensors considered in the above section, the different recovery methods are now applied only relying on the first $1\leq n_{s}\leq 10$ sensors from the full set of 10 sensors. Note that the set of 10 sensors considered so far was determined, from a large set of sensor placement realizations, as that giving the best recovery performance for the K-SVD-based approach. Within this “best set”, they are not sorted or ranked in any way.

The influence of the number of sensors on the recovery accuracy is illustrated in Fig. 4 for different sizes of dictionary. It is important to note that the size of the dictionary here increases with the number of sensors, keeping the ratio $m=n_{{D}}/n_{s}$ constant. For a given dictionary, the performance of the GOBAL approach is seen to significantly improve when more information becomes available (increasing $n_{s}$). In contrast, the performance of the K-SVD-based approach hardly improves and remains rather poor, while the PCA-based performance does not improve and remains around a relative error $\varepsilon\approx 1$. Again, the GOBAL method leads to a situation where the observations are likely to admit a sparse representation in the linear span of the predictor operator, hence allowing an accurate recovery, while no such situation occurs with the other two approaches. It is remarkable that this superior performance of GOBAL is achieved despite the fact that the sensors are located to be “optimal” for another method, namely K-SVD.

4.5 Robustness

An important point in most, if not all, numerical methods is robustness. To be applicable in practice, solution methods should gracefully deteriorate the solution accuracy when input data are affected by a weak deviation from their nominal values. These deviations may come from round-off and truncation errors in the numerics. However, in the present context, the largest (by far) source of error, also termed noise, is expected to originate from the sensor measurements. Characterizing the influence of measurement noise onto the solution accuracy is necessary to demonstrate the benefit and performance of our method.

The noise is assumed to be a Gaussian vector with zero-mean and diagonal variance matrix, as discussed in Sec. 3.5.1, $\boldsymbol{\eta}\sim\mathcal{N}\left(\boldsymbol{0},\delta^{2}\,I_{n_{s}}\right)$. Our GOBAL method determines the posterior density of the estimate under an additive noise assumption, see Eq. (37). However, we here test our method in a context different from the one it is derived for: the accuracy of the solution is studied as the variance $\delta^{2}\,I_{n_{s}}$ of a multiplicative measurement noise varies, that is, the measurement data are actually given by

[TABLE]

with $\odot$ the Hadamard product. Noise affects both the training data and the in situ measurements and a given QoI is then associated with several measurements. The GOBAL algorithm is then run on an extended training set $\left\{{Y},\left\{{S}\left(\theta_{\eta}^{(j)}\right)\right\}_{j}\right\}$. In this work, five realizations, $j=1,\ldots,5$, of each measurement vector are considered.

The Figure 5 shows the evolution of the recovery accuracy $\varepsilon$ as a function of the measurement noise standard deviation $\delta$ in the training set ${S}$, and the noise standard deviation $\delta^{\mathrm{CV}}$ in the actual measurements ${S}^{\mathrm{CV}}$. Several comments are in order. The accuracy of the recovery decreases when the noise variance $\delta^{\mathrm{CV}}$ increases. However, when the noise is strong ($\delta^{\mathrm{CV}}=0.5$), the best performance is achieved with dictionaries trained with strong noise as well, $\delta=0.3$ and $\delta=0.5$. If properly trained, the dictionaries learned from noisy data are hence more robust to noise in the measurements. However, this is at the price of a reduced performance if the measurement noise is low, in which case dictionaries trained without noise perform best. This shows how the learning can be made robust to the expected level of noise of in situ data and that the best performance is obtained when noise in the training step is as similar to that actually affecting the sensors as possible.

Another observation is that, for any noise level, the SBL-based version of GOBAL always performs significantly better than its OMP-based counterpart. There are two reasons for this. First, the OMP-variant of GOBAL is here used with $K=n_{s}$, not allowing any redundancy in the data to inform the non-vanishing mode coefficients, and hence being prone to inaccuracy. Second, OMP is a greedy approach which converges to a local minimum. To address this issue, several initial conditions are usually used, and the best solution is retained. However, in practical situations as mimicked here, the “best” solution cannot be known, nor a posteriori verified. All solutions from the OMP step are hence equivalent, with no mean of discriminating the “good” from the “poor” solutions. We hence use a unique initial condition for OMP, identical for the training and the online step, with the risk for the OMP-based SC step to reach a poor local minimum. The Sparse Bayesian Learning variant of our SC step does not suffer these two limitations and potentially achieves a better overall recovery performance.

4.6 Probabilistic estimate

The probabilistic estimate provided by our GOBAL method is illustrated in Fig. 6. The posterior mean field, $\widehat{\boldsymbol{y}}={D}\,\boldsymbol{\mu}$ (bottom left plot) is seen to compare well with the truth (top left plot). The diagonal part of the posterior standard deviation, $\left[\text{diag}\left({D}\,\varSigma\,{D}^{\mathsf{T}}\right)\right]^{1/2}$, (bottom right plot) illustrates the spatial regions where noise in the measurements, and more generally the limitation introduced by the likelihood model, leads to uncertainty in the estimate. Clearly, regions with the largest estimated standard deviation are those where the flow exhibits the most fluctuations, i.e., the downstream part of the shear layer and the main recirculation vortex within the cavity. In the top right plot, the (absolute value) reconstruction error of the posterior mean, $|{\boldsymbol{y}-\widehat{\boldsymbol{y}}}|$, also exhibits the same general pattern as the posterior standard deviation. That is, regions of the estimated field identified as the most uncertain are indeed roughly those where the error is the largest. It is however important to note that the error is low, even its maximum $\left\|{\boldsymbol{y}-\widehat{\boldsymbol{y}}}\right\|_{\infty}\approx 1.6\,10^{-2}$, compared to the field itself which lies between 0 and about 1.4. This is remarkable since the estimation method targets a reconstruction in the $L^{2}$-sense rather than in the $L^{\infty}$-sense.

5 Closing remarks

In this paper, an estimation technique for inferring high-dimensional quantities is introduced. It relies on learning an approximation basis, the dictionary, in which the quantity of interest is estimated from a set of observations by minimizing the data misfit. The original contribution of the present work is that the dictionary is learned with the specific constraint that it has to be observable from the available sensors. This makes estimation of the associated coefficients more accurate and reliable. An algorithm is presented and allows to derive an observable dictionary from a training set of observations and fields to be inferred. Estimation in the online phase, when in situ, is achieved via an algorithm combining Bayesian estimation and sparse recovery. It results in a probabilistic estimation of the quantity to be inferred. The formulation of the proposed recovery procedure does not involve the size of the high-dimensional quantity of interest but instead scales with the size of the training set, usually much smaller. It results in computationally efficient learning. The method is illustrated in the case of the two-dimensional flow over an open cavity whose velocity field is to be estimated from a small set of point sensors located on the wall. Comparisons with established PCA- and K-SVD-based estimation methods show the superior performance of the present approach.

This study clearly emphasizes the pivotal role of the observations quality onto the recovery performance. While sensor placement is a well-known aspect, the need for a joint derivation of the estimation procedure and determination of an approximation basis is clearly demonstrated in the present work. We strongly advocate learning both the representation basis and the recovery procedure. While not done here, sensor placement could, and actually should, be learned at the same time, as it strongly impacts the resulting observability properties of the dictionary. An integrated approach, accounting for all these aspects as early in the estimation procedure as possible, is key to an efficient and reliable method.

The present work is a first step towards a more general framework. In addition to coupling dictionary learning with sensor placement, future extensions include nonlinear observers where the measurements are not well approximated in the image space of a linear predictor operator. Derivation of the learning step in this more general framework leads to a bi-level optimization problem and will be presented in a future paper.

The estimation procedure introduced here relies on the minimization of the reconstruction residual $L^{2}$-norm. While a common approach, this choice is questionable in some situations. For instance, while suitable for minimizing the estimation error in amplitude, this is a poor choice for minimizing the error in phase. As an example, consider a spatially-located feature well estimated but suffering a slight spatial shift. The associated $L^{2}$-norm estimation error would then be large while the recovery is “good” in the eyeball norm. This motivates alternative norms to be used. Among other choices, Wasserstein distances are commonly used in some other contexts, e.g., image processing, and may be a more appropriate choice. To alleviate the significant computational overhead arising from working with Wasserstein distances, one could instead rely on a framework stable to deformations as, say, under a scattering transform preserving some invariants, to some extent, Mallat (2012).

More generally, an efficient estimation procedure first learns a suitable transformation such that, in the resulting framework, the inverse problem formulates as a problem easier to solve than its original counterpart. Here, transformations promoting sparse solutions where exploited, allowing efficient sparse recovery tools to be used. Reformulating the estimation problem in a reproducing kernel Hilbert space (RKHS), with a suitable kernel and inner product learned from the data, is a more general approach and may allow to “untangle” a complex relationship between, say observations and approximation coefficients, paving the way for an efficient use of tools from the world of linear methods. Learning the manifold onto which the data lie, and a suitable kernel, allows application of techniques presented in this work to a more general context and provides new perspectives for inverse problems. Restricting to only one example, the current work may readily be extended into a filter for observing dynamical systems. Relying on extended observations, as collected over a given time-horizon, the present observable dictionary learning method can be reformulated in a suitable RKHS learned from the data. Since the estimation problem then relies on an extended set of observations, we expect superior performance as compared to the present static estimator.

These developments are the subject of on-going efforts.

Acknowledgement

The first author (LM) gratefully acknowledges stimulating discussions with Youssef M Marzouk and Alessio Spantini.

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.