Computing fixed points of monotone

ℓ p-contractions

Hugo Arbones

NI VER
U S
E

IT
TH

Y
O F

H
G

E
R

D I U
N B

4th Year Project Report

Computer Science and Mathematics
School of Informatics
University of Edinburgh
2026
 Abstract
Fixed points provide a framework for encoding solutions to a wide range of problems,
including equilibria in game theory. Results such as Tarski’s and Banach’s Fixed Point
Theorems guarantee the existence of fixed points for monotone functions and contraction
mappings, respectively. Still, no efficient algorithms are known for computing these
fixed points in general.
This dissertation studies the problem of computing fixed points of functions that are
both monotone and ℓ p -contractive, a setting that is less studied than either property in
isolation.
We first investigate approaches for fixed points of ℓ p -contractions, introduce a decom-
position method that breaks high-dimensional problems into smaller sub-problems,
and develop a heuristic that accelerates existing algorithms when the monotonicity and
ℓ p -contraction properties are combined.

i
 Research Ethics Approval
This project was planned in accordance with the Informatics Research Ethics policy. It
did not involve any aspects that required approval from the Informatics Research Ethics
committee.

Declaration
I declare that this thesis was composed by myself, that the work contained herein is my
own except where explicitly stated otherwise in the text, and that this work has not been
submitted for any other degree or professional qualification except as specified.

(Hugo Arbones)

ii
 Acknowledgements
Thank you to my supervisor, Kousha Etessami, and to Lidl’s instant coffee, which
powered every word of this dissertation.

iii
 Table of Contents

1 Introduction 1

2 Background 2
2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1.1 Contraction Mapping Fixed Points . . . . . . . . . . . . . . . 3
2.1.2 Monotone Functions and Lattice-Theoretic Fixed Points . . . 6
2.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 The ARRIVAL Problem . . . . . . . . . . . . . . . . . . . . 9
2.2.2 Simple Stochastic Games . . . . . . . . . . . . . . . . . . . . 11

3 Theoretical Advancements 14
3.1 Algorithms for ℓ p -contractions . . . . . . . . . . . . . . . . . . . . . 14
3.1.1 Haslebacher’s Algorithm . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Decomposition Algorithm . . . . . . . . . . . . . . . . . . . 28
3.2 Algorithms for monotone ℓ p contractions . . . . . . . . . . . . . . . 41
3.2.1 Bisection Algorithm Description . . . . . . . . . . . . . . . . 41
3.2.2 Applications to Existing Algorithms . . . . . . . . . . . . . . 44

4 Implementation and Results 46

5 Conclusion 49

Bibliography 52

A Figures and Examples 56

iv
 Chapter 1

Introduction

Fixed points provide a way to encode solutions to many problems, from differential
equations [CL55, Lin94] to equilibria in game theory [Sha53]. Two main theorems
guarantee fixed points in different contexts: Tarski’s Fixed Point Theorem ensures
the existence of at least one fixed point of a monotone function on a complete lattice
[Tar55]. Similarly, Banach’s Fixed Point Theorem guarantees a unique fixed point of a
contraction map on a complete metric space [Ban22].
The problem of finding such fixed points when f is given as an explicit Boolean
circuit lies within the complexity class TFNP [FGMS20, EPRY20]. This class contains
total search problems whose solutions are guaranteed to exist and are verifiable in
polynomial time, and is often viewed as the search analogue to NP ∩ coNP. As a result,
these problems are believed unlikely to be NP-hard unless NP = coNP. Despite this,
no polynomial-time algorithms are known for computing fixed points in either of the
above settings.
While fixed points of monotone functions and contractions have been studied extensively
in isolation, their intersection has received significantly less attention. Furthermore,
this intersection is particularly interesting due to its connection to two graph-based
problems that we will later introduce: the ARRIVAL problem [DGK+ 17] and Simple
Stochastic Games [Con92, Sha53]. While the decision versions of these problems are
known to lie in NP ∩ coNP [DGK+ 17, Con92], it is not known whether they lie in P.
Importantly, both can be reduced to the task of finding a fixed point of a monotone
contraction [Has25, Con92, EPRY20].
We make four main contributions. First, we investigate a recent approach for computing
fixed points of contraction maps and provide implementation details that make the
method computable, as the original formulation requires iteration over an infinite
domain. We also explore a sampling-based strategy toward improving efficiency,
establishing partial results and identifying conditions under which such an approach
may be effective. We introduce a decomposition method to break down contraction
fixed-point problems, and show that this method achieves worst-case sub-exponential
time complexity for a certain class of functions. Lastly, we study monotone contractions
and develop a heuristic that can be integrated into many of the existing fixed-point
algorithms to achieve significant speed-ups.

1
 Chapter 2

Background

2.1 Preliminaries
The two problems of interest are the ARRIVAL problem and Simple Stochastic Games.
We now give brief descriptions of both. More detailed definitions are provided in
Section 2.2.
The ARRIVAL problem is defined as follows. It describes a deterministic process on a
directed graph in which each vertex has an out-degree of two, with one of the outgoing
edges labelled as even and the other as odd. Starting from some designated vertex
o, a single token moves through the graph. At each vertex, the token is routed along
the least-used outgoing edge so far, and in case of a tie, it defaults to the even edge.
Assuming that two designated destination vertices t and t are reachable from every
vertex in the graph, the task is to decide which of these two destinations is visited first
by the token. See Figure A.1a for an illustration.
Similarly, (Condon’s) Simple Stochastic Games are defined as follows. We are given a
directed graph, in which each node is controlled by a MAX player, a MIN player or
is a “random node”. Additionally, there are two sink nodes labelled 0 and 1. Starting
from some vertex, the game goes as follows: at MIN nodes, the MIN player selects
the outgoing edge, at MAX nodes, the MAX player selects the outgoing edge, and at
“random nodes” the outgoing edge is selected according to fixed transition probabilities.
The MAX player aims to maximise the probability of reaching the 1 sink, while the
MIN player aims to minimise it. The problem is to determine, for each vertex, the
probability that the MAX player reaches the 1 sink under optimal play by both players.
See Figure A.1b for an illustration.
While simple to describe, these problems have no known polynomial-time solution.
However, both can be reduced to finding the fixed point of a function that is both
monotone and contractive. This connection motivates our focus on algorithms for
computing such fixed points.
To study these problems algorithmically, it is important to clarify the model of access to
the function that describes these fixed points. In the “white-box” model, the function
is given explicitly, for example, as a Boolean circuit, allowing an algorithm to inspect

2
Chapter 2. Background 3

its internal structure. In contrast, in the “black-box” model, the function can only be
accessed through queries. Given an input, the algorithm may evaluate the function,
but it has no access to its internal workings. This gives rise to two ways to measure
computational complexity. We can either measure the number of computational steps
or comparisons in the traditional sense, which is known as time complexity, or we can
measure the query complexity, which is the number of times we call the black-box
oracle f (x).
Throughout this dissertation, we consider both the black-box model and settings in
which the algorithm is given additional information about the function, specifically
whether the i-th component of the function “depends” on the j-th input coordinate for
each pair (i, j). Further details are provided in Section 3.1.2.
We now introduce the definitions and results from the monotone and contraction settings
that will be used throughout this dissertation.

2.1.1 Contraction Mapping Fixed Points

At a high level, contraction maps “pull” points closer together. If we keep applying the
map, the points move progressively closer, eventually reaching a point that no longer
changes under the mapping, i.e., a fixed point. The behaviour of such mappings depends
on how distance is measured. Different choices of metric can affect both the properties
of the function and the efficiency with which fixed points can be computed.
To formalise this, we work in metric spaces, which provide a notion of distance between
elements in a set.
Definition 2.1.1. Let X be a set. A metric on X is a function dX : X × X → R such that
for all x, y, z ∈ X, the following conditions are satisfied:
1. dX (x, y) ≥ 0;
2. dX (x, y) = 0 if and only if x = y;
3. dX (x, y) = dX (y, x);
4. dX (x, z) ≤ dX (x, y) + dX (y, z) .
The pair (X, dX ) is called a metric space.
Definition 2.1.2. A sequence (xn ) is Cauchy on a metric space (X, dX ) if for every
ε > 0, there exists N ∈ N such that dX (xn , xm ) < ε for all n, m > N.
Definition 2.1.3. A metric space (X, dX ) is said to be complete if every Cauchy sequence
in X converges to an element in X.
Definition 2.1.4. Let (X, dX ) be a metric space. A function f : X → X is called a
contraction mapping if there exists a constant 0 ≤ λ < 1 such that for all x, y ∈ X,

dX ( f (x), f (y)) ≤ λ dX (x, y). (2.1)

If we relax the constraint on λ to 0 ≤ λ ≤ 1 we call f non-expansive.

Chapter 2. Background 4

Definition 2.1.5. Let x, y ∈ Rd be vectors where x = (x1 , . . . , xd ) and y = (y1 , . . . , yd ).

For any 1 ≤ p ≤ ∞, the ℓ p distance between x and y is defined as:
1. For 1 ≤ p < ∞:
!1/p
d
dX (x, y) = ||x − y|| p = ∑ |xi − yi| p (2.2)
i=1

2. For p = ∞:
dX (x, y) = ||x − y||∞ = max |xi − yi | (2.3)
1≤i≤d

Throughout this dissertation, we primarily consider the metric space ([0, 1]d , d p ), where
d p (x, y) = ∥x − y∥ p . We refer to contractions under this metric as ℓ p -contractions.
Sometimes, we will assume the contraction coefficient λ > 0 for simplicity; if λ = 0
then f is constant function, and the fixed point is trivial. The reason contractions are
useful is the guarantee of not only existence, but also uniqueness.
Theorem 2.1.6 (Banach’s Fixed Point Theorem [Ban22]). Let (X, dX ) be a complete
metric space, and let f : X → X be a contraction with constant λ < 1. Then there exists
a unique fixed point x∗ ∈ X such that f (x∗ ) = x∗ . Moreover, for any starting point
x0 ∈ X, the iterative sequence xn+1 = f (xn ) converges to x∗ .
Definition 2.1.7. A subset C ⊆ Rd is said to be convex if for all x, y ∈ C and any t ∈ [0, 1],
the point (1 − t)x + ty belongs to C.
Remark 2.1.8. If f is a non-expansive mapping on a general complete metric space, a
fixed point is not guaranteed to exist. However, if the domain is additionally assumed
to be a compact1 and convex subset of Rd (or more generally a Banach space), a fixed
point exists (though it may not be unique) [Bro11].
Definition 2.1.9. Let (X, dX ) be a metric space and f : X → X be a map. For a given
ε > 0, a point x̃ ∈ X is called an ε-approximate fixed point of f if

dX ( f (x̃), x̃) ≤ ε. (2.4)

Lemma 2.1.10. Let (X, dX ) be a complete metric space and let f : X → X be a con-
traction mapping with coefficient λ ∈ (0, 1). Let x∗ ∈ X be the unique fixed point
of f . For any initial point  x0 ∈ X, the sequence defined by xn+1 = f (xn ) requires
log(1/ε)
at most n = O log(1/λ) iterations to guarantee an ε-approximate fixed point, i.e.,
dX (xn , f (xn )) ≤ ε.

Proof. By the definition of a contraction mapping, the distance between xn and xn+1 is
bounded by:
dX (xn , f (xn )) = dX (xn , xn+1 ) ≤ λn dX (x0 , x1 ) (2.5)
1A subset K ⊆ Rd is compact if and only if every sequence in K has a subsequence that converges to
a point in K.
Chapter 2. Background 5

To ensure dX (xn , f (xn )) ≤ ε, it is sufficient to require λn dX (x0 , x1 ) ≤ ε. Assuming

x1 ̸= x0 (otherwise x0 = x∗ and the fixed point is already found), we can rearrange this
to λn ≤ ε/dX (x0 , x1 ). Taking the natural logarithm:
 
ε
n ln(λ) ≤ ln (2.6)
dX (x0 , x1 )

Dividing both sides by ln(λ):

ln(ε) − ln(dX (x0 , x1 ))

n≥ (2.7)
ln(λ)

log ε−1
   
log ε
Treating the initial distance dX (x0 , x1 ) as a constant yields n = O log λ = O log λ−1
.

While Banach’s Fixed Point Theorem provides an iterative method to find the fixed
log(1/ε)
point, this method requires O log(1/λ) steps, which for λ ≈ 1 can become impractical.
As a result, more efficient approaches have been developed.
For Euclidean spaces (i.e. (X, ℓ2 ) with X ⊆ Rd ) more efficient algorithms are possible.
The Inscribed Ellipsoid algorithm by Sikorski et al. [STW93] can find an ε-approximate
fixed point with time complexity O (d log(ε−1 )).
Unfortunately, for other ℓ p -normed metric spaces, we do not have algorithms as efficient
as in the ℓ2 case. The first algorithm for finding fixed points of non-expansive maps
under norms other than ℓ2 was developed for the ℓ∞ norm over [0, 1]2 by Shellman
et al. [SS02], achieving a time complexity of O (log(ε−1 )). The method works by
maintaining a search space that is guaranteed to contain the fixed point. The algorithm
evaluates the function at the geometric centre2 of the current search space, using the
result to discard a significant volume where a fixed point cannot exist. Importantly, when
an area is discarded from the search space, the remaining space maintains geometric
properties that allow quick computation of geometric centres that are used to repeat the
procedure.
This approach was later generalised to [0, 1]d with an O (logd (ε−1 ))-time algorithm by
Shellman et al. [SS03]. The generalisation works by temporarily fixing one coordi-
nate, reducing the problem to d − 1 dimensions to compute a partial solution. It then
substitutes this partial solution back into the function to check the error of that fixed
coordinate. If the error exceeds the allowed tolerance, the algorithm adjusts the fixed
coordinate and repeats the process until the error across all dimensions is small enough.
For the ℓ1 norm, Fearnley et al. [FGMS20] gave an algorithm with the same asymptotic
running time of O (logd (ε−1 )) by following the same general idea as in the ℓ∞ case.
Recently, Chen et al. [CLY25] presented another algorithm for the ℓ∞ case that finds an
ε-approximate fixed point using O (d 2 log(ε−1 )) queries to the function f (x). The algo-
rithm specifically selects points x to query such that each query rules out a large fraction
2The average of all the coordinates in the search area.
Chapter 2. Background 6

of the search space where the fixed point could lie. Haslebacher et al. [HLSW25] ex-
tended this idea to all ℓ p -norms using at most O d 2 log ε−1 + log(1/(1 − λ)) queries

to f (x).
Both of these algorithms achieve polynomial query complexity in d, log(ε−1 ) and
1
log( 1−λ ), however, they require Ω(ε−d ) steps to determine the query points.
Addressing this bottleneck and building on the work of Shellman et al. [SS03] and
Fearnley et al. [FGMS20], Chen et al. [CLY26] established a O (log⌈d/2⌉ (1/ε))-time
algorithm for finding ε-approximate fixed points of non-expansive maps over [0, 1]d for
both the ℓ1 and ℓ∞ norms.
For the ℓ∞ -norm, they demonstrated how to break down a d-dimensional problem into
two lower a and b-dimensional problems. An outer algorithm searches for a solution for
the b-coordinates. Every time this outer algorithm needs to test a guess, it fixes those
b coordinates and runs a separate inner algorithm to find an approximate fixed point
for the remaining a-coordinates. Setting a = 2 and using the result from Shellman et
al. [SS02], they obtained the quadratic speed-up with respect to the result from Fearnley
et al. [FGMS20].
For the ℓ1 -norm, this a, b decomposition does not apply, and instead they singled out
two coordinates at a time; these coordinates are fixed, and just like in the ℓ∞ case, it
runs a recursive inner algorithm to find a solution for the remaining d − 2 coordinates.
However, these d − 2 coordinates are solved to a much smaller ε error; this means that
once this partial solution is found, if the overall error is too large, the algorithm knows
the bulk of the error must be coming from the two fixed coordinates. It turns out only
one of the two coordinates contributes the majority of the error, and pinpointing which
of the two it is allows them to discard half of the search space in that dimension.
The algorithms above summarise the state-of-the-art for the computation of fixed
points in ℓ p metric spaces. Next, we will introduce a complementary order-theoretic
perspective.

2.1.2 Monotone Functions and Lattice-Theoretic Fixed Points

Monotone functions preserve the order of their domain. Whenever x ≤ y, we have
f (x) ≤ f (y). If we consider a finite, ordered set with a least element x0 , and we
repeatedly apply xn+1 = f (xn ) where f is monotone, the resulting sequence is increasing.
Since the set is finite, this sequence must at some point stabilise. When this happens,
we obtain a point x such that f (x) = x, i.e. a fixed point.
The behaviour of these maps depends on the “structure” of this order. In some sets,
every element may be comparable against every other element (e.g. the natural numbers
with their usual order), in others, elements might be incomparable (e.g. in N2 with the
component-wise order, where x ≤ y if and only if xi ≤ yi for all components i)).
To better understand this, we introduce some definitions and results.
Definition 2.1.11. A partially ordered set (or poset) is a pair (P, ≤) consisting of a
set P and a binary relation ≤ on P that satisfies the following three properties for all
Chapter 2. Background 7

x, y, z ∈ P:
1. x ≤ x;
2. If x ≤ y and y ≤ x, then x = y;
3. If x ≤ y and y ≤ z, then x ≤ z.
Definition 2.1.12. Let (P, ≤) be a poset and S ⊆ P. An element u ∈ P is the least upper
bound of S if s ≤ u for all s ∈ S, and for any other upper bound v of S, u ≤ v. An
element l ∈ P is the greatest lower bound of S if l ≤ s for all s ∈ S, and for any other
lower bound m of S, m ≤ l.
Definition 2.1.13. A lattice is a poset (L, ≤) in which every pair of elements x, y ∈ L has
a least upper bound and a greatest lower bound. Furthermore, a poset where every subset
S ⊆ L has both a least upper bound and a greatest lower bound is called a complete
lattice.
A complete lattice has a greatest element (top, denoted ⊤) and a least element (bottom,
denoted ⊥) [Bir40].
The simplest kind of ordered structure is one where all elements can be compared to
one another.
Definition 2.1.14. A totally ordered set (or chain) is a poset (P, ≤) with the additional
property that every pair of elements is comparable. That is, for all x, y ∈ P, either x ≤ y
or y ≤ x.
Definition 2.1.15. A chain C ⊆ P in a poset P is a maximal chain if there exists no
element x ∈ P \C such that C ∪ {x} is also a chain.
Definition 2.1.16. Let (L, ≤L ) and (M, ≤M ) be posets. A function f : L → M is mono-
tone if for all x, y ∈ L, x ≤L y =⇒ f (x) ≤M f (y).
Working directly with abstract posets can be difficult. However, any poset can be
represented as a “product” of totally ordered sets via an “order embedding”.
Definition 2.1.17. Let (P1 , ≤1 ), . . . , (Pd , ≤d ) be partially ordered sets. The product of
these posets is the Cartesian product P1 × · · · × Pd with the coordinate-wise partial order
≤, defined such that for x = (x1 , . . . , xd ) and y = (y1 , . . . , yd ), we have x ≤ y if and only
if xi ≤i yi for all i ∈ {1, . . . , d}.
Definition 2.1.18. Let (L, ≤L ) and (M, ≤M ) be posets. A function f : L → M is an
order embedding if for all x, y ∈ L, x ≤L y ⇐⇒ f (x) ≤M f (y).
A monotone function preserves the order, but it can, for example, map two incomparable
elements to comparable ones or not be injective. An order embedding satisfies x ≤
y ⇐⇒ f (x) ≤ f (y), so the order on L is faithfully represented by its image in M.
Theorem 2.1.19 ( [DM41]). Any poset, and thus any lattice (L, ≤), can be order-
embedded into a Cartesian product of totally ordered sets.
Definition 2.1.20. The dimension of a poset L is defined as the smallest integer d such
that L can be order-embedded into the product of d totally ordered sets.
Chapter 2. Background 8

Theorem 2.1.21 (Tarski’s Fixed Point Theorem [Tar55]). Let (L, ≤) be a complete
lattice and let f : L → L be a monotone function. Then the set of fixed points of f ,
defined as Fix( f ) = {x ∈ L | f (x) = x}, is not empty.
Theorem 2.1.22. Let (L, ≤) be a finite complete lattice of dimension d and |L| elements,
and let f : L → L be a monotone function. The procedure defined by xn+1 = f (xn ) with
x0 = ⊥ converges to the least fixed point of f using O (d · |L|) queries to f in the worst
case.

Proof. We first show that the sequence converges to the least fixed point. The sequence
is monotone non-decreasing: x0 = ⊥ ≤ f (⊥) = x1 by definition of the least element.
Assuming xi−1 ≤ xi , the monotonicity of f implies f (xi−1 ) ≤ f (xi ), giving xi ≤ xi+1 .
By induction, the sequence forms a chain x0 ≤ x1 ≤ x2 ≤ . . . .
Because L is finite, this chain must eventually stabilise at some finite step k, meaning
xk = xk+1 = f (xk ). Thus, x∗ = xk is a fixed point. To see that it is the least fixed point,
let y ∈ L be any fixed point of f (so f (y) = y). We show by induction that xi ≤ y for all
i ≥ 0. For the base case, x0 = ⊥ ≤ y. For the inductive step, assume xi ≤ y. Applying
the monotone function f gives f (xi ) ≤ f (y), which is xi+1 ≤ y. Therefore, the limit
point satisfies x∗ ≤ y, proving x∗ is the least fixed point.
Now, bound the number of queries required to reach x∗ . Since L has dimension d, there
exists an order embedding φ : L → C1 × · · · ×Cd , where each C j is a totally ordered set.
For any x ∈ L, let φ(x) = (φ1 (x), . . . , φd (x)). Because L contains exactly |L| elements,
the image of L under the projection φ j onto any individual chain C j is made up of at
most |L| distinct elements.
If xi is not a fixed point, then xi < xi+1 . Order embeddings preserve strict inequalities,
φ(xi ) < φ(xi+1 ), meaning the “embedded representation” must strictly increase in at
least one of its d coordinates. Since each coordinate can strictly increase a maximum
of |L| − 1 times, the maximum possible number of increases across the entire product
poset is d · (|L| − 1).
Therefore, the sequence must stabilise after at most d · |L| steps. The algorithm evaluates
f at most d ·(|L|−1)+1 times, yielding a worst-case query complexity of O (d ·|L|).

In practice, we are often interested in finding approximate fixed points for functions
defined on a continuous domain (in particular f : [0, 1]d → [0, 1]d ). To use lattice-
theoretic algorithms we can discretise the domain into a finite grid. Under the coordinate-
wise partial order, this grid is a finite complete lattice. Choosing a grid resolution where
adjacent points are separated by a distance proportional to ε (depending on the metric),
we guarantee that the exact fixed point in the continuous domain is close enough to a
grid point that this grid point is an ε-approximate fixed point. Although f may map grid
points off the grid, this can be handled by composing f with a rounding mechanism back
to the grid, which preserves monotonicity and yields an equivalent discrete problem
[EPRY20]. However, this means the lattice contains Ω(ε−d ) points, so in the worst case
the procedure in Theorem 2.1.22 takes exponential time in the dimension.
Chapter 2. Background 9

It is worth noting that by Theorem 2.1.19, any finite lattice of dimension d can be order-
embedded into a product of d totally ordered sets. Therefore, developing algorithms for
this d-dimensional grid solves the problem for any lattice of dimension d, provided the
grid is fine enough.
Since the procedure from Theorem 2.1.22 becomes highly inefficient for our discretisa-
tion, more sophisticated algorithms have been developed.
Chang et al. [CLT08] established a binary-search-based log(|L|)-query algorithm for
fixed points on totally ordered lattices. By decomposing a general lattice into its set of
maximal chains M , this algorithm can be applied to each chain individually. However,
this yields a query complexity of O (|M | log(|L|)), which in practice is intractable.
Dang et al. [DQY24] studied the computational complexity of finding a fixed point
under two different models of access to the monotone function: the black-box model,
in the same sense as stated in Section 2.1.1, and the “polynomial function model”,
where f is given by an explicit algorithm running in time polynomial in log |L|. For the
black-box model, they created an algorithm with a time complexity of O (logd (|L|)).
In the polynomial function model, they introduced a polynomial-time algorithm and
demonstrated that verifying the uniqueness of a fixed point is coNP-hard, even in one
dimension.
In two dimensions, the approach of Dang et al. [DQY24] has a time complexity of
O (log2 (|L|)), which Etessami et al. [EPRY20] later proved to be optimal by obtaining
a corresponding Ω(log2 (|L|)) lower bound
For dimensions greater than two, Fearnley et al. [FPS22] introduced an algorithm that
requires O (log2⌈d/3⌉ (|L|)) queries. Recently, Chen et al. [CL22] improved this bound
to O (log⌈(d+1)/2⌉ (|L|)) queries.
We now introduce two problems that can be solved by finding a fixed point of a
monotone and ℓ p -contracting function.

2.2 Motivation
The ARRIVAL problem and Simple Stochastic Games (SSGs) are particularly relevant
to the class of functions we study. Both (in their decision versions) have been shown
to lie within the complexity class NP ∩ coNP [DGK+ 17, Con92], but neither is known
to have a polynomial-time algorithm. As we will see, recent results show that solving
either problem reduces to finding approximate fixed points of monotone, contracting
functions (in the ℓ1 and ℓ∞ norms, respectively).

2.2.1 The ARRIVAL Problem

The ARRIVAL problem was originally introduced by Dohrau et al. [DGK+ 17] as a
decision problem. Here we consider its total search variant, as formulated by Karthik et
al. [Kar17].
The problem is defined on a switch graph G = (V, E, s0 , s1 ), where each non-terminal
Chapter 2. Background 10

vertex v ∈ V has exactly two outgoing edges, given by “successor functions” s0 (v)
(even) and s1 (v) (odd). The graph has two terminal nodes, t and t, and an initial vertex
o ∈ V.
Let Xt denote the set of vertices from which t is not reachable from o. To remove these
vertices, we overwrite the successor functions so that every v ∈ Xt has both successors
equal to t.
Doing this guarantees that every vertex in G has a directed path to at least one terminal
node, ensuring that the token eventually reaches either t or t¯.
Gärtner et al. [GHH21] developed a sub-exponential-time algorithm for n-vertex graphs
and obtained polynomial-time algorithms for graphs that are close to being acyclic,
specifically graphs that can be made acyclic by removing a constant number of vertices.
They also showed that the problem can be framed as the fixed point of a monotone
function on a finite lattice.
More recently, Haslebacher [Has25] reduced a generalisation of ARRIVAL, called G-
ARRIVAL, to the problem of computing an approximate fixed point of an ℓ1 -contraction
map on [0, 1]n .
We now formalise the fixed point map used in this latter reduction, simplified for the
standard ARRIVAL problem.
Let G = (V, E, s0 , s1 ) be a directed graph with |V | = n, where for each vertex v ∈ V ,
the functions s0 (v) and s1 (v) denote its even and odd successors, respectively. Let
T = {o,t,t} ⊂ V be the set of terminal vertices, and let k = |V \ T | = n − 3 be the
number of non-terminal vertices.
We represent the state of the process by assigning a (scaled) amount of token mass to
each non-terminal vertex.
Let x ∈ [0, 1]k . To ensure we work within [0, 1]k , we divide the mass at each vertex v by
M = 2n . While this is to make the map stay within [0, 1]k , it is unnecessary for the logic
of the problem. The actual token mass is simply M · xv .
To simplify the notation when routing mass through the graph, we extend x to x′ ∈ Rn≥0
defined over all vertices:

1
 if u = o,
′
xu = 0 if u ∈ {t,t}, (2.8)

M · xu if u ∈ V \ T.


This reflects that one unit of mass is created at the origin, and the terminals are sinks.
We now model how mass is split at each vertex. Recall that in the ARRIVAL process,
edges are used in an alternating fashion, sending nearly half of the visits along each
outgoing edge.
To capture this, we define two functions that split an incoming quantity y ≥ 0 into two
nearly equal parts:
 j y k l y m j y k l y m
h0 (y) = min y − , , h1 (y) = max , y− . (2.9)
2 2 2 2
Chapter 2. Background 11

Here, h0 (y) is the token mass routed along the even edge, and h1 (y) to the odd edge.
We now define a map that takes in the current token mass at each vertex and updates it.
About half of the mass goes to the even edge and the remaining to the odd. For each
vertex v, we sum the contributions arriving from predecessors, separating between when
v is reached via an even or odd edge.
Definition 2.2.1. For each non-terminal vertex v ∈ V \ T , define ĝ(λ) : [0, 1]k → [0, 1]k
by !!
1
ĝ(λ) (x)v = min M, λ
M ∑ h0(xu′ ) + ∑ h1(xu′ ) . (2.10)
u:s (u)=v
0 u:s (u)=v
1

This expression puts together all mass sent to v, applies a “discount factor”, and rescales
the result back into the interval [0, 1].
To make sure ĝ(λ) (x) is a contraction, we need to discount it by multiplying it by λ.
 
1
λ ∈ 1− ,1 , (2.11)
1 + n2n

which slightly reduces the total mass at each step.

As shown by Haslebacher [Has25], this map is monotone with respect to the coordinate-
wise partial order and is ℓ1 -contracting with parameter λ.
1−λ
Furthermore, computing an approximate fixed point with an ℓ1 -error of ε < 4 allows
us to recover the solution.
The intuition is as follows. Without the discount factor, the fixed point would correspond
to a conserved flow of token mass through the graph, with exactly one unit reaching the
correct terminal. The discount ensures ĝ(λ) is an ℓ1 -contraction, but causes a small loss
of mass. With the choice of λ in Equation 2.11, the total loss at the terminals is less
than 1, so the terminal receiving more mass at the fixed point is the one reached by the
original process.

2.2.2 Simple Stochastic Games

Simple Stochastic Games (SSGs) are another setting in which fixed-point computations
become useful. They are two-player, zero-sum games played on directed graphs, where
the evolution of the game depends both on the players’ choices and on probabilistic
transitions. Earlier, we described Condon’s Simple Stochastic Games. Now, we give a
more detailed description, introduce Shapley’s Stochastic Games and explain how the
solution to these corresponds to finding a fixed point.
Condon’s SSGs are played on a directed graph G = (V, E) with n = |V | vertices. The
vertices are partitioned into three types: MAX vertices (V1 ), MIN vertices (V2 ), and
random vertices (V0 ). There are also two terminal vertices (a 0-sink and a 1-sink). We
assume at least one terminal is reachable from every vertex. Starting from an initial
vertex, the players move a token through the graph: MAX and MIN choose outgoing
edges at their respective vertices. Transitions at random vertices vi ∈ V0 to neighbours
Chapter 2. Background 12

v j ∈ V follow fixed probabilities pvi ,v j . The objective is to maximise (for MAX) or

minimise (for MIN) the probability of reaching the 1-sink. The problem for us is to
determine, for each vertex, given both players play optimally, this probability. This is
known as the value of the game at that vertex.
A realistic example of an SSG involves a criminal attempting to flee a city by car. As
the criminal navigates the streets, they go through various junctions. The police control
some junctions. Rather than attempting an immediate and potentially dangerous arrest,
the officers instead set up barricades to influence the route, steering the criminal’s
choices. At other intersections, heavy traffic disrupts the criminal’s decision-making,
forcing them to take a direction more or less at random. There are also quieter roads
where there is no traffic or police, and the criminal is free to choose any path.
In this scenario, the police are looking to channel the criminal towards a dead end, where
they can safely arrest them. The criminal, on the other hand, is trying to navigate through
the network to safely leave the city. The question becomes: given both the criminal and
the police work optimally, what is the probability of escape at each intersection?
This example will become relevant in Section [Link].
Shapley’s stochastic games generalise this model. Each vertex vi has action sets for both
players, and each pair of actions determines both an immediate reward and a probability
distribution over successor vertices. Additionally, there is a positive probability that
the game halts at each step, ensuring termination with probability one. The value of a
vertex is defined as the expected total reward under optimal play.
We now connect these games to fixed-point computation. Etessami et al. [EPRY20]
showed that the value of both Condon’s and Shapley’s SSGs can be characterised as
fixed points of monotone, ℓ∞ -contracting maps.
We first make this explicit for Condon’s SSGs. To compute the exact value vector
q∗ ∈ [0, 1]n , we can find the unique fixed point of a map f : [0, 1]n → [0, 1]n .
Definition 2.2.2. For each vertex vi , define


 λ ∑v j ∈V pvi ,v j x j if vi ∈ V0 (Random)

λ max(vi ,v j )∈E x j if vi ∈ V1 (Max)



f (x)i = λ min(vi ,v j )∈E x j if vi ∈ V2 (Min) (2.12)

0 if vi is the 0-sink





if vi is the 1-sink

λ

where λ = 1 − 2−poly(|G|) is a discount factor, and |G| denotes the total bit encoding size
of the game. Because the exact values of the undiscounted game are rational numbers
with denominators bounded by 2poly(|G|) , choosing λ to be this close to 1 ensures that
the map is a contraction while keeping the resulting fixed point close enough to the true
values to allow recovering the exact solution.
The operations defining f (minimum, maximum, and weighted averages) are all order-
preserving and non-expansive in the ℓ∞ norm. The factor λ ensures that f is a contrac-
tion.
Chapter 2. Background 13

A similar formulation applies to Shapley’s stochastic games.

Definition 2.2.3. For a Shapley SSG with n vertices, the value vector r∗ ∈ Rn is the
unique fixed point of a map f : Rn → Rn . For each vertex vi ∈ V , suppose the MAX and
MIN players have mi and ni available actions, respectively. Let j ∈ {1, . . . , mi } denote a
specific action chosen by the MAX player, and k ∈ {1, . . . , ni } denote a specific action
chosen by the MIN player.
Let Aij,k denote the immediate reward and Pij,k (r) denote the probability of transitioning
to state vr when MAX chooses action j, and MIN chooses action k. We define Bi (x) as
an mi × ni matrix with entries
n
Bi (x) j,k = Aij,k + ∑ Pij,k (r)xr (2.13)
r=1

Then the value of the state is given by:

f (x)i = max min yT Bi (x)z (2.14)

y∈Pmi z∈Pni

where Pk represents the set of all valid probability distributions over k actions (i.e., if
y ∈ Pk , then y = (y1 , . . . , yk ) such that yi ≥ 0 for all i ∈ {1, . . . , k} and ∑ki=1 yi = 1).
Intuitively, Equation 2.14 means the MAX player chooses a probability distribution over
their actions to maximise the expected reward, assuming the MIN player will counter
with a distribution chosen to minimise that same reward.
Monotonicity follows from the non-negativity of the transition probabilities, while
contraction arises because due to the halting probability, the total transition probability
is strictly less than 1.
In both ARRIVAL and SSGs, the problem reduces to finding a fixed point of a function
that is both monotone and contracting in some ℓ p -norm. This also raises the question:
can these properties together be exploited algorithmically?
Batziou et al. [BFG+ 25] address this question for monotone ℓ∞ -contractions. They
give an algorithm that computes an ε-approximate fixed point in dimension d using
O ((c · log(1/ε))⌈d/3⌉ ) queries. However, to date, this intersection has not been studied
further.
 Chapter 3

Theoretical Advancements

This chapter investigates algorithmic improvements for finding fixed points of ℓ p -

contractions. We first investigate the query-efficient ℓ p -contraction algorithm of Hasle-
bacher et al. [HLSW25]. In particular, developing implementation details and strategies
for turning the algorithm time efficient via sampling, a direction suggested but not inves-
tigated in their work. While a complete solution is beyond the scope of this dissertation,
we establish a useful stepping stone. If one can uniformly sample from the search space
where the fixed point could lie, and efficiently compute a query point from the sample
that rules out a large, guaranteed fraction of the remaining space upon evaluation, then
an efficient algorithm would follow.
Following this, we will introduce a decomposition algorithm that works particularly
well for a specific subclass of problems by exploiting variable dependencies.
Finally, we will investigate monotone ℓ p -contractions and demonstrate a heuristic that
can discard large portions of the search space under certain conditions and, in turn,
speed up existing algorithms.

3.1 Algorithms for ℓ p -contractions

3.1.1 Haslebacher’s Algorithm
One of the recent algorithms for finding fixed points of ℓ p -contractions is the geometric
framework introduced by Haslebacher et al. [HLSW25]. They generalise an earlier al-
gorithm designed for ℓ∞ -contractions [CLY25] by extending well-established geometric
concepts such as halfspaces and centerpoints from the ℓ2 -norm to the ℓ p -norm (formal
definitions are provided later in this section).
With these tools, they showed that it is possible to find an approximate fixed point, for
any ℓ p -norm, using only a polynomial number of calls to f . However, this guarantee
comes with a caveat: the existence of these query points is non-constructive (via
Brouwer’s Fixed Point Theorem [Bro11]) and computing them requires exponential
time.

14
Chapter 3. Theoretical Advancements 15

Their algorithm is based on the following: the ℓ p -contraction property guarantees that
f (x) will lie closer to the fixed point x∗ than x (under the ℓ p -norm):
|| f (x) − f (y)|| p ≤ λ||x − y|| p =⇒ || f (x) − x∗ || p = || f (x) − f (x∗ )|| p (3.1)
≤ λ||x − x∗ || p < ||x − x∗ || p (3.2)
As a result, we initialise the search space to S0 = [0, 1]d , and every time we query a
point c ∈ [0, 1]d , we can discard from the search space any point closer to c than to f (c).
Haslebacher et al. [HLSW25] show that there always exists a point that guarantees a
1
d+1 fraction of the search space will be discarded with each query.
The challenge is how to compute such “optimal” query points efficiently. Chen et
al. [CLY25] provide a constructive method in the ℓ∞ norm, while Haslebacher et
al. [HLSW25] leave the general ℓ p -norm case open. To develop the practical imple-
mentation of finding these points, we must establish some important definitions and
results.
In Euclidean geometry (i.e. under the ℓ2 -norm), halfspaces are typically defined using
the inner product. Given a point x ∈ Rd and a direction v ∈ Sd−1 , the corresponding
halfspace consists of all points z ∈ Rd satisfying ⟨v, z − x⟩ ≥ 0.
For general ℓ p -norms, however, inner product is not defined, and therefore this charac-
terisation does not directly generalise. To overcome this, Haslebacher et al. [HLSW25]
instead express the same condition in terms of distances. To bypass the need for an
inner product, one can show that in the Euclidean setting,
⟨v, z − x⟩ ≥ 0 ⇐⇒ ∥x − z∥2 ≤ ∥x − εv − z∥2 for all ε > 0. (3.3)
This reformulation no longer depends on the inner product and can therefore be extended
to arbitrary ℓ p -norms. Moreover, it can be interpreted as the limit of the halfspace
obtained by comparing x and x − εv as ε → 0. This motivates the definition of bisector
and limit ℓ p -halfspaces.
Definition 3.1.1 ( [HLSW25]). For fixed p ∈ [1, ∞] and distinct points x, y ∈ Rd , the
p
bisector ℓ p -halfspace Hx,y ⊆ Rd is defined as
p
Hx,y := {z ∈ Rd | ∥x − z∥ p ≤ ∥y − z∥ p }. (3.4)

Definition 3.1.2 ( [HLSW25]). For fixed p ∈ [1, ∞], point x ∈ Rd , and direction v ∈ Rd
p
with v ̸= 0, the limit ℓ p -halfspace Hx,v ⊆ Rd through x in the direction of v is defined as

p
Hx,v := {z ∈ Rd | ∀ε > 0 : ∥x − z∥ p ≤ ∥x − εv − z∥ p }. (3.5)

Observation 3.1.3 ( [HLSW25]). Let p ∈ [1, ∞], x ∈ Rd , and v ∈ Sd−1 be a direction

p
on the unit ℓ2 -sphere, where Sd−1 = {z ∈ Rd : ∥z∥2 = 1}. The limit ℓ p -halfspace Hx,v is
p
the intersection of all bisector ℓ p -halfspaces Hx,x−εv for all ε > 0:
p
\
p
Hx,v = Hx,x−εv . (3.6)
ε>0
p p
Consequently, Hx,v ⊆ Hx,x−εv holds for every ε > 0.
Chapter 3. Theoretical Advancements 16

(a) Bisector ℓ∞ halfspace (b) Limit ℓ∞ halfspace

Figure 3.1: Comparison of ℓ p -halfspaces. (a) The bisector ℓ p -halfspace Hx,y (shaded)
consisting of all points at least as close to x as to y. (b) The limit ℓ p -halfspace originating
at x along the direction v = y − x. The shaded wedge behind x represents the opposite
limit ℓ p -halfspace Hx,−v . Comparing the panels shows that Hx,−v ⊆ Hx,y .

While it is quite easy to verify if a point is contained in a bisector ℓ p -halfspace, it is not

as intuitive to do so for a limit ℓ p -halfspace. Haslebacher provides a useful lemma that
utilises subdifferentials for us to check containment in a limit ℓ p -halfspace.
Definition 3.1.4. Let h : Rd → R be a convex function and x ∈ Rd . A vector u ∈ Rd is
called a subgradient of h at x if and only if

h(x′ ) − h(x) ≥ ⟨u, x′ − x⟩ (3.7)

for all x′ ∈ Rd . The set of all such subgradients is called the subdifferential of h at x,
and is denoted by ∂h(x) ⊆ Rd .
Lemma 3.1.5 ( [HLSW25]). For any p ∈ [1, ∞], a point z ∈ Rd is contained in an
p
ℓ p -halfspace Hx,v if and only if there exists a subgradient u ∈ ∂∥z − x∥ p with ⟨u, v⟩ ≥ 0.
In Euclidean geometry, given a set of points P ∈ Rd there is always some point c ∈ Rd
1
such that every (ℓ2 ) halfspace that contains c, contains at least a d+1 fraction of the
points [Rad38]. Haslebacher’s main result was extending this to ℓ p -normed spaces.
Definition 3.1.6 ( [HLSW25]). For a probability measure µ on Rd with bounded support
and an arbitrary p ∈ [1, ∞], an ℓ p -centerpoint is a point c ∈ Rd such that

p 1 1
µ(Hc,v )≥ µ(Rd ) = (3.8)
d +1 d +1
for all v ∈ Sd−1 .
Definition 3.1.7 ( [HLSW25]). For a finite set of points P ⊆ Rd and an arbitrary
p ∈ [1, ∞], a discrete ℓ p -centerpoint is a point c such that

p |P|
|Hc,v ∩ P| ≥ (3.9)
d +1
for all v ∈ Sd−1 .
 Chapter 3. Theoretical Advancements 17

Theorem 3.1.8 ( [HLSW25]). For any p ∈ [1, ∞] and any probability measure µ on Rd
with bounded support, there always exists an ℓ p -centerpoint c ∈ Rd such that
p 1
µ(Hc,v )≥ µ(Rd ) (3.10)
d +1
for all directions v ∈ Sd−1 .
Furthermore, for any finite set of points P ⊂ Rd , there exists a discrete ℓ p -centerpoint
c ∈ Rd satisfying
p |P|
|Hc,v ∩ P| ≥ (3.11)
d +1
for all v ∈ Sd−1 .

Algorithm 1: Algorithm for Fixed Points of ℓ p -contractions [HLSW25].

Input: Contraction map f : [0, 1]d → [0, 1]d , tolerance ε, ℓ p -norm
Output: Approximate fixed point
1 Initialise search space S0 ← [0, 1]d ;
2 k ← 0;
3 // Compute the ℓ p -centerpoint with respect to the uniform volume measure of the
current search space Sk ;
4 x0 ← Centerpoint(S0 );
5 while ∥ f (xk ) − xk ∥ p > ε do
6 Sk+1 ← Sk \ Hxpk , f (xk ) ;
7 k ← k + 1;
8 xk ← Centerpoint(Sk );
9 end
10 return xk ;

(a) Iteration 1 (b) Iteration 2

Figure 3.2: A query at c creates a bisector ℓ∞ -halfspace Hc,∞f (c) (shaded) that cannot
contain the fixed point x∗ . The search space is then reduced (right), and a new ℓ p -
centerpoint is selected.

Haslebacher et al. prove that a region around the true fixed point x∗ is never discarded
during this process. Because we are guaranteed to never discard the solution, and every
Chapter 3. Theoretical Advancements 18

1
query reduces the search space volume by a fraction of d+1 , the search space shrinks
fast. A full derivation of the complexity bound can be found in their
 original work, their

2 −1
result is that finding an ε-approximate fixed point requires O d log ε + log 1−λ 1

queries to f .

[Link] Implementation on a Discrete Grid

Having established the theoretical foundations of Haslebacher et al. [HLSW25], we

now discuss how to compute these ℓ p -centerpoints.
In practice, finding fixed points requires discretising the continuous space with a grid
Gdb ⊆ [0, 1]d . For an integer b ≥ 1, this grid consists of the set of points x ∈ [0, 1]d with
rational coordinates x1 . . . , xd of the form xi = 2kbi for integers k1 , . . . , kd ∈ {0, 1, . . . , 2b }.
By choosing a sufficiently small grid resolution, in particular b ≥ log2 (d + λd)ε−1 ,

we guarantee there is a point x ∈ Gdb that is an ε-approximate fixed point of f [HLSW25].

The challenge arises because there is no guarantee that the ℓ p -centerpoint will lie in
the grid Gdb . While the problem of solving this issue for general ℓ p is open, rounding
mechanisms have been developed for the ℓ∞ and ℓ1 norms [CLY25, HLSW25]. These
provide a way to transform a given ℓ p -centerpoint c ̸∈ Gdb to a point c′ ∈ Gdb that is a
discrete ℓ p -centerpoint of the discretised search space Gdb ∩ Sk . However, this requires
an even finer grid resolution. For the ℓ1 and ℓ∞ norm we require :
   
2d 1 + λ 1
b ≥ log2 and b ≥ 4 + log2 (3.12)
ε 1−λ ε(1 − λ)
respectively. The lack of a rounding mechanism for other ℓ p -norms means we currently
cannot implement a query-efficient algorithm for general ℓ p -norms.
To compute a centerpoint, we must iterate through each point c in the grid Gdb ∩ Sk ,
where Sk is the search space after k − 1 queries to f (as in Algorithm 1). For each
point c, we must check the discrete centerpoint condition: for every possible direction
p 1
v ∈ Sd−1 , the limit ℓ p -halfspace Hc,v must contain at least a d+1 fraction of the total
points currently in the search space.
By Lemma 3.1.5, checking if a single point z is in a limit ℓ p -halfspace is trivial. The
bottleneck is verifying this condition across every direction v ∈ Sd−1 , which is an
uncountably infinite set.
As we will now prove, for the specific cases of p ∈ {1, ∞}, checking this condition over
an infinite continuous sphere can be reduced to checking a finite, discrete set of test
directions. Chen et al. [CLY25] established this for the ℓ∞ -norm; our work yields the
same test directions through a different derivation and provides an analogous, novel
result for the ℓ1 -norm.
Let P denote the finite set of points in the current search space. For any point c ∈ [0, 1]d
p
and any point z ∈ P, let y = z − c. By Lemma 3.1.5, the membership of z in Hc,v depends
entirely on the sign of the function φy (v):
p
z ∈ Hc,v ⇐⇒ φy (v) = max ⟨u, v⟩ ≥ 0. (3.13)
u∈∂∥y∥ p
Chapter 3. Theoretical Advancements 19

n o
Proposition 3.1.9. Let V∞ = √1 (σ1 , . . . , σd ) σi ∈ {−1, 1} . Then c satisfies Equa-
d
tion 3.13 for p = ∞ for all v ∈ Sd−1 if and only if it satisfies it for all v ∈ V∞ .

Proof. Let z ∈ P and y = z − c. Let I = {i : |yi | = ∥y∥∞ }. The subdifferential of the

ℓ∞ -norm at y is ∂∥y∥∞ = conv({sgn(yi )ei : i ∈ I}).
Because the inner product is linear, its maximum over the convex set ∂∥y∥∞ is achieved
at one of the set’s extreme points [Bau58], giving φy (v) = maxi∈I sgn(yi ) vi .
Thus, z ∈ Hc,v
∞ if and only if max
i∈I sgn(yi )vi ≥ 0. This value depends only on the
signs of the coordinates vi , and can change sign only when some vi = 0, i.e., on the
ℓ2 -hyperplanes vi = 0.
These d ℓ2 -hyperplanes partition Rd (and Sd−1 ) into 2d orthants. Within any orthant,
(sgn(v1 ), . . . , sgn(vd )) is fixed, so φy (v) cannot change sign. Hence, the membership of
z in Hc,v
∞ is constant within each orthant.

Since this holds for every z ∈ P, the quantity |Hc,v

∞ ∩ P| is also constant on each orthant.

Therefore, testing one direction per orthant, v ∈ V∞ , suffices.

Proposition 3.1.10. Consider the set of ℓ2 -hyperplanes defined by vi = 0 for i = 1, . . . , d

and ∑di=1 σi vi = 0 for σi ∈ {−1, 1}. These ℓ2 -hyperplanes split Sd−1 into finitely many
regions. Let V1 ⊂ Sd−1 contain exactly one point from each region. Then c satisfies
Equation 3.13 for p = 1 for all v ∈ Sd−1 if and only if it satisfies it for all v ∈ V1 .

Proof. Fix z ∈ P and define y = z − c. For each coordinate i, the subdifferential of ∥y∥1
is (
{sgn(yi )}, yi ̸= 0,
[∂∥y∥1 ]i = (3.14)
[−1, 1], yi = 0.

Hence, Equation 3.13 with p = 1 can be written as:

φy (v) = ∑ sgn(yi )vi + max ∑ vi ui = ∑ sgn(yi )vi + ∑ |vi | ≥ 0. (3.15)

i:yi ̸=0 u∈[0,1] i:y =0 i:yi ̸=0 i:yi =0
i

Fix τ ∈ {−1, 1}d . In the region Rτ := {v ∈ Sd−1 : sgn(vi ) ∈ {τi , 0} for all i}, we have
|vi | = τi vi . Therefore, on Rτ , φy (v) is a linear function φy (v) = ∑di=1 αi vi , where αi =
sgn(yi ) if yi ̸= 0, and αi = τi if yi = 0. In all cases, αi ∈ {−1, 1}.
So, within any region Rτ , the function φy is a linear function of v.
The only places where this expression can change are when v crosses an ℓ2 -hyperplane
of the form vi = 0. Once the linear function v 7→ ∑di=1 αi vi is fixed, it can only change
from positive to negative (or vice versa) when it crosses the set where it equals zero.
Therefore, the sign of φy (v) depends on which side of the following set of ℓ2 -hyperplanes
v lies on:
d
vi = 0 (i = 1, . . . , d), ∑ σivi = 0 (σi ∈ {−1, 1}). (3.16)
i=1
Chapter 3. Theoretical Advancements 20

These ℓ2 -hyperplanes split Sd−1 into finitely many regions. On any single region F,
none of the expressions vi changes sign, and none of the expressions ∑ σi vi changes
sign. As a result, φy (v) keeps the same sign (or is always zero) throughout F.
1 } is constant on each
Since this is true for every z ∈ P, the indicator function 1{z ∈ Hc,v
region F. Hence, the quantity |Hc,v 1 ∩ P| is also constant on each region.

Let V1 contain one point from each region. Since the centerpoint condition depends
1 ∩ P|, it is enough to check the condition only for each v ∈ V .
only on |Hc,v 1

(a) ℓ∞ sphere partition (b) ℓ1 sphere partition

Figure 3.3: Partition of the unit sphere S2 into regions of “test directions” in R3 . (a) For
the ℓ∞ norm, the planes vi = 0 split the sphere into 8 orthants, requiring one direction
per region. (b) For the ℓ1 norm, additional planes ∑ ±vi = 0 create a finer partition.

With these results, we can write an algorithm that is guaranteed to work for p ∈ {1, ∞}.
However, in practice, this is limited by the required grid resolution in Equation 3.12. To
illustrate this, consider solving an SSG with just five nodes (d = 5) with a contraction
coefficient λ ≈ 1 − 215 to a precision of ε = 1 × 10−4 . We would need a grid with
roughly 2 × 104 points per dimension. This means the full search space will contain
approximately 1027 points. Running the ℓ p -centerpoint procedure over this space would
require approximately 25 (1027 )2 ≈ 1055 operations, making the algorithm completely
intractable.
For this reason, we now consider an alternative approach that replaces the need for
having such a thin grid with sampling.

[Link] A Sampling-Based Algorithm for Approximate Centerpoints

A natural alternative to computing an ℓ p -centerpoint is using sampling, something

suggested as a possible area of future research in Haslebacher et al. [HLSW25].
The approach is as follows: In each step of the algorithm, we replace the true search
space with a uniform random sample drawn from it, then we compute an ℓ p -centerpoint
of this sample, and use it as a proxy for a ℓ p -centerpoint of the search space. The main
issue is whether this substitution is reliable: does a ℓ p -centerpoint computed from a
sample discard a large amount of the search space with high probability?
Chapter 3. Theoretical Advancements 21

In this section, we will investigate that idea by addressing the following three questions:
1. How do we define and compute an approximate ℓ p -centerpoint?
2. How many samples are needed to obtain a good ℓ p -centerpoint approximation?
3. Can we sample uniformly from the search space, and if so, how?

[Link].1 How do we define and compute an approximate centerpoint? Comput-

ing an exact ℓ p -centerpoint is a difficult problem. Even in the ℓ2 case [Ten92], which
motivates working with a relaxed notion.
Recall that a centerpoint c requires that every limit ℓ p -halfspace Hc,v contains at least a
1
d+1 fraction of the total mass of the search space. An α-ℓ p -centerpoint weakens this
guarantee to a fraction α.
Definition 3.1.11. Let µ be a probability measure with bounded support. A point
c ∈ Rd is an α-ℓ p -centerpoint under the ℓ p -norm if, for every direction v ∈ Sd−1 , the
ℓ p -halfspace through c in direction v satisfies:
p


µ Hc,v ≥ α. (3.17)

1
In the ℓ2 norm, there are algorithms that can find a 2d -ℓ2 -centerpoint, with running time
polynomial in the dimension d, the number of samples n and log(q−1 ), where q denotes
the failure probability of the algorithm [Che24]. However, no known algorithm extends
this result to other ℓ p -norms.
Given this constraint, we continue under the assumption that we have access to a black-
box algorithm that returns an α-ℓ p -centerpoint for a given set of n points in d dimensions.
This allows us to focus on the statistical aspect of the problem, namely, how many
samples we would need to make sure, with high probability, that the ℓ p -centerpoint of
the sample is an approximate ℓ p -centerpoint of the measure.

[Link].2 How many samples are needed to obtain a good approximation? We

now investigate how well a finite sample represents the full search space. This heavily
depends on how “complex” the search space is. We will show that, in our case, the
search space is simple enough that, with a small number of samples, for every subset of
the search space, the fraction of sample points that fall into that subset is close to the
true probability mass of that region.
To formalise this, we consider a family of subsets of the space (called ranges), which in
our case will be the family of all bisector ℓ p -halfspaces in Rd . These are the regions
used to discard parts of the search space in each iteration. We require that the sample
approximates the measure of every such set up to a small additive error ε.
Definition 3.1.12. A range space is a pair (X, R ) where X is a ground set and R is a
collection of subsets of X called ranges. A finite subset S̃ ⊆ X is an ε-approximation
for (X, R ) with respect to a probability measure µ on X if, for every R ∈ R :
|S̃ ∩ R|
− µ(R) ≤ ε. (3.18)
|S̃ |
Chapter 3. Theoretical Advancements 22

The number of samples we need to achieve an ε-approximation depends on how complex

the family of sets R is. The complexity of a family of sets is measured by Vapnik–
Chervonenkis (VC) dimension. Specifically, the VC dimension measures the maximum
number of points in X such that the family of sets can shatter them.
Definition 3.1.13. Let (X, R ) be a range space. A subset C ⊆ X is said to be shattered
by R if every possible subset of C can be formed by intersecting C with a set from R ;
formally, if
{C ∩ R | R ∈ R } = P (C), (3.19)
where P (C) denotes the powerset of C.
The VC dimension of R , denoted VCdim(R ), is the supremum of the sizes of all finite
sets C ⊆ X that are shattered by R [VC15]. If R can shatter sets of arbitrarily large
size, then VCdim(R ) = ∞.
An important result from statistical learning theory states that if the VC dimension
of the range space is finite, then a small random sample is sufficient to ensure an
ε-approximation with high probability.
Theorem 3.1.14 ( [Tal94]). Let (X, R ) be a range space with finite VC dimension
V = VCdim(R ). For any q, ε ∈ (0, 1), a random sample S̃ of size
V + log(q−1 )
 
N=O (3.20)
ε2
drawn i.i.d. according to µ is an ε-approximation for (X, R ) with probability at least
1 − p.
To be able to apply this theorem, we need a bound on the VC dimension of the family
of all bisector ℓ p -halfspaces. These are the regions used to define the search space, and
so they are the sets whose measure we need to approximate.
Using a general result on families defined by polynomial inequalities, one can show the
VC dimension of bisector ℓ p -halfspaces is bounded.
Theorem 3.1.15 ( [GJ95]). Let R be a set of ranges over a domain X. If membership
of x ∈ X can be decided by a Boolean formula over s polynomial inequalities of degree
at most deg parametrised by k real variables, then the VC dimension of R is bounded
by V ≤ 2k log(8e · deg · s).
Proposition 3.1.16. For p ∈ Z+ ∪ {∞}, the VC dimension of bisector ℓ p -halfspaces in
Rd is O (d log d).

Proof. We use Theorem 3.1.15. Since any bisector ℓ p -halfspace is parametrised by two
points u, v ∈ Rd , we have a total of k = 2d real parameters.
Case 1: The ℓ∞ norm.
The equation ∥z − u∥∞ ≤ ∥z − v∥∞ can be written as:
!
d d
(zi − ui )2 − (z j − v j )2 ≤ 0
_ ^

Φ∞ (z, u, v) ⇐⇒ . (3.21)
j=1 i=1
Chapter 3. Theoretical Advancements 23

This involves s = d 2 polynomial predicates, each of degree deg = 2. Substituting into

the bound gives V ≤ 4d log(16e · d 2 ) = O(d log d).
Case 2: The ℓ p -norm for integer p.
The condition ∥z − u∥ p ≤ ∥z − v∥ p can be written as:

d d
∑ |zi − ui| p ≤ ∑ |zi − vi| p. (3.22)
i=1 i=1

When p is even, |x| p = x p , so this reduces to a single polynomial inequality of degree

p, giving s = 1 and deg = p. When p is odd, each term |zi − wi | p can be handled by
introducing auxiliary sign variables that total s = O(22d ) predicates [AB99]. However,
since s appears inside a logarithm in Theorem 3.1.15, the 2d term becomes linear in d.
In both cases, the bound gives V ≤ 4d log(8e · p · s) = O(d log d).

Remark 3.1.17. Theorem 3.1.14 provides an asymptotic bound on the number of

samples; a practical implementation of the algorithm requires an explicit, finite sample
size. This exact number is detailed by Talagrand [Tal94]. Likewise, the precise upper
bound on the VC dimension is derived in Goldberg et. al. [GJ95].
We now combine these ideas to analyse how sampling affects the overall algorithm.
In each iteration, we replace the search space with a random sample, compute an
ℓ p -centerpoint of that sample, and use it to discard part of the search space. In each
iteration k, we will draw a new sample from the current search space Sk and compute an
approximate α-ℓ p -centerpoint for that sample. Let T denote the maximum total number
of queries the algorithm needs to take to find an approximate fixed point.
Recall that the ε-approximation guarantee from Theorem 3.1.14 has a failure probability.
Because we want the algorithm to be successful across all T queries, we set the per-
round failure probability to qk = Tq , so the entire procedure fails with a probability of at
most q.
To determine the total number of queries T required to find an ε f p -fixed point of a λ
contraction map under the ℓ p -norm, we consider the volume of the search space. By
querying an α-ℓ p -centerpoint at each step, the algorithm is guaranteed to discard at
least an α fraction of the current search space volume. After T queries, the remaining
volume is bounded by (1 − α)T vol(S0 ).
Haslebacher et al. [HLSW25] established that their algorithm must terminate before the
ε p −ε f p λ
search space volume drops below the volume of an ℓ p ball of radius r = f 2+2λ . Let
Vtarget = vol(B p (0, r)). To guarantee termination, we require:

Vtarget ln(vol(S0 )) − ln(Vtarget )

(1 − α)T ≤ =⇒ T ≥ . (3.23)
vol(S0 ) − ln(1 − α)

1
Lastly, we establish the total number of samples we require to achieve a 2(d+1) -
centerpoint.
Chapter 3. Theoretical Advancements 24

p
Theorem 3.1.18. Let X = Rd , and let R be the family of all bisector ℓ p -halfspaces Hu,v
in Rd . In each iteration k ∈ {1, . . . , T }, let µk be the uniform measure on the search
1
space Sk . Fix εapprox = 2(d+1) and an overall failure probability q ∈ (0, 1). If in each
round we draw a sample S̃k of size
!
VCdim(R ) + log(T /q)
N=O (3.24)
ε2approx
i.i.d. from µk , then with probability at least 1 − q over all T rounds simultaneously,
every discrete ℓ p -centerpoint ck of the sample S̃k is an α-ℓ p -centerpoint of the measure
1
µk with α = 2(d+1) .

Proof. By Proposition 2.2, the VC dimension of the range space (X, R ) of bisector
ℓ p -halfspaces is V = O(d log d). It is easy to see from the definition of VC dimension
that restricting the domain from X to Sk ⊆ X does not increase the VC dimension.
By setting the per-round failure probability to qk = q/T , Theorem 3.1.14 guarantees
that S̃k is an εapprox -approximation for (Sk , RSk ) with probability at least 1 − q/T . All T
samples are valid approximations simultaneously with probability at least 1 − q.
Under this success event, let ck be a discrete ℓ p -centerpoint of S̃k . By definition, for any
p
direction v ∈ Sd−1 , the limit ℓ p -halfspace Hck ,v satisfies:
p
|S̃k ∩ Hck ,v | 1
≥ . (3.25)
|S̃k | d +1
p p
Recall from Observation 3.1.3 that Hck ,v = ε>0 Hck ,ck −εv . This implies that for any
T
p p
ε > 0, the limit ℓ p -halfspace is a subset of the bisector ℓ p -halfspace: Hck ,v ⊆ Hck ,ck −εv .
Consequently:
p
|S̃k ∩ Hck ,ck −εv | |S̃k ∩ Hcpk ,v | 1
≥ ≥ . (3.26)
|S̃k | |S̃k | d +1

Since S̃k is an εapprox -approximation for the range space of bisector ℓ p -halfspaces, we
have for any ε > 0:
p
p |S̃k ∩ Hck ,ck −εv | 1 1 1
µk (Hck ,ck −εv ) ≥ − εapprox ≥ − = . (3.27)
|S̃k | d + 1 2(d + 1) 2(d + 1)

To find the measure of the limit halfspace, consider a decreasing sequence εn → 0+ . As

εn decreases, the bisector halfspaces form a nested sequence:
p p
Hck ,ck −ε1 v ⊇ Hck ,ck −ε2 v ⊇ . . . , (3.28)
p
and their intersection is exactly the limit ℓ p -halfspace Hck ,v .
By applying the continuity of measure from above to this countable sequence, we
obtain:
!
∞
\ p p 1
µk (Hcpk ,v ) = µk Hck ,ck −εn v = lim µk (Hck ,ck −εn v ) ≥ . (3.29)
n=1
n→∞ 2(d + 1)
Chapter 3. Theoretical Advancements 25

1
Thus, ck is an α-ℓ p -centerpoint with α = 2(d+1) for all v ∈ Sd−1 .

Substituting our previous results into Equation 3.24, it is easy to see we only need
a polynomial number of samples on d, log(q−1 ) and log((1 − λ)−1 ) to find a good
approximate centerpoint.

[Link].3 Can we sample uniformly from the search space, and if so, how? Sam-
pling from our search space Sk is non-trivial. While it is possible to sample efficiently
and uniformly from convex sets [DFK91], there is no known method to sample from
general sets. Our search space Sk is the intersection of bisector ℓ p -halfspaces and need
not be convex. In fact, the intersection of ℓ p -halfspaces may not even be connected in
general.
One possible way forward is to relax convexity to a weaker geometric property known
as star-shapedness.
Definition 3.1.19. A nonempty set K ⊂ Rn is star-shaped if there exists a point x0 ∈ K,
called a star center, such that for every x ∈ K and every t ∈ [0, 1], the point (1 −t)x0 +tx
lies in K. Equivalently, K is star-shaped if the entire line segment from x0 to any point
in K stays inside K. (See Figure 3.4a for a visualisation).
Definition 3.1.20. The kernel of a star shaped set K ⊂ Rn , denoted ker(K), is the set of
all star centres:

ker(K) := {x0 ∈ K | (1 − t)x0 + tx ∈ K for all x ∈ K, t ∈ [0, 1]} . (3.30)

Establishing that the search space is star-shaped would enable the use of existing
algorithms for uniform sampling, each with different assumptions and computational
complexity trade-offs:
• Chandrasekaran et al. [CDV10] developed an algorithm for sampling star-shaped
sets whose complexity scales polynomially with the dimension d, the inverse
error tolerance 1/ε, and the inverse kernel fraction 1/η. In this context, ε bounds
the “distance” from the true uniform distribution, while η represents the ratio of
the kernel’s volume to the total volume of the star-shaped set.
• Abbasi-Yadkori et al. [AYBGM17] showed that, provided there exists a well-
behaved1 map from a convex set to the star-shaped set, they can sample uniformly
in O(d 6 ) steps.
The question is whether our specific search space Sk , constructed via Algorithm 1, is
star-shaped.
To understand why this could be the case, we can interpret a bisector ℓ p -halfspace as a
Voronoi cell.
1They require a measure-preserving, bi-Lipschitz map.
Chapter 3. Theoretical Advancements 26

Definition 3.1.21. Let P = {p1 , p2 , . . . , pn } be a set of n distinct points in Rd , referred

to as sites2 . The ℓ p Voronoi cell Vor p (pi , P) associated with the site pi ∈ P is defined as:

Vor p (pi , P) = {x ∈ Rd | ∥x − pi ∥ p ≤ ∥x − p j ∥ p , ∀p j ∈ P} (3.31)

(See Figure 3.4b for a visualisation).

Lemma 3.1.22. For any p ∈ [1, ∞] and any finite set P ⊂ Rd , each ℓ p Voronoi cell
Vor p (pi , P) is star-shaped with respect to its site pi .

Proof. Fix x ∈ P and let y ∈ Vor p (x, P). For any t ∈ [0, 1], define yt := (1 −t)x +ty. We
will show that yt ∈ Vor p (x, P). Note that since yt lies on the line segment connecting x
and y, the distance is additive:

∥y − x∥ p = ∥(y − yt ) + (yt − x)∥ p = ∥y − yt ∥ p + ∥yt − x∥ p . (3.32)

Let z ∈ P be arbitrary. We apply the triangle inequality to the points y, yt , and z:

∥y − z∥ p ≤ ∥y − yt ∥ p + ∥yt − z∥ p . (3.33)

Since y ∈ Vor p (x, P), we know that ∥y − x∥ p ≤ ∥y − z∥ p . Combining this with the
triangle inequality above gives:

∥y − x∥ p ≤ ∥y − yt ∥ p + ∥yt − z∥ p . (3.34)

Now, substituting the identity from Equation 3.32 into the left-hand side:

∥y − yt ∥ p + ∥yt − x∥ p ≤ ∥y − yt ∥ p + ∥yt − z∥ p . (3.35)

Subtracting the term ∥y − yt ∥ p from both sides yields:

∥yt − x∥ p ≤ ∥yt − z∥ p . (3.36)

Since z ∈ P was arbitrary, we conclude that yt ∈ Vor p (x, P).

Generally, the intersection of star-shaped sets is not necessarily star-shaped. However, a

well-known sufficient condition is that if all the sets being intersected share a common
star centre, their intersection is star-shaped [HHMM20].
The fixed point x∗ lies in every bisector ℓ p -halfspace used in the construction of Sk .
Therefore, if x∗ lies in the kernel of each bisector ℓ p -halfspace, then Sk is guaranteed to
be star-shaped.
For certain classes of contraction maps, it can be shown that the exact fixed point x∗
is a star centre. However, this property does not hold in general (see Figure A.2), and
establishing star-shapedness of the search space remains an open question.
2 Note pi is unrelated to the p in the ℓ p -norm.
Chapter 3. Theoretical Advancements 27

(a) star-shaped set and kernel (b) ℓ∞ Voronoi diagram

Figure 3.4: Illustrations of a star-shaped set and an ℓ∞ Voronoi diagram. (a) A star-
shaped set K together with its kernel ker(K), shown as the shaded region. (b) An ℓ∞
Voronoi diagram for a set of sites (black dots). Each Voronoi cell Vor p (pi , P) is star-
shaped with respect to its generating site pi , even though the cells may be non-convex.

Lemma 3.1.23. Let f : Rd → Rd be a monotone ℓ1 -contraction. If there exists x ∈ Rd

with x ≤ f (x) or x ≥ f (x) component-wise, and x∗ is a fixed point of f , then x∗ lies in
the kernel of the bisector ℓ1 -halfspace Hx,1 f (x) .

Proof. Let y = f (x). The condition z ∈ Hx,1 f (x) is equivalent to Φ(z) ≤ 0, where Φ is
defined as:
d
Φ(z) := ∥z − y∥1 − ∥z − x∥1 = ∑ ψi (zi ), (3.37)
i=1
with ψi (t) = |t − yi | − |t − xi |.
Without loss of generality, assume x ≤ y, which means xi ≤ yi for all i = 1, . . . , d.
Consider ψi (t): 
yi − xi
 if t < xi ,
ψi (t) = yi + xi − 2t if xi ≤ t ≤ yi , (3.38)

xi − yi if t > yi .


The function ψi (t) is monotonically non-increasing for all t ∈ R. Furthermore, it

achieves its global minimum value (xi − yi ) for all t ∈ [yi , ∞).
Let z ∈ Hx,1 f (x) (so Φ(z) ≤ 0). Let w = (1 − η)z + ηx∗ for some arbitrary η ∈ [0, 1].
We know that x∗ ≥ y, so for every i, xi∗ ≥ yi . This means xi∗ is in the region where ψi is
at its global minimum. Comparing ψi (wi ) to ψi (zi ):
• Case A: If zi ≤ xi∗ , then zi ≤ wi ≤ xi∗ . Since ψi is non-increasing: ψi (wi ) ≤ ψi (zi ).
• Case B: If zi > xi∗ , then zi > xi∗ ≥ yi . Since both zi and xi∗ are in the region [yi , ∞)
where ψi is constant, ψi (zi ) = ψi (wi ) = ψi (xi∗ ). Thus, ψi (wi ) ≤ ψi (zi ) holds
trivially.
Chapter 3. Theoretical Advancements 28

In both cases, ψi (wi ) ≤ ψi (zi ). Summing over i:

d d
Φ(w) = ∑ ψi (wi ) ≤ ∑ ψi (zi ) = Φ(z). (3.39)
i=1 i=1

Since z ∈ Hx,1 f (x) , we have Φ(z) ≤ 0, which implies Φ(w) ≤ 0. Therefore, w ∈ Hx,1 f (x) .
Since w was arbitrary, the segment [z, x∗ ] ⊆ Hx,1 f (x) .

Remark 3.1.24. Lemma 3.1.23 is not particularly helpful in practice. Though it may be
the case that for some functions there is a theoretical guarantee that a chosen centerpoint
c satisfies c ≤ f (c) or c ≥ f (c) element-wise, in general, an ℓ p -centerpoint that satisfies
either condition may not exist. Furthermore, if we were to enforce x ≤ f (x) or x ≥ f (x)
across the entire domain of f , binary search would efficiently yield the approximate
fixed point [CLT08].
If Sk is shown to be star-shaped in general, algorithms like those by Chandrasekaran
et al. [CDV10] could be directly applied (provided we meet their other kernel-related
requirements). Alternatively, even if Sk is not star-shaped, future research could seek a
well-behaved transformation that maps a convex set to the search space (as required by
Abbasi-Yadkori et al. [AYBGM17]).

3.1.2 Decomposition Algorithm

Several of the existing algorithms for finding fixed points of either monotone functions
or contractions work by partitioning the d coordinates {x1 , . . . , xd } into two sets, which
we call blocks, of sizes a and b [SS03, FGMS20, DQY24, CL22, FPS22, CLY26]. For
a block B ⊆ {1, . . . , d}, we say “solving block B” to mean computing the coordinates
(xi )i∈B of a fixed point of f while keeping the other coordinates fixed.
The algorithms typically proceed by fixing (i.e. treating as constants) the value of
block b to solve block a. However, because block b is fixed to arbitrary values, the
resulting solution is “conditional”: even on block a, it does not need to coincide with
the corresponding coordinates of the fixed point of f , as it depends on the chosen values
of block b. While this partial solution may provide useful information, finding the
global fixed point of f requires further steps. Typically, if T (d, ·) denotes the time it
takes to solve the fixed point problem in d dimensions (ignoring dependence on other
parameters), this decomposition leads to a running time on the order of T (a, ·) T (b, ·).
This will become clearer later in this section.
The reason why solving block a alone does not suffice is that its solution depends on
the fixed values of block b. For this reason, we investigate what happens if we can
split the coordinates in such a way that the coordinate functions of f corresponding to
block a do not depend on block b, specifically in the setting of ℓ p -contractions. If this
is the case, block a can be solved once using arbitrary fixed values for b. No matter
what values are chosen for b, the resulting solution for a will be the same. We can then
proceed to solve block b while keeping block a fixed at its fixed point. We generalise
this idea to decomposing the problem into k blocks; the smaller these blocks are, the
Chapter 3. Theoretical Advancements 29

faster the algorithm runs. Overall, the running time is of the form ∑ki=1 T (|blocki |, ·), as
each block is solved once.
However, it may be the case that the coordinate functions of f corresponding to one
block depend on the coordinates of another block, and vice versa. In this case, the
sequential method above cannot be applied. To address this, we select a subset of
coordinates, which we call a “separator” block, such that once these coordinates are
fixed, the remaining components of f no longer have a mutual dependence. This
allows the remaining coordinates to be decomposed into blocks that can be solved
sequentially, as before. Because the components of f corresponding to the separator
may still depend on the other coordinates, we incur a higher computational cost than
if it was unidirectional. Overall, it results in an overall time complexity of the form
T (|separator|, ·) ∑ki=1 T (|blocki |, ·). Although for general graphs this means the worst-
case runtime does not improve upon that of Chen et al. [CLY26], we will later introduce
certain classes of functions for which this method achieves a worst-case sub-exponential
runtime.
We now formalise these ideas, starting with the notion of a dependency graph. Let d ∈ N
and let f : [0, 1]d → [0, 1]d be a ℓ p -contraction with contraction coefficient λ ∈ (0, 1),
for some p ∈ [0, ∞].
Definition 3.1.25. The dependency graph of f is the directed graph G = (V, E) with
vertex set V = {1, 2, . . . , d}, where there is a directed edge (i, j) ∈ E if and only if i ̸= j
and the i-th coordinate function fi depends on x j . Formally, (i, j) ∈ E if i ̸= j and there
exist x, y ∈ [0, 1]d with xk = yk for all k ̸= j such that fi (x) ̸= fi (y). See Example A.0.1.

[Link] Strongly Connected Component Decomposition

In the dependency graph, vertices may form part of cycles. These correspond to groups
of coordinates that depend on each other in a cyclic way. To properly define these
mutually dependent coordinates, we use the notion of strongly connected components.
Definition 3.1.26. A strongly connected component (SCC) of a directed graph G =
(V, E) is a maximal set of vertices C ⊆ V such that for every pair u, v ∈ C, there exists a
directed path from u to v and from v to u.
SCCs are groups of variables that are interdependent and therefore cannot be separated
as earlier described. If we abstract each SCC, and view it as a single vertex, we obtain
a higher-level view of the dependency graph in which all cycles are contained within
individual “condensed” components (see Figure 3.5).
Definition 3.1.27. A directed acyclic graph (DAG) is a directed graph G = (V, E) that
contains no directed cycles. That is, there is no sequence of vertices v1 , . . . , vk with
k > 1 such that (vi , vi+1 ) ∈ E for all 1 ≤ i < k and v1 = vk .
It is a well-known result that any directed graph G can be uniquely decomposed into
its SCCs C = {C1 , . . . ,Ck }. Furthermore, the condensation graph GSCC is defined as
follows: each Cr ∈ C is represented by a single vertex in GSCC , and we include a
directed edge from Cr to Cs if and only if there exist vertices u ∈ Cr and v ∈ Cs such
Chapter 3. Theoretical Advancements 30

that (u, v) is an edge in G. By construction, all cycles are within individual SCCs, and
hence GSCC is a DAG [CLRS22].
The advantage of a DAG is that we can order its nodes so that all edges go from
earlier nodes to later ones. This allows us to solve the fixed point problem for the
components of f corresponding to each SCC in this order, without recomputing any of
them unnecessarily. This ordering is called a topological ordering.
Definition 3.1.28. A topological ordering of a DAG is an ordering of its vertices such
that for every edge (u, v) ∈ E, vertex u appears before v. Formally, it is a bijection
σ : V → {1, . . . , |V |} satisfying σ(u) < σ(v) for all (u, v) ∈ E. A directed graph admits
a topological ordering if and only if it is a DAG [CLRS22].

(a) Dependency graph (b) Condensation graph

Figure 3.5: SCC decomposition of the system in Example A.0.1. (a) The dependency
graph G. Solid edges denote the dependencies within the same SCC. The dashed
edges represent dependencies between different components. (b) The condensation
graph GSCC obtained by collapsing each SCC into a single node. The resulting graph is
a DAG.

Because GSCC is a DAG, we can find its topological ordering. For our purposes, SCC u
appears before SCC v in the order if v depends on u.
This topological ordering forms the sequence of blocks we must solve one by one, in
order. However, in practice, this introduces a numerical issue. Our algorithms do not
compute exact fixed points, but rather ε-approximate fixed points. Solving one SCC to
a certain tolerance, treating those approximate values as constants, and subsequently
solving the next SCC propagates and magnifies the initial error. To overcome this issue,
we must solve earlier SCCs in the topological order with a smaller tolerance. Before
we go into precisely how much smaller this tolerance needs to be, we must define two
functions.
Definition 3.1.29. Let f : [0, 1]d → [0, 1]d , and let S and C partition V = {1, . . . , d}.
For x ∈ [0, 1]d with sub-vectors xS , xC , and a fixed uS ∈ [0, 1]|S| , we define the mapping
fCuS : [0, 1]|C| → [0, 1]|C| by fCuS (xC ) := ( f (uS , xC ))C (where we optionally reorder the
coordinates so that S precedes C). That is, for each j ∈ C, [ fCuS (xC )] j = [ f (uS , xC )] j .
Chapter 3. Theoretical Advancements 31

Likewise, for a fixed uC ∈ [0, 1]|C| , we define the mapping fSuC : [0, 1]|S| → [0, 1]|S|
analogously, meaning for each i ∈ S, [ fSuC (xS )]i = [ f (xS , uC )]i .
Lemma 3.1.30. Let f : [0, 1]d → [0, 1]d be an ℓ p -contraction with coefficient λ, and let
S and C partition the set V = {1, . . . , d}. For any fixed vector uS ∈ [0, 1]|S| , the mapping
fCuS is an ℓ p -contraction on [0, 1]|C| .

Proof. Let xC and yC be any two vectors in [0, 1]|C| , and define x = (uS , xC ) and y =
(uS , yC ) in [0, 1]d . Because x and y have identical sub-vector uS on the index set S, the
ℓ p -distance between them is only determined by their components in C. For p < ∞:
!1/p
∥x − y∥ p = ∑ |(uS )i − (uS )i| p + ∑ |(xC ) j − (yC ) j | p = ∥xC − yC ∥ p
i∈S j∈C

And for p = ∞:
 
∥x − y∥∞ = max max |(uS )i − (uS )i |, max |(xC ) j − (yC ) j | = ∥xC − yC ∥∞
i∈S j∈C

Since f is an ℓ p -contraction with coefficient λ, we have ∥ f (x) − f (y)∥ p ≤ λ∥x − y∥ p .

The ℓ p -distance between f (x) and f (y) contains the differences though all indices in
V = S ∪C. Since the absolute value raised to the power of p is always non-negative,
discarding the terms corresponding to the index set S yields a lower bound for p < ∞:
!1/p !1/p
∥ f (x) − f (y)∥ p = ∑ | f (x)i − f (y)i| p ≥ ∑ | f (x) j − f (y) j | p
i∈V j∈C

Similarly, for p = ∞, the maximum over a subset of indices yields a lower bound:

∥ f (x) − f (y)∥∞ = max | f (x)i − f (y)i | ≥ max | f (x) j − f (y) j |

i∈V j∈C

The right side of this inequality is the ℓ p -distance ∥ fCuS (xC ) − fCuS (yC )∥ p . Combining
these results yields:

∥ fCuS (xC ) − fCuS (yC )∥ p ≤ ∥ f (x) − f (y)∥ p ≤ λ∥x − y∥ p = λ∥xC − yC ∥ p

Therefore, the mapping fCuS is an ℓ p -contraction.

With the guarantee that each component is a contraction mapping, we can now bound
the “accumulation” of error as we progress through the topological ordering.
Theorem 3.1.31. Let f be an ℓ p -contraction with contraction coefficient λ. To compute
an ε-approximate fixed point of f , split over k SCCs topologically ordered, it is sufficient
ε λ
to solve each r-th SCC to a tolerance of δr ≤ k(K+1) k−r , where K = 1−λ .
Chapter 3. Theoretical Advancements 32

Proof. Since GSCC is a DAG, it can be topologically ordered. Index the components
1, . . . , k according to this order. Let x∗ = (xC∗ 1 , . . . , xC∗ k ) be the exact fixed point and
x̂ = (x̂C1 , . . . , x̂Ck ) our computed approximate fixed points. Writing er = ∥x̂Cr − xC∗ r ∥ p
for the error at step r.
The variables xC1 depend only on themselves, so we solve this self-contained fixed point
problem to tolerance δ1 , giving e1 ≤ δ1 (Lemma 3.1.30 guarantees the mapping fC1 is a
contraction).
Suppose we have computed x̂C1 , . . . , x̂Cr−1 with errors e1 , . . . , er−1 respectively. Let
C<r = ∪r−1
j=1C j denote the earlier components in the order. Because of the topological
ordering, the output of f on Cr does not depend on any later components C>r . Therefore,
we can fix C>r arbitrarily and substitute our approximate earlier fixed points x̂C<r to
x̂
define the mapping fCCr <r on [0, 1]|Cr | .
Let xC∗ r (x̂C<r ) denote the exact fixed point of this mapping, i.e. the unique vector
x̂
satisfying fCCr <r (xC∗ r (x̂C<r )) = xC∗ r (x̂C<r ). We compute an approximation x̂Cr satisfying
∥x̂Cr − xC∗ r (x̂C<r )∥ p ≤ δr .
Applying the triangle inequality:
er = ∥x̂Cr − xC∗ r ∥ p ≤ ∥x̂Cr − xC∗ r (x̂C<r )∥ p + ∥xC∗ r (x̂C<r ) − xC∗ r ∥ p . (3.40)
To bound the second term, we use the fact that both xC∗ r (x̂C<r ) and xC∗ r are fixed points
of their respective mappings, along with the contraction properties of f :
x̂ x∗
∥xC∗ r (x̂C<r ) − xC∗ r ∥ p = ∥ fCCr <r (xC∗ r (x̂C<r )) − fCCr <r (xC∗ r )∥ p (3.41)
r−1
≤ λ ∑ ∥x̂C j − xC∗ j ∥ p + λ∥xC∗ r (x̂C<r ) − xC∗ r ∥ p . (3.42)
j=1
Rearranging:
λ r−1 r−1
∥xC∗ r (x̂C<r ) − xC∗ r ∥ p ≤ ∑ e j = K ∑ e j. (3.43)
1 − λ j=1 j=1

This gives the recurrence er ≤ δr + K ∑r−1 r

j=1 e j . Setting Er = ∑ j=1 e j :

k
Er ≤ δr + (K + 1)Er−1 =⇒ Ek ≤ ∑ δr (K + 1)k−r . (3.44)
r=1

To ensure Ek ≤ ε, it suffices to require each term satisfies δr (K + 1)k−r ≤ ε/k, giving:

ε
δr ≤ . (3.45)
k(K + 1)k−r
Finally, for all p ∈ [1, ∞]:
k
∥x̂ − x∗ ∥∞ ≤ ∥x̂ − x∗ ∥ p ≤ ∑ ∥x̂Cr − xC∗ r ∥ p = Ek . (3.46)
r=1

Thus, Ek ≤ ε is sufficient to guarantee that the overall error is at most ε for any ℓ p -
norm.
Chapter 3. Theoretical Advancements 33

This theorem implies that we can compute the SCC decomposition and its topological
order in O (|V | + |E|) = O (d 2 ) time [Tar72], and subsequently solve each component
to the required accuracy using any existing algorithm.
Let Tbase (n, λ, ε) denote the time required by an existing algorithm to compute an
ε-approximate fixed point of an ℓ p -contraction with contraction coefficient λ. Let
m = ∑kr=1 |Cr | be the total dimension of the entire system, and let Cmax = maxr |Cr |
denote the size of the largest SCC. The total time taken to solve the fixed point of f is
the time taken to find the decomposition plus the sum of the times required to solve
each SCC:
k
(SCC)
Tdecomp (m, λ, ε) ≤ O (m2 ) + ∑ Tbase(|Cr |, λ, δr ). (3.47)
r=1
For the norms of particular interest, i.e. ℓ1 and ℓ∞ , we can use the runtime result from
Chen et al. [CLY26]. To obtain an upper bound on the running time, we substitute k for
k − r in the bound for δr from Theorem 3.1.31 (as k − r < k), and replace the sum over
all components with k instances of the largest component Cmax solved to the tightest
tolerance of δ1 . This yields the bound:
(SCC)

⌈|Cmax |/2⌉ 
Tdecomp (m, λ, ε) ≤ O m2 + k log(ε−1 ) + k log(K + 1) . (3.48)

For graphs where the largest SCC is relatively small, this achieves an improvement over
directly applying the algorithm from Chen et al. [CLY26].

[Link] Separator-Based Decomposition

The SCC decomposition provides an efficient solution when Cmax is small. However,
we need to consider the scenario where the dependency graph contains a large SCC that,
by definition, can’t be decomposed into smaller SCCs. In such cases, we can apply a
separator-based decomposition.
The idea is to partition the dependency graph G into a subset of variables S ⊂ V we call
the “separator”, and blocks C1 ,C2 , . . . ,Ck ⊂ V \ S, such that removing S causes the sub-
graph G[V \ S] to form a DAG over these blocks (see Figure A.3). Let C = kr=1 Cr . We
S

index the blocks in C in topological order so that we can write x = (xS , xC1 , . . . , xCk ) =
(xS , xC ) and decompose f as:
 
fS (xS , xC )
f (x) = . (3.49)
fC (xS , xC )

Because G[V \ S] is a DAG over these blocks, each block Ci depends only on the
separator S and on earlier blocks in the topological order.
To make use of this structure, we temporarily fix the values of the separator variables
to some uS ∈ [0, 1]|S| . Then, to find the fixed point of fCuS we can apply the DAG
decomposition as described before.
By Lemma 3.1.30 (along with Banach’s Fixed Point Theorem), we know that for
any given uS fCuS has a unique fixed point. We denote this exact fixed point, and its
Chapter 3. Theoretical Advancements 34

corresponding δC -approximation output by our solver, as xC∗ (uS ) and x̃C (uS ) respectively.
They satisfy:

xC∗ (uS ) = fCuS (xC∗ (uS )), ∥x̃C (uS ) − xC∗ (uS )∥ p ≤ δC . (3.50)

Knowing that the mapping fCxS has a unique fixed point for any given xS ∈ [0, 1]|S| , we
can express the fixed point for the C coordinates entirely as a function of xS , xC∗ (xS ).
Instead of solving for all d coordinates of f , we can substitute this function xC∗ (xS ) into
fS and treat it as a black-box for now. By doing so, we reduce the problem from d to |S|
dimensions.
To turn this idea into an algorithm, we introduce two functions F and F̃δC .
The first map, F, is defined by passing in the true, exact fixed point of fCxS , denoted
xC∗ (xS ), into fS :
x∗ (xS )
F(xS ) = fS C (xS ). (3.51)
As we will show, F is an ℓ p -contraction. This guarantees it has a unique fixed point. If
we find this fixed point xS∗ and concatenate it with xC∗ (xS∗ ), it yields the exact fixed point
of the full system f .
In practice, we cannot compute this exact fixed point xC∗ (xS ). Therefore, we define the
computable counterpart of F, F̃δC , which instead feeds the δC -approximate fixed point,
x̃C (xS ), into fS . Assuming our chosen approximate solver is deterministic, x̃C (xS ) and
consequently F̃δC are well-defined functions:
x̃ (xS )
F̃δC (xS ) = fS C (xS ). (3.52)

Because it relies on an approximation, F̃δC behaves as an ℓ p -contraction over most of

the space, but as we will soon show, it loses this property within a small radius.
Definition 3.1.32. A map g : [0, 1]d → [0, 1]d is a (λ, γ)-contraction under the ℓ p -norm
if there exist λ ∈ (0, 1) and γ ≥ 0 such that for all x, y ∈ [0, 1]d :

∥g(x) − g(y)∥ p ≤ λ∥x − y∥ p + γ. (3.53)

Proposition 3.1.33. Let g : [0, 1]d → [0, 1]d be a (λ, γ)-contraction under the ℓ p -norm
γ
with λ ∈ (0, 1). Whenever ∥x − y∥ p > 1−λ , the function g is a contraction under the
ℓ p -norm.

Proof. By the definition of a (λ, γ)-contraction under the ℓ p -norm, for any x, y ∈ [0, 1]d :

∥g(x) − g(y)∥ p ≤ λ∥x − y∥ p + γ (3.54)

Given that the ℓ p distance between the points x and y satisfies:

γ
∥x − y∥ p > (3.55)
1−λ
Since λ ∈ (0, 1), we know that 1 − λ > 0. Multiplying both sides by 1 − λ:

(1 − λ)∥x − y∥ p > γ =⇒ ∥x − y∥ p − λ∥x − y∥ p > γ (3.56)

Chapter 3. Theoretical Advancements 35

We can now substitute this upper bound for γ back into Equation 3.54:
∥g(x) − g(y)∥ p < λ∥x − y∥ p + (∥x − y∥ p − λ∥x − y∥ p ) (3.57)
=⇒ ∥g(x) − g(y)∥ p < ∥x − y∥ p (3.58)
γ
Therefore, g is a contraction when ∥x − y∥ p ≥ 1−λ .

Proposition 3.1.34. For any separator S ⊂ V , the function F as defined in Equation 3.51
is an ℓ p -contraction.

Proof. Let uS , wS ∈ [0, 1]|S| . Recall that xC∗ (uS ) and xC∗ (wS ) are the exact fixed points of
fCuS and fCwS i.e. fCuS (xC∗ (uS )) = xC∗ (uS ) and fCwS (xC∗ (wS )) = xC∗ (wS ).
Define the vectors X = (uS , xC∗ (uS )) and Y = (wS , xC∗ (wS )). Applying f to X and Y
yields f (X) = (F(uS ), xC∗ (uS )) and f (Y ) = (F(wS ), xC∗ (wS )).
Because f is an ℓ p -contraction with contraction coefficient λ, we have ∥ f (X) −
f (Y )∥ p ≤ λ∥X −Y ∥ p . We now consider the cases for p < ∞ and p = ∞ separately.
For 1 ≤ p < ∞, raising both sides to the power of p and expanding the expression:
 
p ∗ ∗ p p p ∗ ∗ p
∥F(uS ) − F(wS )∥ p + ∥xC (uS ) − xC (wS )∥ p ≤ λ ∥uS − wS ∥ p + ∥xC (uS ) − xC (wS )∥ p .
(3.59)
p
We distribute λ p on the right side and subtract λ p ∥xC∗ (uS ) − xC∗ (wS )∥ p from both sides:
∥F(uS ) − F(wS )∥ pp + (1 − λ p )∥xC∗ (uS ) − xC∗ (wS )∥ pp ≤ λ p ∥uS − wS ∥ pp . (3.60)
Since λ ∈ (0, 1), we know that (1 − λ p ) > 0. Because norms are non-negative, the entire
p
term (1 − λ p )∥xC∗ (uS ) − xC∗ (wS )∥ p is non-negative. Dropping a non-negative term from
the left side:
∥F(uS ) − F(wS )∥ pp ≤ λ p ∥uS − wS ∥ pp . (3.61)
Taking the p-th root of both sides:
∥F(uS ) − F(wS )∥ p ≤ λ∥uS − wS ∥ p . (3.62)
For p = ∞, we can write:
 
∗ ∗
max ∥F(uS ) − F(wS )∥∞ , ∥xC (uS ) − xC (wS )∥∞ (3.63)
 
≤ λ max ∥uS − wS ∥∞ , ∥xC∗ (uS ) − xC∗ (wS )∥∞ . (3.64)

If the maximum on the right side evaluates to ∥xC∗ (uS ) − xC∗ (wS )∥∞ , the inequality
implies ∥xC∗ (uS ) − xC∗ (wS )∥∞ ≤ λ∥xC∗ (uS ) − xC∗ (wS )∥∞ . Since λ ∈ (0, 1), this requires
∥xC∗ (uS ) − xC∗ (wS )∥∞ = 0, which in turn forces ∥F(uS ) − F(wS )∥∞ = 0 ≤ λ∥uS − wS ∥∞ .
Otherwise, if the maximum on the right side evaluates to ∥uS − wS ∥∞ , meaning:
∥F(uS ) − F(wS )∥∞ ≤ λ∥uS − wS ∥∞ .
In all cases, we obtain ∥F(uS ) − F(wS )∥ p ≤ λ∥uS − wS ∥ p .
Chapter 3. Theoretical Advancements 36

Proposition 3.1.35. For an arbitrary separator S ⊂ V , the function F̃δC as defined in

Equation 3.52, is a (λ, γ)-contraction under the ℓ p -norm with γ = 2λδC .

Proof. Fix xS ∈ [0, 1]|S| . Let Xexact = (xS , xC∗ (xS )) and Xapprox = (xS , x̃C (xS )). Because
their S-coordinates are identical, the ℓ p -distance between these vectors is determined
entirely by the error of the solver on the C-coordinates. By our solver’s guarantee, this
error is bounded by δC :
∥Xexact − Xapprox ∥ p = ∥xC∗ (xS ) − x̃C (xS )∥ p ≤ δC . (3.65)
Next, observe that F(xS ) and F̃δC (xS ) correspond exactly to the S-block components of
the outputs f (Xexact ) and f (Xapprox ). The norm of a full vector is always bounded below
by the norm of any of its constituent blocks. For p < ∞, dropping the non-negative
p p
differences in the C-block yields ∥F(xS ) − F̃δC (xS )∥ p ≤ ∥ f (Xexact ) − f (Xapprox )∥ p . For
p = ∞, the maximum difference over all coordinates is trivially at least the maximum
difference over the S-coordinates. Thus, we can bound the distance between F(xS ) and
F̃δC (xS ):
∥F(xS ) − F̃δC (xS )∥ p ≤ ∥ f (Xexact ) − f (Xapprox )∥ p . (3.66)
Finally, because the map f is an ℓ p -contraction, we can bound the right side of the
previous inequality:
∥ f (Xexact ) − f (Xapprox )∥ p ≤ λ∥Xexact − Xapprox ∥ p ≤ λδC . (3.67)
Putting these inequalities together yields:
∥F(xS ) − F̃δC (xS )∥ p ≤ λδC . (3.68)

Now, let uS , wS ∈ [0, 1]|S| . By the triangle inequality, and using F is an ℓ p -contraction:
∥F̃δC (uS ) − F̃δC (wS )∥ p ≤ ∥F̃δC (uS ) − F(uS )∥ p + ∥F(uS ) − F(wS )∥ p (3.69)
+ ∥F(wS ) − F̃δC (wS )∥ p (3.70)
≤ λδC + λ∥uS − wS ∥ p + λδC (3.71)
= λ∥uS − wS ∥ p + 2λδC .

By partitioning into S and C, we reduce the problem to solving for an approximate fixed
point of F̃δC . Substituting y = xS∗ into ||x − y|| p > 2λδC
1−λ we see that F̃δC is a contraction
on [0, 1]|S| except in an ℓ p -ball centred at the fixed point xS∗ with a radius controlled by
δC .
We now derive how small this δC needs to be for the overall ℓ p -error to be smaller than
ε.
Proposition 3.1.36. Let x∗ ∈ [0, 1]d denote the unique fixed point of f . Suppose x̂S ∈
[0, 1]|S| satisfies ∥x̂S − F̃δC (x̂S )∥ p ≤ δS , and set x̂C = x̃C (x̂S ) such that ∥x̂C − xC∗ (x̂S )∥ p ≤
δC . The combined vector x̂ = (x̂S , x̂C ) satisfies ∥x̂ − x∗ ∥ p ≤ ε provided:
ε(1 − λ) ε(1 − λ)
δS ≤ , δC ≤ , (3.72)
21/p 21/p (1 + λ)
(with the convention that 21/p = 1 when p = ∞).
Chapter 3. Theoretical Advancements 37

Proof. Let f be an ℓ p -contraction, and let x∗ denote its unique fixed point. By the
triangle inequality,

∥x̂ − x∗ ∥ p ≤ ∥x̂ − f (x̂)∥ p + ∥ f (x̂) − f (x∗ )∥ p . (3.73)

Using f (x∗ ) = x∗ and the contraction property,

∥ f (x̂) − f (x∗ )∥ p ≤ λ∥x̂ − x∗ ∥ p , (3.74)

1
=⇒ (1 − λ)∥x̂ − x∗ ∥ p ≤ ∥x̂ − f (x̂)∥ p , hence ∥x̂ − x∗ ∥ p ≤ ∥x̂ − f (x̂)∥ p .
1−λ
(3.75)
It therefore suffices to prove ∥x̂ − f (x̂)∥ p ≤ ε(1 − λ).
For the S-coordinates, since x̂C = x̃C (x̂S ) and F̃δC (x̂S ) = fSx̂C (x̂S ), we have

fS (x̂) = fSx̂C (x̂S ) = F̃δC (x̂S ), (3.76)

and thus
∥x̂S − fS (x̂)∥ p = ∥x̂S − F̃δC (x̂S )∥ p ≤ δS . (3.77)

For the C-coordinates, writing fC (x̂) = fCx̂S (x̂C ) and letting xC∗ (x̂S ) be the fixed point of
fCx̂S , we have

∥x̂C − fC (x̂)∥ p ≤ ∥x̂C − xC∗ (x̂S )∥ p + ∥xC∗ (x̂S ) − fCx̂S (x̂C )∥ p . (3.78)

The first term is bounded by δC . For the second term, using xC∗ (x̂S ) = fCx̂S (xC∗ (x̂S )) and
the contraction property,

∥xC∗ (x̂S ) − fCx̂S (x̂C )∥ p = ∥ fCx̂S (xC∗ (x̂S )) − fCx̂S (x̂C )∥ p ≤ λ∥xC∗ (x̂S ) − x̂C ∥ p ≤ λδC . (3.79)

Hence,
∥x̂C − fC (x̂)∥ p ≤ (1 + λ)δC . (3.80)

If p < ∞, then
p
p
∥x̂ − f (x̂)∥ pp = ∥x̂S − fS (x̂)∥ pp + ∥x̂C − fC (x̂)∥ pp ≤ δS + (1 + λ)δC . (3.81)

Using the bounds in Equation (3.72),

p ε p (1 − λ) p
p ε p (1 − λ) p
δS ≤ , (1 + λ)δC ≤ , (3.82)
2 2
so
∥x̂ − f (x̂)∥ pp ≤ ε p (1 − λ) p , (3.83)
and therefore ∥x̂ − f (x̂)∥ p ≤ ε(1 − λ). If p = ∞, then



∥x̂ − f (x̂)∥∞ = max ∥x̂S − fS (x̂)∥∞ , ∥x̂C − fC (x̂)∥∞ ≤ max δS , (1 + λ)δC ≤ ε(1 − λ),
(3.84)
1/p
using Equation (3.72) with 2 = 1. Substituting into the earlier bound yields ∥x̂ −
x∗ ∥ p ≤ ε, as required.
Chapter 3. Theoretical Advancements 38

So far, we have relied on computing an approximate fixed point of F via F̃δC . As

established, F̃δC is an ℓ p -contraction only on the domain [0, 1]|S| \ B p (xS∗ , 2δCλ
1−λ ) (where
p
B (c, r) denotes the ℓ p ball centred at c with radius r). Because this domain is not
complete, Banach’s Fixed Point Theorem does not apply.
However, any fixed-point algorithm works by querying points in [0, 1]|S| . By setting
the δC sufficiently small, the non-contractive region is contained within our desired
tolerance ball, meaning B p (xS∗ , 2δCλ p ∗
1−λ ) ⊂ B (xS , δS ).
Because of this, if the algorithm queries a point outside our target tolerance, it operates
in a region where F̃δC is an ℓ p -contraction, and it can proceed as usual. If the algorithm
queries a point that falls inside the non-contractive region, it satisfies the condition
∥F̃δC (xS ) − xS ∥ p ≤ δS . The algorithm will immediately halt and return xS as the δS -
approximate fixed point.
Therefore, the lack of a contraction on all of [0, 1]|S| does not affect the algorithm’s
correctness, and any existing algorithm designed for standard ℓ p -contractions can be
applied to F̃δC .
To quantify how the size of a separator impacts the runtime, we introduce the notion of
a balanced vertex separator.
Definition 3.1.37. Let G = (V, E) be a graph with |V | = d. For a constant α ∈ (0, 1),
an α-balanced vertex separator is a subset S ⊆ V such that every SCC of the sub-graph
G[V \ S] has size at most αd.
Intuitively, removing S breaks the graph into pieces, each of which is significantly
smaller than the original graph. Since every evaluation of F̃δC requires computing an
approximate fixed point on the entire C subsystem (which is relatively easy to solve as
it’s made up of smaller components), it is therefore desirable to keep the separator S as
small as possible.
Let M(d, α) denote the worst-case size of the smallest α-balanced separator over all
directed graphs on d vertices. Formally, M(d, α) = maxG:|V |=d minS |S|, where the
minimum is taken over all α-balanced separators S of G. We also let Tsep (d, α) denote
the time required to compute such a separator.
To solve for the approximate fixed point of F, we apply some base fixed-point algorithm
(passing in F̃δC ). Recall that Tbase (n, λ, ε) denotes the time complexity for this base
algorithm on an n-dimensional problem.
In the worst case, finding the approximate fixed point of F̃δC takes at most Tbase (· · · )
queries (though we could also apply this decomposition algorithm recursively to F̃δC ,
which would replace Tbase with Tdecomp ). Because F̃δC relies on the approximate fixed
points of the remaining blocks Ci ⊂ V \ S, every single evaluation of F̃δC requires
sequentially solving these k sub-problems. Maximizing over all graphs on d vertices
and their valid separators yields the following recurrence relation for the total worst-case
Chapter 3. Theoretical Advancements 39

time Tdecomp (d, λ, ε):

" #
k
Tdecomp (d, λ, ε) ≤ Tsep (d, α) + max Tbase (M(d, α), λ, δS ) ∑ Tdecomp (|Cr |, λ, δC ) .
G,S r=1
(3.85)
The base case for a 1-dimensional system is bounded by Tdecomp (1, λ, ε) = O (log(ε−1 )).
Unfortunately, without additional assumptions on the dependency graph, this recursion
yields a worst-case exponential running time.
Proposition 3.1.38. There exist directed graphs on d vertices for which every α-
balanced directed vertex separator must contain at least d(1 − α) vertices. Therefore,
the worst-case minimal separator size is M(d, α) ≥ d(1 − α). Consequently, with-
out additional assumptions, the separator-based decomposition yields a worst-case
exponential running time of T (d) = 2Ω(d) .

Proof. First, we establish a bound on the separator size. Consider a complete directed
graph Kd = (V, E), where |V | = d and bidirectional edges exist between all pairs of
vertices. By definition, any sub-graph of Kd is also a complete directed graph and forms
a single SCC.
Let S be a valid α-balanced separator for Kd . The remaining vertices V \ S form an
sub-graph of size d − |S|, which is a single SCC. To satisfy the α-balanced criteria, the
size of this remaining SCC cannot exceed αd:

d − |S| ≤ αd =⇒ |S| ≥ d(1 − α). (3.86)

Because α ∈ (0, 1) is a constant, 1 − α is strictly positive. Thus, M(d, α) ∈ Ω(d).

Substituting this bound into our runtime recurrence yields T (d, ·) ≥ T ((1−α)d, ·)T (αd, ·).
Let T ′ (d) = log T (d, ·). This transformation yields:

T ′ (d) ≥ T ′ ((1 − α)d) + T ′ (αd). (3.87)

If we expand this recurrence as a recursion tree, the problem of size d splits into two
subproblems at each step. Because the sizes of these subproblems sum exactly to d
((1 − α)d + αd = d), the total “mass” of the input is conserved across every level of the
tree. Consequently, if we expand the tree all the way down to its leaves, the sum of the
sizes of all the leaves will exactly equal d.
Since T ′ (1) ≥ c0 , we have that T ′ (d) = Ω(d). Finally, by exponentiating both sides to
invert our initial transformation, we obtain exponential time complexity for the original
algorithm:
′
T (d) = 2T (d) = 2Ω(d) . (3.88)

Because some dependency graphs make this decomposition slow, we must restrict
the class of graphs we consider. In the following section, we will parametrise the
dependency graph by a property called “treewidth”, which limits the maximum size a
separator can reach, allowing us to escape this exponential worst-case bound.
Chapter 3. Theoretical Advancements 40

[Link] Complexity Bounds via Treewidth

We begin by defining tree decompositions and treewidth.

Definition 3.1.39. A tree decomposition of an undirected graph G = (V, E) is a tree
T = (X, F) whose nodes Xi ∈ X (called bags) are subsets of V satisfying:
1. Every vertex v ∈ V appears in at least one bag.
2. For every edge (u, v) ∈ E, there exists a bag containing both u and v.
3. For every vertex v, the set of bags containing v forms a connected subtree of T .
The width of a decomposition is maxi |Xi | − 1, and the treewidth tw(G) is the minimum
width over all such decompositions. See Figure A.4 for a visualisation.
To use treewidth as a parameter for our directed dependency graph, we must change its
definition. Treewidth is a property of undirected graphs; we consider the undirected
version of the graph, where all directed edges are replaced by undirected ones. Then,
duplicate edges arising from this change are removed. Then we define w = tw(G) to be
the treewidth of this undirected graph.
A possible alternative would be to replace treewidth with other graph properties such as
directed treewidth [JRST01] or DAG-width [BDH+ 12], which may yield tighter bounds
on separator size, as they account for direction. However, we continue with treewidth
as it is a more studied property.
Lemma 3.1.40 ( [CFK+ 15]). Let G be a graph with treewidth w. Then there exists a
1
2 -balanced vertex separator S with |S| ≤ w + 1. Moreover, given a tree decomposition
of G of width w and size d, such a separator can be found in O (dw) time.
For Theorem ??, we have assumed a tree decomposition is given. If this were not the
case, we would need to find such a decomposition. Some algorithms can find a separator
of size at most 2w + 2 (or a certificate that 2w + 2 ≥ d) in time O(d2w ) [Kor23]. This
means if w is relatively small with respect to d (w ≤ ⌊ d2 −1⌋), we can compute separators
recursively on sub-problems, after which Theorem ?? applies.
Finally, this decomposition offers a potential advantage: when multiple contraction
maps share a common dependency graph, we can pre-compute the tree decomposition of
size at most w + 1 using Ω(cd )-time exact algorithms (as opposed to the approximation
mentioned before) [Tod08]. Then, we can reuse this decomposition to calculate the
fixed points more efficiently for different functions fi . One example where this arises
is in the setting of a criminal attempting to flee a city (as mentioned in Section 2.2.2),
where the road network and hence the dependency graph is fixed, even when the police
blockades and traffic jams occur at different junctions.
√ Furthermore, road networks with
d intersections typically have a treewidth of O( d) [RS86], meaning this decomposition
likely works well for this scenario.
Chapter 3. Theoretical Advancements 41

3.2 Algorithms for monotone ℓ p contractions

As we have seen, there is currently no known efficient algorithm for finding approximate
fixed points of ℓ p -contractions. We now add the condition to f that it must be monotone
on [0, 1]d . This section investigates how we can make use of the combination of these
two properties. We develop a method that, in favourable cases, allows us to discard
large portions of the space where the fixed point could lie.

3.2.1 Bisection Algorithm Description

Let L,U ∈ Rd be lower and upper bounds such that L ≤ x∗ ≤ U, where x∗ is the unique
fixed point of f . For a chosen coordinate i ∈ {1, . . . , d} and value mi ∈ [Li ,Ui ], we aim
to determine the position of xi∗ relative to mi .
Proposition 3.2.1. Fix p ∈ [1, ∞]. Let L,U ∈ Rd with L ≤ U (i.e. Li ≤ Ui for all
i), and define B = {z ∈ Rd | L ≤ z ≤ U}. Let f : B → B be a monotone map that is
ℓ p -contracting with coefficient λ ∈ (0, 1), and let x∗ be its unique fixed point in B.
Fix an index i ∈ {1, . . . , d} and a scalar mi ∈ [Li ,Ui ]. Define the vectors x(L) , x(U) ∈ B
componentwise as:
( (
(L) m i if j = i, (U) mi if j = i,
xj = xj = (3.89)
L j if j ̸= i, U j if j ̸= i.
Then:
1. If fi (x(L) ) ≥ mi , then xi∗ ≥ mi . 2. If fi (x(U) ) ≤ mi , then xi∗ ≤ mi .

Proof. We only prove the first case as the second is symmetric. Construct the vector
y ∈ Rd by replacing the i-th component of x∗ with mi :
(
mi if j = i,
yj = ∗ (3.90)
x j if j ̸= i.
Since x∗ ∈ B, we have L j ≤ x∗j ≤ U j for all j ̸= i, and yi = mi ∈ [Li ,Ui ] by assumption,
so y ∈ B. Comparing y with the vectors x(L) and x(U) at each coordinate: for j ̸= i,
(L) (U)
x j = L j ≤ x∗j = y j ≤ U j = x j , and at j = i all three agree. Hence x(L) ≤ y ≤ x(U) .
Since f is monotone: fi (x(L) ) ≤ fi (y) ≤ fi (x(U) ). Since fi (x(L) ) ≥ mi we have fi (y) ≥ mi ,
i.e. fi (y) − mi ≥ 0.
Let δ = xi∗ − mi . Because x∗ = f (x∗ ), we write
δ = fi (x∗ ) − mi = fi (x∗ ) − fi (y) + fi (y) − mi .



(3.91)
Since x∗ and y differ only at coordinate i, ∥x∗ − y∥ p = |δ|. The contraction property
gives
| fi (x∗ ) − fi (y)| ≤ ∥ f (x∗ ) − f (y)∥ p ≤ λ|δ|, (3.92)
so fi (x∗ ) − fi (y) ≥ −λ|δ|. Combined with fi (y) − mi ≥ 0, δ ≥ −λ|δ|. Suppose that we
have δ < 0. Then |δ| = −δ, and the inequality becomes (1 − λ)δ ≥ 0. Since λ < 1 and
δ < 0, this is a contradiction. Hence δ ≥ 0, i.e. xi∗ ≥ mi .
 Chapter 3. Theoretical Advancements 42

This result relies on both the monotonicity and ℓ p -contraction properties. Neither
condition alone suffices (see Example A.0.2).

Algorithm 2: Simple Iteration + Heuristic. See Figure ?? for a visualisation.

Input: Monotone contraction f , bounds L,U, tolerance ε
Output: Approximate fixed point x
1 while ∥U − L∥ p > ε do
2 changed ← false;
3 for i ← 1 to d do
4 Ignore coordinates that have already converged. C is a constant to ensure
the overall error is bounded by ε, it depends on the norm;
5 if Ui − Li > ε/C then
6 mi ← (Li +Ui )/2;
7 Construct x(L) , x(U) as in Proposition 3.2.1;
8 if fi (x(L) ) ≥ mi then
9 Li ← mi ; changed ← true;
10 end
11 else if fi (x(U) ) ≤ mi then
12 Ui ← mi ; changed ← true;
13 end
14 end
15 end
16 if not changed then
17 L ← f (L); U ← f (U);
18 end
19 end
20 return (L +U)/2;

Based on Proposition 3.2.1, we can develop new algorithms that make use of this
property. The simplest algorithm (see Algorithm 2) looks for the fixed point in [0, 1]d
by setting L = 0 and U = 1. Then, an index that satisfies conditions 1. or 2. is chosen.
If all coordinates are “ambiguous,” i.e. fi (x(L) ) < mi < fi (x(U) ) for all i, then we cannot
bisect any of them, and we apply f to U and L, which by the contraction property will
shrink the box slightly and never exclude the fixed point (see Figure ??).
One natural choice for mi is Ui −L i
2 , but it is arbitrary. This raises the question of whether
a better choice of mi could improve performance.
Proposition 3.2.2. There exists a monotone ℓ∞ -contracting f : [0, 1]d → [0, 1]d , with
contraction coefficient λ ∈ (0, 1), for which Algorithm 2 using mi = Ui −L 2
i
requires
1
Ω( 1−λ ) steps to converge. Furthermore, even if a different test value mi is chosen to
avoid ambiguity, the resulting reduction in search space volume is bounded above by
O (1 − λ).
Chapter 3. Theoretical Advancements 43

Proof. Define f : [0, 1]d → [0, 1]d by:

d
λ
f (x)i =
d ∑ x j + (1 − λ)
j=1

This function is monotone and ℓ∞ -contracting. The unique fixed point is x∗ = 1. Set
L = 0 and U = 1. Consider the first iteration for any coordinate i using mi = 0.5.
Consider fi (x(L) ):
λ
fi (x(L) ) = (0.5) + (1 − λ)
d
For large d and λ ≈ 1, both terms approach zero, yielding fi (x(L) ) ≪ 0.5. Thus,
fi (x(L) ) < mi .
Now, consider fi (x(U) ):
 
(U) λ 0.5 0.5λ
fi (x ) = (d − 0.5) + (1 − λ) = λ 1 − +1−λ = 1−
d d d

Because λ < 1 and d ≥ 1, the subtracted term is strictly greater than zero and less than
0.5. Thus, fi (x(U) ) > 0.5, meaning fi (x(U) ) > mi .
Since fi (x(L) ) < mi < fi (x(U) ), the coordinate i is ambiguous. By symmetry, this holds
for all coordinates.
Because all coordinates are ambiguous, we need to rely on the fallback, i.e. L(k+1) =
f (L(k) ). For any coordinate i, because the sequence of lower bounds L(k) monotoni-
cally approaches the fixed point x∗ = 1 from below, resolving the ambiguity requires
fi (L(k) ) > 0.5.
By the definition of the ℓ∞ -contraction, the distance from the lower bound to the fixed
point after k iterations is bounded by ∥1 − L(k) ∥∞ ≤ λk ∥1 − L(0) ∥∞ . Given the initial
(k)
bound L(0) = 0, the worst-case ℓ∞ distance implies Li ≥ 1 − λk for all i. Applying the
contraction property:
fi (L(k) ) ≥ 1 − λk+1
To guarantee fi (L(k) ) > 0.5, it is sufficient that 1 − λk+1 > 0.5, which simplifies to
λk+1 < 0.5. Taking the natural logarithm of both sides yields:
 
− ln(2) 1
k+1 > =Ω
ln(λ) 1−λ

We now investigate whether some other mi ∈ (0, 1) could guarantee a cut while reducing
the search space by at least a 1/poly(d) fraction. For the algorithm to successfully
update a lower bound without ambiguity, Proposition 3.2.1 requires fi (x(L) ) ≥ mi .
Recall that x(L) is constructed by setting the i-th component to mi and all other compo-
nents to L j = 0. The non-ambiguity condition becomes:

λ
fi (x(L) ) = mi + (1 − λ) ≥ mi
d
Chapter 3. Theoretical Advancements 44

Rearranging to derive an upper bound on mi :

 
λ 1−λ 1−λ
1 − λ ≥ mi 1 − =⇒ mi ≤ =d = O (1 − λ)
d λ
1− d d −λ

Because the initial lower bound is Li = 0, moving the lower bound up to mi removes
exactly a length of mi from the search interval. To achieve a polynomial reduction in
1
the search space per step, we would need to choose a cut value mi ≥ poly(d) . However,
avoiding ambiguity requires mi = O (1 − λ). Therefore, any test value mi requires
1
poly(d) ≤ mi = O(1 − λ) . Conversely, choosing an mi small enough to avoid ambiguity
limits the size of the cut to at most O (1 − λ) per step.

It follows that Algorithm 2 has worst-case

 complexity
 no better than simply using
log(1/ε)
xn+1 = f (xn ), which has a runtime of O log(1/λ) . Therefore, we consider alternative
ways in which Proposition 3.2.1 can be used to accelerate existing fixed-point algorithms
for monotone maps.

3.2.2 Applications to Existing Algorithms

While we have not found a way to leverage Proposition 3.2.1 to improve the worst-case
asymptotic complexity of finding fixed points of monotone ℓ p -contractions, it can still
be used to speed up existing algorithms.
We illustrate this by incorporating it into the algorithm of Dang et al. [DQY24] for
monotone functions on [0, 1]d , under an additional ℓ p -contraction assumption on f .
Note that similar modifications apply to the algorithms of Fearnley et al. [FPS22] and
Chen et al. [CL22].
Their algorithm has lower and upper bound vectors L,U ∈ [0, 1]d defining the current
(0)
search region. At each step, it sets xd = 12 (Ld +Ud ), fixes the d-th coordinate to this
value, and recursively computes a fixed point x∗ of the remaining d − 1 coordinates.
(0)
It then evaluates the d-th coordinate of the original function at (x∗ , xd ), i.e. computes
(0) (0)
[ f (x∗ , xd )]d . If this value is larger than xd , then the fixed point must be greater than
(0) (0) (0)
xd (due to monotonicity), so we update Ld ← xd . Otherwise, we set Ud ← xd . If
equality holds (or is ε close), the fixed point has been found.
When d = 1 the problem can be solved using binary search, and the recursion yields
a time-complexity of O(logd (1/ε)). The highest computational cost comes from the
recursive call. To decide how to bisect along coordinate d, the algorithm must solve a
(d − 1)-dimensional problem, which takes O(logd−1 (1/ε)).
Proposition 3.2.1 allows us to sometimes avoid this step. Given bounds L and U, and
(0)
letting md = xd , we construct x(L) = (L1 , . . . , Ld−1 , md ), x(U) = (U1 , . . . ,Ud−1 , md ),
and evaluate f at these two points. If fd (x(L) ) ≥ md , then the proposition guarantees
that xd∗ ≥ md , so we can discard the lower half without the need for a recursive call.
Similarly, if fd (x(U) ) ≤ md , we can discard the upper half. If neither condition holds,
Chapter 3. Theoretical Advancements 45

we fall back to the original recursive step. We can also apply this check across multiple
coordinates to increase the chance of resolving a coordinate without recursion.
In this way, a constant-time test can, in favourable cases, replace an O(logd−1 (1/ε))
recursive subroutine. This does not change the worst-case complexity, but it can
significantly improve the running time in practice by reducing the number of recursive
calls.
 Chapter 4

Implementation and Results

In this chapter, we test multiple algorithms for computing fixed points. We generate
Condon’s SSGs for testing. The functions arising from these are simple, monotone
and ℓ∞ -contracting. To create instances, we create directed graphs with n nodes. These
nodes have a value, initially set to a uniform random value between 0 and 1, and
a function fi (x). This function reads the values of the nodes’ outgoing neighbours
and applies either max, min, or arithmetic mean, multiplied by a factor λ = 1 − 2−n .
Additionally, two nodes are set to have a constant function fi (x) = 1 and f j (x) = 0,
these are the terminal nodes.
We test two graph families.
Family A (random SSGs). For each size n ∈ {1, . . . , 10}, we create a graph with two
terminal nodes and n − 2 non-terminal nodes. Each non-terminal node has a random
out-degree and distinct (also random) outgoing neighbours. We restrict the graphs to
exclude self-loops and duplicate edges, ensuring at least one terminal node is reachable
from any vertex. These SSG instances are meant to represent a randomly chosen
instance where no favourable decomposition is guaranteed.
Family B (separator-friendly SSGs). These graphs are constructed to specifically
favour the decomposition from Section [Link]. We create k SCCs, each of size s. These
SCCs are directed cycles. We then introduce one additional node, a separator, that has
one incoming edge from two nodes, from the last initialised node in that cycle, and
from the ⌊s/2⌋’th node to be initialised. To connect SCCs, we link them in a chain. The
first node of SCC z + 1 has an incoming edge from the separator of SCC z. Finally, we
add a “wrap-around” edge from the separator of SCC k to a node in SCC 1, creating a
“global” cycle across all SCCs.
We compare six algorithms: (1) simple iteration x(t+1) = f (x(t) ), (2) Algorithm 2, (3)
Algorithm from Dang et al. [DQY24], (4) Algorithm from Dang et al. with the heuristic
in Section 3.2.1, (5) recursive separator decomposition with 1D binary search as base
solver, and (6) recursive separator decomposition with 1D simple iteration as base
solver.
For the separator-based algorithms on random graphs, we precompute decompositions

46
Chapter 4. Implementation and Results 47

using a greedy heuristic: repeatedly identifying the largest SCC and removing the
vertex with the largest out-degree in that SCC. This is a practical but not optimal way to
obtain a separator. For separator-friendly graphs, separator nodes are given with the
construction. In both settings, the decomposition preprocessing time is excluded from
the reported runtimes below.
For every run, we record two metrics: query complexity (number of calls to f (x)) and
wall-clock runtime (seconds). All methods use tolerance ε = 2−n . In the current setup,
both setups are averaged over 100 repetitions per graph size, and reported values are
means over repetitions.
We limit the algorithm from Dang et al. to at most 5 vertices to keep the computation
feasible. In plots, we show the measured points and add extrapolated trend lines using a
log-linear fit.
Finally, we do not include Haslebacher et al. [HLSW25] in these experiments, because
its running time is impractical. As seen across all results, the algorithm from Dang

(a) Family A: Function evaluations. (b) Family A: Wall-clock runtime.

(c) Family B: Function evaluations. (d) Family B: Wall-clock runtime.

Figure 4.1: Performance comparison of fixed-point algorithms on random (Family A) and

separator-friendly (Family B) SSGs. Dashed lines indicate extrapolated trends.
Chapter 4. Implementation and Results 48

et al. [DQY24] scales poorly, requiring us to stop testing at n = 5. For random SSGs,
where no good decomposition is guaranteed, simple iteration proves to be the most
efficient method in terms of both time and query count. Decomposition with a simple
base case is also highly competitive, with near-identical query complexity.
However, our approach shows a clear advantage in separator-friendly SSGs. In these
cases, the decomposition (simple base) outperforms all other methods, requiring the
fewest function evaluations. We also observe that the heuristic helps the algorithm from
Dang et al. [DQY24] run much faster initially. Nevertheless, it still slows down rapidly,
which is expected given that its asymptotic performance is not improved.
Something that is interesting is how well simple iteration performs compared to the
other algorithms, and similarly, how the simple iteration base case outperforms the
binary search base case. This is because, unless a function is specifically crafted
in an adversarial manner, simple iteration is typically very fast, making it a strong
general-purpose method.
Finally, Algorithm 2, which can be viewed as simple iteration combined with a heuristic,
performs worse than simple iteration alone. A possible explanation is that the benefit
of quick search space reduction does not outweigh the increased number of function
evaluations required.
There are two caveats to this. First, the decompositions used in the decomposition
algorithm are precomputed; if this were not the case, the runtime results would likely
be very different, as such decompositions are themselves computationally expensive to
obtain. Secondly, for random graphs, the decomposition used is not optimal. As a result,
it is possible that, with a better decomposition, the decomposition-based algorithm
could outperform simple iteration in terms of query complexity.
 Chapter 5

Conclusion

This dissertation has explored the problem of computing fixed points of functions that
are both monotone and ℓ p -contractive, a setting that sits at the intersection of two
well-studied properties that have mostly been considered in isolation. While existence
results are well known, the efficient computation of such fixed points remains, to this
day, an open problem.
Simple algorithms exist, ranging from basic fixed-point iteration xn+1 = f (xn ) to gener-
alisations of binary search. However, in the worst case, these methods require a large
1
number of steps. In particular, most approaches have complexities of the form O 1−λ ,
O(cd ), or O 1ε . For the problems we are interested in, the contraction coefficient λ is

typically close to 1, ε is close to 0, and the dimension d is large. Hence, these methods
become intractable.
A recurring theme throughout this work is that, although there has been progress on
fixed points of ℓ p -contractions, there is still a large gap between query complexity and
truly time-efficient algorithms. Moreover, most existing results focus on the ℓ1 , ℓ∞ , and
ℓ2 norms, which is understandable given their applications, while the general monotone
ℓ p setting remains underdeveloped.
This dissertation makes several contributions towards addressing this gap.
We investigated the algorithm introduced by Haslebacher et al. [HLSW25] for comput-
ing fixed points of ℓ p -contractions. Their method achieves polynomial query complexity;
however, it is not time-efficient, as it requires an exponential-time search over a very fine
discrete grid in order to locate appropriate ℓ p -centre points. In their original formulation,
the ℓ1 case also involves iterating over an infinite set. We resolve this by showing that it
is sufficient to iterate over a large but finite set.
To address the bottleneck of iterating over the full grid, we followed one of their
suggested directions: using sampling. The idea is to sample from the region where the
fixed point may lie and compute sufficiently good query points from these samples. We
show that a polynomial number of samples is enough to approximate the true ℓ p -centre
point of the search space with high probability.

49
Chapter 5. Conclusion 50

However, this approach introduces new issues. In particular, it requires efficient uniform
sampling from the search space, which is non-trivial since the space is not necessarily
convex. We also consider whether the search space might be star-shaped, as this could
make known sampling techniques applicable, but we do not resolve this question. A
second issue is that there is currently no known efficient method for computing ℓ p -centre
points of n points in this general setting, which remains a bottleneck.
The second contribution is a decomposition approach based on the dependency graph
of f . Previous approaches to decomposition typically split coordinates “blindly”,
without using any information about f , which can lead to redundant computation of
subproblems. In contrast, our method uses the dependency graph of the function,
together with strongly connected components and balanced separators, to reduce the
number of times each subproblem needs to be computed.
This idea is inspired by decomposition used in the ARRIVAL problem [GHH21,
DGK+ 17], which we generalise to all ℓ p -contractions. In favourable cases, in par-
ticular when the dependency graph has small balanced separators or bounded treewidth,
this leads to improved performance. One limitation of our analysis is that, for simplicity,
we did not attempt to find tight bounds on the complexity, and it is likely that these
could be improved.
Empirically, the method performs particularly well on “separator-friendly” instances.
Another advantage is that the decomposition is reusable: once computed, the same
separators and SCC order can be used to find the fixed point of different functions that
share the same dependency graph. This becomes particularly useful in problems defined
on road networks, where not only does the dependency graph remain unchanged, but
the treewidth also tends to be small.
That said, the approach has limitations. In the worst case, the runtime remains exponen-
tial. Even for bounded-treewidth graphs, computing good separators is expensive, and
approximate separators can become ineffective when the treewidth is not sufficiently
small relative to the dimension. It would be interesting to understand whether better
graph parameters (instead of treewidth) or tighter analyses could improve this. Another
question is whether monotonicity can be used to improve this decomposition, using
ideas from Dang et al. [DQY24], for example by using the solutions to subproblems to
infer bounds on the separator.
The third contribution is a heuristic for monotone ℓ p -contractions. This heuristic
combines monotonicity and contraction to eliminate large parts of the search space using
only two evaluations of f per coordinate. The idea is simple, and a related statement
already holds for monotone (not necessarily contractive) functions: if f (xL )i > mi , then
there exists a fixed point with xi∗ > mi . Contraction strengthens this by guaranteeing
uniqueness, which makes the heuristic useful for locating where the fixed point lies.
Nevertheless, in many cases the conditions needed to determine where xi∗ lies are not
met, and we must fall back on an alternative algorithm.
In practice, this heuristic can significantly speed up existing algorithms, especially those
designed for monotone functions. For example, in SSGs it works particularly well when
there are many MAX and MIN nodes, since these tend to create more extreme values
Chapter 5. Conclusion 51

and reduce the number of “ambiguous” coordinates. However, as the experiments show,
the overhead of additional function evaluations can sometimes outweigh the gains.
The experimental results suggest that simple iterative methods often outperform more
sophisticated algorithms with stronger theoretical guarantees. This is mainly due to
large constant factors in the latter, which are designed for worst-case behaviour. At
the same time, our evaluation has several limitations: it is restricted to relatively small
instances and focuses entirely on SSGs with fixed distributions of node types (roughly
equal numbers of MIN, MAX, and random nodes). A more complete evaluation could
include other types of functions, such as instances arising from the ARRIVAL problem,
or vary the distribution of node types. It would also be interesting to construct explicit
worst-case instances for each algorithm, although this is not straightforward.
Looking forward, several directions stand out. One of the more interesting directions
towards a fully sampling-based algorithm would be to understand the geometric shape of
the search space in Haslebacher et al. [HLSW25], in particular whether it is star-shaped
or at least connected. Alternatively, new sampling methods specifically designed for
this space could be very useful.
There is also clear potential in combining the different approaches in this disserta-
tion. For example, integrating the heuristic with the decomposition framework could
improve practical performance. Exploring alternative graph parameters, such as di-
rected treewidth, may also lead to better decompositions. It would also be interesting
to identify further problems that can be reduced to finding fixed points of monotone
ℓ p -contractions, especially beyond the ℓ1 , ℓ∞ , and ℓ2 norm settings. This could help
further motivate the study of this class of problems.
In conclusion, this dissertation shows that, while there has been meaningful progress
in understanding fixed points of monotone ℓ p -contractive functions, many challenges
remain. At the same time, the methods explored here suggest several promising direc-
tions.
 Bibliography

[AB99] Martin Anthony and Peter Bartlett. Neural Network Learning: Theoreti-
cal Foundations. Cambridge University Press, 1999.
[AYBGM17] Yasin Abbasi-Yadkori, Peter Bartlett, Victor Gabillon, and Alan Malek.
Hit-and-Run for Sampling and Planning in Non-Convex Spaces. In Pro-
ceedings of the 20th International Conference on Artificial Intelligence
and Statistics, volume 54 of Proceedings of Machine Learning Research,
pages 888–895. PMLR, 2017.
[Ban22] Stefan Banach. Sur les opérations dans les ensembles abstraits et leur
application aux équations intégrales. Fundamenta Mathematicae, 3:133–
181, 1922.
[Bau58] Heinz Bauer. Minimalstellen von Funktionen und Extremalpunkte.
Archiv der Mathematik, 9(4):389–393, 1958.
[BDH+ 12] Dietmar Berwanger, Anuj Dawar, Paul Hunter, Stephan Kreutzer, and Jan
Obdržálek. The dag-width of directed graphs. Journal of Combinatorial
Theory, Series B, 102(4):900–923, 2012.
[BFG+ 25] Eleni Batziou, John Fearnley, Spencer Gordon, Ruta Mehta, and Rahul
Savani. Monotone Contractions. In Proceedings of the 57th Annual
ACM Symposium on Theory of Computing, STOC ’25, pages 507–517.
Association for Computing Machinery, 2025.
[Bir40] Garrett Birkhoff. Lattice Theory. Number v. 25, pt. 2 in American
Mathematical Society colloquium publications. American Mathematical
Society, 1940.
[Bro11] L. E. J. Brouwer. Über Abbildung von Mannigfaltigkeiten. Mathematis-
che Annalen, 71(1):97–115, 1911.
[CDV10] Karthekeyan Chandrasekaran, Daniel Dadush, and Santosh Vempala.
Thin Partitions: Isoperimetric Inequalities and a Sampling Algorithm for
Star Shaped Bodies, pages 1630–1645. SIAM, 2010.
[CFK+ 15] Marek Cygan, Fedor V. Fomin, Łukasz Kowalik, Daniel Lokshtanov,
Dániel Marx, Marcin Pilipczuk, Michał Pilipczuk, and Saket Saurabh.
Parameterized Algorithms. Springer International Publishing, Cham,
2015.

52
Bibliography 53

[Che24] Yeshwanth Cherapanamjeri. Computing Approximate Centerpoints in

Polynomial Time. In 2024 IEEE 65th Annual Symposium on Foundations
of Computer Science (FOCS), pages 1654–1668, 2024.
[CL55] Earl A. Coddington and Norman Levinson. Theory of Ordinary Differ-
ential Equations. International series in pure and applied mathematics.
McGraw-Hill, 1955.
[CL22] Xi Chen and Yuhao Li. Improved Upper Bounds for Finding Tarski Fixed
Points. In Proceedings of the 23rd ACM Conference on Economics and
Computation, EC ’22, pages 1108–1118, New York, NY, USA, 2022.
Association for Computing Machinery.
[CLRS22] Thomas H. Cormen, Charles Eric Leiserson, Ronald L. Rivest, and
Clifford Stein. Introduction to Algorithms. The MIT Press, Cambridge,
Massachusetts, fourth edition, 2022.
[CLT08] Ching-Lueh Chang, Yuh-Dauh Lyuu, and Yen-Wu Ti. The complexity of
Tarski’s fixed point theorem. Theoretical Computer Science, 401(1):228–
235, 2008.
[CLY25] Xi Chen, Yuhao Li, and Mihalis Yannakakis. Computing a Fixed Point of
Contraction Maps in Polynomial Queries. Journal of the ACM, 72(4):1–
20, July 2025.
[CLY26] Xi Chen, Yuhao Li, and Mihalis Yannakakis. Quadratic Speedup for
Computing Contraction Fixed Points, 2026.
[Con92] Anne Condon. The complexity of stochastic games. Information and
Computation, 96(2):203–224, 1992.
[DFK91] Martin Dyer, Alan Frieze, and Ravi Kannan. A random polynomial-time
algorithm for approximating the volume of convex bodies. Journal of
the ACM, 38(1):1–17, 1991.
[DGK+ 17] Jérôme Dohrau, Bernd Gärtner, Manuel Kohler, Jiří Matoušek, and Emo
Welzl. ARRIVAL: A Zero-Player Graph Game in NP ∩ coNP, pages
367–374. Springer International Publishing, Cham, 2017.
[DM41] Ben Dushnik and E. W. Miller. Partially Ordered Sets. American Journal
of Mathematics, 63(3):600–610, 1941.
[DQY24] Chuangyin Dang, Qi Qi, and Yinyu Ye. Computations and Complexities
of Tarski’s Fixed Points and Supermodular Games. In Frontiers of
Algorithmics: 18th International Joint Conference, IJTCS-FAW 2024,
pages 159–174, Berlin, Heidelberg, 2024. Springer-Verlag.
[EPRY20] Kousha Etessami, Christos Papadimitriou, Aviad Rubinstein, and Mihalis
Yannakakis. Tarski’s Theorem, Supermodular Games, and the Complex-
ity of Equilibria. In 11th Innovations in Theoretical Computer Science
Conference (ITCS 2020), volume 151 of Leibniz International Proceed-
Bibliography 54

ings in Informatics (LIPIcs), pages 18:1–18:19, Dagstuhl, Germany,

2020. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.
[FGMS20] John Fearnley, Spencer Gordon, Ruta Mehta, and Rahul Savani. Unique
end of potential line. Journal of Computer and System Sciences, 114:1–
35, 2020.
[FPS22] John Fearnley, Dömötör Pálvölgyi, and Rahul Savani. A Faster Algorithm
for Finding Tarski Fixed Points. ACM Trans. Algorithms, 18(3), October
2022.
[GHH21] Bernd Gärtner, Sebastian Haslebacher, and Hung P. Hoang. A Subex-
ponential Algorithm for ARRIVAL. In 48th International Colloquium
on Automata, Languages, and Programming (ICALP 2021), volume
198 of Leibniz International Proceedings in Informatics (LIPIcs), pages
69:1–69:14. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2021.
[GJ95] Paul W. Goldberg and Mark R. Jerrum. Bounding the Vapnik–
Chervonenkis dimension of concept classes parameterized by real num-
bers. Machine Learning, 18(2):131–148, 1995.
[Has25] Sebastian Haslebacher. ARRIVAL: Recursive framework & ℓ1 contrac-
tion. In 52nd International Colloquium on Automata, Languages, and
Programming (ICALP 2025), volume 334 of Leibniz International Pro-
ceedings in Informatics (LIPIcs), pages 95:1–95:17. Schloss Dagstuhl –
Leibniz-Zentrum für Informatik, 2025.
[HHMM20] G. Hansen, I. Herburt, H. Martini, and M. Moszyńska. Starshaped sets.
Aequationes Mathematicae, 94(6):1001–1092, 2020. Page of interest:
1032.
[HLSW25] Sebastian Haslebacher, Jonas Lill, Patrick Schnider, and Simon Weber.
Query-Efficient Fixpoints of ℓ p -Contractions. In 2025 IEEE 66th Annual
Symposium on Foundations of Computer Science (FOCS), pages 2037–
2053, 2025.
[JRST01] Thor Johnson, Neil Robertson, P. D. Seymour, and Robin Thomas. Di-
rected tree-width. Journal of Combinatorial Theory, Series B, 82(1):138–
154, 2001.
[Kar17] C. S. Karthik. Did the train reach its destination: The complexity of
finding a witness. Information Processing Letters, 121:17–21, 2017.
[Kor23] Tuukka Korhonen. A Single-Exponential Time 2-Approximation Algo-
rithm for Treewidth. SIAM Journal on Computing, pages FOCS21–174–
FOCS21–194, November 2023.
[Lin94] Ernest Lindelöf. Sur l’application des méthodes d’approximations suc-
cessives à l’étude des intégrales réelles des équations différentielles
ordinaires. Journal de Mathématiques Pures et Appliquées, 10:117–128,
1894.
Bibliography 55

[Rad38] R. Rado. A Theorem on General Measure Functions. Proceedings of the

London Mathematical Society, s2-44(1):61–91, 1938.
[RS86] Neil Robertson and P. D. Seymour. Graph minors. II. Algorithmic aspects
of tree-width. Journal of Algorithms, 7(3):309–322, 1986.
[Sha53] L. S. Shapley. Stochastic Games. Proceedings of the National Academy
of Sciences, 39(10):1095–1100, 1953.
[SS02] Spencer Shellman and K. Sikorski. A Two-Dimensional Bisection Enve-
lope Algorithm for Fixed Points. Journal of Complexity, 18(2):641–659,
2002.
[SS03] Spencer Shellman and K. Sikorski. A recursive algorithm for the infinity-
norm fixed point problem. Journal of Complexity, 19(6):799–834, 2003.
[STW93] K. Sikorski, C. W. Tsay, and H. Woźniakowski. An Ellipsoid Algorithm
for the Computation of Fixed Points. Journal of Complexity, 9(1):181–
200, 1993.
[Tal94] Michael Talagrand. Sharper Bounds for Gaussian and Empirical Pro-
cesses. The Annals of Probability, 22(1):28–76, 1994.
[Tar55] Alfred Tarski. A Lattice-Theoretical Fixpoint Theorem and Its Applica-
tions. Pacific Journal of Mathematics, 5(2):285–309, 1955.
[Tar72] Robert Tarjan. Depth-First Search and Linear Graph Algorithms. SIAM
Journal on Computing, 1(2):146–160, 1972.
[Ten92] Shang-Hua Teng. Points, spheres, and separators: a unified geometric
approach to graph partitioning. PhD thesis, Carnegie Mellon University,
USA, 1992. UMI Order No. GAX92-20224.
[Tod08] Ioan Todinca. Exact Algorithms for Treewidth, pages 1–4. Springer US,
Boston, MA, 2008.
[VC15] V. N. Vapnik and A. Ya. Chervonenkis. On the Uniform Convergence
of Relative Frequencies of Events to Their Probabilities, pages 11–30.
Springer International Publishing, Cham, 2015.
 Appendix A

Figures and Examples

(a) An instance of the ARRIVAL problem. The (b) A Condon Simple Stochastic Game (SSG).
token starts at o and moves through (A, B,C) Vertices are MAX (v1 ), MIN (v2 , v3 ), or random
until reaching t or t . Solid and dashed edges (v4 ). Double circles denote the 1- and 0-sinks;
denote the “even” and “odd” successors. random transitions occur with probability 1/2.

Figure A.1: Two motivating examples reducible to fixed-point computation.

56
Appendix A. Figures and Examples 57

Figure A.2: An example where the fixed point x∗ is not a star centre of a bisector ℓ1 -
halfspace. Here, the query point is c = (4, 0), the fixed point is x∗ = (4/3, 4/3) and the
ℓ1 contraction is f (x) = −0.5(x − x∗ ) + x∗ .

Figure A.3: The graph is divided into a separator set S and components C1 and C2 .
Removing the separator S leaves a sub-graph that forms a DAG over SCCs.

(a) The original undirected graph G. (b) A valid tree decomposition of G.

Figure A.4: An example of a graph and its corresponding tree decomposition. Since no
decomposition with smaller maximum bag size exists, the treewidth of G is 3 − 1 = 2.
Appendix A. Figures and Examples 58

Example A.0.1. Let f : R6 → R6 be defined as:

0.1x1 + sin2 (x6 ) + cos2 (x6 ) − 1

 

 0.2x1 + 0.2x3 

 0.4x2 
f (x) = 
 . (A.1)
 0.1x 1 + 0.3x 3


 0.1x2 + 0.1x4 + 0.2x6 
0.5x5

Based on the definition of the dependency graph G = (V, E), the edge set is:

E = {(2, 1), (2, 3), (3, 2), (4, 1), (4, 3), (5, 2), (5, 4), (5, 6), (6, 5)}. (A.2)

Note that while the symbolic expression for f1 contains x6 , the identity sin2 (x6 ) +
cos2 (x6 ) − 1 = 0 ensures f1 (x) = 0.1x1 globally. Because f1 (x) is invariant to changes
in x6 , there exist no x, y ∈ R6 differing only at index 6 such that f1 (x) ̸= f1 (y). Thus,
the edge (1, 6) ∈/ E.
Example A.0.2. Let d = 2 and define the box B = [−1, 1]2 , where L = (−1, −1) and
U = (1, 1).
We will construct two functions f : B → B that satisfy conditions 1 and 2. However,
both will fail to isolate a fixed point.
Case 1: Monotone, but not contractive
Define the map f : B → B as: f (x1 , x2 ) = (x1 , x2 ) This function is monotone since
x ≤ y =⇒ f (x) ≤ f (y). However, it is not a contraction, for any ℓ p -norm, its contrac-
tion coefficient is 1 ∈
/ [0, 1). Because every point in B is a fixed point, applying the
proposition causes the bounds to “track” different points.
Consider bisection steps on the first coordinate. First, test m1 = 0.5 with the lower test
vector:
f1 (x(L) ) = f1 (0.5, −1) = 0.5 ≥ 0.5
Condition 1 holds, so an algorithm would update the lower bound to L1 = 0.5. Now,
test m1 = −0.5 using the upper test vector:

f1 (x(U) ) = f1 (−0.5, 1) = −0.5 ≤ −0.5

Condition 2 holds, so the algorithm updates the upper bound to U1 = −0.5. The
algorithm now concludes the fixed point x∗ must satisfy 0.5 ≤ x1∗ ≤ −0.5. The search
bounds have crossed (L1 > U1 ), resulting in an empty search region.
Case 2: Contractive, but not monotone
Fix the index i = 1 and choose m1 = 0.5. The test vectors are x(L) = (0.5, −1) and
x(U) = (0.5, 1). Define the map f : B → B as:

f (x1 , x2 ) = (0.6|x2 |, 0)

This function is a contraction with respect to any ℓ p -norm with coefficient λ = 0.6 < 1,
because ∥ f (x) − f (y)∥ p = 0.6 |x2 | − |y2 | ≤ 0.6|x2 − y2 | ≤ 0.6∥x − y∥ p . However, it is
Appendix A. Figures and Examples 59

not monotone: for x = (0, −1) and y = (0, 0), we have x ≤ y, but f1 (x) = 0.6 ̸≤ 0 = f1 (y).
The unique fixed point is x∗ = (0, 0). Evaluating f1 at the test vectors gives:

f1 (x(L) ) = 0.6| − 1| = 0.6 ≥ 0.5

f1 (x(U) ) = 0.6|1| = 0.6 ≥ 0.5

The premise for both conditions holds, but the conclusion x1∗ ≥ 0.5 fails for the unique
fixed point.