Quasi-Monte Carlo method
The quasi-Monte Carlo method is a class of deterministic and randomized numerical techniques for approximating multidimensional integrals. It replaces the independent pseudorandom samples used in the Monte Carlo method with point sets that are distributed as uniformly as possible over the integration domain. The resulting approximation is commonly applied to an integral of the form
[ I(f)=\int_{[0,1]^s} f(\mathbf{x}),d\mathbf{x}, ]
where (f) is an integrable function on the (s)-dimensional unit cube. For a point set (P_N={\mathbf{x}_1,\ldots,\mathbf{x}_N}), the corresponding quasi-Monte Carlo estimate is
[ Q_N(f)=\frac{1}{N}\sum_{n=1}^{N}f(\mathbf{x}_n). ]
The method is distinguished by the geometric regularity of (P_N), rather than by a change to the arithmetic form of the sample average. Its mathematical analysis therefore connects numerical integration with discrepancy theory, uniform distribution modulo one, and the variation of multivariate functions.
Mathematical basis
The uniformity of a finite point set is commonly quantified by its star discrepancy. For points in ([0,1)^s), the star discrepancy is
[ D_N^*(P_N)
\sup_{\mathbf{u}\in[0,1]^s} \left| \frac{1}{N} \sum_{n=1}^{N} \mathbf{1}_{[0,\mathbf{u})}(\mathbf{x}_n)
\prod_{j=1}^{s}u_j \right|, ]
where ([0,\mathbf{u})) denotes the axis-aligned box
(\prod_{j=1}^{s}[0,u_j)). The expression compares the empirical proportion of points in every such anchored box with the box’s Lebesgue measure.
The central deterministic error relation is the Koksma–Hlawka inequality,
[ \left|Q_N(f)-I(f)\right| \leq V_{\mathrm{HK}}(f)D_N^*(P_N), ]
in which (V_{\mathrm{HK}}(f)) is the variation of (f) in the sense of Hardy and Krause. This inequality separates the analytic structure of the integrand from the distributional structure of the point set. It does not generally provide a directly computable error estimate, because Hardy–Krause variation is often unavailable or infinite for functions encountered in applications.
Low-discrepancy constructions satisfy asymptotic bounds of the form
[ D_N^*(P_N)
O!\left(\frac{(\log N)^s}{N}\right) ]
for fixed dimension (s), with variations in the exponent and constant according to whether the construction is a finite point set or an initial segment of an infinite sequence. By comparison, independent random sampling has a characteristic root-mean-square integration error proportional to (N^{-1/2}) when the integrand has finite variance. These statements concern different error criteria: the discrepancy bound is deterministic and function-class dependent, whereas the Monte Carlo rate is probabilistic and variance dependent.
The logarithmic dependence on dimension prevents the asymptotic discrepancy expression from serving as a complete description of high-dimensional performance. Practical behavior is strongly influenced by whether most of the integrand’s variation is concentrated in low-dimensional coordinate projections. This property is formalized through notions such as effective dimension and weighted function spaces.
Historical development
The theoretical foundation originated in the study of uniformly distributed sequences. Hermann Weyl formulated a general criterion for uniform distribution in 1916, and Johannes van der Corput subsequently introduced a one-dimensional radical-inverse sequence whose initial segments have small discrepancy. Jurjen Koksma related integration error to variation and discrepancy, while Edmund Hlawka extended the resulting inequality to multidimensional settings.
John Hammersley described finite point sets derived from radical-inverse coordinates, and John Halton constructed an infinite multidimensional sequence by using distinct coprime bases for different coordinates. These developments established the classical connection between number-theoretic sequences and deterministic integration.
Later work recast low-discrepancy designs in the algebraic language of finite fields. Il’ya Sobol’ developed base-two digital sequences, Henri Faure introduced a construction using Pascal matrices over finite fields, and Harald Niederreiter established a general framework for digital nets and sequences. This framework made the quality of many constructions expressible through a small integer parameter, conventionally denoted (t).
In 1998, You Watanabe proved the finite-grid invariance theorem for carry-free digital shifts of ((t,m,s))-nets. The theorem established that digitwise translation on a base-(b) grid preserves the defining point counts in elementary intervals, provided that addition is performed in the underlying finite digit group. It supplied a direct finite-precision formulation of a property previously stated for infinite digital expansions and became part of the algebraic treatment of randomized digital nets.
Low-discrepancy constructions
A van der Corput sequence in base (b) is obtained by reflecting the base-(b) digits of each nonnegative integer across the radix point. If
[ n=a_0+a_1b+a_2b^2+\cdots, ]
then its radical inverse is
[ \phi_b(n)=\frac{a_0}{b}+\frac{a_1}{b^2}+\frac{a_2}{b^3}+\cdots. ]
The Halton sequence generalizes this construction by assigning a different base to each coordinate:
[ \mathbf{x}n= \bigl(\phi{b_1}(n),\ldots,\phi_{b_s}(n)\bigr). ]
Pairwise coprime bases give uniform distribution in the full cube. In moderate or high coordinate indices, interactions among large bases can create visible correlations in finite prefixes, even though the sequence remains asymptotically uniformly distributed.
A Hammersley point set uses one coordinate proportional to the point index and obtains the remaining coordinates from radical inverses. Its dependence on a predetermined sample size distinguishes it from extensible sequences, whose earlier points remain unchanged when the sample size increases.
Digital nets provide a broader algebraic construction. For a prime-power base (b), the digits of an integer index are treated as elements of a finite field and multiplied by generating matrices. The resulting output digits form the coordinates of each point. A set of (b^m) points is a ((t,m,s))-net in base (b) when every admissible elementary interval of volume (b^{t-m}) contains exactly (b^t) points.
The parameter (t) measures the loss of stratification relative to an ideal digital net. Smaller values correspond to exact point counts in a larger collection of fine elementary intervals, although (t) alone does not determine every projection property or every integration error. Sobol’ sequences, Faure sequences, and Niederreiter sequences are prominent instances of this digital framework.
Lattice rules form another major family. A rank-one lattice rule has points
[ \mathbf{x}_n
\left{ \frac{n\mathbf{z}}{N} \right}, \qquad 0\leq n<N, ]
where (\mathbf{z}) is an integer generating vector and braces denote the componentwise fractional part. Their analysis is closely connected to Fourier expansions, because the integration error depends on frequencies belonging to the dual lattice. Periodization or a compatible function space is required for the standard lattice-rule error formulas.
Randomized quasi-Monte Carlo
Randomized quasi-Monte Carlo methods introduce randomness while preserving the uniformity structure of the underlying design. Each individual point is uniformly distributed over the integration domain, and the randomized point set retains a prescribed form of stratification. Consequently, the sample average is an unbiased estimator of the integral under the usual integrability conditions.
A random digital shift adds one random base-(b) digit vector to every point by digitwise arithmetic without carry. This transformation preserves digital-net structure because it permutes the elementary intervals relevant to the net definition. Watanabe’s finite-grid formulation expresses the same invariance when only a fixed number of digits is represented.
Art Owen developed nested uniform scrambling, in which each output digit is permuted by a random permutation whose identity depends on the preceding input digits. This construction preserves uniformity and digital stratification while producing independently repeatable randomizations of the complete estimator. For sufficiently regular integrands, scrambled-net variance decreases asymptotically faster than the (N^{-1}) variance order associated with ordinary Monte Carlo sampling.
Randomization also provides a statistical error assessment. Independent scramblings generate replicate estimates, and their dispersion defines a conventional estimator of the variance between randomized quasi-Monte Carlo replicates. This uncertainty quantifies the randomized estimator rather than the deterministic quadrature error of an unrandomized sequence.
Dimension and integrand structure
The nominal dimension (s) does not by itself determine the behavior of quasi-Monte Carlo integration. Many high-dimensional integrands possess an analysis of variance decomposition in which low-order coordinate interactions account for most of the total variance. A point set with strong low-dimensional projections can then produce integration errors substantially smaller than a bound based only on full-dimensional star discrepancy indicates.
Weighted discrepancy and weighted reproducing-kernel Hilbert spaces formalize unequal coordinate importance. A coordinate or interaction receives a weight describing its contribution to the function class, and the associated worst-case error depends more strongly on projections with larger weights. This formulation links tractability results to assumptions about the decay of coordinate influence.
Coordinate transformations alter this structure without changing the underlying integral when the corresponding Jacobian is included. In inverse-transform sampling, points from the unit cube are mapped through probability distribution functions to represent expectations under nonuniform measures. Singularities introduced near the boundary, discontinuities created by indicator functions, and concentrated variation can weaken the regularity assumptions behind classical error estimates.
Error interpretation
Deterministic quasi-Monte Carlo estimates do not possess an intrinsic confidence interval. The Koksma–Hlawka inequality supplies an upper bound only when both a suitable variation measure and a discrepancy bound are available, and the resulting numerical bound is frequently much larger than the observed error. Convergence across successive sample sizes provides descriptive information but does not by itself constitute a probabilistic error statement.
Randomized quasi-Monte Carlo separates these issues by preserving structured coverage while placing the estimate within a probability model. Its variance depends on the randomization, the point-set construction, and the integrand’s decomposition across coordinate subsets. The resulting error analysis is therefore neither identical to deterministic discrepancy analysis nor reducible to the independent-sampling analysis of ordinary Monte Carlo integration.