On the fast component-by-component algorithm for polynomial lattice rules

We present some simplification and observations concerning the fast component-by-component algorithm for polynomial lattice rules. This algorithm was introduced by Nuyens and Cools in

D. Nuyens and R. Cools. Fast component-by-component construction, a reprise for different kernels. In H. Niederreiter and D. Talay, editors, Monte Carlo and Quasi-Monte Carlo Methods 2004, pages 371–385. Springer-Verlag, 2006. Cited on pp. 82, 153, 178. See also Chapter 6 of Dirk’s PhD Thesis.

Let {\boldsymbol{x} = (x_1,\ldots, x_s)} and {x_i = \xi_{i,1} b^{-1} + \xi_{i,2} b^{-2} + \cdots} be the {b} adic expansion of {x} (unique in the sense that infinitely many of the {\xi_{i,j}} are different from {b-1}). We use the notation as in the book. Assume that {b} is a prime number. For the sake of definiteness, consider the mean square worst-case error

\displaystyle e^2(\boldsymbol{q},p) = - 1 + \frac{1}{b^m} \sum_{n=0}^{b^m-1} \prod_{i=1}^s (1 + \gamma_i \phi_b(\boldsymbol{x}_n)),

where {\boldsymbol{x}_0,\ldots, \boldsymbol{x}_{b^m-1}} is a polynomial lattice rule, {\phi_b(\boldsymbol{x}) = \prod_{i=1}^s \phi_b(x_i)} for {\boldsymbol{x} = (x_1,\ldots, x_s)}, and

\displaystyle \phi_b(x) = \left\{\begin{array}{ll} \frac{1}{6} & \mbox{if } x = 0, \\ \frac{1}{6} - \frac{\xi_{j_0}(b-\xi_{j_0})}{b^{j_0 + 1}} & \mbox{if } \xi_{j_0} \neq 0, \xi_{1} = \ldots = \xi_{j_0-1} = 0. \end{array}\right.

This criterion does not depend on {\xi_{j_0}} if {b = 2} or {3}, only on {j_0}. This does not hold for other values of {b}. In the following we assume that {b \in \{2,3\}}, though the ideas can be generalised to arbitrary bases. We are using the component-by-component algorithm. Assume that {q_1,\ldots, q_{d-1}} are already chosen. Then we define a function of {q \in \{g \in \mathbb{F}_b[x]\setminus \{0\}: \deg(g) < m\}} given by

\displaystyle \begin{array}{lcl} C_1(q) & := & e^2((q_1,\ldots, q_{d-1},q),p) \\ && \\ & = & -1 + \frac{1}{b^m} \sum_{n=0}^{b^m-1} (1+ \gamma_d \phi_b(x_{n,d})) \prod_{i=1}^{d-1} (1 + \gamma_i \phi_b(x_{n,i})). \end{array}

As only {x_{n,d}} depends on {q} for {0 < n < b^m}, we can introduce some simplifications. Let {\mu_{n,d} = \prod_{i=1}^{d-1} (1+ \gamma_i \phi_b(x_{n,i}))}. Then we can write

\displaystyle \begin{array}{lcl} C_1(q) & = & -1 + \frac{1}{b^m} (1+ \gamma_d \phi_b(x_{0,d})) \mu_{0,d} + \frac{1}{b^m} \sum_{n=1}^{b^m-1} \mu_{n,d} \\ && \\ && + \frac{\gamma_d}{b^m} \sum_{n=1}^{b^m-1} \phi_b(x_{n,d}) \mu_{n,d}.  \end{array}

The aim is to find {q} which minimises {C_1(q)}. Note that only {x_{n,d}} depends on {q}, hence we can also find {q} which minimises

\displaystyle C_2(q) = \sum_{n=1}^{b^m-1} \phi_b(x_{n,d}) \mu_{n,d}.

We can simplify {C_2} further. Using the definition of {\phi_b} it follows that the above is equivalent to finding {q} which minimises

\displaystyle C_3(q) = - \sum_{n=1}^{b^m-1} b^{-j_{0,n,q}} \mu_{n,d},

where {j_{0,n,q}} is the smallest integer {j} such that {\xi_{n,i,j} \neq 0} with {x_{n,i} = \xi_{n,i,1} b^{-1} + \xi_{n,i,2} b^{-2} + \cdots}. Note that {x_{n,i} \in \{b^{-m}, 2 b^{-m}, \ldots, (b^{m}-1)b^{-m} \}} for {0 < n < b^m}. Let {\boldsymbol{C}_3 = (C_3(1), \ldots, C_3(b^m-1))^\top}, {\boldsymbol{\mu}_d = (\mu_{1,d},\ldots, \mu_{b^m-1,d})^\top}, and

\displaystyle A_1 = \left(\begin{array}{llll} b^{-j_{0,1,1}} & \ldots & \ldots & b^{-j_{0,1,b^m-1}} \\ b^{-j_{0,2,1}} & \ldots & \ldots & b^{-j_{0,2,b^m-1}} \\ \vdots & & & \vdots \\ b^{-j_{0,b^m-1,1}} & \ldots & \ldots & b^{-j_{0,b^m - 1,b^m-1}} \end{array} \right).

In this notation we have {\boldsymbol{C}_3 = - A_1^\top \boldsymbol{\mu}_d}.

Assume that for an {g \in \mathbb{F}_b[x]} we have

\displaystyle \frac{g(x)}{p(x)} = u_{1,g} x^{-1} + u_{2,g} x^{-2} + \cdots,

where {u_{j,g} \in \mathbb{F}_b}. Then

\displaystyle x_{n,d} = u_{1,n q} b^{-1} + \cdots + u_{m,n q} b^{-m},

where {n, q \in \mathbb{F}_b[x]} are the polynomials associated with {n, q}, and the multiplication is carried out in {\mathbb{F}_b[x]}. Thus {j_{0,n,q}} is the smallest integer {j} such that {u_{j, nq} \neq 0}, where {n,q \in \mathbb{F}[x]}.
For simplicity let us assume that {p} is primitive. Then the monomial {x} is a generating element of the multiplicative group of {\mathbb{F}_b[x]/(\! \pmod p)}. Let {r_{k+l} = j_{0,x^k, x^l}}.

Applying the Rader transform we obtain a circulant matrix

\displaystyle A_2 = \left(\begin{array}{lllll} b^{-r_{0}} & b^{-r_{1}} & \ldots & b^{-r_{b^m-3}} & b^{-r_{b^m-2}} \\ b^{-r_{b^m-2}} & b^{-r_{0}} & \ddots & & b^{-r_{b^m-3}} \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ b^{-r_{2}} & & \ddots & \ddots & b^{-r_{1}} \\ b^{-r_{1}} & b^{-r_{2}} & \ldots & b^{-r_{b^m-2}} & b^{-r_{0}} \end{array} \right).

Now {r_k} is the smallest integer {r} such that {u_{r, x^k} \neq 0}. We have

\displaystyle \deg(x^k \pmod{p}) = \deg(p) - r_k

and hence

\displaystyle r_k = m - \deg(x^k \pmod{p}).

Hence we can simplify the matrix {A_2} further. Let {t_k = \deg(x^k \pmod{p})}, hence {t_k \in \{0,1, \ldots, m-1\}}. Then let

\displaystyle A_3 = \left(\begin{array}{lllll} b^{t_0} & b^{t_{1}} & \ldots & b^{t_{b^m-3}} & b^{t_{b^m-2}} \\ b^{t_{b^m-2}} & b^{t_{0}} & \ddots & & b^{t_{b^m-3}} \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ b^{t_{2}} & & \ddots & \ddots & b^{t_{1}} \\ b^{t_{1}} & b^{t_{2}} & \ldots & b^{t_{b^m-2}} & b^{t_{0}} \end{array} \right).


\displaystyle C_4 = -A_3^\top \boldsymbol{\mu'}_d,

where {C_4 = (c_1,\ldots, c_{b^m-2})^\top} and \boldsymbol{\mu'}_d, is the permuted vector \boldsymbol{\mu}_d. Then we can find {q} as

\displaystyle q = \mathrm{argmin}_{1 \le g \le b^m-2} \; c_g.

The calculation of the matrix {A_3} is simpler because we do not need to calculate the Laurent series expansion of {\frac{g(x)}{p(x)}}.

The same results hold if {p} is irreducible but not primitive. In this case one has to replace {x} by a generating element of the cyclic multiplicative group of {\mathbb{F}_b[x]/(\pmod{p})}.

We add some further observations:

  • The entries in the matrix {A_3} are from the set {1, b, b^2, \ldots, b^{m-1}}. Hence the matrix has only {m} different entries (the total number of entries in {A_3} is {(b^m-1)^2}, which is much larger). This was also pointed out in Dirk’s PhD Thesis, p. 150.
  • In each row and for each {0 \le l < m}, the entry {b^l} occurs {(b-1) b^l} times. The last statement holds since there are {(b-1)b^l} polynomials of degree {l} in {\mathbb{F}_b[x]}.
  • All the entries of the matrix {A_3} are of the form {b^{t_k}}, hence to multiply the vector {\boldsymbol{\mu}_d} with {A_4} we only need to shift the numbers of the vectors {t_k} digits rather instead of multiplying real numbers. Hence in calculating {A_3 \boldsymbol{\mu}_d} no multiplication is required.
  • Let {B = (b^{m-1})_{1 \le k,l \le b^m-2}} be the matrix whose entries are all {b^{m-1}}. Then instead of using {A_3} in the minimisation problem above, one could use {A_4 = B-A_3}. The matrix {A_4} is then a circulant matrix which has {(b-1) b^{m-1} (b^m-1)} entries which are {0}.

The multiplication {A_3 \boldsymbol{\mu}_d} is carried out using the fast Fourier transform, which does not take advantage of these observations.
It is not known whether there is an alternative algorithm which calculates {A_3 \boldsymbol{\mu}_d} in at most {O(b^m m)} operations (and which takes advantage of some of the observations above).


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