Adsii PDF

Algorithms and Data Structures II
(WS 20/21)
Course Notes
Christian Schulz
Heidelberg University
source: xkcd
Last modified: January 4, 2021

Contents
Preface 1
1 Algorithm Engineering 3
1.1 A Brief “History” of Algorithm Engineering . . . . . . . . . . . . 6
1.2 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.7 Algorithm Libraries . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.8 Instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.9 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5 Maximum Flows and Matchings 15

5.1 Algorithms 1956–now . . . . . . . . . . . . . . . . . . . . . . . 16
5.2 Augmenting Path Algorithms . . . . . . . . . . . . . . . . . . . . 17
5.2.1 The Ford Fulkerson Algorithm . . . . . . . . . . . . . . . 18
5.2.2 Two Bad Examples for Ford Fulkerson . . . . . . . . . . 21
5.2.3 Dinic’s Algorithm . . . . . . . . . . . . . . . . . . . . . 21
5.2.4 Disadvantage of Augmenting Path Algorithms . . . . . . 25
5.3 Preflow-Push Algorithms . . . . . . . . . . . . . . . . . . . . . . 26
5.3.1 Heuristic Improvements . . . . . . . . . . . . . . . . . . 34
5.3.2 Experimental Results . . . . . . . . . . . . . . . . . . . . 36
5.4 Minimum Cost Flows . . . . . . . . . . . . . . . . . . . . . . . . 38
5.5 Matching Problems . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.5.1 Maximum Cardinality Bipartite Matching . . . . . . . . . 41
5.6 Book References . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6 Generic Approaches to Optimization 43
i
ii CONTENTS
7 Linear Programming 45
8 Parallel Algorithms 47
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
8.2 Message-Passing Model/Paradigm . . . . . . . . . . . . . . . . . 47
8.3 Communication Cost in Parallel Machines . . . . . . . . . . . . . 49
8.4 Performance metrics . . . . . . . . . . . . . . . . . . . . . . . . 51
8.4.1 Execution Time . . . . . . . . . . . . . . . . . . . . . . . 51
8.4.2 Speedup . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
8.4.3 Superlinear speedup . . . . . . . . . . . . . . . . . . . . 52
8.4.4 Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . 52
8.4.5 Example: Adding n numbers using n processing elements 53
8.4.6 Example: Adding n numbers using p processing elements 54
8.5 Prefix Sums on a Hypercube Network . . . . . . . . . . . . . . . 55
8.6 Quicksort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
8.6.1 Sequential Algorithm . . . . . . . . . . . . . . . . . . . . 57
8.6.2 Parallelizing Quicksort . . . . . . . . . . . . . . . . . . . 58
8.6.3 Parallel Formulation for a Message-passing Platform . . . 59
8.7 Mergesort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
9 External Algorithms 63
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
9.2 The external memory model . . . . . . . . . . . . . . . . . . . . 65
9.3 Stacks and Queues . . . . . . . . . . . . . . . . . . . . . . . . . 66
9.4 Multiway Merge Sort . . . . . . . . . . . . . . . . . . . . . . . . 67
9.5 Internal work of Multiway Mergesort . . . . . . . . . . . . . . . 69
9.6 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
9.7 External Priority Queues . . . . . . . . . . . . . . . . . . . . . . 73
9.8 Semiexternal Kruskal’s Algorithm . . . . . . . . . . . . . . . . . 78
10 Online Algorithms 81
10.1 Competitive Analysis . . . . . . . . . . . . . . . . . . . . . . . . 82
10.1.1 Example: Ski Rental . . . . . . . . . . . . . . . . . . . . 82
10.2 Paging – Deterministic Algorithms . . . . . . . . . . . . . . . . . 83
10.2.1 Longest Forward Distance . . . . . . . . . . . . . . . . . 86
10.2.2 Competitive Paging Algorithms . . . . . . . . . . . . . . 87
10.2.3 Conservative Paging Algorithms . . . . . . . . . . . . . . 89
10.3 Randomized Online Algorithms . . . . . . . . . . . . . . . . . . 91
10.3.1 Randomized Marking Algorithms . . . . . . . . . . . . . 92
10.3.2 Lower Bound for Randomized Paging Algorithms . . . . 94
CONTENTS iii
11 Geometry 99
11.1 Application Domains . . . . . . . . . . . . . . . . . . . . . . . . 99
11.2 Basic Elements and Definitions . . . . . . . . . . . . . . . . . . . 99
11.3 Data Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
11.4 Line Segment Intersection - Case Study . . . . . . . . . . . . . . 101
11.4.1 Orthogonal Line Segment Intersection . . . . . . . . . . . 101
11.4.2 A more Complex Case . . . . . . . . . . . . . . . . . . . 102
11.4.3 The Arbitrary Case . . . . . . . . . . . . . . . . . . . . . 104
11.4.4 Further Improvements . . . . . . . . . . . . . . . . . . . 105
11.5 Convex Hull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
11.5.1 Graham Scan . . . . . . . . . . . . . . . . . . . . . . . . 106
11.6 Smallest Enclosing Circle . . . . . . . . . . . . . . . . . . . . . . 107
11.6.1 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 108
11.7 Range Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
11.7.1 1D Range Search . . . . . . . . . . . . . . . . . . . . . . 109
11.7.2 2D Range Search . . . . . . . . . . . . . . . . . . . . . . 110
11.7.3 Bitvectors . . . . . . . . . . . . . . . . . . . . . . . . . . 114
Bibliography 115
iv CONTENTS
Preface
The present course notes collect and elaborate on material for the course “Algo-
rithmen and Data Structures II” at Heidelberg University. Some parts of the script
were first published in winter term 2010/2011 for a course “Algorithmen II” at
KIT and have been expanded by additional material from an “Algorithms and Data
Structures II” course that I gave in Vienna. The initial version of the writeup is due
to Timo Bingmann, Johannes Fischer, Robert Geisberger, Moritz Kobitzsch, Vitaly
Osipov, Peter Sanders, Dennis Schieferdecker, Christian Schulz (me), Johannes
Singler. Nevertheless, the current state of the course notes is still preliminary –
in software jargon we would call it an α-version. Also, take note that the topics
presented in here sometimes go beyond what was presented during the lectures.
This aditional material is not required to successfully complete the course, it is
given to broaden your algorithmic horizon. Note that little here is original material
and kindly follow the citations to the original sources for more detailed information.
We will extend the chapters of the scriptum while the course is running. Also
NOTE that some chapters go far beyond the material that we covered in the lecture.
The chapters in these notes directly correspond to the chapters on the lecture
slides. For chapters 2–4 on Advanced Data Structures, Shortest Paths and Applica-
tions of DFS we refer to the textbook "Algorithms and Data Structures: The Basic
Toolbox" by Kurt Mehlhorn and Peter Sanders [38], that is also available online
and for free1 .
1
see https://people.mpi-inf.mpg.de/~mehlhorn/Toolbox.html
1
2 PREFACE
Chapter 1
Algorithm Engineering
realistic
algorithm models 3
engineering real
9 10
Inputs
applications
design
4
falsifiable
analysis 5 hypotheses 7 experiments
induction
deduction 6
implementation appl. engineering
perf.−
guarantees 8
algorithm−
libraries
Figure 1.1: Algorithm engineering as a cycle of design, analysis, implementation,

and experimental evaluation driven by falsifiable hypotheses. The numbers refer to
sections.
This is a shortened version of [48]. Algorithms and data structures are at

the heart of every computer application and thus of decisive importance for per-
manently growing areas of engineering, economy, science, and daily life. The
subject of Algorithmics is the systematic development of efficient algorithms
3
4 CHAPTER 1. ALGORITHM ENGINEERING
and therefore has pivotal influence on the effective development of reliable and
resource-conserving technology. We only mention a few spectacular examples.
Fast search in the huge data space of the internet (e.g. using Google) has
changed the way we handle knowledge. This was made possible with full-text
search algorithms that are harvesting matching information out of petabytes of data
within fractions of a second and by ranking algorithms that process graphs with
billions of nodes in order to filter relevant information out of heaps of result data.
Less visible yet similarly important are algorithms for the efficient distribution
and caching of frequently accessed data under massive load fluctuations or even
distributed denial of service attacks.
One of the most far-reaching results of the last years was the ability to read the
human genome. Algorithmics was decisive for the early success of this project
[55]. Rather than just processing the data coming our of the lab, algorithmic con-
siderations shaped the implementation of the applied shotgun sequencing process.
The list of areas where sophisticated algorithms play a key role could be arbi-
trarily continued: computer graphics, image processing, geographic information
systems, cryptography, planning in production, logistics and transportation,. . .
How is algorithmic innovation transferred to applications? Traditionally, algo-
rithmics used the methodology of algorithm theory which stems from mathematics:
algorithms are designed using simple models of problem and machine. Main results
are provable performance guarantees for all possible inputs. This approach often
leads to elegant, timeless solutions that can be adapted to many applications. The
hard performance guarantees lead to reliably high efficiency even for types of inputs
that were unknown at implementation time. From the point of view of algorithm
theory, taking up and implementing an algorithmic idea is part of application devel-
opment. Unfortunately, it can be universally observed that this mode of transferring
results is a slow process. With growing requirements for innovative algorithms,
this causes growing gaps between theory and practice: Realistic hardware with its
parallelism, memory hierarchies etc. is diverging from traditional machine models.
Applications grow more and more complex. At the same time, algorithm theory
develops more and more elaborate algorithms that may contain important ideas but
are usually not directly implementable. Furthermore, real-world inputs are often
far away from the worst case scenarios of the theoretical analysis. In extreme cases,
promising algorithmic approaches are neglected because a mathematical analysis
would be difficult.
Since the early 1990s it therefore became more and more apparent that al-
gorithmics cannot restrict itself to theory. So, what else should algorithmicists
do? Experiments play a pivotal here. Algorithm engineering (AE) is therefore
sometimes equated with experimental algorithmics. However, in this paper we
argue that this view is too limited. First of all, to do experiments, you also have to
implement algorithms. This is often equally interesting and revealing as the experi-
5
ments themselves, needs its own set of techniques, and is an important interface to
software engineering. Furthermore, it makes little sense to view design and analysis
on the one hand and implementation and experimentation on the other hand as
separate activities. Rather, a feedback loop of design, analysis, implementation,
and experimentation that leads to new design ideas materializes as the central
process of algorithmics.
This cycle is quite similar to the cycle of theory building and experimental
validation in Popper’s scientific method [45]. We can learn several things from
this comparison. First, this cycle is driven by falsifiable hypotheses validated
by experiments – an experiment cannot prove a hypothesis but it can support it.
However, such support is only meaningful if there are conceivable outcomes of
experiments that prove the hypothesis wrong. Hypotheses can come from creative
ideas or result from inductive reasoning stemming from previous experiments.
Thus we see a fundamental difference to the deductive reasoning predominant in
algorithm theory. Experiments have to be reproducible, i.e., other researchers have
to be able to repeat an experiment to the extent that they draw the same conclusions
or uncover mistakes in the previous experimental setup.
There are further aspects of AE as a methodology for algorithmics, outside
the main cycle. Design, analysis and evaluation of algorithms are based on some
model of the problem and the underlying machine. Since gaps between theory
and practice often relate to these models, they are an important aspect of AE.
Since we aim at practicality, applications are an important aspect. However
we choose to view applications as being outside the methodology of AE since
it would otherwise become too open ended and because often one algorithm
can be used for quite diverse applications. Also, every new application will
have its own requirements and techniques some of which may be abstracted
away for algorithmic treatment. Still, in order to reduce gaps between theory
and practice, as many interactions as poissible between the application and the
activities of AE should be taken into account: Applications are the basis for realistic
models, they influence the kind of analysis we do, they put constraints on useful
implementations, and they supply realistic inputs and other design parameters for
experiments. On the other hand, the results of analysis and experiments influence
the way an algorithm is used (fast enough for real time or interactive use?,. . . )
and implementations may be the basis for software used in applications. Indeed,
we may view application engineering as a separate process living in both AE
and a concrete application domain where methods from both areas are used to
adapt an algorithm to a particular application. Applications engineering bridges
remaining unavoidable gaps between experimental implementations and production
quality code. Note that there are important differences between these two kinds of
code: fast development, efficiency, and instrumentation for experiments are very
important for AE, while thorough testing, maintainability, simplicity, and tuning for
particular classes of inputs are more important for the applications. Furthermore,
the algorithm engineers may not even know all the applications for which their
algorithms will be used. Hence, algorithm libraries of highly tested codes with
clear simple user interfaces are an important link between AE and applications.
Figure 1.1 summarizes the resulting schema for AE as a methodology for
algorithmics. The following sections will describe the activities in more detail. We
give examples of challenges and results that are a more or less random sample
biased to results we know well.
1.1 A Brief “History” of Algorithm Engineering

The methodology described here is not intended as a revolution but as a description
of observed practices in algorithmic research being compiled into a consistent
methodology. Basically, all the activities in algorithm development described here
have probably been used as long as there are computers. However, in the 1970s
and 1980s algorithm theory had become a subdiscipline of computer science that
was almost exclusively devoted to “paper and pencil” work. Except for a few
papers around D. Johnson, the other activities were mostly visible in application
papers, in operations research, or J. Bentley’s programming pearls column in
Communications of the ACM. In the late 1980s, people within algorithm theory
began to notice increasing gaps between theory and practice leading to important
activities such as the Library of Efficient Data Types and Algorithms (LEDA, since
1988) by K. Mehlhorn and S. Näher and the DIMACS implementation challenges
(http://dimacs.rutgers.edu/Challenges/). It was not before the
end of the 1990s that several workshops series on experimental algorithmics and
algorithm engineering were started.1 There was a Dagstuhl workshop in 2000 [19],
and several overview papers on the subject were published [2, 41, 30, 31, 25].
The term “algorithm engineering” already appears 1986 in the Foreword of [9]
and 1989 in the title of [54]. No discussion of the term is given. At the same time T.
Beth started an initiative to move the CS department of the University of Karlsruhe
more into the direction of an engineering discipline. For example, a new compul-
sory graduate-level course on algorithms was called “Algorithmentechnik” which
can be translated as “algorithm engineering”. Note that the term “engineering”
like in “mechanical engineering” means the application oriented use of science
whereas our current interpretation of algorithm engineering has applications not as
its sole objective but equally strives for general scientific insight as in the natural
sciences. However, in daily work the difference will not matter much.
1
The Workshop on Algorithm Engineering (WAE) is not the engineering track of ESA. The
Alex workshop first held in Italy in 1998 is now the ALENEX workshop held in conjuction with
SODA. WEA, now SEA was first organized in 2002.
1.2. MODELS 7
P. Italiano organized the “Workshop on Algorithm Engineering” in 1997 and

also uses “algorithm engineering” as the title for the algorithms column of EATCS
in 2003 [17] with the following short abstract: “Algorithm Engineering is concerned
with the design, analysis, implementation, tuning, debugging and experimental
evaluation of computer programs for solving algorithmic problems. It provides
methodologies and tools for developing and engineering efficient algorithmic codes
and aims at integrating and reinforcing traditional theoretical approaches for the
design and analysis of algorithms and data structures.” Independently but with the
same basic meaning, the term was used in the influential policy paper [2]. The
present paper basically follows the same line of argumentation attempting to work
out the methodology in more detail and providing a number of hopefully interesting
examples.
1.2 Models
A big difficulty for defining models for problems and machines is that (apparently)
only complex models are adequate images of reality whereas only simple models
lead to simple, widely usable, portable, and analyzable algorithms. Therefore, AE
must simultaneously and carefully abstract from application problems and refine
theoretical models.
A successful example for a machine model is the external memory model (or
I/O model) [1, 56, 39] which is a careful refinement of the von Neumann model
[42]. Instead of a uniform memory, there are two levels of memory. A fast memory
of limited size M and and a slow memory that is accessed in blocks of size B.
While only counting I/O steps in this model can become a highly theoretical game,
we get an abstraction useful for AE if we additionally take internal work into
account and if we are careful to use the right values for the parameters M and
B 2 . Algorithms good in the I/O model are often good in practice although the
model oversimplifies aspects like rotational delays, seek time, disk data density
depending on the track use, cache replacement strategies [35], flexible block sizes,
etc. Sometimes it would even be counterproductive to be too clever. For example,
a program carefully tuned to minimize rotational delays and seek time might
experience severe performance degradation as soon as another application accesses
the disk.
An example for application modelling is the simulation of traffic flows. While
microscopic simulations that take the actual behavior of every car into account
2
A common pitfall when applying the I/O model to disks is to look for a natural physical block
size. This can lead to values (e.g. the size of 512 byte for a decoding unit) that are four orders of
magnitude from the value that should be chosen – a value where data transfer takes about as long as
the average latency for a small block.
are currently limited to fairly small subnetworks, it may soon become possible to
simulate an entire country by only looking at the paths taken by each car.
1.3 Design
As in algorithm theory, we are interested in efficient algorithms. However, in AE,
it is equally important to look for simplicity, implementability, and possibilities for
code reuse. Furthermore, efficiency means not just asymptotic worst case efficiency,
but we also have to look at the constant factors involved and at the performance
for real-world inputs. In particular, some theoretically efficient algorithms have
similar best case and worse case behavior whereas the algorithms used in practice
perform much better on all but contrived examples. An interesting example are
maximum flow algorithms where the asymptotically best algorithm [21] is much
worse than theoretically inferior algorithms [13, 34].
1.4 Analysis
Even simple and proven practical algorithms are often difficult to analyze and
this is one of the main reasons for gaps between theory and practice. Thus, the
analysis of such algorithms is an important aspect of AE. For example, randomized
algorithms are often simpler and faster than their best deterministic competitors
but even simple randomized algorithms are often difficult to analyze.
Many complex optimization problems are attacked using meta heuristics like
(randomized) local search or evolutionary algorithms. Algorithms of this type
are simple and easily adaptable to the problem at hand. However, only very few
such algorithms have been successfully analyzed (e.g. [57]) although performance
guarantees would be of great theoretical and practical value.
An important open problem is partitioning of graphs into approximately equal
sized blocks such that few edges are cut. This problem has many applications,
e.g., in scientific computing. Currently available algorithms with performance
guarantees are too slow for practical use. Practical methods first contract the graph
while preserving its basic structure until only few nodes are left, compute an initial
solution on this coarse representation, and then improve by local search. These
algorithms, e.g., [20] are very fast and yield good solutions in many situations yet
no performance guarantees are known.
An even more famous example for local search is the simplex algorithm for
linear programming. Simple variants of the simplex algorithm need exponential
time for specially constructed inputs. However, in practice, a linear number
of iterations suffices. So far, only subexponential expected runtime bounds are
1.5. IMPLEMENTATION 9
known – for inpracticable variants. However, Spielmann and Teng were able to
show that even small random perturbations of the coefficients of a linear program
suffice to make the expected run time of the simplex algorithm polynomial [53].
This concept of smoothed analysis is a generalization of average case analysis
and an interesting tool of AE also outside the simplex algorithm. Beier and
Vöcking were able to show polynomial smoothed complexity for an important
family of NP-hard problems [7]. For example, this result explains why the knapsack
problem can be efficiently solved in practice and has also helped to improve the
best knapsack solvers. There are interesting interrelations between smoothed
complexity, approximation algorithms, and pseudopolynomial algorithms that is
also an interesting approach to practical solutions of NP-hard problems.
1.5 Implementation
Implementation only appears to be the most clearly prescribed and boring activity in
the cycle of AE. One reason is that there are huge semantic gaps between abstractly
formulated algorithms, imperative programming languages, and real hardware. A
typical example for this semantic gap is the implementation of an O(nm log n)
matching algorithm in [36]. Its abstract description requires a sophisticated data
structure whose efficient implementation only succeeded in [36].
An extreme example for the semantic gap are geometric algorithms which
are often designed assuming exact arithmetics with real numbers and without
considering degenerate cases. The robustness of geometric algorithms has therefore
become an important branch of AE [12, 8, 6].
Even the implementation of relatively simple basic algorithms can be challeng-
ing. You often have to compare several candidates based on small constant factors
in their execution time. Since even small implementation details can make a big
difference, the only reliable way is to highly tune all competitors and run them on
several architectures. It can even be advisable to compare the generated machine
code (e.g. [47, 51],[11]).
Often only implementations give convincing evidence of the correctness and
result quality of an algorithm. For example, an algorithm for planar embedding
[23] was the standard reference for 20 years although this paper only contains
a vague description how an algorithm for planarity testing can be generalized.
Several attempts at a more detailed description contained errors (e.g. [32]). This
was only noticed during the first correct implementation [33]. Similarly, for a long
time nobody suceeded in implementing the famous algorithm for computing three-
connected components from [24]. Only an implementation in 2000 [22] uncovered
and corrected an error. For the related problem of computing a maximal planar
subgraph there was a series of publications in prominent conferences uncovering
errors in the previous paper and introducing new ones – until it turned out that the
proposed underlying data structure is inadequate for the problem [26].
An important consequence for planning AE projects is that important imple-
mentations cannot usually be done as bachelor or master theses but require the
very best students or long term attendance by full time researchers or scientific
programmers.
1.6 Experiments
Meaningful experiments are the key to closing the cycle of the AE process. For
example, experiments on crossing minimization in [27] showed that previous
theoretical results were too optimistic so that new algorithms became interesting.
Experiments can also have a direct influence on the analysis. For example,
reconstructing a curve from a set of measured points is a fundamental variant of an
important family of image processing problems. In [4] an apparently expensive
method based on the travelling salesman problem is investigated. Experiments
indicated that “reasonable” inputs lead to easy instances of the travelling salesman
problem. This observation was later formalized and proven. A quite different
example of the same effect is the astonishing observation that arbitrary access
patterns to data blocks on disk arrays can be almost perfectly balanced when two
redundant copies of each block are placed on random disks [49].
Compared to the natural sciences, AE is in the privileged situation that it can
perform many experiments with relatively little effort. However, the other side of
the coin is highly nontrivial planning, evaluation, archiving, postprocessing, and
interpretation of results. The starting point should always be falsifiable hypotheses
on the behavior of the investigated algorithms which stem from the design, analysis,
implementation, or from previous experiments. The result is a confirmation,
falsification, or refinement of the hypothesis. The results complement the analytic
performance guarantees, lead to a better understanding of the algorithms, and
provide ideas for improved algorithms, more accurate analysis, or more efficient
implementation.
Successful experimentation involves a lot of software engineering. Modular
implementations allow flexible experiments. Clever use of tools simplifies the eval-
uation. Careful documentation and version management help with reproducibility
– a central requirement of scientific experiments, that is challenging due to the
frequent new versions of software and hardware.
1.7. ALGORITHM LIBRARIES 11
1.7 Algorithm Libraries

Algorithm libraries are made by assembling implementations of a number of algo-
rithms using the methods of software engineering. The result should be efficient,
easy to use, well documented, and portable. Algorithm libraries accelerate the
transfer of know-how into applications. Within algorithmics, libraries simplify
comparisons of algorithms and the construction of software that builds on them.
The software engineering involved is particularly challenging, since the applica-
tions to be supported are unknown at library implementation time and because
the separation of interface and (often highly complicated) implementation is very
important. Compared to applications-specific reimplementation, using a library
should save development time without leading to inferior performance. Compared
to simple, easy to implement algorithms, libraries should improve performance.
In particular for basic data structures with their fine-grained coupling between
applications and library this can be very difficult. To summarize, the triangle be-
tween generality, efficiency, and ease of use leads to challenging tradeoffs because
often optimizing one of these aspects will deteriorate the others. It is also worth
mentioning that correctness of algorithm libraries is even more important than for
other software because it is extremely difficult for a user to debug library code that
has not been written by his team. Sometimes it is not even sufficient for a library
to be correct as long as the user does not trust it sufficiently to first look for bugs
outside the library. This is one reason why result checking, certifying algorithms,
or even formal verification are an important aspect of algorithm libraries. All these
difficulties imply that implementing algorithms for use in a library is several times
more difficult / expensive / time consuming / frustrating /· · · than implementations
for experimental evaluation. On the other hand, a good library implementation
might be used orders of magnitude more frequently. Thus, in AE there is a natural
mechanism leading to many exploratory implementations and a few selected library
codes that build on previous experimental experience.
Let us now look at a few successful examples of algorithm libraries. The Library
of Efficient Data Types and Algorithms LEDA [34] has played an important part in
the development of AE. LEDA has an easy to use object-oriented C++ interfaces.
Besides basic algorithms and data structures, LEDA offers a variety of graph
algorithms and geometric algorithms.
Programming languages come with a run-time library that usually offers a few
algorithmic ingredients like sorting and various collection data structures (lists,
queues, sets, . . . ). For example, the C++ standard template library (STL) has a
very flexible interface based on templates. Since so many things are resolved at
compile time, programs that use the STL are often equally efficient as hand-written
C-style code even with the very fine-grained interfaces of collection classes. This
is one of the reasons why our group is looking at implementations of the STL
for advanced models of computation like external computing (STXXL [15]) or

multicore parallelism (MCSTL, GNU C++ standard library [52]). We should also
mention disadvantages of template based libraries: The more flexible their offered
functionality, the more cumbersome it is to use (the upcoming new C++ standard
might slightly improve the situation). Perhaps the worst aspect is coping with extra
long error messages and debugging code with thousands of tiny inlined functions.
Writing the library can be frustrating for an algorithmicist since the code tends
to consist mostly of trivial but lengthy declarations while the algorithm itself is
shredded into many isolated fragments.
The Boost C++ libraries (www.boost.org) are an interesting concept since
they offer a forum for library designers that ensures certain quality standards and
offers the possibility of a library to become part of the C++ standard.
The Computational Geometry Algorithms Library (CGAL) www.cgal.org
that is a joined effort of several AE groups is perhaps one of the most sophisticated
examples of C++ template programming. In particular, it offers many robust
implementations of geometric algorithms that are also efficient. This is achieved
for example by using floating point interval arithmetics most of the time and
switching to exact arithmetics only when a (near)-degenerate situation is detected.
The mechanisms of template programming make it possible to hide much of these
complicated mechanisms behind special number types that can be used in a similar
way as floating point numbers.
1.8 Instances
Collections of realistic problem instances for benchmarking have proven crucial
for improving algorithms. There are interesting collections for a number of NP-
hard problems like the travelling salesman problem 3 , the Steiner tree problem,
satisfiability, set covering, or graph partitioning. In particular for the first three
problems the benchmarks have helped enable astonishing breakthroughs. Using
deep mathematical insights into the structure of the problems one can now compute
optimal solutions even for large, realistic instances of the travelling salesman
problem [5] and of the Steiner tree problem [44]. It is a bit odd that similar
benchmarks for problems that are polynomially solvable are sometimes more
difficult to obtain. For route planning in road networks, realistic inputs have become
available in 2005 [50] enabling a revolution with speedups of up to six orders of
magnitude over Dijkstra’s algorithm and a perspective for many applications [14].
In string algorithms and data compression, real-world data is also no problem.
But for many typical graph problems like flows, random inputs are still common
3
http://www.iwr.uni-heidelberg.de/groups/comopt/software/
TSPLIB95/
1.9. APPLICATIONS 13
practice. We suspect that this often leads to unrealistic results in experimental

studies. Naively generated random instances are likely to be either much easier
or more difficult than realistic inputs. With more care and competition, such as
for the DIMACS implementation challenges, generators emerge that drive naive
algorithms into bad performance. While this process can lead to robust solutions,
it may overemphasize difficult inputs. Another area with lack of realistic input
collections are data structures. Apart from some specialized scenarios like IP
address lookup, few inputs are available for hash tables, search trees, or priority
queues.
1.9 Applications
We could discuss many important applications where algorithms play a major role
and a lot of interesting work remains to be done. Since this would go beyond
the scope of this paper, we only want to mention a few: Bioinformatics (e.g.
sequencing, folding, docking, phylogenetic trees, DNA chip evaluations, reaction
networks); information retrieval (indexing, ranking); algorithmic game theory;
traffic information, simulation and planning for cars, busses, trains, and air traffic;
geographic information systems; communication networks; machine learning; real
time scheduling.
The effort for implementing algorithms for a particular application usually
lies somewhere between the effort for experimental evaluation and for algorithm
libraries depending on the context.
An important goal for AE should be to help shaping the applications (as in
the example for genome sequencing mentioned in the introduction) rather than
act as an ancillary science for other disciplines like physics, biology, mechanical
engineering,. . .
1.10 Conclusions
We hope to have demonstrated that AE is a “round” methodology for the devel-
opment of efficient algorithms which simplifies their practical use. We want to
stress, however, that it is not our intention to abolish algorithm theory. The saying
that there is nothing as practical as good theory remains true for algorithmics
because an algorithm with proven performance guarantees has a degree of gen-
erality, reliability, and predictability that cannot be obtained with any number of
experiments. However, this does not contradict the proposed methodology since
it views algorithm theory as a subset of AE, making it even more rich by asking
additional interesting kinds of questions (e.g. simplicity of algorithms, care for
constant factors, smoothed analysis,. . . ). We also have no intention of criticizing

some highly interesting research in algorithm theory that is less motivated from
applications than by fundamental questions of theoretical computer science such as
computability or complexity theory. However, we do want to criticize those papers
that begin with a vague claim of relevance to some fashionable application area
before diving deep into theoretical constructions that look completely irrelevant
for the claimed application. Often this is not intentionally misleading but more
like a game of “Chinese whispers” where a research area starts as a sensible ab-
straction of an application area but then develops a life of itself, mutating into a
mathematical game with its own rules. Even this can be interesting but researchers
should constantly ask themselves why they are working on an area, whether there
are perhaps other areas where they can have larger impact on the world, and how
false claims for practicality can damage the reputation of algorithmics in practical
computer science.
Chapter 5
Maximum Flows and Matchings
The maximum flow problem is very easy to state: In a capacitated network, we wish
to send as much flow as possible between two special nodes, a source node s and a
sink node t, without exceeding the capacity of any edge. In this chapter we discuss
a number of algorithms for solving the maximum flow problem. These algorithms
are of two types:
1. Augmenting path algorithms which incrementally augment flow along
paths from the source node to the sink node.
2. Preflow-push algorithms that flood the network so that some nodes have
excesses (or buildup of flow). These algorithms incrementally relieve flow
from nodes with excesses by sending flow from the node forward toward the
sink node or backward toward the source node.
Definition 1
We define a network to be a directed weighted graph G = (V, E) with distinct
vertices s and t in G and the restriction that s has no incoming edges and t has
no outgoing edges. We call s a source node and t a sink node. For an edge e, we
call its weight ce the capacity of e (nonnegative!). An (s, t)-flow or simply flow is
a function f : E 7→ R≥0 (fe := f (e)) satisfying the capacity constraints and the
flow conservation constraints:
1. 0 ≤ fe ≤ ce ∀e ∈ E
P P
2. e:source(e)=v fe = e:target(e)=v fe ∀v ∈ V \{s, t}
The capacity constraints state that the flow across any edge is bounded by its
capacity, and the latter constraints state that for each node different from s and
t, the total incoming flow equals the total outgoing flow. We further define the
value of a flow val(f ) = total outgoing flow from s. Our goal is to find a flow with
maximum value. An example can be found in Figure 5.1.
15
16 CHAPTER 5. MAXIMUM FLOWS AND MATCHINGS
10
8 12
10
10 4 12
s t
8 4 2 6
4
10 4 8
Figure 5.1: An example network in which capacities are marked black and the
current flow values on an edge are marked blue. The total flow value is 18 which is
maximal. We also see a cut in the network induced by the colors.
5.1 Algorithms 1956–now

Algorithms for the Maximum Flow Problem have been studied for over five decades
now. Starting with early work in linear programming and spurred by the classic
book of Ford and Fulkerson, the study of this problem has led to continuing im-
provements in the efficiency of network flow algorithms. We give a brief overview
in the following table. Here n is the number of nodes, m the number of arcs/edges
and U the largest capacity in the graph.
Year Author Running time

1956 Ford-Fulkerson O(mnU )
1969 Edmonds-Karp O(m2 n)
1970 Dinic O(mn2 )
1973 Dinic-Gabow O(mn log U )
1974 Karzanov O(n3 )√
1977 Cherkassky O(n2 m)
1980 Galil-Naamad O(mn log2 n)
1983 Sleator-Tarjan O(mn log n)
1986 Goldberg-Tarjan O(mn log(n2 /m))
1987 Ahuja-Orlin O(mn + n2 log U√)
1987 Ahuja-Orlin-Tarjan O(mn log(2 + n log U /m))
1990 Cheriyan-Hagerup-Mehlhorn O(n3 / log n)
1990 Alon O(mn + n8/3 log n)
1992 King-Rao-Tarjan O(mn + n2+e )
1993 Philipps-Westbrook O(mn log n/ log m
n
+ n2 log2+ε n)
m
1994 King-Rao-Tarjan O(mn log n/ log n log n ) if m ≥ 2n log n
1997 Goldberg-Rao O(min{m1/2 , n2/3 }m log(n2 /m) log U )
5.2. AUGMENTING PATH ALGORITHMS 17
5.2 Augmenting Path Algorithms
The basic idea behind augmenting path algorithms is to iteratively identify a path
from s to t such that each edge has some spare capacity. On this path, the edges
with the smallest spare capacity are saturated, i.e. the flow along this path is
incremented by this capacity. This way the flow is incrementally augmented along
paths from the source node to the sink node. To be more precise we need the
definition of a residual graph:
Definition 2
Given a network G = (V, E, c) and a flow f , we define the residual network to be
Gf = (V, Ef , cf ). The residual network Gf with respect to a flow f has the same
node set as G. Every edge of Gf is induced by an edge of G and has so-called
residual capacity. Let e = (u, v) be an arbitrary edge of G. If f (e) < ce then e
is also an edge of Gf . Its residual capacity is cfe = ce − f (e). If f (e) > 0 then
erev = (v, u) is an edge of Gf . Its residual capacity is cferev = f (e). Note that both
can be the case, i.e. e ∈ Ef and erev ∈ Ef . An example is given in Figure 5.2.
In other words, the residual capacity cfe=(u,v) of an edge e is the maximum additional
flow that can be sent from node u to node v using the edges (u, v) and (v, u).
Sending flow from u to v via the edge (v, u) means to decrease the flow that is on
this edge, i.e. decreasing the flow from v to u.
10 10 10
10 12
10 10 12 12
10 10 44
10 40 12 10 12
2
6 4 2 4 6 4 22 44
4 4
10 8 capacity 10 8
4 4
flow 6
4 residual capacity 4 4
Figure 5.2: On the left hand side we see the original graph with a flow f
shown. On the right hand side we see the corresponding residual network. In
this graph the highlighted path is an augmenting path with residual capacity
∆f = min{4, 2, 10, 4, 4} = 2.
5.2.1 The Ford Fulkerson Algorithm
We now describe one of the simples and most intuitive algorithms for solving
the maximum flow problem. This algorithm is known as the augmenting path
algorithm/Ford Fulkerson algorithm. We refer to a directed path from the source
to the sink in the residual network as an augmenting path. We define the residual
capacity ∆f of an augmenting path p as the minimum residual capacity of any
edge in the path. An example for an augmenting path can be found in Figure 5.2.
Observe that, by definition, the capacity of an augmenting path is always positive.
Consequently, whenever the network contains an augmenting path, we can send
additional flow from the source to the sink. The augmenting path algorithm is
essentially based on this simple observation. The algorithm proceeds by identifying
augmenting paths and augmenting flows on these path until the network contains
no such path. Pseudocode for this algorithm can be found in Figure 5.3 and 5.4.
For integer capacities, the algorithm runs in O(mval(f )) time. An example run of
the algorithm can be found in Figure 5.5.
Function FFMaxFlow(G = (V, E), s, t, c : E → N) : E → N

f := 0
while ∃path p = (s, . . . , t) in Gf do
augmentPath(G, s, t, c, f ) // augment f along p
return f
Figure 5.3: The Ford Fulkerson Augmenting Path algorithm.
Procedure augmentPath(G = (V, E), s, t, c : E → N, f : E → N, path p in Gf )

∆f := min cfe
e∈p
foreach (u, v) ∈ p do
if (u, v) ∈ E then f(u,v) += ∆f
else f(v,u) −= ∆f
Figure 5.4: The Ford Fulkerson subroutine for augmenting a path.

10
10 12
4
s t
4
10 8
4
0 0
10 2 10 2
0 0
10 4 10 10 4 10
4 4 4 4
10 8 6 4 4
4 0
+10 +4
0 0
10 0 10 0
0 0
4 12 10 4 12
10
6 2 2 4 6 2 2 4
4 4 4 4 4 4
0 0
+2 +2
Figure 5.5: An example execution of the augmenting path algorithm.
We now proof that the Ford Fulkerson algorithm is correct.
Definition 3
An s-t-cut is a partition of V into S and T = V \S with s ∈ S and t ∈ T . The
capacity of this cut is the sum of the weights of all edges starting in S and ending
in T :
P
c(S) := c(S, T = V \S) := u∈S,v∈T c(u,v)
A cut S is called saturated if f (e) = ce for all e = (u, v), u ∈ S, v ∈ T and

f (e) = 0 for all e = (v, u), v ∈ T, u ∈ S.
The next lemma relates flows and cuts: the capacity of any (s, t)-cut is an upper
bound for the value of any (s, t)-flow. Conversely, the value of any (s, t)-flow is a
lower bound for the capacity of any (s, t)-cut.
Lemma 1
Let f be any (s, t)−flow and let S be any (s, t)-cut. Then
val(f ) ≤ c(S).
If S is saturated then val(f ) = c(S).

Proof.
We have
S→T edges T →S edges
z X}| { z X}| { X X
val(f ) = fe − fe ≤ fe ≤ ce = c(S).
e∈E∩S×T e∈E∩T ×S e∈E∩S×T e∈E∩S×T
For a saturated cut, the inequalities are equalities.

A saturated cut proves the maximality of given flow f . However a saturated
cut is easily extracted from a maximum flow by means of the residual network as
we will see now.
Lemma 2
Let f be a flow, let Gf be the residual network with respect to f , and let S := {v ∈
V : v reachable from s in Gf } be the set of nodes that are reachable from s in Gf .
1. If t ∈ S then f is not maximum.
2. If t 6∈ S then S is a saturated cut and f is maximum.
Proof.
1. Let p be any simple path from s to t in Gf and let ∆f be the minimum

residual capacity of any edge of p. Then ∆f > 0. We construct a flow f 0 of
value val(f )+∆f .

f (e) + ∆f if e is in p

0
f (e) = f (e) − ∆f if erev is in p

f (e) if neither e nor erev belongs to p.

Then f 0 is a flow and val(f 0 ) =val(f ) + ∆f . Note that f 0 is well-defined

because the path is simple.
2. There is no edge (v, w) in Gf with v ∈ S and w ∈ T . Hence, f (e) = c(e)

for any e ∈ E ∩ (S × T ) and f (e) = 0 for any e ∈ E ∩ (T × S), i.e. the cut
S is saturated. Thus f is maximal.
Lemma 3
The algorithm of Ford and Fulkerson is correct.
Proof.
Clearly the algorithm computes a feasible flow f (invariant). At termination
of the algorithm, we have no augmenting path left in Gf . We further have
∀(u, v) ∈ E ∩ (S × T ) : cf(u,v) = 0 and hence f(v,u) = 0. That means that S
is saturated and hence S proofs the maximality of f .
Theorem (Elias/Feinstein/Shannon, Ford/Fulkerson 1956)

The value of an s-t max-flow f equals the minimum capacity of an s-t-cut.
5.2.2 Two Bad Examples for Ford Fulkerson

We mentioned above that the running time for the Ford Fulkerson algorithm
is O(mval(f )) for integers. However the number of augmenting steps is also
dependent of the maximum difference capacity C := maxe∈E ce − mine∈E ce (see
Figure 5.6). It even gets worse since for irrational numbers it is possible that the
algorithm does not terminate at all (see Figure 5.7).
100 100 99 100 99 99

1 1 1
1 1 1 0 1
1 1 1
100 100 100 99 99 99
Figure 5.6: By iteratively augmenting the flow along the two blue paths, we get
200 =val(f ) path augmentations.
5.2.3 Dinic’s Algorithm

There are several maximum flow augmenting path algorithms that send flow along
shortest paths from the source to the sink. Dinic’s algorithm is a popular algorithm
in this class. This algorithm constructs shortest path networks, called layered
networks, and establishes blocking flows (to be defined later) in these networks. In
this section we point out the relationship between layered networks and distance
labels.
a
4 r 4
b
s 4 1 4 t
c
4 1 4
d
√
Figure 5.7: Here r is 5−1 2
. Let p0 = hs, c, b, ti, p1 = hs, a, b, c, d, ti,
p2 = hs, c, b, a, ti, p3 = hs, d, c, b, ti. Then the sequence of augmenting paths
p0 (p1 , p2 , p1 , p3 )∗ is an infinite sequence of positive flow augmentations and the
flow value does not converge to the maximum value 9.
Definition 4 (Distance Labels)

A distance function d : V → N ∪ {0} with respect to the residual capacities cfe is
a function from the set of nodes to the set of nonnegative integers. We say that a
distance function is valid with respect to a flow f if it satisfies the following two
conditions:
1. d(t) = 0
2. d(u) ≤ d(v) + 1 for every edge (u, v) ∈ Gf
We refer to d(u) as the distance label of node u and the conditions above as the
validity conditions. The following properties show why the distance labels might
be of use in designing network flow algorithms.
Property 1
If the distance labels are valid, the distance label d(u) is a lower bound on the
length of the shortest (directed) path from node u to node t in the residual network.
Property 2
If d(s) ≥ n, the residual network contains no directed path from the source node
to the sink node.
We now introduce some additional notation. We say that the distance labels are
exact if for each node u, d(u) equals the length of the shortest path from node u
to node t in the residual network. We can determine exact distance labels for all
nodes in O(m) time by performing a backward breadth-first search (reverseBFS)
of the network starting at the sink node t.
With respect to a given flow f , we can now define the layered network Lf as
follows. First we determine the exact distance labels d in Gf . Then we compute
Function DinitzMaxFlow(G = (V, E), s, t, c : E → N) : E → N

f := 0
while ∃path p = (s, . . . , t) in Gf do
d=Gf .reverseBFS(t) : V → N
Lf = (V, {(u, v) ∈ Ef : d(v) = d(u) − 1}) // layer graph
find a blocking flow fb in Lf
augment f += fb
return f
Figure 5.8: Pseudocode for Dinic’s algorithm.
the layered network which consists of those edges (u, v) in Gf that satisfy the
condition d(u) = d(v) + 1. Observe that by definition every path from the source
to the sink in the layered network Lf is a shortest path in Gf . Observe further that
some edge in Lf might not be contained in any path from the source to the sink.
Dinic’s algorithm now computes a blocking flow fb in Lf .
Definition 5
A blocking flow fb in Lf is a flow such that for each path p = hs, ..., ti in Lf there
exists an saturated edge e ∈ p, i.e. fb (e) = cf (e).
When a blocking flow fb has been established in a network, Dinic’s algorithm adds
fb to f , recomputes the exact distance labels, forms a new layered network, and
repeats these computations. The algorithm terminates, when as it is forming the
new layered networks, it finds that the source is not connected to the sink. It is
possible to show that every time Dinic’s algorithm forms a new layered network,
the distance label of the source node strictly increases. Consequently, Dinic’s
algorithm forms at most n layered networks and runs in O(n2 m) time, since as
shown below a blocking flow can be computed in O(nm). Pseudocode of this
algorithm can be found in Figure 5.8.
We now explain how to compute blocking flows in the layered network Lf .
The main idea is to repeatedly use depth first searches to find augmenting paths in
the layered network. The augmentation of flow along an edge (u, v) reduces the
residual capacity of edge (u, v) and increases the residual capacity of the reversal
edge (v, u). However, each edge of the layered network is admissible, and therefore
Dinic’s algorithm does not add reversal edges to the layered network. Consequently,
the length of every augmenting path is d(s) and in an augmenting path for every
edge (u, v) we have d(u) = d(v) + 1, i.e. the level decreases with each edge. These
facts allow us to determine an augmenting path in the layered network, on average,
in O(n) time. Pseudocode for the algorithm can be found in Figure 5.9. The
Function blockingFlow(H = (V, E)) : E → N

p=hsi : Path v=NodeRef : p.last()
fb := 0
loop // Round
if v = t then // breakthrough
δ:= min {c(e) − fb (e) : e ∈ p}
foreach e ∈ p do
fb (e) += δ
if fb (e) = c(e) then remove e from E
p:= hsi
else if ∃e = (v, w) ∈ E then p.pushBack(w) // extend
else if v = s then return fb // done
else delete the last edge from p in p and E // retreat
Figure 5.9: Pseudocode for computing a blocking flow. The main idea is to
repeatedly use DF S to find augmenting paths.
running time of this algorithm is basically #extends + #retreats + n · #breakthroughs .

We can have at most m breakthroughs since in each breakthrough at least one
edge is saturated. The number of #retreats is also bounded through m since
the edge is deleted from p and E. The number #extends is lower or equal to
#retreats + n · #breakthroughs since a retreat cancels 1 extend and a breakthrough
cancels ≤ n extends. Therefore the overall running time of this algorithm is
O(m + nm) = O(nm). If we have some further assumptions this time can even
get better. For example, if we have unit edge capacities then a breakthrough
saturates all edges on p, i.e. we have amortized constant cost per edge, yielding a
running time of O(m + n) for this subroutine.
Theorem 4
On√a unit network, i.e. unit edge capacities, Dinic’s algorithms makes at most
O( m) rounds.
Proof. √
Consider the layered network after m iterations. It has at
√ most m√edges. The
average number of edges between any two layers √ is m/ m = m. Hence,
there must be at least one layer i with at most m√edges. Now consider the
cut S = {v | d(v) < i}. It has capacity at most m. √ Since each blocking
flow augmentation increases flow by at least one, after m additional blocking
flow
√ augmentations,
√ a maximum flow has been found. Hence, no more than
2 m = O( m) blocking flow augmentations are needed.
2 10 1
a c 2 1
extend a c extend
10 10 12 2
extend extend
3 10 4 10 0 3 4 extend 0
s t breakthrough s t breakthrough
4 4
extend
10 2 1 8 10 1 8
2
b d b d
4 4
2 1 2 1
a c a c
3 4 extend 0 3 4 0
s retreat t breakthrough s t
2
extend extend
8 2 1 8 extend retreat 4
2 1 4
b d
4 extend b d
Figure 5.10: An example run of the blocking flow subroutine.
Putting these things together √

yields a total time for Dinic’s algorithm on unit
capacity networks of O((n + m) m). A more detailed analysis gives a running
time of O m min{(m1/2 , n2/3 } .
Another improvement for general networks was proposed by Sleator and Tarjan
in 1983. They developed the dynamic tree data structure and used it to improve
the worst-case complexity of Dinic’s algorithm from O(n2 m) to O(mn log n).
Dynamic trees are a special type of data structure that permits us to implicitly send
flow on paths of length n in O(log n) steps on average. Since then, researchers
have used this data structure on many occasions to improve the performance of a
range of network flow algorithms.
5.2.4 Disadvantage of Augmenting Path Algorithms

The inherent drawback of the augmenting path algorithms is the computationally
expensive operation of sending flow along a path, which requires O(n) time in the
worst case. Preflow-push algorithms do not suffer from this drawback and obtain
dramatic improvements in the worst-case complexity. To understand this point
better, consider the (artificially extreme) example shown in Figure 5.11.
When applied to this problem, any augmenting path algorithm would discover
l augmenting paths, each of length 6, and would augment 1 unit of flow along each
of these paths. Observe, however, that although all of these paths share the same
first eight 4 edges, each augmentation traverses all of these edges. If we could have
1
1 1
s u 1
8
t
1 1
l edges
Figure 5.11: Drawback of the augmenting path algorithm.
sent l units of flow from s to node u and then sent 1 unit of flow along l different
paths of length 2, we would have saved the repetitive computations in traversing
the common set of edges. This is the essential idea underlying the preflow-push
algorithms.
5.3 Preflow-Push Algorithms

We now study a class of algorithms, known as preflow-push algorithms, for solving
the maximum flow problem. These algorithms are more general, more powerful,
and more flexible than augmenting path algorithms. The best preflow-push algo-
rithms currently outperform the best augmenting path algorithms in theory as well
as in practice. In this section we study the generic preflow-push algorithm.
Augmenting path algorithms send flow by augmenting along a path. The basic
operation further decomposes into the more elementary operation of sending flow
along individual edges. Thus sending a flow of δ units along a path of k edges
decomposes into k basic operations of sending a flow of δ units along each of the
edges of the path. We shall refer to each of these basic operations as a push (see
Figure 5.12). The preflow-push algorithms push flows on individual edges instead
of augmenting paths.
Definition 6
A preflow f is a flow where the flow conservation constraint is relaxed to
inflow outflow
z X}| { z X}| {
excess(v):= f(u,v) − f(v,w) ≥ 0 .
(u,v)∈E (v,w)∈E
This implies for a preflow that it is allowed that the total inflow exceeds the total
outflow. A node v is called active if and only if excess(v) > 0. As an exception,
5.3. PREFLOW-PUSH ALGORITHMS 27
s and t are never active. The residual network Gf with respect to a preflow f is
defined as in the case of a flow.
The preflow-push algorithms maintain a preflow at each intermediate stage. In a

preflow, we have excess(v) ≥ 0 for each v ∈ V \{s, t}. Moreover, because no
edge is leaving from t we have excess(t) ≥ 0 as well. Therefore node s is the only
node with negative excess. The basic operation to manipulate a preflow is a push:
Procedure push(e = (v, w), δ)

assert δ > 0
assert residual capacity of e ≥ δ
assert excess(v) ≥ δ
excess(v) −= δ
excess(w) += δ
if e is reverse edge then f (reverse(e)) −= δ
else f (e) += δ
Figure 5.12: The basic push operation.
A push of δ units from node v to w decreases both excess(v) and cf(v,w) by δ

units and increases both excess(w) and cf(w,v) . We call a push saturating if δ = cfe
and non-saturating if δ < cfe . A non-saturating push at node v reduces excess(v)
to zero. A saturating push across e removes e from the residual network and either
kind of push adds erev to the residual network (if it is not already there).
The question is now which pushes to perform? Goldberg and Tarjan suggested to
put nodes of G (and hence Gf ) onto layers with t on the bottom-most layer and
to perform only pushes which transport excess to a lower layer. We use d(v) to
denote the (number of the) layer containing v. We call an edge e = (v, w) ∈ Gf
eligible if d(w) < d(v).
Let us summarize: a push across an edge e = (v, w) ∈ Gf can be performed if
v is active and e is eligible. It moves δ ≤ min(excess(v), cfe ) units of flow from v
to w. If e is also an edge of G then f (e) is increased by δ, and if e is the reversal
of an edge of G then f (e) is decreased by δ.
What are we going to do when v is active but there is no eligible edge out of v?
In this situation v is relabeled by increasing d(v) by one, i.e. the purpose of the
relabel operation is to create at least one edge on which the algorithm can perform
further pushes. The generic preflow-push algorithm is proposed in Figure 5.13.
Procedure genericPreflowPush(G=(V,E), f)
forall e = (s, v) ∈ E do push(e, cap(e)) // saturate
d(s):= n
d(v):= 0 for all other nodes
while ∃v ∈ V \ {s, t} : excess(v) > 0 do // active node
if ∃e = (v, w) ∈ Ef : d(w) < d(v) then // eligible edge
choose some δ ≤ min excess(v), cfe
push(e, δ) // no new steep edges
else d(v)++ // relabel. No new steep edges
Figure 5.13: The generic preflow-push algorithm.
The obvious choice for δ is min excess(v), cfe . We call an edge e = (v, w) ∈

Gf steep if d(w) < d(v) − 1, i.e., if it reaches down by two or more levels.
It might be instructive to visualize the generic preflow-push algorithm in terms
of a physical network: edges represent flexible water pipes, nodes represent joints
and the distance function measures how far nodes are above the ground. In this
network we wish to send water from the source to the sink. In addition, we visualize
pushed flow as water flowing downhill. Initially, we move the source node upward,
and water flows to its neighbors. In general, water flows downhill towards the sink;
however, occasionally flow becomes trapped locally at a node that has no downhill
neighbors. At this point we move the node upward, and again water flows downhill
toward the sink. As we continue to move nodes upward, the remaining excess
eventually flows back toward the source. The algorithm terminates when all the
water flows either into the sink or flows back to the source.
In the following we show that an arbitrary preflow-push algorithm finds a
maximum flow in time O(n2 m). We divide this into two steps. The first step shows
that if the algorithm terminates then it computed a maximum flow. Then we show
the running time of the algorithm.
Lemma 5
The algorithm maintains a preflow and does not generate steep edges. The nodes s
and t stay on levels 0 and n, respectively.
Proof.
The algorithm clearly maintains a preflow. After the initialization, each edge in
Gf either connects two nodes on level zero or connects a node on level zero to
node s on level n. Thus, there are no steep edges (there are not even any eligible
edges). A relabeling of a node v does not create a steep edge since a node is
only relabeled if there are no eligible edges out of it. A push across an edge
e = (v, w) ∈ Gf may add the edge (w, v) ∈ Gf . However, this edge is not
even eligible. Only active nodes are relabeled and only nodes different from s
and t can be active. Thus s and t stay on layers n and 0, respectively. The
preceding lemma has an interesting interpretation. Since there are no steep edges,
any path from v to t must have length (=number of edges) at least d(v) and any
path from v to s must have length at least d(v) − n. Thus, d(v) is a lower bound
on the distance from v to t and d(v)−n is a lower bound on the distance from v to s.
The next lemma shows that active nodes can always reach s in the residual
network (since they must be able to send their excess back to s). It has the important
consequence that d-labels are bounded by 2n − 1.
Lemma 6
If v is active then there is path from v to s in Gf :
∀ active nodes v : excess(v) > 0 ⇒ ∃ path hv, . . . , si ∈ Gf .
Proof.
The following is an equivalent invariant: There is an s–v path with nonzero flow to
any active node. Suppose now that push((u, w), δ) would destroy this invariant
for a path hs, . . . , w, u, . . . , vi. Then the path hs, . . . , u, . . . , vi reestablishes it.
s w u v
Lemma 7
No distance label ever reaches 2n (∀v ∈ V : d(v) < 2n).
Proof.
Assume that a node v is moved to level 2n = d(v). Since only active nodes are
relabeled, the lemma above implies the existence of a path p (and hence simple
path) in Gf from a node on level 2n to s (which is on level n). This path has at
most n − 1 nodes. We have d(s) = n and since the algorithm does not generate
steep edges d(v) is lower than 2n which is a contradiction.
Lemma 8
When genericPreflowPush terminates f is a maximum flow.
Proof.
When the algorithm terminates, there are no active nodes, i.e. ∀v ∈ V \ {s, t} :
excess(v) = 0, and hence the algorithm terminates with a flow. Call it f . In Gf
there can be no path from s to t since any such path must contain a steep edge
(since s is on level n, t is on level 0). Thus, f is a maximum flow (Max-Flow
Min-Cut Theorem).
Lemma 9
There are at most 2n2 relabels.
Proof.
No distance label ever exceeds 2n (d(v) ≤ 2n). That means v is relabeled at most
2n times. Hence, we have at most |V | · 2n = 2n2 relabel operations.
Lemma 10
There at most nm saturating pushes.
Proof.
A saturating push across an edge e = (v, w) ∈ Gf removes e from Gf . We claim
that v has to be relabeled at least twice before the next push across e and hence
there can be at most n saturating pushes across any edge. To see the claim, observe
that only a push across erev can again add e to Gf . Since pushes occur only across
eligible edges, w must be relabeled at least twice after the saturating push across e
and before the next push across erev . Similarly, it takes two relabels of v before e
becomes eligible again.
w
v lift
w
It is more difficult to bound the number of non-saturating pushes. It depends
heavily on which active node is selected for pushing and on which edge is selected
for pushing. Recall that a non-saturating push reduces excess(v) to zero, i.e.
deactivates a node.
Lemma 11
The number of non-saturating pushes is O(n2 m) if δ = min excess(v), cfe for

arbitrary node and edge selection rules. ( arbitrary preflow-push).

Proof.
The proof makes use of a potential function argument. Consider the potential
function
X
Φ := d(v)
{v:v is active}
We will show:
1. φ ≥ 0, and φ = 0 initially.
2. A non-saturating push decreases φ by at least one.
3. A relabeling increases φ by one.
4. A saturating push increases φ by at most 2n.

Suppose that we have shown (1) to (4). By (3) and (4) the total increase of φ is
at most 2n2 + nm2n, since there are at most 2n2 relabel operations and at most
nm saturating pushes. By (1) the total decrease can be no larger than this. Thus,
the number of non-saturating pushes can be at most O(n2 m).
Operation ∆(Φ) How many times? Total effect
relabel 1 ≤ 2n2 ≤ 2n2
saturating push ≤ 2n ≤ nm ≤ 2n2 m
non-saturating push ≤ −1
It remains to show (1) to (4). (1) is obvious. For (2) we observe that a non-saturating
push deactivates a node. It may or may not activate a node at the level below. In
either case, φ decreases by at least one. For (3) we observe that a relabeling of v
increases d(v) by one and for (4) we observe that a saturating push may activate a
node and that all distance labels are bounded by 2n.
Pushing Excess Out of a Node: Let v be a node with positive excess. We want
to push flow out of v along eligible edges. An edge e ∈ Gf is either also an edge
of G (and then f (e) < ce ) or the reversal of an edge of G (and then f (erev ) > 0)).
We therefore iterate over all edges out of v and all edges into v. For each edge out
of v we push max(excess(v), ce − fe ). If a push decreases the excess of v to zero
we break from the loop. However this implementation yields a degree-factor in the
running time which is easily removed by means of the so called current edge data
structure. We maintain for each node v an edge pointer currentEdge to its sequence
of outgoing edges in Gf such that all edges preceding this pointer in v adjacency
list are not eligible. Recall that an edge (v, w) is eligible if the layer of w is one less
than the layer of v and that no edge goes down more that one layer. Thus relabeling
w cannot make (v, w) eligible, and reversing an edge in an augmentation cannot
make the edge eligible (because all edges in the augmenting path go from lower
layers to higher layers after the augmentation). Only relabeling v can make an
edge out of v eligible. With these observations it is easy to maintain the invariant
that all edges preceding currentEdge in v’s adjacency list are not eligible: When v
is relabel we reset the currentEdge pointer. This yields the following Lemma.
Lemma 12 X
Total cost for searching ≤ 2n · degree(v) = 4nm = O(nm)
v∈V
Theorem 13
The arbitrary preflow-push algorithm finds a maximum flow in time O(n2 m).
Proof.
Lemma 8 gives us partial correctness. The time for initialization is bounded through
O(n + m). To maintain a set of active nodes we use for example a stack or a FIFO
and we use reverse edge pointers to implement push. Lemma 9 bounds the number
of relabel operations to 2n2 . The saturating pushes are less or equal than nm
(Lemma 10) and the non-saturating pushes are in O(n2 m) (Lemma 11). Lemma 12
states the overall search time for eligible edges is in O(nm). Summing up yields a
total time which is in O(n2 m). There are multiple ways to improve this running
time. The first one is the FIFO-rule: The active nodes are kept in a queue and the
first node in the queue is always selected. When a node is relabeled or activated
the node is added to the rear of the queue. The number of non-saturating pushes is
O(n3 ) when the FIFO-rule is uses. This bound is due to Goldberg.
Another variant follows the Highest Level-rule: An active node on the highest
level, i.e., with the maximal dist-value, is selected. Observe that when a maximal
level active node is relabeled, it will be the unique maximal active node after
the relabel. Thus, this rule guarantees that when a node is selected, pushes out
of the node will be performed√until the node becomes inactive. The number of
non-saturating pushes is O(n2 m) when the Highest Level-rule is uses. This rule
can be implemented using a bucket priority queue.
Lemma 14
When√the Highest Level-rule is used, the number of non-saturating pushes is
O(n2 m).
Proof. √
Again we use a potential function argument. Let K = m; this choice of K will
become clear at the end of the proof. For a node v, let
|{w : d(w) ≤ d(v)}|

d0 (v) :=
K
which basically is the scaled number of dominated nodes. We consider
X
Φ := d0 (v)
{v:v is active}
as potential function. We split the execution into phases. We define a phase to

consist of all pushes between two consecutive changes of
d∗ := max {d(v) : v is active}
and call a phase expensive if it contains more than K non-saturating pushes, and
cheap otherwise. We show:
1. ≤ 4n2 phases and ≤ 4n2 K non-saturating pushes in all cheap phases to-
gether
2. Φ ≥ 0 always, Φ ≤ n2 /K initially (obvious)
3. a relabel or saturating push increases Φ by at most n/K.
4. a non-saturating push does not increase Φ.
5. an expensive phase with Q ≥ K non-saturating pushes decreases Φ by at

least Q.
Suppose that we haven shown (1) to (5). (3) and (4) imply that the total increase
of φ is at most (2n2 + mn)n/K and hence the total decrease can be at most this
number plus n2 /K by (2). The number of non-saturating pushes in expensive
phases is therefore bounded by (2n3 + n2 + mn2 )/K. Together with (1) we
conclude that the total number of non-saturating pushes is at most
(2n3 + n2 + mn2 )/K + 4n2 K.

√
Observing that n = O(m) and that the choice K = m balances √ the contributions
2
from expensive and cheap phases, we obtain a bound of O(n m).
It remains to prove (1) to (5). For we observe that d∗ = 0 initially, d∗ ≥ 0

always, and that only a relabel can increase d∗ . Thus, d∗ is increased at most 2n2
times, decreased no more than this, and hence changed at most 4n2 . Therefore the
number of phases is bounded through 4n2 . This directly yields that the number
of non-saturating pushes in all cheap phases together is less or equal than 4n2 K.
(2) is obvious. (3) follows from the observation that d0 (v) ≤ n/K for all v and
at all times. A relabel of v can increase only the d0 -value of v. A saturating
push on (u, w) may activate only w. For (4) observe that a non-saturating push
across an edge (v, u) deactivates v, activates u (if it is not already active), and
that d0 (u) ≤ d0 (v) (we do not push flow away from the sink). For (5) consider
an expensive phase containing Q ≥ K non-saturating pushes. By definition of a
phase, d∗ is constant during a phase, and hence all Q non-saturating pushes must
be out of nodes at level d∗ . The phase is finished either because level d∗ becomes
empty or because a node is moved from level d∗ to level d∗ + 1. In either case,
we conclude that level d∗ contains Q ≥ K nodes at all times during the phase.
Thus, each non-saturating push in the phase decreases φ by at least one (since
d0 (u) ≤ d0 (v) − 1 for a push from v to u).
5.3.1 Heuristic Improvements
What is the best case running time of the algorithm? The running time is Ω(n2 )
if Ω(n) nodes need to be lifted above level n. This is usually the case. The best
case behavior of the other parts of the algorithm is O(m) and hence the cost of
relabeling dominates the best case running time. In this section we will describe
several heuristics that frequently reduce the time spent in relabeling nodes and as a
side-effect reduce the time spent in all other operations. The heuristics will turn the
preflow-push algorithm into a highly effective algorithm for solving flow problems.
The Local Relabeling Heuristic: The local relabeling heuristic applies when-
ever a node is relabeled. It increases the dist-value of v to
d(v):= 1 + min {d(w) : (v, w) ∈ Gf } .
Observe that v is active whenever it is relabeled and that an active node has at least
one outgoing edge in Gf . The expression above is therefore well defined. When v
is relabeled non of the outgoing edges is eligible and hence d(w) ≥ d(v) for all
(v, w) ∈ Gf . Thus the local relabeling heuristic increases d(v) by at least one. It
may increase it by more than one. The correctness of the heuristic follows from
the following alternative description: when a node is relabeled, continue to relabel
it until there is an eligible edge out of it.
The Global Relabeling Heuristic: The global relabeling heuristic updates the
distance values of all nodes. It sets d(v) := Gf .reverseBFS(t) for nodes that can
reached in Gf (nodes that aren’t reachable get another d-value which is not defined
here but in the LEDA book of Mehlhorn/Näher). A breadth-first search requires
time O(m). It should therefore not be applied to frequently. We will apply it
initially and every O(m) edge inspections. In this way Ω(m) time is spent between
applications of the global relabel heuristic and hence the worst case running time
is increased by at most a constant factor. The best case can improve significantly.
10 10 10 10 10
10 12 10 12
10 0 0 10 1 0
4 4
s t s t
6 4 0 6 4 0
d d
cap 10 0 0 8 cap 10 0 0 8
f 10 f 10
excess 10 4 excess 10 4
10 10 10 10
10 12 10 12
10 1 1 10 1 1
4 10 4 10
s t s t
6 4 0 6 4 0
d d
cap 10 0 0 8 cap 10 1 0 8
f 10 f 10
excess 10 4 excess 6 4 4 4
10 10 4 10 10 4
10 12 10 12
10 1 1 10 1 1
4 10 4 10
s t s t
4 0 4 0
6 4 6 4
d d
cap 10 2 0 8 cap 10 7 0 8
f 10 f 8
excess 2 4 4 4 excess 4 4 4
10 10 4 10 10 2
10 12 10 12
10 1 1 10 1 1
4 10 4 12
s t s t
4 0 4 0
6 4 6 4
d d
cap 10
8 7 1 48 cap 10
8 7 1 48
f f
excess 4 4 excess 4 4
2 10 8 10 8
10 10 12
10 1 2 12 10 2 2
4 12 4 12
2
s t s t
4 4 0
6 4 0 6 4
d d
cap 10
8 7 1 48 cap 10
8 7 1 48
f f
excess 4 4 excess 4 4 2
10 8
10 12
10 2 2
4 12
2
s t
4 0
6 4
d
cap 10
8 7 1 68
f
excess 4 4
Figure 5.14: An example run of the arbitrary preflow-push algorithm (12 pushes).
The Gap Heuristic: We come to our last optimization. Consider a relabeling

of a node v in phase one and let dv be the layer of v before the relabeling. If the
layer dv becomes empty by the relabeling of v, then v cannot reach t anymore in
Gf after the relabeling, since any edge crossing the now empty layer would be
steep. If v cannot reach t in Gf then no node reachable from v in Gf can reach t.
We may therefore move v and all nodes reachable from v to layer n whenever the
old layer of v becomes empty by relabeling of v. This is called the gap heuristic.
5.3.2 Experimental Results

In this section we present experimental results for the different algorithms and
configurations we have seen. We use four classes of graphs:
• Random: n nodes, 2n + m edges; all edges (s, v) and (v, t) exist
• Cherkassky and Goldberg (1997) (two graph classes)
• Ahuja, Magnanti, Orlin (1993)
Timings: Random Graphs
Rule BASIC HL LRH GRH GAP LEDA

FF 5.84 6.02 4.75 0.07 0.07 —
33.32 33.88 26.63 0.16 0.17 —
HL 6.12 6.3 4.97 0.41 0.11 0.07
27.03 27.61 22.22 1.14 0.22 0.16
MF 5.36 5.51 4.57 0.06 0.07 —
26.35 27.16 23.65 0.19 0.16 —
Timings: CG1

FF 3.46 3.62 2.87 0.9 1.01 —
15.44 16.08 12.63 3.64 4.07 —
HL 20.43 20.61 20.51 1.19 1.33 0.8
192.8 191.5 193.7 4.87 5.34 3.28
MF 3.01 3.16 2.3 0.89 1.01 —
12.22 12.91 9.52 3.65 4.12 —
Timings: CG2

FF 50.06 47.12 37.58 1.76 1.96 —
239 222.4 177.1 7.18 8 —
HL 42.95 41.5 30.1 0.17 0.14 0.08
173.9 167.9 120.5 0.36 0.28 0.18
MF 45.34 42.73 37.6 0.94 1.07 —
198.2 186.8 165.7 4.11 4.55 —
Timings: AMO

FF 12.61 13.25 1.17 0.06 0.06 —
55.74 58.31 5.01 0.1399 0.1301 —
HL 15.14 15.8 1.49 0.13 0.13 0.07
62.15 65.3 6.99 0.26 0.26 0.14
MF 10.97 11.65 0.04999 0.06 0.06 —
46.74 49.48 0.1099 0.1301 0.1399 —
FF=FIFO node selection, HL=hightest level, MF=modified FIFO

HL= d(v) ≥ n is special,
LRH=local relabeling heuristic, GRH=global relabeling heuristics
Asymptotics, n ∈ {5000, 10000, 20000}

Gen Rule GRH GAP LEDA
rand FF 0.16 0.41 1.16 0.15 0.42 1.05 — — —
HL 1.47 4.67 18.81 0.23 0.57 1.38 0.16 0.45 1.09
MF 0.17 0.36 1.06 0.14 0.37 0.92 — — —
CG1 FF 3.6 16.06 69.3 3.62 16.97 71.29 — — —
HL 4.27 20.4 77.5 4.6 20.54 80.99 2.64 12.13 48.52
MF 3.55 15.97 68.45 3.66 16.5 70.23 — — —
CG2 FF 6.8 29.12 125.3 7.04 29.5 127.6 — — —
HL 0.33 0.65 1.36 0.26 0.52 1.05 0.15 0.3 0.63
MF 3.86 15.96 68.42 3.9 16.14 70.07 — — —
AMO FF 0.12 0.22 0.48 0.11 0.24 0.49 — — —
HL 0.25 0.48 0.99 0.24 0.48 0.99 0.12 0.24 0.52
MF 0.11 0.24 0.5 0.11 0.24 0.48 — — —
5.4 Minimum Cost Flows

The idea of augmenting paths extends to a more general network flow problem.
We define G = (V, E), f , excess, and c as for maximum flows. Let cost : E → R
denote the edge costs Pper unit of flow in addition to a capacity. Consider further
supply : V → R with v∈V supply(v) P = 0. A negative supply is called a demand.
The cost of a flow f is cost(f ) := e∈E f (e)cost(e). A flow is minimum cost if
among all flows of the same value it has minimum costs. The minimum cost flow
problem is that of finding a maximum flow of minimum cost. The min-cost flow
problem is as follows:
X
minimize c(f ) := f (e)cost(e)
e∈E
(
∀v ∈ V : excess(v) = −supply(v) // flow conservation constraints
subject to:
∀e ∈ E : f (e) ≤ c(e) // capacity constraints
Cycle Canceling Algorithm for Min-Cost Flow

To propose the algorithm we first need to define residual cost: For e ∈ E,
we set costf (e) = cost(e) (edge exists in Gf iff f (e) < c(e)), and we set
costf (reverse(e)) = −cost(e) (edge exists in Gf iff f (e) > 0).
Lemma 15 P
A feasible flow is optimal iff 6 ∃ cycle C ∈ Gf : costf (C) = e∈C costf (e) < 0
Proof.
not here (see [3], p. 307)
We now propose the algorithm which has pseudo-polynomial running time:
f := any feasible flow

invariant f is feasible
while ∃ cycle C : costf (C) < 0 do
augment flow around C
Figure 5.15: The cycle canceling algorithm.
But how do we compute a feasible flow? We can do this using maximum

flow algorithms and the following construction: First we set up a maximum flow
network G∗ starting with the min cost flow problem G:
5.4. MINIMUM COST FLOWS 39
• Add a vertex s
• ∀v ∈ V with supply(v) > 0, add edge (s, v) with cap. supply(v)
• Add a vertex t
• ∀v ∈ V with supply(v) < 0, add edge (v, t) with cap. −supply(v)
• find a maximum flow f in G∗
If f saturates the edges leaving s then f is feasible for G. Otherwise there cannot
be a feasible flow f 0 because f 0 could easily be converted into a flow in G∗ with
larger value.
Lemma 16 (Integrality Property)

If all edge capacities are integral then there exists an integral minimum cost flow.
Proof.
Exercise.
There are better algorithms to solve the problem:
Theorem 17
The min-cost flow problem can be solved in time O(mn log n + m2 log maxe∈E c(e)).
Proof.
For details take the courses in optimization or network flows.
Special Cases of Min Cost Flows

Transportation Problem: The transportation problem is a special case of the
minimum cost flow problem with the property that the node set V is partitioned
into two subsets V1 and V2 (of possibly unequal cardinality) so that (1) each node
in V1 is a supply node, (2) each node in V2 is a demand node, (3) each edge starts
in V1 and ends in V2 . We further have ∀e ∈ E : c(e) = ∞. The classical example
of this problem is the distribution of goods from warehouses to customers. In
this context the nodes in V1 represent the warehouses, the nodes in V2 represent
customers (or, more typically, customer zones), and an edge (u.v) represents a
distribution channel from warehouse u to customer v.
Minimum Cost Bipartite Perfect Matching: The transportation problem can
be used to solve the Minimum Cost Bipartite Perfect Matching. We therefore
formulate a transportation problem in a bipartite graph G = (A ∪ B, E ⊆ A × B)
with supply(v) = 1 for v ∈ A, supply(v) = −1 for v ∈ B. An integral flow
defines a matching. Reminder: M ⊆ E is a matching if (V, M ) has maximum
degree one.
5.5 Matching Problems

A matching in a graph G = (V, E) is a set of edges with the property that every
node is incident to at most one edge in the set; thus a matching induces a pairing of
(some of) the nodes in the graph using the edges in E. In a matching, each node
is matched with at most one other node, and some nodes might not be matched
with any other node. The matching problem seeks a matching that optimizes
some criteria. Matching problems on bipartite graphs are called bipartite matching
problems, and those on non-bipartite graphs are called nonbipartite matching
problem . There are two additional ways of categorizing matching problems:
cardinality matching problems, which maximize the number of pairs of nodes
matched, and weighted matching problems which maximize or minimize the
∗ ∗
P
weight of the matching (Find a matching M such that w(M ) := e∈M ∗ w(e)
is maximized). Applications of matching problems arise in matching roommates
to hostels, matching pilots to compatible airplanes, scheduling airline crews for
available flight legs, and assigning duties to bus drivers.
Theorem 18
A maximum weighted matching can be found in time O(nm + n2 log n). [Gabow
1992]
Proof.
Not here.
Theorem 19
There is an O(m) time algorithm that finds a matching of weight at least maxmatching M w(M )/2.
The algorithm is a 1/2-approximation algorithm [18].
Proof.
The algorithm is as follows:
M 0 := ∅
invariant M 0 is a set of simple paths
while E 6= ∅ do // find heavy simple paths
select any v ∈ V with degree(v) > 0 // select a starting node
while degree(v) > 0 do // extend path greedily
(v, w):= heaviest edge leaving v // (*)
0 0
M := M ∪ {(v, w)}
remove v from the graph
v:= w
return any matching M ⊆ M 0 with w(M ) ≥ w(M 0 )/2
// one path at a time, e.g., look at the two ways to take every other edge.
5.5. MATCHING PROBLEMS 41
We now have to proof the approximation ratio: Let M ∗ denote a maximum weight
matching. It suffices to show that w(M 0 ) ≥ w(M ∗ ). Therefore we assign each
edge to that incident node that is deleted first. Now all e∗ in our maximum weight
matching M ∗ are assigned to different nodes (otherwise it wouldn’t be a matching).
Consider any edge e∗ ∈ M ∗ and assume it is assigned to node v. Since e∗ is
assigned to v, it was available in line (*), i.e. it could have been selected in this
line. Hence, there is an edge e ∈ M 0 assigned to v with w(e) ≥ w(e∗ ).
Lemma 20
There is an O m log 1 time algorithm that computes a 23 − -approximation of the

maximum weight matching problem [43].

Proof.
The algorithm is as follows.
M := ∅
repeat Θ n log 1 times:

Let X ∈ V (G) be selected uniformly at random

M := M ⊕ aug(X)
return M
For more details look into the paper.
5.5.1 Maximum Cardinality Bipartite Matching

As introduced in the beginning of this Section a maximum cardinality matching is
a matching of maximum size. We now consider the problem of finding a maximum
cardinality matching in a given bipartite graph G = (L ∪ R, E). It turns out that
this problem can be reduced to a unit network max-flow problem. The construction
is as follows:
G0 = ({s} ∪ L ∪ R ∪ {t} , {(s, u) : u ∈ L} ∪ EL→R ∪ {(v, t) : v ∈ R}).
L R
s t
Here EL→R are all edges E directed from L to R. The capacity of every edge
is set to one. First notice that if G has a matching of size k then G0 has a flow of
value k. To see this let M be a matching of size k. Now we construct a flow that
sends one unit along each edge in M . This flow has value k.
On the other hand if G0 has a flow of value k then G has a matching of size k.
To prove this consider a flow f of value k. First observe that f is an integral flow.
Thus each edge has flow 1 or 0. Consider the set M of edges from L to R that
have flow 1. M has k edges because the value of the flow is equal to the number of
non-zero flow edges crossing the cut (L ∪ {s}, R ∪ {t}). Furthermore, each vertex
has at most one edge in M incident upon it (since all edges have capacity 1).
Thus, to find a maximum cardinality matching in G, we construct G0 , find
the maximum flow in G0 and construct the matching M as above. To solve the
0
unit network√ max-flow problem G we can use Dinitz algorithm which yields a
O((n + m) n) algorithm for the maximum cardinality bipartite matching prob-
lem.
5.6 Book References

This chapter is largely based on Mehlhorn et. al. [37] as well as Ahuja et. al. [3].
Chapter 6
Generic Approaches to Optimization
Please refer to Chapter 12 of the “Data Structures and Algorithms – The Basic
Toolbox” book as well as Chapter 15, 16 and 35 of “Introduction to Algorithms”
by Cormen, Leiserson, Rivest and Stein. These sections include dynamic program-
ming as well as greedy and branch-and-bound algorithms.
43
44 CHAPTER 6. GENERIC APPROACHES TO OPTIMIZATION
Chapter 7
Linear Programming
Please refer to Chapter 12 of the “Data Structures and Algorithms – The Basic
Toolbox” book as well as Chapter 29 of “Introduction to Algorithms” by Cormen,
Leiserson, Rivest and Stein.
45
46 CHAPTER 7. LINEAR PROGRAMMING
Chapter 8
Parallel Algorithms
This chapter is mainly based on [29].
8.1 Introduction
Numerous programming languages and libraries have been developed for explicit
parallel programming. These differ in their view of the address space that they
make available to the programmer, the degree of synchronization imposed on
concurrent activities, and the multiplicity of programs. The message-passing
programming paradigm is one of the oldest and most widely used approaches for
programming parallel computers. Its roots can be traced back in the early days of
parallel processing, and its wide-spread adoption can be attributed to the fact that
it imposes minimal requirements on the underlying hardware.
8.2 Message-Passing Model/Paradigm

There are two key attributes that characterize the message-passing programming
paradigm. The first is that it assumes a partitioned address space, and the second is
that it supports only explicit parallelization.
The logical view of a machine supporting the message-passing paradigm con-
sists of p processes, each with its own exclusive address space. Instances of such
a view come naturally from clustered workstations and non-shared address space
multicomputers. There are two immediate implications of a partitioned address
space. First, each data element must belong to one of the partitions of the space;
hence, data must be explicitly partitioned and placed. This adds complexity to
programming, but encourages locality of access, which is critical for achieving
high performance on non-UMA architecture, since a processor can access its local
data much faster than non-local data on such architectures. The second implication
47
48 CHAPTER 8. PARALLEL ALGORITHMS
is that all interactions (read-only or read/write) require cooperation of two pro-

cesses, the process that has the data and the process that wants to access the data.
This requirement for cooperation adds a great deal of complexity for a number
of reasons. The process that has the data must participate in the interaction even
if it has no logical connection to the events at the requesting process. In certain
circumstances, this requirement leads to unnatural programs. In particular, for dy-
namic and/or unstructured interactions, the complexity of the code written for this
type of paradigm can be very high for this reason. However, a primary advantage
of explicit two-way interactions is that the programmer is fully aware of all the
costs of non-local interactions, and is more likely to think about algorithms (and
mappings) that minimize interactions. Another major advantage of this type of
programming paradigm is that it can be efficiently implemented on a wide variety
of architectures.
The message-passing programming paradigm requires that the parallelism is
coded explicitly by the programmer. That is, the programmer is responsible for
analyzing the underlying serial algorithm/application and identifying ways by
which he or she can decompose the computations and extract concurrency. As a
result, programming using the message-passing paradigm tends to be hard and
intellectually demanding. However, on the other hand, properly written message-
passing programs can often achieve very high performance and scale to a very
large number of processes.
Message-passing programs are often written using the asynchronous or loosely
synchronous paradigms. In the asynchronous paradigm, all concurrent tasks exe-
cute asynchronously. This makes it possible to implement any parallel algorithm.
However, such programs can be harder to reason about, and can have nonde-
terministic behavior due to race conditions. Loosely synchronous programs are
a good compromise between these two extremes. In such programs, tasks or
subsets of tasks synchronize to perform interactions. However, between these
interactions, tasks execute completely asynchronously. Since the interaction hap-
pens synchronously, it is still quite easy to reason about the program. Many of the
known parallel algorithms can be naturally implemented using loosely synchronous
programs.
In its most general form, the message-passing paradigm supports execution of a
different program on each of the p processes. This provides the ultimate flexibility
in parallel programming, but makes the job of writing parallel programs effectively
unscalable. For this reason, most message-passing programs are written using the
single program multiple data (SPMD) approach. In SPMD programs the code
executed by different processes is identical except for a small number of processes
(e.g., the “root” process). This does not mean that the processes work in lock-step.
In an extreme case, even in an SPMD program, each process could execute a
different code (the program contains a large case statement with code for each
8.3. COMMUNICATION COST IN PARALLEL MACHINES 49
process). But except for this degenerate case, most processes execute the same
code. SPMD programs can be loosely synchronous or completely asynchronous.
8.3 Communication Cost in Parallel Machines

The time taken to communicate a message between two nodes in a network is the
sum of the time to prepare a message for transmission and the time taken by the
message to traverse the network to its destination. The principal parameters that
determine the communication latency are as follows:
• Startup time Tstart : The startup time is the time required to handle a message
at the sending and receiving nodes. This includes the time to prepare the
message (adding header, trailer, and error correction information), the time to
execute the routing algorithm, and the time to establish an interface between
the local node and the router. This delay is incurred only once for a single
message transfer.
• Per-hop time Thop : After a message leaves a node, it takes a finite amount
of time to reach the next node in its path. The time taken by the header of a
message to travel between two directly-connected nodes in the network is
called the per-hop time. It is also known as node latency. The per-hop time
is directly related to the latency within the routing switch for determining
which output buffer or channel the message should be forwarded to.
• Per-byte transfer time Tbyte : If the channel bandwidth is r bytes per second,
then each word takes time Tbyte = 1/r to traverse the link. This time is
called the per-byte transfer time. This time includes network as well as
buffering overheads.
This implies that in order to optimize the cost of message transfers, we would need
to:
• Communicate in bulk. That is, instead of sending small messages and paying
a startup costs for each one, we want to aggregate small messages into a
single large message and amortize the startup latency across a larger message.
This is because on typical platforms such as clusters and message-passing
machines, the value of Tstart is much larger than those of Thop or Tbyte .
• Minimize the volume of data. To minimize the overhead paid in terms of

per-byte transfer time Tbyte , it is desirable to reduce the volume of data
communicated as much as possible.
• Minimize distance of data transfer. Minimize the number of hops that a

message must traverse.
While the first two objectives are relatively easy to achieve, the task of minimiz-
ing distance of communicating nodes is difficult, and in many cases an unnecessary
burden on the algorithm designer. This is a direct consequence of the following
characteristics of parallel platforms and paradigms:
• In many message-passing libraries such as MPI, the programmer has lit-

tle control on the mapping of processes onto physical processors. In such
paradigms, while tasks might have well defined topologies and may commu-
nicate only among neighbors in the task topology, the mapping of processes
to nodes might destroy this structure.
• Many architectures rely on randomized (two-step) routing, in which a mes-

sage is first sent to a random node from source and from this intermediate
node to the destination. This alleviates hot-spots and contention on the net-
work. Minimizing number of hops in a randomized routing network yields
no benefits.
• The per-hop time Thop is typically dominated either by the startup latency
Tstart for small messages or by per-byte component `Tbyte for large messages
`. Since the maximum number of hops in most networks is relatively small,
the per-hop time can be ignored with little loss in accuracy.
All of these point to a simpler cost model in which the cost of transferring a
message between two nodes on a network is given by:
Tcomm (`) = Tstart + `Tbyte
This expression has significant implications for architecture-independent al-

gorithm design as well as for the accuracy of runtime predictions. Since this cost
model implies that it takes the same amount of time to communicate between any
pair of nodes, it corresponds to a completely connected network. Instead of de-
signing algorithms for each specific architecture (for example, a mesh, hypercube,
or tree), we can design algorithms with this cost model in mind and port it to any
target parallel computer.
This raises the important issue of loss of accuracy (or fidelity) of predic-
tion when the algorithm is ported from our simplified model (which assumes a
completely connected network) to an actual machine architecture. If our initial
assumption that the Thop term is typically dominated by the Tstart or Tbyte terms is
valid, then the loss in accuracy should be minimal.
8.4. PERFORMANCE METRICS 51
However, it is important to note that our basic cost model is valid only for
uncongested networks. Architectures have varying thresholds for when they get
congested; i.e., a linear array has a much lower threshold for congestion than a
hypercube. Furthermore, different communication patterns congest a given network
to different extents. Consequently, our simplified cost model is valid only as long
as the underlying communication pattern does not congest the network.
8.4 Performance metrics

It is important to study the performance of parallel programs with a view to
determining the best algorithm, evaluating hardware platforms, and examining
the benefits from parallelism. A number of metrics have been used based on the
desired outcome of performance analysis.
8.4.1 Execution Time

The serial runtime of a program is the time elapsed between the beginning and
the end of its execution on a sequential computer. The parallel runtime is the time
that elapses from the moment a parallel computation starts to the moment the last
processing element finishes execution. We denote the serial runtime by Tseq and
the parallel runtime by T (p).
8.4.2 Speedup
When evaluating a parallel system, we are often interested in knowing how much
performance gain is achieved by parallelizing a given application over a sequential
implementation. Speedup is a measure that captures the relative benefit of solving
a problem in parallel. It is defined as the ratio of the time taken to solve a problem
on a single processing element to the time required to solve the same problem on a
parallel computer with p identical processing elements. We denote speedup by the
symbol S.
For a given problem, more than one sequential algorithm may be available, but
all of these may not be equally suitable for parallelization. When a serial computer
is used, it is natural to use the sequential algorithm that solves the problem in the
least amount of time. Given a parallel algorithm, it is fair to judge its performance
with respect to the fastest sequential algorithm for solving the same problem on
a single processing element. Sometimes, the asymptotically fastest sequential
algorithm to solve a problem is not known, or its runtime has a large constant
that makes it impractical to implement. In such cases, we take the fastest known
algorithm that would be a practical choice for a serial computer to be the best
sequential algorithm. We compare the performance of a parallel algorithm to solve

a problem with that of the best sequential algorithm to solve the same problem. We
formally define the speedup S as the ratio of the serial runtime of the best sequential
algorithm for solving a problem to the time taken by the parallel algorithm to solve
the same problem on p processing elements. The p processing elements used by
the parallel algorithm are assumed to be identical to the one used by the sequential
algorithm.
8.4.3 Superlinear speedup

Theoretically, speedup can never exceed the number of processing elements, p. If
the best sequential algorithm takes Tseq units of time to solve a given problem on a
single processing element, then a speedup of p can be obtained on p processing
elements if none of the processing elements spends more than time Tseq /p. A
speedup greater than p is possible only if each processing element spends less
than time Tseq /p solving the problem. In this case, a single processing element
could emulate the p processing elements and solve the problem in fewer than Tseq
units of time. This is a contradiction because speedup, by definition, is computed
with respect to the best sequential algorithm. If Tseq is the serial runtime of the
algorithm, then the problem cannot be solved in less than time Tseq on a single
processing element.
In practice, a speedup greater than p is sometimes observed (a phenomenon
known as superlinear speedup). This usually happens when the work performed
by a serial algorithm is greater than its parallel formulation or due to hardware
features that put the serial implementation at a disadvantage. For example, the data
for a problem might be too large to fit into the cache of a single processing element,
thereby degrading its performance due to the use of slower memory elements. But
when partitioned among several processing elements, the individual data-partitions
would be small enough to fit into their respective processing elements’ caches.
8.4.4 Efficiency
Only an ideal parallel system containing p processing elements can deliver a
speedup equal to p. In practice, ideal behavior is not achieved because while
executing a parallel algorithm, the processing elements cannot devote 100% of
their time to the computations of the algorithm. Efficiency is a measure of the
fraction of time for which a processing element is usefully employed; it is defined
as the ratio of speedup to the number of processing elements. In an ideal parallel
system, speedup is equal to p and efficiency is equal to one. In practice, speedup is
less than p and efficiency is between zero and one, depending on the effectiveness
8.4. PERFORMANCE METRICS 53
with which the processing elements are utilized. We denote efficiency by the
symbol E. Mathematically, it is given by
S
E=
p
8.4.5 Example: Adding n numbers using n processing elements

Consider the problem of adding n numbers by using p = n processing elements.
Initially, each processing element is assigned one of the numbers to be added and,
at the end of the computation, one of the processing elements stores the sum of all
the numbers. Assuming that n is a power of two, we can perform this operation
in log n steps by propagating partial sums up a logical binary tree of processing
elements. Figure 8.1 illustrates the procedure for n = 16. The processing elements
are labeled from 0 to 15. Similarly, the 16 numbers to be added are labeled from
Pj0
to 15. The sum of numbers with consecutive labels from i to j is denoted by i .
P Computing the global sum of 16 partial sums using 16 processing

Figure 8.1:
elements. ji denotes the sum of numbers with consecutive labels from i to j.
Each step shown in Figure 8.1 consists of one addition and the communication
of a single word. The addition can be performed in some constant time, say Tc0 ,
and the communication of a single word can be performed in time Tstart + Tbyte .
Therefore, the addition and communication operations take a constant amount of
time. Thus,
T (p) = Θ(log n)
Since the problem can be solved in Θ(n) time on a single processing element, its
speedup is
n
S = Θ( )
log n
From this equation and the preceding definition, the efficiency of the algorithm for
adding n numbers on n processing elements is
Θ( logn n ) 1
E= = Θ( )
n log n
8.4.6 Example: Adding n numbers using p processing elements
Figure 8.2: Computing the sum of 16 numbers using four processing elements.
An alternate method for adding n numbers using p processing elements is

illustrated in Figure 8.2 for n = 16 and p = 4. In the first step of this algorithm,
each processing element locally adds its n/p numbers in time Tseq = Θ(n/p). Now
the problem is reduced to adding the p partial sums on p processing elements, which
can be done in time Θ(log p) by the method described in the previous example.
The parallel runtime of this algorithm is
T (p) = Θ(n/p + log p)

8.5. PREFIX SUMS ON A HYPERCUBE NETWORK 55
the efficiency is

Tseq (n) 1 p log p
E= = =1−Θ
p(Tseq (n/p) + Θ(log p)) 1 + Θ(p log(p)) /n n
8.5 Prefix Sums on a Hypercube Network

In the following, we describe prefix sum operations and derive expressions for
their time complexity. We assume that the communication time between any
pair of nodes is practically independent of of the number of intermediate nodes
along the paths between them. We also assume that the communication links are
bidirectional; that is, two directly-connected nodes can send messages of size m
to each other simultaneously in time Tstart + `Tbyte . We assume a single-port
communication model, in which a node can send a message on only one of its links
at a time. Similarly, it can receive a message on only one link at a time. However,
a node can receive a message while sending another message at the same time on
the same or a different link.
Algorithm 1 PREFIX_SUMS_HCUBE(my_id, my_number, d, result)

result := my_number
msg := result
for i := 0 to d − 1 do
partner := my_id XOR 2i
send msg to partner
receive number from partner
msg := msg + number
if partner < my_id then
result := result + number
end if
end for
Algorithm 1 gives a procedure for implementing prefix sums on a d-dimensional

hypercube. Communication starts from the lowest dimension of the hypercube and
then proceeds along successively higher dimensions (Line 5). In each iteration,
nodes communicate in pairs so that the labels of the nodes communicating with
each other in the i th iteration differ in the ith least significant bit of their binary
representations (Line 6).The node with label k uses information from only the k-
node subset of those nodes whose labels are less than or equal to k. To accumulate
the correct prefix sum, every node maintains an additional result buffer. This buffer
is denoted by square brackets in Figure 8.3. After an iteration’s communication
steps, each node adds the data it receives during that iteration to its resident data
(denoted by parentheses in the figure) (Line 9). This sum is transmitted in the
following iteration. The content of an incoming message is added to the result
buffer only if the message comes from a node with a smaller label than that of the
recipient node (Line 10).
On a p-node hypercube, the size of each message exchanged in the i th of the
log p steps is `. It takes a pair of nodes time Tstart + `Tbyte to send and receive
messages from each other during the i th step. Hence, the time to complete the
entire procedure is
T = log p(Tstart + `Tbyte )
8.6 Quicksort
This section examines the quicksort algorithm, which has an average complexity of
Θ(n log n). Quicksort is one of the most common sorting algorithms for sequential
computers because of its simplicity, low overhead, and optimal average complexity.
8.6. QUICKSORT 57
8.6.1 Sequential Algorithm

Quicksort is a divide-and-conquer algorithm that sorts a sequence by recursively
dividing it into smaller subsequences. Assume that the n-element sequence to be
sorted is stored in the array A[1 . . . n]. Quicksort consists of two steps: divide and
conquer. During the divide step, a sequence A[q . . . r] is partitioned (rearranged)
into two nonempty subsequences A[q . . . s] and A[s+1 . . . r] such that each element
of the first subsequence is smaller than or equal to each element of the second
subsequence. During the conquer step, the subsequences are sorted by recursively
applying quicksort. Since the subsequences A[q . . . s] and A[s + 1 . . . r] are sorted
and the first subsequence has smaller elements than the second, the entire sequence
is sorted.
How is the sequence A[q . . . r] partitioned into two parts – one with all elements
smaller than the other? This is usually accomplished by selecting one element
x from A[q . . . r] and using this element to partition the sequence A[q . . . r] into
two parts – one with elements less than or equal to x and the other with elements
greater than x. Element x is called the pivot. The quicksort algorithm is presented
in Algorithm 2. This algorithm arbitrarily chooses the first element of the sequence
A[q . . . r] as the pivot. The operation of quicksort is illustrated in Figure 8.4.
Algorithm 2 QUICKSORT(A, q, r)
if q < r then
x := A[q]
s := q
for i := q + 1 to r do
if A[i] ≤ x then
s := s + 1
swap(A[s], A[i])
end if
end for
swap(A[q], A[s])
QUICKSORT(A, q, s)
QUICKSORT(A, s + 1, r)
end if
The complexity of partitioning a sequence of size k is Θ(k). Quicksort’s

performance is greatly affected by the way it partitions a sequence. Consider
the case in which a sequence of size k is split poorly, into two subsequences of
sizes 1 and k − 1. The run time in this case is given by the recurrence relation
T (n) = T (n − 1) + Θ(n), whose solution is T (n) = Θ(n2 ). Alternatively,
consider the case in which the sequence is split well, into two roughly equal-size
.
subsequences of and elements. In this case, the run time is given by the recurrence
relation T (n) = 2T (n/2) + Θ(n), whose solution is T (n) = Θ(n log n). The
second split yields an optimal algorithm. Although quicksort can have O(n2 ) worst-
case complexity, its average complexity is significantly better; the average number
of compare-exchange operations needed by quicksort for sorting a randomly-
ordered input sequence is 1.44n log n, which is asymptotically optimal. There
are several ways to select pivots. For example, the pivot can be the median of
a small number of elements of the sequence, or it can be an element selected at
random. Some pivot selection strategies have advantages over others for certain
input sequences.
8.6.2 Parallelizing Quicksort

Quicksort can be parallelized in a variety of ways. First, consider a naive paral-
lel formulation in the context of recursive decomposition. During each call of
QUICKSORT (see Algorithm 2), the array is partitioned into two parts and each
part is solved recursively. Sorting the smaller arrays represents are two completely
independent subproblems that can be solved in parallel. Therefore, one way to
parallelize quicksort is to execute it initially on a single process; then, when the
algorithm performs its recursive calls, assign one of the subproblems to another
process. Now each of these processes sorts its array by using quicksort and as-
signs one of its subproblems to other processes. The algorithm terminates when
the arrays cannot be further partitioned. Upon termination, each process holds
an element of the array, and the sorted order can be recovered by traversing the
processes as we will describe later. This parallel formulation of quicksort uses
n processes to sort n elements. Its major drawback is that partitioning the array
A[q . . . r] into two smaller arrays, A[q . . . s] and A[s + 1 . . .r], is done by a single
8.6. QUICKSORT 59
process. Since one process must partition the original array A[1 . . . n], the run time
of this formulation is bounded below by Ω(n).
The main limitation of the previous parallel formulation is that it performs the
partitioning step serially. As we will see in subsequent formulations, performing
partitioning in parallel is essential in obtaining an efficient parallel quicksort. To
see why, consider the recurrence equation T (n) = 2T (n/2) + Θ(n), which gives
the complexity of quicksort for optimal pivot selection. The term Θ(n) is due to
the partitioning of the array. Compare this complexity with the overall complexity
of the algorithm, Θ(n log n). From these two complexities, we can think of the
quicksort algorithm as consisting of Θ(log n) steps, each requiring time Θ(n)
that of splitting the array. Therefore, if the partitioning step is performed in time
Θ(1), using Θ(n) processes, it is possible to obtain an overall parallel run time
of Θ(log n). Hence, parallelizing the partitioning step has the potential to yield a
significantly faster parallel formulation.
In the previous paragraph, we hinted that we could partition an array of size
n into two smaller arrays in time Θ(1) by using Θ(n) processes. However, this is
difficult for most parallel computing models. The only known algorithms are for
the abstract PRAM models. Because of communication overhead, the partitioning
step takes longer than Θ(1) on realistic shared-address-space and message-passing
parallel computers. In the following section we present a parallel formulation for a
message-passing platform.
8.6.3 Parallel Formulation for a Message-passing Platform

Let A be an array of n elements that need to be sorted and p = n be the
number of processes. The array A is explicitly distributed among the processes
and each process stores one element d of A. The labels of the processes define
the global order of the sorted sequence. The algorithm starts by selecting a pivot
element v, which is broadcast to all processes. Each process Pi , upon receiving
the pivot, compares its element d with the pivot. In the next phase, the algorithm
first determines which processes will be responsible for recursively sorting the
smaller-than-the-pivot elements (i.e., ≤ v) and which process will be responsible
for recursively sorting the larger-than-the-pivot elements (i.e., > v). Once this is
done, the processes send their elements d to the corresponding processes. After
that, the processes are partitioned into the two groups, and the algorithm proceeds
recursively.
The method to determine the process where the process Pi should send its
element d to, is based on computing the number of processes j = |{Pk | k ≤ i}|
that store elements ≤ v. Thus, is the target process is either j − 1 if d ≤ v, or
p0 + i − j if d > v, where p0 = |{d | d ≤ v}|. Obviously, j and p0 can be computed
using a prefix sum operation and a broadcast.
Algorithm 3 theoQSort(d, i, p)
if p = 1 then
return
end if
r := random element from 0 . . . p − 1{same value in entire partition}
v := d@r{broadcast pivot}
f := Pd ≤ v{1 if d is on left side, 0 otherwise}
j := ik=0 f @k{prefix sum, count elements on left side}
p0 := j@(p − 1){broadcast, result is border index}
if f then
send d to PE j − 1
else
send d to PE p0 + i − j{i − j = ik=0 d@k > v}
P
end if
receive d
if i < p0 then
join left partition; qsort(d, i, p0 )
else
join right partition; qsort(d, i − p0 , p − p0 )
end if
Analysis. The amount of time required to split an array of size n is 2 × Θ(log p)

for broadcasting the pivot element and p0 , Θ(log p) for performing the prefix sum.
Thus, the overall complexity for the split is Θ(log p). This process is repeated for
each of the two subsequences recursively. Assuming that the pivot splits the input
into roughly equal size subsequences, the expected recursion length is Θ(log p).
Thus, the overall complexity of the parallel algorithm is:
T (p) = Θ(log2 p)
8.7 Mergesort
Given k sorted sequences S1 , . . . , Sk and a global rank m, we are asked to find
splitting positions i1 , . . . , ik such that i1 + . . . + ik = m and ∀j, j 0 ∈ 1 . . . k :
Sj [ij − 1] ≤ Sj 0 [ij 0 ]. This is similar to selection, although in this chapter, we are
not interested in the element of global rank m, but rather the splitting positions.
We call finding the splitting positions multiway partitioning. Finding the respective
element is trivial after multiway partitioning, just take the smallest element right
of the split. Merging two sorted sequences S1 and S2 is usually understood as
8.7. MERGESORT 61
efficiently coalescing them into one sorted sequence. Its generalization to multiple
sequences S1 , . . . , Sk is called multiway merging.
It is often used instead of repeated binary merging, since it is more I/O- and
cache-efficient. The high “branching factor“ allows to perform a considerable
amount of work per item. This compensates for the costs of accessing the data on
the (shared or external) memory.
k
{ ······
PE 3
PE 2
PE 1
PE 0
Figure 8.5: Schema of parallel multiway merging, in this case merging a specified
number of elements only.
Pk We parallelize multiway merging using multiway partitioning. Let n =

i=1 |Si |. In parallel, each processing element (PE) i > 1 does multiway par-
titioning on the set of sequences, requesting global rank d i·n p
e. Then, each PE i
merges the part of the input that is determined by the locally computed splitting
positions and by the splitting positions computed by PE i + 1. For i = 0 we use
splitting local positions (0, . . . , 0). The destination offset is also clear from the
ranks. Figure 8.5 shows a schematic view of the algorithm.
The loser tree data structure [28] keeps the next element of each sequence. The
total execution time of the algorithm is O(n/p log k + k log k · log maxi |Si |).
Based on the components explained so far, it is very easy to construct a parallel
sorting algorithm. Each PE sorts about n/p elements sequentially, plus minus 1
due to the rounding. The resulting sorted sequences are merged using parallel
multiway merging afterwards. We assume the sequential algorithm to run in time
O(n log n), i. e. O(n/p log (n/p)) per PE. Substituting k = p and Si = n/p in
the formula for multiway merging results in O(n/p log p + p log p · log (n/p)) The
sequential sorting and the actual merging add up to O(n/p log n), but there is still
the partitioning, so we have in total O(n/p log n + p log p · log (n/p)). Figure 8.6
illustrates the algorithm.
PE 0 PE 1 PE 2 PE 3
sort locally, split exactly
virtually distribute
multiway merge
temporary
memory
copy back
Figure 8.6: Schema of parallel multiway mergesort. For shared memory, the
distribution depicted in the middle row takes place only virtually, i. e. the PEs read
the data from where it is already, there is no copying. The final copying step is
optional, only necessary if the result is expected in the same place as the input.
Chapter 9
External Algorithms
This chapter is mainly based on lecture notes on “Algorithm Engineering“ by Felix

Putze and Peter Sanders.
9.1 Introduction
Massive data sets arise naturally in many domains. Spatial data bases of geographic
information systems like GoogleEarth and NASA’s World Wind store terabytes of
geographically-referenced information that includes the whole Earth. In computer
graphics one has to visualize huge scenes using only a conventional workstation
with limited memory. Billing systems of telecommunication companies evaluate
terabytes of phone call log files. One is interested in analyzing huge network
instances like a web graph or a phone call graph. Search engines like Google and
Yahoo provide fast text search in their data bases indexing billions of web pages.
A precise simulation of the Earth’s climate needs to manipulate with petabytes of
data. These examples are only a sample of numerous applications which have to
process huge amount of data.
For economical reasons, it is not feasible to build all of the computer’s memory
of the fastest type or to extend the fast memory to dimensions that could hold all
relevant data. Instead, modern computer architectures contain a memory hierarchy
of increasing size, decreasing speed and costs from top to bottom: On top, we
have the registers integrated in the CPU, a number of caches, main memory and
finally disks, which are often referenced as external memory as opposed to internal
memory.
The internal memories of computers can keep only a small fraction of the large
data sets. During processing applications need to access the external memory (e. g.
hard disks) very frequently. One such access can be about 106 times slower than a
main memory access. Therefore, disk accesses (I/Os) become the main bottleneck.
63
64 CHAPTER 9. EXTERNAL ALGORITHMS
Figure 9.1: schematic construction of a hard disk
The reason for this high latency is the mechanical nature of the disk access.
Figure 9.1 shows the schematic construction of a hard disk. The time needed
for finding the data position on the disk is called seek time or (seek) latency and
averages to about 3-10 ms for modern disks. The seek time depends on the surface
data density and the rotational speed and can hardly be reduced because of the
mechanical nature of hard disk technology, which still remains the best way to
store massive amounts of data. Note that after finding the required position on the
surface, the data can be transferred at a higher speed which is limited only by the
surface data density and the bandwidth of the interface connecting CPU and hard
disk. This speed is called sustained throughput and achieves up to 80 MByte/s
nowadays. In order to amortize the high seek latency one reads or writes the data
in chunks (blocks). The block size is balanced when the seek latency is a fraction
of the sustained transfer time for the block. Good results show blocks containing a
full track. For older low density disks of the early 90’s the track capacities were
about 16-64 KB. Nowadays, disk tracks have a capacity of several megabytes.
Operating systems implement the so called virtual memory mechanism that
provides an additional working space for an application, mapping an external
memory file (page file) to virtual main memory addresses. This idea supports the
Random Access Machine model in which a program has an infinitely large main
memory with uniform random access cost. Since the memory view is unified in
operating systems supporting virtual memory, the application does not know where
its working space and program code are located: in the main memory or (partially)
swapped out to the page file. For many applications and algorithms with non-linear
access pattern, these remedies are not useful and even counterproductive: the swap
file is accessed very frequently; the data code can be swapped out in favor of data
9.2. THE EXTERNAL MEMORY MODEL 65
blocks; the swap file is highly fragmented and thus many random input/output
operations (I/Os) are needed even for scanning.
9.2 The external memory model

If we bypass the virtual memory mechanism, we cannot apply the RAM model
for analysis anymore since we now have to explicitly handle different levels of
memory hierarchy, while the RAM model uses one large, uniform memory.
Several simple models have been introduced for designing I/O-efficient al-
gorithms and data structures (also called external memory algorithms and data
structures). The most popular and realistic model is the Parallel Disk Model (PDM)
of Vitter and Shriver, see Figure 9.2. In this model, I/Os are handled explicitly by
the application. An I/O operation transfers a block of B consecutive bytes from/to a
disk to amortize the latency. The application tries to transfer D blocks between the
main memory of size M bytes and D independent disks in one I/O step to improve
bandwidth. The input size is N bytes which is (much) larger than M. The main
complexity metrics of an I/O-efficient algorithm in this model are:
• I/O complexity: the number of I/O steps should be minimized (the main
metric),
• CPU work complexity: the number of operations executed by the CPU

should be minimized as well.
The PDM model has become the standard theoretical model for designing and
analyzing I/O-efficient algorithms.
CPU
Memory M
B B
B
Disk 1 Disk i Disk D
Figure 9.2: Vitter’s I/O model with several independent disks

There are some “golden rules” that can guide the process of designing I/O
efficient algorithms: Unstructured memory access is often very expensive as
it comes with 1 I/O per operation whereas we want 1/B I/Os for an efficient
algorithm. Instead, we want to scan the external memory, always loading the
next due block of size B in one step and processing it internally. An optimal
N
scan will only cost scan(N ) := Θ( D·B ) I/Os. If the data is not stored in a
way that allows linear scanning, we can often use sorting to reorder and than
scan it. As we will see in section 9.4, external sorting can be implemented with
N
sort(N ) := Θ( D·B · logM/B N
B
) I/Os.
A simple example of this technique is the following task: We want to reorder
the elements in an array A into an array B using a given “rank” stored in array C.
This should be done in an I/O efficient way.
int[1..N] A,B,C;
for i=1 to N do A[i]:=B[C[i]];
The literal implementation would have worst case costs of Ω(N ) I/Os. For
N = 106 , this would take ≈ T = 10000 seconds ≈ 3 hours. Using the sort-and-
scan technique, we can lower this to sort(N ) and the algorithm would finish in
less than a second:
SCAN C: (C[1]=17,1), (C[2]=5,2), ...

SORT(1st): (C[73]=1,73), (C[12]=2,12), ...
par SCAN : (B[1],73), (B[2],12), ...
SORT(2nd): (B[C[1]],1), (B[C[2]],2), ...
9.3 Stacks and Queues

The contents of this chapter is from the book [40]
Stacks and queues represent dynamic sets of data elements, and support op-
erations for adding and removing elements. They differ in the way elements are
removed. In a stack, a remove operation deletes and returns the set element most
recently inserted (last-in-first-out), whereas in a queue it deletes and returns the set
element that was first inserted (first-in-first-out).
Recall that both stacks and queues for sets of size at most N can be implemented
efficiently in internal memory using an array of length N and a few pointers. Using
this implementation on external memory gives a data structure that, in the worst
case, uses one I/O per insert and delete operation. However, since we can read or
write B elements in one I/O, we could hope to do considerably better. Indeed this
is possible, using the well-known technique of a buffer.
9.4. MULTIWAY MERGE SORT 67
An External Stack
In the case of a stack, the buffer is just an internal memory array of 2B elements
that at any time contains the k set elements most recently inserted, where k ≤ 2B.
Remove operations can now be implemented using no I/Os, except for the case
where the buffer has run empty. In this case a single I/O is used to retrieve the
block of B elements most recently written to external memory.
One way of looking at this is that external memory is used to implement a stack
with blocks as data elements. In other words: The “macroscopic view” in external
memory is the same as the “microscopic view” in internal memory.
Returning to external stacks, the above means that at least B remove operations
are made for each I/O reading a block. Insertions use no I/Os except when the
buffer runs full. In this case a single I/O is used to write the B least recent elements
to a block in external memory. Summing up, both insertions and deletions are done
in 1/B I/O, in the amortized sense. This is the best performance we could hope
for when storing or retrieving a sequence of data items much larger than internal
memory, since no more that B items can be read or written in one I/O. A desired
goal in many external memory data structures is that when reporting a sequence of
elements, only O(1/B) I/O is used per element.
An External Queue
To implement an efficient queue we use two buffers of size B, a read buffer and
a write buffer. Remove operations work on the read buffer until it is empty, in
which case the least recently written external memory block is read into the buffer.
(If there are no elements in external memory, the contents of the write buffer is
transfered to the read buffer.) Insertions are done to the write buffer which when
full is written to external memory. Similar to before, we get at most 1/B I/O per
operation.
9.4 Multiway Merge Sort

Multiway Merge Sort is covered in [16].
We will now study another algorithm based on the concept of Merge Sort which
is especially well suited for external sorting. For external algorithms, an efficient
sorting subroutine is even more important than for main memory algorithms
because one often tries to avoid random disk accesses by ordering the data, allowing
a sequential scan.
Multiway Merge Sort first splits the data into dn/M e runs which fit into main
memory where they are sorted, see Figure 9.3. We merge these runs until only
one is left. Instead of ordinary 2-way-merging, we merge k := M/B runs in a
single pass resulting in a smaller number of merge phases. We only have to keep
one block (containing the currently smallest elements) per run in main memory.
We maintain a priority queue containing the smallest elements of each run in
the current merging step to efficiently keep track of the overall smallest element.
Figure 9.4 exemplifies a 4-way merging.
M
f i=0 i=M i=2M
run sort internal
t i=0 i=M i=2M
Figure 9.3: Run formation
f
make things as simple as possible but no simpler
f __aeghikmnst__aaeilmpsss__bbeilopssu__eilmnoprst
runBuffer st ss ps st
B next internal
M k out ss
t ________aaabbeeeeghiiiiklllmmmnnooppprss
Figure 9.4: Example of 4-way merging with M = 12, B = 2
Every element is read/written twice for forming the runs (in blocks of size B)
and twice for every merging phase. Access granularity is blocks. This leads to the
following (asymptotically optimal) total number of I/Os:
2n 2n l n m
(1 + dlogk #runse) = 1 + logM/B := sort(n) (9.1)
B B M
Let us consider the following realistic parameters: B = 2MB, M = 1GB. For

inputs up to a size of n = 512GB, we get only one merging phase! In general, this
is the case if we can store dn/M e buffers (one for each run) of size B in internal
memory (i.e., n ≤ M 2 /B). Therefore, only one additional level can increase the
I/O volume by 50%.
9.5. INTERNAL WORK OF MULTIWAY MERGESORT 69
9.5 Internal work of Multiway Mergesort

Until now we have only regarded the number of I/Os. In fact, when running
with several disks our sorting algorithm can very well be compute bound, i.e.
prefetching D new blocks requires less time than merging them. We use the
technique of overlapping to minimize wait time for whichever task is bounding our
algorithm in a given environment. Take the following example on run formation (i
denotes a run):
Thread A: Loop { wait-read i; sort i; post-write i};

Thread B: Loop { wait-write i; post-read i+2};
During initialization, runs 1 and 2 are read, i is set to 1. Thread A sorts runs in
memory and writes them to disk. Thread B will wait until run i is finished (and
thread A works on i + 1) and reads the next run i + 2 into the freed space. The
thread doing the more intense work will never wait for the other one.
A similar result can be achieved during the merging step: We have an I/O-
thread (more, if D > 1) and a merging-thread. Merging is I/O bound if and only if
y > DB/2, where y is the number of elements merged during one I/O step1 and
DB/2 is the number of elements that can be moved in and out the merger during
one I/O step.
Consider the I/O bound case. We can show that the I/O thread never blocks:
Consider figure 9.5. The region in which the I/O thread would block is colored
dark grey. Here, there are too few elements in the output buffers (writing them to
disk would be inefficient) and the prefetch and merge buffers are still full. The
blocking region can only be reached from the light grey areas: Merge y elements
into the output buffers and either refill overlap buffers with DB elements2 elements
(left arrow) or write DB items from the output buffers to disk (right arrow). But
these two areas themselves cannot be reached at all by the two possible operations;
therefore, the I/O thread never blocks.
Analogously, the merging thread does not block in the compute bound case.
Here we have two zones of possible blocking: More than 2DB − y elements in
the output buffers or less than kB + y elements in the merge and prefetch buffers 3 .
Again, the constraints on possible transitions in this diagram lead to a single region
from where the blocking areas can be reached which is itself unreachable.
1
Note that y ≤ DB in all interesting cases as otherwise we have an I/O-bound algorithm but
merging is so fast that there is always enough room to fetch D blocks.
2
At the same time y are removed by merging resulting in a total increase of DB − y.
3
y elements are not sufficient in general because these may not be the smallest; but if there are
more than kB elements in merge and prefetch buffers, at least one element can be written to the
output buffer.
r
kB+3DB
#elements in overlap+merge buffers

block I/O
#elements in overlap+merge buffers
kB+2DB
r output
kB+3DB
fetch
kB+2DB+y blocking
kB+DB
kB+2DB output kB+2y
kB+DB+y fetch kB+y
kB block merging
kB+DB
w w
DB−y DB 2DB−y 2DB DB 2DB−y2DB
#elements in output buffers #elements in output buffers
I/O bound case compute bound case
Figure 9.5: The limiting thread never blocks
Overlapping is not always possible. Consider the k runs in figure 9.6. While
merging them, there is no work to do for the prefetching thread. Then, the current
blocks of each run will out of elements simultaneously and merging can only
continue after both buffers are refilled.
1 1B−1 2 3B−1 4 5B−1 6 · · ·

..
.
k 1B−1 2 3B−1 4 5B−1 6 · · ·
Figure 9.6: Overlapping for these k runs is not possible
As internal work influences running time, we need a fast solution for the most
compute intense step during merging: A Tournament Tree (or Loser Tree) is a
specialized data structure for finding the smallest element of all runs, see Figure 9.7.
For k = 2K , it is a complete binary tree with K levels, where each leaf contains
the currently smallest element of one run. Each internal node contains the ’loser’
(i.e., the greater) of the ’competition’ between its two child nodes. Above the root
node, we store the global winner along with a pointer to the corresponding run.
After writing this element to the output buffer, we simply have to move the next
element of its run up until there is a new global winner. Compared to general
purpose data structures like binary heaps, we get exactly log k comparisons (no
hidden constant factors). Similar to the implicit search trees we used for Sample
Sort, Tournament Trees can be implemented as arrays where finding the parent
9.6. EXPERIMENTS 71
node simply maps to an index shift to the right. The inner loop for moving from
leaf to root can be unrolled and contains predictable load instructions and index
computations allowing exploitation of instruction parallelism, see Algorithm 9.8.
1 deleteMin+ 2
insertNext
2 3
4 4 4 4
6 7 9 7 6 7 9 7
4 6 2 7 9 1 4 7 4 6 2 7 9 3 4 7
3
Figure 9.7: A tournament tree
for (int i=(winnerIndex+kReg)>>1; i>0; i>>=1){

currentPos = entry + i;
currentKey = currentPos->key;
if (currentKey < winnerKey) {
currentIndex = currentPos->index;
currentPos->key = winnerKey;
currentPos->index = winnerIndex;
winnerKey = currentKey;
winnerIndex = currentIndex;}}
Figure 9.8: Inner loop of Tournament Tree computation
9.6 Experiments
Experiments on Multiway Merge Sort were performed in 2001 on a 2 × 2GHz
Xeon × 2 threads machine (Intel IV with Netburst) with several 66 MHz PCI-
buses, 4 fast IDE controllers (Promise Ultra100 TX2) and 8 fast IDE disks (IBM
IC35L080AVVA07). This inexpensive (mid 2002) setup gave a high I/O-bandwidth
of 360 MB/s. The keys consisted of 16 GByte random 32 bit integers, run size was
256 MByte, block size B was 2MB (if not otherwise mentioned).
Figure 9.9 shows the running time for different element sizes (for a constant
total data volume of 16 GByte). The smaller the elements, the costlier becomes
internal work, especially during run formation (there are more elements to sort).
With a high I/O throughput and intelligent prefetching algorithms, I/O wait time
400
run formation
350 merging
I/O wait in merge phase
I/O wait in run formation phase
300
time [s] 250
200
150
100
50
0
16 32 64 128 256 512 1024
element size [byte]
Figure 9.9: Multiway Merge Sort with different element sizes
never makes up for more than half of the total running time. This proves the point
that overlapping and tuning of internal work are important.
26
128 GBytes 1x merge
128 GBytes 2x merge
24 16 GBytes
22
sort time [ns/byte]
20
18
16
14
12
128 256 512 1024 2048 4096 8192
block size [KByte]
Figure 9.10: Performance using different block sizes
What is a good block size B? An intuitive approach would link B to the size
of a physical disk block. However, figure 9.10 shows that B is no technology
constant but a tuning parameter: A larger B is better (as it reduces the amortized
costs of O(1/B) I/Os per element), as long as the resulting smaller k still allows
for a single merge phase (see curve for 128GB).
9.7. EXTERNAL PRIORITY QUEUES 73
9.7 External Priority Queues

The material on external priority queues was first published in [46].
We now study a variant of external priority queues which are called sequence
heaps4 : Merging k sorted sequences into one sorted sequence (k-way merging)
is an I/O efficient subroutine used for sorting – we saw this in chapter 9.4. The
basic idea of sequence heaps is to adapt k-way merging to the related but more
dynamical problem of priority queues.
Let us start with the simple case, that at most km insertions take place where
m is the size of a buffer that fits into fast memory. Then the data structure could
consist of k sorted sequences of length up to m. We can use k-way merging for
deleting a batch of the m smallest elements from k sorted sequences.
sorted
1 2 ... k m sequences
B
k−merge insertion
m buffer PQ
PQ
Figure 9.11: A simple external PQ for n < km
A separate binary heap with capacity m allows an arbitrary mix of insertions

and deletions by holding the recently inserted elements. Deletions have to check
whether the smallest element has to come from this insertion buffer. When this
buffer is full, it is sorted, and the resulting sequence becomes one of the sequences
for the k-way merge.
How can we generalize this approach to handle more than km elements? We
cannot increase m beyond M , since the insertion heap would not fit into fast
memory. We cannot arbitrarily increase k, since eventually k-way merging would
start to incur cache faults. Sequence heaps make room by merging all the k
sequences producing a larger sequence of size up to km.
Now the question arises how to handle the larger sequences. Sequence heaps
employ R merge groups G1 , . . . , GR where Gi holds up to k sequences of size
up to mk i−1 . When group Gi overflows, all its sequences are merged, and the
resulting sequence is put into group Gi+1 .
Each group is equipped with a group buffer of size m to allow batched deletion
from the sequences. The smallest elements of these buffers are deleted in batches
of size m0 m. They are stored in the deletion buffer. Fig. 9.12 summarizes the
data structure.
4
By replacing “I/O” with “cache fault”, we can also use this approach one level higher in the
memory hierarchy.
2
1 2 ... k m 1 2 ... k mk 1 2 ... k mk
k-merge k-merge k-merge T3

T1 T2
group group group

buffer 1 m buffer 2 m buffer 3 m
R-merge
deletion buffer m’ m insert heap
Figure 9.12: Overview of the complete data structure for R = 3 merge groups
We now have enough information to understand how deletion and insertion

operations are going to work:
DeleteMin: The smallest elements of the deletion buffer and insertion buffer are
compared, and the smaller one is deleted and returned. If this empties the deletion
buffer, it is refilled from the group buffers using an R-way merge. Before the
refill, group buffers with less than m0 elements are refilled from the sequences in
their group (if the group is nonempty). Figure 9.13 gives an example of a deletion
operation.
DeleteMin works correctly provided the data structure fulfills the heap property,
i.e., elements in the group buffers are not smaller than elements in the deletion
buffer, and in turn, elements in a sorted sequence are not smaller than the elements
in the respective group buffer. Maintaining this invariant is the main difficulty for
implementing insertion.
Insert: New elements are inserted into the insert heap. When its size reaches
m, its elements are sorted (e.g. using merge sort or heap sort). The result is then
merged with the concatenation of the deletion buffer and the group buffer 1. The
smallest resulting elements replace the deletion buffer and group buffer 1. The
remaining elements form a new sequence of length at most m. The new sequence
is finally inserted into a free slot of group G1 . If there is no free slot initially, G1 is
emptied by merging all its sequences into a single sequence of size at most km,
which is then put into G2 . The same strategy is used recursively to free higher level
groups when necessary. When group GR overflows, R is incremented and a new
group is created. When a sequence is moved from one group to the other, the heap
property may be violated. Therefore, when G1 through Gi have been emptied, the
group buffers 1 through i + 1 are merged, and put into G1 . An example insertion
operation is shown in Figures 9.14 and 9.15.
n o p n o p
m q m q
k l r k l r
d j i s w d j i s w
x b g h t u v x b g h t u v
e b f c e b f c
a 3
(a) Deletion of two elements empties insert heap and deletion buffer
n o p n o p
m q m q
k l r k l r
d j i s w d j i s w
x b g h t u v x t u v
e b f c e
b g h
f c
(b) Every Group fills its buffer via k-way-merging, the deletion buffer is filled
from group buffers via M-way-merging
Figure 9.13: Example of a deletion on the sequence
For cached memory, where the speed of internal computation matters, it is also
crucial how to implement the operation of k-way merging. How this can be done
in an efficient way is described in the chapter about Sorting (9.5).
insert(3) v u w v u w
s n t s n t
r m p r m p
q g o q c o
k j f l c h x k j f l h
i e g i
e b d b
d b
b a x a 3
(a) Inserting element 3 leads to overflow of insert heap: it
is merged with deletion buffer and group buffer 1 and then
inserted into group 1
v u w v u w
s n t s n t
r m p k r m p
q c o j q c o
x k j f l h x g l h
e g i e f i
d b d b
b b
a 3 a 3
(b) Overflow in group 1: all old elements are merged and
inserted in next group
v u w n o p
s n t m q
k r m p k l r
j q c o j i s w
x g l h x g h t u v
e f i e f c
d b d b
b b
a 3 a 3
(c) Overflow in group 2: all old elements are merged and inserted
in next group
Figure 9.14: Example of an insertion on the sequence heap (part 1)

n o p n o p
m q m q
k l r k l r
j i s w d j i s w
x g h t u v x b g h t u v
e f c e b f c
d b
b
a 3 a 3
(d) Group buffers are invalid now: merge and inserted them to group 1
Figure 9.15: Example of an insertion on the sequence heap (part 2)
Analysis
We will now give a sketch for the I/O analysis of our priority queues. Let i denote
the number of insertions and an upper bound to the number of deleteMin operations.
First note that Group Gi can overflow at most every m(k i − 1) insertions:
The only complication is the slot in group G1 used for invalid group buffers.
Nevertheless, when groups G1 through Gi contain k sequences each, this can only
happen if
XR
m(k − 1)k j−1 = m(k i − 1)
j=1
insertions have taken place. Therefore, R = logk mI groups suffice.

Now consider the I/Os performed for an element moving on the following
canonical data path: It is first inserted into the insert buffer and then written to
a sequence in group G1 in a batched manner, i.e., 1/B I/Os are charged to the
insertion of this element. Then it is involved in emptying groups until it arrives in
group GR . For each emptying operation, the element is involved into one batched
read and one batched write, i.e., it is charged with 2(R − 1)/B I/Os for tree
emptying operations. Eventually, the element is read into group buffer R yielding a
charge of 1/B I/Os for. All in all, we get a charge of 2R/B I/Os for each insertion.
What remains to be shown (and is ommited here) is that the remaining I/Os
only contribute lower order terms or replace I/Os done on the canonical path. For
example, we save I/Os when an element is extracted before it reaches the last group.
We use the costs charged for this to pay for swapping the group buffers in and out.
Eventually, we have O(sort(I)) I/Os.
In a similar fashion, we can show that I operations inflict I log I key compar-
isons on average. As for sorting, this is a good measure for the internal work,
since in efficient implementations of priority queues for the comparison model, this
number is close to the number of unpredictable branch instructions (whereas loop

control branches are usually well predictable by the hardware or the compiler),
and the number of key comparisons is also proportional to the number of memory
accesses. These two types of operations often have the largest impact on the
execution time, since they are the most severe limit to instruction parallelism in a
super-scalar processor.
Experiments
We now present the results of some experiments conducted to compare our sequence
heap with other priority queue implementations. Random 32 bit integers were used
as keys for another 32 bits of associated information. The operation sequence used
was (Insert − deleteMin − Insert)N (deleteMin − Insert − deleteMin)N .
The choice of this sequence is nontrivial as it can have measurable influence (factor
two and more) on the performance. Figure 9.17 show this: Here we have the se-
quence (Insert (deleteMin Insert)s )N (deleteMin (Insert deleteMin)s )N
for several values of s. For larger s, the performance gets better when N is large
enough. This can be explained with a “locality effect“: New elements tend to be
smaller than most old elements (the smallest of the old elements have long been
removed before). Therefore, many elements never make it into group G1 let alone
the groups for larger sequences. Since most work is performed while emptying
groups, this work is saved. So that these instances should come close to the worst
case. To make clear that sequence heaps are nevertheless still much better than
binary or 4-ary heaps, Figure 9.17 additionally contains their timing for s = 0.
The parameters chosen for the experiments where m0 = 32, m = 256 and
k = 128 on all machines tried. While there were better settings for individual
machines, these global values gave near optimal performance in all cases.
9.8 Semiexternal Kruskal’s Algorithm

A first step for an algorithm that can cope with huge graphs stored on disk is a
semiexternal algorithm: We use Kruskal’s Algorithm but incorporate an external
sorting algorithm. We then just have to scan the edges and maintain the union-find
array in main memory. This only requires one 32 bit word per node to store up to
0..232 − 32 = 4 294 967 264 nodes. There exist asymptotically better algorithms
but these come with discouraging constant factors and significantly larger data
structures.
9.8. SEMIEXTERNAL KRUSKAL’S ALGORITHM 79
bottom up binary heap

200 bottom up aligned 4-ary heap
(T(deleteMin) + T(insert))/log N [ns]
sequence heap
150
100
50
0
1024 4096 16384 65536 218 220 222 223
N
Figure 9.16: Runtime comparison for several PQ implementations (on a 180MHz

MIPS R10000)
s=0, binary heap

s=0, 4-ary heap
(T(deleteMin) + T(insert))/log N [ns]
s=0
s=1
128 s=4
s=16
64
32
256 1024 4096 16384 65536 218 220 222 223
N
Figure 9.17: Runtime comparison for different operation sequences

Chapter 10
Online Algorithms
Online algorithms are a class of algorithms which process their input in a serial
fashion, where future information is only revealed in parts. The algorithms must
fully process each part and produce a partial output before receiving further input.
No information is known about the further input sequence. Therefore, online
algorithms are restricted to making decisions based solely on past events. The
challenge is to design algorithms that perform well under this lack of knowledge
about the future.
In contrast, many algorithms previously developed in these course notes require
the whole input to be available at the beginning and do extensive computation on
the entirety of the problem data. These are called offline algorithms.
However, many practical problems have an intrinsic online nature, for example:
• Paging in a virtual memory system.
• Renting skis vs. buying skis.
• Scheduling jobs in a multi-user system.
• Routing in communication networks.
• Google placing advertisements on websites.
Furthermore, some of the well-known algorithms are online algorithms: insertion-

sort can sort an input sequence that is revealed only one item at a time.
Note that online algorithms do not make particular probabilistic assumptions
on the unknown input. With such assumptions many of the examples above are
typical applications of classic queueing theory.
Most of this chapter is based on the text book [10].
81
82 CHAPTER 10. ONLINE ALGORITHMS
10.1 Competitive Analysis

Prior to focusing on two of these examples we will develop a method to compare
different online algorithms. Their quality is determined by competitive analysis
against an optimal offline algorithm. Thus this analysis quantifies the error an
online algorithm makes due to lack of information.
Definition 7
Given an optimization problem with input I, for which an offline algorithm OPT
can compute an optimal solution using a cost measure OPT(I), then an online
algorithm ALG is called c-competitive if there is a constant α such that for all
finite input sequences I
ALG(I) ≤ c · OPT(I) + α
where ALG(I) is the cost of computing input I using ALG.

If moreover the additive constant α ≤ 0, then one says ALG is strictly c-
competitive
ALG is strictly c-competitive for
ALG(I)
c = sup .
I OPT(I)
The ratio c is also called the competitive ratio of ALG.

Note that c is a worst-case performance measure, depending on the input
sequence with worst relative performance.
An algorithm is simply called competitive if it is c-competitive for some real
constant c. Note that if ALG is c-competitive, then it is also c0 -competitive for all
real constants c0 > c. The goal is to find online algorithms with minimal c.
10.1.1 Example: Ski Rental

We now exemplify this analysis method on a typical online problem: ski rental.
Example Say you are planning to spend your winter holiday at a cost-efficient
hotel in the Alps near Ischgl and do not own a pair of skis. You have the choice of
buying skis for 300 euros from your favorite sports equipment store at home, or to
rent skis for a total of 50 euros in Ischgl. Considering that you do not know how
often you will go skiing, should you rent or buy skis?
As a first simple online algorithm we suggest buying skis on the sixth trip.
Using competitive analysis we can calculate the competitive ratio for this algorithm:
10.2. PAGING – DETERMINISTIC ALGORITHMS 83
• If you make five or fewer skiing trips, you will pay the optimal cost of 50
euro per trip and thus the ratio to an optimal offline algorithm is 5·50
5·50
= 1 for
these inputs.
• If you make six or more skiing trips, you pay 250 + 300 = 550 euros in total.
However, an optimal offline algorithm can take advantage of its information
about the future and pays only 300 euro. For these inputs the ratio is 550
300
= 11
6
,
11
thus our simple algorithm is 6 -competitive.
We now show that no better online algorithm can be devised for this problem.
• Suppose we buy skis earlier, say on the trip x < 6. For fewer than x trips,
the ratio stays 1, but for exactly x trips an optimal algorithm pays 50 x euros
while the online algorithm pays 50(x − 1) + 300. Thus it has a competitive
ratio of
50x + 250 5
=1+ ≥2
50x x
• Suppose we buy skis later, say on the trip x > 6. For exactly x trips an
optimal offline algorithm pays only 300 euros while the online algorithm
pays again 50(x − 1) + 300 euro. The competitive ratio in this case is
50x + 250 5+x
= ≥2
300 6
So both buying times x < 6 and x > 6 yield an algorithm with a competitive ratio
greater than 2. All other algorithms, like renting skis while already possessing a
pair or buying two pairs, are obviously inferior and thus buying after 6 trips is the
optimal algorithm in absence of complete knowledge.
The ski rental example is prototypical of a class of problems in which you
must decide whether to pay smaller repeating costs or a large one-time cost, which
reduces or eliminates the repeating costs. Other examples are whether to buy
a BahnCard or monthly tickets for a public transport system. The idea behind
optimal online algorithms for these problems is to not pay the large one-time cost
until you would have paid the same amount in smaller repeating costs.
10.2 Paging – Deterministic Algorithms

The paging problem is an example of a well-studied set of online algorithms for
which optimal deterministic and randomized algorithms are known. In this section
the paging problem is formalized and five deterministic online algorithms for this
problem are compared via competitive analysis. Randomized online algorithms for
paging are discussed in the next chapter.
page access sequence σ

k pages cache
large memory
Figure 10.1: Paging in a virtual memory system
The Paging Problem Consider a two-level virtual memory system as in fig-

ure 10.1: the system consists of a large main memory storage to which access
is relatively slow. To speed up processing a faster cache is used to temporar-
ily store accessed memory locations. Transfer between both memory stores is
possible only using fixed-sized memory blocks called pages. The slow main mem-
ory P = {p1 , . . . , pN } contains N pages of which any k-subset can be stored
in the cache where k < N . A program specifies a sequence of page accesses
σ : {1, . . . , n} → {1, . . . , N } of finite length |σ| := n. Given a request for a page
pσ(i) , the virtual memory system must make pσ(i) available in the fast memory. If
pσ(i) is already loaded in the fast memory we call this a cache hit, otherwise a page
fault or cache miss occurs and the page must be transferred from main memory to
cache memory. Due to the limited size of fast memory, the paging algorithm must
decide which of the cached pages is to be replaced (evicted from cache). Goal of
the paging algorithm must be to minimize the number of cache misses.
The memory systems of today’s computer architectures are composed of a

multi-level virtual memory system containing CPU registers, three cache levels on
the CPU (L1, L2 and L3), physical RAM and swap partitions on hard-drives. For
our discussion on online algorithms we will use the high-level abstract model and
disregard technical details like k-associative caches or page sizes. We assume a
cache size of k pages, allow cache hits to be free (very fast cache) and penalize a
cache miss with cost 1.
Paging is a natural online problem as future memory access patterns are gen-
erally not known in advance. Table 10.1 contains the most basic and well-known
deterministic paging algorithms.
algorithm strategy: replace page ...

FIFO First In / First Out which was longest in cache
LIFO Last In / First Out which was newest in cache
LFU Least Frequently Used requested least number of times
LRU Least Recently Used with earliest most recent request
FWF Flush When Full evict all pages when cache full
LFD Longest Forward Distance requested farthest in the future
Table 10.1: Deterministic paging algorithms
All algorithms given above except LFD are online algorithms. LFD is an
offline algorithm because it requires knowledge of the future page access order.
Furthermore, all these algorithms except FWF are called demand paging: they
never evict a page unless there is a page fault. One can show that any page
replacement algorithm can be modified to be demand paging without adding a cost
penalty. For FWF this leads to an algorithm which evicts an arbitrary page when
the cache is full.
Figure 10.2 shows a snapshot of LRU and LFD working on an example page
request sequence. The current cache, with k = 3 in both examples, is depicted
in the second row. Conceptually, LRU keeps a backwards reference to the last
time each cached page was used. In the current state shown, the red page will be
evicted to make room for the requested blue page, as the red page has been longest
in cache, relative to the most recent access. LFD, on the over hand, is an offline
algorithm and can therefore make use of the future request sequence. Each page in
the cache references its next access position. When a page must be evicted, LFD
can pick the page which has its next request farthest in the future. In the example,
access to the green page causes a page to be evicted and because the yellow one
will be requested latest of all in the cache, LFD will remove the yellow one.
A further interesting topic is how efficiently the replacement strategies can
be implemented. Simplicity is a key factor for paging algorithms because long
calculations would impact the fast memory’s access performance on cache misses.
However, for competitive analysis of the algorithms with page faults as cost metric,
the implementation complexity is immaterial.
σ σ
LRU LFD
Figure 10.2: Example of LRU and LFD paging algorithm

10.2.1 Longest Forward Distance

Competitive analysis needs an optimal offline algorithm OPT. In the case of
paging algorithms this requirement can not only be met, but the optimal page
eviction sequence can also be calculated easily.
Theorem 21
LFD is an optimal offline algorithm for paging.
Proof.
We prove this theorem by showing that any optimal offline algorithm OPT can be
modified to act like LFD on any given page request sequence σ without increasing
the number of page faults.
This is shown via inductive application of the following claim: given any
paging algorithm ALG and a page request sequence σ, we can construct an offline
algorithm ALG0i for any i = 1, . . . , |σ| so that following three properties are
satisfied:
(i) ALG0i processes the first i − 1 paging requests exactly as ALG does.
(ii) If the i-th request causes a page fault, then ALG0i evicts the page with
longest forward distance from cache memory. So for the i-th request the new
algorithm behaves the same as LFD.
(iii) The cost in page faults of ALG0i does not increase: ALG0i (σ) ≤ ALG(σ).
If this claim holds, the theorem can be proven by applying it n = |σ| times to
any optimal offline algorithm OPT. Start by applying the claim to OPT with
i = 1 to obtain OPT01 , which acts like LFD on the first request. Then apply the
claim to OPT01 with i = 2 to obtain OPT002 and so on. After n steps the algorithm
OPT(n)n acts identically to LFD on the request sequence σ. The number of page
fault of OPT(n)n has not increased due to (iii), nor can it have decreased as OPT
was already optimal.
To finish the proof, it remains to show the claim. Given a paging algorithm
ALG, we first construct ALG0i . Figure 10.3 shows an illustration of this construc-
tion.
For the first i − 1 requests, both algorithms must behave identically. Therefore,
after the (i − 1)-th request the cached page sets are identical. Only after handling
the i-th request for page x these sets might differ, as ALG0i will pick the page a
with longest forward distance to evict. We denote the cached page set of ALG after
i requests by X ∪ {a} and of ALG0i by X ∪ {b}, where X is the common page
set of size k − 1. If a = b, then there was no page fault or ALG already behaves
like ALG0i for the i-th page request and all further page requests can be handled
1st i-th
ALG ... ... .

σ ...
ALG0i ... ... .
LFD
x= a= b= r=
Figure 10.3: Example showing cached page sets of ALG and ALG0i
identically. We therefore assume a 6= b. For the requests following the i-th, ALG0i
shall mimic ALG, except when either page a or b is requested. This strategy is
possible, because both algorithms’ page sets share k − 1 pages. If at some point
ALG evicts a, then ALG0i can evict b instead, and from there on both are identical.
Otherwise b will be requested before a, as a has longest forward distance. When b
is requested, ALG must evict a page r to make room for b, while ALG0i still has
b in cache. If r = a, then both page sets are equal from this point on. Otherwise
ALG now has X̃ ∪ {a} cached while ALG0i has X̃ ∪ {r} with a common page set
X̃. Again ALG0i can continue like ALG, until eventually a is requested. At this
time ALG0i can evict r to make room for a, after which both page sets are identical.
The constructed algorithm ALG0i satisfies conditions (i) and (ii). The difficulty
lies in showing (iii): ALG0i may not suffer more page faults than ALG. An
additional page fault can only occur before both page sets are identical again. Both
algorithms behave the same except for requests to pages a and b, so both incur the
same page faults for access to other pages. If a 6= b, then ALG suffers a page fault
when b is accessed before a. ALG then evicts r, while ALG0i still holds r. The
page r then takes the role of b. This can happen multiple times further increasing
ALG’s cost. Eventually a is requested, for which ALG0i suffers a page fault and
can then remove r. This page fault is covered by the previous one when r was
evicted by ALG. Thus ALG0i never suffers more page faults than ALG.
10.2.2 Competitive Paging Algorithms

We now have an optimal paging algorithm OPT = LFD. Of course, LFD is an
offline algorithm and we now try to find the best online algorithm by competitive
analysis against LFD.
Of the previously represented online algorithms, two can immediately be ex-
cluded: LIFO and LFU are not competitive, as we can show that their competitive
ratio c can grow infinitely large.
Let k be the cache size and P = {p1 , . . . , pk , pk+1 } the pages in the slow main
memory. When processing the request sequence
σ = p1 , p2 , . . . , pk , (pk+1 , pk )m
for some integer m, LIFO will first cache pages p1 , . . . , pk . For the (k + 1)-st
request LIFO will evict pk to make room for pk+1 , which it will remove again
on the (k + 2)-nd request and so on. Thus LIFO will incur page faults on every
request after the k-th one, in total 2m, whereas LFD only suffers k + 1 page faults
in total. We can therefore provide for every constant c a request sequence of length
|σ| > k + c(k + 1), which violates definition 7.
Again let k be the cache size and P = {p1 , . . . , pk , pk+1 } the pages in the slow
main memory. For every integer m we can consider the request sequence
σ = p1m , p2m , . . . , pk−1

m
, (pk , pk+1 )m
LFU will load the pages p1 , . . . , pk−1 for the first k − 1 page faults. After the
first (k − 1)m requests however, LFU will suffer a page fault for each of the 2m
following requests. For the first page fault pk , the cache still has one empty slot, but
following this each request forces LFU to evict a page: pk when pk+1 is requested
and vice versa. Since m can be made arbitrarily large and LFD requires only k + 1
page faults, the competitive ratio of LFU is unbounded.
Following this sobering result, we show another theorem, which establishes a
lower bound on the competitive ratio for any online paging algorithm.
Theorem 22
No deterministic online paging algorithm can achieve a better competitive ratio
than k, where k is the cache size.
Proof.
Let ALG be an online paging algorithm working on a cache size of k. We show
there is an arbitrarily long request sequence σ for which
|σ| = ALG(σ) ≥ k · LFD(σ)
Assume there are k + 1 pages in slow memory and that the cache already holds k
pages. Thus there is always a page not in the cache. Since ALG is deterministic,
a malicious page request sequence can be constructed for which each request
results in a page fault: in each step a page not in cache is requested. For such a
sequence ALG(σ) = |σ|, whereas LFD incurs at most one page fault on every k
requests: i.e. LFD(σ) ≤ |σ|k
. This is due to LFD being able to have the next k page
requests in cache for each position in the sequence. Thus ALG(σ) ≥ k · LFD(σ)
or ALG(σ)
LFD(σ)
≥ k.
10.2.3 Conservative Paging Algorithms

For the continued analysis of online paging algorithms in this chapter, we now
refine the paging problem to allow different sizes of cache memory allotted to the
online and optimal offline algorithm. In the (h, k)-paging problem we compare
an optimal offline algorithm working on h cache pages to an online algorithm with
cache k ≥ h. With h = k we still include the equal cache size scenario, while
for h < k we give the online algorithm the advantage of more cache memory as
compensation for its lack of knowledge. This approach may lead to more realistic
competitive ratios.
Note that some paging algorithms incur what is known as Bélády’s anomaly:
on some request sequences these algorithms suffer fewer page faults with a smaller
cache size. FIFO for example is prone to this anomaly. It shows that the effects of
different cache sizes can be subtle.
To conduct competitive analysis of a broader range of paging algorithms we
define the following class:
Definition 8
A paging algorithm is called conservative if it incurs at most k page faults on any
request sequence containing at most k distinct pages.
LRU, FIFO and even FWF are conservative. LFU and LIFO are not con-
servative, they are not even competitive. The page requests sequences in the
proof of non-competitiveness are good starting points for the reader to show this
observation.
We now show that all conservative paging algorithms have at least the following
competitive ratio.
Theorem 23
Any conservative online algorithm ALG, which works on a cache of size k, is
k
k−h+1
-competitive to any optimal offline algorithm OPT with h ≤ k pages of
cache.
Proof.
Let σ be a fixed page request sequence. We first partition this sequence into phases:
phase 0 is the empty sequence and for i > 0 let phase i be the longest sequence of
page requests following phase i − 1 that contains at most k distinct pages. Such a
σ
1 2 3 4 5 6
Figure 10.4: The 4-phase partition of an example sequence

phase i
f g
k distinct pages
OPT has h pages in cache including f
Figure 10.5: Illustration of the page faults of OPT
partition is well-defined and called a k-phase partition of σ. It does not depend on

any particular algorithm. (see figure 10.4 for an example)
The k-phase partition is constructed so that the definition of a conservative
paging algorithm can easily be applied: in each sub-sequence a conservative
algorithm will have at most k page faults.
Now consider how many page faults an optimal offline paging algorithm OPT
working with h ≤ k cache pages can have in a k-phase (Figure 10.5 shows an
illustration). Take any phase i ≥ 1, then name the first page f and regard the
input sequence following f up to the first request of phase i + 1, which is labeled
g if it exists. An optimal paging algorithm using h cache pages will have h − 1
pages other than f in fast memory. However, since the input sequence following
f contains k distinct pages, at least k − (h − 1) pages are missing from cache
memory. Therefore at least k − h + 1 unavoidable page faults will occur even to
OPT in the sequence following f up to and including g. If the last phase contains
only k 0 < k distinct page requests, the argument shows that at least k 0 − h + 1 page
faults occur for OPT. For competitive analysis this can be ignored, as it is only an
additive constant cost.
So during each k-phase the competitive online algorithm ALG suffers at most
k page faults, while for OPT at least k − h + 1 page faults occur. Therefore, we
showed the following inequality for every input sequence σ
k
ALG(σ) ≤ · OPT(σ) + r
k−h+1
where r ≤ k is the maximum number of page faults of ALG for the last phase.
We can apply this theorem to some example problems: for (k, k)-paging every
conservative online algorithm, like FIFO, LRU or FWF, has a competitive ratio
of k, which we showed in theorem 22 to be a lower bound for this problem. These
online algorithms are therefore optimal in this context.
Applying theorem 23 with h = k2 shows that conservative online algorithms are
2-competitive for the (h, 2h)-paging problem. The lower bound of k only applies
to the special case (k, k)-paging.
10.3. RANDOMIZED ONLINE ALGORITHMS 91
10.3 Randomized Online Algorithms

In this chapter we will explore how to use randomization to avoid the lower bound
of k for (k, k)-paging.
Randomized algorithms are allowed to make random choices in their online
decision process. For analysis of these algorithms we view the problem as a
game between an online player, the online algorithm, and a malicious adversary.
This adversary attempts to pick a page request sequence which maximizes the
competitive ratio. In the deterministic case there was just one natural adversary
model, which knows both the algorithm and its current cache, and from these it
could easily determine which request would lead to a page fault (theorem 22).
For analysis of randomized algorithms we must formalize which information is
available to an adversary. All adversary models assume that the algorithm itself,
the current page set and the probability distributions from which the algorithm
draws random choices are known. The following models distinguish between how
much of the drawn random variates is known in advance and when the page request
sequence must be fixed.
• An oblivious adversary does not know the random choices and must choose
the entire request sequence in advance. This is the weakest adversary as it
can only exploit a-priori knowledge about the algorithm.
• An adaptive-online adversary is allowed to determine the next page request
by inspecting all random choices previously made by the online algorithm.
This models an input-output scenario where the adversary receives the output
and can then devise new malicious inputs. The input requests are determined
online, without knowing the future random choices.
• An adaptive-offline adversary knows the outcome of all random choices the
online algorithm is going to make in advance and can determine the page
request sequence offline beforehand.
These three adversary models are also sometimes called weak, medium and
strong models. Randomization does not help against adaptive-offline adver-
saries: one can show that if there exists a c-competitive randomized online al-
gorithm against an adaptive-offline adversary, then there also exists a deterministic
c-competitive online algorithm for the same problem.
We will focus on analysis against an oblivious adversary. Since this adversary
must define the page request sequence σ in advance, we can assign it the same cost
OPT(σ) as an optimal deterministic paging algorithm requires on the sequence.
To rank a randomized online algorithm ALG against an oblivious adversary, we
need to regard ALG(σ) as a random variable of the possible random choices that
the algorithm makes. We therefore provide
Definition 9
A randomized online algorithm ALG is called c-competitive against an oblivious
adversary if there is a constant α such that for all finite input sequences σ
E(ALG(σ)) ≤ c · OPT(σ) + α
where E(ALG(σ)) is the expected value of the random variable ALG(σ) over all
possible random choices in ALG.
Note that the randomized competitive ratio is not defined as a worst-case

measure over all random choices in ALG, but is defined as the worst case over
all input sequences. For the randomized competitive ratio one needs to finds the
input sequence with the worst expected cost of ALG over all possible random
choices. For the oblivious adversary, the input sequence does not depend on the
random choices of ALG and therefore OPT(σ) can be fixed independently of
ALG’s choices.
10.3.1 Randomized Marking Algorithms

We now distinguish a special class of online paging algorithms called marking
algorithms, which includes a specific randomized paging algorithm that we will
analyze against an oblivious adversary.
Definition 10
We call a paging algorithm a marking algorithm, if its operations can formally
be described using a one bit mark for each page of the slow memory. Initially, with
an empty cache, all pages are unmarked. On request of a page, this page is marked
and loaded into cache memory, if it is not already there. Only unmarked pages may
be evicted if the algorithm needs to removed a page from cache. If all pages are
marked and there is a page fault, then all marked pages are unmarked and one of
them is evicted.
The marking/unmarking phases of this definition are similar but not identical
to the k-phase partition used in the proof of theorem 23. Marking algorithms
can only differ in the decision, which of the unmarked pages is evicted. FWF
is maybe the most primitive marking algorithm, as it just evicts all unmarked
pages on a page fault requiring page eviction. More elaborately, LRU is also a
marking algorithm and it picks the unmarked page least recently used. FIFO is
not a marking algorithm.
For competitive analysis we target RMARK, which is the marking algorithm
that evicts an unmarked page chosen uniformly at random when required.
phase i − 1 phase i (k distinct pages)

k distinct pages new old
mi pages j-th old page
Figure 10.6: Worst-case for RMARK: new and old pages in a k-phase partition
In context of the following analysis we need a common definition: for any

k ∈ N1 the k-th harmonic number Hk is defined as
k
X 1 1 1 1
Hk = =1+ + + ··· +
i=1
i 2 3 k
The harmonic number can be bounded by ln k < Hk ≤ 1 + ln k for all k ≥ 1.
Theorem 24
The randomized paging algorithm RMARK is (2Hk )-competitive against an obliv-
ious adversary.
Proof.
We again look at the k-phase partition of a fixed request sequence σ and analyze
the expected number of page faults in phase i for RMARK. At the beginning of
phase i there are k pages in cache as requested in phase i − 1. We call these pages
old, whereas pages that are requested in phase i, but not in phase i − 1, are labeled
as new. Every page in phase i is either old or new, let mi be the number of new
pages. Access to each of the mi new pages must result in a page fault, whereas
old pages may still be in the cache from the previous phase. Consecutive accesses
to the same old or new pages do not incur any further page faults and we can
therefore concentrate on distinct page requests. The worst possible order of page
requests for the online algorithm is first all new pages followed by the old pages
(see figure 10.6), because each new page request evicts an old page which might
be used later. The online algorithm, however, cannot take into account which of
the old pages will be requested.
So there are mi unavoidable page faults for the new pages, which are marked
and therefore cannot be evicted in this phase. We need to determine the expected
value of page faults for the k − mi old pages. For this consider the probability that
the j-th old page is still in the cache when it is first requested: since it is unmarked
it must be one of the k − mi − (j − 1) unmarked old pages in the cache (there
are only k − mi slots for old pages and j − 1 are already taken). Dividing this
by the total number of unmarked old pages k − (j − 1) in this request sequence
(including those not requested) yields the probability of a cache hit for the j-th old
page. So, the probability for the j-th old page to cause a page fault is
k − mi − (j − 1) mi
1− =
k − (j − 1) k − (j − 1)
To get the expected number of faults in phase i, we sum over all old pages:
k−m
Xi
mi 1 1 1
mi + = mi + mi + + ··· +
j=1
k − j + 1 k k − 1 mi + 1
= mi + mi (Hk − Hmi ) = mi (Hk − Hmi + 1) ≤ mi Hk
If we say the request sequence σ has N phases, then we get
N
X
E(RMARK(σ)) ≤ Hk mi
i=1
We now have an upper bound for the cost of RMARK, to complete the competitive
analysis we require a lower bound on OPT(σ). Here we need to look at two
adjacent phases i − 1 and i. Since there are mi new pages in phase i, there must be
at least k + mi distinct pages in the request sequence consisting of phases i − 1 and
i. Thus in phase i − 1 and i even OPT must incur at least mi page faults. Hence,
we have the following lower bound.
N
1X
OPT(σ) ≥ mi
2 i=1
By dividing the two bounds we can calculate a randomized competitive ratio of

RMARK:
Hk N
P
E(RMARK(σ)) mi
≤ 1 PNi=1 = 2Hk
OPT(σ) 2 i=1 m i

So with RMARK we have a simple randomized algorithm with a competitive
ratio of 2Hk . Note that Hk is usually a lot smaller than k (for example H8 ≈ 2.72),
therefore even this crude randomized algorithm can be better in practice than LRU
or FIFO, and probably is even easier to implement.
10.3.2 Lower Bound for Randomized Paging Algorithms

We now consider whether it is possible to design a randomized paging algorithm,
that has a better competitive ratio than RMARK. This is indeed possible but more
complicated. Instead we show a lower bound for the competitive ratio for any
randomized paging algorithm.
σ U M U M U M M U M U M M M U U M M U
M Marked page
1
U Unmarked page E(ALG) ≥ k−j+1
Figure 10.7: Eight sub-phases of a (super-)phase each containing some marked

pages and one new unmarked page
Theorem 25
No randomized online paging algorithm working with cache size k and slow
memory size N ≥ k + 1 can have a better competitive ratio than Hk against an
oblivious adversary.
Proof.
Let ALG be a randomized online algorithm with k cache pages and N = k + 1
pages in slow memory. Since an oblivious adversary must construct the request
sequence in advance, we must devise a method which yields a sequence that
results in a competitive ratio of Hk for any ALG, then the competitive ratio, as
a worst-case measure over all sequences, is greater for all possible algorithms.
The adversary knows about the algorithm’s randomization mechanism and can
therefore maintain a probability vector (p1 , . . . , pN ), where pi is the probability for
page i of the slow memory to not be in the current cache set. This vector can be
calculated by the adversary offline using the probabilityPdistributions in ALG for
each position in the request sequence iteratively, while N i=1 pi = 1 holds in each
step. Thus the expected cost for a request to page i is pi .
We continue to construct a malicious sequence based on the k-phase partition
schema used in the proof to theorem 23. In each phase ALG will suffer at least
Hk page faults, while OPT can have k − 1 pages in cache and incurs only 1 page
fault. This then results in a competitive ratio of at least H1k . To construct the phases,
the adversary uses a “reverse” marking algorithm and thus must store the set of
currently marked pages.
Each phase starts with all pages unmarked. During the phase all k pages will
be accessed at least once and thereby marked. The trick is to use the probability
vector to increase the expected number of page faults between two accesses to
unmarked pages. This is possible because there are k + 1 pages in total.
Each k-phase is composed of k sub-phases, which in turn are composed of zero
or more requests to marked pages already in cache and terminate with one request
to an unmarked page, which may or may not be in cache (see figure 10.7). At the
start of the j-th sub-phase there are exactly k − (j − 1) unmarked pages. As each
sub-phase increases the number of marked pages by one, thus at the end of the
k-th sub-phase all will become unmarked again and a new (super-)phase begins
with one marked and k unmarked pages (N = k + 1). If we can ensure that the
1
expected cost of ALG for sub-phase j is k−j+1 , then the total expected cost for
phase i would be
k
X 1 1 1 1
= + + · · · + = Hk (10.1)
j=1
k−j+1 k k−1 1
At the start of the j-th sub-phase the set of marked pages M has cardinalityP
j, while
the set of unmarked pages U contains (k + 1) − j items. Consider γ := i∈M pi ,
which is the total probability of any marked page not being in the cache. If γ = 0,
1
there is an unmarked page a ∈ U with pa ≥ k+1−j as |U | = k+1−j. By requesting
this page the expected cost can be met and the sub-phase ended. Otherwise with
γ > 0, there exists a marked page m with pm > 0. Let ε := pm and request this
page as the first in the sub-phase. Then generate subsequent requests to marked
pages using following loop:
1
While the expected cost for ALG is less than k+1−j and γ > ε, request
the marked page x with maximal probability px .
The expected cost for ALG increases with each iteration of the loop as px > 0,
so the loop must terminate. (Actually, each additional request contributes at least
γ
|M |
> |Mε | > 0 cost, because γ > ε.) If the loop terminates with expected cost
1
≥ k+1−j , then the adversary can pick any unmarked page to end the sub-phase. If,
1
however, the loop terminates with expected cost < k+1−j , then γ ≤ ε must have
been the reason. In this case the adversary must request the unmarked page b with
maximal pb , for which pb ≥ 1−γ|U |
1−γ
= k+1−j . The total expected cost for ALG is in
this case ε + pb , which is
1−γ 1−ε 1
ε + pb ≥ ε + ≥ε+ ≥
k+1−j k+1−j k+1−j
1
Thus the total expected cost for any sub-phase is ≥ k+1−j and therefore for the
whole phase at least Hk (due to eq. 10.1). The variable ε is used as an arbitrary
bound to the number of marked pages requested.
In each phase the expected cost for ALG is at least Hk , while the cost for OPT
is at most 1. Thus the competitive ratio of ALG is at least Hk . Since we did not
make any assumption about how ALG works on the page request sequence, this
lower bound holds for all randomized online paging algorithms.
This closes these chapters on competitive analysis of online algorithms. We
have shown many analytic results that tend to be very pessimistic about competitive
performance. But keep in mind that this performance is always against a much
more powerful, optimal offline algorithm. The results of these optimal algorithms
are often also not easy to come by. In practice the best paging algorithms are LRU
and RMARK (see figure 8.1 in [10]).
In spite of these disadvantages, competitive analysis is a solid foundation to
give guarantees on the worst-case performance of online algorithms.
Chapter 11
Geometry
Geometrical algorithms cover a wide range of algorithms defined for sets of

geometrical objects. These objects may include points, line segments, planes,
curves and many other. In this chapter we will discuss a small selection of these
elements, present fitting data structures and some algorithms that perform on these
elements.
11.1 Application Domains

Geometric algorithms influence a lot of fields. The most prominent include com-
puter graphics, robotics, geographic information systems (GIS) or computer aided
design (CAD). Geometric algorithms enable us, to handle such amounts of data as
is necessary for modern GIS or in the latest PC-game. Whether we have to perform
hidden object removal in complex 3D-scenes or efficiently retrieve data about gas
pipes within a given region, efficient data structures and algorithms are necessary.
Typical questions for geometrical algorithms may include intersections between
line segments, to calculate the convex hull of a set of points or localise a point in
an partition of the d-dimensional space.
In this chapter we will study three distinct problems in detail. Those are line
segment intersection, convex hull and the smallest enclosing circle problem.
11.2 Basic Elements and Definitions

Our algorithms will focus on two very basic elements. The elements we use are
the points and line segments. For completeness we give following definitions:
Definition 11
A point p is defined as a vector ∈ Rd . We use p.i to refer to the i-th element of the
99
100 CHAPTER 11. GEOMETRY
vector. In two- as well as three-dimensional space, we also use p.x, p.y and p.z to
refer to the elements coordinates.
Definition 12
Given two points a and b, we define a line segment ab as {α · a + (1 − α) · b|α ∈ [0, 1]}.
Definition 13
A polygon in an d-dimensional space is a set of line segments given as a list of
points P = p1 , . . . , pn with pn = p1 and pi , pi+1 for i ∈ [1, n − 1] defining the
outline of the polygon.
We call a polygon convex if for any two points a and b in P the line segment
ab is contained in P .
11.3 Data Structures

We will not study data structures in this place. We will only give a small example for
geometric data structures as presented in Figure 11.1.
The Figure shows the combination
of the Voronoi diagram (dashed
lines) and the Delaunay triangu-
lation for a given set of points.
Voronoi diagrams define the re-
gions that are closest to one of
the given points and therefore may
be interesting whenever the trad-
ing area of a supermarket chain or
something similar has to be com-
puted. Triangulations become inter-
esting whenever a plane has to be
represented, as graphic cards need
triangles to draw surfaces. The
Delaunay triangulation of a set of Figure 11.1: Voronoi Diagram and Delaunay
points is a triangulation directly Triangulation
connected to the Voronoi diagram.
The triangulation introduces one edge for any two vertices that share a common
boundary in the Voronoi diagram. A Delaunay triangulation yields a “natural
looking” triangulation which can be utilized for terrain rendering purposes for
example.
11.4. LINE SEGMENT INTERSECTION - CASE STUDY 101
11.4 Line Segment Intersection - Case Study

The line segment intersection problem can be defined as follows: Given S a set
S = {s1 , . . . , sn } of line segments, calculate the set of all intersections s,t∈S s∩t.
First we will discuss the orthogonal line segment intersection problem. In this
problem the line segments are layed out on an orthogonal grid.
11.4.1 Orthogonal Line Segment Intersection

This problem has its application domains in circuit development (where do con-
ductors intersect) or in GIS (road crossings, bridges).
First we have a look at a naive implementation of the line segment intersection
algorithm.
foreach {s, t} ⊆ S do
if s ∩ t 6= ∅ then
output {s, t}
The algorithm obviously performs in Θ(n2 ) which is far to slow for large data
sets. Naturally the questions arises: what is the best we can do?
Since we have to look at each edge
at least once and have to output
every intersection, we get a mini-
...
mal complexity of Ω(n + k) with
k denoting the number of existing
...
intersections. Without proof we

give the lower bound for compari-
son based line segment intersection
as Ω(n · log(n) + k). Obviously Figure 11.2: Layout with Θ(n2 ) intersections
we can construct cases containing
Θ(n2 ) (see Figure 11.2) intersections, thus resulting in a running time of Θ(n2 ) for
any algorithm. Still most realistic inputs have k = O(n).
To improve our algorithm, we can take into account the spacial locality. Ob-
viously, pairs of line segments, whose y-coordinates are far apart, cannot have a
common point. Imagine all segments to be projected onto the y-axis: we only need
to test line segments for intersection, whose y-intervals interlap. To find those pairs
of lines, we imagine a horizontal line ` sweeping downwards over the plane. This
line is called sweep-line and algorithms using this principle we call plane-sweep
algorithms. While we sweep down the plane, we keep track of all vertical segments
intersecting `. This yields the following algorithm:
T =hi : SortedSequence of Segment

invariant T stores the vertical segments intersecting `
Q:= sort(h(y, s) : ∃hor. segment s at y or ∃vert. segment s starting/ending at yi)
// tie breaking: vert. starting events first, vert. finishing events last
foreach (y, s) ∈ Q in descending order do
if s is a vertical segment and starts at y then T.insert(s)
else if s is a vertical segment and ends at y then T.remove(s)
else // we have a horizontal segment s = (x1 , y)(x2 , y)
foreach t = (x, y1 )(x, y2 ) ∈ T with x ∈ [x1 , x2 ] do
output {s, t}
handle horizontal segments on ` // interval intersection problem
Keeping a sorted range T , the at most n insert and remove operations can be
performed within O(log(n)) each. Finding the boundaries of the interval containing
the intersections can be done in O(log(n)), the output of the intersections in O(ks )
with ks denoting the number of intersections Pon the horizontal segment s. Overall
this yields a running time of O(n · log(n) + s ks ) = O(n · log(n) + k).
11.4.2 A more Complex Case

Now we will focus on a more
complex case. We will now al-
low a more arbitrary positioning of
the line segments while excluding
some special cases. For now we do
not allow horizontal line segments,
overlapping segments and we do
not allow one intersection to be a
common intersection of more than
two line segments.
To extend our basic sweep-line
algorithm for orthogonal line seg-
ment intersection, we introduce so Figure 11.3: Events during the execution of
called events. These events are the sweep-line algorithm
illustrated in Figure 11.3. These
events include all changes that affect our sweep-line algorithm in any way, to be
exact they consist of segment starting points, segment end points and intersections.
Using this definition, we only have to test for intersection of elements that are next
to each other in T at some event.
11.4. LINE SEGMENT INTERSECTION - CASE STUDY 103
Lemma 26
Let si and sj be two non horizontal segments whose interiors intersect in a singe
point p, and assume there is no third segment passing through p. Then there is an
event point above p where si and sj become adjacent and are tested for intersection.
Proof.
Choose ` in such a way that there is no event located on ` and no event between `
and the horizontal line through p. Then si and sj are adjacent in T for this position
of `. Since T = hi at the beginning of the algorithm, there exists an event where si
and sj become adjacent in T and are tested for intersection.
This observations yields an improved version of our algorithm. We define two
functions to aid our algorithm description. The first is called findNewEvent and
calculates whether there is an intersection event between two given line segments.
Function findNewEvent(s, t) // O(1 + log n)

if ∗s and ∗t intersect at y 0 < y then
Q.insert((y 0 , intersection, (s, t)))
Due to the insert into Q, findNewEvent runs in O(1 + log(n)). The second
function performs the event processing step. We distinguish between three version
of handleEvent for the three different types of events.
Function handleEvent(y, start, s, T, Q) // n×

h:= T.insert(s) // O(log n)
prev:= pred(h) // O(1)
next:= succ(h) // O(1)
findNewEvent(prev, h)
findNewEvent(h, next)
Function handleEvent(y, finish, s, T, Q) // n×

h:= T.locate(s) // O(log n)
prev:= pred(h) // O(1)
next:= succ(h) // O(1)
T.remove(s) // O(log n)
findNewEvent(prev, next)
Function handleEvent(y, intersection, (a, b), T, Q)// k×

output (∗s ∩ ∗t) // O(1)
T.swap(a, b) // O(1)
prev:= pred(b) // O(1) a b
next:= succ(a) // O(1)
findNewEvent(prev, b)
findNewEvent(a, next)
Due to the n line segments, we have to handle n start and n finish events. The
number of intersection events is linear in the number of intersections as can be
proven by an interpretation as a planar graph and by using eulers formula. For a
detailed analysis we refer to [? ].
T =hi : SortedSequence of Segment

invariant T stores the relative order of the seqments intersecting `
Q : MaxPriorityQueue
n o
Q:= Q ∪ (y, start, s) : s = (x, y)(x0 , y 0 ) ∈ S, y ≤ y 0 // O(n log n)
n o
Q:= Q ∪ (y, finish, s) : s = (x, y)(x0 , y 0 ) ∈ S, y ≥ y 0 // O(n log n)
while Q 6= ∅ do
(y, type, s):= Q.deleteMin // O((n + k) log n)
handleEvent(y, type, s, T, Q)
The initialization of our algorithm performs in O(n · log(n)), the event loop in
O((n + k) · log(n)). Thus the overall running time is in O((n + k) · log(n)).
11.4.3 The Arbitrary Case

To allow horizontal lines, we have to introduce following order on Q:
(x, y) ≺ (x0 , y 0 ) ⇔ y > y 0 ∨ y = y 0 ∧ x < x0
This can be interpreted as a sweep-line with an infinitely small slope. With this we
can once more adapt our algorithm:
To finally also allow overlapping segments, we can compute the line g(s) for
any segment s ∈ S that contains s. Sorting S along g(s) enables us to solve a
1-dimensional overlapping problem for each segment.
11.5. CONVEX HULL 105
Function handleEvent(p = (x, y))

U := segments starting at p // from Q
C:= segments with p in their interior // from T
L:= segments finishing at p // from Q
if |U | + |C| + |L| ≥ 2 then report intersection @ p
T.remove(L ∪ C)
T.insert(C ∪ U ) such that order just below p is correct
if U ∪ C = ∅ then
findNewEvent(T.findPred(p), T.findSucc(p), p)
else
findNewEvent(T.findPred(p), T.findLeftmost(p), p)
findNewEvent(T.findRightmost(p), T.findSucc(p), p)
Function findNewEvent(s, t, p)
if s and t intersect at a point p0 p then
if p0 6∈ Q then Q.insert(p0 )
11.4.4 Further Improvements

The algorithm as specified has one problem. The space consumption can become
quite large, as any intersection is kept in the queue, even though the segments may
not be adjacent at the moment. This yields a space consumption of Θ(n + k) which
can be quite large for a large k. But we can improve this to O(n) by a simple
modification. If we delete the intersection events from the queue for segments not
longer adjacent, the space consumption obviously falls to O(n). The algorithm
also does still report all intersections correctly, as all segments become adjacent
again just before they intersect.
11.5 Convex Hull

In this section we will discuss the problem of
the convex hull. For simplicity we focus on the
2-dimensional space R2 . pn
The formal definition of the Convex-Hull- p1

Problem in R2 is given as follows:
Given a set of points P = p1 , . . . , pn ∈
R2 . Calculate a convex polygon C, using only Figure 11.4: Upper/Lower Hull
points in P with ∀p ∈ P : p contained in C.
We discuss an algorithm running in O(sort(n)) time. First we assume the points
are given in lexicographical order (i.e. a point pi comes before pj if pi .x < pj .x or
pi .x = pj .x and pi .y < pj .y. If the points are not given in this order, we can sort
them in O(sort(n)) additional time.
For this algorithm, we first split the set of points along the line segment p1 pn .
Now the algorithm can be implemented as two convex hull calculations that process
the lower-hull, which represents the hull of points below p1 pn , and the upper-hull
separately (illustrated in Figure 11.4). In the following we will focus on the
upper-hull. The algorithm for the lower-hull can be formulated equivalently.
11.5.1 Graham Scan

Before we can actually discuss the algorithm, we have to give two definitions.
Definition 14
For a set of consecutive line segments, given as a set of points P = p1 , . . . , pn , we
speak of a right turn at i if pi+1 is to the right of the line formed by pi−1 and pi . In
the same matter we speak of a left turn, if pi+1 is to the left of the line formed by
pi−1 and pi .
Looking only at the upper-hull, we observe that the convex-hull consists of a
set of line segments performing only right-turns from p1 to pn and the line segment
p1 pn . Thus, we can give the following algorithm to compute the upper-hull:
Function upperHull(p1 , . . . , pn )
L=hpn , p1 , p2 i : Stack of Point
invariant L is the upper hull of hpn , p1 , . . . , pi i
for i := 3 to n do
while ¬rightTurn(L.secondButlast, L.last, pi ) do
L.pop
L:= L ◦ hpi i
return L
To prove our algorithm correct, we will prove the invariant.

Proof.
Before the for-loop, the set L contains only three vertices. Thus, L is a convex
hull by definition. Now suppose that L forms the convex-hull of hp1 , . . . , pi−1 i
and consider the addition of pi . After the execution of the while-loop, the points
in L form a chain of right turns. To verify our assumption, we have to show that
all points ∈ {p1 , . . . , pi−1 } not contained in L are below this chain. Due to the
induction hypothesis, we know that a point below the new chain could only exist
in the vertical slab between pi and pi − 1 since the new chain lies above the old
chain. Because of the lexicographical order, no such point exists.
11.6. SMALLEST ENCLOSING CIRCLE 107
The running time of the algorithm upper-hull itself lies in O(n). This is, as the
for loop obviously runs in O(n) time. The additional work done by the while loop
only adds an additional O(n) in total, as each point can only be deleted once from
L. Thus the overall running time is dominated by the sorting and the algorithm
performs in O(n · log(n)).
In three-dimensional space, the number of vertices does not directly imply a linear
complexity of the convex hull. Still one can show that the number of edges and the
number of faces are both in O(n). We refer to [? ] for a proof. For even higher
dimensions, the output complexity rises to O nbd/2c .
11.6 Smallest Enclosing Circle

In this part we will go on to discuss the smallest enclosing circle problem which is
defined as follows:
Given a set of points P = {p1 , . . . , pn } ∈ Rd , calculate a sphere K with
minimal radius so that P ⊆ K. We will discuss a simple algorithm with expected
running time O(n) as introduced by Welzl in 1991. An illustration of the problem
can be found in Figure 11.5.
Figure 11.5: Enclosing Balls in 2-dimensional space
The algorithm is defined by a recursive function:
Function smallestEnclosingBallWithPoints(P, Q)
if |P | = 1 ∨ |Q| = d + 1 then return ball(Q)
pick random x ∈ P
B:= smallestEnclosingBallWithPoints(P \ {x} , Q)
if x ∈ B then return B
return smallestEnclosingBallWithPoints(P \ {x} , Q ∪ {x})
Note the base of our recursion: a sphere for only one point can be given as
the point itself without any radius. Also a set of d + 1 vertices uniquely defines a
sphere in d-dimensional space. In Q we store a set of border nodes that form the
current sEB.
To prove our algorithm correct we have to prove two claims.
The first states:
Lemma 27
Smallest Enclosing Ball With Points (sEB) are non-ambiguous.
Proof.
Assume: ∃ sEB’s B1 6= B2 . Since B1 and B2 are both sEB’s P ⊆ B1 ∧ P ⊆ B2
holds. This directly induces P ⊆ B1 ∩ B2 ⊆ sEB( B1 ∩ B2 ) = B and radius(B)
< radius(B1 ). This cannot hold, since B1 is a sEB.
Lemma 28
If a point x is not in the sEB of P \ {x} , Q, then x is on the border of sEB(P ).
Proof.
We show the contraposition: Assume x is not on the border of sEB(P ). Then
sEB(P ) = sEB(P \ {x}) = B, due to the non-ambiguous sEBs.
11.6.1 Analysis
We have given a probabilistic algorithm for the smallest enclosing ball problem
and proven it correct. In this section we want to analyse the expected running time
for the algorithm.
We start by defining the base of our recurrence. As already mentioned above,
our recursion T (p, q) ends for T (p, d + 1) = T (1, p) = 0.
For T (p, q) following equation holds:
T (p, q) ≤ 1 + T (p − 1, q) + P[x ∈ / B] · T (p, q + 1) as we only have to execute the
second recursive call, whenever x ∈ / B. Now let us inspect P[x ∈ / B] for a fixed
subset of points p1 , . . . , pi and its sEB Q. If we remove a point pj , when does Q
change? This happens only if we remove one of the d+1 points that uniquely define
Q. If there are more than d + 1 points on the boundary of Q, the smallest enclosing
ball cannot change after the removal of one point. Thus, since we already know q
points to be on the boundary of Q, the number of points that could change the sEB
is given by d + 1 − q. Therefore T (p, q) ≤ 1 + T (p − 1, q) + d+1−q p
· T (p, q + 1)
holds.
Using this equation, we can calculate T (p, d + 1) for different values of d. For
example d = 1 would yield an expected running time of 3 · n whereas d = 2 would
result in a running time of 10 · n. In general T (p, 0) ≤ d! · n holds.
11.7. RANGE SEARCH 109
11.7 Range Search

The most general version of range searching problem is given as follows: Given
a set of geometric objects P and a geometric range Q, determine which objects
of P Q contains. In short we will write P ∩ Q for such a range reporting query.
Furthermore we might be interested only in the number of elements contained in
Q. We call this a range counting query and write |P ∩ Q|. In most applications
of range search, we get a lot of querys for a given set of points. Thus we exploit
precomputation to answer the queries as fast as possible. Intuitively we would
expect the precomputation the be possible in O(n log n) and consume only O(n)
space. Counting queries should be possible to answer in O(log n) while we
would expect range reporting queries to be possible in O(k + log n) or at least in
O(k log n) with k denoting the number of contained elements.
11.7.1 1D Range Search

In the 1-dimensional case, we consider a list of numbers. To answer range reporting
and range counting queries, we utilize a binary range tree, as depicted in Figure 11.6.
The tree consists of a doubly linked list containing the given elements in sorted
order. On top we keep a tree structure, containing the splitting element for the
underlying range as key.
The construction is obviously bounded by sorting the elements in O(n log n)
time. Range queries traverse the tree for both ends of the range, giving the start
and the end of the list, containing all elements contained in the range. This can
be done in O(log n). Afterwards, the range can then be reported in O(k). Range
counting queries can be answered in O(log n) if we keep the number of elements
contained in the range covered by a given tree node up-to-date.
[6 ,15]? 17
3 13
2 5 11 19
2 3 5 7 11 13 17 19 00
Figure 11.6: A 1D range tree.

11.7.2 2D Range Search

Range queries in 2-dimensional space have their application in geocoding or
collision detection. A possible representation for answering those kind of queries
fast is the Wavelet Tree introduced by Chazelle in 1988. The wavelet tree utilizes
bitvectors to represent itself.
Be reminded that in the following, coordinates are represented by their ranks
in each dimension respectively, i. e. we have only integers. On each level of the
tree, these ranks are newly calculated, i. e. compacted. Futhermore, the sequence
of points is always kept sorted in y-direction.
Definition 15
For a given range of elements X = hx1 , . . . , xn i, a Wavelet Tree is recursively
defined as follows. The tree keeps a bitvector b with b[i] = 1 iff xi > bn/2c.
Furthermore, it contains two children, each recursivly defined as the Wavelet Trees
of hxi : xi ≤ bn/2ci as the left child and hxi − bn/2c : xi > bn/2ci as the right
child respectively. This recursion is stopped as soon as n ≤ n0 , where n0 is a
tuning parameter. In this case, the sequence is stored directly. A sample wavelet
tree is shown in Figure 11.7.
0
1
1
1
0
0
1
0
1 1
1 0
0 0
0 1
21,12 21,12 22,11 22,11
Figure 11.7: Example Wavelet Tree
Definition 16
The Dominance Counting Problem is defined as follows: Given a sequence of
elements X = hx1 , . . . , xn i as well as an pair (x, y),
find |{xi : xi ∈ X ∧ xi ≥ x ∧ i ≥ y}|. (dominance criterion)
Using Wavelet Trees, we can solve this problem with Algorithm 11.8. In every
recursion step, the number of elements yr less or equal to y is calculated using the
rank function of the bitvector. How this can be accomplished will be discussed later.
If the desired value of x is found in the left child, we know of at least dn/2e − yr
elements to fulfill the dominance criterion. The rest can be calculated recursively
with (x, y − yr ) as input for the left child, as all (xi , yi ) with yi < y ∧ xi > x are
not considered in the left child. The recursion for x > bn/2c can be explained in a
similar matter. For an illustration, the reader can turn to Figure 11.9.
If we implement the rank to work in constant time (see below), the algorithm
runs in O(log n) as we only follow one recursive path and have a recursion depth
of O(log n).
Function intDominanceCount(x, y) // |[x, n] × [y, n] ∩ P |

if n ≤ n0 then return |[x, n] × [y, n] ∩ P | // brute force
yr := b.rank(y) // Number of els ≤ y in right half
if x ≤ bn/2c then
return `.intDominanceCount(x, y − yr ) + dn/2e − yr
else
return r.intDominanceCount(x − bn/2c , yr )
Figure 11.8: Algorithm intDominanceCount
0
1
1
1
0 4−1=3
3 0
1 1
0
2
1 2−0=2 1
1 0
2 0 0
00 1
2
21,12 21,12 22,11 22,11
Figure 11.9: Example query for dominance counting using a Wavelet Tree
Using Algorithm 11.10, we can also perform range counting queries on integer
values. The algorithm is illustrated in Figure 11.11.
Function intRangeCount([x, x0 ] × [y, y 0 ])

return
intDominanceCount(x, y)−
intDominanceCount(x0 , y)−
intDominanceCount(x, y 0 )+
intDominanceCount(x0 , y 0 )
Figure 11.10: Algorithm: integer range counting

n 11111111111111
00000000000000
0000000
1111111
00000000000000
11111111111111
0000000
1111111
11111111111111
00000000000000
0000000
1111111
00000000000000
11111111111111
0000000
1111111
00000000000000
11111111111111
0000000
1111111
y’ 11111111111111
00000000000000
0000000
1111111
00000000000000
11111111111111
0000000
1111111
00000000000000
11111111111111
0000000
1111111
00000000000000
11111111111111
0000000
1111111
00000000000000
11111111111111
0000000
1111111
00000000000000
11111111111111
0000000
1111111
y 11111111111111
00000000000000
0000000
1111111
x11111111111111
00000000000000
0000000
1111111
x’ n
Figure 11.11: Integer Range Counting using Integer Dominance Counting
Wavelet Tree Dominance Reporting Query

To implement the Dominance Reporting Query, we look at a subproblem first.
Definition 17
The One Sided Reporting problem is given as: Given a set P = {(x1 , y1 ), . . . (xn , yn )}
and a value y, find P ∩ [1, n] × [y, n].
The One Sided Reporting problem can be solved by Algorithm 11.12. It is

based on the same principles as Algorithm 11.8. The only difference is that we
have to consider all elements above the given y value. Thus we have to perform a
recursion step in both children.
Function oneSidedReporting(y) // |[1, n] × [y, n] ∩ P |

if n ≤ n0 then return |[1, n] × [y, n] ∩ P | // brute force
R:= ∅
if yr < n2 then R:= R ∪ r.oneSidedReporting(yr )
if y − yr < n2 then R:= R ∪ `.oneSidedReporting(y − yr )
return R
Figure 11.12: Algorithm for One Sided Reporting in a Wavelet Tree
Having solved the One Sided Reporting problem, we are ready to report all
elements dominating an (x, y) pair with Algorithm 11.13:
Function intDominanceReporting(x, y) // |[x, n] × [y, n] ∩ P |

if n ≤ n0 then return |[x, n] × [y, n] ∩ P | // brute force
R:= ∅ // Result
if x ≤ bn/2c then // Both halves interesting
n
if y − yr < 2 then R:= R ∪ `.intDominanceReporting(x, y − yr )
if yr < n2 then R:= R ∪ r.oneSidedReporting(yr )
else if yr < n2 then R:= R ∪ r.intDominanceReporting(x − bn/2c , yr )
return R
Figure 11.13: Algorithm for Dominance Reporting using a Wavelet Tree
A generalization to the range reporting query can be achieved. The range

reporting query, as all other queries before, is defined recursively. We report
elements in the moment that all elements above a given y value are dominant. The
rest of the algorithm is identical to the range counting algorithm.
We can adapt the one sided reporting to work with a minimal and a maximal
rank and report the elements within a range. The rest of the adaption is trivial
following the Illustration in Figure 11.14.
2D Range Search on floating point numbers

Now consider the following problem:
Definition 18
The 2D-Range Search problem is defined as follows.
Given a set of data P = {p1 , . . . , pn } ⊆ R2 , and an axis parallel rectangle
Q = [x, x0 ] × [y, y 0 ]
1. find P ∩ Q (range query)
2. find |P ∩ Q| (range counting)
We can reduce the general problem onto the integer problem on [1, n] × [1, n].
This can be achieved by replacing each pair (x, y) ∈ P with (rank(x, xi ), rank(y, yi )),
where rank(x, xi ) denotes the rank of x in a sorted range of the {xi |(xi , y) ∈
P for y ∈ R}. This replacement can be achieved by sorting the xi and yi respec-
tivly im combination with a binary search for every xi , yi in the resulting ranges.
Thus we can achieve a general range query in O(n log n), using the integer range
counting methods discussed above.
...
...
Figure 11.14: Illustration of general range reporting Query
11.7.3 Bitvectors
Bitvectors, as used in the range query algorithms discussed above, store a set of
bits and offer a set of distinct operations. For our use case, the interesting operation
is the rank(j) operation.
To support the rank operation in constant time, we precompute the ranks of
every iB-th element, were B can be seen as a tuning parameter. Than rank(j)
can be computed as the sum of the precomputed value j mod B and rank(v[j
mod B], j]). If we chose B = Θ(log n), we use O(n) time and O(n) bits. This
process is illustrated in Figure 11.15.
The constant lookup can be either realized using a lookup table or by the
instruction population count which on some systems (e. g. SSE 4.2) counts the
number of set bits in a machine word.
In a similar manner, we can calculate the position of the i-th set bit b.select(i)
in constant time. Furthermore, it is possible to reach the information theoretical
optimal space of n + o(n) bits.
0 4 8 11
10010011 10001110 00001110 10100010
rank=8+2=10
Figure 11.15: Illustration of the rank function on a bit vector

Bibliography
[1] A. Aggarwal and J. S. Vitter. The input/output complexity of sorting and

related problems. Communications of the ACM, 31(9):1116–1127, 1988.
[2] A. V. Aho, D. S. Johnson, R. M. Karp, S. Rao Kosaraju, C. C. McGeoch,

C. H. Papadimitriou, and P. Pevzner. Emerging opportunities for theoretical
computer science. SIGACT News, 28(3):65–74, 1997.
[3] Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin. Network

Flows: Theory, Algorithms, and Applications. Prentice Hall, 1993.
[4] E. Althaus and K. Mehlhorn. Traveling salesman-based curve reconstruction

in polynomial time. SIAM Journal on Computing, 31(1):27–66, 2002.
[5] D. Applegate, R. Bixby, V. Chvátal, and W. Cook. Implementing the Dantzig-

Fulkerson-Johnson algorithm for large traveling salesman problems. Math.
Programming, 97(1–2):91–153, 2003.
[6] H. Bast, S. Funke, P. Sanders, and D. Schultes. Fast routing in road networks
with transit nodes. Science, 316(5824):566, 2007.
[7] R. Beier and B. Vöcking. Typical properties of winners and losers in discrete
optimization. In 36th ACM Symposium on the Theory of Computing, pages
343–352, 2004.
[8] E. Berberich, A. Eigenwillig, M. Hemmer, S. Hert, L. Kettner, K. Mehlhorn,

J. Reichel, S. Schmitt, E. Schömer, and N. Wolpert. EXACUS—efficient and
exact algorithms for curves and surfaces. In 13th ESA, volume 3669 of LNCS,
pages 155–166, 2005.
[9] T. Beth and M. Clausen, editors. Applied Algebra, Algebraic Algorithms and
Error-Correcting Codes, 4th International Conference, AAECC-4, Karlsruhe,
FRG, September 23-26, 1986, Proceedings, volume 307 of Lecture Notes in
Computer Science. Springer, 1988.
115
116 BIBLIOGRAPHY
[10] A. Borodin and R. El-Yaniv. Online Computation and Competitive Analysis.

Cambridge University Press, Cambridge, 1998.
[11] G. S. Brodal, R. Fagerberg, and K. Vinther. Engineering a cache-oblivious

sorting algorithm. In 6th Workshop on Algorithm Engineering and Experi-
ments, 2004.
[12] C. Burnikel, J. Könemann, K. Mehlhorn, S. Näher, S. Schirra, and C. Uhrig.

Exact geometric computation in leda. In SCG ’95: 11th annual symposium
on Computational geometry, pages 418–419, New York, NY, USA, 1995.
ACM.
[13] B. V. Cherkassky and A. V. Goldberg. On implementing push-relabel method

for the maximum flow problem. In IPCO: 4th Integer Programming and
Combinatorial Optimization Conference, 1995.
[14] D. Delling, P. Sanders, D. Schultes, and D. Wagner. Engineering route plan-

ning algorithms. In Algorithmics of Large and Complex Networks, volume
5515 of LNCS State-of-the-Art Survey, pages 117–139. Springer, 2009.
[15] R. Dementiev, L. Kettner, and P. Sanders. STXXL: Standard Template Library

for XXL data sets. Software Practice & Experience, 38(6):589–637, 2008.
[16] R. Dementiev and P. Sanders. Asynchronous parallel disk sorting. In 15th

ACM Symposium on Parallelism in Algorithms and Architectures, pages
138–148, San Diego, 2003.
[17] C. Demetrescu, I. Finocchi, G. F., and Italiano. Algorithm engineering.

Bulletin of the EATCS, 79:48–63, 2003.
[18] D.E. Drake and S. Hougardy. A simple approximation algorithm for the
weighted matching problem. Information Processing Letters, 85(4):211–213,
2003.
[19] R. Fleischer, B. Moret, and E. Meineche Schmidt, editors. Experimental Al-

gorithmics From Algorithm Design to Robust and Efficient Software, volume
2547 of LNCS. Springer, 2002.
[20] V. Kumar G. Karypis. Multilevel k-way partitioning scheme for irregular

graph. J. Parallel Distrib. Comput., 48(1), 1998.
[21] A. V. Goldberg and S. Rao. Beyond the flow decomposition barrier. Journal
of the ACM, 45(5):1–15, September 1998.
BIBLIOGRAPHY 117
[22] C. Gutwenger and P. Mutzel. A linear time implementation of SPQR trees.

In J. Marks, editor, Graph Drawing, volume 1984 of LNCS, pages 77–90.
Springer-Verlag, 2001.
[23] J. Hopcroft and R. E. Tarjan. Efficient planarity testing. J. of the ACM,

21(4):549–568, 1974.
[24] R. E. Tarjan J. E. Hopcroft. Dividing a graph into triconnected components.

SIAM J. Comput., 2(3):135–158, 1973.
[25] D. S. Johnson. A theoretician’s guide to the experimental analysis of al-

gorithms. In M. Goldwasser, D. S. Johnson, and C. C. McGeoch, editors,
Proceedings of the 5th and 6th DIMACS Implementation Challenges. Ameri-
can Mathematical Society, 2002.
[26] M. Jünger, S. Leipert, and P. Mutzel. A note on computing a maximal planar

subgraph using PQ-trees. IEEE Transactions on Computer-Aided Design,
17(7):609–612, 1998.
[27] M. Jünger and P. Mutzel. 2-layer straightline crossing minimization: Perfor-

mance of exact and heuristic algorithms. Journal of Graph Algorithms and
Applications (JGAA), 1(1):1–25, 1997.
[28] D. E. Knuth. The Art of Computer Programming—Sorting and Searching,

volume 3. Addison Wesley, 2nd edition, 1998.
[29] V. Kumar, A. Grama, A. Gupta, and G. Karypis. Introduction to Parallel

Computing. Addison Wesley, 2003.
[30] C. C. McGeoch, D. Precup, and P. R. Cohen. How to find big-oh in your data
set (and how not to). In Advances in Intelligent Data Analysis, number 1280
in LNCS, pages 41–52, 1997.
[31] C.C. McGeoch and B. M. E. Moret. How to present a paper on experimental

work with algorithms. SIGACT News, 30(4):85–90, 1999.
[32] K. Mehlhorn. Data Structures and Algorithms, Vol. I — Sorting and Searching.
EATCS Monographs on Theoretical CS. Springer-Verlag, Berlin/Heidelberg,
Germany, 1984.
[33] K. Mehlhorn and P. Mutzel. On the embedding phase of the Hopcroft and
Tarjan planarity testing algorithm. Algorithmica, 16(2):233–242, 1996.
[34] K. Mehlhorn and S. Näher. The LEDA Platform of Combinatorial and

Geometric Computing. Cambridge University Press, 1999.
118 BIBLIOGRAPHY
[35] K. Mehlhorn and P. Sanders. Scanning multiple sequences via cache memory.
Algorithmica, 35(1):75–93, 2003.
[36] K. Mehlhorn and G. Schäfer. Implementation of weighted matchings in
general graphs: The power of data structure. ACM Journal of Experimental
Algorithmics, 7, 2002.
[37] Kurt Mehlhorn and Stefan Näher. LEDA, A platform for Combinatorial and
Geometric Computing. Cambridge University Press, 1999.
[38] Kurt Mehlhorn and Peter Sanders. Algorithms and Data Structures: The
Basic Toolbox. Springer, 2008.
[39] U. Meyer, P. Sanders, and J. Sibeyn, editors. Algorithms for Memory Hierar-
chies, volume 2625 of LNCS Tutorial. Springer, 2003.
[40] U. Meyer, P. Sanders, and J. Sibeyn, editors. Algorithms for Memory Hierar-
chies, volume 2625 of LNCS Tutorial. Springer, 2003.
[41] B. M. E. Moret. Towards a discipline of experimental algorithmics. In 5th
DIMACS Challenge, DIMACS Monograph Series, 2000.
[42] J. von Neumann. First draft of a report on the EDVAC. Technical report,
University of Pennsylvania, 1945.
[43] S. Pettie and P. Sanders. A simpler linear time 2/3-approximation for max-
imum weight matching. Information Processing Letters, 91(6):271–276,
2004.
[44] T. Polzin and S. V. Daneshmand. Extending reduction techniques for the
Steiner tree problem. In 10th Eur. Symp. on Algorithms, volume 2461 of
LNCS, pages 795–807. Springer, 2002.
[45] K. R. Popper. Logik der Forschung. Springer, 1934. English Translation: The
Logic of Scientific Discovery, Hutchinson, 1959.
[46] P. Sanders. Fast priority queues for cached memory. In ALENEX ’99, Work-
shop on Algorithm Engineering and Experimentation, number 1619 in LNCS,
pages 312–327. Springer, 1999.
[47] P. Sanders. Fast priority queues for cached memory. ACM Journal of Experi-
mental Algorithmics, 5, 2000.
[48] P. Sanders. Algorithm engineering – an attempt at a definition. In Efficient
Algorithms, volume 5760 of Lecture Notes in Computer Science, pages 321–
340. Springer, 2009.
BIBLIOGRAPHY 119
[49] P. Sanders, S. Egner, and J. Korst. Fast concurrent access to parallel disks. In
11th ACM-SIAM Symposium on Discrete Algorithms, pages 849–858, 2000.
[50] P. Sanders and D. Schultes. Highway hierarchies hasten exact shortest path
queries. In 13th European Symposium on Algorithms (ESA), volume 3669 of
LNCS, pages 568–597. Springer, 2005.
[51] P. Sanders and S. Winkel. Super scalar sample sort. In 12th European
Symposium on Algorithms, volume 3221 of LNCS, pages 784–796. Springer,
2004.
[52] J. Singler, P. Sanders, and F. Putze. MCSTL: The multi-core standard template
library. In 13th International Euro-Par Conference, volume 4641 of LNCS,
pages 682–694. Springer, 2007.
[53] D. Spielman and S.-H. Teng. Smoothed analysis of algorithms: why the
simplex algorithm usually takes polynomial time. In 33rd ACM Symposium
on Theory of Computing, pages 296–305, 2001.
[54] D. Gollman T. Beth. Algorithm engineering for public key algorithms. IEEE
Journal on Selected Areas in Communications, 7(4):458–466, 1989.
[55] J. C. Venter, M. D. Adams, and E. W. Myers et al. The sequence of the human
genome. Science, 291(5507):1304–1351, 2001.
[56] J. S. Vitter and E. A. M. Shriver. Algorithms for parallel memory, I: Two

level memories. Algorithmica, 12(2/3):110–147, 1994.
[57] I. Wegener. Simulated annealing beats metropolis in combinatorial opti-

mization. In 32nd International Colloquium on Automata, Languages and
Programming, pages 589–601, 2005.

Adsii PDF

Uploaded by

Copyright:

Available Formats

Adsii PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Adsii PDF

Uploaded by

Copyright:

Available Formats

Algorithms and Data Structures II

Last modified: January 4, 2021

5 Maximum Flows and Matchings 15

6 Generic Approaches to Optimization 43

Figure 1.1: Algorithm engineering as a cycle of design, analysis, implementation,

This is a shortened version of [48]. Algorithms and data structures are at

1.1 A Brief “History” of Algorithm Engineering

P. Italiano organized the “Workshop on Algorithm Engineering” in 1997 and

1.7 Algorithm Libraries

for advanced models of computation like external computing (STXXL [15]) or

practice. We suspect that this often leads to unrealistic results in experimental

constant factors, smoothed analysis,. . . ). We also have no intention of criticizing

Maximum Flows and Matchings

5.1 Algorithms 1956–now

Year Author Running time

5.2 Augmenting Path Algorithms

5.2.1 The Ford Fulkerson Algorithm

Function FFMaxFlow(G = (V, E), s, t, c : E → N) : E → N

Figure 5.3: The Ford Fulkerson Augmenting Path algorithm.

Procedure augmentPath(G = (V, E), s, t, c : E → N, f : E → N, path p in Gf )

Figure 5.4: The Ford Fulkerson subroutine for augmenting a path.

Figure 5.5: An example execution of the augmenting path algorithm.

We now proof that the Ford Fulkerson algorithm is correct.

A cut S is called saturated if f (e) = ce for all e = (u, v), u ∈ S, v ∈ T and

If S is saturated then val(f ) = c(S).

For a saturated cut, the inequalities are equalities.

1. If t ∈ S then f is not maximum.

2. If t 6∈ S then S is a saturated cut and f is maximum.

1. Let p be any simple path from s to t in Gf and let ∆f be the minimum

Then f 0 is a flow and val(f 0 ) =val(f ) + ∆f . Note that f 0 is well-defined

2. There is no edge (v, w) in Gf with v ∈ S and w ∈ T . Hence, f (e) = c(e)

Theorem (Elias/Feinstein/Shannon, Ford/Fulkerson 1956)

5.2.2 Two Bad Examples for Ford Fulkerson

100 100 99 100 99 99

5.2.3 Dinic’s Algorithm

Definition 4 (Distance Labels)

Function DinitzMaxFlow(G = (V, E), s, t, c : E → N) : E → N

Figure 5.8: Pseudocode for Dinic’s algorithm.

Function blockingFlow(H = (V, E)) : E → N

running time of this algorithm is basically #extends + #retreats + n · #breakthroughs .

Figure 5.10: An example run of the blocking flow subroutine.

Putting these things together √

5.2.4 Disadvantage of Augmenting Path Algorithms

Figure 5.11: Drawback of the augmenting path algorithm.

5.3 Preflow-Push Algorithms

The preflow-push algorithms maintain a preflow at each intermediate stage. In a

Procedure push(e = (v, w), δ)

Figure 5.12: The basic push operation.

A push of δ units from node v to w decreases both excess(v) and cf(v,w) by δ

Figure 5.13: The generic preflow-push algorithm.

∀ active nodes v : excess(v) > 0 ⇒ ∃ path hv, . . . , si ∈ Gf .

arbitrary node and edge selection rules. ( arbitrary preflow-push).

2. A non-saturating push decreases φ by at least one.

3. A relabeling increases φ by one.

4. A saturating push increases φ by at most 2n.

|{w : d(w) ≤ d(v)}|

as potential function. We split the execution into phases. We define a phase to

2. Φ ≥ 0 always, Φ ≤ n2 /K initially (obvious)

3. a relabel or saturating push increases Φ by at most n/K.