quasi-Monte Carlo Integration Research Articles

Differentiable Owen Scrambling

Quasi-Monte Carlo integration is at the core of rendering. This technique estimates the value of an integral by evaluating the integrand at well-chosen sample locations. These sample points are designed to cover the domain as uniformly as possible to achieve better convergence rates than purely random points. Deterministic low-discrepancy sequences have been shown to outperform many competitors by guaranteeing good uniformity as measured by the so-called discrepancy metric, and, indirectly, by an integer t value relating the number of points falling into each domain stratum with the stratum area (lower t is better). To achieve randomness, scrambling techniques produce multiple realizations preserving the t value, making the construction stochastic. Among them, Owen scrambling is a popular approach that recursively permutes intervals for each dimension. However, relying on permutation trees makes it incompatible with smooth optimization frameworks. We present a differentiable Owen scrambling that regularizes permutations. We show that it can effectively be used with automatic differentiation tools for optimizing low-discrepancy sequences to improve metrics such as optimal transport uniformity, integration error, designed power spectra or projective properties, while maintaining their initial t -value as guaranteed by Owen scrambling. In some rendering settings, we show that our optimized sequences improve the rendering error.

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.

Analyzing the Robustness of Complex Networks with Attack Success Rate.

Analyzing the robustness of networks against random failures or malicious attacks is a critical research issue in network science, as it contributes to enhancing the robustness of beneficial networks and effectively dismantling harmful ones. Most studies commonly neglect the impact of the attack success rate (ASR) and assume that attacks on the network will always be successful. However, in real-world scenarios, attacks may not always succeed. This paper proposes a novel robustness measure called Robustness-ASR (RASR), which utilizes mathematical expectations to assess network robustness when considering the ASR of each node. To efficiently compute the RASR for large-scale networks, a parallel algorithm named PRQMC is presented, which leverages randomized quasi-Monte Carlo integration to approximate the RASR with a faster convergence rate. Additionally, a new attack strategy named HBnnsAGP is introduced to better assess the lower bound of network RASR. Finally, the experimental results on six representative real-world complex networks demonstrate the effectiveness of the proposed methods compared with the state-of-the-art baselines.

Quasi-Random Discrete Ordinates Method to Radiative Transfer Equation with Linear Anisotropic Scattering

The modeling of energy transport via radiative transfer is important to many practical high temperature engineering applications. Furthermore, the computation of solutions to the radiative transfer equation (RTE) plays a fundamental part in it. The quasi-random discrete ordinates method was developed as an alternative to mitigate the ray effect found in the classical discrete ordinates method solutions. The former method was originally developed for transport problems with isotropic scattering and it is here extended and tested to problems with linear anisotropic scattering. Its main idea is to approximate the integral term of the RTE by a quasi-Monte Carlo integration. The discrete system of differential equations arising from it can be solved by a variety of classical discretization methods, here computed with a SUPG finite element scheme. The novel developments are tested for selected manufactured solutions. The achieved good results indicate the potential of the novel method to be applied to the solution of radiative transfer problems with anisotropic scattering.

On a Reduced Component-by-Component Digit-by-Digit Construction of Lattice Point Sets

Abstract In this paper, we study an efficient algorithm for constructing point sets underlying quasi-Monte Carlo integration rules for weighted Korobov classes. The algorithm presented is a reduced fast component-by-component digit-by-digit (CBC-DBD) algorithm, which is useful for situations where the weights in the function space show a sufficiently fast decay. The advantage of the algorithm presented here is that the computational effort can be independent of the dimension of the integration problem to be treated if suitable assumptions on the integrand are met. By considering a reduced digit-by-digit construction, we allow an integration algorithm to be less precise with respect to the number of bits in those components of the problem that are considered less important. The new reduced CBC-DBD algorithm is designed to work for the construction of lattice point sets, and the corresponding integration rules (so-called lattice rules) can be used to treat functions in different kinds of function spaces. We show that the integration rules constructed by our algorithm satisfy error bounds of almost optimal convergence order. Furthermore, we give details on an efficient implementation such that we obtain a considerable speed-up of a previously known CBC-DBD algorithm that has been studied in the paper Digit-by-digit and component-by-component constructions of lattice rules for periodic functions with unknown smoothness by Ebert, Kritzer, Nuyens, and Osisiogu, published in the Journal of Complexity in 2021. This improvement is illustrated by numerical results.

Random Effects Multinomial Processing Tree Models: A Maximum Likelihood Approach

The present article proposes and evaluates marginal maximum likelihood (ML) estimation methods for hierarchical multinomial processing tree (MPT) models with random and fixed effects. We assume that an identifiable MPT model with S parameters holds for each participant. Of these S parameters, R parameters are assumed to vary randomly between participants, and the remaining S-R parameters are assumed to be fixed. We also propose an extended version of the model that includes effects of covariates on MPT model parameters. Because the likelihood functions of both versions of the model are too complex to be tractable, we propose three numerical methods to approximate the integrals that occur in the likelihood function, namely, the Laplace approximation (LA), adaptive Gauss–Hermite quadrature (AGHQ), and Quasi Monte Carlo (QMC) integration. We compare these three methods in a simulation study and show that AGHQ performs well in terms of both bias and coverage rate. QMC also performs well but the number of responses per participant must be sufficiently large. In contrast, LA fails quite often due to undefined standard errors. We also suggest ML-based methods to test the goodness of fit and to compare models taking model complexity into account. The article closes with an illustrative empirical application and an outlook on possible extensions and future applications of the proposed ML approach.

A Mixed-Effects Model in Which the Parameters of the Autocorrelated Error Structure Can Differ between Individuals

Research in psychology has seen a rapid increase in the usage of experience sampling methods and daily diary methods. The data that result from using these methods are typically analyzed with a mixed-effects or a multilevel model because it allows testing hypotheses about the time course of the longitudinally assessed variable or the influence of time-varying predictors in a simple way. Here, we describe an extension of this model that does not only allow to include random effects for the mean structure but also for the residual variance, for the parameter of an autoregressive process of order 1 and/or the parameter of a moving average process of order 1. After we have introduced this extension, we show how to estimate the parameters with maximum likelihood. Because the likelihood function contains complex integrals, we suggest using adaptive Gauss-Hermite quadrature and Quasi-Monte Carlo integration to approximate it. We illustrate the models using a real data example and also report the results of a small simulation study in which the two integral approximation methods are compared.

Analysis of Preintegration Followed by Quasi–Monte Carlo Integration for Distribution Functions and Densities

Analysis of Preintegration Followed by Quasi–Monte Carlo Integration for Distribution Functions and Densities

Parton distributions need representative sampling

In global QCD fits of parton distribution functions (PDFs), a large part of the estimated uncertainty on the PDFs originates from the choices of parametric functional forms and fitting methodology. We argue that these types of uncertainties can be underestimated with common PDF ensembles in high-stake measurements at the Large Hadron Collider and Tevatron. A fruitful approach to quantify these uncertainties is to view them as arising from sampling of allowed PDF solutions in a multidimensional parametric space. This approach applies powerful insights gained in recent statistical studies of large-scale population surveys and quasi-Monte Carlo integration methods. In particular, PDF fits may be affected by the big data paradox, which stipulates that more experimental data do not automatically raise the accuracy of PDFs---close attention to the data quality and sampling of possible PDF solutions is as essential. To test if the sampling of the PDF uncertainty of an experimental observable is truly representative of all acceptable solutions, we introduce a technique (``a hopscotch scan'') based on a combination of parameter scans and stochastic sampling. With this technique, we examine the PDF uncertainty on key LHC cross sections at 13 TeV obtained with the public NNPDF4.0 fitting code, while accounting for the likelihood distribution. We show that the uncertainties on the charm distribution at a large momentum fraction $x$ and gluon PDF at small $x$ are enlarged. In PDF ensembles obtained in the analytic minimization (Hessian) formalism, the tolerance on the PDF uncertainty must be based on sufficiently complete sampling of PDF functional forms and choices of the experiments.

Improved bounds on the gain coefficients for digital nets in prime power base

We study randomized quasi-Monte Carlo integration by scrambled nets. The scrambled net quadrature has long gained its popularity because it is an unbiased estimator of the true integral, allows for a practical error estimation, achieves a high order decay of the variance for smooth functions, and works even for Lp-functions with any p≥1. The variance of the scrambled net quadrature for L2-functions can be evaluated through the set of the so-called gain coefficients.In this paper, based on the system of Walsh functions and the concept of dual nets, we provide improved upper bounds on the gain coefficients for digital nets in general prime power base. Our results explain the known bound by Owen (1997) for Faure sequences, the recently improved bound by Pan and Owen (2022) for digital nets in base 2 (including Sobol' sequences as a special case), and their finding that all the nonzero gain coefficients for digital nets in base 2 must be powers of two, all in a unified way.

Improving the accuracy of estimators for the two-point correlation function

Aims. We show how to increase the accuracy of estimates of the two-point correlation function without sacrificing efficiency. Methods. We quantify the error of the pair-counts and of the Landy & Szalay estimator by comparing them with exact reference values. The standard method, using random point sets, is compared to geometrically motivated estimators and estimators using quasi-Monte Carlo integration. Results. In the standard method, the error scales proportionally to 1/√Nr, with Nr being the number of random points. In our improved methods, the error scales almost proportionally to 1/Nq, where Nq is the number of points from a low-discrepancy sequence. We study the run times of the new estimator in comparison to those of the standard estimator, keeping the same level of accuracy. For the considered case, we always see a speedup ranging from 50% up to a factor of several thousand. We also discuss how to apply these improved estimators to incompletely sampled galaxy catalogues.

Construction-Free Median Quasi-Monte Carlo Rules for Function Spaces with Unspecified Smoothness and General Weights

We study quasi-Monte Carlo (QMC) integration of smooth functions defined over the multi-dimensional unit cube. Inspired by a recent work of Pan and Owen, we study a new construction-free median QMC rule which can exploit the smoothness and the weights of function spaces adaptively. For weighted Korobov spaces, we draw a sample of $r$ independent generating vectors of rank-1 lattice rules, compute the integral estimate for each, and approximate the true integral by the median of these $r$ estimates. For weighted Sobolev spaces, we use the same approach but with the rank-1 lattice rules replaced by high-order polynomial lattice rules. A major advantage over the existing approaches is that we do not need to construct good generating vectors by a computer search algorithm, while our median QMC rule achieves almost the optimal worst-case error rate for the respective function space with any smoothness and weights, with a probability that converges to 1 exponentially fast as $r$ increases. Numerical experiments illustrate and support our theoretical findings.

MatBuilder

Many applications ranging from quasi-Monte Carlo integration over optimal control to neural networks benefit from high-dimensional, highly uniform samples. In the case of computer graphics, and more particularly in rendering, despite the need for uniformity, several sub-problems expose a low-dimensional structure. In this context, mastering sampling uniformity over projections while preserving high-dimensional uniformity has been intrinsically challenging. This difficulty may explain the relatively small number of mathematical constructions for such samplers. We propose a novel approach by showing that uniformity constraints can be expressed as an integer linear program that results in a sampler with the desired properties. As it turns out, complex constraints are easy to describe by means of stratification and sequence properties of digital nets. Formalized using generator matrix determinants, our new MatBuilder software solves the set of constraints by iterating the linear integer program solver in a greedy fashion to compute a problem-specific set of generator matrices that can be used as a drop-in replacement in the popular digital net samplers. The samplers created by MatBuilder achieve the uniformity of classic low discrepancy sequences. More importantly, we demonstrate the benefit of the unprecedented versatility of our constraint approach with respect to low-dimensional problem structure for several applications.

Adaptive Quasi-Monte Carlo Finite Element Methods for Parametric Elliptic PDEs

We introduce novel adaptive methods to approximate moments of solutions of partial differential Equations (PDEs) with uncertain parametric inputs. A typical problem in Uncertainty Quantification is the approximation of the expected values of quantities of interest of the solution, which requires the efficient numerical approximation of high-dimensional integrals. We perform this task by a class of deterministic quasi-Monte Carlo integration rules derived from Polynomial lattices, that allows to control a-posteriori the integration error without querying the governing PDE and does not incur the curse of dimensionality. Based on an abstract formulation of adaptive finite element methods (AFEM) for deterministic problems, we infer convergence of the combined adaptive algorithms in the parameter and physical space. We propose a selection of examples of PDEs admissible for these algorithms. Finally, we present numerical evidence of convergence for a model diffusion PDE.

Cascaded Sobol' sampling

Rendering quality is largely influenced by the samplers used in Monte Carlo integration. Important factors include sample uniformity (e.g., low discrepancy) in the high-dimensional integration domain, sample uniformity in lower-dimensional projections, and lack of dominant structures that could result in aliasing artifacts. A widely used and successful construction is the Sobol' sequence that guarantees good high-dimensional uniformity and consequently results in faster convergence of quasi-Monte Carlo integration. We show that this sequence exhibits low uniformity and dominant structures in low-dimensional projections. These structures impair quality in the context of rendering, as they precisely occur in the 2-dimensional projections used for sampling light sources, reflectance functions, or the camera lens or sensor. We propose a new cascaded construction, which, despite dropping the sequential aspect of Sobol' samples, produces point sets exhibiting provably perfect dyadic partitioning (and therefore, excellent uniformity) in consecutive 2-dimensional projections, while preserving good high-dimensional uniformity. By optimizing the initialization parameters and performing Owen scrambling at finer levels of binary representations, we further improve over Sobol's integration convergence rate. Our method does not incur any overhead as compared to the generation of the Sobol' sequence, is compatible with Owen scrambling and can be used in rendering applications.

Weighted integration over a hyperrectangle based on digital nets and sequences

Quasi-Monte Carlo (QMC) methods are equal weight quadrature rules to approximate integrals over the unit cube with respect to the uniform measure. In this paper we discuss QMC integration with respect to general product measures defined on an arbitrary hyperrectangle. We only require that the cumulative distribution function is invertible. We develop a worst-case error bound and study the dependence of the error on the number of points and the dimension for digital nets and sequences as well as polynomial lattice point sets, which are mapped to the domain using the inverse cumulative distribution function. We do not require any smoothness properties of the probability density function and the worst-case error does not depend on the particular choice of density function and its smoothness. The component-by-component construction of polynomial lattice rules is based on a criterion which depends only on the size of the hyperrectangle but is otherwise independent of the product measure.

Walsh functions, scrambled [formula omitted]-nets, and negative covariance: Applying symbolic computation to quasi-Monte Carlo integration

We investigate base b Walsh functions for which the variance of the integral estimator based on a scrambled (0,m,s)-net in base b is less than or equal to that of the Monte-Carlo estimator based on the same number of points. First we compute the Walsh decomposition for the joint probability density function of two distinct points randomly chosen from a scrambled (t,m,s)-net in base b in terms of certain counting numbers and simplify it in the special case t is zero. Using this, we obtain an expression for the covariance of the integral estimator in terms of the Walsh coefficients of the function. Finally, we prove that the covariance of the integral estimator is negative when the Walsh coefficients of the function satisfy a certain decay condition. To do this, we use creative telescoping and recurrence solving algorithms from symbolic computation to find a sign equivalent closed form expression for the covariance term.

QMC integration errors and quasi-asymptotics

Abstract A crude Monte Carlo (MC) method allows to calculate integrals over a d-dimensional cube. As the number N of integration nodes becomes large, the rate of probable error of the MC method decreases as O ⁢ ( 1 / N ) {O(1/\sqrt{N})} . The use of quasi-random points instead of random points in the MC algorithm converts it to the quasi-Monte Carlo (QMC) method. The asymptotic error estimate of QMC integration of d-dimensional functions contains a multiplier 1 / N {1/N} . However, the multiplier ( ln ⁡ N ) d {(\ln N)^{d}} is also a part of the error estimate, which makes it virtually useless. We have proved that, in the general case, the QMC error estimate is not limited to the factor 1 / N {1/N} . However, our numerical experiments show that using quasi-random points of Sobol sequences with N = 2 m {N=2^{m}} with natural m makes the integration error approximately proportional to 1 / N {1/N} . In our numerical experiments, d ≤ 15 {d\leq 15} , and we used N ≤ 2 40 {N\leq 2^{40}} points generated by the SOBOLSEQ16384 code published in 2011. In this code, d ≤ 2 14 {d\leq 2^{14}} and N ≤ 2 63 {N\leq 2^{63}} .

Faster Monte Carlo estimation of joint models for time-to-event and multivariate longitudinal data

Quasi-Monte Carlo (QMC) methods using quasi-random sequences, as opposed to pseudo-random samples, are proposed for use in the joint modelling of time-to-event and multivariate longitudinal data. The QMC integration framework extends the Monte Carlo Expectation Maximisation approaches that are commonly adopted, namely using ordinary and antithetic variates. The motivation of QMC integration is to increase the convergence speed by using nodes that are scattered more uniformly. Through simulation, estimates and computational times are compared and this is followed with an application to a clinical dataset. There is a distinct speed advantage in using QMC methods for small sample sizes and QMC is comparable to the antithetic MC method for moderate sample sizes. The new method is available in an updated version of the R package joineRML.

Domain Uncertainty Quantification in Computational Electromagnetics

We study the numerical approximation of time-harmonic, electromagnetic fields inside a lossy cavity of uncertain geometry. Key assumptions are a possibly high-dimensional parametrization of the uncertain geometry along with a suitable transformation to a fixed, nominal domain. This uncertainty parametrization results in families of countably parametric, Maxwell-like cavity problems that are posed in a single domain, with inhomogeneous coefficients that possess finite, possibly low spatial regularity, but exhibit holomorphic parametric dependence in the differential operator. Our computational scheme is composed of a sparse grid interpolation in the high-dimensional parameter domain and an ${H}{curl}$ -conforming edge element discretization of the parametric problem in the nominal domain. As a stepping-stone in the analysis, we derive a novel Strang-type lemma for Maxwell-like problems in the nominal domain, which is of independent interest. Moreover, we accommodate arbitrary small Sobolev regularity of the electric field and also cover uncertain isotropic constitutive or material laws. The shape holomorphy and edge-element consistency error analysis for the nominal problem are shown to imply convergence rates for multilevel Monte Carlo and for quasi-Monte Carlo integration, as well as sparse grid approximations, in uncertainty quantification for computational electromagnetics. They also imply expression rate estimates for deep ReLU networks of shape-to-solution maps in this setting. Finally, our computational experiments confirm the presented theoretical results.