projet_optim_m2/main.tex

% ------------------------------------------------------------------------
\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
  simmulated annealing and genetic $(5+5)$ search with a ad-hoc 
  crossover function. 
\end{abstract}


\section{Introduction}

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

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 corresponging 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 1600 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}

%# TODO sa n'est pas evolutionnaire
\subsection{Evolutionary heuristic: Simulated annealing and local search}
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 simultations), 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 simulated annealing local search 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.

\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 (5+5) search}


\begin{figure}
 \begin{mdframed}
\includegraphics[scale=0.4]{crossover_flow.pdf}
\caption{Flowchart of the crossover algorithm.}
  \label{cross_flow}
\end{mdframed}
\end{figure}

\subsubsection{Dependence on the parameter p}
First experiments were made to select the value for the crossover parameter
p. 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.

\begin{figure}
\includegraphics[scale=0.3]{Results/res_gen_2.png}
\caption{Dependence on parameter p: D=2}
  \label{res_gen2}
\end{figure}

\begin{figure}
\includegraphics[scale=0.3]{Results/res_gen_2_zoom.png}
\caption{Dependence on parameter p (zoom): D=2}
\label{res_gen2z}
\end{figure}
\begin{figure}
\includegraphics[scale=0.3]{Results/res_gen_3_zoom.png}
\caption{Dependence on parameter p: D=3}
  \label{res_gen3}
\end{figure}

\begin{figure}
\includegraphics[scale=0.3]{Results/res_gen_4_zoom.png}
\caption{Dependence on parameter p: 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.
An interesting phenomena can be observed: the error rates are somehow 
invariant w.r.t.\ the number of iteration 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.


\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*{Acknoledgments}
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}