Subgroup-based Rank-1 Lattice Quasi-Monte Carlo

A $d$-dimensional lattice $\Lambda$ is a set of points that contains no limit points and satisfies [22]

A widely known lattice is the unit lattice $\mathbb{Z}^{d}$ whose components are all integers. A general lattice is constructed by a generator matrix. Given a generator matrix $\boldsymbol{B}\in\mathbb{R}^{d\times d}$, a $d$-dimensional lattice $\Lambda$ can be constructed as

A generator matrix is not unique to a lattice $\Lambda$, namely, a lattice $\Lambda$ can be obtained from a different generator matrices.

A lattice point set for integration is constructed as $\Lambda\cap[0,1)^{d}$. This step may require an additional search (or check) for all the points inside the unit cube.

A rank-1 lattice is a special case of the general lattice, which has a simple operation for point set construction instead of directly using Eq.(2). A rank-1 lattice point set can be constructed as

where $\boldsymbol{z}\in\mathbb{Z}^{d}$ is the so-called generating vector, and the big $\langle\cdot\rangle$ denotes the operation of taking the fractional part of the input number elementwise. Compared with the general lattice rule, the construction form of the rank-1 lattice already ensures the constructed points to be inside the unit cube without the need for any further checks.

Given a rank-1 lattice set $X$ in the unit cube, we can also construct a randomized point set. Sample a random variable $\boldsymbol{\Delta}\sim Uniform[0,1]^{d}$, we can construct a point set $\widetilde{X}$ by random shift as [8]

Several criteria have been studied in the literature for good lattice construction through computer search. Worst case error is one of the most widely used criteria for functions in a reproducing kernel Hilbert space (RKHS) [8]. However, this criterion requires the prior knowledge of functions and the assumption of the RKHS. Without assumptions of the functions, it is reasonable to construct a good lattice by designing an evenly spaced point set. Minimizing the covering radius is a good way for evenly spaced point set construction.

As minimizing the covering radius of the lattice is equivalent to maximizing the packing radius [7], we can construct the point set through maximizing the packing radius (separating distance) of the lattice. Define the covering radius and packing radius of a set of points $X=\{x_{1},...,x_{N}\}$ as Eq.(5) and Eq.(6), respectively:

The $l_{p}$-norm-based toroidal distance [11] between two lattice points ${\bf{x}}\in[0,1]^{d}$ and ${\bf{x}^{\prime}}\in[0,1]^{d}$ can be defined as:

Because the difference (subtraction) between two lattice points is still a lattice point, and a rank-1 lattice has a period 1, the packing radius $\rho_{X}$ of a rank-1 lattice can be calculated as

where $\|{\bf{x}}\|_{T_{2}}$ denotes the $l_{2}$-norm-based toroidal distance between $\bf{x}$ and $\bf{0}$, symbol $X\setminus{\bf{0}}$ denotes the set $X$ excludes the point $\boldsymbol{0}$. This formulation calculates the packing radius with a time complexity of $\mathcal{O}(Nd)$ rather than $\mathcal{O}(N^{2}d)$ in pairwise comparison. However, the computation of the packing radius is not easily accelerated by fast Fourier transform due to the minimum operation in Eq.(8).

From the definition of rank-1 lattice, we know the packing radius of rank-1 lattice with $n$ points can be reformulated as

where

Then, we have

where $\text{min}(\cdot,\cdot)$ denotes the elementwise min operation between two inputs.

Suppose $n$ is a prime number, from number theory, we know that for a primitive root $g$, the residue of $\{g^{0},g^{1},\cdots,g^{n-2}\}$ modulo $n$ forms a cyclic group under multiplication, and $g^{n-1}\equiv 1\;\text{mod}\;n$. Since $(g^{\frac{n-1}{2}})^{2}=g^{n-1}\equiv 1\;\text{mod}\;n$, we know that $g^{\frac{n-1}{2}}\equiv-1\;\text{mod}\;n$.

Because of the one-to-one correspondence between the residue of $\{g^{0},g^{1},\cdots,g^{n-2}\}$ modulo $n$ and the set $\{1,2,\cdots,n-1\}$, we can construct the generating vector as

without loss of generality, where $m_{1},\cdots,m_{d}$ are integer components to be designed. Denote $\boldsymbol{\bar{z}}=[g^{\frac{n-1}{2}+m_{1}},g^{\frac{n-1}{2}+m_{2}},\cdots,g^{\frac{n-1}{2}+m_{d}}]\;\text{mod}\;n$, maximizing the separating distance $\rho_{X}$ is equivalent to maximizing

Suppose $2d$ divides $n-1$, i.e., $2d|(n-1)$, by setting $m_{i}=g^{\frac{(i-1)(n-1)}{2d}}$ for $i\in\{1,\cdots,d\}$, we know that $H=\{g^{m_{1}},g^{m_{2}},\cdots,g^{m_{d}},g^{\frac{n-1}{2}+m_{1}},g^{\frac{n-1}{2}+m_{2}},\cdots,g^{\frac{n-1}{2}+m_{d}}\}$ is equivalent to setting $\{g^{0},g^{\frac{n-1}{2d}},\cdots,g^{\frac{(2d-1)(n-1)}{2d}}\}$ mod $n$ , and it forms a subgroup of the group $\{g^{0},g^{1},\cdots,g^{n-2}\}$ mod $n$. From Lagrange’s theorem in group theory [10], we know that the cosets of the subgroup $H$ partition the entire group $\{g^{0},g^{1},\cdots,g^{n-2}\}$ into equal-size, non-overlapping sets, and the number of cosets of $H$ is $\frac{n-1}{2d}$. Thus, we know that distance $\text{min}(g^{k}\boldsymbol{z}\;\text{mod}\;n,g^{k}\boldsymbol{\bar{z}}\;\text{mod}\;n)$ for $k\in\{0,\cdots,n-2\}$ has $\frac{n-1}{2d}$ different values, and there are the same numbers of items for each value.

Thus, we can construct the generating vector as

In this way, the constructed rank-1 lattice is more regular as it has few different distinct pairwise distance values, and for each distance, the same number of items obtain this value. Usually, the constructed regular lattice is more evenly spaced, and it has a large minimum pairwise distance. We confirm this empirically through extensive experiments in Section 5.

We summarize our construction method and the properties of the constructed rank-1 lattice in Theorem 1.

As shown in Theorem 1, our method can construct regular rank-1 lattice through a very simple closed-form construction, which does not require any exhaustive computer search.

We show the regular property of rank-1 lattices in terms of $l_{p}$-norm-based toroidal distance.

Theorem 2 gives the upper and lower bounds of the minimum pairwise distance of any non-degenerate rank-1 lattice. The term ‘non-degenerate’ means that the elements in the generating vector are not equal, i.e., $z_{i}\neq z_{j}$.

We now show that our subgroup-based rank-1 lattice can achieve the optimal minimum pairwise distance when $n=2d+1$ is a prime number.

Corollary 1 shows a case when our subgroup rank-1 lattice obtains the maximum minimum pairwise toroidal distance. It is useful for expensive black-box functions, where the number of function queries is small. Empirically, we find that our subgroup rank-1 lattice can achieve near-optimal pairwise toroidal distance in many other cases.

We evaluate the minimum distance of our subgroup rank-1 lattice by comparing with Hua’s method [13] and the Korobov [17] searching method. We denote ‘SubGroup’ as our subgroup rank-1 lattice, ‘Hua’ as rank-1 lattice constructed by Hua’s method [13], and ‘Korobov’ as rank-1 lattice constructed by exhaustive computer search in Korobov form [17].

We set the dimension $d$ as in $\{50,100,200,500\}$. For each dimension $d$, we set the number of points $n$ as the first ten prime numbers such that $2d$ divides $n\!-\!1$, i.e., $2d\big{|}(n\!-\!1)$. The minimum $l_{1}$-norm-based toroidal distance and the minimum $l_{2}$-norm-based toroidal distance for each dimension are reported in Table 5.1 and Table 2, respectively. The larger the distance, the better.

We can observe that our subgroup rank-1 lattice achieves consistently better (larger) minimum distances than Hua’s method in all the cases. Moreover, we see that subgroup rank-1 lattice obtains, in 20 out of 40 cases, the same $l_{2}$-norm-based toroidal distance and in 24 out of 40 cases the same $l_{1}$-norm-based toroidal distance compared with the exhaustive computer search in Korobov form. The experiments show that our subgroup rank-1 lattice achieves the optimal toroidal distance in exhaustive computer searches in Korobov form in over half of all the cases. Furthermore, the experimental result shows that our subgroup rank-1 lattice obtains a competitive distance compared with the exhaustive Korobov search in the remaining cases. Note that our subgroup rank-1 lattice is a closed-form construction which does not require computer search, making our method more appealing and simple to use.

Time Comparison of Korobov searching and our sub-group rank-1 lattice. The table below shows the time cost (seconds) for lattice construction. The run time for Korobov searching grows fast to hours. Our method can run in less than one second, achieving a $10^{4}\times$ to $10^{5}\times$ speed-up. The speed-up increases when $n$ and $d$ becomes larger.

We evaluate our subgroup rank-1 lattice on the integration test problem

We compare with i.i.d. Monte Carlo, a Hua’s rank-1 lattice [13], Korobov searching rank-1 lattice [16], Halton sequence, and Sobol sequence  [8]. For both Halton sequence and Sobol sequence, we use the scrambling technique suggested in [8]. For all the QMC methods, we use the random shift technique as in Eq.(4).

We fix $b=2$ and $c=1$ in all the experiments. We set dimension $d=100$ and $d=500$, respectively. We set the number of points $n$ as the first ten prime numbers such that $2d$ divides $n\!-\!1$, i.e., $2d\big{|}(n\!-\!1)$.

The mean approximation error ($\frac{|Q(f)-I(f)|}{|I(f)|}$) with error bars over 50 independent runs for each dimension $d$ is presented in Figure 2. We can see that Hua’s method obtains a smaller error than i.i.d Monte Carlo on the 50-d problem, however, it becomes worse than MC on 500-d and 1000-d problems. Moreover, our subgroup rank-1 lattice obtains a consistent smaller error on all the tested problems than Hua and MC. In addition, our subgroup rank-1 lattice achieves a slightly better performance than Halton, Sobol and Korobov searching method.

We evaluate the performance of subgroup rank-1 lattice on kernel approximation tasks by comparing with other QMC baseline methods. We test the kernel approximation of the Gaussian kernel, the zeroth-order arc-cosine kernel, and the first-order arc-cosine kernel as in [6].

We compare subgroup rank-1 lattice with a Hua’s rank-1 lattice [13], Halton sequence, Sobol sequence  [8] and standard i.i.d. Monte Carlo sampling. For both the Halton sequence and Sobol sequence, we use the scrambling technique suggested in [8]. For both subgroup rank-1 lattice and Hua’s rank-1 lattice, we use the random shift as in Eq.(4). We evaluate the methods on the DNA [28] and the SIFT1M [14] dataset over 50 independent runs. Each run contains 2000 random samples to construct the Gram matrix. The bandwidth parameter of Gaussian kernel is set to 15 in all the experiments.

The mean Frobenius norm approximation error (${\|{\widetilde{K}\!-\!K}\|}_{F}/{\left\|K\right\|}_{F}$) and maximum norm approximation error (${\|{\widetilde{K}\!-\!K}\|}_{\infty}/{\|K\|}_{\infty}$) with error bars on DNA [28] dataset are plotted in Figure 3. The results on SIFT1M [14] is given in Figure 6 in the supplement. The experimental result shows that subgroup rank-1 lattice consistently obtains a smaller approximation error compared with other baselines.

For general Boltzmann machines with continuous state in $[0,1]$, the energy function of $\boldsymbol{x}\in[0,1]^{d}$ is defined as $E(\boldsymbol{x})=-(\boldsymbol{x}^{\top}\boldsymbol{W}\boldsymbol{x}+\boldsymbol{b}^{\top}\boldsymbol{x})/d$. The normalization constant is $Z=\int_{[0,1]^{d}}{\exp{(-E(\boldsymbol{x})})d\boldsymbol{x}}$. For inference, the marginal likelihood of observation $\boldsymbol{v}\in\mathbb{R}^{d}$ is $\mathcal{L}(\boldsymbol{v})=\int_{[0,1]^{d}}{\exp{(-f(\boldsymbol{v}))}\exp{(-E(\boldsymbol{h})})/Z\text{d}\boldsymbol{h}}$ with function $f(\boldsymbol{v})=-(\boldsymbol{v}^{\top}\boldsymbol{W}_{v}\boldsymbol{v}+2\boldsymbol{v}^{\top}\boldsymbol{W}_{h}\boldsymbol{h}+\boldsymbol{b}_{v}^{\top}\boldsymbol{v})/d$, where $\boldsymbol{h}\in\mathbb{R}^{d}$ denotes the hidden states.

We evaluate our method on approximation of the normalization constant and inference by comparing with i.i.d. Monte Carlo (MC), slice sampling (SS) and Hamiltonian Monte Carlo (HMC). We generate the elements of $\boldsymbol{W}$, $\boldsymbol{W}_{v}$, $\boldsymbol{W}_{h}$, $\boldsymbol{b}$ and $\boldsymbol{b}_{v}$ by sampling from standard Gaussian $\mathcal{N}(0,1)$. These parameters are fixed and kept the same for all the methods in comparison. For inference, we generate an observation $\boldsymbol{v}\in[0,1]^{d}$ by uniformly sampling and keep it fixed and same for all the methods. For SS and HMC, we use the slicesample function and hmcSampler function in MATLAB, respectively. We use the approximation of i.i.d. MC with $10^{7}$ samples as the pseudo ground-truth. The approximation errors $|\widehat{Z}-Z|/Z$ and $|\widehat{\mathcal{L}}-\mathcal{L}|/\mathcal{L}$ are shown in Fig.4(a)-4(d) and Fig.4(e)-4(h), respectively. our method consistently outperforms MC, HMC and SS on all cases. Moreover, our method is much cheaper than SS and HMC.

Comparison to sequential Monte Carlo. When the positive density region takes a large fraction of the entire domain, our method is very competitive. When it is only inside a small part of a large domain, our method may not be better than sequential adaptive sampling. In this case, it is interesting to take advantage of both lattice and adaptive sampling. E.g., one can employ our subgroup rank-1 lattice as a rough partition of the domain to find high mass regions, then take sequential adaptive sampling on the promising regions with the lattice points as the start points. Also, it is interesting to consider recursively apply our subgroup rank-1 lattice to refine the partition. Moreover, our subgroup-based rank-1 lattice enables black-box evaluation without the need for gradient information. In contrast, several sequential sampling methods, e.g., HMC, need a gradient of density function for sampling.