\lecturetitle{16}{N-body Simulations 2}{22 Oct 97}{Steve Gribble}
\begin{itemize} 
  \item N-body continued: Callahan \& Kosaraju Algorithm

  \item Fast Multipole Method (Review)
    \begin{enumerate}
       \item The Four Ideas
       \item Interaction Lists
    \end{enumerate}
  \item Callahan and Kosaraju's Algorithm
    \begin{enumerate}
       \item Well Separated Pair Decompositions
       \item Definitions
       \item Algorithm: Tree Decomposition
       \item Algorithm: Well-Separated Realization
       \item Bounding the size of the realization
    \end{enumerate}
   \item N-body Summary
\end{itemize}

%%%% Body of notes

\section{Fast Multipole Method (Review)}

\subsection{The Four Ideas}

The four main ideas being the fast multipole method (FMM) were:
\begin{itemize}
\item {\bf Multipole expansion} - a Taylor-series like expansion 
      of the potential due to a cluster of particles around an origin.
      The expansion converges at points outside the cluster
      so that the potential can be approximated by only keeping the 
      first $p$ terms.
      We call this the outer or multipole expansion.
\item {\bf Error bounds on dropping terms} - the number of terms in the
      expansion dictates the accuracy of the expansion
\item {\bf Dual local expansion} - an expansion of the potential due to
      particles outside a cluster which converges for points inside 
      the cluster.
      We call this the inner or local expansion.
\item {\bf Translation of expansions} - being able to translate a inner
      or outer expansion to a new origin, or translate an outer to
      an inner expansion at a different origin.
\end{itemize}

\begin{figure}[htbp]
\centerline{
  \psfig{file=../lecture16/tree.eps,width=300pt}
}
\caption{Fifth idea - using a tree for translations}\label{treetrans}
\end{figure}

A fifth idea (Figure~\ref{treetrans}) is to use a tree-structure to
organize the particles and to determine a strategy for grouping the
particles and performing the outer/outer, outer/inner, and inner/inner
translations.  We call each node of the tree a {\em cell}.

\subsection{Interaction Lists}

In the fast multipole method, the interaction list
(Figure~\ref{interaction}) dictates the particles that a given cell
will interact with (cross edges in Figure~\ref{treetrans}).  The {\em
interaction list} for a cell $c$ consists of all the cells at the same
level of the tree which (1) are well-separated from $c$ and (2) are
children of cells that are not well-separated from the parent of $c$.
In Figure~\ref{interaction}, for example, we assume a cell is well
separated if it is 2 cells away.  Clearly all the outer shaded cells
are 2 away from the center cell $c$ and are therefore well-separated.
Also the parents of these cells are not well-separated from the parent
of $c$ (note that parent cells are marked with the thicker lines).
One issue is how to find particles in the interaction list.  Here
are two possible methods:

\begin{enumerate}

\item {\bf Uniform distributions:} a hierarchy of meshes can be used in
this situation.  Simply look at neighbors at an appropriate resolution in
the mesh hierarchy.

\item {\bf Non-uniform distributions:} a mesh hierarchy may be very
sparsely populated, so a Barnes-Hut tree-like structure may be more
appropriate in this case.

\end{enumerate}

\begin{figure}[htbp]
\centerline{
  \psfig{file=../lecture16/interaction.eps,width=250pt}
}
\caption{The Interaction List}\label{interaction}
\end{figure}

\section{Callahan and Kosaraju's Algorithm}
\index{Callahan and Kosaraju}

We now discuss an algorithm due to Callahan and Kosaraju for
generating a tree and interactions-lists for non-uniform distributions
(``A Decomposition of Multi-Dimensional Point-Sets with Applications
to $k$-Nearest-Neighbors and $n$-Body Potential Fields.'', JACM, 42(1),
January 1995).  The algorithm has the nice property that it is very simple,
although the proof of its bounds is a little involved.  The particular
tree-structure it creates is called a well-separated pair
decomposition.  It uses a somewhat more general definition of
interaction lists than used in the previous section (for example,
it allows interactions between cells on different levels).  The
algorithm works in any dimension, requires at most $O(n \log n)$ time
to build the tree for $n$ particles and guarantees that there will be
a total of $O(n)$ interactions (size of the interaction lists summed
over all cells).  

\subsection{Well Separated Pair Decompositions}

The tree structure used by the algorithm is based on the notion of a
tree decomposition of points.  Given the set of points \{a, b, c, d,
e, f, g\}, a {\em tree decomposition} of these points repeatedly
splits this set into subsets, as in Figure~\ref{treedecomp}.

\begin{figure}[htbp]
\centerline{
  \psfig{file=../lecture16/treedecomp.eps,width=300pt}
}
\caption{A tree decomposition $D$ of points.}
\label{treedecomp}
\end{figure}

A {\em well-separated realization} of a tree decomposition is a set of
interaction edges between vertices (cells) in the tree (edges
represent outer/inner translations) such that (1) the cells at the two
ends of the edge are well-separated (we will define well-separated
later in these notes) and (2) there is exactly one path between any
two leaves that goes up the tree any number of levels (including 0),
across a single interaction-edge, and back down the tree any number of
levels.  The second property guarantees that each particle will
interact exactly once with every other particle.
Figure~\ref{wellsepreal} shows an example of a well-separated
realization for the tree decomposition given in
Figure~\ref{treedecomp}.  You can verify that there is exactly one
path between every leaf.  We will define 
well-separated realizations more formally in 
Section~\ref{sec:definitions}.

\begin{figure}[htbp]
\centerline{
  \psfig{file=../lecture16/wellsepreal.eps,width=400pt}
}
\caption{A well-separated ``realization'' of $D$.}
\label{wellsepreal}
\end{figure}

Together the tree-decomposition and the well-separated realization on
that tree are called a {\em well separated pair decomposition}.
The Callahan-Kosaraju algorithm consists of generating the tree and
the realization in two steps.  We will see that we can build the
tree-decomposition in $O(n \log n)$ time and find the well-separated
realization given the tree in $O(n)$ time.  Furthermore we can bound
the size of the well-separated realization (total number of
interaction edges) by $O(n)$.  

The well-separated pair decomposition has applications both to the
N-body problem and to computational geometry.  For the N-body problem
the tree is used for outer/outer translations (going up the tree) and
inner/inner translations (going down the tree), and the well-separated
realization is used for outer/inner translations.  In computational
geometry the well-separated pair decomposition can be used to find the
$k$-nearest-neighbors for each point.  This requires some augmentation
to the tree and is discussed in the Callahan and Kosaraju paper.

For the N-body problem the time taken by one step of building the
tree, finding the interaction edges, and doing all the translations to
calculate the forces is $k_1 n \log n + k_2 n$, where the first term
is for building the tree and the second linear term is for finding the
interaction edges and translating the expansions.  In practice the
constant $k_2$ is very much larger than $k_1$.  Furthermore in many
cases it is possible to use the same tree (or only slightly modified
tree) for multiple steps.  Therefore from a practical standpoint the
algorithm runs in almost linear time.

\subsection{Definitions}
\label{sec:definitions}

\begin{itemize}

\item {\bf Bounding Rectangle $R(P)$}: the smallest rectangle that
contains points P and is aligned with the axes (Figure~\ref{boundrect}).
\begin{figure}[ht]
\centerline{
  \psfig{file=../lecture16/boundrect.eps,width=175pt}
}
\caption{Bounding Rectangle}\label{boundrect}
\end{figure}

\item 
{\bf $l_{\max}$}: maximum length over the dimensions of a rectangle.  
$l_{\max}(P)$ is the maximum length of $R(P)$.

\item 
{\bf Well-separated}: Given two rectangles, let $r$ be the smallest
radius that contains them both (see Figure~\ref{wellsep}).
For a given {\em separation constant} $s$ two rectangles
are well separated if $d > s r$.
\begin{figure}[ht]
\centerline{
  \psfig{file=../lecture16/wellsep.eps,width=200pt}
}
\caption{Well Separated}\label{wellsep}
\end{figure}

\item {\bf Interaction Product $A \otimes B$ }: 
\[ A \otimes B = \{ \{p, p'\} | p \in A, p' \in B, p \neq p'\}\]

\item {\bf Realization of $A \otimes B$ } is a set $\{\{A_1, B_1\}, \ldots
\{A_k, B_k\}\}$ such that:
\begin{enumerate}
\item $A_i \subseteq A, B_i \subseteq B$ for all $i = 1 \ldots k$
\item $A_i \cap B_i = \emptyset$
\item $ (A_i \otimes B_i) \cap (A_j \otimes B_j) = \emptyset$ for all $i
\neq j$, {\em i.e.}, there are no repeated interactions
\item $A \otimes B = \cup_{i=1}^{k} A_i \otimes B_i $, {\em i.e.}, all
interactions are included.
\end{enumerate}

This formalizes the interaction edges in the tree and guarantees that
there is only one path between two leaves.  In particular rule 2
guarantees that no leaf interacts with itself, rule 3 guarantees that
no leaf interacts more than once with another leaf and rule 4
guarantees that every leaf interacts with every other leaf.

\item {\bf A well separated realization} is a realization such that for all $i$, $A_i$
and $B_i$ are well separated.

\item {\bf A well-separated pair decomposition} is a tree decomposition
plus a well-separated realization $P \otimes P$ where subsets are nodes of
the tree (Figure~\ref{subsetnode}).
\begin{figure}[ht]
\centerline{
  \psfig{file=../lecture16/wellseppair_node.eps,width=250pt}
}
\caption{Subsets are nodes of the tree}\label{subsetnode}
\end{figure}

\end{itemize}

\subsection{Algorithm: Tree Decomposition}

The algorithm partitions space by recursively
cutting bounding rectangles in half along their longest dimension.
It can be defined as follows:
\begin{quote}
\begin{tabbing}
{\bf function} Build\_Tree$(P)$\\
~~~~// $P$ is a set of points\\[.08in]
~~~~{\bf if} $|P| = 1$ {\bf then return} Leaf($P$)\\
~~~~{\bf else}\\
~~~~~~~~$S$ = plane that splits $R(P)$ perpendicular to its longest dimension\\
~~~~~~~~$P_1, P_2$ = Split $P$ by the plane $S$\\
~~~~~~~~{\bf return} Tree\_Node(Build\_Tree($P_1$), Build\_Tree($P_2$), $R(P)$)\\
{\bf end}
\end{tabbing}
\end{quote}
An example of the algorithm on a set of 8 points 
is given in (Figure~\ref{tree_algo}).
\begin{figure}[ht]
\centerline{\psfig{file=../lecture16/tree_algo.eps,width=350pt}}
\caption{Building the tree}\label{tree_algo}
\end{figure}
The algorithm returns a tree with the bounding rectangle stored
at each internal node and the original points at the leaves. 
The bounding rectangles are needed for building the well-separated
realization ({\em e.g.} to determine when nodes are well separated).

For uniformly distributed sets of points the Build\_Tree algorithm
will run in $O(n \log n)$ time without any special tricks since
finding the bounding rectangle and splitting $P$ take $O(n)$ time and
the two recursive calls will be on sets of approximately the same size
({\em i.e.}, we have a standard $T(n) = 2T(n/2) + O(n)$ style
recurrence).  Furthermore the algorithm is highly parallel since the
two recursive calls can be made in parallel and the data can be split
in parallel In fact for most ``practical'' sets of points this will be
true.

If the cut does not evenly divide the points, however, the naive
algorithm could take $O(n^2)$ time.  For example, if there is a
sequence of points that are exponentially approaching a fixed point.
In this case the algorithm can be slightly modified so that it still
runs in $O(n \log n)$ time.  This modified version is outlined as
follows.

\begin{enumerate}
\item Keep the points in linked lists sorted by each dimension
\item In the selected dimension, come in from both sides until the cut is
found
\item Remove cut elements and put aside
\item Repeat making cuts until the size of the largest subset is less than
$(2/3)n$
\item Create subsets and make recursive calls.
\end{enumerate}

The runtime for building the tree is thus specified by

\[T(n) = \sum_{i=1}^{k} T(n_i) + O(n)\] such that $\sum_{i=1}^{k} n_i = n$
and $\mbox{max}_{i=1}^{k}(n_i) = \frac{2}{3}n$.

This solves to $O(n \log n)$ independent of how inbalanced the tree
ends up being.  A problem is that it does not parallelize easily.
Callahan and Kosaraju show how it can be parallelized, but the
parallel version is even more complicated and is not work-efficient
(the processor-time product is $O(n \log^2 n)$).
We will not cover this complicated parallel version since in practice
parallelizing the simple version is both more efficient and likely
to be highly effective.

\newcommand{\wsr}{WSR}
\newcommand{\wsrp}{WSR\_Pair}
\subsection{Algorithm: Well-Separated Realization}

Once we have a decomposition tree we can use it to build the
well-separated realization.  The following algorithm returns the
set of interaction edges in the well-separated realization.
The basic idea is to first generate interaction edges within
the left and right children (recursive calls to \wsr) and then generate
interaction edges between the two children (call to \wsrp).

\begin{quote}
\begin{tabbing}
{\bf function} \wsr($T$)\\
~~~~// $T$ is a decomposition tree\\[.08in]
~~~~{\bf if} Is\_Leaf($T$) {\bf return} $\emptyset$\\
~~~~{\bf else}\\
~~~~~~~~{\bf return} \=\wsr(Left($T$)) $\cup$ \\
                     \>\wsr(Right($T$)) $\cup$ \\
                     \>\wsrp(Left($T$),Right($T$))\\
{\bf end}
\end{tabbing}
\end{quote}

\begin{quote}
\begin{tabbing}
{\bf function} \wsrp($T_1,T_2$)\\
~~~~// $T_1$ and $T_2$ are two disjoint decomposition trees\\[.08in]
~~~~{\bf if} Well\_Separated($T_1$, $T_2$) {\bf return} \{Edge($T_1,T_2$)\}\\
~~~~{\bf else}\\
~~~~~~~~{\bf if} $(l_{\max}(T_1) > l_{\max}(T_2))$\\
~~~~~~~~~~~~{\bf return} \wsrp(Left($T_1$), $T_2$) 
                    $\cup$ \wsrp(Right($T_1$), $T_2$)\\
~~~~~~~~{\bf else}\\
~~~~~~~~~~~~{\bf return} \wsrp($T_1$, Left($T_2$)) 
                    $\cup$ \wsrp($T_1$, Right($T_2$))\\
{\bf end}
\end{tabbing}
\end{quote}
The idea of \wsrp($T_1$, $T_2$) is to return an interaction edge if
$T_1$ and $T_2$ are well separated.  Otherwise it picks the node with
a larger $l_{\max}$ and splits that node and calls itself recursively
on the smaller and the two halves of the larger.

Figure~\ref{real_algo} gives an example of the algorithm along with
call tree for \wsrp$(T_1,T_2)$.  In the example the algorithm
determines that $l_{\max}(T_2) > l_{\max}(T_1)$ and therefore splits
$T_2$ in the recursive calls.  At the next level it splits $T_1$, at
which point it determines that $T_{11}$ is well separated from
$T_{21}$ and therefore makes no more recursive calls on that branch.


\begin{figure}[ht]
\centerline{\psfig{file=../lecture16/real_algo.eps,width=300pt}}
\caption{Generating the Realization.  The tree represents
the call tree for \wsrp.}
\label{real_algo}
\end{figure}

\subsection{Bounding the size of the realization}

By bounding the size of the realization we both bound the number of
interaction edges and also bound the runtime for building the
realization since the runtime of \wsr{} is proportional to the
number of edges it creates.

We will use the following to bound the total number of interaction edges
generated by \wsr{} (size of realization):

\begin{figure}
\centerline{
  \psfig{file=../lecture16/cube_touch.eps,width=100pt}
}
\caption{Non-intersecting rectangles}\label{cube_touch}
\end{figure}

\begin{itemize}
\item Properties of the tree nodes - show bounding rectangles 
      are not ``too thin''
\item Bound \# of non-intersecting rectangles that can touch a cube of
fixed size (Figure~\ref{cube_touch}).  This assumes the rectangles
are not too thin.
\item Bound \# of calls to \wsrp($T_1$, $T_2$) for each cell
$T_1$ or $T_2$.
\end{itemize}

\begin{figure}
\centerline{
  \psfig{file=../lecture16/outer_rect.eps,width=200pt}
}
\caption{Outer Rectangle}\label{outer_rect}
\end{figure}

We define an {\em outer rectangle} for each bounding rectangle such that
these outer rectangles fill space ({\em i.e.}, the union of all rectangles
at a given level of a tree fill the initial space).
The outer rectangle for each node in the decomposition tree is
defined as follows and illustrated in Figure~\ref{outer_rect}.

\[ 
   \hat{R}(P) = \left\{ 
      \begin{array}{ll}
         \mbox{Smallest enclosing cube} \hspace*{.3in} & \mbox{at root} \\
            & \\
         \mbox{Use hyperplane that} & \\
         \mbox{divides $R(P)$ to} & \mbox{internally} \\
         \mbox{divide } \hat{R}(P) &
       \end{array}
    \right. \]

\begin{lemma} 
$l_{\min}(\hat{R}(A)) \geq \frac{l_{\max}(p(A))}{2}$, where
$p(A) = $ parent of A in the tree.
\end{lemma}

\proof{
By induction.

\noindent
{\bf Base:} true when $p(A) = $ root, since
      $l_{\min}(\hat{R}(A)) = \frac{l_{\max}(p(A))}{2}$.

\noindent
{\bf Induction:} if true for $p(A)$ then true for $A$, 

\begin{description}
\item[Case 1:] $l_{\min}(\hat{R}(A)) = l_{\min}(\hat{R}(p(A)))$. 
True, since $l_{\max}(p(A)) \leq l_{\max}(p(p(A)))$ and by induction
$l_{\min}(\hat{R}(p(A))) = \frac{l_{\max}(p(p(A)))}{2}$.

\item[Case 2:] $l_{\min}(\hat{R}(A)) < l_{\min}(\hat{R}(p(A)))$.
First we note that $l_{\max}(p(A))$ and
$l_{\min}(\hat{R}(A))$ must be along the same dimension since
Build\_Tree cut perpendicular to $l_{\max}(p(A))$ and if
$l_{\min}(\hat{R}(A))$ has decreased it must have been cut.
Since $\hat{R}(p(A))$ contains $R(p(A))$ cutting
$R(p(A))$ in half will mean that
$l_{\min}(\hat{R}(A)) \geq \frac{l_{\max}(p(A))}{2}$.
\end{description}
}

\begin{figure}
\centerline{\psfig{file=../lecture16/lemmaproof.eps,width=200pt}}
\caption{Proof of Lemma~\protect\ref{lemma:cubes}}
\label{lemmaproof}
\end{figure}

\begin{lemma}
Given a d-cube $C$ and a set $S$ of disjoint d-cubes of length
$l'$ that overlap $C$, then \[ |S| \leq {\left( \frac{l(C)}{l'} + 2
\right)}^d \].
\label{lemma:cubes}
\end{lemma}

\proof{By picture (Figure~\ref{lemmaproof}).}

\begin{figure}
\centerline{\psfig{file=../lecture16/moreproof.eps,width=300pt}}
\caption{Well separatedness}
\label{moreproof}
\end{figure}

It's not hard to show that if $l_{\max}(P_1) \geq l_{\max}(P_2)$, then for
$P_2$ not to be well separated from $P_1$ it must intersect a d-cube of
length $(s\sqrt{d} + 2\sqrt{d} + 1) l_{\max}(P_1)$ centered at the center of
$P_1$ (Figure~\ref{moreproof}).

Note: we will only split $P_1$ in \wsrp($P_1$, $P_2$) if $P_1$ and
$P_2$ are not well separated.

\subsubsection*{Putting everything together}

Consider for a given node P of the decomposition tree, all invocations
\wsrp($P$,$P'$) such that $l_{\max}(P) \geq l_{\max}(P')$. 
\begin{enumerate}
\item all $P'$ must be independent (consider recursive structure of
      the algorithm \wsrp)
\item we know $l_{\min}(\hat{R}(P')) \geq \frac{l_{\max}(p(P'))}{2}$ so we
can bound the number of $P'$ to \[ |P'| < {(2(s\sqrt{d} + 2\sqrt{d} + 1) + 2)^d}
\].
\end{enumerate}

Therefore, the total \# of calls to \wsrp($P_1$,$P_2$) is bound by \[ 2n
{(2(s\sqrt{d} + 2\sqrt{d} + 1) + 2)}^d = O(n)\] 
The factor of two
at the front comes from the fact that we are counting non-leaf invocations of
\wsrp($T_1$,$T_2$) while the interactions are created by leaf invocations.
However, the total number of leafs is at most twice the number of
internal nodes---see Figure~\ref{leaf_only}.

We note that this bound is very conservative and in practice the algorithm
actually uses many fewer interactions.

\begin{figure}
\centerline{
  \psfig{file=../lecture16/leaf_only.eps,width=300pt}
}
\caption{Counting leaves only}\label{leaf_only}
\end{figure}

\section{N-body Summary}

\begin{itemize}

\item {\bf Algorithms}
\begin{itemize}
\item Particle-Particle
\item Particle-mesh (several variants)
\item Particle-cell (Appel, Barnes-Hut, Press)
\item Cell-cell (FMM, Callahan-Kosaraju)
\end{itemize}

\item {\bf Applications}
\begin{itemize}
\item Astronomy
\item Molecular Dynamics
\item Fluid Dynamics (elementless methods)
\item Plasma Physics
\end{itemize}

\end{itemize}

