This work proposes and analyzes a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random coefficients and …
This work proposes and analyzes a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random coefficients and forcing terms (input data of the model). To compute solution statistics, the sparse grid stochastic collocation method uses approximate solutions, produced here by finite elements, corresponding to a deterministic set of points in the random input space. This naturally requires solving uncoupled deterministic problems as in the Monte Carlo method. If the number of random variables needed to describe the input data is moderately large, full tensor product spaces are computationally expensive to use due to the curse of dimensionality. In this case the sparse grid approach is still expected to be competitive with the classical Monte Carlo method. Therefore, it is of major practical relevance to understand in which situations the sparse grid stochastic collocation method is more efficient than Monte Carlo. This work provides error estimates for the fully discrete solution using $L^q$ norms and analyzes the computational efficiency of the proposed method. In particular, it demonstrates algebraic convergence with respect to the total number of collocation points and quantifies the effect of the dimension of the problem (number of input random variables) in the final estimates. The derived estimates are then used to compare the method with Monte Carlo, indicating for which problems the former is more efficient than the latter. Computational evidence complements the present theory and shows the effectiveness of the sparse grid stochastic collocation method compared to full tensor and Monte Carlo approaches.
This work proposes and analyzes a stochastic collocation method for solving elliptic partial differential equations with random coefficients and forcing terms. These input data are assumed to depend on a …
This work proposes and analyzes a stochastic collocation method for solving elliptic partial differential equations with random coefficients and forcing terms. These input data are assumed to depend on a finite number of random variables. The method consists of a Galerkin approximation in space and a collocation in the zeros of suitable tensor product orthogonal polynomials (Gauss points) in the probability space, and naturally leads to the solution of uncoupled deterministic problems as in the Monte Carlo approach. It treats easily a wide range of situations, such as input data that depend nonlinearly on the random variables, diffusivity coefficients with unbounded second moments, and random variables that are correlated or even unbounded. We provide a rigorous convergence analysis and demonstrate exponential convergence of the "probability error" with respect to the number of Gauss points in each direction of the probability space, under some regularity assumptions on the random input data. Numerical examples show the effectiveness of the method. Finally, we include a section with developments posterior to the original publication of this work. There we review sparse grid stochastic collocation methods, which are effective collocation strategies for problems that depend on a moderately large number of random variables.
In this paper we present some theoretical results on the Arbitrary Lagrangian Eulerian (ALE) formulation. This formulation may be used when dealing with moving domains and consists in recasting the …
In this paper we present some theoretical results on the Arbitrary Lagrangian Eulerian (ALE) formulation. This formulation may be used when dealing with moving domains and consists in recasting the governing differential equation and the related weak formulation in a frame of reference moving with the domain. The ALE technique is first presented in the whole generality for conservative equations and a result on the regularity of the underlying mapping is proven. In a second part of the work, the stability property of two types of finite element ALE schemes for parabolic evolution problems are analyzed and its relation with the so-called Geometric Conservation Laws is addressed.
In this work we focus on the numerical approximation of the solution u of a linear elliptic PDE with stochastic coefficients. The problem is rewritten as a parametric PDE and …
In this work we focus on the numerical approximation of the solution u of a linear elliptic PDE with stochastic coefficients. The problem is rewritten as a parametric PDE and the functional dependence of the solution on the parameters is approximated by multivariate polynomials. We first consider the stochastic Galerkin method, and rely on sharp estimates for the decay of the Fourier coefficients of the spectral expansion of u on an orthogonal polynomial basis to build a sequence of polynomial subspaces that features better convergence properties, in terms of error versus number of degrees of freedom, than standard choices such as Total Degree or Tensor Product subspaces. We consider then the Stochastic Collocation method, and use the previous estimates to introduce a new class of Sparse Grids, based on the idea of selecting a priori the most profitable hierarchical surpluses, that, again, features better convergence properties compared to standard Smolyak or tensor product grids. Numerical results show the effectiveness of the newly introduced polynomial spaces and sparse grids.
In the solution of fluid-structure interaction (FSI) problems, partitioned procedures are modular algorithms that involve separate fluid and structure solvers that interact in an iterative framework through the exchange of …
In the solution of fluid-structure interaction (FSI) problems, partitioned procedures are modular algorithms that involve separate fluid and structure solvers that interact in an iterative framework through the exchange of suitable transmission conditions at the FS interface. In this work we study, using Fourier analysis, the convergence of partitioned algorithms based on Robin transmission conditions. We derive, for different models of the fluid and the structure, a frequency-dependent reduction factor at each iteration of the partitioned algorithm, which is minimized by choosing optimal values of the coefficients in the Robin transmission conditions. Two-dimensional numerical results are also reported, which highlight the effectiveness of the optimization procedure.
We propose an innovative method for the accurate estimation of surfaces and spatial fields when prior knowledge of the phenomenon under study is available. The prior knowledge included in the …
We propose an innovative method for the accurate estimation of surfaces and spatial fields when prior knowledge of the phenomenon under study is available. The prior knowledge included in the model derives from physics, physiology, or mechanics of the problem at hand, and is formalized in terms of a partial differential equation governing the phenomenon behavior, as well as conditions that the phenomenon has to satisfy at the boundary of the problem domain. The proposed models exploit advanced scientific computing techniques and specifically make use of the finite element method. The estimators have a penalized regression form and the usual inferential tools are derived. Both the pointwise and the areal data frameworks are considered. The driving application concerns the estimation of the blood flow velocity field in a section of a carotid artery, using data provided by echo-color Doppler. This applied problem arises within a research project that aims at studying atherosclerosis pathogenesis. Supplementary materials for this article are available online.
In the framework of stochastic processes, the connection between the dynamic programming scheme given by the Hamilton-Jacobi-Bellman equation and a recently proposed control approach based on the Fokker-Planck equation is …
In the framework of stochastic processes, the connection between the dynamic programming scheme given by the Hamilton-Jacobi-Bellman equation and a recently proposed control approach based on the Fokker-Planck equation is discussed. Under appropriate assumptions it is shown that the two strategies are equivalent in the case of expected cost functionals, while the Fokker-Planck formalism allows considering a larger classof objectives. To illustratethe connection between the two control strategies, the cases of an Itō stochastic process and of a piecewise-deterministic process are considered.
We study a class of models at the interface between statistics and numerical analysis. Specifically, we consider nonparametric regression models for the estimation of spatial fields from pointwise and noisy …
We study a class of models at the interface between statistics and numerical analysis. Specifically, we consider nonparametric regression models for the estimation of spatial fields from pointwise and noisy observations, which account for problem-specific prior information, described in terms of a partial differential equation governing the phenomenon under study. The prior information is incorporated in the model via a roughness term using a penalized regression framework. We prove the well-posedness of the estimation problem, and we resort to a mixed equal order finite element method for its discretization. Moreover, we prove the well-posedness and the optimal convergence rate of the proposed discretization method. Finally the smoothing technique is extended to the case of areal data, particularly interesting in many applications.
Abstract We design and analyse the performance of a multilevel ensemble Kalman filter method (MLEnKF) for filtering settings where the underlying state-space model is an infinite-dimensional spatio-temporal process. We consider …
Abstract We design and analyse the performance of a multilevel ensemble Kalman filter method (MLEnKF) for filtering settings where the underlying state-space model is an infinite-dimensional spatio-temporal process. We consider underlying models that needs to be simulated by numerical methods, with discretization in both space and time. The multilevel Monte Carlo sampling strategy, achieving variance reduction through pairwise coupling of ensemble particles on neighboring resolutions, is used in the sample-moment step of MLEnKF to produce an efficent hierarchical filtering method for spatio-temporal models. Under sufficent regularity, MLEnKF is proven to be more efficient for weak approximations than EnKF, asymptotically in the large-ensemble and fine-numerical-resolution limit. Numerical examples support our theoretical findings.
We present a theoretical analysis of the CORSING (COmpRessed SolvING) method for the numerical approximation of partial differential equations based on compressed sensing. In particular, we show that the best …
We present a theoretical analysis of the CORSING (COmpRessed SolvING) method for the numerical approximation of partial differential equations based on compressed sensing. In particular, we show that the best $s$-term approximation of the weak solution of a PDE with respect to a system of $N$ trial functions, can be recovered via a Petrov-Galerkin approach using $m \ll N$ test functions. This recovery is guaranteed if the local $a$-coherence associated with the bilinear form and the selected trial and test bases fulfills suitable decay properties. The fundamental tool of this analysis is the restricted inf-sup property, i.e., a combination of the classical inf-sup condition and the well-known restricted isometry property of compressed sensing.
In this paper, we present a multilevel Monte Carlo (MLMC) version of the Stochastic Gradient (SG) method for optimization under uncertainty, in order to tackle Optimal Control Problems (OCP) where …
In this paper, we present a multilevel Monte Carlo (MLMC) version of the Stochastic Gradient (SG) method for optimization under uncertainty, in order to tackle Optimal Control Problems (OCP) where the constraints are described in the form of PDEs with random parameters. The deterministic control acts as a distributed forcing term in the random PDE and the objective function is an expected quadratic loss. We use a Stochastic Gradient approach to compute the optimal control, where the steepest descent direction of the expected loss, at each iteration, is replaced by independent MLMC estimators with increasing accuracy and computational cost. The refinement strategy is chosen a-priori such that the bias and, possibly, the variance of the MLMC estimator decays as a function of the iteration counter. Detailed convergence and complexity analyses of the proposed strategy are presented and asymptotically optimal decay rates are identified such that the total computational work is minimized. We also present and analyze an alternative version of the multilevel SG algorithm that uses a randomized MLMC estimator at each iteration. Our methodology is validated through a numerical example.
We consider rank-$1$ lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the …
We consider rank-$1$ lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the non-periodic cosine space and Chebyshev space, via tent transform and then cosine transform, to transfer known results from the periodic setting into new insights for the non-periodic settings. Fast discrete cosine transform can be applied for the reconstruction phase. To reduce the size of the auxiliary index set in the associated component-by-component (CBC) construction for the lattice generating vectors, we work with a bi-orthonormal set of basis functions, leading to three methods for function reconstruction in the non-periodic settings. We provide new theory and efficient algorithmic strategies for the CBC construction. We also interpret our results in the context of general function approximation and discrete least-squares approximation.
Gaussian random fields are widely used as building blocks for modeling stochastic processes. This paper is concerned with the efficient representation of d-point correlations for such fields, which in turn …
Gaussian random fields are widely used as building blocks for modeling stochastic processes. This paper is concerned with the efficient representation of d-point correlations for such fields, which in turn enables the representation of more general stochastic processes that can be expressed as a function of one (or several) Gaussian random fields. Our representation consists of two ingredients. In the first step, we replace the random field by a truncated Karhunen--Loève expansion and analyze the resulting error. The parameters describing the d-point correlation can be arranged in a tensor, but its storage grows exponentially in d. To avoid this, the second step consists of approximating the tensor in a low-rank tensor format, the so-called tensor train decomposition. By exploiting the particular structure of the tensor, an approximation algorithm is derived that does not need to form this tensor explicitly and allows processing correlations of order as high as d=20. The resulting representation is very compact, and its use is illustrated for elliptic partial differential equations with random Gaussian forcing terms.
This work embeds a multilevel Monte Carlo (MLMC) sampling strategy into the Monte Carlo step of the ensemble Kalman filter (EnKF), thereby yielding a multilevel ensemble Kalman filter (MLEnKF) which …
This work embeds a multilevel Monte Carlo (MLMC) sampling strategy into the Monte Carlo step of the ensemble Kalman filter (EnKF), thereby yielding a multilevel ensemble Kalman filter (MLEnKF) which has provably superior asymptotic cost to a given accuracy level. The development of MLEnKF for finite-dimensional state-spaces in the work [20] is here extended to models with infinite-dimensional state- spaces in the form of spatial fields. A concrete example is given to illustrate the results.
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency …
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency only, built at few selected values of the parameters. This, in particular, requires matching the respective poles by solving an optimization problem. If the frequency surrogates are constructed by a suitable rational interpolation strategy, frequency and parameters can both be sampled in an adaptive fashion. This, in general, yields frequency surrogates with different numbers of poles, a situation addressed by our proposed algorithm. Moreover, we explain how our method can be applied even in high-dimensional settings, by employing locally-refined sparse grids in parameter space to weaken the curse of dimensionality. Numerical examples are used to showcase the effectiveness of the method, and to highlight some of its limitations in dealing with unbalanced pole matching, as well as with a large number of parameters.
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are …
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are shown to satisfy a discrete variational formulation. By exploiting this property, we establish stability of our schemes: we show that our explicit and semi-implicit versions are conditionally stable under a parabolic type CFL condition which does not depend on the smallest singular value of the DLR solution; whereas our implicit scheme is unconditionally stable. Moreover, we show that, in certain cases, the semi-implicit scheme can be unconditionally stable if the randomness in the system is sufficiently small. Furthermore, we show that these schemes can be interpreted as projector-splitting integrators and are strongly related to the scheme proposed by Lubich et al. [BIT Num. Math., 54:171-188, 2014; SIAM J. on Num. Anal., 53:917-941, 2015], to which our stability analysis applies as well. The analysis is supported by numerical results showing the sharpness of the obtained stability conditions.
This work describes the convergence analysis of a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random …
This work describes the convergence analysis of a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random coefficients and forcing terms (input data of the model). To compute solution statistics, the sparse grid stochastic collocation method uses approximate solutions, produced here by finite elements, corresponding to a deterministic set of points in the random input space. This naturally requires solving uncoupled deterministic problems and, as such, the derived strong error estimates for the fully discrete solution are used to compare the computational efficiency of the proposed method with the Monte Carlo method. Numerical examples illustrate the theoretical results and are used to compare this approach with several others, including the standard Monte Carlo.
The present work concerns the approximation of the solution map S associated to the parametric Helmholtz boundary value problem, i.e. , the map which associates to each (real) wavenumber belonging …
The present work concerns the approximation of the solution map S associated to the parametric Helmholtz boundary value problem, i.e. , the map which associates to each (real) wavenumber belonging to a given interval of interest the corresponding solution of the Helmholtz equation. We introduce a least squares rational Padé-type approximation technique applicable to any meromorphic Hilbert space-valued univariate map, and we prove the uniform convergence of the Padé approximation error on any compact subset of the interval of interest that excludes any pole. This general result is then applied to the Helmholtz solution map S , which is proven to be meromorphic in ℂ, with a pole of order one in every (single or multiple) eigenvalue of the Laplace operator with the considered boundary conditions. Numerical tests are provided that confirm the theoretical upper bound on the Padé approximation error for the Helmholtz solution map.
In this work, we consider the approximation of Hilbert space-valued meromorphic functions that arise as solution maps of parametric PDEs whose operator is the shift of an operator with normal …
In this work, we consider the approximation of Hilbert space-valued meromorphic functions that arise as solution maps of parametric PDEs whose operator is the shift of an operator with normal and compact resolvent, e.g., the Helmholtz equation. In this restrictive setting, we propose a simplified version of the Least-Squares Padé approximation technique studied in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284] following [J. Approx. Theory 95 (1998), pp. 203–2124]. In particular, the estimation of the poles of the target function reduces to a low-dimensional eigenproblem for a Gramian matrix, allowing for a robust and efficient numerical implementation (hence the "fast" in the name). Moreover, we prove several theoretical results that improve and extend those in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284], including the exponential decay of the error in the approximation of the poles, and the convergence in measure of the approximant to the target function. The latter result extends the classical one for scalar Padé approximation to our functional framework. We provide numerical results that confirm the improved accuracy of the proposed method with respect to the one introduced in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284] for differential operators with normal and compact resolvent.
We consider an optimal control problem (OCP) for a partial differential equation (PDE) with random coefficients. The optimal control function is a deterministic, distributed forcing term that minimizes an expected …
We consider an optimal control problem (OCP) for a partial differential equation (PDE) with random coefficients. The optimal control function is a deterministic, distributed forcing term that minimizes an expected quadratic regularized loss functional. For the numerical approximation of this PDE-constrained OCP, we replace the expectation in the objective functional by a suitable quadrature formula and, eventually, discretize the PDE by a Galerkin method. To practically solve such an approximate OCP, we propose an importance sampling version the SAGA algorithm [A. Defazio, F. Bach, and S. Lacoste-Julien, Advances in Neural Information Processing Systems 27, Curran Associates, Red Hook, NY, 2014, pp. 1646--1654], a type of stochastic gradient algorithm with a fixed-length memory term, which computes at each iteration the gradient of the loss functional in only one quadrature point, randomly chosen from a possibly nonuniform distribution. We provide a full error and complexity analysis of the proposed numerical scheme. In particular we compare the complexity of the generalized SAGA algorithm with importance sampling, with that of the stochastic gradient and the conjugate gradient (CG) algorithms, applied to the same discretized OCP. We show that SAGA converges exponentially in the number of iterations as for a CG algorithm and has a similar asymptotic computational complexity, in terms of computational cost versus accuracy. Moreover, it features good preasymptotic properties, as shown by our numerical experiments, which makes it appealing in a limited budget context.
Abstract This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice—a setting that, as pointed out by Zeng et …
Abstract This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice—a setting that, as pointed out by Zeng et al. (Monte Carlo and Quasi-Monte Carlo Methods 2004, Springer, New York, 2006) and Zeng et al. (Constr. Approx. 30: 529–555, 2009), allows fast evaluation by fast Fourier transform, so avoiding the need for a linear solver. The main contribution of the paper is the application to the approximation problem for uncertainty quantification of elliptic partial differential equations, with the diffusion coefficient given by a random field that is periodic in the stochastic variables, in the model proposed recently by Kaarnioja et al. (SIAM J Numer Anal 58(2): 1068–1091, 2020). The paper gives a full error analysis, and full details of the construction of lattices needed to ensure a good (but inevitably not optimal) rate of convergence and an error bound independent of dimension. Numerical experiments support the theory.
In this work, we propose high-order basis-update & Galerkin (BUG) integrators based on explicit Runge-Kutta methods for large-scale matrix differential equations. These dynamical low-rank integrators are high-order extensions of the …
In this work, we propose high-order basis-update & Galerkin (BUG) integrators based on explicit Runge-Kutta methods for large-scale matrix differential equations. These dynamical low-rank integrators are high-order extensions of the BUG integrator and are constructed by performing one time-step of the first-order BUG integrator at each stage of the Runge-Kutta method. In this way, the resulting Runge-Kutta BUG integrator is robust to the presence of small singular values and does not involve backward time-integration steps. We provide an error bound, which shows that the Runge-Kutta BUG integrator retains the order of convergence of the associated Runge-Kutta method until the error reaches a plateau corresponding to the low-rank truncation and which vanishes as the rank increases. This error bound is finally validated numerically on three different test cases. The results demonstrate the high-order convergence of the Runge-Kutta BUG integrator and its superior accuracy compared to other low-rank integrators proposed in the literature.
The Active Subspace (AS) method is a widely used technique for identifying the most influential directions in high-dimensional input spaces that affect the output of a computational model. The standard …
The Active Subspace (AS) method is a widely used technique for identifying the most influential directions in high-dimensional input spaces that affect the output of a computational model. The standard AS algorithm requires a sufficient number of gradient evaluations (samples) of the input output map to achieve quasi-optimal reconstruction of the active subspace, which can lead to a significant computational cost if the samples include numerical discretization errors which have to be kept sufficiently small. To address this issue, we propose a multilevel version of the Active Subspace method (MLAS) that utilizes samples computed with different accuracies and yields different active subspaces across accuracy levels, which can match the accuracy of single-level AS with reduced computational cost, making it suitable for downstream tasks such as function approximation. In particular, we propose to perform the latter via optimally-weighted least-squares polynomial approximation in the different active subspaces, and we present an adaptive algorithm to choose dynamically the dimensions of the active subspaces and polynomial spaces. We demonstrate the practical viability of the MLAS method with polynomial approximation through numerical experiments based on random partial differential equations (PDEs).
Abstract In this manuscript, we present a collective multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty, and develop a …
Abstract In this manuscript, we present a collective multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty, and develop a novel convergence analysis of collective smoothers and collective two-level methods. The multigrid algorithm is based on a collective smoother that at each iteration sweeps over the nodes of the computational mesh, and solves a reduced saddle-point system whose size is proportional to the number N of samples used to discretized the probability space. We show that this reduced system can be solved with optimal O ( N ) complexity. The multigrid method is tested both as a stationary method and as a preconditioner for GMRES on three problems: a linear-quadratic problem, possibly with a local or a boundary control, for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and $$L^1$$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msup> <mml:mi>L</mml:mi> <mml:mn>1</mml:mn> </mml:msup> </mml:math> -norm penalization on the control, in which the multigrid scheme is used as an inner solver within a semismooth Newton iteration; a risk-averse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits excellent performances and robustness with respect to the parameters of interest.
This manuscript presents a framework for using multilevel quadrature formulae to compute the solution of optimal control problems constrained by random partial differential equations. Our approach consists in solving a …
This manuscript presents a framework for using multilevel quadrature formulae to compute the solution of optimal control problems constrained by random partial differential equations. Our approach consists in solving a sequence of optimal control problems discretized with different levels of accuracy of the physical and probability discretizations. The final approximation of the control is then obtained in a postprocessing step, by suitably combining the adjoint variables computed on the different levels. We present a convergence analysis for an unconstrained linear quadratic problem, and detail our framework for the specific case of a Multilevel Monte Carlo quadrature formula. Numerical experiments confirm the better computational complexity of our MLMC approach compared to a standard Monte Carlo sample average approximation, even beyond the theoretical assumptions.
We propose a novel framework of generalised Petrov-Galerkin Dynamical Low Rank Approximations (DLR) in the context of random PDEs. It builds on the standard Dynamical Low Rank Approximations in their …
We propose a novel framework of generalised Petrov-Galerkin Dynamical Low Rank Approximations (DLR) in the context of random PDEs. It builds on the standard Dynamical Low Rank Approximations in their Dynamically Orthogonal formulation. It allows to seamlessly build-in many standard and well-studied stabilisation techniques that can be framed as either generalised Galerkin methods, or Petrov-Galerkin methods. The framework is subsequently applied to the case of Streamine Upwind/Petrov Galerkin (SUPG) stabilisation of advection-dominated problems with small stochastic perturbations of the transport field. The norm-stability properties of two time discretisations are analysed. Numerical experiments confirm that the stabilising properties of the SUPG method naturally carry over to the DLR framework.
We perform an error analysis of a fully discretised Streamline Upwind Petrov Galerkin Dynamical Low Rank (SUPG-DLR) method for random time-dependent advection-dominated problems. The time integration scheme has a splitting-like …
We perform an error analysis of a fully discretised Streamline Upwind Petrov Galerkin Dynamical Low Rank (SUPG-DLR) method for random time-dependent advection-dominated problems. The time integration scheme has a splitting-like nature, allowing for potentially efficient computations of the factors characterising the discretised random field. The method allows to efficiently compute a low-rank approximation of the true solution, while naturally "inbuilding" the SUPG stabilisation. Standard error rates in the L2 and SUPG-norms are recovered. Numerical experiments validate the predicted rates.
A kernel method for estimating a probability density function from an independent and identically distributed sample drawn from such density is presented. Our estimator is a linear combination of kernel …
A kernel method for estimating a probability density function from an independent and identically distributed sample drawn from such density is presented. Our estimator is a linear combination of kernel functions, the coefficients of which are determined by a linear equation. An error analysis for the mean integrated squared error is established in a general reproducing kernel Hilbert space setting. The theory developed is then applied to estimate probability density functions belonging to weighted Korobov spaces, for which a dimension-independent convergence rate is established. Under a suitable smoothness assumption, our method attains a rate arbitrarily close to the optimal rate. Numerical results support our theory.
.In this work, we present, analyze, and implement a class of multilevel Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis–Hastings proposals for Bayesian inverse problems. In this context, …
.In this work, we present, analyze, and implement a class of multilevel Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis–Hastings proposals for Bayesian inverse problems. In this context, the likelihood function involves solving a complex differential model, which is then approximated on a sequence of increasingly accurate discretizations. The key point of this algorithm is to construct highly coupled Markov chains together with the standard multilevel Monte Carlo argument to obtain a better cost-tolerance complexity than a single-level MCMC algorithm. Our method extends the ideas of Dodwell et al., [SIAM/ASA J. Uncertain. Quantif., 3 (2015), pp. 1075–1108] to a wider range of proposal distributions. We present a thorough convergence analysis of the ML-MCMC method proposed, and show, in particular, that (i) under some mild conditions on the (independent) proposals and the family of posteriors, there exists a unique invariant probability measure for the coupled chains generated by our method, and (ii) that such coupled chains are uniformly ergodic. We also generalize the cost-tolerance theorem of Dodwell et al. to our wider class of ML-MCMC algorithms. Finally, we propose a self-tuning continuation-type ML-MCMC algorithm. The presented method is tested on an array of academic examples, where some of our theoretical results are numerically verified. These numerical experiments evidence how our extended ML-MCMC method is robust when targeting some pathological posteriors, for which some of the previously proposed ML-MCMC algorithms fail.KeywordsBayesian inversionmultilevel Monte CarloMarkov chain Monte Carlouncertainty quantificationMSC codes35R3062F1562P3065M3260J0560J22
We present a multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty. The algorithm is based on a collective smoother …
We present a multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty. The algorithm is based on a collective smoother that at each iteration sweeps over the nodes of the computational mesh, and solves a reduced saddle-point system whose size depends on the number $N$ of samples used to discretized the probability space. We show that this reduced system can be solved with optimal $O(N)$ complexity. We test the multigrid method on three problems: a linear-quadratic problem for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and $L^1$-norm penalization on the control, in which the multigrid scheme is used within a semismooth Newton iteration; a risk-adverse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits very good performances and robustness with respect to all parameters of interest.
In this paper we propose an unbiased Monte Carlo maximum likelihood estimator for discretely observed Wright-Fisher diffusions. Our approach is based on exact simulation techniques that are of special interest …
In this paper we propose an unbiased Monte Carlo maximum likelihood estimator for discretely observed Wright-Fisher diffusions. Our approach is based on exact simulation techniques that are of special interest for diffusion processes defined on a bounded domain, where numerical methods typically fail to remain within the required boundaries. We start by building unbiased maximum likelihood estimators for scalar diffusions and later present an extension to the multidimensional case. Consistency results of our proposed estimator are also presented and the performance of our method is illustrated through a numerical example.
In this work, we consider the problem of estimating the probability distribution, the quantile or the conditional expectation above the quantile, the so called conditional-value-at-risk (CVaR), of output quantities of …
In this work, we consider the problem of estimating the probability distribution, the quantile or the conditional expectation above the quantile, the so called conditional-value-at-risk (CVaR), of output quantities of complex random differential models by the Multilevel Monte Carlo (MLMC) method. We follow an approach that recasts the estimation of the above quantities to the computation of suitable parametric expectations. In this work, we present novel computable error estimators for the estimation of such quantities, which are then used to optimally tune the MLMC hierarchy in a continuation type adaptive algorithm. We demonstrate the efficiency and robustness of our adaptive continuation-MLMC in an array of numerical test cases.
In this paper, we set the mathematical foundations of the Dynamical Low-Rank Approximation (DLRA) method for stochastic differential equations (SDEs). DLRA aims at approximating the solution as a linear combination …
In this paper, we set the mathematical foundations of the Dynamical Low-Rank Approximation (DLRA) method for stochastic differential equations (SDEs). DLRA aims at approximating the solution as a linear combination of a small number of basis vectors with random coefficients (low rank format) with the peculiarity that both the basis vectors and the random coefficients vary in time. While the formulation and properties of DLRA are now well understood for random/parametric equations, the same cannot be said for SDEs and this work aims to fill this gap. We start by rigorously formulating a Dynamically Orthogonal (DO) approximation (an instance of DLRA successfully used in applications) for SDEs, which we then generalize to define a parametrization independent DLRA for SDEs. We show local well-posedness of the DO equations and their equivalence with the DLRA formulation. We also characterize the explosion time of the DO solution by a loss of linear independence of the random coefficients defining the solution expansion and give sufficient conditions for global existence.
Accurate and reliable electricity load forecasts are becoming increasingly important as the share of intermittent resources in the system increases. Distribution System Operators (DSOs) are called to accurately forecast their …
Accurate and reliable electricity load forecasts are becoming increasingly important as the share of intermittent resources in the system increases. Distribution System Operators (DSOs) are called to accurately forecast their production and consumption to place optimal bids in the day-ahead market. Violations of their dispatch-plan requires activation of reserve-power which has a direct cost for the DSO, and also necessitate available reserve-capacity. Forecasts must account for the volatility of weather-parameters that impacts both the production and consumption of electricity. If DSO-loads are small or lower-granularity forecasts are needed, traditional statistical methods may fail to provide reliable performance since they rely on a priori statistical distributions of the variables to forecast. In this paper we introduce a probabilistic load forecast (PLF) method based on empirical copulas. Our model is data-driven, does not need a priori assumption on parametric distribution for variables, nor the dependence structure (copula), but employs a kernel density estimate of the underlying distribution using beta kernels that have bounded support on the unit hypercube. The method naturally supports variables with widely different distributions, such as weather data (including forecasted ones) and historic electricity consumption, and produces a conditional probability distribution for every time step in the forecast, which allows inferring the quantiles of interest. The proposed non-parametric approach is highly flexible and can produce meaningful forecasts even at very low aggregated levels (e.g. neighborhoods). We present results from an open dataset and showcase the strength of the model with respect to Quantile Regression using standard probabilistic evaluation metrics.
Abstract The discretization of robust quadratic optimal control problems under uncertainty using the finite element method and the stochastic collocation method leads to large saddle‐point systems, which are fully coupled …
Abstract The discretization of robust quadratic optimal control problems under uncertainty using the finite element method and the stochastic collocation method leads to large saddle‐point systems, which are fully coupled across the random realizations. Despite its relevance for numerous engineering problems, the solution of such systems is notoriously challenging. In this manuscript, we study efficient preconditioners for all‐at‐once approaches using both an algebraic and an operator preconditioning framework. We show in particular that for values of the regularization parameter not too small, the saddle‐point system can be efficiently solved by preconditioning in parallel all the state and adjoint equations. For small values of the regularization parameter, robustness can be recovered by the additional solution of a small linear system, which however couples all realizations. A mean approximation and a Chebyshev semi‐iterative method are proposed to solve this reduced system. We consider a random elliptic partial differential equation whose diffusion coefficient is modeled as an almost surely continuous and positive random field, though not necessarily uniformly bounded and coercive. We further provide estimates of the dependence of the spectrum of the preconditioned system matrix on the statistical properties of the random field and on the discretization of the probability space. Such estimates involve either the first or second moment of the random variables and , where is the spatial domain. The theoretical results are confirmed by numerical experiments, and implementation details are further addressed.
In this work, we consider the problem of estimating the probability distribution, the quantile or the conditional expectation above the quantile, the so called conditional-value-at-risk, of output quantities of complex …
In this work, we consider the problem of estimating the probability distribution, the quantile or the conditional expectation above the quantile, the so called conditional-value-at-risk, of output quantities of complex random differential models by the MLMC method. We follow the approach of (reference), which recasts the estimation of the above quantities to the computation of suitable parametric expectations. In this work, we present novel computable error estimators for the estimation of such quantities, which are then used to optimally tune the MLMC hierarchy in a continuation type adaptive algorithm. We demonstrate the efficiency and robustness of our adaptive continuation-MLMC in an array of numerical test cases.
In this work, we tackle the problem of minimising the Conditional-Value-at-Risk (CVaR) of output quantities of complex differential models with random input data, using gradient-based approaches in combination with the …
In this work, we tackle the problem of minimising the Conditional-Value-at-Risk (CVaR) of output quantities of complex differential models with random input data, using gradient-based approaches in combination with the Multi-Level Monte Carlo (MLMC) method. In particular, we consider the framework of multi-level Monte Carlo for parametric expectations and propose modifications of the MLMC estimator, error estimation procedure, and adaptive MLMC parameter selection to ensure the estimation of the CVaR and sensitivities for a given design with a prescribed accuracy. We then propose combining the MLMC framework with an alternating inexact minimisation-gradient descent algorithm, for which we prove exponential convergence in the optimisation iterations under the assumptions of strong convexity and Lipschitz continuity of the gradient of the objective function. We demonstrate the performance of our approach on two numerical examples of practical relevance, which evidence the same optimal asymptotic cost-tolerance behaviour as standard MLMC methods for fixed design computations of output expectations.
We present a combination technique based on mixed differences of both spatial approximations and quadrature formulae for the stochastic variables to solve efficiently a class of Optimal Control Problems (OCPs) …
We present a combination technique based on mixed differences of both spatial approximations and quadrature formulae for the stochastic variables to solve efficiently a class of Optimal Control Problems (OCPs) constrained by random partial differential equations. The method requires to solve the OCP for several low-fidelity spatial grids and quadrature formulae for the objective functional. All the computed solutions are then linearly combined to get a final approximation which, under suitable regularity assumptions, preserves the same accuracy of fine tensor product approximations, while drastically reducing the computational cost. The combination technique involves only tensor product quadrature formulae, thus the discretized OCPs preserve the (possible) convexity of the continuous OCP. Hence, the combination technique avoids the inconveniences of Multilevel Monte Carlo and/or sparse grids approaches, but remains suitable for high dimensional problems. The manuscript presents an a-priori procedure to choose the most important mixed differences and an asymptotic complexity analysis, which states that the asymptotic complexity is exclusively determined by the spatial solver. Numerical experiments validate the results.
Abstract This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice—a setting that, as pointed out by Zeng et …
Abstract This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice—a setting that, as pointed out by Zeng et al. (Monte Carlo and Quasi-Monte Carlo Methods 2004, Springer, New York, 2006) and Zeng et al. (Constr. Approx. 30: 529–555, 2009), allows fast evaluation by fast Fourier transform, so avoiding the need for a linear solver. The main contribution of the paper is the application to the approximation problem for uncertainty quantification of elliptic partial differential equations, with the diffusion coefficient given by a random field that is periodic in the stochastic variables, in the model proposed recently by Kaarnioja et al. (SIAM J Numer Anal 58(2): 1068–1091, 2020). The paper gives a full error analysis, and full details of the construction of lattices needed to ensure a good (but inevitably not optimal) rate of convergence and an error bound independent of dimension. Numerical experiments support the theory.
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are …
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are shown to satisfy a discrete variational formulation. By exploiting this property, we establish stability of our schemes: we show that our explicit and semi-implicit versions are conditionally stable under a parabolic type CFL condition which does not depend on the smallest singular value of the DLR solution; whereas our implicit scheme is unconditionally stable. Moreover, we show that, in certain cases, the semi-implicit scheme can be unconditionally stable if the randomness in the system is sufficiently small. Furthermore, we show that these schemes can be interpreted as projector-splitting integrators and are strongly related to the scheme proposed by Lubich et al. [BIT Num. Math., 54:171-188, 2014; SIAM J. on Num. Anal., 53:917-941, 2015], to which our stability analysis applies as well. The analysis is supported by numerical results showing the sharpness of the obtained stability conditions.
In the current work we present two generalizations of the Parallel Tempering algorithm, inspired by the so-called continuous-time Infinite Swapping algorithm. Such a method, found its origins in the molecular …
In the current work we present two generalizations of the Parallel Tempering algorithm, inspired by the so-called continuous-time Infinite Swapping algorithm. Such a method, found its origins in the molecular dynamics community, and can be understood as the limit case of the continuous-time Parallel Tempering algorithm, where the (random) time between swaps of states between two parallel chains goes to zero. Thus, swapping states between chains occurs continuously. In the current work, we extend this idea to the context of time-discrete Markov chains and present two Markov chain Monte Carlo algorithms that follow the same paradigm as the continuous-time infinite swapping procedure. We analyze the convergence properties of such discrete-time algorithms in terms of their spectral gap, and implement them to sample from different target distributions. Numerical results show that the proposed methods significantly improve over more traditional sampling algorithms such as Random Walk Metropolis and (traditional) Parallel Tempering.
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency …
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency only, built at few selected values of the parameters. This, in particular, requires matching the respective poles by solving an optimization problem. If the frequency surrogates are constructed by a suitable rational interpolation strategy, frequency and parameters can both be sampled in an adaptive fashion. This, in general, yields frequency surrogates with different numbers of poles, a situation addressed by our proposed algorithm. Moreover, we explain how our method can be applied even in high-dimensional settings, by employing locally-refined sparse grids in parameter space to weaken the curse of dimensionality. Numerical examples are used to showcase the effectiveness of the method, and to highlight some of its limitations in dealing with unbalanced pole matching, as well as with a large number of parameters.
Because of their robustness, efficiency and non-intrusiveness, Monte Carlo methods are probably the most popular approach in uncertainty quantification to computing expected values of quantities of interest (QoIs). Multilevel Monte …
Because of their robustness, efficiency and non-intrusiveness, Monte Carlo methods are probably the most popular approach in uncertainty quantification to computing expected values of quantities of interest (QoIs). Multilevel Monte Carlo (MLMC) methods significantly reduce the computational cost by distributing the sampling across a hierarchy of discretizations and allocating most samples to the coarser grids. For time dependent problems, spatial coarsening typically entails an increased time-step. Geometric constraints, however, may impede uniform coarsening thereby forcing some elements to remain small across all levels. If explicit time-stepping is used, the time-step will then be dictated by the smallest element on each level for numerical stability. Hence, the increasingly stringent CFL condition on the time-step on coarser levels significantly reduces the advantages of the multilevel approach. By adapting the time-step to the locally refined elements on each level, local time-stepping (LTS) methods permit to restore the efficiency of MLMC methods even in the presence of complex geometry without sacrificing the explicitness and inherent parallelism.
We consider rank-$1$ lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the …
We consider rank-$1$ lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the non-periodic cosine space and Chebyshev space, via tent transform and then cosine transform, to transfer known results from the periodic setting into new insights for the non-periodic settings. Fast discrete cosine transform can be applied for the reconstruction phase. To reduce the size of the auxiliary index set in the associated component-by-component (CBC) construction for the lattice generating vectors, we work with a bi-orthonormal set of basis functions, leading to three methods for function reconstruction in the non-periodic settings. We provide new theory and efficient algorithmic strategies for the CBC construction. We also interpret our results in the context of general function approximation and discrete least-squares approximation.
We consider an optimal control problem (OCP) for a partial differential equation (PDE) with random coefficients. The optimal control function is a deterministic, distributed forcing term that minimizes an expected …
We consider an optimal control problem (OCP) for a partial differential equation (PDE) with random coefficients. The optimal control function is a deterministic, distributed forcing term that minimizes an expected quadratic regularized loss functional. For the numerical approximation of this PDE-constrained OCP, we replace the expectation in the objective functional by a suitable quadrature formula and, eventually, discretize the PDE by a Galerkin method. To practically solve such an approximate OCP, we propose an importance sampling version the SAGA algorithm [A. Defazio, F. Bach, and S. Lacoste-Julien, Advances in Neural Information Processing Systems 27, Curran Associates, Red Hook, NY, 2014, pp. 1646--1654], a type of stochastic gradient algorithm with a fixed-length memory term, which computes at each iteration the gradient of the loss functional in only one quadrature point, randomly chosen from a possibly nonuniform distribution. We provide a full error and complexity analysis of the proposed numerical scheme. In particular we compare the complexity of the generalized SAGA algorithm with importance sampling, with that of the stochastic gradient and the conjugate gradient (CG) algorithms, applied to the same discretized OCP. We show that SAGA converges exponentially in the number of iterations as for a CG algorithm and has a similar asymptotic computational complexity, in terms of computational cost versus accuracy. Moreover, it features good preasymptotic properties, as shown by our numerical experiments, which makes it appealing in a limited budget context.
A kernel method for estimating a probability density function (pdf) from an i.i.d. sample drawn from such density is presented. Our estimator is a linear combination of kernel functions, the …
A kernel method for estimating a probability density function (pdf) from an i.i.d. sample drawn from such density is presented. Our estimator is a linear combination of kernel functions, the coefficients of which are determined by a linear equation. An error analysis for the mean integrated squared error is established in a general reproducing kernel Hilbert space setting. The theory developed is then applied to estimate pdfs belonging to weighted Korobov spaces, for which a dimension independent convergence rate is established. Under a suitable smoothness assumption, our method attains a rate arbitrarily close to the optimal rate. Numerical results support our theory.
Because of their robustness, efficiency and non-intrusiveness, Monte Carlo methods are probably the most popular approach in uncertainty quantification to computing expected values of quantities of interest (QoIs). Multilevel Monte …
Because of their robustness, efficiency and non-intrusiveness, Monte Carlo methods are probably the most popular approach in uncertainty quantification to computing expected values of quantities of interest (QoIs). Multilevel Monte Carlo (MLMC) methods significantly reduce the computational cost by distributing the sampling across a hierarchy of discretizations and allocating most samples to the coarser grids. For time dependent problems, spatial coarsening typically entails an increased time-step. Geometric constraints, however, may impede uniform coarsening thereby forcing some elements to remain small across all levels. If explicit time-stepping is used, the time-step will then be dictated by the smallest element on each level for numerical stability. Hence, the increasingly stringent CFL condition on the time-step on coarser levels significantly reduces the advantages of the multilevel approach. To overcome that bottleneck we propose to combine the multilevel approach of MLMC with local time-stepping (LTS). By adapting the time-step to the locally refined elements on each level, the efficiency of MLMC methods is restored even in the presence of complex geometry without sacrificing the explicitness and inherent parallelism. In a careful cost comparison, we quantify the reduction in computational cost for local refinement either inside a small fixed region or towards a reentrant corner.
In this work, we present, analyze, and implement a class of Multi-Level Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis-Hastings proposals for Bayesian inverse problems. In this context, …
In this work, we present, analyze, and implement a class of Multi-Level Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis-Hastings proposals for Bayesian inverse problems. In this context, the likelihood function involves solving a complex differential model, which is then approximated on a sequence of increasingly accurate discretizations. The key point of this algorithm is to construct highly coupled Markov chains together with the standard Multi-level Monte Carlo argument to obtain a better cost-tolerance complexity than a single-level MCMC algorithm. Our method extends the ideas of Dodwell, et al. "A hierarchical multilevel Markov chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow," \textit{SIAM/ASA Journal on Uncertainty Quantification 3.1 (2015): 1075-1108,} to a wider range of proposal distributions. We present a thorough convergence analysis of the ML-MCMC method proposed, and show, in particular, that (i) under some mild conditions on the (independent) proposals and the family of posteriors, there exists a unique invariant probability measure for the coupled chains generated by our method, and (ii) that such coupled chains are uniformly ergodic. We also generalize the cost-tolerance theorem of Dodwell et al., to our wider class of ML-MCMC algorithms. Finally, we propose a self-tuning continuation-type ML-MCMC algorithm (C-ML-MCMC). The presented method is tested on an array of academic examples, where some of our theoretical results are numerically verified. These numerical experiments evidence how our extended ML-MCMC method is robust when targeting some \emph{pathological} posteriors, for which some of the previously proposed ML-MCMC algorithms fail.
The discretization of robust quadratic optimal control problems under uncertainty using the finite element method and the stochastic collocation method leads to large saddle-point systems, which are fully coupled across …
The discretization of robust quadratic optimal control problems under uncertainty using the finite element method and the stochastic collocation method leads to large saddle-point systems, which are fully coupled across the random realizations. Despite its relevance for numerous engineering problems, the solution of such systems is notoriusly challenging. In this manuscript, we study efficient preconditioners for all-at-once approaches using both an algebraic and an operator preconditioning framework. We show in particular that for values of the regularization parameter not too small, the saddle-point system can be efficiently solved by preconditioning in parallel all the state and adjoint equations. For small values of the regularization parameter, robustness can be recovered by the additional solution of a small linear system, which however couples all realizations. A mean approximation and a Chebyshev semi-iterative method are investigated to solve this reduced system. Our analysis considers a random elliptic partial differential equation whose diffusion coefficient $\kappa(x,\omega)$ is modeled as an almost surely continuous and positive random field, though not necessarily uniformly bounded and coercive. We further provide estimates on the dependence of the preconditioned system on the variance of the random field. Such estimates involve either the first or second moment of the random variables $1/\min_{x\in \overline{D}} \kappa(x,\omega)$ and $\max_{x\in \overline{D}}\kappa(x,\omega)$, where $D$ is the spatial domain. The theoretical results are confirmed by numerical experiments, and implementation details are further addressed.
Abstract We design and analyse the performance of a multilevel ensemble Kalman filter method (MLEnKF) for filtering settings where the underlying state-space model is an infinite-dimensional spatio-temporal process. We consider …
Abstract We design and analyse the performance of a multilevel ensemble Kalman filter method (MLEnKF) for filtering settings where the underlying state-space model is an infinite-dimensional spatio-temporal process. We consider underlying models that needs to be simulated by numerical methods, with discretization in both space and time. The multilevel Monte Carlo sampling strategy, achieving variance reduction through pairwise coupling of ensemble particles on neighboring resolutions, is used in the sample-moment step of MLEnKF to produce an efficent hierarchical filtering method for spatio-temporal models. Under sufficent regularity, MLEnKF is proven to be more efficient for weak approximations than EnKF, asymptotically in the large-ensemble and fine-numerical-resolution limit. Numerical examples support our theoretical findings.
We study the Darcy boundary value problem with lognormal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and …
We study the Darcy boundary value problem with lognormal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and approximating the expected value of the stochastic solution of the PDE by the expected value of its Taylor polynomial. The recursive deterministic equation satisfied by the expected value of the Taylor polynomial (first moment equation) is formally derived. Well-posedness and regularity results for the recursion are proved to hold in Sobolev space-valued Hölder spaces with mixed regularity. The recursive first moment equation is then discretized by means of a sparse approximation technique, and the convergence rates are derived.
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency …
We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency only, built at few selected values of the parameters. This, in particular, requires matching the respective poles by solving an optimization problem. If the frequency surrogates are constructed by a suitable rational interpolation strategy, frequency and parameters can both be sampled in an adaptive fashion. This, in general, yields frequency surrogates with different numbers of poles, a situation addressed by our proposed algorithm. Moreover, we explain how our method can be applied even in high-dimensional settings, by employing locally-refined sparse grids in parameter space to weaken the curse of dimensionality. Numerical examples are used to showcase the effectiveness of the method, and to highlight some of its limitations in dealing with unbalanced pole matching, as well as with a large number of parameters.
This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice---a setting that, as pointed out by Zeng, Leung, Hickernell …
This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice---a setting that, as pointed out by Zeng, Leung, Hickernell (MCQMC2004, 2006) and Zeng, Kritzer, Hickernell (Constr. Approx., 2009), allows fast evaluation by fast Fourier transform, so avoiding the need for a linear solver. The main contribution of the paper is the application to the approximation problem for uncertainty quantification of elliptic partial differential equations, with the diffusion coefficient given by a random field that is periodic in the stochastic variables, in the model proposed recently by Kaarnioja, Kuo, Sloan (SIAM J. Numer. Anal., 2020). The paper gives a full error analysis, and full details of the construction of lattices needed to ensure a good (but inevitably not optimal) rate of convergence and an error bound independent of dimension. Numerical experiments support the theory.
Abstract We present and analyze a novel wavelet–Fourier technique for the numerical treatment of multidimensional advection–diffusion–reaction equations based on the COmpRessed SolvING (CORSING) paradigm. Combining the Petrov–Galerkin technique with the …
Abstract We present and analyze a novel wavelet–Fourier technique for the numerical treatment of multidimensional advection–diffusion–reaction equations based on the COmpRessed SolvING (CORSING) paradigm. Combining the Petrov–Galerkin technique with the compressed sensing approach the proposed method is able to approximate the largest coefficients of the solution with respect to a biorthogonal wavelet basis. Namely, we assemble a compressed discretization based on randomized subsampling of the Fourier test space and we employ sparse recovery techniques to approximate the solution to the partial differential equation (PDE). In this paper we provide the first rigorous recovery error bounds and effective recipes for the implementation of the CORSING technique in the multidimensional setting. Our theoretical analysis relies on new estimates for the local $a$-coherence, which measures interferences between wavelet and Fourier basis functions with respect to the metric induced by the PDE operator. The stability and robustness of the proposed scheme are shown by numerical illustrations in the one-, two- and three-dimensional cases.
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are …
We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are shown to satisfy a discrete variational formulation. By exploiting this property, we establish stability of our schemes: we show that our explicit and semi-implicit versions are conditionally stable under a parabolic type CFL condition which does not depend on the smallest singular value of the DLR solution; whereas our implicit scheme is unconditionally stable. Moreover, we show that, in certain cases, the semi-implicit scheme can be unconditionally stable if the randomness in the system is sufficiently small. Furthermore, we show that these schemes can be interpreted as projector-splitting integrators and are strongly related to the scheme proposed by Lubich et al. [BIT Num. Math., 54:171-188, 2014; SIAM J. on Num. Anal., 53:917-941, 2015], to which our stability analysis applies as well. The analysis is supported by numerical results showing the sharpness of the obtained stability conditions.
We study the Darcy boundary value problem with log-normal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and …
We study the Darcy boundary value problem with log-normal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and approximating the expected value of the stochastic solution of the PDE by the expected value of its Taylor polynomial. The recursive deterministic equation satisfied by the expected value of the Taylor polynomial (first moment equation) is formally derived. Well-posedness and regularity results for the recursion are proved to hold in Sobolev space-valued Holder spaces with mixed regularity. The recursive first moment equation is then discretized by means of a sparse approximation technique, and the convergence rates are derived.
An existence result is presented for the dynamical low rank (DLR) approximation for random semi-linear evolutionary equations. The DLR solution approximates the true solution at each time instant by a …
An existence result is presented for the dynamical low rank (DLR) approximation for random semi-linear evolutionary equations. The DLR solution approximates the true solution at each time instant by a linear combination of products of deterministic and stochastic basis functions, both of which evolve over time. A key to our proof is to find a suitable equivalent formulation of the original problem. The so-called Dual Dynamically Orthogonal formulation turns out to be convenient. Based on this formulation, the DLR approximation is recast to an abstract Cauchy problem in a suitable linear space, for which existence and uniqueness of the solution in the maximal interval are established.
Weighted least squares polynomial approximation uses random samples to determine projections of functions onto spaces of polynomials. It has been shown that, using an optimal distribution of sample locations, the …
Weighted least squares polynomial approximation uses random samples to determine projections of functions onto spaces of polynomials. It has been shown that, using an optimal distribution of sample locations, the number of samples required to achieve quasi-optimal approximation in a given polynomial subspace scales, up to a logarithmic factor, linearly in the dimension of this space. However, in many applications, the computation of samples includes a numerical discretization error. Thus, obtaining polynomial approximations with a single level method can become prohibitively expensive, as it requires a sufficiently large number of samples, each computed with a sufficiently small discretization error. As a solution to this problem, we propose a multilevel method that utilizes samples computed with different accuracies and is able to match the accuracy of single-level approximations with reduced computational cost. We derive complexity bounds under certain assumptions about polynomial approximability and sample work. Furthermore, we propose an adaptive algorithm for situations where such assumptions cannot be verified a priori . Finally, we provide an efficient algorithm for the sampling from optimal distributions and an analysis of computationally favorable alternative distributions. Numerical experiments underscore the practical applicability of our method.
We numerically study the bias and the mean square error of the estimator in Spatial Regression with Partial Differential Equation (SR-PDE) regularization. SR-PDE is a novel smoothing technique for data …
We numerically study the bias and the mean square error of the estimator in Spatial Regression with Partial Differential Equation (SR-PDE) regularization. SR-PDE is a novel smoothing technique for data distributed over two-dimensional domains, which allows to incorporate prior information formalized in termof a partial differential equation. This technique also enables an accurate estimation when the shape of the domain is complex and it strongly influences the phenomenon under study.
We study the Darcy boundary value problem with log-normal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and …
We study the Darcy boundary value problem with log-normal permeability field. We adopt a perturbation approach, expanding the solution in Taylor series around the nominal value of the coefficient, and approximating the expected value of the stochastic solution of the PDE by the expected value of its Taylor polynomial. The recursive deterministic equation satisfied by the expected value of the Taylor polynomial (first moment equation) is formally derived. Well-posedness and regularity results for the recursion are proved to hold in Sobolev space-valued H\"older spaces with mixed regularity. The recursive first moment equation is then discretized by means of a sparse approximation technique, and the convergence rates are derived.
We show that multigrid ideas can be used to reduce the computational complexity of estimating an expected value arising from a stochastic differential equation using Monte Carlo path simulations. In …
We show that multigrid ideas can be used to reduce the computational complexity of estimating an expected value arising from a stochastic differential equation using Monte Carlo path simulations. In the simplest case of a Lipschitz payoff and a Euler discretisation, the computational cost to achieve an accuracy of O(ϵ) is reduced from O(ϵ −3 ) to O(ϵ −2 (log ϵ) 2 ). The analysis is supported by numerical results showing significant computational savings.
This work proposes and analyzes a stochastic collocation method for solving elliptic partial differential equations with random coefficients and forcing terms. These input data are assumed to depend on a …
This work proposes and analyzes a stochastic collocation method for solving elliptic partial differential equations with random coefficients and forcing terms. These input data are assumed to depend on a finite number of random variables. The method consists of a Galerkin approximation in space and a collocation in the zeros of suitable tensor product orthogonal polynomials (Gauss points) in the probability space, and naturally leads to the solution of uncoupled deterministic problems as in the Monte Carlo approach. It treats easily a wide range of situations, such as input data that depend nonlinearly on the random variables, diffusivity coefficients with unbounded second moments, and random variables that are correlated or even unbounded. We provide a rigorous convergence analysis and demonstrate exponential convergence of the "probability error" with respect to the number of Gauss points in each direction of the probability space, under some regularity assumptions on the random input data. Numerical examples show the effectiveness of the method. Finally, we include a section with developments posterior to the original publication of this work. There we review sparse grid stochastic collocation methods, which are effective collocation strategies for problems that depend on a moderately large number of random variables.
This work proposes and analyzes a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random coefficients and …
This work proposes and analyzes a Smolyak-type sparse grid stochastic collocation method for the approximation of statistical quantities related to the solution of partial differential equations with random coefficients and forcing terms (input data of the model). To compute solution statistics, the sparse grid stochastic collocation method uses approximate solutions, produced here by finite elements, corresponding to a deterministic set of points in the random input space. This naturally requires solving uncoupled deterministic problems as in the Monte Carlo method. If the number of random variables needed to describe the input data is moderately large, full tensor product spaces are computationally expensive to use due to the curse of dimensionality. In this case the sparse grid approach is still expected to be competitive with the classical Monte Carlo method. Therefore, it is of major practical relevance to understand in which situations the sparse grid stochastic collocation method is more efficient than Monte Carlo. This work provides error estimates for the fully discrete solution using $L^q$ norms and analyzes the computational efficiency of the proposed method. In particular, it demonstrates algebraic convergence with respect to the total number of collocation points and quantifies the effect of the dimension of the problem (number of input random variables) in the final estimates. The derived estimates are then used to compare the method with Monte Carlo, indicating for which problems the former is more efficient than the latter. Computational evidence complements the present theory and shows the effectiveness of the sparse grid stochastic collocation method compared to full tensor and Monte Carlo approaches.
Parametric partial differential equations are commonly used to model physical systems. They also arise when Wiener chaos expansions are used as an alternative to Monte Carlo when solving stochastic elliptic …
Parametric partial differential equations are commonly used to model physical systems. They also arise when Wiener chaos expansions are used as an alternative to Monte Carlo when solving stochastic elliptic problems. This paper considers a model class of second order, linear, parametric, elliptic PDE's in a bounded domain D with coefficients depending on possibly countably many parameters. It shows that the dependence of the solution on the parameters in the diffusion coefficient is analytically smooth. This analyticity is then exploited to prove that under very weak assumptions on the diffusion coefficients, the entire family of solutions to such equations can be simultaneously approximated by multivariate polynomials (in the parameters) with coefficients taking values in the Hilbert space [Formula: see text] of weak solutions of the elliptic problem with a controlled number of terms N. The convergence rate in terms of N does not depend on the number of parameters in V which may be countable, therefore breaking the curse of dimensionality. The discretization of the coefficients from a family of continuous, piecewise linear finite element functions in D is shown to yield finite dimensional approximations whose convergence rate in terms of the overall number N dof of degrees of freedom is the minimum of the convergence rates afforded by the best N-term sequence approximations in the parameter space and the rate of finite element approximations in D for a single instance of the parametric problem.
We consider a finite element approximation of elliptic partial differential equations with random coefficients. Such equations arise, for example, in uncertainty quantification in subsurface flow modeling. Models for random coefficients …
We consider a finite element approximation of elliptic partial differential equations with random coefficients. Such equations arise, for example, in uncertainty quantification in subsurface flow modeling. Models for random coefficients frequently used in these applications, such as log-normal random fields with exponential covariance, have only very limited spatial regularity and lead to variational problems that lack uniform coercivity and boundedness with respect to the random parameter. In our analysis we overcome these challenges by a careful treatment of the model problem almost surely in the random parameter, which then enables us to prove uniform bounds on the finite element error in standard Bochner spaces. These new bounds can then be used to perform a rigorous analysis of the multilevel Monte Carlo method for these elliptic problems that lack full regularity and uniform coercivity and boundedness. To conclude, we give some numerical results that confirm the new bounds.
Recently there has been a growing interest in designing efficient methods for the solution of ordinary/partial differential equations with random inputs. To this end, stochastic Galerkin methods appear to be …
Recently there has been a growing interest in designing efficient methods for the solution of ordinary/partial differential equations with random inputs. To this end, stochastic Galerkin methods appear to be superior to other nonsampling methods and, in many cases, to several sampling methods. However, when the governing equations take complicated forms, numerical implementations of stochastic Galerkin methods can become nontrivial and care is needed to design robust and efficient solvers for the resulting equations. On the other hand, the traditional sampling methods, e.g., Monte Carlo methods, are straightforward to implement, but they do not offer convergence as fast as stochastic Galerkin methods. In this paper, a high-order stochastic collocation approach is proposed. Similar to stochastic Galerkin methods, the collocation methods take advantage of an assumption of smoothness of the solution in random space to achieve fast convergence. However, the numerical implementation of stochastic collocation is trivial, as it requires only repetitive runs of an existing deterministic solver, similar to Monte Carlo methods. The computational cost of the collocation methods depends on the choice of the collocation points, and we present several feasible constructions. One particular choice, basedon sparse grids, depends weakly on the dimensionality of the random space and is more suitable for highly accurate computations of practical applications with large dimensional random inputs. Numerical examples are presented to demonstrate the accuracy and efficiency of the stochastic collocation methods.
We compare the convergence behavior of Gauss quadrature with that of its younger brother, Clenshaw–Curtis. Seven-line MATLAB codes are presented that implement both methods, and experiments show that the supposed …
We compare the convergence behavior of Gauss quadrature with that of its younger brother, Clenshaw–Curtis. Seven-line MATLAB codes are presented that implement both methods, and experiments show that the supposed factor-of-2 advantage of Gauss quadrature is rarely realized. Theorems are given to explain this effect. First, following O'Hara and Smith in the 1960s, the phenomenon is explained as a consequence of aliasing of coefficients in Chebyshev expansions. Then another explanation is offered based on the interpretation of a quadrature formula as a rational approximation of $\log((z+1)/(z-1))$ in the complex plane. Gauss quadrature corresponds to Padé approximation at $z=\infty$. Clenshaw–Curtis quadrature corresponds to an approximation whose order of accuracy at $z=\infty$ is only half as high, but which is nevertheless equally accurate near $[-1,1]$.
We study multivariate tenser product problems in the worst case and average case settings. They are defined on functions of d variables. For arbitrary d, we provide explicit upper bounds …
We study multivariate tenser product problems in the worst case and average case settings. They are defined on functions of d variables. For arbitrary d, we provide explicit upper bounds on the costs of algorithms which compute an ϵ-approximation to the solution. The cost bounds are of the form (c(d) + 2)β1(β2 + β3(ln 1/ϵ)/(d − 1))β4(d − 1)(1/ϵ)β5. Here c(d) is the cost of one function evaluation (or one linear functional evaluation), and βi′s do not depend on d; they are determined by the properties of the problem for d = 1. For certain tensor product problems, these cost bounds do not exceed c(d)Kϵ−p for some numbers K and p, both independent of d. However, the exponents p which we obtain are too large. We apply these general estimates to certain integration and approximation problems in the worst and average case settings. We also obtain an upper bound, which is independent of d, for the number, n(ϵ, d), of points for which discrepancy (with unequal weights) is at most ϵ, n(ϵ, d) ≤ 7.26ϵ−2.454, ∀d, ϵ ≤ 1.
In this work we focus on the numerical approximation of the solution u of a linear elliptic PDE with stochastic coefficients. The problem is rewritten as a parametric PDE and …
In this work we focus on the numerical approximation of the solution u of a linear elliptic PDE with stochastic coefficients. The problem is rewritten as a parametric PDE and the functional dependence of the solution on the parameters is approximated by multivariate polynomials. We first consider the stochastic Galerkin method, and rely on sharp estimates for the decay of the Fourier coefficients of the spectral expansion of u on an orthogonal polynomial basis to build a sequence of polynomial subspaces that features better convergence properties, in terms of error versus number of degrees of freedom, than standard choices such as Total Degree or Tensor Product subspaces. We consider then the Stochastic Collocation method, and use the previous estimates to introduce a new class of Sparse Grids, based on the idea of selecting a priori the most profitable hierarchical surpluses, that, again, features better convergence properties compared to standard Smolyak or tensor product grids. Numerical results show the effectiveness of the newly introduced polynomial spaces and sparse grids.
Based on the parametric deterministic formulation of Bayesian inverse problems with unknown input parameter from infinite-dimensional, separable Banach spaces proposed in Schwab and Stuart (2012 Inverse Problems 28 045003), we …
Based on the parametric deterministic formulation of Bayesian inverse problems with unknown input parameter from infinite-dimensional, separable Banach spaces proposed in Schwab and Stuart (2012 Inverse Problems 28 045003), we develop a practical computational algorithm whose convergence rates are provably higher than those of Monte Carlo (MC) and Markov chain Monte Carlo methods, in terms of the number of solutions of the forward problem. In the formulation of Schwab and Stuart, the forward problems are parametric, deterministic elliptic partial differential equations, and the inverse problem is to determine the unknown diffusion coefficients from noisy observations comprising linear functionals of the system's response. The sparsity of the generalized polynomial chaos representation of the posterior density being implied by sparsity assumptions on the class of the prior (Schwab and Stuart 2012), we design, analyze and implement a class of adaptive, deterministic sparse tensor Smolyak quadrature schemes for the efficient approximate numerical evaluation of expectations under the posterior, given data. The proposed, deterministic quadrature algorithm is based on a greedy, iterative identification of finite sets of most significant, 'active' chaos polynomials in the posterior density analogous to recently proposed algorithms for adaptive interpolation (Chkifa et al 2012 Report 2012-NN, 2013 Math. Modelling Numer. Anal. 47 253–80). Convergence rates for the quadrature approximation are shown, both theoretically and computationally, to depend only on the sparsity class of the unknown, but are bounded independently of the number of random variables activated by the adaptive algorithm. Numerical results for a model problem of coefficient identification with point measurements in a diffusion problem confirm the theoretical results.
We consider the problem of numerically approximating the solution of an elliptic partial differential equation with random coefficients and homogeneous Dirichlet boundary conditions. We focus on the case of a …
We consider the problem of numerically approximating the solution of an elliptic partial differential equation with random coefficients and homogeneous Dirichlet boundary conditions. We focus on the case of a lognormal coefficient and deal with the lack of uniform coercivity and uniform boundedness with respect to the randomness. This model is frequently used in hydrogeology. We approximate this coefficient by a finite dimensional noise using a truncated Karhunen–Loève expansion. We give estimates of the corresponding error on the solution, both a strong error estimate and a weak error estimate, that is, an estimate of the error commited on the law of the solution. We obtain a weak rate of convergence which is twice the strong one. In addition, we give a complete error estimate for the stochastic collocation method in this case, where neither coercivity nor boundedness is stochastically uniform. To conclude, we apply these results of strong and weak convergence to two classical cases of covariance kernel choices, the case of an exponential covariance kernel on a box and the case of an analytic covariance kernel, yielding explicit weak and strong convergence rates.
This paper is concerned with the construction of optimized sparse grid approximation spaces for elliptic pseudodifferential operators of arbitrary order. Based on the framework of tensor-product biorthogonal wavelet bases and …
This paper is concerned with the construction of optimized sparse grid approximation spaces for elliptic pseudodifferential operators of arbitrary order. Based on the framework of tensor-product biorthogonal wavelet bases and stable subspace splittings, we construct operator-adapted subspaces with a dimension smaller than that of the standard full grid spaces but which have the same approximation order as the standard full grid spaces, provided that certain additional regularity assumptions on the solution are fulfilled. Specifically for operators of positive order, their dimension is<inline-formula content-type="math/mathml"><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" alttext="upper O left-parenthesis 2 Superscript upper J Baseline right-parenthesis"><mml:semantics><mml:mrow><mml:mi>O</mml:mi><mml:mo stretchy="false">(</mml:mo><mml:msup><mml:mn>2</mml:mn><mml:mrow class="MJX-TeXAtom-ORD"><mml:mi>J</mml:mi></mml:mrow></mml:msup><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:annotation encoding="application/x-tex">O(2^{J})</mml:annotation></mml:semantics></mml:math></inline-formula>independent of the dimension<inline-formula content-type="math/mathml"><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" alttext="n"><mml:semantics><mml:mi>n</mml:mi><mml:annotation encoding="application/x-tex">n</mml:annotation></mml:semantics></mml:math></inline-formula>of the problem, compared to<inline-formula content-type="math/mathml"><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" alttext="upper O left-parenthesis 2 Superscript upper J n Baseline right-parenthesis"><mml:semantics><mml:mrow><mml:mi>O</mml:mi><mml:mo stretchy="false">(</mml:mo><mml:msup><mml:mn>2</mml:mn><mml:mrow class="MJX-TeXAtom-ORD"><mml:mi>J</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msup><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:annotation encoding="application/x-tex">O(2^{Jn})</mml:annotation></mml:semantics></mml:math></inline-formula>for the full grid space. Also, for operators of negative order the overall cost is significantly in favor of the new approximation spaces. We give cost estimates for the case of continuous linear information. We show these results in a constructive manner by proposing a Galerkin method together with optimal preconditioning. The theory covers elliptic boundary value problems as well as boundary integral equations.
Monte Carlo methods are a very general and useful approach for the estimation of expectations arising from stochastic simulation. However, they can be computationally expensive, particularly when the cost of …
Monte Carlo methods are a very general and useful approach for the estimation of expectations arising from stochastic simulation. However, they can be computationally expensive, particularly when the cost of generating individual stochastic samples is very high, as in the case of stochastic PDEs. Multilevel Monte Carlo is a recently developed approach which greatly reduces the computational cost by performing most simulations with low accuracy at a correspondingly low cost, with relatively few simulations being performed at high accuracy and a high cost. In this article, we review the ideas behind the multilevel Monte Carlo method, and various recent generalizations and extensions, and discuss a number of applications which illustrate the flexibility and generality of the approach and the challenges in developing more efficient implementations with a faster rate of convergence of the multilevel correction variance.
Stochastic collocation methods for approximating the solution of partial differential equations with random input data (e.g., coefficients and forcing terms) suffer from the curse of dimensionality whereby increases in the …
Stochastic collocation methods for approximating the solution of partial differential equations with random input data (e.g., coefficients and forcing terms) suffer from the curse of dimensionality whereby increases in the stochastic dimension cause an explosion of the computational effort. We propose and analyze a multilevel version of the stochastic collocation method that, as is the case for multilevel Monte Carlo (MLMC) methods, uses hierarchies of spatial approximations to reduce the overall computational complexity. In addition, our proposed approach utilizes, for approximation in stochastic space, a sequence of multidimensional interpolants of increasing fidelity which can then be used for approximating statistics of the solution as well as for building high-order surrogates featuring faster convergence rates. A rigorous convergence and computational cost analysis of the new multilevel stochastic collocation method is provided in the case of elliptic equations, demonstrating its advantages compared to standard single-level stochastic collocation approximations as well as MLMC methods. Numerical results are provided that illustrate the theory and the effectiveness of the new multilevel method.
We analyze the convergence of compressive sensing based sampling techniques for the efficient evaluation of functionals of solutions for a class of high-dimensional, affine-parametric, linear operator equations which depend on …
We analyze the convergence of compressive sensing based sampling techniques for the efficient evaluation of functionals of solutions for a class of high-dimensional, affine-parametric, linear operator equations which depend on possibly infinitely many parameters. The proposed algorithms are based on so-called "non-intrusive" sampling of the high-dimensional parameter space, reminiscent of Monte Carlo sampling. In contrast to Monte Carlo, however, a functional of the parametric solution is then computed via compressive sensing methods from samples of functionals of the solution. A key ingredient in our analysis of independent interest consists of a generalization of recent results on the approximate sparsity of generalized polynomial chaos representations (gpc) of the parametric solution families in terms of the gpc series with respect to tensorized Chebyshev polynomials. In particular, we establish sufficient conditions on the parametric inputs to the parametric operator equation such that the Chebyshev coefficients of gpc expansion are contained in certain weighted $\ell _p$-spaces for $0<p\leq 1$. Based on this we show that reconstructions of the parametric solutions computed from the sampled problems converge, with high probability, at the $L_2$, resp. $L_\infty$, convergence rates afforded by best $s$-term approximations of the parametric solution up to logarithmic factors.
Parametrized families of PDEs arise in various contexts such as inverse problems, control and optimization, risk assessment, and uncertainty quantification. In most of these applications, the number of parameters is …
Parametrized families of PDEs arise in various contexts such as inverse problems, control and optimization, risk assessment, and uncertainty quantification. In most of these applications, the number of parameters is large or perhaps even infinite. Thus, the development of numerical methods for these parametric problems is faced with the possible curse of dimensionality. This article is directed at (i) identifying and understanding which properties of parametric equations allow one to avoid this curse and (ii) developing and analysing effective numerical methods which fully exploit these properties and, in turn, are immune to the growth in dimensionality. Part I of this article studies the smoothness and approximability of the solution map, that is, the map $a\mapsto u(a)$ , where $a$ is the parameter value and $u(a)$ is the corresponding solution to the PDE. It is shown that for many relevant parametric PDEs, the parametric smoothness of this map is typically holomorphic and also highly anisotropic, in that the relevant parameters are of widely varying importance in describing the solution. These two properties are then exploited to establish convergence rates of $n$ -term approximations to the solution map, for which each term is separable in the parametric and physical variables. These results reveal that, at least on a theoretical level, the solution map can be well approximated by discretizations of moderate complexity, thereby showing how the curse of dimensionality is broken. This theoretical analysis is carried out through concepts of approximation theory such as best $n$ -term approximation, sparsity, and $n$ -widths. These notions determine a priori the best possible performance of numerical methods and thus serve as a benchmark for concrete algorithms. Part II of this article turns to the development of numerical algorithms based on the theoretically established sparse separable approximations. The numerical methods studied fall into two general categories. The first uses polynomial expansions in terms of the parameters to approximate the solution map. The second one searches for suitable low-dimensional spaces for simultaneously approximating all members of the parametric family. The numerical implementation of these approaches is carried out through adaptive and greedy algorithms. An a priori analysis of the performance of these algorithms establishes how well they meet the theoretical benchmarks.
This paper addresses optimization problems constrained by partial differential equations with uncertain coefficients. In particular, the robust control problem and the average control problem are considered for a tracking type …
This paper addresses optimization problems constrained by partial differential equations with uncertain coefficients. In particular, the robust control problem and the average control problem are considered for a tracking type cost functional with an additional penalty on the variance of the state. The expressions for the gradient and Hessian corresponding to either problem contain expected value operators. Due to the large number of uncertainties considered in our model, we suggest evaluating these expectations using a multilevel Monte Carlo (MLMC) method. Under mild assumptions, it is shown that this results in the gradient and Hessian corresponding to the MLMC estimator of the original cost functional. Furthermore, we show that the use of certain correlated samples yields a reduction in the total number of samples required. Two optimization methods are investigated: the nonlinear conjugate gradient method and the Newton method. For both, a specific algorithm is provided that dynamically decides which and how many samples should be taken in each iteration. The cost of the optimization up to some specified tolerance $\tau$ is shown to be proportional to the cost of a gradient evaluation with requested root mean square error $\tau$. The algorithms are tested on a model elliptic diffusion problem with lognormal diffusion coefficient. An additional nonlinear term is also considered.
We propose an adaptive sparse grid stochastic collocation approach based upon Leja interpolation sequences for approximation of parameterized functions with high-dimensional parameters. Leja sequences are arbitrarily granular (any number of …
We propose an adaptive sparse grid stochastic collocation approach based upon Leja interpolation sequences for approximation of parameterized functions with high-dimensional parameters. Leja sequences are arbitrarily granular (any number of nodes may be added to a current sequence, producing a new sequence) and thus are a good choice for the univariate composite rule used to construct adaptive sparse grids in high dimensions. When undertaking stochastic collocation one is often interested in constructing weighted approximation where the weights are determined by the probability densities of the random variables. This paper establishes that a certain weighted formulation of one-dimensional Leja sequences produces a sequence of nodes whose empirical distribution converges to the corresponding limiting distribution of the Gauss quadrature nodes associated with the weight function. This property is true even for unbounded domains. We apply the Leja sparse grid approach to several high-dimensional problems and demonstrate that Leja sequences are often superior to more standard sparse grid constructions (e.g., Clenshaw--Curtis), at least for interpolatory metrics.
Stochastic sampling methods are arguably the most direct and least intrusive means of incorporating parametric uncertainty into numerical simulations of partial differential equations with random inputs. However, to achieve an …
Stochastic sampling methods are arguably the most direct and least intrusive means of incorporating parametric uncertainty into numerical simulations of partial differential equations with random inputs. However, to achieve an overall error that is within a desired tolerance, a large number of sample simulations may be required (to control the sampling error), each of which may need to be run at high levels of spatial fidelity (to control the spatial error). Multilevel sampling methods aim to achieve the same accuracy as traditional sampling methods, but at a reduced computational cost, through the use of a hierarchy of spatial discretization models. Multilevel algorithms coordinate the number of samples needed at each discretization level by minimizing the computational cost, subject to a given error tolerance. They can be applied to a variety of sampling schemes, exploit nesting when available, can be implemented in parallel and can be used to inform adaptive spatial refinement strategies. We extend the multilevel sampling algorithm to sparse grid stochastic collocation methods, discuss its numerical implementation and demonstrate its efficiency both theoretically and by means of numerical examples.
For the low‐rank approximation of time‐dependent data matrices and of solutions to matrix differential equations, an increment‐based computational approach is proposed and analyzed. In this method, the derivative is projected …
For the low‐rank approximation of time‐dependent data matrices and of solutions to matrix differential equations, an increment‐based computational approach is proposed and analyzed. In this method, the derivative is projected onto the tangent space of the manifold of rank‐r matrices at the current approximation. With an appropriate decomposition of rank‐r matrices and their tangent matrices, this yields nonlinear differential equations that are well suited for numerical integration. The error analysis compares the result with the pointwise best approximation in the Frobenius norm. It is shown that the approach gives locally quasi‐optimal low‐rank approximations. Numerical experiments illustrate the theoretical results.
We discuss in this paper the numerical approximation of fluid-structure interaction (FSI) problems dealing with strong added-mass effect. We propose new semi-implicit algorithms based on inexact block-$LU$ factorization of the …
We discuss in this paper the numerical approximation of fluid-structure interaction (FSI) problems dealing with strong added-mass effect. We propose new semi-implicit algorithms based on inexact block-$LU$ factorization of the linear system obtained after the space-time discretization and linearization of the FSI problem. As a result, the fluid velocity is computed separately from the coupled pressure-structure velocity system at each iteration, reducing the computational cost. We investigate explicit-implicit decomposition through algebraic splitting techniques originally designed for the FSI problem. This approach leads to two different families of methods which extend to FSI the algebraic pressure correction method and the Yosida method, two schemes that were previously adopted for pure fluid problems. Furthermore, we have considered the inexact factorization of the fluid-structure system as a preconditioner. The numerical properties of these methods have been tested on a model problem representing a blood-vessel system.
In this paper quasi-Monte Carlo (QMC) methods are applied to a class of elliptic partial differential equations (PDEs) with random coefficients, where the random coefficient is parametrized by a countably …
In this paper quasi-Monte Carlo (QMC) methods are applied to a class of elliptic partial differential equations (PDEs) with random coefficients, where the random coefficient is parametrized by a countably infinite number of terms in a Karhunen--Loève expansion. Models of this kind appear frequently in numerical models of physical systems, and in uncertainty quantification. The method uses a QMC method to estimate expected values of linear functionals of the exact or approximate solution of the PDE, with the expected value considered as an infinite dimensional integral in the parameter space corresponding to the randomness induced by the random coefficient. The error analysis, arguably the first rigorous application of the modern theory of QMC in weighted spaces, exploits the regularity with respect to both the physical variables (the variables in the physical domain) and the parametric variables (the parameters corresponding to randomness). In the weighted-space theory of QMC methods, “weights,” describing the varying difficulty of different subsets of the variables, are introduced in order to make sure that the high-dimensional integration problem is tractable. It turns out that the weights arising from the present analysis are of a nonstandard kind, being of neither product nor order dependent form, but instead a hybrid of the two---we refer to these as “product and order dependent weights,” or “POD weights” for short. These POD weights are of a simple enough form to permit a component-by-component construction of a randomly shifted lattice rule that has optimal convergence properties for the given weighted space setting. If the terms in the expansion for the random coefficient have an appropriate decay property, and if we choose POD weights that minimize a certain upper bound on the error, then the solution of the PDE belongs to the joint function space needed for the analysis, and the QMC error (in the sense of a root-mean-square error averaged over shifts) is of order $\mathcal{O}(N^{-1+\delta})$ for arbitrary $\delta>0$, where $N$ denotes the number of sampling points in the parameter space. Moreover, this convergence rate and slower convergence rates under more relaxed conditions are achieved under conditions similar to those found recently by Cohen, De Vore, and Schwab [Found. Comput. Math., 10 (2010), pp. 615--646] for the same model by best $N$-term approximation. We analyze the impact of a finite element (FE) discretization on the overall efficiency of the scheme, in terms of accuracy versus overall cost, with results that are comparable to those of the best $N$-term approximation.
Low-rank approximations to large time-dependent matrices and tensors are the subject of this paper. These matrices and tensors either are given explicitly or are the unknown solutions of matrix and …
Low-rank approximations to large time-dependent matrices and tensors are the subject of this paper. These matrices and tensors either are given explicitly or are the unknown solutions of matrix and tensor differential equations. Based on splitting the orthogonal projection onto the tangent space of the low-rank manifold, novel time integrators for obtaining approximations by low-rank matrices and low-rank tensor trains were recently proposed. By standard theory, the Lie--Trotter and Strang projector-splitting methods are first and second order accurate, respectively, but the usual error bounds break down when the low-rank approximation has small singular values. This happens when the singular values of the solution decay without a distinct gap or when the effective rank of the solution is overestimated. On the other hand, the integrators are exact when given time-dependent matrices or tensors are already of the prescribed rank. We provide an error analysis which unifies these properties. We show that in cases where the exact solution is an $\varepsilon$-perturbation of a low-rank matrix or tensor train, the error of the projector-splitting integrator is favorably bounded in terms of $\varepsilon$ and the stepsize, independently of the smallness of the singular values. Such a result does not hold for any standard integrator. Numerical experiments illustrate the theory.