Faculty of Electrical Engineering, Mathematics, and Computer Science,
University of Twente, The Netherlands and
https://people.utwente.nl/b.mantheyb.manthey@utwente.nlhttps://orcid.org/0000-0001-6278-5059
Faculty of Electrical Engineering, Mathematics, and Computer Science,
University of Twente, The Netherlands and
https://people.utwente.nl/j.vanrhijnj.vanrhijn@utwente.nlhttps://orcid.org/0000-0002-3416-7672Supported
by NWO grant OCENW.KLEIN.176. Bodo Manthey and Jesse van Rhijn [100]Theory of
computation Randomness, geometry and discrete structures [100]Theory of
computation Approximation algorithms analysis [100]Theory of computation
Discrete optimization [#1][Restated]

# Worst-Case and Smoothed Analysis of Hartigan’s Method for k-Means Clustering

Bodo Manthey  Jesse van Rhijn

###### Abstract

We analyze the running time of Hartigan’s method, an old algorithm for the
$k$-means clustering problem. First, we construct an instance on the line on
which the method can take $2^{\Omega(n)}$ steps to converge, demonstrating
that Hartigan’s method has exponential worst-case running time even when
$k$-means is easy to solve. As this is in contrast to the empirical
performance of the algorithm, we also analyze the running time in the
framework of smoothed analysis. In particular, given an instance of $n$ points
in $d$ dimensions, we prove that the expected number of iterations needed for
Hartigan’s method to terminate is bounded by
$k^{12kd}\cdot\operatorname{poly}(n,k,d,1/\sigma)$ when the points in the
instance are perturbed by independent $d$-dimensional Gaussian random
variables of mean $0$ and standard deviation $\sigma$.

###### keywords:

k-means clustering, smoothed analysis, probabilistic analysis, local search,
heuristics

###### category:

##  1 Introduction

Clustering is an important problem in computer science, from both a practical
and a theoretical perspective. On the practical side, identifying clusters of
similar points in large data sets has relevance to fields ranging from physics
to biology to sociology. Recent advances in machine learning and big data have
made the need for efficient clustering algorithms even more apparent. On the
theoretical side, clustering problems continue to be a topic of research from
the perspective of approximation algorithms, heuristics, and computational
geometry.

Perhaps the best-studied clustering problem is that of $k$-means clustering.
In this problem, one is given a finite set of points
$\mathcal{X}\subseteq\mathbb{R}^{d}$ and an integer $k$. The goal is to
partition the points into $k$ subsets, such that the sum of squared distances
of each point to the centroid of its assigned cluster, also called its cluster
center, is minimized.

Despite great effort to devise approximation algorithms for $k$-means
clustering, the method of choice remains Lloyd’s method [11]. This method
starts with an arbitrary choice of centers, and assigns each point to its
closest center. The centers are then moved to the centroids of each cluster.
On the next iteration, each point is again reassigned to its closest center,
and the process repeats.

It is not hard to show that this process strictly decreases the objective
function whenever either a cluster center changes position, or a point is
reassigned. Hence, no clustering can show up twice during an execution of this
algorithm. Since the number of partitions of $n$ points into $k$ sets is at
most $k^{n}$, the process must eventually terminate.

Although Lloyd’s method has poor approximation performance both in theory and
in practice [2], its speed has kept it relevant to practitioners. This is in
startling contrast to its worst-case running time, which is exponential in the
number of points [15].

To close the gap between theory and practice, Arthur et al. have shown that
Lloyd’s method terminates in expected polynomial time on perturbed point sets,
by means of a smoothed analysis [1]. This provides some theoretical
justification for the use of Lloyd’s method in practice.

Another, less well-known heuristic for clustering is Hartigan’s method [8]. In
this method, one proceeds point-by-point. Given an arbitrary clustering, one
checks whether there exists a point that can be reassigned to a different
cluster, such that the objective function decreases. If such a point exists,
it is reassigned to this new cluster. If no such points exist, the algorithm
terminates and the clustering is declared locally optimal.

Although at first sight Hartigan’s method might seem like a simpler version of
Lloyd’s method, it is qualitatively different. If Lloyd’s method reassigns a
point $x$ from cluster $i$ to cluster $j$, then $x$ must be closer to the
center of cluster $j$ than to that of cluster $i$. In Hartigan’s method, this
is not true; $x$ may be reassigned even when there are no cluster centers
closer to $x$ than its current center. This can be beneficial, as Telgarsky &
Vattani showed that Hartigan’s method is more powerful than Lloyd’s method
[14].

To be precise, every local optimum of Hartigan’s method is also a local
optimum of LLoyd’s method, while the converse does not hold. Telgarsky &
Vattani moreover performed computational experiments, which show that
Hartigan’s method not only tends to find better clusterings than Lloyd’s, but
also has a similar running time on practical instances. Despite these
promising results, theoretical knowledge of Hartigan’s method is lacking.

In this paper, we aim to advance our understanding of this heuristic. Our
contributions are twofold. First, we construct an instance on the line on
which Hartigan’s method can take $2^{\Omega(n)}$ iterations to terminate.
Considering that $k$-means clustering can be solved exactly in polynomial time
in $d=1$, this shows that the worst-case running time of Hartigan’s method is
very poor even on easy instances. This is in contrast to Lloyd’s method, where
all known non-trivial lower bounds require $d\geq 2$.

######  Theorem 1.1 (restate=lowerbound).

For each $m\in\mathbb{N}_{\geq 2}$ there exists an instance of $k$-means
clustering on the line with $n=4m-3$ points and $k=2m-1$ clusters on which
Hartigan’s method can take $2^{\Omega(n)}$ iterations to converge to a local
optimum.

Second, we attempt to reconcile Theorem 1.1 with the observed practical
performance of Hartigan’s method. We perform a smoothed analysis of its
running time, in which each point in an arbitrary instance is independently
perturbed by a Gaussian random variable of variance $\sigma^{2}$.

######  Theorem 1.2 (restate=smoothedcomplexity).

Let $n,k,d\in\mathbb{N}$, and assume $4kd\leq n$. Fix a set of $n$ points
$\mathcal{Y}\subseteq[0,1]^{d}$, and assume that each point in $\mathcal{Y}$
is independently perturbed by a $d$-dimensional Gaussian random variable with
mean $0$ and standard deviation $\sigma$, yielding a new set of points
$\mathcal{X}$. Then the expected running time of Hartigan’s method on
$\mathcal{X}$ is bounded by

$O\mathopen{}\mathclose{{}\left(\frac{k^{12kd+5}d^{11}n^{12.5+\frac{1}{d}}\ln^{4.5}(nkd)}{\sigma^{8}}}\right)=k^{12kd}\cdot\operatorname{poly}(n,k,d,1/\sigma).$

This is a first step to settling a conjecture by Telgarsky & Vattani that
Hartigan’s method, like Lloyd’s, should have polynomial smoothed running time.

##  2 Preliminaries and Notation

Given vectors $x,y\in\mathbb{R}^{d}$, we write $\langle x,y\rangle$ for the
standard Euclidean inner product on $\mathbb{R}^{d}$, and $\|x\|=\sqrt{\langle
x,x\rangle}$ for the standard norm.

Given a set of $k$ clusters
$\mathcal{C}=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{k}\\}$, a configuration of
a cluster $\mathcal{C}_{i}\in\mathcal{C}$ is an assignment of a set of points
to $\mathcal{C}_{i}$. We will denote the clusters by calligraphic letters, and
their configurations by regular letters; i.e., the configuration of
$\mathcal{C}_{i}$ will be denoted $C_{i}$. This distinction is sometimes
useful. For the majority of this paper, however, we will not make this
distinction explicitly, and will refer to both a cluster and its configuration
interchangeably by regular letters.

Given a finite set of points $S\subseteq\mathbb{R}^{d}$, we define the center
of mass of $S$ as

$\operatorname{cm}(S)=\frac{1}{|S|}\sum_{x\in S}x.$

With this definition, we can formally define the objective function of
$k$-means. Let $C=\\{C_{i}\\}_{i=1}^{k}$ be a partition of a finite set of
points $\mathcal{X}\subseteq\mathbb{R}^{d}$. Then the objective function of
$k$-means is

$\Phi(C)=\sum_{i=1}^{k}\sum_{x\in
C_{i}}\|x-\operatorname{cm}(C_{i})\|^{2}=\sum_{i=1}^{k}\Phi(C_{i}),$

where we define $\Phi(C_{i})=\sum_{x\in
C_{i}}\|x-\operatorname{cm}(C_{i})\|^{2}$. We will also refer to $\Phi(C)$ as
the potential function.

For both the worst-case and smoothed complexity bounds, be need to analyze the
improvement of a single iteration. Thus, we need a simple expression for this
quantity. Lemmas 2.1 and 2.2 allow us to obtain such an expression. These
results were already obtained by Telgarsky & Vattani [14].

######  Lemma 2.1.

Let $S$ and $T$ be two disjoint nonempty sets of points in $\mathbb{R}^{d}$.
Then

$\Phi(S\cup
T)-\Phi(S)-\Phi(T)=\frac{|S|\cdot|T|}{|S|+|T|}\cdot\|\operatorname{cm}(S)-\operatorname{cm}(T)\|^{2}.$

######  Lemma 2.2.

Let $S$ and $T$ be two disjoint nonempty sets of points in $\mathbb{R}^{d}$
with $|S|>1$. Suppose we move a point $x\in S$ from $S$ to $T$. Then

$\Phi(S\setminus\\{x\\})+\Phi(T\cup\\{x\\})-\Phi(T)-\Phi(S)=\frac{|T|}{|T|+1}\|\operatorname{cm}(T)-x\|^{2}-\frac{|S|}{|S|-1}\|\operatorname{cm}(S)-x\|^{2}.$

Let $C$ be some clustering of $\mathcal{X}$. Suppose in some iteration of
Hartigan’s, we move $x\in C_{i}$ to $C_{j}$. Let the gain of this iteration be
denoted $\Delta_{x}(C_{i},C_{j})$. Then Lemma 2.2 tells us that

$\Delta_{x}(C_{i},C_{j})=\frac{|C_{i}|}{|C_{i}|-1}\|x-\operatorname{cm}(C_{i})\|^{2}-\frac{|C_{j}|}{|C_{j}|+1}\|x-\operatorname{cm}(C_{j})\|^{2}.$

At first sight, it seems like Lemma 2.2 leaves open the possibility that a
cluster is left empty. The following lemma shows that this can never happen.

######  Lemma 2.3.

No iteration can leave a cluster empty.

######  Proof 2.4.

Suppose before an iteration, $C_{i}=\\{x\\}$ for some $x\in X$, and after the
iteration $C_{i}^{\prime}=\emptyset$ and $C_{j}^{\prime}=C_{j}\cup\\{x\\}$,
i.e. $x$ is moved from cluster $i$ to cluster $j$. The gain of this iteration
is then (Lemma 2.1)

$\Phi(C_{i})+\Phi(C_{j})-\Phi(\emptyset)-\Phi(C_{j}\cup\\{x\\})=\Phi(C_{j})-\Phi(C_{j}\cup\\{x\\})=-\frac{|C_{j}|}{|C_{j}|+1}\|x-\operatorname{cm}(C_{j})\|^{2}\leq
0,$

since $\operatorname{cm}(C_{i})=x$ and $\Phi(\emptyset)=0$. Since every
iteration must improve the clustering, this concludes the proof.

##  3 Exponential Lower Bound

In this section, we construct a family of $k$-means instances on the line on
which Hartigan’s method can take an exponential number of iterations before
reaching a local optimum. To be precise, we prove the following theorem.

*

The construction we employ is similar to the construction used by Vattani for
Lloyd’s method [15]. However, Hartigan’s method only reassigns a single point
in each iteration, and we are free to choose which point we reassign.
Moreover, we are even free to choose which cluster we move a point to if there
are multiple options. This allows us to simplify the construction and embed it
in a single dimension, rather than the plane used by Vattani.

We define a set of $m$ gadgets $G_{i}$, $i\in\\{0,\ldots,m-1\\}$. Each gadget
except for the “leaf” gadget $G_{0}$ consists of four points, and has two
clusters $G_{i}(\mathcal{C}_{0})$ and $G_{i}(\mathcal{C}_{1})$ associated with
it. Moreover, each gadget except $G_{0}$ has three distinguished states,
called “morning”, “afternoon”, and “asleep”. The leaf gadget only has two
states, “awake” and “asleep ”.

During the morning state, a gadget $G_{i}$ watches $G_{i-1}$. If $G_{i-1}$
falls asleep, then it is awoken by $G_{i}$; this is achieved by moving a point
of $G_{i}$ to one of the clusters of $G_{i-1}$. This allows $G_{i-1}$ to
perform a sequence of iterations, which ends with $G_{i-1}$ back in its
morning state.

Meanwhile, $G_{i}$ performs a sequence of iterations that transition it to its
afternoon state. During the afternoon state, it once more watches $G_{i-1}$.
When the latter falls asleep, $G_{i}$ once again wakes $G_{i-1}$, and
transitions itself to its asleep state.

The leaf gadget $G_{0}$, as it does not watch any gadgets, only ever awakens
and immediately falls asleep again.

We end the sequence of iterations once gadget $m-1$ falls asleep. Observe that
with this construction, $G_{i}$ falls asleep twice as often as $G_{i+1}$. With
the condition that $G_{m-1}$ falls asleep once, we obtain a sequence of at
least $2^{m-1}$ iterations. With $n=4m-3$, this yields Theorem 1.1.

For space reasons, we only describe the instance and the exponential-length
sequence here. The proof that this sequence is improving, which completes the
proof of Theorem 1.1, is deferred to the full version.

###  3.1 Formal Construction

We now give a detailed construction of a unit gadget, $G$. All gadgets except
for $G_{0}$ are scaled and translated versions of $G$. The unit gadget is a
tuple $G=(S,\mathcal{C}_{0},\mathcal{C}_{1})$, where
$S=\\{a,b,p,q\\}\subseteq\mathbb{R}$, and $\mathcal{C}_{0}$ and
$\mathcal{C}_{1}$ are two clusters. The positions of the points in $S$ are
given in Table 1. In addition, the gadget is depicted schematically in Figures
1 and 2. Note that the relative positions of the points in these figures do
not correspond to Table 1, but are chosen for visual clarity.

Point | $a$ | $b$ | $p$ | $q$ | $f$ | $t_{0}$
---|---|---|---|---|---|---
Position | 9 | 6 | 5 | 13 | 0 | 8
Table 1: Positions of the points in $S(G)$, the leaf point $f$, and the
translation vector $t_{0}$ between gadgets $G_{1}$ and $G_{2}$.

We remark that the points in Table 1 are not simply chosen by trial-and-error.
As will be explained shortly, we can obtain from our construction a series of
inequalities that must be satisfied by the points in $S$. We then obtained
these points by solving the model

$\displaystyle\min$ $\displaystyle\quad
a^{2}+b^{2}+p^{2}+q^{2}+f^{2}+t_{0}^{2}$ s.t. each move decreases the
clustering cost, $\displaystyle\quad a,b,p,q,f,t_{0}\in\mathbb{Z}$

using Gurobi [7]. The first constraint here amounts to satisfying a series of
inequalities of the form $\Delta_{x}(A,B)>0$ for $x\in S(G)$ and $A,B$ subsets
of the points in a gadget and its neighboring gadgets. For space reasons, we
defer their derivation and verification to the full version. The objective
function here is purely chosen so that Gurobi prefers to choose small integers
in the solution.

To construct $G_{i}$ from the unit gadget (for $i\geq 1)$, we scale the unit
gadget by a factor $5^{i-1}$, and translate it by
$t_{i}=\sum_{j=0}^{i-1}5^{j}t_{0}$, where $t_{0}=8$. Since each gadget only
ever exchanges points with its neighbors in the sequence we are about to
construct, it will suffice in proving Theorem 1.1 to consider only iterations
involving $G_{i}$, $G_{i-1}$ and $G_{i+1}$ for some fixed $i>2$. For the leaf
gadget, we simply have $G_{0}=(S_{0},\mathcal{C}_{0})$, where
$S_{0}=\\{f\\}=\\{0\\}$.

Before we go on to construct an improving sequence of exponential length, we
define the earlier-mentioned states. For ease of notation, we will refer to
the points of $G_{i}$ as $a_{i}$, $b_{i}$, and so on, and to the clusters of
$G_{i}$ as $\mathcal{C}_{0}(G_{i})$ and $\mathcal{C}_{1}(G_{i})$. Then we say
the state of $G_{i>0}$ is:

  * •

asleep, if $C_{0}(G_{i})=\\{b_{i}\\}$ and $C_{1}(G_{i})=\\{a_{i},q_{i}\\}$ (in
this state, $p_{i}$ is in some cluster of $G_{i-1}$);

  * •

morning, if $C_{0}(G_{i})=\\{p_{i},q_{i},b_{i}\\}$ and
$C_{1}(G_{i})=\\{a_{i}\\}$;

  * •

afternoon, if $C_{0}(G_{i})=\\{b_{i}\\}$ and
$C_{1}(G_{i})=\\{p_{i},q_{i},a_{i}\\}$.

For the leaf gadget, we say its state is:

  * •

asleep, if $C_{0}(G_{0})=\\{f\\}$;

  * •

awake, otherwise.

We now explicitly determine a sequence of iterations of exponential length. In
the proof of Theorem 1.1, we show that this sequence is improving. To analyze
the sequence, we consider the perspective of $G_{i}$ as it wakes up $G_{i-1}$
and falls asleep; and then as it is awoken by $G_{i+1}$. We first consider
only the case that $G_{i-1}\neq G_{0}$. See Figure 1 and Figure 2 for a
schematic depiction of the sequence described below.

Morning$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$

Waking up $G_{i-1}$ (1)$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$

Waking up $G_{i-1}$ (2)$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$

Afternoon$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$

Waking up $G_{i-1}$ / asleep$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$

Figure 1: Schematic depiction of the interactions between $G_{i}$ and
$G_{i-1}$ during the morning and afternoon phases of $G_{i}$.

#### Morning.

We start with $G_{i}$ in the morning state, and $G_{i-1}$ asleep. To wake up
$G_{i-1}$, the point $p_{i}$ moves to $\mathcal{C}_{1}(G_{i-1})$, which
currently contains $a_{i-1}$ and $q_{i-1}$. This triggers the wakeup phase of
$G_{i-1}$; we will analyze this phase later from the perspective of $G_{i}$.
When the wakeup phase completes, $\mathcal{C}_{1}(G_{i-1})$ contains $a_{i-1}$
and $p_{i}$, and $p_{i}$ moves to $\mathcal{C}_{1}(G_{i})$. Subsequently
$q_{i}$ moves from $\mathcal{C}_{0}(G_{i})$ to $\mathcal{C}_{1}(G_{i})$.
Observe that this puts $G_{i}$ into the afternoon state.

#### Afternoon.

In this state, $G_{i}$ is once again watching $G_{i-1}$. Once the latter falls
asleep, $p_{i}$ moves from $\mathcal{C}_{1}(G_{i})$ to
$\mathcal{C}_{1}(G_{i-1})$, which triggers another wakeup phase of $G_{i-1}$.
Additionally, this move causes $G_{i}$ to fall asleep. Thus, at the end of the
wakeup phase of $G_{i-1}$, we have $G_{i+1}$ wake up $G_{i}$.

#### Waking up.

First, the point $p_{i+1}$ joins $\mathcal{C}_{1}(G_{i})$. Next, $p_{i}$ moves
from $\mathcal{C}_{1}(G_{i-1})$ to $\mathcal{C}_{0}(G_{i})$. Then, $q_{i}$
moves from $\mathcal{C}_{1}(G_{i})$ to $\mathcal{C}_{0}(G_{i})$, and finally,
$p_{i+1}$ leaves $\mathcal{C}_{1}(G_{i})$, and joins either
$\mathcal{C}_{1}(G_{i+1})$ (if $G_{i+1}$ was in the morning state when waking
up $G_{i}$) or $\mathcal{C}_{0}(G_{i+1})$ (if $G_{i+1}$ was in the afternoon
state; in this case, the move of $p_{i+1}$ occurs during the wakeup phase of
$G_{i+1}$).

#### Leaf gadget.

The leaf gadget does not watch or wake up any other gadgets. It only wakes up
when $p_{1}$ moves into $\mathcal{C}_{0}(G_{0})$, and falls asleep again when
$p_{1}$ moves back to a cluster of $G_{1}$.

#### Initialization.

The sequence starts with all gadgets in the morning state, except for $G_{0}$,
which is asleep.

At every step, we have the gadget with the smallest index that is not asleep
wake up the gadget that it is watching. From this sequence of iterations, we
can retrieve a series of inequalities, each of which encodes the condition
that the gain of every iteration must be positive. To prove Theorem 1.1, we
must show that the points in Table 1 satisfy these inequalities.

An implementation of the sequence described above is provided in the following
link: https://pastebin.com/raw/RhZyut4X.

Asleep$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$$p_{i+1}$

$G_{i+1}$ wakes up $G_{i}$
(1)$a_{i}$$b_{i}$$p_{i}$$q_{i}$$a_{i-1}$$q_{i-1}$$p_{i+1}$

$G_{i+1}$ wakes up $G_{i}$ (2)$a_{i}$$b_{i}$$p_{i}$$q_{i}$$p_{i+1}$

$G_{i+1}$ wakes up $G_{i}$ (3)$a_{i}$$b_{i}$$p_{i}$$q_{i}$$p_{i+1}$

Morning$a_{i}$$b_{i}$$p_{i}$$q_{i}$$p_{i+1}$

Figure 2: Schematic depiction of the interactions between $G_{i}$, $G_{i-1}$
and $G_{i+1}$ during the wakeup phase of $G_{i}$. Note that the final state of
$G_{i}$ corresponds to the first state depicted in Figure 1.

##  4 Smoothed Analysis

For a smoothed analysis, the first hope might be to straightforwardly adapt a
smoothed analysis of Lloyd’s algorithm, e.g. that of Arthur, Manthey and
Röglin [1]. On closer inspection, however, such analyses strongly rely on a
couple of properties of Lloyd’s method that are not valid in Hartigan’s.

First, in Lloyd’s algorithm the hyperplane that bisects two cluster centers
also separates their corresponding clusters, since every point is always
assigned to the cluster center closest to itself. Second, the two stages of
Lloyd’s algorithm, moving the cluster centers and reassigning points, both
decrease the potential. Neither of these properties are satisfied by
iterations of Hartigan’s method. Hence, any analysis that relies on either
property cannot be easily repurposed.

Instead, we will use a different technique, more closely related to the
analysis of the Flip heuristic for Max-Cut with squared Euclidean distances by
Etscheid and Röglin [6]. The main result we will work towards in this section
is stated in Theorem 1.2.

###  4.1 Technical Preliminaries

Let $\mathcal{Y}\subseteq[0,1]^{d}$ be a set of $n$ points. Throughout the
remainder, we will denote by $\mathcal{X}$ the set of points obtained by
perturbing each point in $\mathcal{Y}$ independently by a $d$-dimensional
Gaussian vector of mean $0$ and standard deviation $\sigma\leq 1$. Note that
this last assumption is not actually a restriction. If $\sigma>1$, we scale
down the set $\mathcal{Y}$ so that $\mathcal{Y}\subseteq[0,1/\sigma]^{d}$, and
subsequently perturb the points by Gaussian variables with $\sigma=1$. Since
the number of iterations required to terminate is invariant under scaling of
the input point set, this is equivalent to the original instance.

Our analysis is based on the standard technique of proving that it is unlikely
that a sequence of iterations decreases the potential function by a small
amount. For this technique to work, we additionally require the potential
function to be bounded from above and from below with sufficiently high
probability. Since it is obvious that the potential is non-negative for any
clustering, it is enough to guarantee that the perturbed point set
$\mathcal{X}$ lies within the hypercube $[-D/2,D/2]^{d}$ for some finite $D$.
To that end, we have the following lemma.

######  Lemma 4.1.

Let $D=\sqrt{2n\ln(nkd)}$. Then .

Similar results to Lemma 4.1 can be found in previous works on the smoothed
analysis of algorithms on Gaussian-perturbed point sets [13, 1]. The only
difference in our version is the value of $D$. Hence, we omit the proof.

Lemma 4.1 allows us to assume that all points lie within $[-D/2,D/2]^{d}$
after the perturbation. Formally, we must take into account the failure event
that any point lies outside this hypercube. However, since the probability of
this event is at most $k^{-n}$, this adds only a negligible $+1$ to the
smoothed complexity bound which we prove in Theorem 1.2. We therefore ignore
the failure event in the sequel.

We need to show that we can approximate the gain of an iteration if we have a
good approximation to the cluster centers. Recall that
$\Delta_{x}(C_{i},C_{j})$ is the gain of moving a point $x$ from $C_{i}$ to
$C_{j}$. Since we wish to use approximations to the centers of $C_{i}$ and
$C_{j}$, it is convenient to define the variable

$\Delta_{x}^{|C_{i}|,|C_{j}|}(a,b)=\frac{|C_{i}|}{|C_{i}|-1}\|x-a\|^{2}-\frac{|C_{j}|}{|C_{j}|+1}\|x-b\|^{2}.$

This variable is the gain that would be incurred if the centers of $C_{i}$ and
$C_{j}$, with fixed sizes $|C_{i}|$ and $|C_{j}|$, were $a$ and $b$. Indeed,
note that
$\Delta_{x}^{|C_{i}|,|C_{j}|}(\operatorname{cm}(C_{i}),\operatorname{cm}(C_{j}))=\Delta_{x}(C_{i},C_{j})$.
When their intended values are clear from context, we will often omit the
superscripts $|C_{i}|$ and $|C_{j}|$ from $\Delta_{x}^{|C_{i}|,|C_{j}|}(a,b)$.

###  4.2 Approximating Iterations

Before we begin with the analysis proper, we provide a rough outline of our
analysis. Suppose we tile the hypercube $[-D/2,D/2]^{d}$ with a rectangular
grid of spacing $\epsilon$. Then any point in $[-D/2,D/2]^{d}$ is at a
distance of at most $\sqrt{d}\epsilon$ from some grid point. Since we need the
positions of the cluster centers $c_{i}=\operatorname{cm}(C_{i})$ for
$i\in[k]$, we guess $k$ grid points $c^{\prime}_{i}$ for their positions. If
we guess correctly, meaning $c^{\prime}_{i}$ is the grid point closest to
$c_{i}$ for each $i\in[k]$, then we can approximate the gain $\Delta$ of an
iteration by replacing the cluster centers with these grid points in the
formula for $\Delta$ (Lemma 4.2).

The price for this approximation is a union bound over all choices of the grid
points. However, we can compensate for this by noticing that, when we move a
point between clusters, we know exactly how the cluster centers move. Thus, if
the guessed grid points are good approximations, we can obtain new good
approximations by moving them the same amount. Thus, we only need to guess
once, and can use this guess for a sequence of iterations. Then we can bound
the probability that all iterations in this sequence yield a small
improvement.

######  Lemma 4.2.

Suppose the point $x$ moves from cluster $i$ to cluster $j$. Let $C_{i}$ and
$C_{j}$ denote the configurations of these clusters before this move, and let
$c_{i}=\operatorname{cm}(C_{i})$ and $c_{j}=\operatorname{cm}(C_{j})$. Let
$c_{i}^{\prime}$ and $c_{j}^{\prime}$ be two points such that
$\|c_{i}-c_{i}^{\prime}\|,\|c_{j}-c_{j}^{\prime}\|\leq\epsilon$ for some
$0\leq\epsilon\leq\sqrt{d}D$. Then

$|\Delta_{x}(C_{i},C_{j})-\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})|\leq
9\sqrt{d}D\epsilon,$

In particular, $\Delta_{x}(C_{i},C_{j})\in(0,\epsilon]$ implies
$|\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})|\leq 10\sqrt{d}D\epsilon$.

######  Proof 4.3.

Observe that

$\|x-c_{i}\|^{2}=\|x-c_{i}^{\prime}+c_{i}^{\prime}-c_{i}\|^{2}=\|x-c_{i}^{\prime}\|^{2}+\|c_{i}-c_{i}^{\prime}\|^{2}+2\langle
c_{i}^{\prime}-c_{i},x\rangle.$

Thus,

$\Delta_{x}(C_{i},C_{j})=\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})+\frac{|C_{i}|}{|C_{i}|-1}\mathopen{}\mathclose{{}\left(\|c_{i}-c_{i}^{\prime}\|^{2}+2\langle
c_{i}^{\prime}-c_{i},x\rangle}\right)\\\
-\frac{|C_{j}|}{|C_{j}|+1}\mathopen{}\mathclose{{}\left(\|c_{j}-c_{j}^{\prime}\|^{2}+2\langle
c_{j}^{\prime}-c_{j},x\rangle}\right).$

By the Cauchy-Schwarz inequality, $|\langle
c_{i}^{\prime}-c_{i},x\rangle|\leq\epsilon\cdot\|x\|$. Moreover, since
$x\in[-D/2,D/2]^{d}$, we have $\|x\|\leq\sqrt{d}D$. Similarly,
$\|c_{i}^{\prime}-c_{i}\|\leq\epsilon\leq\sqrt{d}D$ by assumption.

Moving $\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})$ to the left and taking an
absolute value, we then obtain

$|\Delta_{x}(C_{i},C_{j})-\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})|\leq\mathopen{}\mathclose{{}\left(\frac{|C_{i}|}{|C_{i}|-1}+\frac{|C_{j}|}{|C_{j}|+1}}\right)\cdot
3\sqrt{d}D\epsilon.$

To finish the proof, observe that by Lemma 2.3 the first term inside the
parentheses is at most $2$, while the second term is bounded by $1$. We then
have that $\Delta_{x}(C_{i},C_{j})\in(0,\epsilon]$ implies
$\Delta_{x}(c_{i}^{\prime},c_{j}^{\prime})\in(-9\sqrt{d}D\epsilon,(9\sqrt{d}D+1)\epsilon]$,
which yields the lemma.

In the following, we fix a set $A\subseteq\mathcal{X}$ of active points which
will move during a sequence of Hartigan’s method. We also fix the
configuration of the active points, the sizes of the clusters $|C_{1}|$ and
$|C_{2}|$, and the order $\pi:A\to[|A|]$ in which the points move. Observe
that these data also fix the sizes of the clusters whenever a new point moves.

While performing a sequence of iterations, the cluster centers move. Hence,
even if we have a good approximation to a cluster center, it may not remain a
good approximation after the iteration. However, if we know which points are
gained and lost by each cluster, then we can compute new good approximations
to the cluster centers from the old approximations. The following lemma
captures this intuition.

######  Lemma 4.4.

Let $t_{1}$, $t_{2}$ be two iterations of Hartigan’s method in a sequence in
which the points $A\subseteq\mathcal{X}$ move, with $t_{1}<t_{2}$. Suppose
between $t_{1}$ and $t_{2}$, cluster $i$ loses the points $S_{-}$ and gains
the points $S_{+}$. Let $c_{i}(t)$ denote the cluster center of cluster $i$
before $t$ takes place, and let $C_{i}^{t}$ denote its configuration before
$t$. Let $c_{i}^{\prime}(t_{1})\in\mathbb{R}^{d}$, and
$c_{i}^{\prime}(t_{2})=\frac{|C_{i}^{t_{1}}|}{|C_{i}^{t_{2}}|}c_{i}^{\prime}(t_{1})+\frac{1}{|C_{i}^{t_{2}}|}\mathopen{}\mathclose{{}\left(\sum_{x\in
S_{+}}x-\sum_{x\in S_{-}}x}\right)$. Then

$\|c_{i}^{\prime}(t_{2})-c_{i}(t_{2})\|=\frac{|C_{i}^{t_{1}}|}{|C_{i}^{t_{2}}|}\cdot\|c_{i}^{\prime}(t_{1})-c_{i}(t_{2})\|.$

Moreover, if $\|c_{i}^{\prime}(0)-c_{i}(0)\|\leq\epsilon$, then
$\|c_{i}^{\prime}(t_{j})-c_{i}(t_{j})\|\leq 2|A|\epsilon$ for all $j\in[|A|]$.

######  Proof 4.5.

Since the center of a cluster is defined as its center of mass, we can write

$|C_{i}^{t_{2}}|\operatorname{cm}(C_{i}^{t_{2}})=\sum_{x\in C_{i}^{t_{1}}\cup
S_{+}\setminus
S_{-}}x=|C_{i}^{t_{1}}|\operatorname{cm}(C_{i}^{t_{1}})+\sum_{x\in
S_{+}}x-\sum_{x\in S_{-}}x.$

Thus,

$|C_{i}^{t_{2}}|c_{i}(t_{2})=|C_{i}^{t_{1}}|c_{i}(t_{1})+\sum_{x\in
S_{+}}x-\sum_{x\in S_{-}}x.$

Observe then that

$\|c_{i}^{\prime}(t_{2})-c_{i}(t_{2})\|=\frac{|C_{i}^{t_{1}}|}{|C_{i}^{t_{2}}|}\cdot\|c_{i}^{\prime}(t_{1})-c_{i}(t_{1})\|.$

This proves the first claim. To prove the second claim, we obtain by
telescoping

$\|c_{i}(t_{j})-c_{i}^{\prime}(t_{j})\|=\frac{|C_{i}^{0}|}{|C_{i}^{t_{j}}|}\cdot\|c_{i}(0)-c_{i}^{\prime}(0)\|\leq(|A|+1)\epsilon\leq
2|A|\epsilon,$

since at most $|A|$ points are active during any subsequence.

###  4.3 Analyzing Sequences

We now know that we can closely approximate the gain of a sequence of
iterations, provided that we have good approximations to the cluster centers
at the start of the sequence. The next step is then to show that there is only
a small probability that such an approximate sequence improves the potential
by a small amount. For that, we first require the following technical lemma.

######  Lemma 4.6.

Let $X$ be a $d$-dimensional Gaussian random variable with arbitrary mean and
standard deviation $\sigma\leq 1$, and let $Z=a\|X\|^{2}+\langle v,X\rangle$
for fixed $a\in\mathbb{R}\setminus\\{0\\}$ and $v\in\mathbb{R}^{d}$. Then the
probability that $Z$ falls in an interval of size $\epsilon\leq 1$ is bounded
from above by
$O\mathopen{}\mathclose{{}\left(\frac{1}{|a|\sigma^{2}}\sqrt{\frac{\epsilon}{d}}}\right)$.

######  Proof 4.7.

Let $Z_{i}=aX_{i}^{2}+v_{i}X_{i}$, so that $Z=\sum_{i=1}^{d}Z_{i}$. We define
the auxiliary variable $\bar{Z}_{i}=Z_{i}+v_{i}^{2}/(4a)$. Since $a$ and $v$
are fixed, the densities of $Z$ and $\bar{Z}$ are identical up to translation,
and so we can analyze $\bar{Z}$ instead. Observe that
$\bar{Z}_{i}/a=\mathopen{}\mathclose{{}\left(X_{i}+\frac{v_{i}}{2a}}\right)^{2}$.
Thus, $\bar{Z}/a$ is equal in distribution to $\|Y\|^{2}$, where $Y$ is a
$d$-dimensional Gaussian variable with mean $\mu-v/(4a)$ and variance
$\sigma^{2}$. We see then that $\bar{Z}/a$ has the density of a non-central
chi-squared distribution.

For $\lambda\geq 0$, denote by $f(x,\lambda,d)$ the non-central
$d$-dimensional chi-squared density with non-centrality parameter $\lambda$
and standard deviation $\sigma$. Then [10]

$\displaystyle f(x,\lambda,d)$
$\displaystyle=\sum_{i=0}^{\infty}\frac{e^{-\lambda/2}(\lambda/2)^{i}}{i!}f(x,0,d+2i).$

Now observe that $f(x,0,d)$ is bounded from above by
$O(1/(\sqrt{d}\sigma^{2}))$ for $d\geq 2$. Adding in the scaling factor of
$1/|a|$ then yields the lemma for $d\geq 2$, since
$\epsilon\leq\sqrt{\epsilon}$ for $\epsilon\leq 1$.

For $d=1$, we have

$f(x,0,1)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\cdot\frac{e^{-\frac{x}{2\sigma^{2}}}}{\sqrt{x/\sigma^{2}}}.$

Let $I$ be an interval of size $\epsilon$. Then The first term is bounded by
$O(\epsilon/\sigma^{2})$ by the same argument we used for $d\geq 2$. For the
second term, we use the expression for $f(x,0,1)$ above to bound the integral
as

$\int_{I}f(x,0,1)\mathop{}\\!\mathrm{d}x\leq\frac{1}{\sqrt{2\pi\sigma^{2}}}\int_{0}^{\epsilon}\frac{e^{-\frac{x}{2\sigma^{2}}}}{\sqrt{x/\sigma}}\mathop{}\\!\mathrm{d}x=O(\sqrt{\epsilon}/\sigma)=O(\sqrt{\epsilon}/\sigma^{2}),$

where the final bound follows since $\sigma\leq 1$. This proves the lemma for
$d=1$ when we again add in the scaling factor $1/|a|$.

With Lemma 4.6, we can show that a single fixed approximate iteration is
unlikely to yield a small improvement.

######  Lemma 4.8.

Let $a,b\in\mathbb{R}^{d}$ be fixed. Let $\Delta_{x}(a,b)$ be the improvement
of the first move of $x$ in $S$, if the cluster centers in this iteration are
located at $a$ and $b$. Let $I$ be an interval of size $\epsilon\leq 1$. Then

######  Proof 4.9.

By Lemma 2.2, we have

$\Delta_{x}(a,b)=\frac{|C_{i}|}{|C_{i}|-1}\|x-a\|^{2}-\frac{|C_{j}|}{|C_{j}|+1}\|x-b\|^{2}\\\
=\mathopen{}\mathclose{{}\left(\frac{|C_{i}|}{|C_{i}|-1}-\frac{|C_{j}|}{|C_{j}|+1}}\right)\|x\|^{2}+\mathopen{}\mathclose{{}\left\langle
2\mathopen{}\mathclose{{}\left(\frac{|C_{j}|}{|C_{j}|+1}b-\frac{|C_{i}|}{|C_{i}|-1}a}\right),x}\right\rangle\\\
+\frac{|C_{i}|}{|C_{i}|-1}\|a\|^{2}-\frac{|C_{j}|}{|C_{j}|+1}\|b\|^{2},$

where $|C_{i}|$ and $|C_{j}|$ denote the sizes of clusters $i$ and $j$ before
the iteration, and we assume $x$ moves from cluster $i$ to cluster $j$.

Since the sizes of the clusters as well as $a$ and $b$ are fixed, the last
term in the above is fixed, and hence we may disregard it when analyzing .
Since $x$ is a Gaussian random variable, we can apply Lemma 4.6 to find

It remains to bound quantity in the inner brackets from below. Since each
cluster is bounded in size by $n$, we have

$\frac{|C_{i}|}{|C_{i}|-1}-\frac{|C_{j}|}{|C_{j}|+1}\geq\frac{n}{n-1}-\frac{n}{n+1}=\frac{2n}{(n-1)(n+1)}\geq\frac{1}{n},$

and we are done.

As stated at the start of the analysis, analyzing a single iteration is not
enough to prove Theorem 1.2. The following lemma extends Lemma 4.8 to a
sequence of iterations, given a fixed point set $A\subseteq\mathcal{X}$ that
moves in the sequence.

######  Lemma 4.10.

Fix an active set $A$ and starting cluster sizes $|C_{i}|$ for $i\in[k]$.
Moreover, fix an order $\pi:A\to[|A|]$ in which the points in $A$ move, i.e.,
$\pi(x)<\pi(y)$ means $x$ moves for the first time before $y$ moves for the
first time. Let $\Delta$ denote the minimum improvement of a sequence
satisfying these hypotheses over all possible configurations of
$\mathcal{X}\setminus A$. Then for $\epsilon\leq 1$,

######  Proof 4.11.

For $x\in A$, let $\Delta_{x}$ denote the improvement of the first move of
$x\in A$. We label the points in $A$ as $(x_{1},\ldots,x_{|A|})$ according to
$\pi$. Let $\Delta=(\Delta_{i})_{i=1}^{|A|}$.

To compute the vector $\Delta$, we would need to know the configuration and
positions of the points $P=\mathcal{X}\setminus A$, since these are required
to compute the $k$ cluster centers. However, if we had approximations to the
cluster centers in every iteration corresponding to the entries of $\Delta$,
then then we could compute an approximation to $\Delta$ by Lemma 4.2.

Since the cluster centers are convex combinations of points in
$[-D/2,D/2]^{d}$, we know that the cluster centers at the start of $S$ must
also lie in $[-D/2,D/2]^{d}$. Thus, there exist grid points $c_{i}^{\prime}$
($i\in[k]$) within a distance $\sqrt{d}\epsilon$ of the initial cluster
centers.

Knowing these grid points, we would like to apply Lemma 4.4 in order to update
the approximate cluster centers whenever a new point moves. We then need to
know the points gained and lost by each cluster between first moves of each
$x\in A$. Observe that to obtain this information, it suffices to know the
configuration of the active points before the first move of each $x\in A$.
Thus, we fix these configurations.

We collect the gain of each first move of a point in $A$, where we replace the
cluster centers by these approximations, into a vector $\Delta^{\prime}$. By
the reasoning above and by Lemma 4.2, if there exist initial cluster centers
$c_{i}$ ($i\in[k]$) such that $\Delta_{x}\in(0,\epsilon]$ for all $x\in A$,
then there exist grid points $c_{i}^{\prime}$, such that
$|\Delta_{x}^{\prime}|\leq 20|A|dD\epsilon$ for all $x\in A$.

By this reasoning, it suffices to obtain a bound on . We can then take a union
bound over these events for all $(D/\epsilon+1)^{kd}\leq(2D/\epsilon)^{kd}$
choices of $c_{i}^{\prime}$ for $i\in[k]$, and a union bound over the
configuration of $A$ before the first move of each $x\in A$.

To show that is bounded as desired, we consider the following algorithm.

  1. 1.

Set $t=1$.

  2. 2.

Reveal $x_{t}$, and compute
$\Delta_{x_{t}}(c_{i_{t}}^{\prime},c_{j_{t}}^{\prime})$, where $x_{t}$ moves
from $C_{i_{t}}$ to $C_{j_{t}}$.

  3. 3.

If $|\Delta_{x_{t}}(c_{i_{t}}^{\prime},c_{j_{t}}^{\prime})|>20|A|dD\epsilon$,
then return $\sf{GOOD}$ and halt.

  4. 4.

If $t=|A|$, return $\sf{BAD}$.

  5. 5.

Update the positions of the approximate cluster centers using Lemma 4.4.

  6. 6.

Continue executing moves in the sequence until we encounter the first move of
$x_{t+1}$. Observe that the information we fixed before executing this
algorithm suffices to compute approximations to the cluster centers whenever a
new point moves.

  7. 7.

Set $t\leftarrow t+1$ and go to step 2.

The sequence of iterations improves the potential by at most $\epsilon$ only
if the above algorithm returns $\sf{BAD}$. We now argue that

Let $\mathsf{BAD}_{t}$ be the event that the above algorithm loops for at
least $t$ iterations. Then . Since , we can immediately conclude that for all
$t\in\\{2,\ldots,|A|\\}$,

By Lemma 4.8, we have . Thus, is bounded as claimed.

Taking a union bound over all choices of the approximate grid points at the
start of the sequence yields the factor $(2D/\epsilon)^{kd}$. Finally, we must
take a union bound over the configuration of $A$ before the first move of each
$x\in A$, yielding a factor $k^{|A|^{2}}$, which concludes the proof.

Armed with Lemma 4.10, we can bound the probability that there exists a
sequence in which a fixed number of points moves, which improves the potential
by at most $\epsilon$.

######  Lemma 4.12.

Let $\Delta_{\mathrm{min}}$ denote the minimum improvement of any sequence of
moves in which exactly $4kd$ distinct points switch clusters. Then for
$\epsilon\leq 1$,

######  Proof 4.13.

Fix an active set $A$ of $4kd$ distinct points, an order $\pi:A\to[|A|]$ in
which the points in $A$ move, and the sizes of the two clusters at the start
of the sequence.

We have by Lemma 4.10

We conclude the proof by a union bound over the choices of $A$, $\pi$, and the
sizes of the clusters at the start of the sequence, which yields a factor of
at most $(4kd)^{4kd}\cdot n^{4kd+1}$.

With Lemma 4.12, we are in a position to prove the main result of this
section. The proof is essentially mechanical, following techniques used in
many previous smoothed analyses [1, 3, 4, 5, 6, 13].

*

######  Proof 4.14.

First, we recall that the point set $\mathcal{X}$ is contained in
$[-D/2,D/2]^{d}$. This yields an upper bound for the value of the potential
function for the initial clustering $C$,

$\Phi(C)=\sum_{i=1}^{k}\sum_{x\in C_{i}}\|x-\operatorname{cm}(C_{i})\|^{2}\leq
kndD^{2}.$

We divide the sequence of iterations executed by Hartigan’s method into
contiguous disjoint blocks during which exactly $4kd$ distinct points move. By
Lemma 4.12, we know that the probability that any such block yields a bad
improvement is small.

Let $T$ be the number of such blocks traversed by the heuristic before we
reach a local optimum. Then

This probability becomes nontrivial when

$t>\mathopen{}\mathclose{{}\left\lceil\frac{O(1)\cdot
k^{8kd+5}d^{11}D^{7}n^{9+\frac{1}{d}}}{\sigma^{8}}}\right\rceil=:t^{\prime}.$

Observe that $t^{\prime}=\Omega(kndD^{2})$, justifying our use of Lemma 4.12
above. Thus, we find

$\displaystyle\mathbb{E}(T)=\sum_{t=1}^{k^{n}}$

The upper limit of $k^{n}$ to the sum is simply the number of possible
clusterings of $n$ points into $k$ sets, which is a trivial upper bound to the
number of iterations. To conclude, we observe that any block in which exactly
$4kd$ distinct points move has a length of at most $k^{4kd}$, as otherwise
some clustering would show up twice. Thus, we multiply $\mathbb{E}(T)$ by
$k^{4kd}$ to obtain a bound for the smoothed complexity. Finally, we insert
the value of $D=\sqrt{2n\ln(nkd)}$.

##  5 Discussion

Theorems 1.1 and 1.2 provide some of the first rigorous theoretical results
concerning Hartigan’s method that have been found since Telgarsky & Vattani
explored the heuristic in 2010 [14]. Of course, many interesting open
questions still remain.

Theorem 1.1 establishes the existence of exponential-length sequences on the
line, but leaves open the possibility that a local optimum may be reachable
more efficiently by a different improving sequence. To be precise: given an
instance of $k$-means clustering on the line and an initial clustering, does
there always exist a sequence of iterations of Hartigan’s method of length
$\operatorname{poly}(n,k)$ starting from this clustering and ending in a local
optimum? Although the $d=1$ case appears very restricted at first sight, this
question seems surprisingly difficult to answer.

In addition, the construction we use in Theorem 1.1 requires $k=\Theta(n)$
clusters. This opens up the question whether similar worst-case constructions
can be made using fewer, perhaps even $O(1)$, clusters. Note that this is not
true for Lloyd’s method, since the number of iterations of Lloyd’s method is
bounded by $n^{O(kd)}$ [9], which is polynomial for $k,d\in O(1)$.

Theorem 1.2 entails, to our knowledge, the first step towards settling the
conjecture by Telgarsky & Vattani [14] that Hartigan’s method has polynomial
smoothed complexity. Our result is reminiscent of the smoothed complexity
bound of Lloyd’s method obtained in 2009 by Manthey & Röglin [12], which is
$k^{kd}\cdot\operatorname{poly}(n,1/\sigma)$. In the case of Lloyd’s method,
the smoothed complexity was later settled to
$\operatorname{poly}(n,k,d,1/\sigma)$ [1].

Observe that our bound is polynomial for constant $k$ and $d$, and even for
$kd\log k\in O(\log n)$. While this is certainly an improvement over the
trivial upper bound of $k^{n}$, it falls short of a true polynomial bound. We
hope that our result can function as a first step to a
$\operatorname{poly}(n,k,d,1/\sigma)$ smoothed complexity bound of Hartigan’s
method.

We remark that the exponents in the bound in Theorem 1.2 can be easily
improved by a constant factor for $d\geq 2$. The reason is that in Lemma 4.6,
the factor $\sqrt{\epsilon}$ emerges from the $d=1$ case, while for $d\geq 2$
this would be $\epsilon$. We chose to combine these cases for the sake of
keeping the analysis simple, as we expect the bound in Theorem 1.2 would be
far from optimal regardless.

The main challenges to improving this bound are twofold. First, we must take a
union bound over the configuration of the active points each time we apply
Lemma 4.4, yielding factors of $k^{O(kd)}$. Second, we must analyze sequences
in which $\Theta(kd)$ points move in order to guarantee a significant
potential decrease. This incurs a factor of the length of such a sequence,
which is another source of a factor $k^{O(kd)}$

One avenue for resolving this problem might be to analyze shorter sequences in
which a significant number of points move. Angel et al. used such an approach
in their analysis of the Flip heuristic for Max-Cut. They identify in any
sequence $L$ of moves a shorter subsequence $B$, such that the number of
unique vertices that flip in $B$ is linear in the length of $B$. The major
challenge is then to find sufficient independence in such a short subsequence,
which in our case seems challenging, as we need to compensate for a factor
$\epsilon^{-kd}$ in Lemma 4.10.

## References

  * [1] David Arthur, Bodo Manthey, and Heiko Röglin.  Smoothed Analysis of the k-Means Method.  Journal of the ACM, 58(5):19:1–19:31, October 2011.  doi:10.1145/2027216.2027217.
  * [2] David Arthur and Sergei Vassilvitskii.  K-means++: The advantages of careful seeding.  In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’07, pages 1027–1035, USA, January 2007\. Society for Industrial and Applied Mathematics.
  * [3] Matthias Englert, Heiko Röglin, and Berthold Vöcking.  Worst Case and Probabilistic Analysis of the 2-Opt Algorithm for the TSP.  Algorithmica, 68(1):190–264, January 2014.  doi:10.1007/s00453-013-9801-4.
  * [4] Matthias Englert, Heiko Röglin, and Berthold Vöcking.  Smoothed Analysis of the 2-Opt Algorithm for the General TSP.  ACM Transactions on Algorithms, 13(1):10:1–10:15, September 2016\.  doi:10.1145/2972953.
  * [5] Michael Etscheid and Heiko Röglin.  Smoothed Analysis of the Squared Euclidean Maximum-Cut Problem.  In Nikhil Bansal and Irene Finocchi, editors, Algorithms - ESA 2015, Lecture Notes in Computer Science, pages 509–520, Berlin, Heidelberg, 2015. Springer.  doi:10.1007/978-3-662-48350-3_43.
  * [6] Michael Etscheid and Heiko Röglin.  Smoothed Analysis of Local Search for the Maximum-Cut Problem.  ACM Transactions on Algorithms, 13(2):25:1–25:12, March 2017.  doi:10.1145/3011870.
  * [7] Gurobi Optimization LLC.  Gurobi Optimizer Reference Manual.  Gurobi Optimization, LLC, 2023.
  * [8] J. A. Hartigan and M. A. Wong.  Algorithm AS 136: A K-Means Clustering Algorithm.  Journal of the Royal Statistical Society. Series C (Applied Statistics), 28(1):100–108, 1979.  arXiv:2346830, doi:10.2307/2346830.
  * [9] M Inaba, Naoki Katoh, and Hiroshi Imai.  Variance-based k-clustering algorithms by Voronoi diagrams and randomization.  IEICE Transactions on Information and Systems, E83D, June 2000.
  * [10] Norman L. Johnson, Samuel Kotz, and Narayanaswamy Balakrishnan.  Continuous Univariate Distributions, Volume 2.  John Wiley & Sons, May 1995.
  * [11] S. Lloyd.  Least squares quantization in PCM.  IEEE Transactions on Information Theory, 28(2):129–137, March 1982\.  doi:10.1109/TIT.1982.1056489.
  * [12] Bodo Manthey and Heiko Röglin.  Improved Smoothed Analysis of the k-Means Method.  In Proceedings of the 2009 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), Proceedings, pages 461–470. Society for Industrial and Applied Mathematics, January 2009.  doi:10.1137/1.9781611973068.51.
  * [13] Bodo Manthey and Rianne Veenstra.  Smoothed Analysis of the 2-Opt Heuristic for the TSP: Polynomial Bounds for Gaussian Noise.  In Leizhen Cai, Siu-Wing Cheng, and Tak-Wah Lam, editors, Algorithms and Computation, Lecture Notes in Computer Science, pages 579–589, Berlin, Heidelberg, 2013. Springer.  doi:10.1007/978-3-642-45030-3_54.
  * [14] Matus Telgarsky and Andrea Vattani.  Hartigan’s Method: K-means Clustering without Voronoi.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 820–827. JMLR Workshop and Conference Proceedings, March 2010.
  * [15] Andrea Vattani.  K-means Requires Exponentially Many Iterations Even in the Plane.  Discrete & Computational Geometry, 45(4):596–616, June 2011.  doi:10.1007/s00454-011-9340-1.

##  Appendix A Proof of Theorem 1.1

*

######  Proof A.1.

Inspecting the sequence detailed in Section 3.1, we observe that the morning
phase consists of three iterations, the afternoon phase of one iteration, and
the waking up phase of six iterations. Each iteration yields a single
inequality, and thus we must verify ten inequalities in total.

Recall the notation

$\Delta_{x}(S,T)=\frac{|S|}{|S|-1}\|x-\operatorname{cm}(S)\|^{2}-\frac{|T|}{|T|+1}\|x-\operatorname{cm}(T)\|^{2},$

which is the gain obtained by moving the point $x$ from $S$ to $T$. Then each
iteration will yield an inequality of the form $\Delta_{x}(S,T)>0$.

Note that we may assume w.l.o.g. that $G_{i}$ is the unit gadget, and as such
that for each triple of corresponding points $x_{i-1},x_{i},x_{i+1}$ it holds
that $x_{i-1}=(x_{i}-t_{0})/5$ and $x_{i+1}=5x_{i}+t_{0}$.

#### Morning.

This phase yields the inequalities:

$\displaystyle\Delta_{p_{i}}(\\{p_{i},q_{i},b_{i}\\},\\{a_{i-1},q_{i-1}\\})>0,$
(1) $\displaystyle\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{a_{i}\\})>0,$ (2)
$\displaystyle\Delta_{q_{i}}(\\{b_{i},q_{i}\\},\\{a_{i},p_{i}\\})>0.$ (3)

For the second of these inequalities, we refer to the description of the
wakeup phase, which shows that $\mathcal{C}_{1}(G_{i-1})$ takes the value
$\\{a_{i-1},p_{i}\\}$ at the moment $p_{i}$ moves back to a cluster of
$G_{i}$.

#### Afternoon.

This phase only yields

$\displaystyle\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{a_{i-1},q_{i-1}\\})>0.$
(4)

#### Waking up.

Here, we have

$\displaystyle\Delta_{p_{i+1}}(\\{p_{i+1},q_{i+1},b_{i+1}\\},\\{a_{i},q_{i}\\})>0,$
$\displaystyle\Delta_{p_{i+1}}(\\{p_{i+1},q_{i+1},a_{i+1}\\},\\{a_{i},q_{i}\\})>0,$
$\displaystyle\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{b_{i}\\})>0,$
$\displaystyle\Delta_{q_{i}}(\\{a_{i},q_{i},p_{i+1}\\},\\{b_{i}\\})>0,$
$\displaystyle\Delta_{p_{i+1}}(\\{a_{i},p_{i+1}\\},\\{a_{i+1}\\})>0,$
$\displaystyle\Delta_{p_{i+1}}(\\{a_{i},p_{i+1}\\},\\{b_{i+1}\\})>0.$

Observe that some of these inequalities are equivalent to one another, or to
the inequalities specified in the morning and afternoon phases. The only
unique inequalities are

$\displaystyle\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{b_{i}\\})>0,$ (5)
$\displaystyle\Delta_{q_{i}}(\\{a_{i},q_{i},p_{i+1}\\},\\{p_{i},b_{i}\\})>0.$
(6)

#### Leaf gadget.

Finally, we consider the case where $G_{i-1}$ is the leaf gadget. This yields
the additional inequalities

$\displaystyle\Delta_{p_{i}}(\\{p_{i},q_{i},b_{i}\\},\\{f\\})>0,$ (7)
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{a_{i}\\})>0,$ (8)
$\displaystyle\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{f\\})>0,$ (9)
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{b_{i}\\})>0.$ (10)

Verifying the inequalities is now simply a matter of plugging in the values
from Table 1.

#### Morning.

We start by computing some necessary quantities:

$\displaystyle q_{i-1}$ $\displaystyle=(q_{i}-t_{0})/5=(13-8)/5=1,$
$\displaystyle a_{i-1}$ $\displaystyle=(a_{i}-t_{0})/5=(9-8)/5=0.2,$
$\displaystyle\operatorname{cm}(\\{p_{i},q_{i},b_{i}\\})$
$\displaystyle=\frac{1}{3}(5+13+6)/3=8$
$\displaystyle\operatorname{cm}(\\{a_{i-1},q_{i-1}\\})$
$\displaystyle=\frac{1}{2}(0.2+1)=0.6,$
$\displaystyle\operatorname{cm}(\\{a_{i-1},p_{i}\\})$
$\displaystyle=\frac{1}{2}(0.2+5)=2.6,$
$\displaystyle\operatorname{cm}(\\{a_{i}\\})$ $\displaystyle=a_{i}=9,$
$\displaystyle\operatorname{cm}(\\{b_{i},q_{i}\\})$
$\displaystyle=\frac{1}{2}(6+13)=9.5,$
$\displaystyle\operatorname{cm}(\\{a_{i},p_{i}\\})$
$\displaystyle=\frac{1}{2}(9+5)=7.$

Using these quantities, we directly compute:

$\displaystyle\Delta_{p_{i}}(\\{p_{i},q_{i},b_{i}\\},\\{a_{i-1},q_{i-1}\\})$
$\displaystyle=\frac{3}{2}\mathopen{}\mathclose{{}\left(5-8}\right)^{2}-\frac{2}{3}\mathopen{}\mathclose{{}\left(5-0.6}\right)^{2}\approx
0.5933>0,$ $\displaystyle\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{a_{i}\\})$
$\displaystyle=2\cdot\mathopen{}\mathclose{{}\left(5-2.6}\right)^{2}-\frac{1}{2}\cdot\mathopen{}\mathclose{{}\left(5-9}\right)^{2}=3.52>0,$
$\displaystyle\Delta_{q_{i}}(\\{b_{i},q_{i}\\},\\{a_{i},p_{i}\\})$
$\displaystyle=2\cdot\mathopen{}\mathclose{{}\left(13-9.5}\right)^{2}-\frac{1}{2}\cdot\mathopen{}\mathclose{{}\left(13-7}\right)^{2}=6.5>0.$

Thus, Equations 1, 2 and 3 are satisfied by the points given in Table 1.

#### Afternoon.

The inequality we must verify for this phase is

$\displaystyle\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{a_{i-1},q_{i-1}\\})>0.$

We already know that $\operatorname{cm}(\\{a_{i-1},q_{i-1})=0.6$. Moreover,

$\operatorname{cm}(\\{a_{i},p_{i},q_{i}\\})=\frac{1}{3}(9+5+13)=9.$

Thus,

$\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{a_{i-1},q_{i-1}\\})=\frac{3}{2}(5-9)^{2}-\frac{2}{3}(5-0.6)^{2}\approx
11.093>0$

as claimed.

#### Waking up.

We must verify

$\displaystyle\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{b_{i}\\})>0,$
$\displaystyle\Delta_{q_{i}}(\\{a_{i},q_{i},p_{i+1}\\},\\{p_{i},b_{i}\\})>0.$

We already know that $\operatorname{cm}(\\{a_{i-1},p_{i}\\})=2.6$. Thus,

$\Delta_{p_{i}}(\\{a_{i-1},p_{i}\\},\\{b_{i}\\})=2\cdot(5-2.6)^{2}-0.5\cdot(5-6)^{2}=11.02>0.$

Next, we have $p_{i+1}=5p_{i}+t_{0}=33$. Thus,

$\displaystyle\operatorname{cm}(\\{a_{i},q_{i},p_{i+1}\\})$
$\displaystyle=\frac{1}{3}(9+13+33)=55/3,$
$\displaystyle\operatorname{cm}(\\{p_{i},b_{i}\\})$
$\displaystyle=\frac{1}{2}(5+6)=5.5,$

enabling us to verify that

$\Delta_{q_{i}}(\\{a_{i},q_{i},p_{i+1}\\},\\{p_{i},b_{i}\\})=\frac{3}{2}\cdot(13-55/3)^{2}-\frac{2}{3}\cdot(13-5.5)^{2}\approx
5.167>0.$

#### Leaf gadget.

We must verify

$\displaystyle\Delta_{p_{i}}(\\{p_{i},q_{i},b_{i}\\},\\{f\\})>0,$
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{a_{i}\\})>0,$
$\displaystyle\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{f\\})>0,$
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{b_{i}\\})>0.$

We first compute

$\displaystyle\operatorname{cm}(\\{p_{i},q_{i},b_{i}\\})$
$\displaystyle=\frac{1}{3}(5+13+6)=8,$
$\displaystyle\operatorname{cm}(\\{f,p_{i}\\})$
$\displaystyle=\frac{1}{2}(0+5)=2.5.$

Finally,

$\displaystyle\Delta_{p_{i}}(\\{p_{i},q_{i},b_{i}\\},\\{f\\})$
$\displaystyle=\frac{3}{2}(5-8)^{2}-0.5\cdot(5-0)^{2}=1>0,$
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{a_{i}\\})$
$\displaystyle=2\cdot(5-2.5)^{2}-0.5\cdot(5-9)^{2}=4.5>0,$
$\displaystyle\Delta_{p_{i}}(\\{a_{i},p_{i},q_{i}\\},\\{f\\})$
$\displaystyle=\frac{3}{2}\cdot(5-9)^{2}-0.5\cdot(5-0)^{2}=11.5>0,$
$\displaystyle\Delta_{p_{i}}(\\{f,p_{i}\\},\\{b_{i}\\})$
$\displaystyle=2\cdot(5-2.5)^{2}-0.5\cdot(5-6)^{2}=12>0.$

This concludes the verification of Equations 1 to A and therefore the proof of
Theorem 1.1.


World impact of kernel European Union 9 countries
from Google matrix analysis of the world trade network

Justin Loye1,2, Leonardo Ermann3,4, Dima L. Shepelyansky1,∗

1 Laboratoire de Physique Théorique du CNRS, IRSAMC, Université de Toulouse,
CNRS, UPS, 31062 Toulouse, France

2 Institut de Recherche en Informatique de Toulouse, Université de Toulouse,
UPS, 31062 Toulouse, France

3 Departamento de Física Teórica, GIyA, Comisión Nacional de Energía Atómica.
Av. del Libertador 8250, 1429 Buenos Aires, Argentina

4 Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Buenos
Aires, Argentina

$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

We use the United Nations COMTRADE database for analysis of the multiproduct
world trade network. With this data, considered for years 2012-2018, we
determined the world trade impact of the Kernel of EU 9 countries (KEU9),
being Austria, Belgium, France, Germany, Italy, Luxembourg, Netherlands,
Portugal, Spain, considered as one united country. We apply the advanced
Google matrix analysis for investigation of the influence of KEU9 and show
that KEU9 takes the top trade network rank positions thus becoming the main
player of the world trade being ahead of USA and China. Our network analysis
provides additional mathematical grounds in favor of the recent proposal [2]
of KEU9 super-union which is based only on historical, political and economy
basis.

Dated: October 21, 2020

## Introduction

The economy of European Union (EU) is considered as the second world largest
economy after United States (US) [3] even if there are also other opinions
placing China (CN) on the first position [4]. At present EU includes 27 member
states and about 447 million population [5]. While the global EU economy and
population are really huge the political action of member states [6] is not
always coherent pushing in some cases in different directions. Due to this
reason there is a proposal, pushed forward by Christian Saint-Etienne, to
consider the Kernel EU 9 (KEU9) states (or countries), which are tightly
linked by historical, political and economic relations, as a strongly united
kernel group of EU that would allow to perform coherent actions of these KEU9
states [2] (follow also discussion of this proposal at [7]). These 9 kernel
states include Austria (AT), Belgium (BE), France (FR), Germany (DE), Italy
(IT), Luxembourg (LU), Netherlands (NL), Portugal (PT), Spain (ES) [2] with
the total population of about 305 millions [5].

A variety of arguments in favor of possible coherent political and economic
actions of KEU9 group is presented and analyzed in [2]. However, this analysis
is not based on detailed mathematical grounds pushing forward arguments of
historical, political and economic heuristic reasons. Here, we put forward the
mathematical foundations for this KEU9 proposal presenting the mathematical
and statistical analysis of the world trade database of UN COMTRADE [8]. This
database presents an exceptional variety of data on trade exchange between all
world UN registered countries on a scale of more than 50 years with more than
$10^{4}$ of trade commodities (products). The transactions are expressed in
their dollar (USD) values of a given year. The World Trade Organization (WTO)
Statistical Review 2018 [9] demonstrates the vital importance of the
international trade between countries for their development and progress. Also
the whole world economy is deeply influenced by the world trade [10]. Hence,
this database is well appropriate for verification of how strong and important
is the world influence of KEU9 group on the world economy. Thus here we use
the UN COMTRADE database [8] for mathematical and statistical analysis of
heuristic arguments presented in favor of KEU9 in [2].

The trade transfer between countries represents the multiproduct World Trade
Network (WTN). The modern methods of Google matrix approach [11, 12, 13] are
well suited for the analysis of transactions on the WTN. The detailed
description of Google matrix applications to WTN are described in [14, 15, 16,
17]. Here we apply these methods considering KEU9 countries as one country
thus excluding trade transfers between them and keeping only ingoing and
outgoing trade flows to this group from other countries.

We point that various research groups investigated the statistical properties
of WTN (see e.g. [18, 19, 20, 21, 22, 23, 24]). However, as discussed in [14,
15, 16] the Google matrix approach has significant advantages for analysis of
weighted directed trade networks since it takes into account multiple
iterative transactions and thus provides a new and more detailed analysis of
trade influence propagation compared to the usual approach based on export and
import flows.

## Materials and Methods

### Google matrix construction of WTN

We consider the trade exchange between $N_{c}=186$ (185 countries + KEU9)
world countries and $N_{p}=10$ products given by 1 digit from the the Standard
International Trade Classification (SITC) Rev. 1, and for years 2012, 2014,
2016, 2018 taken from UN COMTRADE [8]. These 10 products contain all smaller
subdivided specific products which number goes up to $\sim 10^{4}$. The list
of these 10 products is given in Table 1. The list of world countries is
available at [14, 15]. Following the approach developed in [14, 15] we obtain
$N_{p}$ money matrices $M^{p}_{c,c^{\prime}}$ which give product $p$ transfer
(in USD) from country $c^{\prime}$ to country $c$. The Google matrices $G$ for
the direct trade flow and $G^{*}$ for the inverted trade flow have the size of
$N=N_{c}N_{p}=1860$ nodes. They are constructed by normalization of all column
of outgoing weighted links to unity. There is also the part with a damping
factor $\alpha=0.5$ describing random trade-surfer jumps to all nodes with a
certain personalized vector taking into account the weight of each product in
the global trade volume. The construction procedure of $G$ and $G^{*}$ is
described in detail in [15, 16]. The general properties and various examples
of Google matrices of various networks are given in [11, 12, 13].

The stationary probability distribution of Markov transitions described by the
Google matrix $G$ is given by the PageRank vector $P$ with maximal eigenvalue
$\lambda=1$: $GP=\lambda P=P$ [11, 12]. For the inverted flow described by
$G^{*}$ matrix we have similarly the CheiRank vector $P^{*}$, being the
eigenvector of $G^{*}P^{*}=P^{*}$. The importance and detailed statistical
analysis of the CheiRank vector were demonstrated in [25] (see also [14, 15,
26]). We define PageRank $K$ and CheiRank $K^{*}$ indexes by monotonic
ordering of probabilities of PageRank vector $P$ and of CheiRank vector
$P^{*}$ as $P(K)\geq P(K+1)$ and $P^{*}(K^{*})\geq P^{*}(K^{*}+1)$ with
$K,K^{*}=1,\ldots,N$. By taking a sum over all products $p$ we obtain the
PageRank and CheiRank probabilities of a given country as
$P_{c}=\sum_{p}P_{cp}$ and ${P^{*}}_{c}=\sum_{p}{P^{*}}_{cp}$ (and in a
similar way product probabilities $P_{p},{P^{*}}_{p}$) [15, 16]). From these
probabilities we obtain the related indexes $K_{c},{K^{*}}_{c}$. In a similar
way we define from import and export trade volume the probabilities
$\hat{P}_{p}$, $\hat{P}^{*}_{p}$, $\hat{P}_{c}$, $\hat{P}^{*}_{c}$,
$\hat{P}_{pc}$, $\hat{P}^{*}_{pc}$ and corresponding indexes $\hat{K}_{p}$,
$\hat{K}^{*}_{p}$, $\hat{K}_{c}$, $\hat{K}^{*}_{c}$, $\hat{K}$, $\hat{K}^{*}$
(the import and export probabilities are normalized to unity via the total
import and export volumes, see details in [15, 16]). We note that
qualitatively PageRank probability is proportional to the volume of ingoing
trade flow and CheiRank respectively to outgoing flow. Thus, approximately we
can say that the high import gives a high PageRank $P$ and a high export a
high CheiRank $P^{*}$ probabilities.

### Reduced Google matrix

We also use the REGOMAX algorithm described in detail in [27, 28]. This
algorithm allows to compute efficiently a reduced Google matrix $G_{R}$ of
size $N_{r}\times N_{r}$ that accounts all transitions of direct and indirect
pathways happening in the full Google matrix $G$ between $N_{r}$ nodes of
interest. For the selected $N_{r}$ nodes their PageRank probabilities are the
same as for the global network with $N$ nodes (up to a constant multiplicative
factor which takes into account that the sum of PageRank probabilities over
$N_{r}$ nodes is unity). The matrix $G_{R}$ can be presented as as a sum of
three matrix components that clearly distinguish direct and indirect
interactions: $G_{\mathrm{R}}=G_{rr}+G_{\mathrm{pr}}+G_{\mathrm{qr}}$ [28].
Thus $G_{rr}$ is produced by the direct links between selected $N_{r}$ nodes
in the global network of $N\gg N_{r}$ nodes. The component $G_{pr}$ is rather
close to the matrix in which each column is given by the PageRank vector
$P_{r}$ (up to a constant multiplier). Due to that $G_{pr}$ does not give much
information about direct and indirect links between selected $N_{r}$ nodes.
The most interesting and nontrivial role is played by the component $G_{qr}$,
which accumulates the contribution of all indirect links between selected
$N_{r}$ nodes appearing due to multiple pathways via the global network of $N$
nodes. The exact formulas for these three components of $G_{R}$ are given in
[27, 28].

### Sensitivity of trade balance

Following [15, 16], we use the trade balance of a given country with PageRank
and CheiRank probabilities defined as
$B_{c}=(P^{*}_{c}-P_{c})/(P^{*}_{c}+P_{c})$. In a similar way we have from
ImportRank and ExportRank probabilities as
$\hat{B}_{c}=({\hat{P}^{*}}_{c}-\hat{P}_{c})/({\hat{P}^{*}}_{c}+\hat{P}_{c})$.
The sensitivity of trade balance $B_{c}$ to the price of energy or machinery
can be obtained from the change of corresponding money volume flow related to
SITC Rev.1 code $p=3$ (mineral fuels) or $p=7$ (machinery) by multiplying it
by $(1+\delta)$, then computing all rank probabilities and the derivative
$dB_{c}/d\delta$.

The efficiency of the above Google matrix methods has been demonstrated not
only for the WTN but also for variety of other directed networks including
Wikipedia networks [26, 29, 30, 31] and biological networks of protein-protein
interactions [32, 33].

## Results

### CheiRank and PageRank of countries

We start the presentation of obtained results from showing the distribution of
world countries on the plane of CheiRank-PageRank indexes $(K,K^{*})$ given in
Fig. 1 (left panel). Here, for a better visibility, we show only countries
with $K,K^{*}\leq 60$, each country is marked by a circle with its flag. For a
comparison we also present in Fig. 1 (right panel) the distribution of
countries on the plane of ExportRank-ImportRank $\hat{K}$, $\hat{K}^{*}$ (in
both panels, for compactness, we keep index $K$ which in fact corresponds to
$K_{c}$ index of a country obtained by a summation over all products). The top
20 countries with their indexes are given in Table 2.

The main feature of Fig. 1 and Table 2 is that KEU9 takes the top leading
position in PageRank and CheiRank indexes $K,K^{*}$ in 2018 (this leadership
is also present in other studied years 2012, 2014, 206 as it is shown in
Supporting Information (SupInf) Fig.S1). This result is significantly
different from the Import-Export volume ranking where in 2018 China is leading
in export and USA in import. We argue that the Google matrix analysis via
PageRank and CheiRank treats in a deeper way the multiplicity of trade
relations between world countries compared to the standard Import-Export
approach which takes into account only one step trade links.

Another important feature of Google matrix analysis is a significant
improvement of positions of certain countries compared to their usual Import-
Export ranking (see Fig. 1, Table 2). Thus Russia moves to the fourth CheiRank
position $K^{*}=6$ compared to its ExportRank $\hat{K}^{*}=7$. Also India has
strong CheiRank-PageRank position $K^{*}=7,K=5$ compared to Export-ImportRanks
$\hat{K}^{*}=11$, $\hat{K}=7$. Also there is a significant reduction of
positions of Switzerland from $\hat{K}^{*}=12$, $\hat{K}=11$ to
$K^{*}=18,K=14$. In our opinion these results demonstrate a significant hidden
power or weakness of trade relations of certain countries due to the
multiplicity and variety of their trade relations which are not visible in a
standard Export-Import approach.

The main message of the results of this part is the world top leading position
of KEU9 in CheiRank and PageRank trade that gives confirmation of the strength
and importance of KEU9 countries discussed in [2].

### Trade balance of countries

We present the world map of trade balance $B_{c}$ of countries obtained from
CheiRank-PageRank and ExportRank-ImportRank probabilities in Fig. 2 for year
2018 (other years 2012, 2014, 2016 are given in SupInf FigS2; the
distributions of import and export of countries for all years are shown in
Figs.S3,S4). The comparison of two ways of balance computation shows that
Export-Import approach does not capture the influence of Russia and China on
the world trade exchange. In contrast the CheiRank-PageRank approach directly
highlights the multistep network influence of Russia and China on the world
trade flows and their balance. We also see a strong positive CheiRank-PageRank
balance for Japan. In both approaches the balance of US is close slightly
negative. There is a relative increase of KEU9 balance in CheiRank-PageRank
description compared to the standard Export-Import one. We attribute this to
the fact that CheiRank-PageRank description takes into account the
multiplicity of trade links which better describes a broad variety of KEU9
trade exchange.

We note that by definition we have the balance bounds $-1\leq B_{c}\leq 1$.
The actual obtained bounds are $(-0.94,0.73)$ and $(-0.25,0.31)$ for Export-
Import and CheiRank-PageRank descriptions (see Fig. 2). We attribute a
reduction of bounds for the latter case to a multiplicity of network links
that reduce fluctuations in trade exchange.

### Sensitivity of trade balance to specific products

As described above we determine the sensitivity of trade balance of countries
$dB_{c}/d\delta_{s}$ to specific products using the sensitivity definition
from CheiRank-PageRank and Export-Import probabilities. The sensitivity
results for $s=3$ (mineral fuels) are given in Fig. 3. The CheiRank-PageRank
approach shows that the most profitable countries with the highest values of
$dB_{c}/d\delta_{3}$ are Saudi Arabia and Russia (Kazakhstan also has high
sensitivity). This is rather natural since these countries are the highest
petroleum producers. The strongly negative impact is well visible for
Australia, China and countries of Latin America. USA and KEU9 sensitivities
being close to zero.

In contrast the sensitivity from Export-Import approach gives of the top
position Algeria (followed by Brunei). Among countries with strongly negative
sensitivities we have India, Pakistan and China while Australia is slightly
positive. In this Export-Import approach USA is slightly positive and KEU9 is
slightly negative.

This shows a significant difference between the usual Export-Import analysis
and the Google matrix approach. We argue that the latter approach takes into
account the multiplicity of trade links and flows thus highlighting in a
better way the multistep trade relations between countries.

The sensitivities of countries to the product $s=7$ (machinery) is shown in
Fig. 4. Here both approaches give the most positive countries being Japan,
S.Korea and China. In the Export-Import approach KEU9 has a bit higher
positive sensitivity compared to the CheiRank-PageRank method. Thus we have
for both methods of CheiRank-PageRank and Export-Import: KEU9
$dB_{c}/d\delta_{7}=0.015$, $d\hat{B}_{c}/d\delta_{7}=0.043$; slightly
negative values for USA $dB_{c}/d\delta_{7}=-0.019$,
$d\hat{B}_{c}/d\delta_{7}=-0.027$; Russia has strongly negative values
$dB_{c}/d\delta_{7}=-0.145$, $d\hat{B}_{c}/d\delta_{7}=-0.169$.

In the above Figs. 3, 4 we presented results for year 2018. The same type of
data for years 2012, 2014, 2016 are given in SupInf Figs. S5,S6.

Above we considered the sensitivity of trade balance to a global price
variation of a given product applicable to the whole world with a homogeneous
price increase of product for all countries. It is also interesting to
consider the sensitivity of country trade balance when the product price is
changed only by one country. In this way we obtain the sensitivity
$dB_{c}/d\delta_{cs}$ of countries to a product of a given country. This
specific sensitivity is shown in Fig. 5 in respect to price variation of $s=7$
(machinery) from KEU9, USA, and China in year 2018. The Export-Import approach
gives strongly positive sensitivity for the country which increased price of
machinery (respectively KEU9, USA, China). All other countries have
sensitivity close to zero or negative. The result from CheiRank-PageRank
analysis is different. For KEU9 machinery price increase the positive
sensitivity is obtained for Czechia, Slovakia, Hungary (with the sensitivity
values $0.028,0.017,0.015$ respectively). For USA case the positive
sensitivity is obtained for Mexico, Canada with repective values
$0.045,0.014$. For China case the positive sensitivity is obtained for Korea,
Philippines, Malaysia (with sensitivity respective values
$0.031,0.023,0.014$). This increase is related to strong network links between
these countries well captured by the Google matrix analysis.

In Fig. 6 we show the sensitivity $dB_{c}/d\delta_{cs}$ from both approaches
for $s=3$ (mineral fuels) of Russia. Again as for $s=7$ we see that the
Export-Import approach gives the strong positive sensitivity only for Russia.
In contrast the CheiRank-PageRank approach shows that Uzbekistan, Kazakhstan,
Ukraine (with values $dB_{c}/d\delta_{cs}=0.032,0.021,0.012$ respectively)
also gain the positive sensitivity in the case of price increase of $s=3$ of
Russia. This also confirms the strength of Google matrix analysis which
captures multiple trade links between countries.

### Sensitivity of trade balance to labor cost

It is interesting to analyze the sensitivity of a country trade balance
$dB_{c}/d\sigma_{c^{\prime}}$ to a labor cost variation in a given country.
This analysis is done by increasing the price of all products of a given
country by a factor $1+\sigma_{c^{\prime}}$ followed by a renormalization of
sum all column elements to unity. Such an approach has been developed and
studied in [34] for the world economic activities from World Trade
Organization data. Here, at the difference of price shock of one product, the
price increase affects all product flows from a given country corresponding to
a global increase of labor cost in a given country. Of course, the price
increase is considered to be very small corresponding to the linear response
regime. The labor cost sensitivity $dB_{c}/d\sigma_{c^{\prime}}$ is computed
numerically in the same manner as the product sensitivity
$dB_{c}/d\delta_{c^{\prime}}$ discussed above.

As discussed in [34] the most strong labor cost sensitivity
$dB_{c}/d\delta_{c^{\prime}}$ is naturally obtained for the country itself
with $c=c^{\prime}$. Therefore, below in Fig. 7 we present the diagonal
results for $dB_{c}/d\delta_{c^{\prime}}$ at $c\neq c^{\prime}$ by a separate
magenta color while all other countries sensitivity are characterized by color
bar.

For KEU9 the strongest sensitivity values $dB_{c}/d\delta_{c^{\prime}}$ are
obtained for countries: Czechia, Tunisia, Morocco with positive values
$0.027,0.014,0.013$ and S.Korea, Canada, Russia with negative values
$-0.074,-0.072,-0.062$. For USA case these are: Canada, Mexico, Venezuela with
positive values $0.086,0.063,0.008$ and S.Korea, Japan, United Kingdom with
negative values $-0.045,-0.044,-0.044$; for China we find: Korea, Congo,
Philippines with positive values $0.033,0.015,0.014$ and Mexico, Canada,
Poland with negative values $-0.090,-0.076,-0.074$. For Russia we obtain:
Uzbekistan, Belarus, Kazakhstan with positive values $0.028,0.026,0.020$ and
Norway, Sweden, United Kingdom with negative values $-0.019,-0.016,-0.013$.

This shows that the Google matrix analysis allows to determine on pure
mathematical grounds the mutual trade dependences of world countries.

### Network structure of trade from reduced Google matrix

We use the REGOMAX algorithm described above to obtain the reduced Google
matrix of trade flows between certain selected countries. We choose the case
of 4 countries which has a strong world trade influence KEU9, USA, China,
Russia with 10 trade products of Table 1. In this way the size of $G_{R}$
(Import or PageRank direct flow direction) and ${G^{*}}_{R}$ (Export or
CheiRank inverted flow direction) is equal to $40$. For clarity we show only 4
directed links corresponding to the most strong matrix elements from a given
node (only non-diagonal terms are shown). The obtained networks are shown in
Fig. 8 and Fig. 9 respectively.

The network structure shown in Fig. 8 from $G_{R}$ shows that the main
importing nodes are machinery ($s=7$) of KEU9 and USA. In a similar way the
main exporting nodes of ${G^{*}}_{R}$ are again machinery product of KEU9,
China and USA (Fig. 9). This clearly show the importance of machinery product
for the world trade.

## Discussion

Above we considered the trade influence of kernel EU 9 countries (KEU9)
considered as a one united state following the proposal pushed forward by
Christian Saint-Etienne in [2]. The analysis is done on the bases of
multiproduct trade data provided by UN COMTRADE [8]. Our results are based on
the advanced Google matrix analysis of multiproduct world trade network flows
for years 2012-2018 between all world countries registered at UN. They clearly
show that KEU9 takes the world leading position in PageRank and CheiRank
probabilities being ahead of USA and China. This mathematical network analysis
demonstrates that KEU9 becomes the main player in the international trade.
This provides additional mathematical foundation for the historical,
economical and political arguments presented in [2] in the favor of coherent
strong impact of KEU9 (if united) on the world development.

We also show that the Google matrix analysis allows to obtained significantly
deeper information about world trade comparing to the Import-Export analysis
usually used in economy studies.

## Acknowledgments

We thank Katia Jaffres-Runser (INP ENSEEIHT Toulouse) for useful discussions.
This research has been partially supported through the grant NANOX $N^{\circ}$
ANR-17-EURE-0009 (project MTDINA) in the frame of the Programme des
Investissements d’Avenir, France and in part by APR 2019 call of University of
Toulouse and by Region Occitanie (project GoIA). We thank UN COMTRADE for
providing us a friendly access to their detailed database.

## Supporting Information

Supporting Information file presents Figures S1-S6 extending information
presented in Figures of the main part to years 2012, 2014, 2016 and more
details.

## References

  * 1.
  * 2. C. Saint-Etienne  (2018)  Osons l’Europe des nations  Editions de l’Observatoire/Humensis
  * 3. Wikipedia:Economy of the European Union
https://en.wikipedia.org/wiki/Economy_of_the_European_Union (Retrieved
October, 2020)

  * 4. Largest Economies in the World
https://www.thebalance.com/world-s-largest-economy-3306044 (Retrieved October,
2020)

  * 5. Wikipedia: European Union
https://en.wikipedia.org/wiki/European_Union (Retrieved October, 2020)

  * 6. Wikipedia: Politics of the European Union
https://en.wikipedia.org/wiki/Politics_of_the_European_Union (Retrieved
October, 2020)

  * 7. CNEWS: Face a l’info avec Eric Zemmour Saison 1, 14 Nov 2019
https://www.cnews.fr/emission/2019-11-14/face-linfo-du-14112019-899420
(Retrieved October 2020)

  * 8. United Nations Commodity Trade Statistics Database
Available: http://comtrade.un.org/db/. (Retrieved October 2020).

  * 9. World Trade Organization (2018)  World Trade Statistical Review 2018
Available:
https://www.wto.org/english/res_e/statis_e/wts2018_e/wts18_toc_e.htm
(Retrieved October 2020).

  * 10. Krugman PR, Obstfeld M, Melitz M  (2011)  International economics: theory & policy  Prentic Hall, New Jersey
  * 11. Brin S, Page L  (1998)  The anatomy of a large-scale hypertextual Web search engine  Computer Networks and ISDN Systems 30, 107
  * 12. Langville AM, Meyer CD  (2006)  Google’s PageRank and beyond: the science of search engine rankings  Princeton University Press, Princeton
  * 13. Ermann L, Frahm KM, Shepelyansky DL  (2015)  Google matrix analysis of directed networks  Rev. Mod. Phys. 87, 1261
  * 14. Ermann L, Shepelyansky DL  (2011)  Google matrix of the world trade network  Acta Physica Polonica A 120, A158
  * 15. Ermann L, Shepelyansky DL  (2015)  Google matrix analysis of the multiproduct world trade network  Eur. Phys. J. B 88, 84
  * 16. Coquide C, Ermann L, Lages J, Shepelyansky DL  (2019)  Influence of petroleum and gas trade on EU economies from the reduced Google matrix analysis of UN COMTRADE data  Eur. Phys. J. B 92, 171
  * 17. Coquide C, Lages J, Shepelyansky DL  (2020)  Crisis contagion in the world trade network  Appl. Network Science (Springer) 5, 67
  * 18. Serrano MA, Boguna M, Vespignani A  (2007)  Patterns of dominant flows in the world trade web,  Journal of Economic Interaction and Coordination 2(2), 111
  * 19. Fagiolo G, Reyes J, Schiavo S  (2009)  World-trade web: Topological properties, dynamics, and evolution,  Phys. Rev. E 79, 036115
  * 20. He J, Deem MW  (2010)  Structure and response in the world trade network  Phys. Rev. Lett. 105, 198701
  * 21. Fagiolo G, Reyes J, Schiavo S  (2010)  The evolution of the world trade web: a weighted-network analysis  Journal of Evolutionary Economics 20, 479
  * 22. Barigozzi M, Fagiolo G, Garlaschelli D  (2010)  Multinetwork of international trade: a commodity-specific analysis  Phys. Rev. E 81, 046104
  * 23. De Benedictis L, Tajoli L  (2011)  The world trade network,  The World Economy 34(8), 1417
  * 24. (2014) Deguchi T, Takahashi K, Takayasu H, Takayasu M  Hubs and authorities in the world trade network using a weighted hits algorithm,  PLoS ONE 9(7), 1
  * 25. (2010) Chepelianskii AD  Towards physical laws for software architecture,  arXiv:1003.5455 [cs.SE]
  * 26. Zhirov AO, Zhirov OV, Shepelyansky DL  (2010)  Two-dimensional ranking of Wikipedia articles,  Eur. Phys. J. B 77, 523
  * 27. Frahm KM, Shepelyansky DL  (2016)  Reduced Google matrix,  arXiv:1602.02394[physics.soc]
  * 28. Frahm KM, Jaffres-Runser K, Shepelyansky DL,  (2016)  Wikipedia mining of hidden links between political leaders,  Eur. Phys. J. B 89, 269
  * 29. El Zant S, Jaffres-Runser K, Shepelyansky DL,  (2018)  Capturing the influence of geopolitical ties from Wikipedia with reduced Google matrix,  PLoS ONE 13(8), e0201397
  * 30. Coquide C, Lages J, Shepelyansky DL  (2019)  World influence and interactions of universities from Wikipedia networks,  Eur. Phys. J. B 92, 3
  * 31. Rollin G, Lages J, Serebriyskaya TS, Shepelyansky DL  (2019)  Interactions of pharmaceutical companies with world countries, cancers and rare diseases from Wikipedia network analysis,  PLoS ONE 14(12), e0225500
  * 32. Lages J, Shepelyansky DL, Zinovyev A  (2018)  Inferring hidden causal relations between pathway members using reduced Google matrix of directed biological networks  PLoS ONE 13(1), e0190812
  * 33. Zinovyev A, Czerwinska U, Cantini L, Barillot E, Frahm KM, Shepelyansky DL  (2020)  Collective intelligence defines biological functions in Wikipedia as communities in the hidden protein connection network,  PLoS Comput Biol 16(2), e1007652
  * 34. Kandiah V, Escaith H, Shepelyansky DL  (2015)  Google matrix of the world network of economic activities,  Eur. Phys. J. B 88, 186

Table 1: Product code and name for SITC classification level 1 Code | Name
---|---
0 | Food and live animals
1 | Beverages and tobacco
2 | Crude materials,inedible,except fuels
3 | Mineral fuels etc
4 | Animal and vegetable oils and fats
5 | Chemicals and related products,n.e.s.
6 | Basic manufactures
7 | Machinery,transport equipment
8 | Miscellaneous manufactured articles
9 | Goods not classified elsewhere
Table 2: Top 20 Ranking of PageRank ($K$), CheiRank ($K^{*}$), ImportRank and ExportRank for the year 2018. The table lists country ranks of both panels of Fig. 1. Rank | PageRank ($K$) | CheiRank ($K^{*}$) | ImportRank | ExportRank
---|---|---|---|---
1 | KEU9 | KEU9 | USA | China
2 | USA | China | KEU9 | KEU9
3 | China | USA | China | USA
4 | United Kingdom | Japan | Japan | Japan
5 | India | Republic of Korea | United Kingdom | Republic of Korea
6 | U. Arab Emirates | Russian Federation | Republic of Korea | United Kingdom
7 | Japan | India | India | Russian Federation
8 | Mexico | U. Arab Emirates | Canada | Mexico
9 | Republic of Korea | United Kingdom | Mexico | Canada
10 | Canada | Singapore | Singapore | Singapore
11 | Singapore | South Africa | Switzerland | India
12 | Russian Federation | Thailand | Poland | Switzerland
13 | Turkey | Malaysia | U. Arab Emirates | Malaysia
14 | Switzerland | Canada | Russian Federation | Australia
15 | Poland | Turkey | Thailand | U. Arab Emirates
16 | Australia | Mexico | Vietnam | Saudi Arabia
17 | South Africa | Australia | Australia | Thailand
18 | Saudi Arabia | Switzerland | Turkey | Vietnam
19 | Thailand | Poland | Malaysia | Brazil
20 | Brazil | Brazil | Czechia | Poland
Figure 1: Circles with county flags show positions of countries on the plane
of PageRank-CheiRank indexes $(K,K^{*})$ (summation is done over all products)
(left panel) and on the plane of ImportRank-ExportRank $\hat{K}$,
$\hat{K}^{*}$ from trade volume (right panel); data are shown only for index
values less than $61$ for year 2018; data for years 2012, 2014, 2016 are shown
in SupInf Fig.S1; KEU9 is marked by EU flag. Figure 2: World map of
probabilities balance $B_{c}=({P_{c}}^{*}-P_{c})/({P_{c}}^{*}+P_{c})$ of
countries. Top panel: probabilities are computed from the trade volume of
Export-Import. Bottom panel: probabilities are computed from PageRank and
CheiRank vectors. Data are shown for year 2018 (similar data for years 2012,
2014, 2016 are given in SupInf Fig.S2). Balance values are marked by color
with the corresponding color bar ($j$); countries absent in the UN COMTRADE
report are marked by black color (here and in other Figs.)  Figure 3:
Sensitivity of country balance $dB_{c}/d\delta_{s}$ for product $s=3$ (mineral
fuels). Top panel: probabilities are computed from the trade volume of Export-
Import. Bottom panel: probabilities are computed from PageRank and CheiRank
vectors. Sensitivity are marked by color with the corresponding color bar
($j$). Data are shown for year 2018 Figure 4: Same as in Fig. 3 but for
product $s=7$ (machinery). Data are shown for year 2018 (similar data for
years 2012, 2014, 2016 are given in SupInf Fig.S6). Figure 5: Sensitivity of
country balance $dB_{c}/d\delta_{cs}$ for product price $s=7$ (machinery) of
KEU9 (top row), USA (middle row), China (bottom row). Left column:
probabilities are computed from PageRank and CheiRank vectors. Right column:
probabilities are computed from the trade volume of Export-Import. Data are
shown for year 2018. Balance values are marked by color with the corresponding
color bar ($j$). Sensitivity values for KEU9, USA, China are
$dB_{c}/d\delta_{cs}=0.12,0.12,0.15$ (left column) and $0.21,0.18,0.24$ (right
column) respectively, these values are marked by separate magenta color to
highlight sensitivity of other countries in a better way. Figure 6:
Sensitivity of country balance $dB_{c}/d\delta_{cs}$ for product price $s=3$
(mineral fuels) of Russia. Left panel: probabilities are computed from
PageRank and CheiRank vectors; Right panel: probabilities are computed from
the trade volume of Export-Import. Data are shown for year 2018. Sensitivity
values for Russia are $dB_{c}/d\delta_{cs}=0.13$ (left panel) and $0.24$
(right panel); these values are marked by separate magenta color to highlight
sensitivity of other countries in a better way. Figure 7: Labor cost
sensitivity of country balance $dB_{c}/d\sigma_{c^{\prime}}$ for KEU9 (top
left panel), USA (top right panel), China (bottom left panel), Russia (bottom
right panel) for year 2018 (sensitivity is obtained from PageRank-CheiRank
probabilities). The diagonal term $dB_{c}/d\sigma_{c}$ is shown by a fixed
magenta color with numerical values being
$dB_{c}/d\sigma_{c}=0.30,,0.31,0.30,0.29$ respectively for KEU9, USA, China,
Russia (this allows to highlight the effect for other countries in a better
way)). Figure 8: Network trade structure between KEU9 (marked as EU), USA
(US), China (CN), Russia (RU), Japan (JP) with 10 products of Table 1 in year
2018. Network is obtained from the reduced Google matrix $G_{R}$ by tracing
four strongest outgoing links. Countries are shown by circles with two letters
of country and product index from Table 1. The arrow direction from node $A$
to node $B$ means that $B$ imports from $A$. All $40$ nodes are shown. Figure
9: Same in Fig. 8 but for the reduced Google matrix ${G^{*}}_{R}$
corresponding of inverted trade flows (CheiRank or Export direction). The
arrow direction from node $A$ to node $B$ means that $B$ exports to $A$.

SUPPORTING INFORMATION FOR:
World impact of kernel EU 9 countries
from Google matrix analysis of UN COMTRADE network

Justin Loye1,2, Leonardo Ermann3,4, Dima L. Shepelyansky1,∗

1 Laboratoire de Physique Théorique du CNRS, IRSAMC, Université de Toulouse,
CNRS, UPS, 31062 Toulouse, France

2 Institut de Recherche en Informatique de Toulouse, Université de Toulouse,
UPS, 31062 Toulouse, France

3 Departamento de Física Teórica, GIyA, Comisión Nacional de Energía Atómica.
Av. del Libertador 8250, 1429 Buenos Aires, Argentina

4 Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Buenos
Aires, Argentina

$\ast$ Webpage: www.quantware.ups-tlse.fr/dima

## Additional data

Here we present additional Figures related to the material presented in the
main part of the article.

Figure S1: Same as Fig.1 for years 2012, 2014, 2016 (from top to bottom).
Figure S2: Same as Fig.2, left column shows data from CheiRank and PageRank
vectors, right column shows data from Export and Import volumes; data are
given for years 2012, 2014, 2016 (from top to bottom). Figure S3: World map of
countries Import volume in USD for years 2012, 2014, 2016 (from top to
bottom); values lower then $10^{8}$ are clipped for better visibility. Figure
S4: World map of countries Export volume in USD for years 2012, 2014, 2016
(from top to bottom); values lower then $10^{8}$ are clipped for better
visibility. Figure S5: Same as Fig.3, left column shows data from CheiRank and
PageRank vectors, right column shows data from Export and Import volumes; data
are given for years 2012, 2014, 2016 (from top to bottom). Figure S6: Same as
Fig.4, left column shows data from CheiRank and PageRank vectors, right column
shows data from Export and Import volumes; data are given for years 2012,
2014, 2016 (from top to bottom).