Recent advances in higher order quasi-Monte Carlo methods

TL;DR

This paper reviews recent theoretical and practical advances in higher order quasi-Monte Carlo methods, highlighting their development, analysis, and applications in complex computational problems.

Contribution

It provides a unified overview of how Walsh analysis has advanced the theory and applications of HoQMC methods since their inception.

Findings

Significant progress in discrepancy and multivariate integration theory.

Successful applications to PDEs with random coefficients.

Enhanced understanding of Walsh analysis in HoQMC development.

Abstract

In this article we review some of recent results on higher order quasi-Monte Carlo (HoQMC) methods. After a seminal work by Dick (2007, 2008) who originally introduced the concept of HoQMC, there have been significant theoretical progresses on HoQMC in terms of discrepancy as well as multivariate numerical integration. Moreover, several successful and promising applications of HoQMC to partial differential equations with random coefficients and Bayesian estimation/inversion problems have been reported recently. In this article we start with standard quasi-Monte Carlo methods based on digital nets and sequences in the sense of Niederreiter, and then move onto their higher order version due to Dick. The Walsh analysis of smooth functions plays a crucial role in developing the theory of HoQMC, and the aim of this article is to provide a unified picture on how the Walsh analysis enables recent developments of HoQMC both for discrepancy and numerical integration.

Equations

$$I(f)={\int }_{[0,1{]}^s}f(\mathbf{x})\mathrm{d}\mathbf{x}.$$

$$I(f;P)=\frac{1}{|P|}\sum _{\mathbf{x}\in P}f(\mathbf{x}),$$

$$\mathbb{E}[I(f;P)]=I(f)\text{and}\mathbb{V}[I(f;P)]=\frac{\mathbb{V}[f]}{|P|},$$

$$|I(f;P)-I(f)|\le V_{\mathrm{HK}}(f)D^{\ast }(P),$$

$$x\oplus y=\sum _{i=1}^{\infty }{\zeta }_i{b}^{-i},\text{where ζi=ξi+ηi(modb)}.$$

$$x_{h,j}=\frac{{\xi }_{1,h,j}}{b}+\frac{{\xi }_{2,h,j}}{b^2}+\cdots +\frac{{\xi }_{n,h,j}}{b^n},$$

$$({\xi }_{1,h,j},{\xi }_{2,h,j},\dots ,{\xi }_{n,h,j})=({\eta }_0,{\eta }_1,\dots ,{\eta }_{m-1})\cdot C_j^{\top }.$$

$$x_{h,j}=\frac{{\xi }_{1,h,j}}{b}+\frac{{\xi }_{2,h,j}}{b^2}+\cdots ,$$

$$({\xi }_{1,h,j},{\xi }_{2,h,j},\dots )=({\eta }_0,{\eta }_1,\dots )\cdot C_j^{\top }.$$

$$n=\underset{1\le j\le s}{\max }\underset{1\le l\le m}{\max }{K}_j(l),$$

$$P^{\perp }:=\left\{\mathbf{k}=(k_1,\dots ,k_s)\in {\mathbb{N}}_0^s\mid C_1^{\top }{\nu }_n(k_1)\oplus \cdots \oplus C_s^{\top }{\nu }_n(k_s)=0\in {\mathbb{F}}_b^m\right\},$$

$${\nu }_n(k)=({\kappa }_0,\dots ,{\kappa }_{n-1}{)}^{\top }\in {\mathbb{F}}_b^n$$

$$k={\kappa }_1{b}^{c_1-1}+{\kappa }_2{b}^{c_2-1}+\cdots +{\kappa }_v{b}^{c_v-1}$$

$${\mu }_1(k_1,\dots ,k_s)=\sum _{j=1}^s{\mu }_1(k_j).$$

$${\mu }_1(P^{\perp }):=\underset{\mathbf{k}\in P^{\perp }\setminus \{0\}}{\min }{\mu }_1(\mathbf{k}).$$

$$t:=m-{\mu }_1(P^{\perp })+1,$$

$$E=\prod _{j=1}^s\left[\frac{a_j}{b^{c_j}},\frac{a_j+1}{b^{c_j}}\right)$$

$$D^{\ast }(P)\le B_{s,b,t}^{(1)}\frac{m^{s-1}}{b^m}.$$

$$D^{\ast }(\{{\mathbf{x}}_0,\dots ,{\mathbf{x}}_{N-1}\})\le B_{s,b,t}^{(2)}\frac{(\log N{)}^s}{N},$$

$$\left|\left\{\mathbf{k}\in P^{\perp }\setminus \{0\}\mid {\mu }_1(\mathbf{k})=z\right\}\right|\le \{\begin{array}{ll}0& \text{if z<μ1(P⊥),}\\ b^{z-{\mu }_1(P^{\perp })+1}(z+1{)}^{s-1}& \text{otherwise.}\end{array}$$

$$\sum _{\mathbf{k}\in P^{\perp }\setminus \{0\}}{b}^{-\lambda {\mu }_1(\mathbf{k})}\le 2^{s-1}{b}^{\lambda }\frac{({\mu }_1(P^{\perp }){)}^{s-1}}{b^{\lambda {\mu }_1(P^{\perp })}}\sum _{z=1}^{\infty }{b}^{(1-\lambda )z}{z}^{s-1}.$$

$$A_z=\sum _{\begin{array}{c}{z}_1,\dots ,z_s\in {\mathbb{N}}_0\\ z_1+\cdots +z_s=z\end{array}}\left|\left\{\mathbf{k}\in P^{\perp }\setminus \{0\}\mid {\mu }_1(k_1)=z_1,\dots ,{\mu }_1(k_s)=z_s\right\}\right|.$$

$$\left|\left\{\mathbf{k}\in P^{\perp }\setminus \{0\}\mid {\mu }_1(k_1)=z_1,\dots ,{\mu }_1(k_s)=z_s\right\}\right|\le \{\begin{array}{ll}0& \text{if z<μ1(P⊥),}\\ b^{z-{\mu }_1(P^{\perp })+1}& \text{otherwise.}\end{array}$$

$$A_z$$

$$=b^{z-{\mu }_1(P^{\perp })+1}\prod _{j=1}^{s-1}\frac{z+j}{j}\le b^{z-{\mu }_1(P^{\perp })+1}(z+1{)}^{s-1},$$

$$\sum _{\mathbf{k}\in P^{\perp }\setminus \{0\}}{b}^{-\lambda {\mu }_1(\mathbf{k})}$$

$$=b^{-\lambda ({\mu }_1(P^{\perp })-1)}\sum _{z=1}^{\infty }{b}^{(1-\lambda )z}(z+{\mu }_1(P^{\perp }){)}^{s-1}$$

$$\le 2^{s-1}({\mu }_1(P^{\perp }){)}^{s-1}{b}^{-\lambda ({\mu }_1(P^{\perp })-1)}\sum _{z=1}^{\infty }{b}^{(1-\lambda )z}{z}^{s-1},$$

$$\frac{x^{e_j-z-1}}{(p_j(x){)}^i}=\sum _{l=1}^{\infty }\frac{a^{(j)}(i,z,l)}{x^l}\in {\mathbb{F}}_b((x^{-1}))$$

$$c_{k,l}^{(j)}=a^{(j)}\left(\lfloor \frac{k-1}{e_j}\rfloor +1,(k-1)\mathrm{mod}{e}_j,l\right).$$

$$\frac{x^{e_j-1}}{p_j(x)},\dots ,\frac{1}{p_j(x)},\frac{x^{e_j-1}}{(p_j(x){)}^2},\dots ,\frac{1}{(p_j(x){)}^2},\dots .$$

$$t\le \sum _{j=1}^s(e_j-1).$$

Full text

Keywords: Higher order quasi-Monte Carlo, digital nets and sequences, Walsh analysis, discrepancy, numerical integration

Mathematics Subject Classification: 11K38, 41A55, 42C10, 65C05 (primary), 65D30, 65D32

1 Introduction

For an integrable function $f\colon[0,1]^{s}\to\mathbb{R}$, we denote the integral of $f$ by

[TABLE]

Monte Carlo/Quasi-Monte Carlo (QMC) methods are a class of numerical algorithms for approximating $I(f)$ based on pointwise function evaluations. Let $P\subset[0,1]^{s}$ be a finite multiset, that is, if an element occurs multiple times, it is counted according to its multiplicity. Then $I(f)$ is approximated by

[TABLE]

where $|P|$ denotes the cardinality of $P$. It is obvious that this algorithm is exact for any choice of $P$ if $f$ is a constant function, but except for such a trivial case, a careful design of $P$ and the accompanying theoretical analysis are required to show that the algorithm works well for various functions $f$.

One fairly easy but sensible approach is to choose each point $\boldsymbol{x}$ independently and uniformly from $[0,1]^{s}$. This is widely known under the name of Monte Carlo methods [39]. Many fundamental results in probability theory, including the law of large numbers and the central limit theorem, apply to this approach. Looking at $I(f;P)$ as a random variable (with $P$ being the underlying stochastic variable), we have

[TABLE]

for any function $f\in L_{2}([0,1]^{s})$, where $\mathbb{V}[f]$ on the right-hand side of the second equality denotes the variance of $f$. This means, Monte Carlo methods work for any square-integrable functions, but the approximation error converges only probabilistically at the notorious “one over square root of $N$” rate. Thus we have a trade-off between versatility and efficiency.

In some applications where the Monte Carlo convergence is considered too slow, one needs to improve efficiency while discarding versatility of Monte Carlo methods to some extent. QMC methods are one of the standard choices for this purpose. The classical but still central result in QMC methods is the celebrated Koksma-Hlawka inequality:

[TABLE]

where $V_{\mathrm{HK}}(f)$ denotes the total variation of $f$ in the sense of Hardy and Krause, and $D^{*}(P)$ denotes the star-discrepancy of $P$ (we shall give a precise definition of $D^{*}(P)$ later in Section 5). Although a class of functions we can deal with is restricted to functions with bounded total variations (i.e., we discard versatility to some extent), through a clever design of $P$ such that $D^{*}(P)$ is of order better than $|P|^{-1/2}$, the convergence rate can be improved (i.e., we improve efficiency). In fact, there are many explicit constructions of so-called digital $(t,m,s)$-nets and digital $(t,s)$-sequences achieving the star-discrepancy of order $(\log N)^{s-1}/N$ and $(\log N)^{s}/N$, respectively,111To be precise, for an infinite sequence of points $\mathcal{S}$, this means that there exists a constant $C_{s}>0$ depending only on $s$ such that the star-discrepancy of the first $N$ elements of $\mathcal{S}$ is bounded by $C_{s}(\log N)^{s}/N$ uniformly for all $N$. see [49, 24]. Hence, it follows from the Koksma-Hlawka inequality that the integration error decays faster than the “one over square root of $N$” rate.

One natural question in this line is then “Can we improve efficiency further while sacrificing versatility to more extent?” Higher order quasi-Monte Carlo (HoQMC) methods due to Dick [7, 8] provide an affirmative solution to this question. Now let us focus on functions $f$ having square-integrable partial mixed derivatives up to order $\alpha>1$ in each variable, which obviously means that we discard versatility to more extent than standard QMC methods. In return for this drawback, however, the order of convergence of the integration error can be improved to $(\log N)^{c(\alpha,s)}/N^{\alpha}$ with some exponent $c(\alpha,s)>0$ by employing so-called higher order digital nets and sequences as quadrature nodes $P$.222For higher order digital sequences, this order of convergence does not hold uniformly for all $N$, but holds for a geometric spacing of $N$. It is known that this cannot be improved [53]. Hence, when the considered integrand is smooth enough, HoQMC methods can be much more efficient than standard QMC methods, not to mention Monte Carlo methods. Of course, one may doubt if there is any chance of encountering with such smooth functions in practice. Fortunately, there have been several successful and promising applications of HoQMC methods reported already in the literature. These include [18, 17, 20, 30, 31] on applications to partial differential equations with random coefficients (see also the review article [41]), and [11, 32, 12] on applications to Bayesian estimation/inversion problems.

Recently there have been significant theoretical progresses on HoQMC methods. The first major step was made by Dick and Pillichshammer [25]. They proved that order 5 digital sequences achieve the best possible order of $L_{2}$-discrepancy, which is $(\log N)^{s/2}/N$, uniformly for all $N$, and moreover, they proved that order 3 digital nets of $N$ points achieve the best possible order of $L_{2}$-discrepancy, which is $(\log N)^{(s-1)/2}/N$ (here again, we shall give a precise definition of $L_{2}$-discrepancy later in Section 5). Prior to their work, there had been only one explicit construction of finite point sets (for arbitrarily fixed dimension $s$) with the best possible order of $L_{2}$-discrepancy due to Chen and Skriganov [5, 55]. Therefore, higher order digital nets (resp. sequences) are now recognized as the second (resp. first) explicit construction of optimal order $L_{2}$-discrepancy point sets (resp. sequences). More recently, several refined analyses for generalizing or extending the work of Dick and Pillichshammer have been conducted [10, 42, 15, 16, 4].

Another major step was made in a series of papers [36, 37, 38], where the authors refined the integration error analysis for smooth functions due to Dick [7, 8] and proved that order $(2\alpha+1)$ digital nets and sequences achieve the best possible order of the worst-case error for a reproducing kernel Hilbert space with dominating mixed smoothness $\alpha$, which is $(\log N)^{(s-1)/2}/N^{\alpha}$. Note that the original work by Dick [8] proves the worst-case error of order $(\log N)^{s\alpha}/N^{\alpha}$ for order $\alpha$ digital nets and sequences, see also [3]. Other than higher order digital nets and sequences, only the Frolov lattice rule in conjunction with periodization of integrands has been proven to achieve the same, best possible order of the worst-case error so far [29, 64, 45].

There is a common source for obtaining the result of [25] and that of [37, 38], which is the Walsh analysis. To introduce the concept of HoQMC methods originally, Dick managed to prove the decay of the Walsh coefficients of smooth functions [7, 8]. In fact, his digit interlacing construction of higher order digital nets and sequences, which shall be described in Section 3.2, is carefully designed to exploit the decay of the Walsh coefficients. In order to improve his seminal results, it may be sensible to attempt to exploit some further aspect of the Walsh coefficients. Both the result of [25] and that of [37, 38] rely not only on the decay but also on the sparsity of the Walsh coefficients.

In this article we mainly focus on the papers [25, 38] and provide a unified picture on how the Walsh analysis enables recent developments of HoQMC methods both for discrepancy and numerical integration. The rest of this article is organized as follows. In Section 2, we explain about standard QMC methods based on digital nets and sequences in the sense of Niederreiter [49]. Although integer lattices are another important class of QMC point sets, see for instance [56] and [19, Section 5], we do not cover them in this article. In Section 3, we introduce the definitions of higher order digital nets and sequences, and provide an explicit construction algorithm due to Dick [8]. In Section 4, we introduce the definition of the Walsh functions and give some key connection to digital nets. Thereafter, recent advances in HoQMC methods for discrepancy are described in Section 5, while those for numerical integration are in Section 6. We shall highlight an analogy between the approach by [25] and that of [38], where exploiting both the decay and the sparsity of the Walsh coefficients plays a crucial role. We conclude the article with some future research directions.

Notation. Throughout this article, we shall use the following notation. Let $\mathbb{N}$ be the set of positive integers and we write $\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$. For a prime $b$, let $\mathbb{F}_{b}$ be the finite field with $b$ elements and we identify $\mathbb{F}_{b}$ with the set of integers $\{0,1,\ldots,b-1\}$ equipped with addition and multiplication modulo $b$. For $x\in[0,1]$, its $b$-adic expansion $x=\sum_{i=1}^{\infty}\xi_{i}b^{-i}$ with $\xi_{i}\in\mathbb{F}_{b}$ is understood to be unique in the sense that infinitely many of the $\xi_{i}$’s are different from $b-1$ if $x\neq 1$ and that all of the $\xi_{i}$’s are equal to $b-1$. Note that for $k=1\in\mathbb{N}$ we use the $b$-adic expansion $1b^{0}$, whereas for $x=1\in[0,1]$ we use $(b-1)b^{-1}+(b-1)b^{-2}+\cdots$. It will be clear from the context which expansion we use. The operator $\oplus$ denotes the digitwise addition modulo $b$, that is, for $x,y\in[0,1]$ with $b$-adic expansions given by $x=\sum_{i=1}^{\infty}\xi_{i}b^{-i}$ and $y=\sum_{i=1}^{\infty}\eta_{i}b^{-i}$, respectively, $\oplus$ is defined as

[TABLE]

Similarly we use $\oplus$ for digitwise addition for non-negative integers based on the $b$-adic expansions. In case of vectors in $[0,1]^{s}$ or $\mathbb{N}_{0}^{s}$, the operator $\oplus$ is applied componentwise.

2 Standard quasi-Monte Carlo

2.1 Digital nets and sequences

We start with a general construction scheme for a class of QMC point sets called digital nets due to Niederreiter [49]. Note that both of digital $(t,m,s)$-nets and higher order digital nets can be regarded as special subclasses of them.

Definition 1** (Digital nets).**

Let $m,n\in\mathbb{N}$ and let $C_{1},\ldots,C_{s}$ be $n\times m$ matrices over $\mathbb{F}_{b}$. For an integer $0\leq h<b^{m}$ with $b$-adic expansion $h=\eta_{0}+\eta_{1}b+\cdots+\eta_{m-1}b^{m-1}$, define the point $\boldsymbol{x}_{h}=(x_{h,1},\ldots,x_{h,s})\in[0,1]^{s}$ by

[TABLE]

where

[TABLE]

The set $P=\{\boldsymbol{x}_{h}\mid 0\leq h<b^{m}\}\subset[0,1]^{s}$ is called a digital net over $\mathbb{F}_{b}$ (with generating matrices $C_{1},\ldots,C_{s}$).

It is obvious from the definition that the parameter $m$ determines the total number of points, which is $b^{m}$, while the parameter $n$ determines the precision of points. We can extend this definition to construct infinite sequences of points called digital sequences. Again, both of digital $(t,s)$-sequences and higher order digital sequences can be regarded as special subclasses of them.

Definition 2** (Digital sequences).**

Let $C_{1},\ldots,C_{s}$ be $\mathbb{N}\times\mathbb{N}$ matrices over $\mathbb{F}_{b}$. For an integer $h\in\mathbb{N}_{0}$ with $b$-adic expansion $h=\eta_{0}+\eta_{1}b+\cdots$, where all but a finite number of $\eta_{i}$ are 0, define the point $\boldsymbol{x}_{h}=(x_{h,1},\ldots,x_{h,s})\in[0,1]^{s}$ by

[TABLE]

where

[TABLE]

The sequence of points $\mathcal{S}=\{\boldsymbol{x}_{h}\mid h\in\mathbb{N}_{0}\}\subset[0,1]^{s}$ is called a digital sequence over $\mathbb{F}_{b}$ (with generating matrices $C_{1},\ldots,C_{s}$).

Remark 1**.**

Assume that for each $1\leq j\leq s$ there exists a function $K_{j}\colon\mathbb{N}\to\mathbb{N}$ such that $C_{j}=(c_{k,l}^{(j)})_{k,l\in\mathbb{N}}$ satisfies $c_{k,l}^{(j)}=0$ whenever $k>K_{j}(l)$. Then for any $m\in\mathbb{N}$ the first $b^{m}$ elements of a digital sequence over $\mathbb{F}_{b}$ with generating matrices $C_{1},\ldots,C_{s}$ can be identified with a digital net over $\mathbb{F}_{b}$ with generating matrices $C_{1}^{[n\times m]},\ldots,C_{s}^{[n\times m]}$ with the precision parameter

[TABLE]

where $C_{j}^{[n\times m]}$ denotes the upper-left $n\times m$ submatrix of $C_{j}$.

2.2 Quality measure

In order to generate point sets or sequences from the above construction scheme such that the star-discrepancy is small, we need to design generating matrices $C_{1},\ldots,C_{s}$ properly. In this subsection, we introduce the widely-used quality measure called $t$-value, which is based on the Niederreiter-Rosenbloom-Tsfasman (NRT) weight function [46, 54].

Definition 3** (Dual nets).**

Let $m,n\in\mathbb{N}$ and let $P$ be a digital net over $\mathbb{F}_{b}$ with generating matrices $C_{1},\ldots,C_{s}\in\mathbb{F}_{b}^{n\times m}$. The dual net of $P$, denoted by $P^{\perp}$, is defined by

[TABLE]

where

[TABLE]

for $k\in\mathbb{N}_{0}$ with $b$-adic expansion $k=\kappa_{0}+\kappa_{1}b+\cdots$, where all but a finite number of $\kappa_{i}$ are 0.

Remark 2**.**

Let $\mathcal{S}$ be a digital sequence over $\mathbb{F}_{b}$ for which there exist functions $K_{j}\colon\mathbb{N}\to\mathbb{N}$ such that $C_{j}=(c_{k,l}^{(j)})_{k,l\in\mathbb{N}}$ satisfies $c_{k,l}^{(j)}=0$ whenever $k>K_{j}(l)$. Then the dual net can be defined for the first $b^{m}$ elements of $\mathcal{S}$ for any $m\in\mathbb{N}$, since they can be identified with a digital net as discussed in Remark 1.

Definition 4** (NRT weight function).**

For $k\in\mathbb{N}$, we denote the $b$-adic expansion of $k$ by

[TABLE]

with $\kappa_{1},\ldots,\kappa_{v}\in\mathbb{F}_{b}\setminus\{0\}$ and $c_{1}>\cdots>c_{v}>0$. Then the NRT weight function $\mu_{1}\colon\mathbb{N}_{0}\to\mathbb{N}_{0}$ is defined by $\mu_{1}(0)=0$ and $\mu_{1}(k)=c_{1}$. In case of vectors in $\mathbb{N}_{0}^{s}$, we define

[TABLE]

We are ready to introduce the definition of $t$-value.

Definition 5** ($t$-value).**

Let $m,n\in\mathbb{N}$ and let $P$ be a digital net over $\mathbb{F}_{b}$ with generating matrices $C_{1},\ldots,C_{s}\in\mathbb{F}_{b}^{n\times m}$. We write

[TABLE]

Then the $t$-value of $P$ is defined by

[TABLE]

and $P$ is called a digital $(t,m,s)$-net over $\mathbb{F}_{b}$.

This definition of $t$-value is based on the concept of duality theory of digital nets as originally studied in [50]. There is another but equivalent definition of $t$-value: let $\rho$ be the largest integer such that, for any choice $d_{1},\ldots,d_{s}\in\mathbb{N}_{0}$ with $d_{1}+\cdots+d_{s}=\rho$,

the first $d_{1}$ row vectors of $C_{1}$

the first $d_{2}$ row vectors of $C_{2}$

$\vdots$

the first $d_{s}$ row vectors of $C_{s}$

are linearly independent over $\mathbb{F}_{b}$. Then the $t$-value can be also defined by $m-\rho$.

Because of the linear independence of the row vectors of generating matrices, any digital $(t,m,s)$-net over $\mathbb{F}_{b}$ has the following equi-distribution property: every $b$-adic elementary box of the form

[TABLE]

with $c_{1},\ldots,c_{s}\geq 0$, $c_{1}+\cdots+c_{s}=m-t$ and $0\leq a_{j}<b^{c_{j}}$ for all $j$, whose volume is $b^{t-m}$, contains $b^{t}$ points exactly. Hence, as the $t$-value is smaller, digital nets are more equi-distributed over $[0,1]^{s}$. This is why the $t$-value works as a quality measure of digital nets.

For digital sequences, the $t$-value is defined as follows.

Definition 6**.**

Let $\mathcal{S}$ be a digital sequence over $\mathbb{F}_{b}$. $\mathcal{S}$ is called a digital $(t,s)$-sequence over $\mathbb{F}_{b}$ if there exists $t\in\mathbb{N}_{0}$ such that the first $b^{m}$ points of $\mathcal{S}$ is a digital $(t,m,s)$-net for any $m\geq t$.

The following result states that the star-discrepancy of digital $(t,m,s)$-nets and digital $(t,s)$-sequences are of order $(\log N)^{s-1}/N$ and $(\log N)^{s}/N$, respectively, as mentioned in the first section. We refer to [49, Theorems 4.10 and 4.17] for the proof.

Theorem 1**.**

The following holds true:

Let $P$ be a digital $(t,m,s)$-net over $\mathbb{F}_{b}$. There exists a constant $B^{(1)}_{s,b,t}$ such that the star-discrepancy of $P$ is bounded by

[TABLE] 2. 2.

Let $\mathcal{S}=\{\boldsymbol{x}_{h}\mid h\in\mathbb{N}_{0}\}$ be a digital $(t,s)$-sequence over $\mathbb{F}_{b}$. There exists a constant $B^{(2)}_{s,b,t}$ such that the star-discrepancy of the first $N$ points of $\mathcal{S}$ is bounded by

[TABLE]

for any $N\geq 2$.

We end this subsection by providing one useful result in analyzing the integration error of QMC rules using digital $(t,m,s)$-nets.

Lemma 1**.**

Let $P$ be a digital $(t,m,s)$-net over $\mathbb{F}_{b}$. The following holds true:

For $z\in\mathbb{N}_{0}$,

[TABLE] 2. 2.

For any real $\lambda>1$,

[TABLE]

Proof.

In this proof we put $A_{z}=\left|\left\{\boldsymbol{k}\in P^{\perp}\setminus\{\boldsymbol{0}\}\mid\mu_{1}(\boldsymbol{k})=z\right\}\right|$. It holds that

[TABLE]

Following [55, Lemma 2.2], the summand is bounded above by

[TABLE]

Thus we have $A_{z}=0$ if $z<\mu_{1}(P^{\perp})$, since each summand is 0. For $z\geq\mu_{1}(P^{\perp})$, this bound gives

[TABLE]

which proves the first assertion of the lemma.

Using the result of the first assertion and then applying the change of variables $z\mapsto z+\mu_{1}(P)-1$, we have

[TABLE]

where the last sum over $z$ is finite since $\lambda>1$. Hence we complete the proof. ∎

2.3 Explicit constructions

Theorem 1 together with the Koksma-Hlawka inequality (1) gives a motivation to construct digital $(t,m,s)$-nets or digital $(t,s)$-sequences with small $t$-value. In fact, many explicit constructions of digital $(t,s)$-sequences with small $t$-value are already known. Examples are given by Sobol’ [58], Faure [26], Niederreiter [47], Tezuka [63], Niederreiter and Xing [51] as well as many others. Here we give one example from [47, 63].

Let $p_{1},p_{2},\ldots\in\mathbb{F}_{b}[x]$ be a sequence of distinct monic irreducible polynomials over $\mathbb{F}_{b}$ with $\deg(p_{1})\leq\deg(p_{2})\leq\cdots$. For each $j\in\mathbb{N}$, let $e_{j}=\deg(p_{j})$ and consider the following Laurent series expansion

[TABLE]

for integers $i\geq 1$ and $0\leq z<e_{j}$. Define the matrix $C_{j}=(c_{k,l}^{(j)})_{k,l\in\mathbb{N}}$ by

[TABLE]

Here we see that the rows of $C_{j}$ (from upper to lower) correspond to the Laurent series expansions of

[TABLE]

Hence, we have $K_{j}(l)=l$ for any $j,l$ in the light of Remark 1, that is, $c_{k,l}^{(j)}=0$ whenever $k>l$, meaning that $C_{j}$ is an upper triangular matrix.

It is straightforward from the definition that this explicit construction of digital sequences is extensible in dimension. The first $s$ matrices $C_{1},\ldots,C_{s}$ generate a digital $(t,s)$-sequence over $\mathbb{F}_{b}$ with

[TABLE]

We refer to [24, Theorem 8.2] for the proof of this fact.

Remark 3**.**

Some comments are in order:

Replacing the numerator $x^{e_{j}-z-1}$ of (2) by $x^{z}$, this construction algorithm is the same as the one originally introduced in **[47]**, which is nowadays known as Niederreiter sequences. For the original Niederreiter sequence, the $t$-value is strictly equal to $\sum_{j=1}^{s}(e_{j}-1)$ **[21]**. 2. 2.

A generalization of Niederreiter sequences by Tezuka **[63]** is to use the set of linearly independent polynomials $\{y_{j,i,z}(x)\mid 0\leq z<e_{j}\}$ over $\mathbb{F}_{b}$ for the numerator of (2) rather than the simplest set $\{x^{z}\mid 0\leq z<e_{j}\}$. 3. 3.

The Sobol’ sequences due to **[58]** are a subclass of the generalized Niederreiter sequences over $\mathbb{F}_{2}$, where the primitive polynomials are used for $p_{1},p_{2},\ldots$. Recently, Faure and Lemieux **[27]** gave some precise connections between the Sobol’ sequences and the generalized Niederreiter sequences.

2.4 Polynomial lattice point sets

We end this section by introducing another important class of digital nets called polynomial lattice point sets introduced by Niederreiter [48].

Definition 7**.**

Let $m\in\mathbb{N}$, and let $p\in\mathbb{F}_{b}[x]$ and $\boldsymbol{q}=(q_{1},\ldots,q_{s})\in(\mathbb{F}_{b}[x])^{s}$ such that $\deg(p)=m$ and $\deg(q_{j})<m$. For $1\leq j\leq s$, consider the Laurent series expansion

[TABLE]

and define the Hankel matrix $C_{j}=(c_{k,l}^{(j)})_{1\leq k,l\leq m}\in\mathbb{F}_{b}^{m\times m}$ by

[TABLE]

Then a digital net over $\mathbb{F}_{b}$ with these generating matrices $C_{1},\ldots,C_{s}$ is called a polynomial lattice point set (with modulus $p$ and generating vector $\boldsymbol{q}$).

Indeed, polynomial lattice point sets can be constructed without using generating matrices explicitly. Define the map $v_{m}\colon\mathbb{F}_{b}((x^{-1}))\to[0,1]$ by

[TABLE]

We identify $h\in\mathbb{N}_{0}$, whose finite $b$-adic expansion is given by $h=\eta_{0}+\eta_{1}b+\cdots$, with the polynomial over $\mathbb{F}_{b}$ given by $h(x)=\eta_{0}+\eta_{1}x+\cdots$. Put

[TABLE]

Then a point set $\{\boldsymbol{x}_{h}\mid 0\leq h<b^{m}\}$ is nothing but the polynomial lattice point set as defined above.

The modulus $p$ is often chosen to be either the monomial $p(x)=x^{m}$ or irreducible. The difficulty is in how to choose the generating vector $\boldsymbol{q}$. In particular, for $s\geq 3$, no explicit way for this choice has been known yet. Currently one of the most standard approaches is to recursively choose one component $q_{j}$ from the set $\{q\in\mathbb{F}_{b}[x]\mid\deg(q)<m\}$ which minimizes a considered criterion while the earlier ones $q_{1},\ldots,q_{j-1}$ kept unchanged. This greedy algorithm is known as component-by-component construction [57]. Another well-known approach is to restrict ourselves to the form

[TABLE]

for $q\in\mathbb{F}_{b}[x]$ with $\deg(q)<m$, and then to choose one optimal $q$ with respect to a considered criterion. This algorithm is known as Korobov construction. Compared to Korobov construction, component-by-component construction has the advantages that it is extensible in dimension and that the fast algorithm using the fast Fourier transform has been well-established [52]. However, neither of both is extensible in the number of points.

3 Higher order quasi-Monte Carlo

3.1 Quality measure

To introduce the definitions of higher order digital nets and sequences, we start with generalizing the NRT weight function.

Definition 8** (Dick weight function).**

Let $\alpha\in\mathbb{N}$. For $k\in\mathbb{N}$ we denote the $b$-adic expansion of $k$ by

[TABLE]

with $\kappa_{1},\ldots,\kappa_{v}\in\mathbb{F}_{b}\setminus\{0\}$ and $c_{1}>\cdots>c_{v}>0$. Then the Dick weight function $\mu_{\alpha}\colon\mathbb{N}_{0}\to\mathbb{N}_{0}$ is defined by $\mu_{\alpha}(0)=0$ and

[TABLE]

In case of vectors in $\mathbb{N}_{0}^{s}$ we define

[TABLE]

It is obvious that the Dick weight function coincides with the NRT weight function when $\alpha=1$. As a natural generalization of digital $(t,m,s)$-nets and $(t,s)$-sequences based on the NRT weight function, higher order digital nets and sequences due to Dick [7, 8] are defined by using the Dick weight function as follows:

Definition 9** (Higher order digital nets).**

Let $\alpha\in\mathbb{N}$. Let $m,n\in\mathbb{N}$ and let $P$ be a digital net over $\mathbb{F}_{b}$ with generating matrices $C_{1},\ldots,C_{s}\in\mathbb{F}_{b}^{n\times m}$. We write

[TABLE]

Then $P$ is called an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$ with

[TABLE]

Definition 10** (Higher order digital sequences).**

Let $\mathcal{S}$ be a digital sequence over $\mathbb{F}_{b}$. $\mathcal{S}$ is called an order $\alpha$ digital $(t_{\alpha},s)$-sequence over $\mathbb{F}_{b}$ if there exists $t_{\alpha}\in\mathbb{N}_{0}$ such that the first $b^{m}$ points of $\mathcal{S}$ is an order $\alpha$ digital $(t_{\alpha},m,s)$-net for any $m\geq t_{\alpha}/\alpha$.

Obviously, for fixed $\alpha$, the $t_{\alpha}$-value works as a quality measure of order $\alpha$ digital nets and sequences.

Our definition of higher order digital nets is again based on the concept of duality theory of digital nets, and looks different from the original definition by Dick which has an additional parameter $\beta$. When $n\geq\alpha m$, however, by setting $\beta=1$, our definition becomes equivalent to his original definition based on the linear independence of row vectors of generating matrices described below: let $\rho$ be the largest integer such that, for any choice $1\leq d_{j,\nu_{j}}<\cdots<d_{j,1}\leq n$, where $0\leq\nu_{j}\leq m$ for all $1\leq j\leq s$ with

[TABLE]

the $d_{1,\nu_{1}},\ldots,d_{1,1}$-th row vectors of $C_{1}$

the $d_{2,\nu_{2}},\ldots,d_{2,1}$-th row vectors of $C_{2}$

$\vdots$

the $d_{s,\nu_{s}},\ldots,d_{s,1}$-th row vectors of $C_{s}$

are linearly independent over $\mathbb{F}_{b}$. Then a digital net with generating matrices $C_{1},\ldots,C_{s}$ is an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$ with $t_{\alpha}=\alpha m-\rho$. This linear independence of the row vectors of generating matrices ensures that a higher order digital net has a similar geometric equi-distribution property to what is described in Subsection 2.2, see [24, Chapter 15.3] for more details.

Here we provide one useful property called propagation rule of higher order digital nets and sequences shown in [7, Theorem 3.3]. We give a different proof which does not rely on the linear independence of generating matrices.

Lemma 2**.**

For $\beta\in\mathbb{N}$, $\beta\geq 2$, let $P$ and $\mathcal{S}$ be an order $\beta$ digital $(t_{\beta},m,s)$-net and an order $\beta$ digital $(t_{\beta},s)$-sequence over $\mathbb{F}_{b}$, respectively. Then, for any $1\leq\alpha<\beta$, $P$ and $\mathcal{S}$ are also an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$ and an order $\alpha$ digital $(t_{\alpha},s)$-sequence over $\mathbb{F}_{b}$, respectively, both with $t_{\alpha}\leq\lceil t_{\beta}\alpha/\beta\rceil$.

Proof.

First we prove that the inequality

[TABLE]

holds for any $\boldsymbol{k}\in\mathbb{N}_{0}^{s}$ and $1\leq\alpha\leq\beta$. Since the weight function for vector $\boldsymbol{k}$ is defined as the sum of the weight function for each coordinate, it suffices to prove the one-dimensional case. Since the result for $k=0$ follows trivially, let us consider $k>0$. Denote the $b$-adic expansion of $k$ by

[TABLE]

with $\kappa_{1},\ldots,\kappa_{v}\in\mathbb{F}_{b}\setminus\{0\}$ and $c_{1}>\cdots>c_{v}>0$, and write $c_{v+1}=c_{v+2}=\cdots=0$. Then we have

[TABLE]

which proves the assertion.

Now let us consider an order $\beta$ digital $(t_{\beta},m,s)$-net over $\mathbb{F}_{b}$. Using the above inequality we obtain

[TABLE]

Thus $P$ is an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$ with

[TABLE]

Given that $t_{\alpha}$ is a non-negative integer and that the fraction $(\beta-\alpha)/\beta$ is less than 1, the $t_{\alpha}$-value can be bounded above by $\lceil t_{\beta}\alpha/\beta\rceil$. The result for an order $\beta$ digital $(t_{\beta},s)$-sequence follows immediately. ∎

Importantly this result implies that any higher order digital nets and sequences are also digital $(t,m,s)$-nets and digital $(t,s)$-sequences, respectively.

We end this subsection by providing two useful results in analyzing the integration error of QMC rules using higher order digital nets. The first lemma is a higher order version of Lemma 1.

Lemma 3**.**

For $\alpha\geq 2$, let $P$ be an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$. The following holds true:

For $z\in\mathbb{N}_{0}$,

[TABLE] 2. 2.

For any real $\lambda>1/\alpha$,

[TABLE]

Proof.

Similarly to the proof of Lemma 1, we put $A_{\alpha,z}=\left|\left\{\boldsymbol{k}\in P^{\perp}\setminus\{\boldsymbol{0}\}\mid\mu_{\alpha}(\boldsymbol{k})=z\right\}\right|$. For $1\leq j\leq s$, denote the $b$-adic expansion of $k_{j}\in\mathbb{N}_{0}$ by

[TABLE]

with $\kappa_{1,j},\ldots,\kappa_{v_{j},j}\in\mathbb{F}_{b}\setminus\{0\}$ and $c_{1,j}>\cdots>c_{v_{j},j}>0$, and write $c_{v_{j}+1,j}=c_{v_{j}+2,j}=\cdots=0$. If $k_{j}=0$, let $c_{1,j}=c_{2,j}=\cdots=0$. Then we have

[TABLE]

It can be inferred from [24, Proof of Lemma 15.20] that the summand is bounded above by

[TABLE]

Thus, we have $A_{\alpha,z}=0$ if $z<\mu_{\alpha}(P^{\perp})$, since each summand is 0. For $z\geq\mu_{\alpha}(P^{\perp})$, this bound gives

[TABLE]

which proves the first assertion of the lemma.

Using the result of the first assertion and then applying the change of variables $z\mapsto z+\mu_{\alpha}(P)$, we have

[TABLE]

where the last sum over $z$ is finite since $\lambda>1/\alpha$. Hence we complete the proof. ∎

Before stating the second useful lemma, we need to recall the notion of “type $(p,q)$” introduced in [38].

Definition 11**.**

For $k,l\in\mathbb{N}_{0}$, we denote the $b$-adic expansions of $k$ and $l$ by

[TABLE]

respectively, where $\kappa_{1},\ldots,\kappa_{v},\lambda_{1},\ldots,\lambda_{w}\in\{1,\ldots,b-1\}$, $c_{1}>c_{2}>\cdots>c_{v}>0$ and $d_{1}>d_{2}>\cdots>d_{w}>0$. For $k=0$ ($l=0$, resp.), we assume that $v=0$ and $\kappa_{0}b^{c_{0}-1}=0$ ($w=0$ and $\lambda_{0}b^{d_{0}-1}=0$, resp.). For $p,q\in\mathbb{N}_{0}$, we write

[TABLE]

where the empty sum equals 0. Then we say that $(k,l)$ is of type $(p,q)$ if $k^{(p)}=l^{(q)}$ and $\kappa_{p}b^{c_{p}-1}\neq\lambda_{q}b^{d_{q}-1}$, where we set $\kappa_{0}b^{c_{0}-1}=\lambda_{0}b^{d_{0}-1}=0$, except the case $k=l$ where we say that $(k,l)$ is of type $(0,0)$.

In what follows, we write $(k,l)\in T_{\geq\alpha}$ if $(k,l)$ is of type $(p,q)$ with $p+q\geq\alpha$. In case of vectors in $\mathbb{N}_{0}^{s}$ we write $(\boldsymbol{k},\boldsymbol{l})\in T_{\geq\alpha}$ if there exists at least one index $1\leq j\leq s$ such that $(k_{j},l_{j})\in T_{\geq\alpha}$.

Now the following result, which can be regarded as a generalization of the result shown in [25, Lemma 3.7], is proven in [38, Lemma 8]. Here we state the result in a slightly more general form.

Lemma 4**.**

For $\alpha\in\mathbb{N}$, let $P$ be an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$. For $z\in\mathbb{N}_{0}$, there exists a constant $B_{\alpha,b,s,t_{\alpha}}>0$ such that the following holds:

[TABLE]

3.2 Digit interlacing construction

Here we give an explicit construction of higher order digital nets and sequences based on the digit interlacing function due to Dick [7, 8]:

Definition 12**.**

Let $\alpha\in\mathbb{N}$ and let $\boldsymbol{x}=(x_{1},\ldots,x_{\alpha})\in[0,1]^{\alpha}$. For $1\leq j\leq\alpha$, we denote the $b$-adic expansion of $x_{j}$ by $x_{j}=\sum_{i=1}^{\infty}\xi_{i,j}b^{-i}$. The digit interlacing function (of factor $\alpha)$ $\mathcal{D}_{\alpha}:[0,1]^{\alpha}\to[0,1]$ is defined by

[TABLE]

In case of vectors in $[0,1]^{\alpha s}$, we apply $\mathcal{D}_{\alpha}$ to every non-overlapping consecutive $\alpha$ components, i.e.,

[TABLE]

Lemma 5**.**

The following holds true:

Let $P$ be a digital $(t,m,\alpha s)$-net over $\mathbb{F}_{b}$. The set $\mathcal{D}_{\alpha}(P)=\{\mathcal{D}_{\alpha}(\boldsymbol{x})\mid\boldsymbol{x}\in P\}$ is an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$ with

[TABLE] 2. 2.

Let $\mathcal{S}$ be a digital $(t,\alpha s)$-sequence over $\mathbb{F}_{b}$. The sequence $\mathcal{D}_{\alpha}(\mathcal{S})=\{\mathcal{D}_{\alpha}(\boldsymbol{x})\mid\boldsymbol{x}\in\mathcal{S}\}$ is an order $\alpha$ digital $(t_{\alpha},s)$-sequence over $\mathbb{F}_{b}$ with

[TABLE]

Since there are many explicit constructions of digital $(t,m,s)$-nets and $(t,s)$-sequences with small $t$-value for an arbitrarily dimension $s$, as described in Subsection 2.3, the above lemma from [8, Theorems 4.11 and 4.12] directly implies that higher order digital nets and sequences can be explicitly constructed.

Remark 4**.**

Let $P$ be a digital $(t,m,\alpha s)$-net over $\mathbb{F}_{b}$ with generating matrices $C_{1},\ldots,C_{\alpha s}\in\mathbb{F}_{b}^{m\times m}$. Let $\boldsymbol{c}_{i}^{(j)}$ denote the $i$-th row vector of $C_{j}$. For each $1\leq j\leq s$, construct the matrix $D_{j}\in\mathbb{F}_{b}^{\alpha m\times m}$, whose $i$-th row vector is denoted by $\boldsymbol{d}_{i}^{(j)}$, from the matrices $C_{\alpha(j-1)+1},\ldots,C_{\alpha j}$ as

[TABLE]

Then the set $\mathcal{D}_{\alpha}(P)$ is a digital net over $\mathbb{F}_{b}$ with generating matrices $D_{1},\ldots,D_{s}$.

Similarly, the sequence $\mathcal{D}_{\alpha}(\mathcal{S})$ can be identified with a digital sequence over $\mathbb{F}_{b}$ with generating matrices $D_{1},\ldots,D_{s}\in\mathbb{F}_{b}^{\mathbb{N}\times\mathbb{N}}$ which are constructed from the generating matrices $C_{1},\ldots,C_{\alpha s}\in\mathbb{F}_{b}^{\mathbb{N}\times\mathbb{N}}$ of $\mathcal{S}$.

To construct an interlaced finite point set $\mathcal{D}_{\alpha}(P)$, one can use polynomial lattice point sets in dimension $\alpha s$ instead of digital $(t,m,\alpha s)$-nets. The resulting point set $\mathcal{D}_{\alpha}(P)$ is called an interlaced polynomial lattice point set, and has been often used in applications of HoQMC methods, see [18, 17, 20, 11, 30, 31, 32, 12]. Here we need to find good generating vectors $\boldsymbol{q}=(q_{1},\ldots,q_{\alpha s})$, but the digit interlacing composition makes it nontrivial whether each component $q_{j}$ can be searched for one-by-one or consecutive $\alpha$ components $q_{\alpha(j-1)+1},\ldots,q_{\alpha j}$ should be searched for simultaneously. The papers [35, 33] originally gave a justification for employing the former approach, i.e., component-by-component construction.

4 Walsh functions

4.1 Definitions

The Walsh functions were originally introduced by Walsh [66] and have been studied thereafter, for instance, in [28, 6]. In what follows, let $\omega_{b}$ denote the primitive $b$-th root of unity $\exp(2\pi\sqrt{-1}/b)$. The one-dimensional Walsh functions are defined as follows.

Definition 13**.**

Let $k\in\mathbb{N}_{0}$ with $b$-adic expansion $k=\sum_{i=0}^{\infty}\kappa_{i}b^{i}$, where all but a finite number of $\kappa_{i}$ are 0. The $k$-th $b$-adic Walsh function ${}_{b}\mathrm{wal}_{k}\colon[0,1]\to\{1,\omega_{b},\ldots,\omega_{b}^{b-1}\}$ is defined by

[TABLE]

where the unique $b$-adic expansion of $x\in[0,1]$ is denoted by $x=\sum_{i=1}^{\infty}\xi_{i}b^{-i}$.

It is clear from the definition that every Walsh function is piecewise constant since it depends only on some finite number of the digits $\xi_{i}$. Multi-dimensional Walsh functions are given by generalizing the one-dimensional Walsh functions:

Definition 14**.**

Let $s\geq 1$. For $\boldsymbol{k}=(k_{1},\ldots,k_{s})\in\mathbb{N}_{0}^{s}$, the $\boldsymbol{k}$-th $b$-adic Walsh function ${}_{b}\mathrm{wal}_{\boldsymbol{k}}\colon[0,1]^{s}\to\{1,\omega_{b},\ldots,\omega_{b}^{b-1}\}$ is defined by

[TABLE]

Several important properties of the Walsh functions are listed below, see [24, Appendix A.2] for the proof.

Lemma 6**.**

The following holds true:

For $\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}$ and $\boldsymbol{x},\boldsymbol{y}\in[0,1]^{s}$,

[TABLE] 2. 2.

For $\boldsymbol{k}\in\mathbb{N}_{0}^{s}$,

[TABLE] 3. 3.

For $\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}$,

[TABLE] 4. 4.

For any $s\in\mathbb{N}$, the Walsh system $\{{}_{b}\mathrm{wal}_{\boldsymbol{k}}\mid\boldsymbol{k}\in\mathbb{N}_{0}^{s}\}$ is a complete orthonormal system in $L_{2}([0,1]^{s})$

It follows from the fourth assertion of Lemma 6 that we can define the Walsh series of $f\in L_{2}([0,1]^{s})$:

[TABLE]

where $\hat{f}(\boldsymbol{k})$ denotes the $\boldsymbol{k}$-th Walsh coefficient of $f$ defined by

[TABLE]

For any continuous function $f$ which satisfies $\sum_{\boldsymbol{k}\in\mathbb{N}_{0}^{s}}|\hat{f}(\boldsymbol{k})|<\infty$, the above Walsh series of $f$ equals pointwise to $f$ itself, see [24, Theorem A.20].

4.2 Connection to digital nets

It follows from the first assertion of Lemma 6 that the Walsh functions hold the following important character property. The proof can be found, for instance, in [24, Lemma 4.75].

Lemma 7**.**

Let $P\subset[0,1]^{s}$ be a digital net over $\mathbb{F}_{b}$. For $\boldsymbol{k}\in\mathbb{N}_{0}^{s}$ we have

[TABLE]

In what follows, by using this lemma, we show how the Walsh functions play a crucial role in analyzing the QMC integration error.

As preparation, let us consider a reproducing kernel Hilbert space $H$ equipped with reproducing kernel $K\colon[0,1]^{s}\times[0,1]^{s}\to\mathbb{R}$ and inner product $\langle\cdot,\cdot\rangle_{K}$. The norm of $f\in H$ is simply given by $\|f\|_{K}=\sqrt{\langle f,f\rangle_{K}}$. The worst-case error of QMC integration using a point set $P$ is defined by

[TABLE]

while the initial error is defined as reference by

[TABLE]

Both of the initial error and the worst-case error have explicit formulas relying only on $K$ and $P$ as follows, see [24, Chapter 2.3.3] for the proof.

Proposition 1**.**

For a reproducing kernel Hilbert space $H$ whose reproducing kernel satisfies $\int_{[0,1]^{s}}\sqrt{K(\boldsymbol{x},\boldsymbol{x})}\,\mathrm{d}\boldsymbol{x}<\infty$, the squared initial error is given by

[TABLE]

The squared worst-case error of QMC integration using a point set $P$ is given by

[TABLE]

Now let us consider the Walsh series of a reproducing kernel $K$:

[TABLE]

where the $(\boldsymbol{k},\boldsymbol{l})$-th Walsh coefficient $\hat{K}(\boldsymbol{k},\boldsymbol{l})$ is defined by

[TABLE]

Again the pointwise equality holds between the above Walsh series and $K$ itself if $K$ is continuous and $\sum_{\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}}|\hat{K}(\boldsymbol{k},\boldsymbol{l})|<\infty$.

The following proposition provides a simple expression of the squared worst-case error when the point set is a digitally shifted digital net

[TABLE]

where $P$ is a digital net over $\mathbb{F}_{b}$ and $\boldsymbol{\delta}\in[0,1]^{s}$. As far as the authors know, this result is not available in the literature in this full generality, so that we provide a proof for the sake of completeness.

Proposition 2**.**

Let $P$ be a digital net over $\mathbb{F}_{b}$ and $\boldsymbol{\delta}\in[0,1]^{s}$. For a reproducing kernel Hilbert space $H$ whose reproducing kernel $K$ is continuous and satisfies $\int_{[0,1]^{s}}\sqrt{K(\boldsymbol{x},\boldsymbol{x})}\,\mathrm{d}\boldsymbol{x}<\infty$ and $\sum_{\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}}|\hat{K}(\boldsymbol{k},\boldsymbol{l})|<\infty$, we have

[TABLE]

Proof.

It is trivial from the definition of the Walsh functions that the first term on the right-hand side of (1) is equal to $\hat{K}(\boldsymbol{0},\boldsymbol{0})$. For the second term on the right-hand side of (1), by using the symmetry of $K$, the first and second assertions of Lemma 6 and Lemma 7, we have

[TABLE]

Finally, for the third term on the right-hand side of (1), by using the first assertion of Lemma 6 and Lemma 7, we have

[TABLE]

Altogether we obtain

[TABLE]

which completes the proof. ∎

Using Proposition 2 we obtain the following result.

Corollary 1**.**

Let $P$ be a digital net over $\mathbb{F}_{b}$ and $H$ be a reproducing kernel Hilbert space whose reproducing kernel $K$ is continuous and satisfies $\int_{[0,1]^{s}}\sqrt{K(\boldsymbol{x},\boldsymbol{x})}\,\mathrm{d}\boldsymbol{x}<\infty$ and $\sum_{\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}}|\hat{K}(\boldsymbol{k},\boldsymbol{l})|<\infty$. Then we have

[TABLE]

and

[TABLE]

Proof.

The first assertion follows immediately from Proposition 2 by considering the case $\boldsymbol{\delta}=\boldsymbol{0}$. For the second assertion, it follows from Proposition 2 and the third assertion of Lemma 6 that

[TABLE]

Thus we are done. ∎

We would like to emphasize that the worst-case error of QMC rule using a digital net is given by the double sum of the Walsh coefficients of reproducing kernel, whereas the shift-averaged worst-case error is simply given by the single sum of their diagonal elements. This way applying a random digital shift has an effect on vanishing all non-diagonal terms, which sometimes makes an error analysis much easier as we shall see in the subsequent sections.

5 Discrepancy

5.1 Definitions

As represented by the Koksma-Hlawka inequality (1), the star-discrepancy, or more generally speaking, the $L_{p}$-discrepancy is an extremely important quantitative measure of how uniformly a point set is distributed over $[0,1]^{s}$.

Definition 15**.**

Let $P$ be a point set in $[0,1]^{s}$. For $\boldsymbol{y}=(y_{1},\ldots,y_{s})\in[0,1]^{s}$, we write $[\boldsymbol{0},\boldsymbol{y})=[0,y_{1})\times\cdots\times[0,y_{s})$ and define the so-called local discrepancy function

[TABLE]

where $\boldsymbol{1}_{[\boldsymbol{0},\boldsymbol{y})}$ denotes the indicator function which is equal to 1 if $\boldsymbol{x}\in[\boldsymbol{0},\boldsymbol{y})$, and 0 otherwise. For $1\leq p\leq\infty$, the $L_{p}$-discrepancy of $P$ is defined as the $L_{p}$-norm of $\Delta_{P}(\boldsymbol{y})$, i.e.,

[TABLE]

with the obvious modification for $p=\infty$.

We speak of the star-discrepancy if $p=\infty$ and the different notions such as $D^{*}(P)$ have been often used in the literature.

5.2 $L_{2}$-discrepancy and worst-case error

In this section, we are particularly concerned with the $L_{2}$-discrepancy. Here we follow the exposition of [24, Chapter 2.4] to give a connection between the $L_{2}$-discrepancy and the worst-case error for some reproducing kernel Hilbert space.

Let $H^{\downarrow}$ denote the reproducing kernel Hilbert space whose reproducing kernel is given by

[TABLE]

It is known from [1, Section 8] that the function space $H^{\downarrow}$ is the $s$-fold tensor product of the univariate reproducing kernel Hilbert space with reproducing kernel $K^{\downarrow}(x,y)=\min(1-x,1-y)$. Moreover, $H^{\downarrow}$ contains all absolutely continuous functions which satisfy

[TABLE]

where we write $\boldsymbol{x}_{u}=(x_{j})_{j\in u}$ and $(\boldsymbol{x}_{u},\boldsymbol{1})$ is the vector whose $j$-th component is equal to $x_{j}$ if $j\in u$, and 1 otherwise. Note that the symbol $\downarrow$ is used to represent this “anchor” property of the space. The inner product is given by

[TABLE]

for $f,g\in H^{\downarrow}$. Then the following result is known.

Lemma 8**.**

For any point set $P\subset[0,1]^{s}$, we have

[TABLE]

Therefore, when $P$ is a digital net over $\mathbb{F}_{b}$, combining this identity with Corollary 1 gives an expression of the squared $L_{2}$-discrepancy:

[TABLE]

and also an expression of the shift-averaged squared $L_{2}$-discrepancy:

[TABLE]

In [25, Lemma 2.2], Dick and Pillichshammer conducted an exact evaluation of the Walsh coefficients $\widehat{K^{\downarrow}}(\boldsymbol{k},\boldsymbol{l})$ for all $\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}$ when $b=2$. Here we simplify their result in a way that will be sufficient for readers to understand the proof of optimal order $L_{2}$-discrepancy bounds.

Lemma 9**.**

For $\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}$, we have

$|\widehat{K^{\downarrow}}(\boldsymbol{k},\boldsymbol{l})|\leq 3^{-s}\cdot 2^{-\mu_{1}(\boldsymbol{k})-\mu_{1}(\boldsymbol{l})}$. 2. 2.

$\widehat{K^{\downarrow}}(\boldsymbol{k},\boldsymbol{l})=0$* if $(\boldsymbol{k},\boldsymbol{l})\in T_{\geq 3}$.*

Here we recall that the notion $T_{\geq\alpha}$ was introduced in Subsection 3.1. The first assertion of the lemma shows the decay of the Walsh coefficients, while the second assertion shows the sparsity of the Walsh coefficients.

5.3 Optimal order $L_{2}$-discrepancy bounds

Before stating the result by Dick and Pillichshammer [25] on the $L_{2}$-discrepancy bound for higher order digital nets, we start from another result by the same authors shown in [23] by using Lemma 1.

Theorem 2**.**

Let $P$ be a digital $(t,m,s)$-net over $\mathbb{F}_{2}$. Then there exists a constant $D^{(1)}_{s,t}$ such that the following holds:

[TABLE]

Proof.

Using Lemma 9 and the second assertion of Lemma 1 (with $\lambda=2$) we have

[TABLE]

Since we have $\mu_{1}(P^{\perp})=m-t+1$, the result of the theorem follows. ∎

This theorem directly implies the existence of a digital shift $\boldsymbol{\delta}\in[0,1]^{s}$ such that the digitally shifted digital $(t,m,s)$-net satisfies the best possible order of $L_{2}$-discrepancy:

[TABLE]

However, this is a probabilistic result since we do not know how to find such $\boldsymbol{\delta}$ explicitly.

We note that the optimal exponent $(s-1)$ of the numerator in Theorem 2 comes from the inequality

[TABLE]

given in the first assertion of Lemma 1. As is clear from the proof, since the expression of the shift-averaged squared $L_{2}$-discrepancy has only the diagonal terms of the Walsh coefficients, there is no necessity of exploiting the sparsity of the Walsh coefficients. In order to obtain a deterministic counterpart of Theorem 2, however, it seems insufficient to exploit the decay of the Walsh coefficients only, and one approach is to exploit both the decay and the sparsity of the Walsh coefficients simultaneously. To do this, we rely on Lemma 4. The following theorem is from [25, Theorem 4.1].

Theorem 3**.**

Let $P$ be an order 3 digital $(t_{3},m,s)$-net over $\mathbb{F}_{2}$. Then there exists a constant $D_{s,t_{3}}^{(2)}$ such that the following holds:

[TABLE]

Proof.

Using Lemmas 9 and 4 (with $\alpha=3$) and then applying the change of variables $z\to z+2\mu_{1}(P^{\perp})$, we have

[TABLE]

The last sum over $z$ is trivially finite. It follows from Lemma 2 that an order 3 digital $(t_{3},m,s)$-net over $\mathbb{F}_{2}$ is also a digital $(t,m,s)$-net over $\mathbb{F}_{2}$ with $t\leq\lceil t_{3}/3\rceil$, i.e., $\mu_{1}(P^{\perp})\geq m-\lceil t_{3}/3\rceil+1$. Thus, the result of the theorem follows. ∎

Remark 5**.**

If we do not take the sparsity of the Walsh coefficients into account, we might proceed like

[TABLE]

where we used Lemma 1 in the last inequality. Since the last sum over $z$ obviously diverges, this argument ends up with a trivial upper bound.

We recall that Theorem 3 is for order 3 digital nets over $\mathbb{F}_{2}$, a class of finite point sets. By considering order 5 digital sequences, as the other main result of the paper [25], Dick and Pillichshammer proved the optimal order $L_{2}$-discrepancy bound which holds uniformly for all $N$.

6 Numerical integration

6.1 Sobolev spaces

Let us move onto multivariate numerical integration problem. The function space of our interest in this section is defined as follows. First let us consider the one-dimensional case. For $\alpha\in\mathbb{N}$, the Sobolev space with smoothness $\alpha$ is

[TABLE]

where $f^{(r)}$ denotes the $r$-th derivative of $f$. According to [65, Chapter 10.2], the space $H_{\alpha}$ is a reproducing kernel Hilbert space with the reproducing kernel

[TABLE]

for $x,y\in[0,1]$, where $B_{r}$ denotes the Bernoulli polynomial of degree $r$, and with the inner product

[TABLE]

for $f,g\in H_{\alpha}$.

For the $s$-dimensional case, we consider the $s$-fold tensor product space of the one-dimensional space introduced above. Thus the Sobolev space $H_{\alpha,s}$ which we consider is simply given by $H_{\alpha,s}=\bigotimes_{j=1}^{s}H_{\alpha}$. Again it is known from [1, Section 8] that the reproducing kernel of the space $H_{\alpha,s}$ is the product of the reproducing kernels for the one-dimensional space $H_{\alpha}$. Therefore, $H_{\alpha,s}$ is the reproducing kernel Hilbert space with the reproducing kernel

[TABLE]

and with the inner product

[TABLE]

for $f,g\in H_{\alpha,s}$. In the above, we used the following notation: For $u\subseteq\{1,\ldots,s\}$ and $\boldsymbol{x}\in[0,1]^{s}$, we write $\boldsymbol{x}_{u}=(x_{j})_{j\in u}$. Moreover, for $\boldsymbol{r}_{u}=(r_{j})_{j\in u}\in\{0,\ldots,\alpha-1\}^{|u|}$, $(\boldsymbol{r}_{u},\boldsymbol{\alpha})$ denotes the $s$-dimensional vector whose $j$-th component equals $r_{j}$ if $j\in u$, and $\alpha$ otherwise. Note that an integral and sum over the empty set is defined to be the identity operator. Since the dimension $s$ is fixed, we shall simply write $H_{\alpha}$ and $K_{\alpha}$ instead of $H_{\alpha,s}$ and $K_{\alpha,s}$, respectively.

If $P$ is a digital net over $\mathbb{F}_{b}$, Corollary 1 gives

[TABLE]

and

[TABLE]

Similarly to Lemma 9, the following result is known.

Lemma 10**.**

For $\boldsymbol{k},\boldsymbol{l}\in\mathbb{N}_{0}^{s}$, we have

$|\widehat{K_{\alpha}}(\boldsymbol{k},\boldsymbol{l})|\leq C_{\alpha,b}^{s}b^{-\mu_{\alpha}(\boldsymbol{k})-\mu_{\alpha}(\boldsymbol{l})}$* with*

[TABLE]

where

[TABLE] 2. 2.

$\widehat{K_{\alpha}}(\boldsymbol{k},\boldsymbol{l})=0$* if $(\boldsymbol{k},\boldsymbol{l})\in T_{\geq 2\alpha+1}$.*

The first assertion on the decay of the Walsh coefficients was shown in [3], while the second assertion of the sparsity of the Walsh coefficients was shown in [38]. Regarding the first assertion, we also refer to more recent works [62, 67] which introduce different approaches from the one by Dick [7, 8, 9] for evaluating the Walsh coefficients. In many cases, one may obtain smaller constants $C_{\alpha,b}$.

Remark 6**.**

A lower bound on the worst-case error of order $(\log N)^{(s-1)/2}/N^{\alpha}$, which holds for any (non-linear/adaptive) quadrature rule based on $N$ function evaluations, can be proven by adapting the bump function technique from [2]. We also refer to [22, Theorem 4 and Appendix] whose result directly applies to the present problem. As we shall show, higher order digital nets and sequences achieve this best possible order exactly.

6.2 Optimal order error bounds

First we discuss what happens if we exploit only the decay of the Walsh coefficients. A similar result has been already proven in [3, Theorem 30] but with the slightly worse exponent of the $\log N$ term.

Theorem 4**.**

Let $\alpha\geq 2$. Let $P$ be an order $\alpha$ digital $(t_{\alpha},m,s)$-net over $\mathbb{F}_{b}$. Then there exists a constant $E_{\alpha,b,s,t_{\alpha}}^{(1)}>0$ such that the following holds:

[TABLE]

Proof.

Using the first assertion of Lemma 10 and the second assertion of Lemma 3 (with $\lambda=2$), we have

[TABLE]

By considering the equality $\mu_{\alpha}(P^{\perp})=\alpha m-t_{\alpha}+1$, we prove the theorem. ∎

Again the exponent $(s\alpha-1)$ of the numerator of the theorem stems from the inequality

[TABLE]

given in the first assertion of Lemma 3. It seems hard to fundamentally improve this bound. Thus, even before exploiting the sparsity of the Walsh coefficients, there is a difficulty in obtaining the best possible order of the shift-averaged worst-case error.

To overcome this issue, let us go back to the bound given in the first assertion of Lemma 1:

[TABLE]

As we discussed in Subsection 5.3, this gives the optimal order of the shift-averaged $L_{2}$-discrepancy. Considering that the best possible exponent of the $\log N$ term for the present integration problem is the same as that of the $L_{2}$-discrepancy, which is $(s-1)/2$, one idea is to switch the weight function from $\mu_{\alpha}$ to $\mu_{1}$ in the error analysis. The following interpolation inequality was shown in [36] to realize this:

Lemma 11**.**

Let $\alpha,\beta\in\mathbb{N}$ with $1<\alpha<\beta$. For any $\boldsymbol{k}\in\mathbb{N}_{0}^{s}$ we have

[TABLE]

By using this inequality together with the propagation rule of higher order digital nets (Lemma 2), we can improve the order of the shift-averaged worst-case error to best possible. The following result is from [36].

Theorem 5**.**

Let $\alpha\geq 2$. Let $P$ be an order $2\alpha$ digital $(t_{2\alpha},m,s)$-net over $\mathbb{F}_{b}$. Then there exists a constant $E_{\alpha,b,s,t_{2\alpha}}^{(2)}>0$ such that the following holds:

[TABLE]

Proof.

We write $A=(\alpha-1)/(2\alpha-1)$ and $B=\alpha/(2\alpha-1)$. Using the first assertion of Lemma 10, Lemma 11 (with $\beta=2\alpha$) and the second assertion of Lemma 1 (with $\lambda=2B$), we have

[TABLE]

Since $2B-1=1/(2\alpha-1)>0$, the last sum over $z$ is finite. Using Lemma 2 we have

[TABLE]

and so

[TABLE]

from which the result of the theorem follows. ∎

It is important to recall that we consider order $2\alpha$ digital nets in this theorem, instead of order $\alpha$ digital nets as in Theorem 4. The order $2\alpha$ is required here to ensure the finiteness of the sum

[TABLE]

Finally, combining the idea of switching the weight function (Theorem 5) with the idea of exploiting both the decay and the sparsity of the Walsh coefficients simultaneously (Theorem 3), we arrive at a deterministic counterpart of Theorem 5, proven in [38].

Theorem 6**.**

Let $\alpha\geq 2$. Let $P$ be an order $(2\alpha+1)$ digital $(t_{2\alpha+1},m,s)$-net over $\mathbb{F}_{b}$. Then there exists a constant $E_{\alpha,b,s,t_{2\alpha+1}}^{(3)}>0$ such that the following holds:

[TABLE]

Proof.

We write $A=(\alpha-1)/(2\alpha)$ and $B=(\alpha+1)/(2\alpha)$. Using Lemmas 10, 11 (with $\beta=2\alpha+1$) and 4 (with $\alpha$ replaced by $2\alpha+1$) in order and then applying the change of variables $z\to z+2\mu_{1}(P^{\perp})$, we have

[TABLE]

Since $B-1/2=1/(2\alpha)>0$, the last sum over $z$ is finite. Using Lemma 2 we have

[TABLE]

and so

[TABLE]

from which the result of the theorem follows. ∎

The error bound shown in Theorem 6 also applies to the first $b^{m}$ points of an order $(2\alpha+1)$ digital $(t_{2\alpha+1},s)$-sequence over $\mathbb{F}_{b}$ if $m\geq t_{2\alpha+1}/(2\alpha+1)$, since they can be identified with an order $(2\alpha+1)$ digital $(t_{2\alpha+1},m,s)$-net over $\mathbb{F}_{b}$.

7 Conclusions and outlook

In this article we have reviewed some of recent results on HoQMC methods with the particular aim to provide a unified picture on how the Walsh analysis enables these developments. The challenge in analyzing either the $L_{2}$-discrepancy or the worst-case error in a reproducing kernel Hilbert space is that we have to deal with the double sum of the Walsh coefficients. Considering the shift-averaged worst-case error instead, the problem becomes much easier in many cases, but of course, the outcome will remain probabilistic. To obtain a deterministic counterpart of such probabilistic result, exploiting the decay and the sparsity of the Walsh coefficients simultaneously seems to be a reasonable strategy to attack the problem. In fact, as we have seen, both the optimal order $L_{2}$-discrepancy bound in [25] and the optimal order quadrature error bound in [38] are obtained by employing this strategy.

Looking into the future, there are some possible directions for further research as raised below.

Choice of an orthonormal basis: The system of Walsh functions fits quite well with digital nets as emphasized in this article, but is not the only choice. Indeed, recent papers [42, 15, 16, 4] use the system of Haar functions and succeed in generalizing or extending the result of [25]. For instance, it was proven in [15] that order 2 (instead of order 5) digital sequences achieve the best possible order of $L_{2}$-discrepancy. Also, prior to the works of [36, 37, 38], Hinrichs et al. used the system of Faber functions to analyze the worst-case error of order 2 digital nets for different function spaces [40]. As natural questions from the current status, one may ask “Can we get better results in numerical integration problems by using the system of Haar/Faber functions? Can we lower the necessary order of digital nets and sequences from $2\alpha+1$ to achieve the best possible error rate?” We do not have any progress on these questions so far. 2. 2.

Alternative construction scheme: Higher order digital nets and sequences can be explicitly constructed through the digit interlacing function (Definition 12 and Lemma 5). Quite recently, it has been shown in [14] that Richardson extrapolation plays an alternative role to the digit interlacing when the class of underlying point sets are restricted to polynomial lattice point sets. Also, the paper [34] proposed a different usage of Richardson extrapolation in the context of HoQMC methods. It is desirable to have more different options for explicit construction of higher order digital nets and sequences. 3. 3.

Universality for various function classes: One major drawback of higher order digital nets and sequences is that we need to construct point sets or sequences depending on dominating mixed smoothness of the considered function space. Propagation rule (Lemma 2) says that, once we construct point sets or sequences which work for a certain smoothness $\beta$, they also work for any smaller smoothness $1\leq\alpha<\beta$ but not for larger smoothness. Ideally what we want in practice is point sets or sequences which work for all ranges of smoothness. One straightforward idea is to construct infinite order digital nets and sequences and then study their propagation rule. In this line of research, we refer to [62, 67] for the Walsh analysis of infinitely many times differentiable functions, and furthermore, to [44, 59, 60, 61, 13, 43] for the relevant literature.

Acknowledgments

This article is based on the talk which the first author gave during the discrepancy workshop of RICAM special semester “Multivariate Algorithms and their Foundations in Number Theory”. He would like to thank the organizers of the workshop, Dmitriy Bilyk, Josef Dick, and Friedrich Pillichshammer, for their kind invitation. The authors sincerely acknowledge their colleague Takehito Yoshiki. Some of the outcomes brought from the collaboration with him play a central part of this article. The work of T. G. is supported by JSPS Grant-in-Aid for Young Scientists No. 15K20964. The work of K. S. is supported by Grant-in-Aid for JSPS Fellows No. 17J00466.