# 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 diﬀerent 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).$

Let

$\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})}$.

• 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).