STAR DISCREPANCY BOUNDS BASED ON HILBERT SPACE FILLING CURVE STRATIFIED SAMPLING AND ITS APPLICATIONS

Many problems in industry, statistics, physics and finance need to introduce multivariable integrals. In order to effectively represent, analyze and process these integrals on the computer, it is necessary to discretize the sample points. A natural problem is whether these discrete sampling values can achieve the optimal fast convergence approximation of multivariable integrals. This is investigated by the discrepancy function method and sampling theory. And which also makes it a popular research topic in applied mathematics and engineering science.

For the multivariate function integral $\int_{[0,1]^{d}}f(x)dx$ defined on a $d-$dimensional unit cube, if a discrete sampling point set $P_{N,d}=\{t_{1},t_{2},\ldots,t_{N}\}$ is selected, it can be approximated by the sample mean $\frac{1}{N}\sum_{t\in P_{N,d}}f(t)$. Then the approximation error can be calculated by:

Among various techniques to estimate the error of $E_{N}(f)$, star discrepancy has proven to be among the most efficient approaches. Star discrepancy was first used to measure the irregularity of the point distribution. For a point set $P_{N,d}=\{t_{1},t_{2},\ldots,t_{N}\}$, the star discrepancy is defined as follows:

where $\lambda$ denotes the $d$-dimensional Lebesgue measure, and $I_{[0,x]}$ denotes the characteristic function defined on the $d$-dimensional rectangle $[0,x]$.

A well-known result in discrepancy theory is the Koksma-Hlawka inequality [1, 2]. It shows that the integral approximation error of any multivariate function can be separated into two parts, one is the star discrepancy of the defined point set $P_{N,d}$, and the other is the total variation of the function $f$ (under the definition of Hardy and Krause), i.e,:

The equation (1.2) indicates that when $f$ is a function of bounded variation, a smaller upper bound on the star discrepancy yields a more accurate approximation error of the integration. Therefore, we focus on the research of the star discrepancy in the following.

There are many kinds of constructive sampling sets, also known as low discrepancy point sets [3], which obtain a satisfactory convergence rate of star discrepancy bound. For a fixed dimension $d$, it can reach $O((\ln N)^{\alpha_{d}}/N)$, where $\alpha_{d}\geq 0$ is a constant that depends only on the dimension $d$. This kind of point set is widely used in many fields, such as colored noise, learning theory and finance  [4, 5, 6, 7], but this kind of point set also has some limitations. First of all, the convergence rate is strongly dependent on the dimensions. For high dimensional space, the convergence rate becomes slower, even worse than the convergence rate $O(N^{-\frac{1}{2}})$ obtained by classical Monte Carlo in some cases. While the classical Monte Carlo sampling method is independent of the dimensions, so it has strong applicability for especially large-dimensional data. Secondly, because the low discrepancy point set is a constructive sampling, it is difficult to deal with a large number of point sets generated by random sampling in real life, such as the processing of noise-contaminated signals. Therefore, it is necessary to explore a new sampling method that inherits the fast convergence characteristics of low discrepancy point sets in moderate dimensions and the strong applicability of Monte Carlo point sets to high dimensions.

In order to achieve a better approximation effect for multivariate integrals in higher dimensions, a large number of points must be sampled if we use the low discrepancy point sets. Which loses information about the moderate sample size. In other words, if it is required to find the optimal upper bound of the star discrepancy in the higher dimensional space under the moderate sample size, the low discrepancy point set is not applicable. Therefore, in recent years, many scholars have carried out a series of studies on this problem and tried to use the appropriate sampling point set. This study is called the search for the optimal pre-asymptotic star discrepancy bound. In 2001, Heinrich, Novak [8] et al. proved that the optimal pre-asymptotic upper bound can be obtained as long as the number of sampling points linearly depends on the dimension, and the upper bound is $c\cdot\frac{\sqrt{d}}{{\sqrt{N}}}$, where $c$ is an unknown positive constant. In 2011, Aistleitner [9] calculated that $c$ is at most 10. This constant has been improved to 2.53  [10] so far. From a practical point of view, the following probabilistic pre-asymptotic upper bound  [11] can be obtained, that is,

holds with probability at least $q$ for the Monte Carlo point set $Y$.

Combined with the K-H inequality, it is easy to achieve the optimal probability approximation of the bounded variation function on the $d-$dimensional unit cube. Especially for high-dimensional data, this approximation overcomes some disadvantages of the traditional deterministic low discrepancy point sets, and gives a good approximation with a moderate sample size.

By optimizing the cover technique, the probability bound of (1.3) is improved by  [10], which gives

holds with probability at least $q$ for the Monte Carlo point set $Y$.

It is worth mentioning that the above pre-asymptotic star discrepancy bounds and the optimal probability approximation bounds are all obtained by using the Monte Carlo point set. Although the Monte Carlo method can get better approximation for high dimensional data, the convergence speed is slow, and the convergence speed is $O(N^{-\frac{1}{2}})$. How to improve this convergence speed and make the original approximation still applied to moderate sample size is a concern of scholars in recent years. Some scholars have noticed that the main reason for the slow convergence of the Monte Carlo method is the large variance of the sample mean function. Therefore, they have focused on the variance reduction technique and tried to use some new sampling methods to improve the approximation effect. For example, the probabilistic star discrepancy bound using Latin hypercube sampling  [10] and the expected bound  [12, 13] obtained by using jittered sampling. Jittered sampling is formed by isometric grid partition, and $[0,1]^{d}$ is divided into $m^{d}$ small sub-cubes $Q_{i},1\leq i\leq N,$ each with side length $\frac{1}{m},$ Figures 1 and 2 give $2$-dimensional and $3$-dimensional illustrations.

Moreover, although the classical jittered sampling inherits some advantages of stratified sampling, it is more demanding for the number of sampling points, i.e., $N=m^{d}$. The exponential dependence of the sampling number on $d$, which will inevitably lead to the curse of dimensionality, so it is difficult to obtain the pre-asymptotic star discrepancy bound. However, jittered sampling retains the value of theoretical and numerical research, such as, it can achieve a faster convergence rate than the Monte Carlo sampling in a relatively low-dimensional space. While Latin hypercube sampling overcomes the dimensionality problem and reduces the variance of the sample mean, the star discrepancy convergence order is still $O(N^{-\frac{1}{2}})$, see  [10].

If the problem is to reduce the error of the (1.2) under suitable sampling number, then for a class of bounded variation function spaces, there is

holds with probability at least $q$.

So how to improve the convergence speed of the approximation error in (1.5) and preserve approximation information under the appropriate sampling number, based on the above introduction and discussion, there are two main ideas: one is to use the low discrepancy point sets, but note that the convergence order is $O((\ln N)^{\alpha_{d}}/N)$, where $\alpha_{d}\geq 0$ is a constant that depends on the dimension $d$. Generally for different point sets, $\alpha_{d}$ is at least $O(d)$, so for high-dimensional space, the number of samples required is particularly large; the other is to use the idea of stratified sampling, but for classical jittered sampling, the disadvantage is that $N=m^{d}$ sampling points are required, and for high-dimensional space, it faces the curse of dimensionality. Therefore, it is necessary to return to the general stratified sampling method. We try to use a sampling method that combines the advantages of low discrepancy point sets and stratified sampling, and to overcome the dimension problem to a certain extent.

The Hilbert space filling curve sampling method coincidentally inherits many advantages of general stratified sampling, and we describe it in detail below. This method was developed by A. B. Owen, Z. He et al. to study the variance of the sample mean function  [14, 15]. In our main result, we use Hilbert space filling curve sampling to obtain an improved upper bound of probability star discrepancy, and we give

holds with probability at least $q$ for the Hilbert space filling curve sampling set $X$, where $c(d,q)$ is an explicit and computable constant depending on $d$ and $q$.

We can easily derive the probability star discrepancy bounds to the weighted form, and give the uniform probability approximation error in the weighted function space, that is, for a function in the weighted Sobolev space $F_{d,1,\gamma}$, the following

holds with probability at least $\epsilon$ for the Hilbert space filling curve sampling set $X$.

The strong partition principle for the star discrepancy version is an open question mentioned in [12, 20], we use the probabilistic results to prove the following conclusion:

where $\Omega=\{\Omega_{1},\Omega_{2},\ldots,\Omega_{N}\}$ is an arbitrary equivolume partition and $P_{\Omega}$ is the corresponding stratified sampling set, and $Y=\{Y_{1},Y_{2},Y_{3},\ldots,Y_{N}\}$ is simple random $d-$dimensional samples.

We also apply the main result to give explicit bounded expressions for a set of sampling inequalities in the general Hilbert space.

We prove:

holds with probability at least $1-\epsilon$($\epsilon$ the probability that $D_{N}^{*}(X)\leq\frac{c(d,\epsilon)}{N^{\frac{1}{2}+\frac{1}{2d}}}$ must hold), where $0<\gamma\leq 1$, and $f$ belongs to Hilbert space.

In the end, we use the variant K-H inequality, combined with the obtained probabilistic star discrepancy bound, to give an approximation error for the piecewise smooth function $f$ defined on a Borel convex subset of $[0,1]^{d}$.

We prove:

holds with probability at least $q$, where $D_{N}^{*}(X_{1},X_{2},\ldots,X_{N})$ is the star discrepancy of Hilbert space filling curve sampling $X$ and $q$ denotes the probability that the upper bound $D_{N}^{*}(X)\leq\frac{c(d,q)}{N^{\frac{1}{2}+\frac{1}{2d}}}$ must hold. Besides,

The symbol $\frac{\partial^{|u|}}{\partial x_{n}}f(x)$ is the partial derivative of $f$ with respect to those components of $x$ with index in $u$.

The remainder of this paper is organised as follows. Section 2 presents the preliminaries. Section 3 presents our main result. Section 4 gives some applications of the main result. Finally, in section 5, we conclude the paper with a short summary.

1. Hilbert space filling curve sampling

We first introduce Hilbert space filling curve sampling, which is abbreviated as HSFC sampling. We mainly adopt the definitions and symbols in [14, 15]. HSFC sampling is actually a stratified sampling formed by a special partition method.

Let $a_{i}$ be the first $N=b^{m}$ points of the van der Corput sequence(van der Corput 1935) in base $b\geq 2,m=0,1,\ldots$. The integer $i-1\geq 0$ is written in base $b$ as

for $a_{ij}\in\{0,\ldots,b-1\}$. Then, $a_{i}$ is defined by

The scrambled version of $a_{1},a_{2},\ldots,a_{N}$ is $x_{1},x_{2},\ldots,x_{N}$ written as

where $x_{ij}$ are defined through random permutations of the $a_{ij}$. These permutations depend on $a_{ik},$ for $k<j.$ More precisely, $x_{i1}=\pi(a_{i1}),$ $x_{i2}=\pi_{a_{i1}}(a_{i2})$ and generally for $j\geq 2,$

Each random permutation is uniformly distributed over the $b!$ permutation of $\{0,\ldots,b-1\}$, and is mutually independent of the others. Thanks to the nice property of nested uniform scrambling, the data values in the scrambled sequence can be reordered such that

independently with

for $i=1,2,\ldots,N=b^{m}$.

Let

where $H$ is a mapping.

Then

is the corresponding stratified samples. Set $r_{i}$ be the diameter of $E_{i}$, according to the property of HSFC sampling, the following estimation holds,

2. Minkowski content

We use the definition of Minkowski content in [16], which provides convenience for analysing the boundary properties of the test set $B$ in (1.1). For a set $\Omega\subset[0,1]^{d}$, define

where $(\partial\Omega)_{\epsilon}=\{X\in\mathbb{R}^{d}|dist(x,\partial\Omega)\leq\epsilon\}$. If $\mathscr{M}(\partial\Omega)$ exists and is finite, then $\partial\Omega$ is said to admit $(d-1)-$dimensional Minkowski content. If $\Omega$ is a convex set, then it is easy to see that $\partial\Omega$ admits $(d-1)-$dimensional Minkowski content. Furthermore, $\mathscr{M}(\partial\Omega)\leq 2d$ as the outer surface area of a convex set in $[0,1]^{d}$ is bounded by the surface area of the unit cube $[0,1]^{d}$, which is $2d$.

3. $\delta$-covers

To discretise the star discrepancy, we use the definition of $\delta-$covers as in[17].

From [18], combined with Stirling’s formula, the following estimate holds for $\mathcal{N}(d,\delta)$, that is, for any $d\geq 1$ and $\delta\in(0,1)$,

Furthermore, the following lemma provides a convenient way to estimate the star discrepancy with $\delta-$covers.

4. Bernstein inequality

At the end of this section, we will restate the Bernstein inequality which will be used in the estimation of star discrepancy bounds.

The following will give three applications of the theorem  3.1. One is the application of uniform integral approximation in weighted function space, and the other is to prove the strong partition principle of star discrepancy version. The strong partition principle is the problem proposed by  [20, 21], which has been perfectly solved for the $L_{2}-$discrepancy, but remains to be solved for the star discrepancy version. The third is the sampling theorem in general Hilbert space. Finally, it is applied to the integral approximation problem on the Borel convex subset in $[0,1]^{d}$.