quasi-Monte Carlo Research Articles

The Representative Points of Generalized Alpha Skew-t Distribution and Applications

Assuming the underlying statistical distribution of data is critical in information theory, as it impacts the accuracy and efficiency of communication and the definition of entropy. The real-world data are widely assumed to follow the normal distribution. To better comprehend the skewness of the data, many models more flexible than the normal distribution have been proposed, such as the generalized alpha skew-t (GAST) distribution. This paper studies some properties of the GAST distribution, including the calculation of the moments, and the relationship between the number of peaks and the GAST parameters with some proofs. For complex probability distributions, representative points (RPs) are useful due to the convenience of manipulation, computation and analysis. The relative entropy of two probability distributions could have been a good criterion for the purpose of generating RPs of a specific distribution but is not popularly used due to computational complexity. Hence, this paper only provides three ways to obtain RPs of the GAST distribution, Monte Carlo (MC), quasi-Monte Carlo (QMC), and mean square error (MSE). The three types of RPs are utilized in estimating moments and densities of the GAST distribution with known and unknown parameters. The MSE representative points perform the best among all case studies. For unknown parameter cases, a revised maximum likelihood estimation (MLE) method of parameter estimation is compared with the plain MLE method. It indicates that the revised MLE method is suitable for the GAST distribution having a unimodal or unobvious bimodal pattern. This paper includes two real-data applications in which the GAST model appears adaptable to various types of data.

Stabilized weighted reduced order methods for parametrized advection-dominated optimal control problems governed by partial differential equations with random inputs

Abstract In this work, we analyze Parametrized Advection-Dominated distributed Optimal Control Problems with random inputs in a Reduced Order Model (ROM) context. All the simulations are initially based on a finite element method (FEM) discretization; moreover, a space-time approach is considered when dealing with unsteady cases. To overcome numerical instabilities that can occur in the optimality system for high values of the Péclet number, we consider a Streamline Upwind Petrov–Galerkin technique applied in an optimize-then-discretize approach. We combine this method with the ROM framework in order to consider two possibilities of stabilization: Offline-Only stabilization and Offline-Online stabilization. Moreover we consider random parameters and we use a weighted Proper Orthogonal Decomposition algorithm in a partitioned approach to deal with the issue of uncertainty quantification. Several quadrature techniques are used to derive weighted ROMs: tensor rules, isotropic sparse grids, Monte-Carlo and quasi Monte-Carlo methods. We compare all the approaches analyzing relative errors between the FEM and ROM solutions and the computational efficiency based on the speedup-index.

Dimension reduction for Quasi-Monte Carlo methods via quadratic regression

Quasi-Monte Carlo (QMC) methods have been gaining popularity in computational finance as they are competitive alternatives to Monte Carlo methods that can accelerate numerical accuracy. This paper develops a new approach for reducing the effective dimension combined with a randomized QMC method. A distinctive feature of the proposed approach is its sample-based transformation that enables us to choose a flexible manipulation via regression. In the proposed approach, the first step is to perform a regression using the samples to estimate the parameters of the regression model. An optimal transformation is proposed based on the regression result to minimize the effective dimension. An advantage of this approach is that adopting a statistical approach allows greater flexibility in selecting the regression model. In addition to a linear model, this paper proposes a dimension reduction method based on a linear-quadratic model for regression. In numerical experiments, we focus on pricing different types of exotic options to test the effectiveness of the proposed approach. The numerical results show that different regression models are chosen depending on the underlying risk process and the type of derivative securities. In particular, we show several examples where the proposed method works while existing dimension reductions are ineffective.

Seismic retrofit optimization of water distribution systems based on the reduction of uncertain damage scenarios

Pre-event interventions, such as retrofit of critical components, are effective strategies to improve the seismic performance and resilience of infrastructure systems. This paper proposes a damage scenario reduction method aims to reduce the computation burden within the seismic retrofit optimization for water distribution systems (WDSs) considering the uncertain damages of pipelines. The optimization framework consists of a multi-objective optimization model for WDS pipeline retrofit and a sequential optimization model for the reduction of uncertain damage scenarios. The multi-objective optimization model is established to minimize the retrofit cost and maximize the seismic performance of WDS. The sequential optimization model, namely weighted optimization-based probabilistic scenarios (OPS + W), was designed to mitigate the computational burden of the multi-objective optimization model, which arises from considering pipeline damage uncertainties through quasi-Monte Carlo (QMC) sampling that generates numerous damage scenarios. The OPS + W selects a few representative damage scenarios from numerous QMC scenarios, aiming to minimize the error of pipeline damage probabilities counted by the representative scenarios, where the error is weighted by the importance of pipeline. The effectiveness of OPS + W is verified through its application in seismic retrofit optimization of two WDS cases and compared with other scenario reduction methods. Application results reveal that the retrofit strategies derived from OPS + W have the highest similarity to those obtained from QMC scenarios while requiring only 3 %–10 % of the computation time compared to using QMC scenarios.

Groundwater LNAPL Contamination Source Identification Based on Stacking Ensemble Surrogate Model

Groundwater LNAPL (Light Non-Aqueous Phase Liquid) contamination source identification (GLCSI) is essential for effective remediation and risk assessment. Addressing the GLCSI problem often involves numerous repetitive forward simulations, which are computationally expensive and time-consuming. Establishing a surrogate model for the simulation model is an effective way to overcome this challenge. However, how to obtain high-quality samples for training the surrogate model and which method should be used to develop the surrogate model with higher accuracy remain important questions to explore. To this end, this paper innovatively adopted the quasi-Monte Carlo (QMC) method to sample from the prior space of unknown variables. Then, this paper established a variety of individual machine learning surrogate models, respectively, and screened three with higher training accuracy among them as the base-learning models (BLMs). The Stacking ensemble framework was utilized to integrate the three BLMs to establish the ensemble surrogate model for the groundwater LNAPL multiphase flow numerical simulation model. Finally, a hypothetical case of groundwater LNAPL contamination was designed. After evaluating the accuracy of the Stacking ensemble surrogate model, the differential evolution Markov chain (DE-MC) algorithm was applied to jointly identify information on groundwater LNAPL contamination source and key hydrogeological parameters. The results of this study demonstrated the following: (1) Employing the QMC method to sample from the prior space resulted in more uniformly distributed and representative samples, which improved the quality of the training data. (2) The developed Stacking ensemble surrogate model had a higher accuracy than any individual surrogate model, with an average R2 of 0.995, and reduced the computational burden by 99.56% compared to the inversion process based on the simulation model. (3) The application of the DE-MC algorithm effectively solved the GLCSI problem, and the mean relative error of the identification results of unknown variables was less than 5%.

Markov chain quasi-Monte Carlo method for forecasting fire hotspots in Sarawak, Malaysia.

Stochastic modeling approaches have attracted many researchers to the field. However, fire hotspot detection suffers from not using a Markov chain quasi-Monte Carlo (MCQMC) as a forecasting methodology. This paper proposes improvements to the computational time by combining the strengths of the Markov chain Monte Carlo (MCMC) and quasi-Monte Carlo (QMC) methods. The proposed method can lead to more precise and stable results, particularly in problems with high-dimensional integration or complex probability distributions. The proposed method is applied to a case study of fire hotspot detection in Sarawak, Malaysia. The outcome of this study reveals that the MCQMC method is more computationally efficient, taking only 0.0746seconds compared to MCMC's 0.0914seconds and QMC's 0.0994seconds. It is shown that the best option derived by the proposed method is effective in predicting fire hotspots and providing quick solutions to protect the environment and communities from forest fires.

Multilevel Monte Carlo Methods for Stochastic Convection–Diffusion Eigenvalue Problems

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection–diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov–Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

Efficient Methods for Reliability Sensitivity Analysis of Distribution Parameters and Their Application

To efficiently evaluate the influence of the distribution parameters of the input variables on the failure probability of engineering structures and improve the reliability and safety of engineering structures in a targeted manner, new methods for the global reliability sensitivity analysis (RSA) of distribution parameters are proposed in this study based on the cubature formula (CF), surrogate sampling probability density function (SSPDF), and quasi-Monte Carlo (QMC) method. By introducing CF, the proposed methods can effectively improve the computational efficiency of the nested expectation and variance operators in the reliability sensitivity indices of the distribution parameters. Based on the concept of SSPDF, a surrogate importance sampling probability density function was developed. This not only overcomes the problem of the computational effort of propagating parameter uncertainty to the failure probability function (FPF), which depends on the dimensionality of the parameters; it also further improves the efficiency of the RSA of the parameters in the case of a small failure probability. Finally, by incorporating the idea of the QMC method, the process of calculating the reliability sensitivity indices of the parameters is reduced from a double-loop to a single-loop one. Three engineering examples are used in this study to demonstrate the efficiency and accuracy of the new algorithms.

Parabolic PDE-constrained optimal control under uncertainty with entropic risk measure using quasi-Monte Carlo integration

We study the application of a tailored quasi-Monte Carlo (QMC) method to a class of optimal control problems subject to parabolic partial differential equation (PDE) constraints under uncertainty: the state in our setting is the solution of a parabolic PDE with a random thermal diffusion coefficient, steered by a control function. To account for the presence of uncertainty in the optimal control problem, the objective function is composed with a risk measure. We focus on two risk measures, both involving high-dimensional integrals over the stochastic variables: the expected value and the (nonlinear) entropic risk measure. The high-dimensional integrals are computed numerically using specially designed QMC methods and, under moderate assumptions on the input random field, the error rate is shown to be essentially linear, independently of the stochastic dimension of the problem—and thereby superior to ordinary Monte Carlo methods. Numerical results demonstrate the effectiveness of our method.

Probabilistic analysis of bolted joint anti-loosening under cyclic transverse load

The loosening of bolted joints under cyclic transverse load is a crucial failure mode, and the uncertainty of the relevant parameters of the bolted joint will affect the loosening state. However, there are few bolt-loosening probability analysis and evaluation methods based on parameter uncertainty in previous studies. In this paper, different loosening states are analyzed based on the relationship between pitch torque, thread friction torque, and bearing friction torque. Considering the uncertainty of process parameters, a probabilistic analysis method of the bolted joint anti-loosening under cyclic transverse load is proposed. In the proposed method, the finite element model verified by the fastening experiment is used to verify the accuracy of the bolted joints’ loosening criterion. Then, the loosening criterion is used as the Limit State Function (LSF). According to the LSF, an Adaptive Kriging-Kernel Density Importance Sampling (AK-KDIS) method is proposed for solving the anti-loosening reliability. This method combines adaptive Kriging model with the importance sampling method, effectively reducing the number of sampling samples and the number of LSF calls. The classical and robust Quasi-Monte Carlo (QMC) method is also used to evaluate the reliability of anti-loosening to verify the accuracy of the proposed AK-KDIS method. Finally, the practical application of probabilistic analysis of anti-loosening is presented, and the effectiveness of the proposed method is analyzed in detail. The results indicate that the AK-KDIS method has higher efficiency while ensuring calculation accuracy compared with the QMC method.

A search for short-period Tausworthe generators over F b with application to Markov chain quasi-Monte Carlo

A one-dimensional sequence u 0 , u 1 , u 2 , … ∈ [ 0 , 1 ) is said to be completely uniformly distributed (CUD) if overlapping s-blocks ( u i , u i + 1 , … , u i + s − 1 ) , i = 0 , 1 , 2 , … , are uniformly distributed for every dimension s ≥ 1 . This concept naturally arises in Markov chain quasi-Monte Carlo (QMC). However, the definition of CUD sequences is not constructive, and thus there remains the problem of how to implement the Markov chain QMC algorithm in practice. Harase [A table of short-period Tausworthe generators for Markov chain quasi-Monte Carlo. J Comput Appl Math. 2021;384:Paper No. 113136, 12.] focussed on the t-value, which is a measure of uniformity widely used in the study of QMC, and implemented short-period Tausworthe generators (i.e. linear feedback shift register generators) over the two-element field F 2 that approximate CUD sequences by running for the entire period. In this paper, we generalize a search algorithm over F 2 to that over arbitrary finite fields F b with b elements and conduct a search for Tausworthe generators over F b with t-values zero (i.e. optimal) for dimension s = 3 and small for s ≥ 4 , especially in the case where b = 3, 4, and 5. We provide a parameter table of Tausworthe generators over F 4 , and report a comparison between our new generators over F 4 and existing generators over F 2 in numerical examples using Markov chain QMC.

Stratified sample tiling

The paper introduces a practical method for the construction of large-scale point sets for analysis of computer models. The constructed experimental design is useful for (probabilistic) integration, construction of approximation or a screening. The essence of the presented approach is the stratification of the design domain into an orthogonal grid of substrata and a subsequent tiling with tiles of points. If optimized, such tiles experience a major reduction of the number of degrees of freedom in the optimization process. That way, optimal or near-optimal point patterns can be feasibly identified and are further utilized for construction of larger point sets, thanks to the idea of self-similarity and structured space stratification The space-filling properties of the resulting point sets may be further enhanced by various “scrambling” strategies, which may remove the undesired sample collapsibility achieved via regular tiling. The performance of the constructed point sets is compared to Quasi Monte Carlo (QMC), Randomized Quasi Monte Carlo (RQMC) sequences, which are still today considered by engineers and even scientists as choices for variance reduction of numerical integration Further, the mentioned sampling strategies are compared in the terms of robustness when integrating a multivariate function with a localized feature. It is concluded that the proposed sampling approach reaches a superior performance in numerical integration and identification of function extremes as compared to sampling methods used by practicing researchers and engineers. Additionally, the reader is supplied with the open-access, ready-to-use implementation of the presented algorithm named SampleTiler.

Uncertainty quantification for random domains using periodic random variables

We consider uncertainty quantification for the Poisson problem subject to domain uncertainty. For the stochastic parameterization of the random domain, we use the model recently introduced by Kaarnioja et al. (SIAM J. Numer. Anal., 2020) in which a countably infinite number of independent random variables enter the random field as periodic functions. We develop lattice quasi-Monte Carlo (QMC) cubature rules for computing the expected value of the solution to the Poisson problem subject to domain uncertainty. These QMC rules can be shown to exhibit higher order cubature convergence rates permitted by the periodic setting independently of the stochastic dimension of the problem. In addition, we present a complete error analysis for the problem by taking into account the approximation errors incurred by truncating the input random field to a finite number of terms and discretizing the spatial domain using finite elements. The paper concludes with numerical experiments demonstrating the theoretical error estimates.

Efficient risk estimation via nested multilevel quasi-Monte Carlo simulation

We consider the problem of estimating the probability of a large loss from a financial portfolio, where the future loss is expressed as a conditional expectation. Since the conditional expectation is intractable in most cases, one may resort to nested simulation. To reduce the complexity of nested simulation, we present an improved multilevel Monte Carlo (MLMC) method by using quasi-Monte Carlo (QMC) to estimate the portfolio loss in each financial scenario generated via Monte Carlo. We prove that using QMC can accelerate the convergence rates in both the crude nested simulation and the multilevel nested simulation. Under certain conditions, the complexity of the proposed MLMC method can be reduced to O(ϵ−2(logϵ)2). On the other hand, we find that using QMC in MLMC encounters a high-kurtosis phenomenon due to the existence of indicator functions. To remedy this, we propose a smoothed method which uses logistic sigmoid functions to approximate indicator functions. Numerical results show that the optimal MLMC complexity O(ϵ−2) is almost attained even in moderate high dimensions.

Uncertainties of band sound levels when estimated from monochromatic outdoor sound propagation calculations

Octave band and third octave bands sound levels are widely used in environmental acoustics, and most numerical models for outdoor sound propagation estimate them through the combination of several monochromatic calculations. As some of these models can be very time consuming when long range distance calculations are needed, the number of monochromatic calculations performed for each band is often necessarily very limited. That induces uncertainties on the estimate of the noise level of the band that can be large when the spectrum amplitude varies widely inside a band. This is the case for example when interferences due to ground effects occur. The paper presents estimates of these uncertainties depending on frequency band and on source-receiver locations. The uncertainties are evaluated thanks to a quasi Monte-Carlo technique from many calculations using a Rudnick propagation model and a Latin Hypercube Sampling numerical experimental design for input parameters (source-receiver locations). The results provide practical indications on the minimum number of monochromatic calculations to do per band, in order not to exceed a given uncertainty level.

Application of high-credible statistical results calculation scheme based on least squares Quasi-Monte Carlo method in multimodal stochastic problems

Solving stochastic problems based on multimodal distributions with high accuracy and efficiency is the focus of current research in the stochastic field. However, since exact statistical results in practical engineering problems are unknown, how to efficiently obtain statistical results with high credibility is of great importance. Therefore, a high-credible statistical results calculation scheme is proposed in this paper. First, the multimodal model is transformed into a weighted sum form of multiple unimodal models by Gaussian mixture model. Then, a high-credible statistical results calculation scheme is designed based on the reusability of low-discrepancy sequences in the quasi-Monte Carlo (QMC) method. The scheme calculates the coefficient of variation for the random response statistical results of the unimodal model across various sample ranges. When the coefficient of variation is less than the tolerance error, which indicates low fluctuation, the statistical results are deemed to have converged to the exact solution. In order to further accelerate the calculation efficiency of the proposed scheme, this paper proposes a least squares quasi-Monte Carlo (LSQMC) method to yield high accuracy statistical moments. Compared with the averageness weights in the QMC method, the non-averageness weights obtained by LSQMC method can more effectively reflect the numerical Characteristics of the random variables and the location distribution of the sampling points in the probability space. Finally, numerical examples have demonstrated that the high-credible statistical results calculation scheme based on the LSQMC method can obtain the high accuracy statistical moments with high efficiency in multimodal stochastic problems.

The fast reduced QMC matrix–vector product

We study the approximation of integrals of the form ∫Df(x⊤A)dμ(x), where the measure μ is of product form and A is a matrix, by quasi-Monte Carlo (QMC) rules N−1∑k=0N−1f(xk⊤A). We are interested in cases where the main computational cost in the approximation arises from calculating the products xk⊤A. We design QMC rules for which the computation of xk⊤A, k=0,1,…,N−1, can be done in a fast way, and for which the approximation error of the QMC rule is similar to the standard QMC error. We do not require that the matrix A has any particular structure.Problems of this form arise in some important applications in statistics and uncertainty quantification. For instance, this approach can be used when approximating the expected value of some function with a multivariate normal random variable with some given covariance matrix, or when approximating the expected value of the solution of a PDE with random coefficients.The speed-up of the computation time of our approach is sometimes better and sometimes worse than the fast QMC matrix–vector product from [Josef Dick, Frances Y. Kuo, Quoc T. Le Gia, and Christoph Schwab, Fast QMC Matrix-Vector Multiplication, SIAM J. Sci. Comput. 37 (2015), no. 3, A1436–A1450]. As in that paper, our approach applies to lattice point sets and polynomial lattice point sets, but also applies to digital nets (we are currently not aware of any approach which allows one to apply the fast QMC matrix–vector paper from the aforementioned paper of Dick, Kuo, Le Gia, and Schwab to digital nets).The method in this paper does not make use of the fast Fourier transform, instead we use repeated values in the quadrature points to derive a significant reduction in the computation time. Such a situation naturally arises from the reduced CBC construction of lattice rules and polynomial lattice rules. The reduced CBC construction has been shown to reduce the computation time for the CBC construction. Here we show that it can additionally be used to also reduce the computation time of the underlying QMC rule. One advantage of the present approach is that it can be combined with random (digital) shifts, whereas this does not apply to the fast QMC matrix–vector product from the earlier paper of Dick, Kuo, Le Gia, and Schwab.

On the reconstruction of functions from values at subsampled quadrature points

This paper is concerned with function reconstruction from samples. The sampling points used in several approaches are (1) structured points connected with fast algorithms or (2) unstructured points coming from, e.g., an initial random draw to achieve an improved information complexity. We connect both approaches and propose a subsampling of structured points in an offline step. In particular, we start with quasi-Monte Carlo (QMC) points with inherent structure and stable L 2 L_2 reconstruction properties. The subsampling procedure consists of a computationally inexpensive random step followed by a deterministic procedure to further reduce the number of points while keeping its information. In these points functions (belonging to a reproducing kernel Hilbert space of bounded functions) will be sampled and reconstructed from whilst achieving state of the art error decay. Our method is dimension-independent and is applicable as soon as we know some initial quadrature points. We apply our general findings on the d d -dimensional torus to subsample rank-1 lattices, where it is known that full rank-1 lattices lose half the optimal order of convergence (expressed in terms of the size of the lattice). In contrast to that, our subsampled version regains the optimal rate since many of the lattice points are not needed. Moreover, we utilize fast and memory efficient Fourier algorithms in order to compute the approximation. Numerical experiments in several dimensions support our findings.

An efficient quasi-Monte Carlo method with forced fixed detection for photon scatter simulation in CT.

Detected scattered photons can cause cupping and streak artifacts, significantly degrading the quality of CT images. For fast and accurate estimation of scatter intensities resulting from photon interactions with a phantom, we first transform the path probability of photons interacting with the phantom into a high-dimensional integral. Secondly, we develope a new efficient algorithm called gQMCFFD, which combines graphics processing unit(GPU)-based quasi-Monte Carlo (QMC) with forced fixed detection to approximate this integral. QMC uses low discrepancy sequences for simulation and is deterministic versions of Monte Carlo. Numerical experiments show that the results are in excellent agreement and the efficiency improvement factors are 4 ∼ 46 times in all simulations by gQMCFFD with comparison to GPU-based Monte Carlo methods. And by combining gQMCFFD with sparse matrix method, the simulation time is reduced to 2 seconds in a single projection angle and the relative difference is 3.53%.

Bypassing the quadrature exactness assumption of hyperinterpolation on the sphere

This paper focuses on the approximation of continuous functions on the unit sphere by spherical polynomials of degree n via hyperinterpolation. Hyperinterpolation of degree n is a discrete approximation of the L2-orthogonal projection of the same degree with its Fourier coefficients evaluated by a positive-weight quadrature rule that exactly integrates all spherical polynomials of degree at most 2n. This paper aims to bypass this quadrature exactness assumption by replacing it with the Marcinkiewicz–Zygmund property proposed in a previous paper. Consequently, hyperinterpolation can be constructed by a positive-weight quadrature rule (not necessarily with quadrature exactness). This scheme is referred to as unfettered hyperinterpolation. This paper provides a reasonable error estimate for unfettered hyperinterpolation. The error estimate generally consists of two terms: a term representing the error estimate of the original hyperinterpolation of full quadrature exactness and another introduced as compensation for the loss of exactness degrees. A guide to controlling the newly introduced term in practice is provided. In particular, if the quadrature points form a quasi-Monte Carlo (QMC) design, then there is a refined error estimate. Numerical experiments verify the error estimates and the practical guide.