Construction Algorithms for Higher Order Polynomial Lattice Rules

We recently submitted the manuscript

This paper fits into the work on higher order quasi-Monte Carlo rules which started with [22] and [31]. It can be viewed as the higher order extension of [7], where classical polynomial lattice rules were considered.

1. Construction Algorithms

As stated in the title of the manuscript, we present construction algorithms for higher order polynomial lattice rules. The construction is based on the worst-case error rather than the quality parameter {t}. This allows us to find good higher order polynomial lattice rules for weighted function spaces. It also presents a feasible alternative to the direct construction introduced in [31] (and [22] for the periodic case).

We use a variety of approaches to achieve our results. One is the component-by-component (cbc) construction. This is a greedy algorithm where one chooses the component of the generating vector {\boldsymbol{q} = (q_1,\ldots, q_d)} starting with {d = 1} and increasing the dimension by one in each step.

One of the differences to the cbc algorithm for classical lattice rules and polynomial lattice rules is that one cannot simply set {q_1 = 1}. The reason for this is that even in the one dimensional case the construction of a quasi-Monte Carlo (qMC) rule which achieves higher order convergence of the worst-case error is not trivial (in the classical setting one simply chooses the points {0, 1/N, 2/N, \ldots, (N-1)/N}). Hence one needs to search also for a good one-dimensional rule.

Further, one also wants to obtain higher order polynomial lattice rules which can be used if the smoothness of the integrand is not known exactly. To find such vectors we use a sieve principle. In its simplest form, the idea is the following: Let {A \subset \mathbb{N}} be a finite set with {a} elements and let {B, C \subseteq A}, where {B} has {b} elements and {C} has {c} elements. The sieve principle then just states that if {b, c > a/2}, then {B \cap C \neq \emptyset}. This method can be generalised to more than {2} sets, which is what we use in the paper.

Note that there is a propagation rule (Propagation Rule II in [44]) which roughly states that a higher order net which is optimal for a smoothness {\alpha} is also automatically optimal for all functions with smoothness {\alpha'} with {\alpha' \le \alpha}. Unfortunately this propagation rule applies only to the quality parameter {t}. An analogous propagation rule for the worst-case error is not known. Since the cbc construction is not based on the quality parameter {t} (we do not have any information about the {t} value of the higher order polynomial lattice rules constructed using the cbc approach in our manuscript), we resort to the sieve principle to find higher order polynomial lattice rules which are optimal for a range of smoothness classes.

The construction cost of the fast cbc algorithm is {O(n b^n s)} operations. Since {n} is ideally of size {\alpha m}, where {b^m} is the number of points, one quickly reaches the limits of computational feasibility. For instance, assume {b = 2} and we want to have, say, {2^{13} \approx 10^4} points. If we want quadratic convergence ({\alpha = 2}) then {b^n = 2^{2 \cdot 13} = 2^{26} \approx 6 \cdot 10^7}. For {\alpha = 3} we obtain {2^{39} \approx 10^{12}}. Hence it is difficult to find point sets with more than say {10^4} points, even for {\alpha = 2}.

One way to reduce the construction cost further would be to construct classical polynomial lattice rules (i.e. {\alpha = 1}) and then use the method from [31] to obtain higher order polynomial lattice rules. In this case the {\alpha} does not appear in the exponent as in {b^{n} \approx b^{\alpha m}}, but rather one has to increase the dimension by a factor of {\alpha}, that is, one needs to construct classical polynomial lattice rules with {b^m} points in dimension {s \alpha}. Hence one can expect a construction cost of {O(m b^m \alpha s)}. How well such an approach would work is currently not known.

2. Higher order nets: how well do they measure up?

Since the direct construction method of [31] is currently the only method to yield applicable higher order digital nets, it is not know how well it measures up. That is, it is not known whether there is still significant room for improvement of the direct construction method. Numerical investigations of the method presented in this manuscript will shed more light on whether the direct construction method works well for point sets whose size is within the range of the present construction method.

In order to be better able to judge how well higher order polynomial lattice rules work, one could also perform a numerical comparison in the periodic case. In Korobov spaces, the natural candidate for numerical integration is a lattice rule, which can achieve arbitrarily high rates of convergence. The construction of higher order polynomial lattice rules for periodic functions is somewhat simpler, since in this case one can always choose {n = m}. Hence it should be possible to construct periodic higher order polynomial lattice rules for a much larger number of points in much higher dimension. Until now, the higher order polynomial lattice rules for periodic functions have not been studied (though the higher order polynomial lattice rules which work for non-periodic functions also work for periodic ones).

Numerical comparisons with sparse grids would also be interesting, but such comparisons have not been performed as of now.

3. An Ansatz

It is interesting to note that higher order polynomial lattices are constructed differently from the construction of higher order nets from [31]. The construction from [31] uses classical digital {(t,m,\alpha s)}-nets to obtain digital {(t,\alpha,1, \alpha m \times m, s)}-nets.

On the other hand, a higher order polynomial lattice is defined in the following way: Take a classical polynomial lattice {\{\boldsymbol{x}_0,\ldots, \boldsymbol{x}_{b^{\alpha m}-1}\}} with {b^{\alpha m}} points. Then the first {b^m} points {\{\boldsymbol{x}_0,\ldots, \boldsymbol{x}_{b^m-1}\}} form a higher order polynomial lattice. In the construction algorithms though, one uses a different search criterion in the classical case than in the present paper. Hence the construction is different in this sense. That is to say, one cannot expect that the first {b^m} points of a classical {(t,\alpha m, s)}-net will achieve higher order of convergence.

It would be interesting though, to have a construction for higher order nets which uses a classical digital {(t,\alpha m, s)}-net to construct a digital {(t',\alpha, 1, \alpha m \times m, s)}-net. This would be a variation of the propagation rule described in the post on duality theory and propagation rules, which can be found here. The difficulty here is that the linear independence condition for the {(t,\alpha m, s)}-net is different from the linear independence condition of the {(t',\alpha, 1, \alpha m \times m, s)}-net and that each row has only {m} entries in the higher order case (compared to {\alpha m} in the classical case). How these difficulties can be overcome is currently not known.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s