pLab



Home

Spectral-Test Home



LLL-Spectral Test

C-Source

Math-Source

Maintained by  Karl.Entacher@fh-sbg.ac.at  



Spectral Test Server

LLL - Version


The sections below shows how to perform efficient lattice parameter searches for Monte Carlo - and quasi Monte Carlo node sets. The classical measure for such parameter searches is the spectral test which is based on a calculation of the shortest vector in a lattice. Instead of the shortest vector we discuss an application of an approximation given by the LLL algorithm for lattice basis reduction. We empirically demonstrate the speed up and the quality loss obtained by the LLL reduction. For an implementation of the LLL-spectral test in C choose the menue item C - Source.


Spectral Test Approximation with the LLL-Algorithm

Table of Contents:


Application and Analysis of Lattices for Monte Carlo - and Quasi Monte Carlo Methods


Monte Carlo (MC)

Classical node sets for Monte Carlo (MC) methods are obtained by linear congruential random number generators (LCGs). LCGs have extensively been applied for a long time and they are the most common random number generators. Although we have to mention that recent versions use implementations based on a combination of LCGs, or multiple recursive generators (MRGs) to get improved quality and huge periods. For many classical and recent examples, see [8,16]. The definitions and basic properties of linear random number generators are contained in [11,14,15,22]. LCGs are generated by means of the recursion $y_{n+1}\equiv a y_n + b\!\!\pmod m$, $n\ge 0$, and by an initial seed $y_0$, $a\not=1$, $b$, $y_0 \in \mathbb{Z}_m$ (the least residue system modulo $m$). Normalized PRNs in $[0,1[$ are obtained by the transformation $x_n:=y_n / m$. We consider only multiplicative where the modulus $m$ is prime LCGs ($b = 0$) and the multiplier $a$ is a primitive root modulo $m$. Therefore the recursion above, with seed $y_0\not= 0$, produces a sequence of integers in $\mathbb{Z}_m$ with maximal period $m-1$.

A central property of linear congruential generators in general (this holds also for combined LCGs and for MRGs), is that arbitrary $s$-dimensional vectors

\begin{displaymath} {\vec{x}}_i := (x_i, x_{i+j_1}, \ldots , x_{i+j_{s-1}}), \end{displaymath} (1)

with fixed lags $j_1$, $\ldots$ , $j_{s-1}$, are contained in certain grid structures [18].

For $j_1 = 1$, $\ldots$, $j_{s-1} = s-1$, the case which has been studied in detail, these vectors are called overlapping $s$-tuples. In our case, for multiplicative LCGs, the latter $s$-tuples produce intersections of a lattice $L_s(a,m)$ with the $s$-dimensional unit cube $I^s$, where

\begin{displaymath} L_s(a,m) := \left\{ \sum_{i=1}^s k_i \cdot {\vec{b}}_i \; : \; k_i \in \mathbb{Z}\right\}, \end{displaymath} (2)

denotes a $s$-dimensional lattice with lattice basis
\begin{displaymath} {\vec{b}}_1 = (1,a,\ldots,a^{s-1})/m, \; {\vec{b}}_2 = (0,1,0,\ldots,0), \ldots , {\vec{b}}_s = (0,0,\ldots,0,1), \end{displaymath} (3)

see [11,14,15,21,22,24].

In practice usually non-overlapping $s$-tuples

\begin{displaymath} {\vec{x}}_i := (x_{is}, x_{is+1}, \ldots , x_{is+s-1}), \quad i\ge 0, \end{displaymath} (4)

are used to produce ``independent'' random points in $I^s$. If $gcd(s,m-1) = 1$, hence the dimension $s$ does not divide the period $m-1$ of the generator, then all possible non-overlapping tuples ${\vec{x}}_i$ originate in the same lattice as above. Note that for the computation of all non-overlapping $s$-tuples with an LCG one has to generate at least $s$ times the period. If $gcd(s,m-1) > 1$, these vectors produce true subsets of the lattice $L_s(a,N)$ which must not have a pure lattice structure in general [1]. A lattice is obtained if one considers the union of all possible subsets produced by such tuples.

Essentially different lattices $L_s(a',m')$ are obtained for vectors (1) with arbitrary lags which occur if, for example, lagged subsequences from the output of an LCG are used [9,18].


Quasi-Monte Carlo (QMC)

MC - and also QMC methods can efficiently be applied for an approximate calculation of high dimensional integrals.

Consider the standard domain $I^s:=[0,1[^s$ in dimension $s\ge 2$, a point (node) set $P = \{ \vec{x}_1,\ldots , \vec{x}_N\}$ in $I^s$, $N\in \mathbb{N}$, and a function $f:I^s\longrightarrow \mathbb{R}$. The (quasi) Monte Carlo approximation of an integral $E(f):=\int_{I^s} f(\vec{x})\;d\vec{x}$ is computed by the average of the integrand over the point set $P$

\begin{displaymath} S_N(f,P):=\frac{1}{N} \sum_{n = 1}^N f(\vec{x}_n). \end{displaymath} (5)

A classical method for quasi Monte Carlo integration uses the intersection of the full lattice $L_s(a,m)$, $m = N$, with the $s$-dimensional unit cube $I^s$ as node set for the calculation of the latter approximation.

Note that the difference of MC and QMC integration in our case lies in the fact that for MC one may use a large modulus $m$ and only a ``random'' part of the lattice $L_s(a,m)$, and for QMC, one may apply the sample size $N$ as modulus and therefore the full lattice $L_s(a,N)$.

The application of $L_s(a,m)$ is called rank-1 lattice rule or Korobov lattice rule [13,19,22,26]. The quality of $L_s(a,m)$ depends on the generating vector $\vec{b}_1$ in (3). For well chosen vectors $\vec{b}_1$, the Korobov lattice rule is also called method of ``good lattice points'' (GLPs). More general lattice rules are for example obtained for different vectors in basis (3).


Lattice Assessment

The coarseness of the lattice $L_s(a,m)$ may change dramatically if either the dimension $s$ or the multiplier $a$ are varied. To get reliable linear random number generators for a large class of applications, it is necessary to assess the quality of the several lattices produced by various tuple constructions described above. Furthermore, the selection of ``good'' lattices $L_s(a,m)$ provides excellent QMC node sets as well.

The spectral test (geometric version) is a classical measure for the quality of $s$-dimensional lattices $L_s$. Specifically, this test determines the maximum distance $d_s$ between adjacent hyper-planes, taken over all families of parallel hyper-planes which contain all points of the lattice. The smaller $d_s$ the more regular is the point structure.

Widely used is a normalized spectral test $S_s := d^*_s / d_s$, $2 \le s \le 8$, for which $0 \le S_s \le 1$ (values near 1 imply a good lattice structure). The constants $d^*_s$ are absolute lower bounds on $d_s$, see [14, p. 105] and [11, Sect. 7.7]. L'Ecuyer [17] used also certain lower bounds $d^*_s$ for dimensions $s > 8$ in order to compute $S_s$ for arbitrary dimensions.

The algorithm to calculate the spectral test is based on the dual lattice of $L_s$, since the maximal distance of adjacent hyper-planes $d_s$ is equal to one over the length of the shortest vector of the dual lattice [6,14]. Historically this test is due to Coveyou and MacPherson [5] who used multivariate Fourier analysis to study the quality of LCGs. An efficient implementation of the spectral test for arbitrary multiple recursive generators is given in [18].


Spectral Test Approximation with the LLL Algorithm

In this section we want to demonstrate the effects which appear if the shortest vector in the spectral test calculation is replaced by an approximation obtained by the Lenstra Lenstra Lovász (LLL) basis reduction algorithm [4,20,23]. The calculation of the shortest vectors in a lattice is performed by variants of the Fincke-Pohst algorithm which are in general not polynomial in dimension [10]. Applying the LLL algorithm, an approximation of the shortest vector can be calculated in polynomial time [4,23,27]. For recent discussions on the complexity of lattice problems see [2,3], and also other papers from [7].

For a given basis $B = \{\vec{b}_1,\ldots ,\vec{b}_s\}$ of a lattice $L_s$, the LLL algorithm finds a new basis $B' = \{\vec{b}'_1,\ldots ,\vec{b}'_s\}$, with

\begin{displaymath} \vert\vert\vec{b}'_1\vert\vert\le 2^{\frac{s-1}{2}}\cdot \vert\vert\vec{v}\vert\vert, \end{displaymath} (6)

where $\vec{b}'_1$ denotes the shortest vector in $B'$, $\vert\vert.\vert\vert$ the euclidian norm, and $\vec{v}$ an arbitrary non-zero vector in the lattice $L_s$. Therefore the latter inequality also holds for a shortest vector $\vec{v}_1$ in $L_s$. Further properties of LLL reduced bases are given in the book of Cohen [4]. The latter author quotes: ``We see that the vector $\vec{b}'_1$ in a reduced basis is, in a very precise sense, not too far from being the shortest non-zero vector of $L_s$. In fact, it often is the shortest, and when it is not, one can, most of the time, work with $\vec{b}'_1$ instead of the actual shortest vector''. Moreover, Pohst [23] supplements: ``Examples show that there exist lattices with LLL reduced bases $\{\vec{b}'_1, \ldots , \vec{b'}_s\}$ with $\vert\vert\vec{b}'_1\vert\vert^2 \ge (4/3)^{s-2} \vert\vert\vec{v}_1\vert\vert^2$. These observations certainly do not favour applications of LLL reduced bases. However, the results in practice are in general much better than the worst case estimates.''

In the following we demonstrate that the above statements on the good behavior in practice apply also for the spectral test for parameter selection of MC and QMC node sets. Therefore, in most cases, it is sufficient to apply the LLL reduction instead of the calculation of the shortest vector.

Consider our lattices $L_s:=L_s(a,p)$, $p$ prime, as defined in Section 1. The LLL spectral test $d'_s$ of $L_s$ will obviously be defined as one over the length of the first vector of the LLL reduced dual lattice basis of $L_s$. Similar as in Subsection 1.3 we may also use a normalized LLL spectral test $S'_s := d^*_s / d'_s$, $s\ge 2$. Further we will use the notation $M'_j:=\min_{2\le s \le j} S'_s$ which in the original form $M_j:=\min_{2\le s \le j} S_s$ has often been applied for LCG parameter searches [11,17]. For fixed prime numbers $p$ we will write $d'_s(a,p)$, $S'_s(a,p)$ and $M'_j(a,p)$ for the corresponding figures of merit for lattice $L_s(a,p)$.

To exhibit different behavior of the spectral test and its LLL version for our lattices we considered the nearest prime numbers $p_j$ to $2^j$, $9 \le j \le 28$ which are for example $\{2^9 -3, 2^{10}-3, 2^{11}+5, \ldots \}$ and, for each of these primes, the set of primitive roots $a$ modulo $p_j$. Note that assessments of lattices $L_s(a,p)$ for MC node set selection are in general restricted to the set $A$ of all primitive roots $a$ modulo $p$ since for such primitve roots the corresponding LCG guarantees maximal period. There are $\phi (p-1)$ primitive roots where $\phi$ denotes the Euler totient function. Further there are certain lattices $L_s(a,p)$ for $a \in A$ which are equivalent with respect to the spectral test and therefore it suffices to assess $L_s(a,p)$ for a number of $\phi (p-1)/2$ primitive roots $a$ [11,17].

Figure 1 shows relative frequences of the occurrence of different values of $M_8(a,p)$ and $M'_8(a,p)$ (left graphics) and $d_s(a,p)$ and $d'_s(a,p)$, $2 \le s \le 8$ (right graphics). For prime numbers $p_j$, $j \le 18$ we considered all relevant primitive roots modulo $p_j$, and for $p_j$ with $19 \le j \le 28$, we have randomly chosen a set of $2^{15}$ primitive roots $a$. The relative number of different values $d_s(a,p)$ and $d'_s(a,p)$, $2 \le s \le 8$ is very low (maximum about three percent) but increasing with the dimension. However for the measures $M_8$ and $M'_8$ the frequences are significantly lower. We have chosen $M_8$ since the normalization constants $d^*_s$ used for the normalized spectral test are best possible for dimensions $2 \le s \le 8$, which is not the case for larger dimensions [17].

In Figure 2 we exhibit the magnitude of the differences of the measures. The left graphics shows the maximal values $max_a:=\max_a ( d_s(a,p) / d'_s(a,p) )^2$ for each prime number $p_j$, $12\le j \le 28$ and dimension $2 \le s \le 8$. Note that almost all values are between one and two. There are some outliers, the largest one equals 7 which for dimension $s=7$ is still clearly lower than the theoritical bound $2^{(s-1)}$ given in (6). For dimension $s=2$ all values are equal to one which means that in all cases in dimension two the LLL reduction already provided the shortest vector i.e. $d_s = d'_s$. The latter property can also be seen from the right graphics in Figure 1 and 2. The right graphics shows the mean values of the absolute differences between $S_s(a,p)$ and $S'_s(a,p)$.

Figure 1: Relative frequences of the occurrence of different values of $M_8(a,p)$ and $M'_8(a,p)$ (left graphics) and $d_s(a,p)$ and $d'_s(a,p)$, $2 \le s \le 8$ (right graphics).
\begin{figure}\par\begin{center} \begin{tabular}{cc} \epsfxsize =6.5cm\epsfbox{e... ...psfxsize =6.5cm\epsfbox{eps/freqOfDiff.eps}\end{tabular}\end{center}\end{figure}

Figure 2: The maximal values $max_a$ for each prime number $p_j$, $12\le j \le 28$ and dimension $2 \le s \le 8$ (left graphics) and the mean values of the absolute differences of $S_s(a)$ and $S'_s(a)$ (right graphics).
\begin{figure}\par\begin{center} \begin{tabular}{cc} \epsfxsize =6.5cm\epsfbox{e... ...size =6.5cm\epsfbox{eps/meanOfDiffNorm.eps}\end{tabular}\end{center}\end{figure}

From these empirical tests we can confirm the statements of Cohen and Pohst above, i.e. for our applications the vector $\vec{b}'_1$ of the LLL reduced dual basis of $L_s(a,b)$ is very often the shortest vector and if not the results are much better than the theoretical bounds. Our comparisons have been calculated using a Mathematica implementation of the Fincke-Pohst algorithm by Wilberd van der Kallen, University of Utrecht NL, which is based on a previous LLL reduction. Figure 3 exhibits the time performance of our calculations. The figure displays the mean values of the time used for the calculation of $d_s(a,p)$ divided by the means for $d'_s(a,p)$, the means are taken over 64 different primitive roots $a$ for each prime $p$. The LLL reduction is for small dimensions about a factor 5 faster than the calculation of the shortest vector, but the speedup in our Mathematica implementation decreases for increasing dimension and prime size.

Figure 3: Mean values of the time spent for $d_s(a,p)$ divided by the means for $d'_s(a,p)$, the means taken over 64 different primitive roots $a$ for each prime $p$.
\begin{figure}\begin{center} \begin{tabular}{c} \epsfxsize =7.5cm\epsfbox{eps/meanTimeFac.eps}\end{tabular}\end{center}\end{figure}

We further used Victor Shoup's C++ implementation of the LLL algorithm [25] and performed the same parameter searches for optimal multipliers with respect to $M_8$ which were carried out in [17, Table 2]. For the exhaustive searches in the latter paper which were computed for prime numbers $p \in \{2^8-5, 2^9-3, 2^{10}-3, \ldots , 2^{26}-5\}$ we got exactly the same multipliers with respect to $M'_8$ and equal values for the measures $M'_8$ and $M_8$. For the random searches for primes $p \in \{2^{27}-39, 2^{28}-57, \ldots , 2^{37}-25\}$, for example, we easily obtained improved results with our algorithm. Table 1 shows the results given in [17, Table 2] and our multipliers selected with LLL for the latter primes.


Table 1: Multipliers selected via the spectral test from [17] and our results obtained by a random search using LLL reduction.
$p$ $a$ [17] $M_8$ [17] $a$ $M'_8=M_8$
$2^{27}-39 = 134217689$ 45576512 0.75874 45576512 0.75874
$2^{28}-57 = 268435399$ 150873839 0.74215 150873839 0.74215
$2^{29}-3 = 536870909$ 520332806 0.75238 435136037 0.75356

$2^{30}-35 = 1073741789$

 771645345 0.74881 325079677 0.75432

$2^{31}-1 = 2147483647$

 1583458089 0.72771 598753959 0.73435
    117879879 0.74309
    629824009 0.74880
    1355089539 0.74972

$2^{32} -5 = 4294967291$

 1588635695 0.74530 3265168268 0.74870

$2^{33}-9 = 8589934583 $

 7425194315 0.73666 8137022074 0.75316

$2^{34}-41 = 17179869143$

 5295517759 0.73607 10771374442 0.73899
    1491142424 0.75157

$2^{35}-31 = 34359738337$

 3124199165 0.74740 23314821278 0.75022

$2^{36}-5 = 68719476731$

 49865143810 0.72011 24365995562 0.75969
    46865245638 0.76825

$2^{37}-25 = 137438953447$

 76886758244 0.73284 64192466011 0.73997
       

Bibliography

1
L. Afflerbach.
The sub-lattice structure of linear congruential random number generators.
Manuscripta Mathematica, 55:455-465, 1986.

2
J Ajtai.
Generating Hard Instances of Lattice Problems.
Report TR96-007, Electronic Colloquium on Computational Complexity ECCC, 1996.
Web-page: [7].

3
J Cai.
Some Recent Progress on the Complexity of Lattice Problems.
Report TR99-006, Electronic Colloquium on Computational Complexity ECCC, 1999.
Web-page: [7].

4
H. Cohen.
A Course in Computational Algebraic Number Theory, volume 138 of Graduate Texts in Mathematics.
Springer-Verlag, 1993.

5
R.R. Coveyou and R.D. MacPherson.
Fourier analysis of uniform random number generators.
J. Assoc. Comput. Mach., 14:100-119, 1967.

6
U. Dieter.
How to Calculate Shortest Vectors in a Lattice.
Math. Comp., 29(131):827-833, 1975.

7
ECCC.
Electronic Colloquium on Computational Complexity ECCC.
http://www.eccc.uni-trier.de/eccc/.

8
K. Entacher.
A collection of selected pseudorandom number generators with linear structures - advanced version.
Technical report, Dept. of Mathematics, University Salzburg, Austria, available at: http://random.mat.sbg.ac.at/, 1998.
The previous version is published as technical report 97-1, ACPC-Austrian Center for Parallel Computation, University of Vienna, Austria, 1997.

9
K. Entacher.
Parallel Streams of Linear Random Numbers in the Spectral Test.
ACM Transactions on Modeling and Computer Simulation, 9(1):31-44, 1999.

10
U. Fincke and M. Pohst.
Improved methods for calculating vectors of short length in a lattice, including a complexity analysis.
Math. Comp., 44:463-471, 1985.

11
G.S. Fishman.
Monte Carlo: Concepts, Algorithms, and Applications, volume 1 of Springer Series in Operations Research.
Springer, New York, 1996.

12
P. Hellekalek and G. Larcher (eds.).
Random and Quasi-Random Point Sets, volume 138 of Lecture Notes in Statistics.
Springer, Berlin, 1998.

13
F.J. Hickernell.
Lattice Rules: How Well Do They Measure Up.
In [12] pp 109-166.

14
D.E. Knuth.
The Art of Computer Programming, volume 2: Seminumerical Algorithms.
Addison-Wesley, Reading, MA, 2nd edition, 1981.

15
P. L'Ecuyer.
Uniform random number generation.
Ann. Oper. Res., 53:77-120, 1994.

16
P. L'Ecuyer.
Random Number Generation.
In Jerry Banks, editor, Handbook of Simulation, Chapter 4. Wiley, 1998.

17
P. L'Ecuyer.
Tables of Linear Congruential Generators of Different Sizes and Good Lattice Structure.
Mathematics of Computation, 68(225):249-260, 1999.

18
P. L'Ecuyer and R. Couture.
An Implementation of the Lattice and Spectral Tests for Multiple Recursive Linear Random Number Generators.
INFORMS Journal on Computing, 9(2):209-217, 1997.

19
P. L'Ecuyer and C. Lemieux.
Variance Reduction via Lattice Rules.
Management Science, to appear, 2000.

20
A.K. Lenstra, H.W. Lenstra, and L. Lovász.
Factoring polynomials with rational coefficients.
Mathematische Annalen, 261:515-534, 1982.

21
G. Marsaglia.
The structure of linear congruential sequences.
In S. K. Zaremba, editor, Applications of Number Theory to Numerical Analysis, pages 248-285. Academic Press, New York, 1972.

22
H. Niederreiter.
Random Number Generation and Quasi-Monte Carlo Methods.
SIAM, Philadelphia, 1992.

23
M.E. Pohst.
Computational Algebraic Number Theory.
DMV Seminar Band 21. Birkhäuser Verlag, 1993.

24
B.D. Ripley.
The lattice structure of pseudo-random number generators.
Proc. Roy. Soc. London Ser. A, 389:197-204, 1983.

25
V. Shoup.
NTL: A Library for doing Number Theory.
http://www.shoup.net/.

26
I.H. Sloan and S. Joe.
Lattice Methods for Multiple Integration.
Oxford Univ. Press, New York, 1994.

27
A. Storjohann.
Faster Algorithms for Integer Lattice Basis Reduction.
Technical report, Institute of Scientific Computing, ETH Zürich, 1996.
Available at http://www.inf.ethz.ch/research/wr/.