123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502 |
- % ------------------------------------------------------------------------
- \documentclass{llncs}
- \input{prelude}
- \begin{document}
- \title{$\mbox{\EightStarBold}$ Discrepancies for generalized Halton points\\
- Comparison of three heuristics for generating points set}
- % \titlerunning{} % abbreviated title (for running head)
- % also used for the TOC unless
- % \toctitle is used
- \author{Thomas Espitau$\mbox{}^{\mbox{\SnowflakeChevron}}$ \and Olivier Marty$\mbox{}^{\mbox{\SnowflakeChevron}}$}
- %
- % \authorrunning{} % abbreviated author list (for running head)
- %
- %%%% list of authors for the TOC (use if author list has to be modified)
- % \tocauthor{}
- %
- \institute{$\mbox{}^{\mbox{\SnowflakeChevron}}$ ENS Cachan}
- \maketitle
- \makeatletter
- \renewcommand\bibsection%
- {
- \section*{\refname
- \@mkboth{\MakeUppercase{\refname}}{\MakeUppercase{\refname}}}
- }
- \makeatother
- \begin{abstract}
- Geometric discrepancies are standard measures to quantify the irregularity of
- distributions. They are an important notion in numerical integration.
- One of the most important discrepancy notions is the so-called star
- discrepancy. Roughly speaking, a point set of low star discrepancy value
- allows for a small approximation error in quasi-Monte Carlo integration.
- In this work we present a tool realizing the implantation of three
- basics heuristics for construction low discrepancy points sets
- in the generalized Halton model: fully random search, local search with
- simulated annealing and genetic $(5+5)$ search with a ad-hoc
- crossover function.
- \end{abstract}
- \section{General architecture of the tool}
- The testing tool is aimed to be modular: it is made of independents blocks that
- are interfaced trough a scheduler. More precisely a master wrapper is written
- in Python that calls a first layer which performs the chosen heuristic. This
- layer is written in C++ for performances. The given discrepancy algorithm
- --- written in C --- is called when evaluations of a state is needed.
- The wrapper dispatch the computations on the multi-core architecture of
- modern computers\footnote{for us, between 2 and 4 physical cores and 4 or 8
- virtual cores}. This basic architecture is described in figure~\ref{main_flow}.
- More precisely the class diagram for the unitary test layer is
- presented in figure~\ref{class_flow}.
- Experiments were conducted on two machines:
- \begin{itemize}
- \item 2.4 GHz Intel Dual Core i5 hyper-threaded to 2.8GHz, 8 Go 1600 MHz DDR3.
- \item 2.8 GHz Intel Quad Core i7 hyper-threaded to 3.1GHz, 8 Go 1600 MHz DDR3.
- \end{itemize}
- \begin{figure}
- \includegraphics[scale=0.6]{main_flow.pdf}
- \caption{Tool overview}
- \label{main_flow}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=1]{graph_classes.pdf}
- \caption{Class dependencies}
- \label{class_flow}
- \end{figure}
- On these machines, some basic profiling has make clear that
- the main bottleneck of the computations is hiding in the \emph{computation
- of the discrepancy}. The chosen algorithm and implantation of this
- cost function is the DEM-algorithm~\cite{Dobkin} of
- \emph{Magnus Wahlstr\o m}~\cite{Magnus}.\medskip
- All the experiments has been conducted on dimension 2,3,4
- --- with a fixed Halton basis 7, 13, 29, 3 ---. Some minor tests have
- been made in order to discuss the dependency of the discrepancy and
- efficiency of the heuristics with regards to the values chosen for the
- prime base. The average results remains roughly identical when taking
- changing these primes and taking them in the range [2, 100]. For such
- a reason we decided to pursue the full computations with a fixed
- basis.
- \subsection{Algorithmic insights}
- To perform an experiment we made up a
- loop above the main algorithm which calls the chosen heuristic multiple
- times in order to smooth up the results and obtain more exploitable datas.
- Then an arithmetic mean of the result is performed on the values. In addition
- extremal values are also given in order to construct error bands graphs.
- A flowchart of the conduct of one experiment is described in the
- flowchart~\ref{insight_flow}. The number of iteration of the heuristic is
- I and the number of full restart is N. The function Heuristic() correspond to
- a single step of the chosen heuristic. We now present an in-depth view of
- the implemented heuristics.
- \begin{figure}
- \begin{mdframed}
- \includegraphics[scale=0.4]{insight.pdf}
- \caption{Flowchart of a single experiment}
- \label{insight_flow}
- \end{mdframed}
- \end{figure}
- Graph are presented not with the usual box plot to show the
- error bounds, but in a more graphical way with error bands. The graph
- of the mean result is included inside a band of the same color which
- represents the incertitude with regards to the values obtained.
- \section{Heuristics developed}
- \subsection{Fully random search (Test case)}
- The first heuristic implemented is the random search. We generate
- random permutations, compute the corresponding sets of Halton points
- and select the best set with regard to its discrepancy.
- The process is wrapped up in the
- flowchart~\ref{random_flow}. In order to generate at each step random
- permutations, we transform them directly from the previous ones.
- More precisely the permutation is a singleton object which have method
- random, built on the Knuth Fisher Yates shuffle. This algorithm allows
- us to generate an uniformly chosen permutation at each step. We recall
- this fact and detail the algorithm in the following section.
- \begin{figure}
- \begin{mdframed}
- \includegraphics[scale=0.4]{flow_rand.pdf}
- \caption{Flowchart of the random search}
- \label{random_flow}
- \end{mdframed}
- \end{figure}
- \subsubsection{The Knuth-Fisher-Yates shuffle}
- The Fisher–Yates shuffle is an algorithm for generating a random permutation
- of a finite sets. The Fisher–Yates shuffle is unbiased, so that every
- permutation is equally likely. We present here the Durstenfeld variant of
- the algorithm, presented by Knuth in \emph{The Art of Computer programming}
- vol. 2~\cite{Knuth}.
- The algorithm's time complexity is here $O(n)$, compared to $O(n^2)$ of
- the naive implementation.
- \begin{algorithm}[H]
- \SetAlgoLined
- \SetKwFunction{Rand}{Rand}
- \SetKwFunction{Swap}{Swap}
- \KwData{A table T[1..n]}
- \KwResult{Same table T, shuffled}
- \For{$i\leftarrow 1$ \KwTo $n-1$}{
- $j \leftarrow$ \Rand{$[1,n-i]$}\;
- \Swap{$T[i], T[i+j]$}\;
- }
- \caption{KFY algorithm}
- \end{algorithm}
- \begin{lemma}
- The resulting permutation of KFY is unbiased.
- \end{lemma}
- \begin{proof}
- Let consider the set $[1,\ldots n]$ as the vertices of a random graph
- constructed as the trace of the execution of the algorithm:
- an edge $(i,j)$ exists in the graph if and only if the swap of $T[i]$ and
- $T[j]$ had been executed. This graph encodes the permutation represented by
- $T$. To be able to encode any permutation the considered graph must be
- connected --- in order to allow any pairs of points to be swapped ---.
- Since by construction every points is reached by an edge, and that there
- exists exactly $n-1$ edges, we can conclude directly that any permutation can
- be reached by the algorithm. Since the probability of getting a fixed graph
- of $n-1$ edges with every edges of degree at least one is $n!^{-1}$, the
- algorithm is thus unbiased.
- \end{proof}
- \subsubsection{Results and stability}
- We first want to analyze the dependence of the results on the number of
- iterations of the heuristic, in order to discuss its stability.
- The results are compiled in the figures~\ref{rand_iter2},~\ref{rand_iter3},
- restricted to a number of points between 80 and 180.
- We emphasize on the fact the lots of datas appears on the graphs,
- and the error bands representation make them a bit messy. These graphs
- were made for extensive internal experiments and parameters researches.
- The final wrap up graphs are much more lighter and only presents the best
- results obtained.
- As expected from a fully random search, the error bands are very large for
- low number of iterations ($15\%$ of the value for 400 iterations) and tends
- to shrink with a bigger number of iterations (around $5\%$ for 1500 iterations).
- This shrinkage is a direct consequence of well known concentrations bounds
- (Chernoff and Asuma-Hoeffding).
- The average results are quite stable, they decrease progressively with
- the growing number of iterations, but seem to get to a limits after 1000
- iterations. This value acts as a threshold for the interesting number of iterations.
- As such interesting results can be conducted with \emph{only} 1000 iterations,
- without altering too much the quality of the set with regards to its
- discrepancy and this heuristic.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/random_iter.png}
- \caption{Dependence on iterations, dimension 2}
- \label{rand_iter2}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/random_iter_3.png}
- \caption{Dependence on iterations, dimension 3}
- \label{rand_iter3}
- \end{figure}
- \subsection{Local search with simulated annealing}
- The second heuristic implemented is a randomized local search with
- simulated annealing. This heuristic is inspired by the physical
- process of annealing in metallurgy.
- Simulated annealing interprets the physical slow cooling as a
- slow decrease in the probability of accepting worse solutions as it
- explores the solution space.
- More precisely a state is a $d$-tuple of permutations, one for each dimension,
- and the neighbourhood is the set of $d$-tuple of permutations which can be obtained
- by application of exactly one transposition of one of the permutation of
- the current state.
- The selection phase is dependant on the current temperature:
- after applying a random transposition on one of the current permutations, either
- the discrepancy of the corresponding Halton set is decreased and the
- evolution is kept, either it does not but is still kept with
- a probability $e^{\frac{\delta}{T}}$ where $\delta$ is the difference
- between the old and new discrepancy, and $T$ the current temperature.
- If the de discrepancy has decreased, the temperature $T$ is multiplied
- by a factor $\lambda$ (fixed to $0.992$ in all our simulations), hence
- is decreased.
- The whole algorithm is described in the flowchart~\ref{flow_rec}.
- \begin{figure}
- \begin{mdframed}
- \includegraphics[scale=0.4]{flow_recuit.pdf}
- \caption{Flowchart of the local search with simulated annealing heuristic}
- \label{flow_rec}
- \end{mdframed}
- \end{figure}
- \subsubsection{Dependence on the temperature}
- First experiments were made to select the best initial temperature.
- Results are compiled in graphs~\ref{temp_2},~\ref{temp3},\ref{temp3_z}.
- Graphs~\ref{temp_2},~\ref{temp3} represent the results obtained respectively
- in dimension 2 and 3 between 10 and 500 points. The curve obtained is
- characteristic of the average evolution of the discrepancy optimization
- algorithms for Halton points sets: a very fast decrease for low number of
- points --- roughly up to 80 points --- and then a very slow one
- after~\cite{Doerr}.
- The most interesting part of these results are concentrated between 80 and 160
- points were the different curves splits. The graph~\ref{temp3_z} is a zoom
- of~\ref{temp3} in this window. We remark on that graph that the lower the
- temperature is, the best the results are, with a threshold at $10^{-3}$.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/resu_2_temp.png}
- \caption{Dependence on initial temperature: D=2}
- \label{temp_2}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/resu_temp3.png}
- \caption{Dependence on initial temperature: D=3}
- \label{temp3}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/resu_temp3_zoom.png}
- \caption{Dependence on initial temperature (zoom): D=3}
- \label{temp3_z}
- \end{figure}
- \subsubsection{Stability with regards to the number of iterations}
- As for the fully random search heuristic we investigated the stability
- of the algorithm with regards to the number of iterations. We present here
- the result in dimension 3 in the graph~\ref{iter_sa}. Once again we
- restricted the window between 80 and 180 points were curves are split.
- An interesting phenomena can be observed: the error rates are somehow
- invariant w.r.t.\ the number of iterations and once again the 1000 iterations
- threshold seems to appear --- point 145 is a light split between iteration
- 1600 and the others, but excepted for that point, getting more than 1000
- iterations tends be be a waste of time. The error rate is for 80 points the
- biggest and is about $15\%$ of the value, which is similar to the error
- rates for fully random search with 400 iterations.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/sa_iter.png}
- \caption{Dependence on iterations number for simulated annealing : D=3}
- \label{iter_sa}
- \end{figure}
- \subsection{Genetic $(\mu+\lambda)$ search}
- The third heuristic implemented is the $(\mu+\lambda)$ genetic search.
- This heuristic is inspired from the evolution of species: a family
- of $\mu$ genes is known (they are generated randomly at the beginning),
- from which $\lambda$ new genes are derived. A gene is the set of parameters
- we are optimizing, i.e. the permutations.
- Each one is derived either from one gene applying a mutation
- (here a transposition of one of the permutations), or from two
- genes applying a crossover : a blending of both genes (the
- algorithm is described in details further). The probability of
- making a mutation is $c$, the third parameter of the algorithm,
- among $\mu$ and $\lambda$. After that, only the $\mu$ best
- genes are kept, according to their fitness, and the process
- can start again.
- Because making variations over $\mu$ or $\lambda$ does not change fundamentally
- the algorithm, we have chosen to fix $\mu=\lambda=5$ once and for all,
- which seemed to be a good trade-off between the running time of
- each iteration and the size of the family.
- \subsubsection{Crossover algorithm}
- We designed a crossover for permutations. The idea is simple: given two
- permutations $A$ and $B$ of $\{1..n\}$, it constructs a new permutations
- $C$ value after value, in a random order (we use our class permutation
- for it). For each $i$, we take either $A_i$ or $B_i$. If exactly
- one of those values is available (understand it was not already chosen)
- we choose it. If both are available, we choose randomly and we remember
- the
- second. If both are unavailable, we choose a remembered value. The
- flowchart of this algorithm is described in figure~\ref{cross_flow}.
- \begin{figure}
- \begin{mdframed}
- \includegraphics[scale=0.4]{crossover_flow.pdf}
- \caption{Flowchart of the crossover algorithm.}
- \label{cross_flow}
- \end{mdframed}
- \end{figure}
- The benefits of this method are that it keeps common values of $A$ and $B$,
- the values $C_i$ are often among $\{A_i, B_i\}$ (hence $C$ is close to $A$
- and $B$), and it does not favor either the beginning or the ending of
- permutations.
- \begin{algorithm}[H]
- \SetAlgoLined
- \SetKwFunction{Rand}{Rand}
- \SetKwFunction{RPO}{a random permutation of }
- \SetKwFunction{RVI}{a random value in }
- \KwData{Two permutations A[1..n], B[1..n]}
- \KwResult{A permutation C[1..n]}
- $pi \leftarrow \RPO \{1,\dots,n\}$\;
- $available \leftarrow \{\}$\;
- $got \leftarrow \{\}$\;
- \For{$i\leftarrow 1$ \KwTo $n$}{
- $j \leftarrow pi_i$\;
- $a \leftarrow A_j$\;
- $b \leftarrow B_j$\;
- \If{$a \in got \land b \in got$}{
- $v \leftarrow \RVI available$\;
- }
- \ElseIf{$a \in got$}{
- $v \leftarrow b$\;
- }
- \ElseIf{$b \in got$}{
- $v \leftarrow a$\;
- }
- \Else{
- \If{$\Rand(\{0,1\}) = 1$}{
- $v \leftarrow a$\;
- $available \leftarrow available \cup \{b\}$\;
- }
- \Else{
- $v \leftarrow b$\;
- $available \leftarrow available \cup \{a\}$\;
- }
- }
- $C_j \leftarrow v$\;
- $got \leftarrow got \cup \{v\}$\;
- $available \leftarrow available \setminus \{v\}$\;
- }
- \caption{Permutations crossover}
- \end{algorithm}
- \subsubsection{Dependence on the parameter $c$}
- First experiments were made to select the value for the crossover parameter
- $c$. Results are compiled in graphs~\ref{res_gen2},~\ref{res_gen2z},\ref{res_gen3}
- and~\ref{res_gen4}.
- Graph~\ref{res_gen2}, represents the results obtained
- in dimension 2 between 10 and 500 points. The curve obtained is, with no
- surprise again,
- the characteristic curve of the average evolution of the discrepancy we already
- saw with the previous experiments.
- The most interesting part of these results are concentrated --- once again ---
- between 80 and 160 points were the different curves splits.
- The graph~\ref{res_gen2z} is a zoom of~\ref{res_gen2} in this window, and
- graphs~\ref{res_gen3} and~\ref{res_gen4} are focused directly into it too.
- We remark that in dimension 2, the results are better for $c$ close to $0.5$
- whereas for dimension 3 and 4 the best results are obtained for $c$ closer to
- $0.1$.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/res_gen_2.png}
- \caption{Dependence on parameter $c$: D=2}
- \label{res_gen2}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/res_gen_2_zoom.png}
- \caption{Dependence on parameter $c$ (zoom): D=2}
- \label{res_gen2z}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/res_gen_3_zoom.png}
- \caption{Dependence on parameter $c$: D=3}
- \label{res_gen3}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/res_gen_4_zoom.png}
- \caption{Dependence on parameter $c$: D=4}
- \label{res_gen4}
- \end{figure}
- Once again we investigated the stability
- of the algorithm with regards to the number of iterations. Once again we
- restricted the window between 80 and 180 points were curves are split.
- The results are compiled in graph~\ref{gen_iter}.
- An interesting phenomena can be observed: the error rates are getting really
- big for 1400 iterations at very low points (up to 120),
- even if the average results are stables after the threshold 1000 iterations,
- like we get before.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/gen_iter.png}
- \caption{Stability w.r.t.\ number of iterations: D=2}
- \label{gen_iter}
- \end{figure}
- \section{Results}
- Eventually we made extensive experiments to compare the three previously
- presented heuristics. The parameters chosen for the heuristics have been
- chosen using the experiments conducted in the previous sections
- Results are compiled in the last
- figures~\ref{wrap2},~\ref{wrap2z},~\ref{wrap3z},~\ref{wrap4z}. The
- recognizable curve of decrease
- of the discrepancy is still clearly recognizable in the graph~\ref{wrap2},
- made for points ranged between 10 and 600. We then present the result
- in the --- now classic --- window 80 points - 180 points ---.
- For all dimensions, the superiority of non-trivial algorithms --- simulated
- annealing and genetic search --- is clear over fully random search.
- Both curves for these heuristics are way below the error band of random
- search. As a result \emph{worse average results of non trivial heuristics are
- better than best average results when sampling points at random}.
- In dimension 2~\ref{wrap2z}, the best results are given by the simulated annealing,
- whereas in dimension 3 and 4~\ref{wrap3z},~\ref{wrap4z}, best results are
- given by genetic search. It is also noticeable that in that range
- of points the error rates are roughly the same for all heuristics:
- \emph{for 1000 iteration, the stability of the results is globally the
- same for each heuristic}.
- \begin{figure}
- \includegraphics[scale=0.3]{Results/wrap_2.png}
- \caption{Comparison of all heuristics: D=2}
- \label{wrap2}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/wrap_2_zoom.png}
- \caption{Comparison of all heuristics (zoom): D=2}
- \label{wrap2z}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/wrap_3.png}
- \caption{Comparison of all heuristics: D=3}
- \label{wrap3z}
- \end{figure}
- \begin{figure}
- \includegraphics[scale=0.3]{Results/wrap_4.png}
- \caption{Comparison of all heuristics: D=4}
- \label{wrap4z}
- \end{figure}
- \section{Conclusion}
- \section*{Acknowledgments}
- We would like to thank Magnus Wahlstrom from the Max Planck Institute for Informatics
- for providing an implementation of the DEM algorithm [DEM96].
- We would also like to thank Christoff Durr and Carola Doerr
- for several very helpful talks on the topic of this work.
- Both Thomas Espitau and Olivier Marty supported by the French Ministry for
- Research and Higher Education, trough the Ecole Normale Supérieure.
- \bibliographystyle{alpha}
- \bibliography{bi}
- \end{document}
|