# Fast Primal-Dual Update against Local Weight Update
in Linear Assignment Problem and Its Application

Kohei Morita Indeed Japan, Tokyo, Japan. Email<EMAIL_ADDRESS>Shinya
Shiroshita Preferred Networks, Tokyo, Japan. Email<EMAIL_ADDRESS>Yutaro Yamaguchi Osaka University, Osaka, Japan. Email:
<EMAIL_ADDRESS>Yu Yokoi Tokyo Institute of Technology,
Tokyo, Japan. Email<EMAIL_ADDRESS>
###### Abstract

We consider a dynamic situation in the weighted bipartite matching problem:
edge weights in the input graph are repeatedly updated and we are asked to
maintain an optimal matching at any moment. A trivial approach is to compute
an optimal matching from scratch each time an update occurs. In this paper, we
show that if each update occurs locally around a single vertex, then a single
execution of Dijkstra’s algorithm is sufficient to preserve optimality with
the aid of a dual solution. As an application of our result, we provide a
faster implementation of the envy-cycle procedure for finding an envy-free
allocation of indivisible items. Our algorithm runs in $\mathrm{O}(mn^{2})$
time, while the known bound of the original one is $\mathrm{O}(mn^{3})$, where
$n$ and $m$ denote the numbers of agents and items, respectively.

#### Keywords

Assignment problem, Bipartite matching, Envy-free allocation of indivisible
items (EF1/EFX).

##  1 Introduction

The weighted matching problem in bipartite graphs is a fundamental
combinatorial optimization problem. Initiated by the so-called Hungarian
method [11], a variety of efficient algorithms for the problem have been
developed [3, 6, 9, 10, 16].

We consider a dynamic situation in this problem: we are given an edge-weighted
bipartite graph in which edge weights are repeatedly updated, and asked to
maintain an optimal matching at any moment. While many papers on dynamic
matching consider fast approximation for the situation in which edges
arbitrarily appear and disappear (one by one), we focus on maintaining an
exact optimal solution when the underlying graph is fixed and edges change
their weights (possibly simultaneously).111In a similar context, Sankowski
[13] proposed pseudopolynomial-time exact maintenance of the maximum weight of
a matching in a dynamic bipartite graph. For the literature of usual dynamic
graph algorithms including those for matchings, see, e.g., [4, 7, 8, 15].

A trivial approach is to compute an optimal matching from scratch each time
the weights are updated. This simple method is applicable to any update; using
a typical efficient algorithm, it requires $\mathrm{O}(n(m+n\log n))$ time per
update, where $n$ and $m$ denote the numbers of vertices and edges,
respectively.

In this paper, we show that if each update occurs around a single vertex, a
more efficient algorithm can be designed. Specifically, in such a situation,
by maintaining a dual optimal solution simultaneously, it suffices to perform
Dijkstra’s algorithm only once for updating both primal and dual
solutions.222The readers familiar with the Hungarian method may think that the
idea is obvious and moreover it is a folklore. We, however, could not find any
description of this fact, and we believe that it is worth reporting with a
clearly described proof. This requires $\mathrm{O}(m+n\log n)$ time per
update.

We first describe our idea in the setting of the assignment problem, i.e., a
special case in which the input bipartite graph is balanced and we are asked
to find an optimal perfect matching. Suppose that the edge weights are updated
around a vertex $s$. On the primal side, we can obtain a new optimal matching
from an old one by exchanging the edges along a minimum-weight cycle
intersecting $s$ in the auxiliary graph. Such a cycle can be found by
computing shortest paths from $s$ in the modified auxiliary graph obtained by
flipping the matching edge incident to $s$. On the dual side, from the
shortest path length, we can obtain a new potential with respect to the new
optimal matching. Since the modified auxiliary graph has no negative cycles,
these primal/dual updates can be done in $\mathrm{O}(nm)$ time by the
Bellman–Ford algorithm. Our main observation is that when the weights are
updated around a single vertex, the old potential can be easily adjusted to
the modified auxiliary graph; this enables us to apply Dijkstra’s algorithm,
which requires $\mathrm{O}(m+n\log n)$ time.

We then discuss cases in which the input bipartite graph is unbalanced or the
perfectness of a matching is not imposed. Although such cases can be reduced
to the assignment problem by adding dummy vertices and edges, a naïve
reduction may essentially increase the input size. We show that by avoiding an
explicit reduction, the same computational time bound can be obtained.

We also demonstrate that our technique is useful for obtaining a faster
implementation of the envy-cycle procedure for finding an _EF1 (envy-free up
to one item)_ allocation of indivisible items [12]. We observe that the envy-
cycle procedure can be modified so that the part of cancelling envy cycles is
formulated as the assignment problem, whose input graph is updated around a
single vertex in each iteration. Our primal-dual update algorithm fits this
situation. While the well-known running time bound of the original envy-cycle
procedure is $\mathrm{O}(mn^{3})$, our variant runs in $\mathrm{O}(mn^{2})$
time, where $n$ and $m$ denote the numbers of agents and items, respectively.

The rest of the paper is organized as follows. In Section 2, we provide the
necessary definitions and basic facts on the assignment problem. In Section 3,
we describe our algorithm for updating primal/dual solutions in response to
weight update around a vertex in the assignment problem, and discuss
equivalent formulations of the weighted bipartite matching problem. In Section
4, we present an improved variant of the envy-cycle procedure.

##  2 Preliminaries

In the _(linear) assignment problem_ , we are given an edge-weighted bipartite
graph and asked to find a perfect matching of minimum total weight. We briefly
review the theory behind typical primal-dual approaches for this problem (see
[14] for more details). Throughout the paper, except for Section 3.2, we
assume that the input graph is simple and has a perfect matching.

The input is a pair $(G,w)$ of a bipartite graph $G=(V^{+},V^{-};E)$ with
vertex set $V$ partitioned into $V^{+}$ and $V^{-}$ and edge set $E$, where
$|V^{+}|=|V^{-}|=n$ and $|E|=m$, and edge weights $w\colon
E\to\mathbb{R}$.333It is often assumed that the input graph is a complete
bipartite graph. The present formulation reduces to this setting by adding all
the absent edges with sufficiently large weight. This reduction, however, does
not preserve the input size if the original graph is not dense, i.e.,
$m=\mathrm{o}(n^{2})$. In this sense, we strictly evaluate the computational
complexity. We also regard $G$ as a directed graph in which each edge is
directed from $V^{+}$ to $V^{-}$.

A _matching_ $\mu\subseteq E$ in $G$ is an edge set in which every vertex
appears as an endpoint at most once, and $\mu$ is called _perfect_ if all the
vertices appear exactly once. We also regard a perfect matching $\mu$ as a
bijection from $V^{+}$ to $V^{-}$ and vice versa. That is, we write $\mu(u)=v$
and $\mu^{-1}(v)=u$ if $e=(u,v)\in\mu$. We define the weight of a matching
$\mu$ as $w(\mu)=\sum_{e\in\mu}w(e)$.

Let $\mu$ be a matching in $G=(V,E)=(V^{+},V^{-};E)$. The _auxiliary graph_
for $\mu$, denoted by $G_{\mu}=(V,E_{\mu})$, is a directed graph that is
obtained from $G$ by flipping the direction of the edges in $\mu$, i.e.,

$\displaystyle E_{\mu}\coloneqq(E\setminus\mu)\cup\bar{\mu}=\\{\,e\mid e\in
E\setminus\mu\,\\}\cup\\{\,\bar{e}\mid e\in\mu\,\\},$

where $\bar{e}$ denotes the reverse edge of $e$. We also define the _auxiliary
weight_ $w_{\mu}\colon E_{\mu}\to\mathbb{R}$ by

$w_{\mu}(e)\coloneqq\begin{cases}w(e)&(e\in E\setminus\mu),\\\
-w(\bar{e})&(e\in\bar{\mu}).\end{cases}$

We define the weight of a path or a cycle in $G_{\mu}$, say $P$, as
$w_{\mu}(P)=\sum_{e\in P}w_{\mu}(e)$. We call a path/cycle _shortest_ if it is
of minimum weight (under some specified condition) and _negative_ if it is of
negative weight.

Let $V_{\mu}^{+}$ and $V_{\mu}^{-}$ denote the sets of vertices not matched by
$\mu$ in $V^{+}$ and in $V^{-}$, respectively. An _augmenting path_ in
$G_{\mu}$ is a directed path starting from a vertex in $V_{\mu}^{+}$ and
ending at a vertex in $V_{\mu}^{-}$. Let $P$ be an augmenting path or a cycle
in $G_{\mu}$. We then obtain another matching, denoted by $\mu\triangle P$,
from $\mu$ by exchanging the edges in $E\setminus\mu$ and in $\mu$ along $P$,
i.e., by adding the edges in $(E\setminus\mu)\cap P$ and by removing the edges
in $\mu\cap P$. Note that $w(\mu\triangle P)=w(\mu)+w_{\mu}(P)$ by definition.

A celebrated feature of a minimum-weight matching, which is utilized in many
algorithms, is described below.

######  Lemma 2.1 (cf. [14, Section 17.2]).

Let $\mu$ be a matching of size $k$.

  * •

If $\mu$ is a minimum-weight matching in $G$ subject to $|\mu|=k$, then
$(G_{\mu},w_{\mu})$ has no negative cycles; these conditions are equivalent
when $k=n$.

  * •

If $\mu$ is a minimum-weight matching in $G$ subject to $|\mu|=k$ and $P$ is a
shortest augmenting path in $(G_{\mu},w_{\mu})$, then
$\mu^{\prime}\coloneqq\mu\triangle P$ is a minimum-weight matching in $G$
subject to $|\mu^{\prime}|=k+1$.

By Lemma 2.1, a minimum-weight matching $\mu$ of size $k$ can be augmented to
a minimum-weight matching $\mu^{\prime}$ of size $k+1$ by finding a shortest
path in a directed graph with no negative cycles. This is done in
$\mathrm{O}(nm)$ time by the Bellman–Ford algorithm. Furthermore, this is done
in $\mathrm{O}(m+n\log n)$ time by Dijkstra’s algorithm (with the aid of the
Fibonacci heap) if we simultaneously maintain a potential (i.e., a dual
solution), which is defined as follows.

For a matching $\mu$ and a function $p\colon V\to\mathbb{R}$, we define the
_reduced weight_ $w_{\mu,p}\colon E_{\mu}\to\mathbb{R}$ by

$w_{\mu,p}(e)\coloneqq w_{\mu}(e)+p(u)-p(v)\quad(e=(u,v)\in E_{\mu}).$

We say that $p$ is a _potential for $(G_{\mu},w_{\mu})$_ if $w_{\mu,p}(e)\geq
0$ for every $e\in E_{\mu}$.

We here summarize the Hungarian method with this speed-up. First, we set
$\mu_{0}\leftarrow\emptyset$ and

$p_{0}(v)\leftarrow\begin{cases}0&(v\in V^{+}),\\\ \min\\{\,w(e)\mid
e\in\delta_{G}(v)\,\\}&(v\in V^{-}),\end{cases}$

where $\delta_{G}(v)$ denotes the set of edges incident to $v$ in $G$. Then,
for each $k=0,1,\dots,n-1$, we successively compute a minimum-weight perfect
matching $\mu_{k+1}$ of size $k+1$ and a potential $p_{k+1}$ for
$(G_{\mu_{k+1}},w_{\mu_{k+1}})$ using $\mu_{k}$ and $p_{k}$ as follows.

In each iteration, we compute a shortest path from $V_{\mu_{k}}^{+}$ to each
vertex $v\in V$ in $(G_{\mu_{k}},w_{\mu_{k}})$ by Dijkstra’s algorithm with
the aid of the potential $p_{k}$. On the primal side, let $P$ be a shortest
augmenting path, and set $\mu_{k+1}\leftarrow\mu_{k}\triangle P$. On the dual
side, for each vertex $v\in V$, let $d(v)$ denote the weight of a shortest
path from $V_{\mu_{k}}^{+}$ to $v$, and set $p_{k+1}(v)\leftarrow d(v)$ (where
we set $p_{k+1}(v)\leftarrow p_{k}(v)$ if $v$ is not reachable from
$V_{\mu_{k}}^{+}$).

The bottleneck is executing Dijkstra’s algorithm $n$ times; hence, the
computational time of the whole algorithm is bounded by $\mathrm{O}(n(m+n\log
n))$.

##  3 Primal/Dual Update against Weight Update around a Vertex

In this section, we propose an $\mathrm{O}(m+n\log n)$-time algorithm for
updating an optimal matching and a corresponding potential when the edge
weights are updated around a single vertex. We show this for the assignment
problem in Section 3.1, and discuss other polynomially equivalent formulations
of the weighted bipartite matching problem in Section 3.2.

###  3.1 Assignment Problem Case

The problem under consideration is described below. Suppose we have a pair of
primal and dual optimal solutions for an instance $(G,w)$ of the assignment
problem, i.e., we have a minimum-weight perfect matching $\mu$ and a potential
$p$ for $(G_{\mu},w_{\mu})$. Now, the edge weights are updated arbitrarily but
only around a single vertex $s\in V$, and we are asked to find a pair of
primal and dual optimal solutions for the new instance $(G,{\tilde{w}})$. By
symmetry, we can assume $s\in V^{+}$. In addition, for the sake of simplicity,
we assume that $\delta_{G}(s)=\\{s\\}\times V^{-}$ by adding at most $n$
absent edges with sufficiently large weights.

Our algorithm is shown in Algorithm 1, which works as follows. Let
$e^{\prime}=(s,t)$ be the edge in $\mu$ incident to $s$ and let
$\mu^{\prime}\coloneqq\mu\setminus\\{e^{\prime}\\}$. Then, in
$G_{\mu^{\prime}}$, all the edges in $\delta_{G}(s)=\\{s\\}\times V^{-}$ are
directed from $s$ to $V^{-}$. Our objective is to find shortest paths from $s$
in $G_{\mu^{\prime}}$ with respect to the new weights ${\tilde{w}}$. Note that
$p$ is a potential for $(G_{\mu},w_{\mu})$ (i.e., $w_{\mu,p}(e)\geq 0$ for
every $e\in E_{\mu}$), but it may not be a potential for the current graph
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$. However, since
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$ differs from
$(G_{\mu},w_{\mu})$ only around $s$, we can recover the nonnegativity
condition only by modifying the value of $p$ on $s$. Specifically, we define
$p^{\prime}$ as $p^{\prime}(s)\coloneqq-\min\\{\,{\tilde{w}}(e)-p(v)\mid
e=(s,v)\in\delta_{G}(s)\,\\}$ and $p^{\prime}(v)\coloneqq p(v)$ for every
$v\in V\setminus\\{s\\}$. Then, $p^{\prime}$ is indeed a potential for
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$, and we can efficiently
compute a shortest path from $s$ to each vertex by Dijkstra’s algorithm. Using
the obtained paths and their weights, we update the primal and dual solutions.

Algorithm 1 Primal Dual Update

Input:

A bipartite graph $G=(V^{+},V^{-};E)$, edge weights $w\colon E\to\mathbb{R}$,
a minimum-weight perfect matching $\mu$ in $(G,w)$, a potential $p$ for
$(G_{\mu},w_{\mu})$, a vertex $s\in V^{+}$ with $\delta_{G}(s)=\\{s\\}\times
V^{-}$, and new edge weights ${\tilde{w}}\colon E\to\mathbb{R}$ such that
${\tilde{w}}(e)=w(e)$ for every edge $e\in E\setminus\delta_{G}(s)$.

Output:

A minimum-weight perfect matching ${\tilde{\mu}}$ in $(G,{\tilde{w}})$ and a
potential ${\tilde{p}}$ for $(G_{{\tilde{\mu}}},{\tilde{w}}_{{\tilde{\mu}}})$.

  1. 1.

Let $e^{\prime}=(s,t)$ be the edge in $\mu$ incident to $s$ and let
$\mu^{\prime}=\mu\setminus\\{e^{\prime}\\}$.

  2. 2.

Set $p^{\prime}(s)=-\min\\{\,{\tilde{w}}(e)-p(v)\mid
e=(s,v)\in\delta_{G}(s)\,\\}$ and $p^{\prime}(v)=p(v)$ for every $v\in
V\setminus\\{s\\}$. Compute a shortest path $P_{v}$ from $s$ to each $v\in V$
in $(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$ by Dijkstra’s algorithm
with the aid of potential $p^{\prime}$.

  3. 3.

Set ${\tilde{\mu}}\leftarrow\mu^{\prime}\triangle P_{t}$. For each $v\in V$,
set ${\tilde{p}}(v)\leftarrow\tilde{w}_{\mu^{\prime}}(P_{v})$.

  4. 4.

Return ${\tilde{\mu}}$ and ${\tilde{p}}$.

######  Theorem 3.1.

Algorithm 1 correctly finds ${\tilde{\mu}}$ and ${\tilde{p}}$ in
$\mathrm{O}(m+n\log n)$ time.

###### Proof.

The algorithm clearly runs in $\mathrm{O}(m+n\log n)$ time. We prove its
correctness.

We first confirm that $p^{\prime}$ in Step 2 is indeed a potential for
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$. Since $p$ is a potential for
$(G_{\mu},w_{\mu})$, we have
${\tilde{w}}_{\mu^{\prime},p^{\prime}}(e)=w_{\mu,p}(e)\geq 0$ for every edge
$e\in E_{\mu^{\prime}}\setminus\delta_{G}(s)\subseteq E_{\mu}$. Also, as every
edge $e=(s,v)\in\delta_{G}(s)$ is directed from $s$ to $v$ in
$G_{\mu^{\prime}}$, we have

$\displaystyle{\tilde{w}}_{\mu^{\prime},p^{\prime}}(e)$
$\displaystyle={\tilde{w}}(e)+p^{\prime}(s)-p^{\prime}(v)$
$\displaystyle={\tilde{w}}(e)-p(v)-\min\left\\{\,{\tilde{w}}(e^{*})-p(v^{*})\mid
e^{*}=(s,v^{*})\in\delta_{G}(s)\,\right\\}$ $\displaystyle\geq 0.$

We next show the correctness of the primal update. Since $e^{\prime}=(s,t)$
itself forms an $s$–$t$ path in $G_{\mu^{\prime}}$, the weight of $P_{t}$ is
at most ${\tilde{w}}(e^{\prime})$; hence,
${\tilde{w}}({\tilde{\mu}})\leq{\tilde{w}}(\mu^{\prime})+{\tilde{w}}(e^{\prime})={\tilde{w}}(\mu)$.
Suppose to the contrary that there is a perfect matching $\mu^{*}$ in $G$ with
${\tilde{w}}(\mu^{*})<{\tilde{w}}({\tilde{\mu}})\leq{\tilde{w}}(\mu)$. Then,
$\mu^{*}\triangle\mu$ contains a negative cycle $C$ in
$(G_{\mu},{\tilde{w}}_{\mu})$ (where each edge in $\mu$ is flipped). Since
$(G_{\mu},w_{\mu})$ has no negative cycles (by Lemma 2.1) and ${\tilde{w}}$
differs from $w$ only around $s$, the cycle $C$ must intersect $s$. This
implies that $C$ is the only negative cycle in $\mu^{*}\triangle\mu$ and hence
${\tilde{w}}_{\mu}(C)\leq{\tilde{w}}(\mu^{*})-{\tilde{w}}(\mu)$. In addition,
as $\bar{e}^{\prime}$ is the only edge entering $s$ in $G_{\mu}$, it is
contained in $C$. By removing $\bar{e}^{\prime}$ from $C$, we obtain an
$s$–$t$ path $P^{*}_{t}$ in $(G_{\mu},{\tilde{w}}_{\mu^{\prime}})$ such that
${\tilde{w}}_{\mu^{\prime}}(P^{*}_{t})={\tilde{w}}_{\mu}(C)-{\tilde{w}}_{\mu}(\bar{e}^{\prime})\leq{\tilde{w}}(\mu^{*})-{\tilde{w}}(\mu)+{\tilde{w}}(e^{\prime})<{\tilde{w}}({\tilde{\mu}})-{\tilde{w}}(\mu^{\prime})={\tilde{w}}_{\mu^{\prime}}(P_{t})$,
which contradicts the fact that $P_{t}$ is a shortest path.

Finally, we show that ${\tilde{p}}$ is a potential for
$(G_{{\tilde{\mu}}},{\tilde{w}}_{{\tilde{\mu}}})$. Since ${\tilde{p}}(v)$ is
the weight of the shortest $s$–$v$ path $P_{v}$ in
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$ for any $v\in V$, the function
${\tilde{p}}$ is clearly a potential for
$(G_{\mu^{\prime}},{\tilde{w}}_{\mu^{\prime}})$. That is, we have
${\tilde{w}}_{\mu^{\prime},{\tilde{p}}}(e)={\tilde{w}}_{\mu^{\prime}}(e)+{\tilde{p}}(u)-{\tilde{p}}(v)\geq
0$ for each $e=(u,v)\in E_{\mu^{\prime}}$. In particular, the equality holds
for every edge $e$ used in the shortest $s$–$v$ path $P_{v}$ for some $v\in
V$. Note that $G_{{\tilde{\mu}}}$ is obtained from $G_{\mu^{\prime}}$ by
flipping the direction of the edges on $P_{t}$, whose
${\tilde{w}}_{\mu^{\prime},{\tilde{p}}}$ values are $0$. Therefore,
${\tilde{w}}_{{\tilde{\mu}},{\tilde{p}}}(e)\geq 0$ holds for every edge $e\in
E_{{\tilde{\mu}}}$, and hence ${\tilde{p}}$ is a potential for
$(G_{{\tilde{\mu}}},{\tilde{w}}_{{\tilde{\mu}}})$. ∎

Using Algorithm 1, we can design a simple $\mathrm{O}(n(m+n\log n))$-time
algorithm for the assignment problem as follows. We first set
$w^{\prime}(e)\leftarrow 0$ for each $e\in E$, find any perfect matching $\mu$
in $G$, and initialize $p(v)\leftarrow 0$ for each $v\in V$. This part is done
in $\mathrm{O}(nm)$ time by a naïve augmenting path algorithm (or faster by
sophisticated algorithms, but this is not a bottleneck). Clearly, $\mu$ is a
minimum-weight perfect matching in $(G,w^{\prime})$ and $p$ is a potential for
$(G_{\mu},w^{\prime}_{\mu})$. We then repeatedly make $w^{\prime}$ close to
the actual input weights $w$ by updating it around a single vertex $s\in
V^{+}$ in any order.

In each iteration, we set

${\tilde{w}}(e)\leftarrow\begin{cases}w(e)&(e\in\delta_{G}(s)),\\\
w^{\prime}(e)&(\text{otherwise}),\\\ \end{cases}$

and find a minimum-weight perfect matching ${\tilde{\mu}}$ in
$(G,{\tilde{w}})$ and a potential ${\tilde{p}}$ for
$(G_{\tilde{\mu}},{\tilde{w}}_{\tilde{\mu}})$ by applying Algorithm 1 to
$(G,w^{\prime},\mu,p,s,{\tilde{w}})$, where we temporarily add the absent
edges around $s$ with sufficiently large weights if necessary. We then update
$w^{\prime}\leftarrow{\tilde{w}}$, $\mu\leftarrow{\tilde{\mu}}$, and
$p\leftarrow{\tilde{p}}$.

After the iterations are completed, we have $w^{\prime}=w$. The correctness
and the computational time bound directly follow from Theorem 3.1.

###  3.2 Other Formulations

In this section, we discuss other polynomially equivalent formulations of our
problem so that we can deal with various settings. Consider the situation in
which we are required to find a maximum-weight (not necessarily perfect)
matching in an edge-weighted (not necessarily balanced) bipartite graph.
First, by negating all the edge weights, the objective turns into minimizing
the total weight of a matching.

Let $G=(V^{+},V^{-};E)$ be a bipartite graph with $|V^{+}|\geq|V^{-}|>0$, and
let $w\colon E\to\mathbb{R}$. Without loss of generality, we assume that $w$
is nonpositive, since any positive-weight edge can be excluded from an optimal
matching. We say that a matching $\mu$ in $G$ is _right-perfect_ if all the
vertices in $V^{-}$ appear exactly once. By adding a vertex $v^{\prime}$ to
$V^{+}$ with an edge $(v^{\prime},v)$ with weight $0$ for each $v\in V^{-}$,
any instance can be transformed into an instance having a right-perfect
optimal matching. Note that the input size does not essentially increase. In
the following, let $(G,w)$ be the graph after this transformation, so we have
$|V^{+}|>|V^{-}|$. Let $n=|V^{+}|$ and $m=|E|$.

We now consider the following analogous problem: we have a minimum-weight
right-perfect matching $\mu$ in an edge-weighted unbalanced bipartite graph
$(G,w)$ and a potential $p$ for $(G_{\mu},w_{\mu})$, the edge weights $w$
change arbitrarily but only around a single vertex $s\in V$, and we are asked
to find a minimum-weight right-perfect matching ${\tilde{\mu}}$ in the new
graph $(G,{\tilde{w}})$ and a potential ${\tilde{p}}$ for
$(G_{\tilde{\mu}},{\tilde{w}}_{\tilde{\mu}})$.

By adding $|V^{+}|-|V^{-}|$ dummy vertices to $V^{-}$ that are adjacent to all
the vertices in $V^{+}$ by edges with weight $0$, this problem reduces to the
problem solved in Section 3.1. Let $(G^{\prime},w^{\prime})$ be the edge-
weighted balanced bipartite graph obtained by this reduction, and let
$V^{\prime}$ and $E^{\prime}$ denote its vertex set and edge set,
respectively. We then have $|V^{\prime}|=2n=\mathrm{O}(n)$, but $|E^{\prime}|$
may not be bounded by $\mathrm{O}(m+n\log n)$ if $G$ is sparse and far from
balanced. Hence, this naïve reduction may worsen the computational time bound
in Theorem 3.1.

We show that this problem can be solved in $\mathrm{O}(m+n\log n)$ time by
emulating Algorithm 1 for $(G^{\prime},w^{\prime})$ with avoiding its explicit
construction. Let $\mu^{\prime}$ and $p^{\prime}$ be the corresponding optimal
solutions and $(G^{\prime},{\tilde{w}}^{\prime})$ be the corresponding new
graph in the reduced problem.

######  Lemma 3.2.

If $(G^{\prime}_{\mu^{\prime}},{\tilde{w}}^{\prime}_{\mu^{\prime}})$ contains
a negative cycle, then there exists a shortest cycle that intersects at most
one dummy vertex.

###### Proof.

Let $C$ be a shortest cycle in
$(G^{\prime}_{\mu^{\prime}},{\tilde{w}}^{\prime}_{\mu^{\prime}})$ that
consists of as few edges as possible. We show that $C$ intersects at most one
dummy vertex. Suppose to the contrary that $C$ intersects two distinct dummy
vertices $v_{1}$ and $v_{2}$ in the extended $V^{-}$. Then, $C$ is divided
into two paths, $P_{1}$ from $v_{1}$ to $v_{2}$ and $P_{2}$ from $v_{2}$ to
$v_{1}$. By Lemma 2.1, the negative cycle $C$ must intersect $s$; as $s$ is
not a dummy vertex, we assume that $P_{1}$ intersects $s$ and $P_{2}$ does not
without loss of generality.

Let $u\in V^{+}$ be the vertex just before $v_{1}$ on $P_{2}$. If
${\mu^{\prime}}(u)=v_{2}$, then
${\tilde{w}}^{\prime}_{\mu^{\prime}}(P_{2})={\tilde{w}}(u,v_{2})+{\tilde{w}}(u,v_{1})=0$
by definition. Otherwise, $G^{\prime}_{\mu^{\prime}}$ contains a directed edge
$e=(u,v_{2})$, and one can obtain a cycle $C^{\prime}$ from $P_{2}$ by
replacing the last edge $(u,v_{1})$ with $e$. This $C^{\prime}$ does not
contain $s$, and hence
${\tilde{w}}^{\prime}_{\mu^{\prime}}(P_{2})={\tilde{w}}^{\prime}_{\mu^{\prime}}(C^{\prime})\geq
0$ by Lemma 2.1. Thus, in either case, we have
${\tilde{w}}^{\prime}_{\mu^{\prime}}(P_{2})\geq 0$, which implies
${\tilde{w}}^{\prime}_{\mu^{\prime}}(P_{1})<0$.

Let $u^{\prime}\in V^{+}$ be the vertex just before $v_{2}$ on $P_{1}$. If
${\mu^{\prime}}(u^{\prime})=v_{1}$, then
${\tilde{w}}^{\prime}_{\mu^{\prime}}(P_{1})={\tilde{w}}(u^{\prime},v_{1})+{\tilde{w}}(u^{\prime},v_{2})=0$,
a contradiction. Thus, there exists a directed edge
$e^{\prime}=(u^{\prime},v_{1})$ in $G^{\prime}_{\mu^{\prime}}$, and a cycle
$C^{\prime\prime}$ with
${\tilde{w}}^{\prime}_{\mu^{\prime}}(C^{\prime\prime})\leq{\tilde{w}}^{\prime}_{\mu^{\prime}}(C)$
and $|C^{\prime\prime}|<|C|$ can be obtained from $C$ by replacing the subpath
consisting of $(u^{\prime},v_{2})$ and $P_{2}$ with $e^{\prime}$. This
contradicts the choice of $C$; thus, we are done. ∎

######  Lemma 3.3.

For any vertex $v\in V$, there exists a shortest $s$–$v$ path in
$(G^{\prime}_{\mu^{\prime}},{\tilde{w}}^{\prime}_{\mu^{\prime}})$ that
intersects at most one dummy vertex. Moreover, for all the dummy vertices, the
weights of the shortest paths are the same.

###### Proof.

The first claim is shown by almost the same argument as the proof of Lemma
3.2. Pick a shortest $s$–$v$ path $P$ in
$(G^{\prime}_{\mu^{\prime}},{\tilde{w}}^{\prime}_{\mu^{\prime}})$ that
consists of as few edges as possible. If $P$ intersects two distinct dummy
vertices $v_{1}$ and $v_{2}$, then the subpath between $v_{1}$ and $v_{2}$ can
be skipped without increasing the weight, which contradicts the minimality of
$P$.

To see the second claim, for two distinct dummy vertices $v_{1}$ and $v_{2}$,
let $u_{1}$, $u_{2}$ be the vertices such that
$(u_{1},v_{1}),(u_{2},v_{2})\in\mu^{\prime}$. Then, by definition,
$(G^{\prime}_{\mu^{\prime}},{\tilde{w}}^{\prime}_{\mu^{\prime}})$ has four
edges $(v_{1},u_{1})$, $(v_{2},u_{2})$, $(u_{1},v_{2})$, and $(u_{2},v_{1})$
of weight $0$, and hence we have $d(v_{2})\leq d(v_{1})$ and $d(v_{1})\leq
d(v_{2})$. ∎

Lemmas 3.2 and 3.3 imply that, in order to emulate Algorithm 1 for
$(G^{\prime},w^{\prime})$, it suffices to introduce only one dummy vertex $x$
that is adjacent to all the vertices in $V^{+}$ with edges of weight $0$ and
to modify the perfect matching constraint so that $x$ appears exactly
$|V^{+}|-|V^{-}|$ times. Thus, we have completed the proof.

##  4 An Improved Variant of Envy-Cycle Procedure

In this section, as an application of our algorithm, we present an improved
variant of the _envy-cycle procedure_ for finding an envy-free allocation of
indivisible items. We describe the problem setting and the envy-cycle
procedure in Section 4.1, and demonstrate how to improve its computational
time via the assignment problem in Section 4.2.

###  4.1 Envy-Free Allocation of Indivisible Items

Consider a situation in which $n$ agents share $m$ indivisible items. Let
$N=[n]\coloneqq\\{1,2,\dots,n\\}$ and $M=[m]$ be the sets of agents and items,
respectively. Each item subset $X\subseteq M$ is called a _bundle_. Each agent
$i\in N$ has a _valuation_ $v_{i}\colon 2^{M}\to\mathbb{R}_{\geq 0}$ meaning
that agent $i$ evaluates a bundle $X\subseteq M$ by a nonnegative real
$v_{i}(X)$. We assume that $v_{i}$ is _monotone_ , i.e., $X\subseteq
Y\subseteq M\implies v_{i}(X)\leq v_{i}(Y)$. In addition, $v_{i}$ is said to
be _additive_ if $v_{i}(X)=\sum_{k\in X}v_{i,k}$ for every bundle $X\subseteq
M$, where $v_{i,k}\coloneqq v_{i}(\\{k\\})$ $(k\in M)$.

An _allocation_ is an $N$-indexed subpartition $(X_{i})_{i\in N}$ of $M$,
i.e., $X_{i}\subseteq M$ for each $i\in N$ and $X_{i}\cap X_{j}=\emptyset$ for
distinct $i,j\in N$. In an allocation $(X_{i})_{i\in N}$, an agent $i\in N$
_envies_ another agent $j\in N$ if $v_{i}(X_{i})<v_{i}(X_{j})$. Since each
item is indivisible, it is usually impossible to achieve a completely envy-
free allocation of the whole item set $M$. In the following, we describe two
well-studied relaxations of envy-freeness [5, 12].

######  Definition 4.1.

Let $\mathcal{X}=(X_{i})_{i\in N}$ be an allocation.

  * •

$\mathcal{X}$ is _envy-free up to one item_ (_EF1_ for short) if, for every
pair of agents $i,j\in N$, we have $v_{i}(X_{i})\geq v_{i}(X_{j})$ or
$v_{i}(X_{i})\geq v_{i}(X_{j}\setminus\\{k\\})$ for some $k\in X_{j}$.

  * •

$\mathcal{X}$ is _envy-free up to any item_ (_EFX_ for short) if, for every
pair of agents $i,j\in N$, we have $v_{i}(X_{i})\geq
v_{i}(X_{j}\setminus\\{k\\})$ for every $k\in X_{j}$.

Clearly, an EFX allocation is EF1. It is known that an EF1 allocation always
exists and can be found by a rather simple algorithm, the so-called _envy-
cycle procedure_ [12], whereas the existence of an EFX allocation is known for
very restricted situations and is widely open in general (see, e.g., [1]).

The envy-cycle procedure is summarized as follows. We initialize
$\mathcal{X}\leftarrow(X_{i})_{i\in N}$ with $X_{i}=\emptyset$ for each $i\in
N$, and distribute the items one-by-one in any order. In each iteration, we
consider the _envy graph_ $H(\mathcal{X})=(N,F(\mathcal{X}))$ for the current
allocation $\mathcal{X}$, which is defined by

$F(\mathcal{X})\coloneqq\\{\,(i,j)\mid i,j\in
N,~{}v_{i}(X_{i})<v_{i}(X_{j})\,\\}.$

While $H(\mathcal{X})$ contains a cycle, we update $\mathcal{X}$ by exchanging
the bundles along it; that is, we pick a cycle $C$ and update $X_{i}\leftarrow
X_{j}$ for all the edges $(i,j)\in C$ simultaneously.444This exchange is
usually repeated until some agent has no incoming edges, and it suffices. We
employ the present form to make its behavior analogous to the improved variant
of the envy-cycle procedure described in the next section. Subsequently, as
$H(\mathcal{X})$ has no cycles, there exists an agent $i\in N$ having no
incoming edges. To such an agent $i$, we distribute an unallocated item, say
$k\in M\setminus\bigcup_{j\in N}X_{j}$; that is, we set $X_{i}\leftarrow
X_{i}\cup\\{k\\}$.

It is easy to observe that at any point of this algorithm, $\mathcal{X}$ is an
EF1 allocation of the set of items distributed so far. The computational time
is bounded by counting the number of edges in $H(\mathcal{X})$, which can
increase only when an item is distributed and must decrease when the bundles
are exchanged along a cycle.

######  Theorem 4.2 (Lipton et al. [12, Theorem 2.1]).

The envy-cycle procedure correctly finds an EF1 allocation in
$\mathrm{O}(mn^{3})$ time.

In [2], it was pointed out that if the valuations are additive and the values
of items are in the same order for all the agents, this algorithm can be
employed to find an EFX allocation.

######  Theorem 4.3 (Barman and Krishnamurthy [2, Lemma 3.5]).

Suppose that the valuations of all agents are additive, and $v_{i,1}\geq
v_{i,2}\geq\dots\geq v_{i,m}$ for every agent $i\in N$. Then, if the items are
distributed in ascending order of $k=1,2,\dots,m$, the envy-cycle procedure
correctly finds an EFX allocation in $\mathrm{O}(mn^{3})$ time.

######  Remark 4.4.

When the valuations are additive, an EF1 allocation can be computed much
faster by the _round-robin protocol_ [5]: based on any fixed-in-advance order,
the agents select the most preferable remaining item one by one until no item
is left. This procedure can be implemented to run in $\mathrm{O}(nm\log m)$
time, where the bottleneck is sorting the items for each agent.

This procedure, however, cannot be easily extended to nonadditive valuations.
If the agents select items one by one in a round-robin fashion so that each
agent selects an item achieving the maximum increase of the value of his/her
bundle, then the resultant allocation is not necessarily EF1; such an instance
is easily constructed, e.g., by considering the situation when there are
complementary items (i.e., the valuations are superadditive for some pair of
items).

###  4.2 Improvement via Assignment Problem

First, we slightly modify the definition of an allocation by regarding it as a
combination of a “division” and an “assignment.” Hereafter, an _allocation_ is
a pair $(\mathcal{X},\mu)$ of an $N$-indexed subpartition
$\mathcal{X}=(X_{i})_{i\in N}$ of $M$ and a bijection $\mu\colon
N\to\mathcal{X}$, which means that bundle $\mu(i)\in\mathcal{X}$ is assigned
to agent $i\in N$. We say that $(\mathcal{X},\mu)$ is _EF1_ or _EFX_ if this
is true for subpartition $(X^{\prime}_{i})_{i\in N}$ with
$X^{\prime}_{i}=\mu(i)$ $(i\in N)$ in the sense of the original definition
(see Definition 4.1). The _welfare_ of an allocation $(\mathcal{X},\mu)$ is
defined as $W(\mathcal{X},\mu)=\sum_{i\in N}v_{i}(\mu(i))$.

######  Definition 4.5.

For an allocation $(\mathcal{X},\mu)$ with the bundle set
$\mathcal{X}=(X_{i})_{i\in N}$, the _assign-envy graph_
$G(\mathcal{X},\mu)=(N,\mathcal{X};E(\mathcal{X},\mu))$ is a bipartite graph
with edge weights $w\colon E(\mathcal{X},\mu)\to{\mathbb{R}}$ defined by

$\displaystyle E(\mathcal{X},\mu)$ $\displaystyle\coloneqq\\{\,(i,\mu(i))\mid
i\in N\,\\}\cup\\{\,(i,X_{j})\mid i\in
N,~{}X_{j}\in\mathcal{X},~{}v_{i}(\mu(i))<v_{i}(X_{j})\,\\},$ $\displaystyle
w(e)$ $\displaystyle\coloneqq-v_{i}(X_{j})\quad(e=(i,j)\in
E(\mathcal{X},\mu)).$

Note that the bijection $\mu\colon N\to\mathcal{X}$ is regarded as a perfect
matching in $G(\mathcal{X},\mu)$.555The negative sign appearing in the
definition of edge weights is used for a clear correspondence to the
assignment problem, in which the objective is to find a minimum-weight perfect
matching whereas in the current problem an allocation with a large welfare is
preferable. Also note that if we contract each edge $(i,\mu(i))$ in
$G(\mathcal{X},\mu)$, the resulting graph will coincide with the envy graph
$H(\mathcal{X}^{\prime})$ for $\mathcal{X}^{\prime}=(X^{\prime}_{i})_{i\in N}$
with $X^{\prime}_{i}=\mu(i)$ $(i\in N)$, which has been defined in Section
4.1. Thus, a cycle in $H(\mathcal{X}^{\prime})$ corresponds to a cycle in
$G(\mathcal{X},\mu)_{\mu}$ and vice versa.

######  Lemma 4.6.

Let $(G,w)$ be the assign-envy graph for an allocation $(\mathcal{X},\mu)$,
and $(G_{\mu},w_{\mu})$ be the auxiliary weighted graph for $\mu$. If
$G_{\mu}$ has a cycle $C$, then the weight $w_{\mu}(C)$ is negative.

###### Proof.

If an allocation is updated by exchanging bundles along a cycle in the envy
graph, each agent in the cycle receives a more valuable bundle; hence, the
welfare strictly increases. This is also true for the assign-envy graph. By
the definitions of $w$ and $(G_{\mu},w_{\mu})$, the weight $w_{\mu}(C)$ is
equal to $W(\mathcal{X},\mu)-W(\mathcal{X},\mu^{\prime})<0$, where
$\mu^{\prime}$ is the assignment after the exchange along $C$. ∎

######  Lemma 4.7.

Let $(G,w)$ be the assign-envy graph for an allocation $(\mathcal{X},\mu)$,
and $\mu^{*}$ be a minimum-weight perfect matching in $(G,w)$. Let
$(G^{*},w^{*})\coloneqq G(\mathcal{X},\mu^{*})$ and let
$(G^{*}_{\mu^{*}},w^{*}_{\mu^{*}})$ be the auxiliary weighted graph for
$\mu^{*}$. Then, the directed graph $G^{*}_{\mu^{*}}$ has no cycles.

###### Proof.

Suppose to the contrary that $G^{*}_{\mu^{*}}$ has a cycle $C$. We then have
$w^{*}_{\mu^{*}}(C)<0$ by Lemma 4.6. Let $\mu^{\prime}$ be the assignment
obtained from $\mu^{*}$ by the exchange along $C$. We then have
$W(\mathcal{X},\mu^{\prime})>W(\mathcal{X},\mu^{*})$ and
$v_{i}(\mu^{\prime}(i))>v_{i}(\mu^{*}(i))\geq v_{i}(\mu(i))$ for every $i\in
N$ with $\mu^{\prime}(i)\neq\mu^{*}(i)$. The latter implies that all the edges
in $\mu^{\prime}$ exist in $G$, and then the former contradicts that $\mu^{*}$
is a minimum-weight perfect matching in $(G,w)$. ∎

By Lemma 4.7, instead of repeatedly exchanging the bundles along cycles in the
envy graph, it suffices to find a minimum-weight perfect matching in the
assign-envy graph once in each iteration of the envy-cycle procedure. Note
that if $(\mathcal{X},\mu)$ is an EF1 allocation, then $(\mathcal{X},\mu^{*})$
is also EF1 since $v_{i}(\mu^{*}(i))\geq v_{i}(\mu(i))$ for every $i\in N$ by
the definition of the assign-envy graph.666Even if we start with the same EF1
allocation, the EF1 allocation obtained from a minimum-weight perfect matching
may be different from that obtained by repeated exchanges in the original
procedure, since the envy graph will change after each exchange. Thus, we
obtain a variant of the envy-cycle procedure as shown in Algorithm 2.

Algorithm 2 Envy-Cycle Procedure via Assignment Problem

Input:

A set $N$ of agents, a set $M$ of items, and valuations $v_{i}$ $(i\in N)$.

Output:

An EF1 allocation $(\mathcal{X},\mu)$ of $M$.

  1. 1.

Initialize $\mathcal{X}\leftarrow(X_{i})_{i\in N}$ with $X_{i}=\emptyset$ for
each $i\in N$, and $\mu(i)\leftarrow X_{i}$ for each $i\in N$.

  2. 2.

For each item $k\in M$ (in any order), do the following.

    1. (a)

Find a minimum-weight perfect matching $\mu^{*}$ in $(G(\mathcal{X},\mu),w)$,
and update $\mu\leftarrow\mu^{*}$.

    2. (b)

Pick a bundle $X_{i}\in\mathcal{X}$ that has only one incident edge in
$G(\mathcal{X},\mu)$ (which is $(\mu^{-1}(X_{i}),X_{i})\in\mu$), and update
$X_{i}\leftarrow X_{i}\cup\\{k\\}$.

  3. 3.

Return $\mathcal{X}$.

The number of edges in $G(\mathcal{X},\mu)$ is always bounded by
$\mathrm{O}(n^{2})$. Hence, using the Hungarian method in Step 2(a), this
algorithm correctly finds an EF1 allocation in $\mathrm{O}(mn^{3})$ time,
which is the same bound as Theorem 4.2.

We improve this bound with the aid of Algorithm 1 given in Section 3.1. We
deal with $G(\mathcal{X},\mu)$ as the complete bipartite graph by adding all
the absent edges with sufficiently large weights. Then, in each iteration of
Step 2, only the weights of edges around bundle $X_{i}\in\mathcal{X}$ picked
in Step 2(b) change; hence, we can correctly update the current assignment
$\mu$ and the corresponding potential in $\mathrm{O}(n^{2})$ time by Algorithm
1. Thus, the total computational time is bounded by $\mathrm{O}(mn^{2})$,
leading to the following theorems, which are analogous to Theorems 4.2 and
4.3.

######  Theorem 4.8.

Algorithm 2 correctly finds an EF1 allocation in $\mathrm{O}(mn^{2})$ time.

######  Theorem 4.9.

Suppose that the valuations of all agents are additive, and $v_{i,1}\geq
v_{i,2}\geq\dots\geq v_{i,m}$ for every agent $i\in N$. Then, if the items are
chosen in ascending order of $k=1,2,\dots,m$ in Step 2, Algorithm 2 correctly
finds an EFX allocation in $\mathrm{O}(mn^{2})$ time.

######  Remark 4.10.

The computational time is actually bounded by using the number of edges in
$G(\mathcal{X},\mu)$, which may be $\Omega(n^{2})$ in the worst case. In
practice, however, the algorithm runs much faster than this theoretical bound,
which is also true for the original envy-cycle procedure, because an acyclic
graph tends to be sparse. In addition, in this special situation, without
keeping a potential, we can update the current assignment $\mu$ by solving the
shortest path problem in directed acyclic graphs as follows; this slightly
improves the running time when $G(\mathcal{X},\mu)$ is almost always sparse,
in particular, with $\mathrm{o}(n\log n)$ edges.

Let $(G,w)$ and $(G^{\prime},w^{\prime})$ be the assign-envy graphs for the
allocations $(\mathcal{X},\mu)$ just before and after Step 2(b), respectively.
As $\mu$ is a minimum-weight perfect matching in $(G,w)$, the auxiliary
weighted graph $(G_{\mu},w_{\mu})$ has no cycles by Lemma 4.6. In addition,
since $(G^{\prime},w^{\prime})$ coincides with $(G,w)$ except around
$X_{i}\in\mathcal{X}$ picked in Step 2(b), every cycle in $G^{\prime}_{\mu}$
intersects $X_{i}$. This means that the auxiliary weighted graph
$(G_{\mu^{\prime}},w_{\mu^{\prime}})$ to which we apply Dijkstra’s algorithm
in Step 2 of Algorithm 1 is always acyclic, where
$\mu^{\prime}=\mu\setminus\\{(\mu^{-1}(X_{i}),X_{i})\\}$. Thus, Dijkstra’s
algorithm with a potential is not necessary, and the elementary, linear-time
dynamic programming is sufficient for computing a shortest path from
$\mu^{-1}(X_{i})$ to $X_{i}$.

## Acknowledgments

The proposed algorithms were developed via the preparation of Problem K of the
International Collegiate Programming Contest (ICPC) 2021 Asia Yokohama
Regional Contest [17]. We are grateful to all of the contest organizers and
participants.

This work was supported by JSPS KAKENHI Grant Numbers 20K19743 and 20H00605,
and JST PRESTO Grant Number JPMJPR212B.

## References

  * [1] H. Aziz, B. Li, H. Moulin, and X. Wu. Algorithmic fair allocation of indivisible items: a survey and new questions. arXiv:2202.08713, 2022.
  * [2] S. Barman and S. K. Krishnamurthy. Approximation algorithms for maximin fair division. ACM Transactions on Economics and Computation, 8:1 (2020), No.5, 28pp.
  * [3] R. Barr, F. Glover, and D. Klingman. The alternating path basis algorithm for assignment problems. Mathematical Programming, 13 (1977), pp. 1–13.
  * [4] A. Bernstein, A. Dudeja, and Z. Langley. A framework for dynamic matching in weighted graphs. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing (STOC 2021), pp. 668–681, 2021.
  * [5] I. Caragiannis, D. Kurokawa, H. Moulin, A. D. Procaccia, N. Shah, and J. Wang. The unreasonable fairness of maximum Nash welfare. ACM Transactions on Economics and Computation, 7:3 (2019), No. 12, 32pp.
  * [6] J. Edmonds and R. M. Karp. Theoretical improvements in algorithmic efficiency for network flow problems. Journal of the ACM, 19 (1972), pp. 248–264.
  * [7] K. Hanauer, M. Henzinger, and C. Schulz. Recent advances in fully dynamic graph algorithms — quick reference guide. arXiv:2102.11169, 2021.
  * [8] M. Henzinger. The state of the art in dynamic graph algorithms. In Proceedings of 44th International Conference on Current Trends in Theory and Practice of Computer Science (SOFSEM 2018), pp. 40–44, 2018.
  * [9] M. S. Hung and W. O. Rom. Solving the assignment problem by relaxation. Operations Research, 28:4 (1980), pp. 969–982.
  * [10] M. Iri. A new method of solving transportation-network problems. Journal of the Operations Research Society of Japan, 3 (1960), pp. 27–87.
  * [11] H. W. Kuhn. The Hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2 (1955), pp. 83–97.
  * [12] R. J. Lipton, E. Markakis, E. Mossel, and A. Saberi. On approximately fair allocations of indivisible goods. In Proceedings of the 5th ACM Conference on Electronic Commerce (EC 2004), pp. 125–131, 2004.
  * [13] P. Sankowski. Faster dynamic matchings and vertex connectivity. In Proceedings of the 19th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2007), pp. 118–126, 2007.
  * [14] A. Schrijver. Combinatorial Optimization: Polyhedra and Efficiency, Springer, 2003.
  * [15] D. Stubbs and V. V. Williams. Metatheorems for dynamic weighted matching. In Proceedings of 8th Innovations in Theoretical Computer Science Conference (ITCS 2017), No. 58, 2017.
  * [16] N. Tomizawa. On some techniques useful for solution of transportation network problems. Networks, 1 (1971), pp. 173–194.
  * [17] https://icpc.iisf.or.jp/2021-yokohama/regional/


# Affine-Invariant Robust Training

Oriol Barbany Mayor

###### Abstract

The field of adversarial robustness has attracted significant attention in
machine learning. Contrary to the common approach of training models that are
accurate in average case, it aims at training models that are accurate for
worst case inputs, hence it yields more robust and reliable models. Put
differently, it tries to prevent an adversary from fooling a model. The study
of adversarial robustness is largely focused on $\ell_{p}-$bounded adversarial
perturbations, i.e. modifications of the inputs, bounded in some $\ell_{p}$
norm. Nevertheless, it has been shown that state-of-the-art models are also
vulnerable to other more natural perturbations such as affine transformations,
which were already considered in machine learning within data augmentation.
This project reviews previous work in spatial robustness methods and proposes
evolution strategies as zeroth order optimization algorithms to find the worst
affine transforms for each input. The proposed method effectively yields
robust models and allows introducing non-parametric adversarial perturbations.

##  1 Introduction

Adversarial attacks are the modifications, also known as perturbations in the
literature, that transform an input with the purpose of causing the model to
misclassify the modified inputs. Not any perturbation is allowed, but rather
the modified inputs must be similar to the original input. For this reason,
adversarial attacks specially play a key role in scenarios where reliability
and security are important. The major work on robust training has been
centered on perturbations constrained within an $\ell_{p}$ ball (Goodfellow et
al., 2014; Madry et al., 2017; Shafahi et al., 2019; Cohen et al., 2019). Even
if these perturbations significantly drop the performance of state-of-the-art
models, they are unlikely to exist in the real world. On the other hand,
spatial transforms considered in this work can arise in real scenarios, e.g.
by simply taking the same picture from a different position in the case of
images. It has been shown that neural network-based image classifiers are also
vulnerable to inputs modified with spatial transforms, in particular affine
transforms (Fawzi and Frossard, 2015), which are the core of this project.

The idea of incorporating spatial transforms in the training of neural
networks is not novel, since it was already studied in the context of data
augmentation more than a decade ago. Data augmentation is based on randomly
perturbing the samples from the training dataset and considering both the
perturbed and the original sample in the training procedure. This has been
considered as a measure to extend the cardinality of a dataset and to make the
model invariant to the transformations considered in the process of
augmentation. Moreover, data augmentation can act as a regularizer in
preventing overfitting in neural networks (Simard et al., 2003) and improve
performance in imbalanced class problems (Bowyer et al., 2011). In fact, data
augmentation is such a common practice that deep learning frameworks usually
have built-in data augmentation utilities that not only allow to do affine
transforms but more general projection transforms (Abadi et al., 2015).

Models trained with data augmentation offers relatively small robustness in
front of adversarial attacks than those trained with the non-augmented dataset
(Engstrom et al., 2017). In fact, in the presented experiments (see section
4), they yield worse results than with standard training. Nevertheless, using
carefully crafted perturbations during training instead of random
perturbations before training as in data augmentation can improve the
robustness of a classifier (Engstrom et al., 2017). This project focuses on
implementing algorithms that can find adversarial examples with low
computational cost on arbitrary transformations in order to perform robust
training with them. Yet another goal is to accurately estimate the robustness
against these arbitrary transformations.

##  2 Problem formulation

An adversarial example for an input tensor $\mathbf{x}$ and a classifier $C$,
is defined as a tensor $\mathbf{x}^{\prime}$ that causes $C$ to output a
different label, i.e. $C(\mathbf{x})\neq C(\mathbf{x}^{\prime})$, while the
adversarial example is similar (for a human) to the original input. In this
work we will focus on the case where the inputs are images.

To the best of our knowledge, most of the previous works in adversarial
robustness for images, uses error metrics that are based on pixel-wise signal
differences. Concretely, they assume that $\mathbf{x}^{\prime}$ is a valid
adversarial example if
$\left\lVert\mathbf{x}-\mathbf{x}^{\prime}\right\rVert_{p}\leq\varepsilon$ for
some $p\in[0,\infty]$, typically $p=\\{1,2,\infty\\}$, and a small enough
$\varepsilon$ (Goodfellow et al., 2014; Madry et al., 2017; Shafahi et al.,
2019; Cohen et al., 2019).

Nevertheless, pixel-wise metrics do not reflect the similarity between an
image and its version after applying an affine transform. Hence, small
variations such as rotations and translations can lead to large values in all
$p-$norms while being almost the same for a human. In fact, the concept of
image similarity is very loosely defined in the adversarial robustness
literature, since formally capturing the notion of human perception is
extremely difficult (Sharif et al., 2018). We address this issue in Section 5.

Figure 1 depicts the $\ell_{2}$ and $\ell_{\infty}$ norms when applying
translations to a random image drawn from the MNIST dataset which can be seen
in Figure 2. With the previous notation, that means $\mathbf{x}^{\prime}$ is
the image resulting from translating $\mathbf{x}$, and the graph plots
$\left\lVert\mathbf{x}-\mathbf{x}^{\prime}\right\rVert_{p}$ for
$p\in\\{2,\infty\\}$ against the horizontal and vertical translation in
pixels.

An ideal metric for images, would be invariant against affine transforms, or
at least for small translations (since as for a human, the images of Figure 2
certainly represent the same number), and give a nearly flat loss landscape
with almost zero value. This is not satisfied in Figure 4. Moreover, $p-$norms
are not necessarily symmetric in the space of translations, which would be
another desirable property that can be easily seen to not hold for the
$\ell_{\infty}$ case.

Figure 1: $\ell_{2}$ and $\ell_{\infty}$ norms of difference signal
$\mathbf{x}-\mathbf{x}^{\prime}$ in the space of vertical and horizontal
translations. This means that $\mathbf{x}^{\prime}$ is obtained by translating
the image $\mathbf{x}$ by the quantities (in pixels) of the x and y axis.
Figure 2: Original (left) and transformed (right) image with translation
$(3,-2)$.

##  3 Attack methods

We extend the attack space described in (Engstrom et al., 2017) by allowing
general affine transforms, which include translations and rotations but also
scale changes and shears. For a translation
$(\delta_{u},\delta_{v})\in\mathbb{R}^{2}$, a rotation $\theta\in[0,2\pi)$, a
shear $\phi\in[0,2\pi)$, and a scale $(s_{u},s_{v})\in\mathbb{R}_{+}^{2}$, the
affine transform matrix for parameters
$\mathbf{t}:=\begin{bmatrix}\delta_{u}&\delta_{v}&\theta&\phi&s_{u}&s_{v}\end{bmatrix}^{T}$
is defined as

$\displaystyle
A(\mathbf{t}):=\begin{bmatrix}s_{u}\cos(\theta)&-s_{v}\sin(\theta+\phi)&\delta_{u}\\\
s_{u}\sin(\theta)&s_{v}\cos(\theta+\phi)&\delta_{v}\\\ 0&0&1\\\ \end{bmatrix}$
(1)

where $\delta_{u}=\delta_{v}=\theta=\phi=0$ and $s_{u}=s_{v}=1$ recovers the
identity matrix.

We formulate the robust training procedure as the following saddle point
problem:

$\displaystyle\min_{W}\operatorname{{\mathbb{E}}}_{(\mathbf{x},y)\sim\mathcal{D}}\\!\\!\left[\max_{\mathbf{t}\in\mathcal{S}}\mathcal{L}(W,f(\mathbf{x},A(\mathbf{t})),y)\right]$
(2)

where $W$ are the weights of the neural network, $\mathcal{L}$ is the loss
function, and $f:\mathbb{R}^{M\times N}\times\mathbb{R}^{3\times
3}\to\mathbb{R}^{M\times N}$ performs the affine transform to the input image
$\mathbf{x}\in\mathbb{R}^{M\times N}$. The function $f$ transforms the
coordinates of a point $(u,v)$ on the original image to a new coordinate
$(u^{\prime},v^{\prime})$ in the adversarial image as follows:

$\displaystyle\begin{bmatrix}u^{\prime}\\\ v^{\prime}\\\
1\end{bmatrix}=A(\mathbf{t})\begin{bmatrix}u\\\ v\\\ 1\end{bmatrix}$ (3)

where $A(\mathbf{t})$ is the transformation matrix defined in (1). Note that
since $A(\mathbf{t})$ is a $3\times 3$ matrix, input and output vectors are of
dimension 3. Thus, the addition of a dummy third component is translated into
a proper normalization of the coordinates of the transformed image. Put
differently, $(u^{\prime},v^{\prime},w)$ represents a point at location
$(u^{\prime}/w,v^{\prime}/w)$.

If a pixel gets mapped to a non-integer coordinate, i.e.
$(u^{\prime},v^{\prime})\notin\mathbb{N}^{2}$, we compute its value by
bilinear interpolation. Affine transforms sometimes require the value of
pixels outside the image, hence padding is required. We use zero padding,
which is specially meaningful for MNIST since the background and the padding
are both black in this case (see s 2 and 3).

The set $\mathcal{S}\subseteq\mathbb{R}^{d}$ in (2) is a nonempty compact
space of allowed perturbations containing the parameters associated to the
identity operator. In this project, $\mathcal{S}$ is a hyper-rectangle as the
one defined in (Engstrom et al., 2017), that is, each parameter has
handcrafted limits independent of other parameters as it can be seen in Table
1.

Parameter | Range
---|---
$(\delta_{u},\delta_{v})$ | $[-3,3]^{2}$
$\theta$ (in degrees) | $[-30,30]$
$\phi$ (in degrees) | $[-20,20]$
$(s_{u},s_{v})$ | $[0.6,1.4]^{2}$
Table 1: Parameter ranges for MNIST and CIFAR10.

Moreover, $\mathcal{S}$ can be defined as a ball in some other metric as
discussed in Section 5.

Sticking to the definition of adversarial examples, we want the transformed
image to be similar to the original one. Directly applying the mapping (3) can
give an image that doesn’t satisfy this similarity criterion as depicted in
Figure 3. This is due to the implicit translation that a rotation, scaling and
shearing imply.

Figure 3: Difference between centered and non-centered rotation (directly
using (1)) with $\theta=50^{\circ}$.

In order to avoid the latter, we enforce that the pixel in the center is
mapped to its same position, i.e. we enforce that for
$(u,v)=(\frac{M}{2},\frac{N}{2})$, we have
$(u^{\prime},v^{\prime})=(\frac{M}{2},\frac{N}{2})$ in (3)111Note that for odd
M or N, we have to round $M/2$ or $N/2$ since the coordinates of a pixel are
integral.. Moreover, to simplify the re-centering, we sequentially apply
translation, shearing, scaling and finally a rotation to the input image. When
we say that we apply e.g. a translation, we mean that the only free parameters
are $(\delta_{u},\delta_{v})$ and all the others are set to their value in the
identity mapping.

The following sections introduce the baseline methods described in (Engstrom
et al., 2017) and propose two zeroth order optimization methods that find
adversarial transformations.

###  3.1 Baseline methods

The first baseline method is Grid search, which is based on discretizing
$\mathcal{S}$ and exhaustively examining every possible transformation in the
discretized space. These transformations are evaluated sequentially, and for
each image we pick the worst transformation according to Assumption 1. The
Grid search method requires a lot of forward passes even when the number of
free parameters is low. Therefore, this method is only used to test the
robustness of trained models but not in the training procedure.

######  Assumption 1.

Let $\mathbf{x}_{j}^{\prime}:=f(\mathbf{x},A(\mathbf{t}_{j}))$ be the image
transformed with transform $\mathbf{t}_{j}$. When evaluating a set of
transforms $\\{\mathbf{t}_{1},\dots,\mathbf{t}_{i},\dots,\mathbf{t}_{N}\\}$
for a given image $\mathbf{x}$, transformation $\mathbf{t}_{i}$ is considered
worse than the transformations $\\{\mathbf{t}_{1},\dots,\mathbf{t}_{i-1}\\}$,
if it either misclassifies the input while the previous transforms don’t, i.e.
$C(\mathbf{x}_{i}^{\prime})\neq y$, $C(\mathbf{x}_{j}^{\prime})=y$ $\forall
j\in\\{1,\dots,i-1\\}$, or achieves the maximum loss in case of a tie. That
is, when either $\textbf{t}_{i}$ doesn’t misclassify neither do the previous
ones, i.e. $C(\mathbf{x}_{j}^{\prime})=y$ $\forall j\in\\{1,\dots,i\\}$ but
$i=\operatorname*{arg\,max}_{j\in\\{1,\dots,i\\}}\mathcal{L}(W,\mathbf{x}_{j}^{\prime},y)$,
or $\textbf{t}_{i}$ misclassifies as some of the previous ones do and
$i=\operatorname*{arg\,max}_{j\in\\{1,\dots,i\\}:C(\mathbf{x}_{j}^{\prime})\neq
y}\mathcal{L}(W,\mathbf{x}_{j}^{\prime},y)$.

The second baseline method, which is used both in training and evaluation, is
Worst-of-$k$. This method is based on sampling $k$ different transformations
uniformly and independently sampled from $\mathcal{S}$ and using the worst one
(following Assumption 1) as an adversarial example. As we increase $k$, this
attack interpolates between a random choice and Grid search.

Note that both baseline methods require specific numerical bounds for each
parameter giving a feasible transform, i.e. we need to explicitly specify the
support both to discretize and to define the uniform distribution over
$\mathcal{S}$. These methods have been adapted from (Engstrom et al., 2017)
and tested in the setup described in Section 4.

###  3.2 Evolution Strategies

As discussed in (Engstrom et al., 2017), due to the highly non-concave loss
landscape in the space of the parameters defining the spatial transform,
first-order methods lack any guarantee of optimality and cannot find good
adversarial examples in practice. The loss landscape of some random samples
from MNIST can be seen in Figure 4.

Algorithm 1 provides the pseudo-code for a general evolution strategy adapted
from (Akimoto and Hansen, 2018) for constrained optimization. The two
algorithms proposed in this section are variants that differ on how to update
the previous parameters.

Algorithm 1 EvolutionStrategy$(F,\mathcal{S})$

1: Input: $F:\mathbb{R}^{d}\to\mathbb{R}$, the function to minimize and
$\mathcal{S}\subseteq\mathbb{R}^{d}$, the constraint set implicitly given by a
membership oracle

2: Initialize: Mean vector $\mathbf{m}\in\mathcal{S}$, step-size
$\sigma\in\mathbb{R}^{+}$ and covariance matrix $C\in\mathbb{R}^{d\times d}$

3: while Stopping criterion not met do

4: repeat

5: $\mathbf{t}\sim\mathbf{m}+\sigma\mathcal{N}(\mathbf{0},C)$

6: until $\mathbf{t}\notin\mathcal{S}$, given by membership oracle

7: Update$(\mathbf{t},\mathbf{m},\sigma,C,F)$

8: end while

9: return $\mathbf{t}$

In Algorithm 1, we need to initialize $\mathbf{m}\in\mathcal{S}$, yet we claim
that only a membership oracle is needed. Note that $\mathcal{S}$ represents
the set of allowed transformations, so one could initialize $\mathbf{m}$ to be
the parameters that give the identity operator. In case of having the
aforementioned affine transform, this translates into setting
$\mathbf{m}=\begin{bmatrix}\delta_{u}&\delta_{v}&\theta&\phi&s_{u}&s_{v}\end{bmatrix}^{T}=\begin{bmatrix}0&0&0&0&1&1\end{bmatrix}^{T}$,
which is contained in $\mathcal{S}$. Following the notation in (2), we can
take $F:\mathbf{t}\mapsto-\mathcal{L}(W,f(\mathbf{x},A(\mathbf{t})),y)$ for a
given image $\mathbf{x}$ with associated ground-truth label $y$, model weights
$W$ and transformation $A(\mathbf{t})$. Nevertheless, in practice we choose
the worst transforms following Assumption 1. This is because achieving a
misclassification and a maximal loss are not the same even though both are
related and the goal of an adversarial attack is to misclassify the
transformed inputs. For the sake of notation, we represent this criteria with
$F$. Note that reversing the sign is needed since Algorithm 1 requires a
function to minimize and we are trying to find the transformation that
maximizes the loss. Nevertheless, in the experiments we take $F$ such that the
number of misclassified samples have more importance than the loss.

Unlike the sampling from a uniform distribution or the Grid search used in
(Engstrom et al., 2017) to find adversarial examples, given that the Gaussian
has infinite support, we could get a transform that is not in the set of
feasible transformations $\mathcal{S}$. If this is the case, we resample the
transformations such that $\mathbf{t}\notin\mathcal{S}$, for which we only
need the membership oracle. The number of times that we will have to resample
depends on the definition of $\mathcal{S}$. In the following section, we will
further discuss the implications on the current setup.

In the current experiments, the stopping criterion of Algorithm 1 is simply a
maximum number of iterations in the outer loop, which is equivalent to the
number of calls to Update. As further explained in the next sections, in the
case of (1+1)-ES, only one evaluation of $F(\mathbf{t})$ and hence one forward
pass is needed to update the parameters. However, in the case of CMA-ES one
needs several evaluations inside the Update function.

The key point of using evolution strategies is that, instead of naively
sampling points in $\mathcal{S}$, we can learn the worst transforms for each
image separately at a given iteration. Thus, this can lead to a more robust
model by incorporating the worst transforms during each training step, which
by assumption change smoothly over the transformation space from one iteration
to the next one. For this reason, the parameters $\mathbf{m},\sigma,C$ are
stored for each image in the training set and when we call Algorithm 1, we
initialize $\mathbf{m},\sigma,C$ to the previous parameter values
corresponding to the same image obtained on the previous epoch. The mean
corresponds to the worst transform (according to Assumption 1) seen so far for
each image. Note that in evaluation, Algorithm 1 is only run once per image
and hence the parameters $\mathbf{m},\sigma,C$ are not stored for each image
once the algorithm terminates.

Figure 4: Loss landscape of three random example drawn from the MNIST dataset
when performing horizontal translations and rotations restricted according to
Table 1. The model and its loss function are described in Section 4. In
particular this model has been fully-trained using the standard training
method.

####  3.2.1 (1+1)-ES

The simplest evolution strategy used to solve (2) is (1+1)-ES, which is based
on sampling an isometric Gaussian with varying mean and variance. This means
that in Algorithm 1, we enforce $C=I$, the identity matrix. To effectively
explore all the allowed transformation space $\mathcal{S}$, we normalize all
the parameters in such a way that its range is $[-1,1]^{d}$. Note that this
needs $\mathcal{S}$ to be symmetric around zero in every dimension, so we have
to not only re-scale all the parameters but also shift the scale for each axis
$(s_{u},s_{v})$ to be centered around 1 (see Table 1). So in this case, we
initialize $\mathbf{m}=\mathbf{0}$.

Note that for the current setup, the probability of having the new sample in
$\mathcal{S}$ is at least $\frac{1}{2d}$, where $d$ is the number of free
parameters of the transformation (1), given that the previous iterate was in
$\mathcal{S}=[-1,1]^{d}$. In practice, we don’t have to resample $\mathbf{t}$
too often given that the worst case is attained when sampling from a Gaussian
centered in one corner of $\mathcal{S}$. Moreover, the sampling of a Gaussian
doesn’t incur much overhead cost and the computation of the attack is
dominated by the number of forward passes.

To avoid resampling, especially in the first iterates, we set the initial
variance $\sigma_{0}$ such that the probability to resample in the first
iterate, namely $\epsilon$ is low. We can analytically compute the initial
sigma for a given $\epsilon$ as follows:

$\displaystyle\epsilon$
$\displaystyle:=\Pr[\mathbf{X}\notin[-1,1]^{d}]=\Pr[\lor_{i=1}^{d}X_{i}\notin[-1,1]]=1-\Pr[\land_{i=1}^{d}X_{i}\in[-1,1]]$
(4)
$\displaystyle=1-\Pr[X_{1}\in[-1,1]]^{d}=1-\left[\int_{-1}^{1}\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}e^{-\frac{x^{2}}{2\sigma_{0}^{2}}}dx\right]^{d}$
(5)

where $X_{i}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma_{0}^{2})$. The
choice of $\epsilon\approx 0.5$ gives the $\sigma_{0}=0.4$ that we use along
the experiments with $d=3$. However, when using this method in evaluation, we
found better results setting $\sigma_{0}=0.75$ since Algorithm 1 only runs for
one times in this case. There are several variants of the (1+1)-ES algorithm.
In this case, we follow the one-fifth success rule (Auger, 2009), resulting in
the update described in Algorithm 2. The variance increases if there is
improvement with respect to the previous iteration and decreases for
unsuccessful iterations. The factors of increase and decrease are theoretical
to reach the one-fifth success rule described in (Auger, 2009).

Algorithm 2 Update$(\mathbf{t},\mathbf{m},\sigma,C,F)$ for (1+1)-ES with
1/5-th success rule

1: Input: Transform $\mathbf{t}$, mean $\mathbf{m}$, step-size $\sigma$,
Covariance $C$, and $F:\mathbb{R}^{d}\to\mathbb{R}$, the function to minimize

2: if $F(\mathbf{t})\leq F(\mathbf{m})$ then

3: $\mathbf{m}\leftarrow\mathbf{t}$

4: $\sigma\leftarrow 1.5\sigma$

5: else

6: $\sigma\leftarrow 1.5^{-1/4}\sigma$

7: end if

Given the nature of Algorithm 2 along with the memory of the evolution
strategy in training, $\sigma$ may vanish. Empirically, we observed that
indeed the step-size tends to vanish and even if theoretically possible, it
never explodes. When $\sigma\approx 0$ we get a deterministic behavior, i.e.
$\mathbf{t}$ is no longer updated and takes the value
$\mathbf{t}\approx\mathbf{m}$. This happens for example when $\mathbf{m}$
corresponds to a local maximum. In this case, all the transformations that are
close to $\mathbf{m}$, which are the ones sampled with higher probability are
worst in terms of the objective value. Thus in such situation it is virtually
impossible to improve the current solution if $\sigma$ is small. To tackle
this, we introduce a restart condition that sets the variance back to its
initial value $\sigma_{0}$. This restart condition is based on comparing the
current variance to the equivalent of having 5 decreases in a row from the
initial iteration. If the current $\sigma$ is lower than this threshold, we
restart it. This translates into adding

$\displaystyle\sigma\leftarrow\begin{cases}\sigma_{0}&\text{if
}\sigma\leq\sigma_{0}(1.5^{-1/4})^{5}\\\ \sigma&\text{otherwise}\end{cases}$
(6)

at the end of Algorithm 2 for (1+1)-ES.

####  3.2.2 CMA-ES with diagonal decoding and unitary determinant

Especially for larger values of $d$, the number of free parameters in (1), it
is interesting to learn relations among parameters. In other words, it could
be that there is a correlation between parameters that yield a large loss and
the independence assumption of (1+1)-ES would be to simplistic in this case.

In order to see if indeed the (1+1)-ES algorithm was too simple, we tested
CMA-ES, which allows full covariance matrices (Akimoto and Hansen, 2018) and
not just a fixed identity matrix. In particular we use CMA-ES with diagonal
decoding. which defines the Update function used in this case (Akimoto, 2019).

Avoiding a deterministic behavior is not as easy as comparing a scalar with a
threshold as in (6) for the case of (1+1)-ES. The determinant, which can be
interpreted as the oriented volume of the hyper-parallelepiped defined by the
column vectors of $C$, can be used to spot such deterministic behavior. This
latter is associated to the case when the determinant, and hence the volume in
the former interpretation, is too small.

Instead of restarting the covariance matrix, we fix our determinant to be
unitary, which translates into adding the step

$\displaystyle C\leftarrow\frac{C}{\det(C)}$ (7)

at the end of the Update procedure described in (Akimoto, 2019) to prevent
shrinking.

Even if in this case there is no need to normalize the parameters since the
Gaussian in CMA-ES is not necessarily isometric, it is simpler to do so in the
case of $\mathcal{S}$ being an hyper-rectangle. This simplifies the
implementation of the membership oracle and allows using the same
initialization as for the (1+1)-ES algorithm. The caveat of this algorithm, is
that in order to update the full covariance matrix, we need to perform several
evaluations of the objective function with the same parameters. Such number of
evaluations with fixed parameters is known as the population size. In
particular, this population size has to be of at least 3, and by default it is
set to $4+\lfloor 3\log(d)\rfloor$ (Akimoto, 2019). If comparing algorithms
with respect to the number of forward passes, one should take this into
account and reduce the number of iterations in the outer loop of Algorithm 1,
which is the stopping criterion of the latter.

As stated before, we use this method for larger values of $d$, implying that
the probability of resampling is higher. Differing from (1+1)-ES, we don’t
have independence of coordinates. Moreover, the second equality of (4) doesn’t
hold in this case due to the dependence of coordinates. Thus, it could happen
that given the value in some coordinates in $\mathcal{S}$, it is very unlikely
to get feasible values, i.e.
$\Pr[\land_{i\in\mathcal{X}}X_{i}\in[-1,1]|\land_{i\notin\mathcal{X}}X_{i}\in[-1,1]]$
may be very low for some $\mathcal{X}\subseteq[d]$.

To avoid resampling, we introduced two variants that take advantage of the
normalization of the transform space into the hyper-cube
$\mathcal{S}=[-1,1]^{d}$. These are projection onto $\mathcal{S}$ and applying
$\tanh(\cdot)$ element-wise to the sampled transformation (Meunier et al.,
2019). Applying the hyperbolic tangent to the samples of the Gaussian is
compatible with the infinite support of such distribution yet avoids
concentration in the extreme points as with the projection approach.

##  4 Experiments

The set of allowed transformations ($\mathcal{S}$ as defined in (2)) in all
the presented experiments is defined as the Cartesian product of all the sets
of Table 1. Note that in the latter, $\mathcal{S}$ is defined for specific
datasets, which is because parameters such as the translation and the scaling
are more meaningful if they depend on the image size or the nature of the
dataset. We only allow a translation of $10\%$ of the image size222MNIST
images are of size $M=N=28$ pixels and CIFAR10 images of size $M=N=32$ pixels.
After rounding to the nearest integer we get the $\pm 3$ of Table 1. as in
(Engstrom et al., 2017). The other limits were chosen so as to give visually
similar transformed images in the limit cases. Put differently, even in limit
cases, it should be clear for a human that the label of the original and
transformed images are the same.

Model architecture: For MNIST, we use a convolutional neural network
consisting of two convolutional layers with 32 and 64 filters each. Both
layers are followed by a max pooling layer. At the end of the last
convolutional layer, two fully connected layers follow. ReLU is chosen as
activation function for all cases, and Cross Entropy is used as the loss
function $\mathcal{L}$. This model corresponds to the CNN of the MNIST
Challenge (Madry et al., 2017), which is different from the model for the same
dataset used in (Engstrom et al., 2017). This means that the results of the
experiments won’t necessarily be the same even though the Worst-of-$k$ and
Grid search methods are virtually the same as in (Engstrom et al., 2017). We
used a different model since the description of the MNIST model used in
(Engstrom et al., 2017) only explains the differences with respect to a
baseline model given in the TensorFlow tutorial (Abadi et al., 2015), which is
no longer existent. For all methods, we used 10k training iterations, that is
roughly 21 epochs.

Data augmentation: The data augmentation procedure used in the following
experiments is the same used in (Engstrom et al., 2017). This consists in
doing random flips in the horizontal dimension of an image. In particular,
with probability 0.5 we get the same image and with probability 0.5 we get the
flipped version. Note that a flip in the horizontal dimension can be modelled
as an affine transform with $\theta=\phi=\pi$. Clearly this transform is not
included in $\mathcal{S}$ (see Table 1) yet it yields more robust models in
(Engstrom et al., 2017). As a further clarification, adversarial perturbations
are computed at each training iteration and the transformed versions of the
inputs are used as the training dataset for such iteration. Nevertheless, for
this implementation of data augmentation, perturbations are random and
computed once before training starts. Such transforms and the original inputs
are used as the training dataset, which remains static in all the training
procedure.

Table 2 shows the accuracy of different classifiers on MNIST with affine
transforms restricted to only rotation and translation (both vertical and
horizontal), hence $d=3$. Note that the Grid search method is clearly the most
effective one333When evaluating adversarial attacks, the lower the accuracy in
evaluation, the better. When the attack is incorporated in training, we expect
the accuracy drop to be significantly damped with respect to (at least)
standard training., but with the current setup and the granularity described
in (Engstrom et al., 2017), it takes 775 forward passes to find each attack.
Thus, it is infeasible to incorporate this method in training, even with only
3 dimensions. When $d$ increases, the number of forward passes for Grid search
explodes, since it is $\mathcal{O}(\exp(d))$.

Table 2: Accuracy averaged over 3 classifiers on MNIST trained with different seeds. The rows and the columns represent the training and evaluation methods respectively. For each pair of training method and evaluation, the mean and standard deviation of the accuracies is provided. The best accuracies for each attack and the overall best methods for robust training and attacking are in bold font. Model | Nat. | W-10 | Grid | (1+1)-ES-10 | (1+1)-ES-100
---|---|---|---|---|---
Standard | $\boldsymbol{99.05\pm 0.06}$ | $27.19\pm 0.15$ | $0.00\pm 0.00$ | $37.61\pm 0.40$ | $37.87\pm 0.52$
Aug. 30 | $98.40\pm 0.13$ | $14.18\pm 0.67$ | $0.00\pm 0.00$ | $21.93\pm 0.80$ | $21.55\pm 0.90$
W-10 | $98.91\pm 0.00$ | $\boldsymbol{96.81\pm 0.01}$ | $87.16\pm 0.00$ | $96.82\pm 0.03$ | $96.82\pm 0.03$
(1+1)-ES-10 | $98.88\pm 0.00$ | $95.97\pm 0.05$ | $86.24\pm 0.00$ | $96.09\pm 0.05$ | $96.08\pm 0.07$
(1+1)-ES-20 | $98.96\pm 0.02$ | $96.77\pm 0.01$ | $87.05\pm 0.00$ | $96.17\pm 0.05$ | $96.17\pm 0.05$
(1+1)-ES-50 | $98.92\pm 0.08$ | $96.76\pm 0.07$ | $\boldsymbol{90.08\pm 0.00}$ | $\boldsymbol{96.97\pm 0.06}$ | $\boldsymbol{96.97\pm 0.00}$

As seen in Table 2, the evaluation with Grid search always yields lower
accuracy than the other methods, thus being the best procedure to find
adversarial examples for trained models. Nevertheless, even evaluating with
Grid search can take a long time and thus other methods of evaluation can be
preferred in some setups. Also note than Worst-of-10 can find adversarial
examples yielding lower accuracies than (1+1)-ES for most of the models. This
is because in this case, the key point of learning the parameters for each
individual image along several epochs cannot be used. Furthermore, in this
case having to sample from Gaussians may impair that adversarial examples
yielding lower accuracies in classification, usually located near the extreme
points of $\mathcal{S}$, are explored.

Given that the (1+1)-ES algorithm only updates the transform if it is better,
again regarding the same criterion described in Assumption 1, adding (1+1)-ES
iterations can only improve. However, when evaluating with (1+1)-ES, only mild
improvements are observed when augmenting the number of forward passes from 10
to 100.

Overall, one can observe that data augmentation is the worst approach for
robustness in all cases, yielding even worse results than with standard
training. This is attributed to the fact that the transformation used in data
augmentation is not representative of the adversarial examples seen when
evaluating. We tried this method of data augmentation since it yielded better
results than standard training in (Engstrom et al., 2017) even if the
transform is also not in the attack space. Intuitively, incorporating the same
affine transforms in data augmentation should increase the robustness of the
model, but we wanted to take the setting of (Engstrom et al., 2017) as a
baseline.

W-10 and (1+1)-ES obtain very similar accuracies, and in particular the
variant with 50 iterations, (1+1)-ES-50, achieves better results according to
the strongest attack (Grid search) for $d=3$. Nonetheless, having 5 times as
many forward passes per training iteration doesn’t increase much the accuracy
in general as seen in Table 2, but it roughly multiplies the training time by
a factor of 5. Thus, in some setups, alternatives with less forward passes may
be preferred. Finally, the best performance with natural evaluation was
achieved by the standard model in both experiments as expected. Thus, in
comparison to standard models, robust training yields models that have
slightly worst accuracy for the test dataset but much better accuracy for the
same inputs after being adversarially transformed.

Table 3 shows the accuracies for different models on MNIST in the case of
general affine transforms. As aforementioned, the focus of this project is
building robust models rather than proposing attack methods. Moreover the best
attack, as seen in Table 2 can be achieved by using Grid Search and, as seen
in Table 2, ES-(1+1) is not better than Worst-of-$k$ in evaluation. Thus, in
Table 3 only natural and Worst-of-$k$ evaluations are provided.

Regarding the CMA-ES method, the hypothesis of a possible non termination of
the inner loop of Algorithm 1 discussed in Section 3.2.2 is confirmed. Indeed,
in all the cases it happened that eventually, given the values of some
parameters of the transform, it is nearly impossible to get feasible values
for the rest of parameters. Therefore, the experiments for CMA-ES were only
conducted for the projected and the hyperbolic tangent approaches. For this
methods we use 12 forward passes, that is 2 calls to Update for each training
iteration. This is because in this case the population size is $4+\lfloor
3\log 6\rfloor=6$. In Table 3, we also test whether the determinant
normalization (7) helps or not. The suffix FreeDet in the models of Table 3
means that the determinant is free in such case, hence no determinant
normalization.

Table 3: Accuracy averaged over 3 classifiers on MNIST trained with different seeds along with its standard deviation. In this case, we allow general affine transformations, that is $d=6$. Model | Nat. | W-10 | W-100
---|---|---|---
Standard | $\mathbf{98.91\pm 0.00}$ | $0.16\pm 0.00$ | $0.00\pm 0.00$
Augmented | $98.61\pm 0.04$ | $0.06\pm 0.01$ | $0.00\pm 0.00$
W-10 | $97.05\pm 0.00$ | $\mathbf{75.91\pm 0.00}$ | $\mathbf{39.58\pm 0.00}$
(1+1)-ES-10 | $92.29\pm 9.95$ | $54.03\pm 30.10$ | $24.14\pm 15.62$
CMA-ES-12-Proj | $9.80\pm 0.00$ | $9.80\pm 0.00$ | $9.80\pm 0.00$
CMA-ES-12-Tanh | $9.80\pm 0.00$ | $9.80\pm 0.00$ | $9.80\pm 0.00$
CMA-ES-12-Proj-FreeDet | $34.48\pm 0.00$ | $0.30\pm 0.00$ | $0.00\pm 0.00$
CMA-ES-12-Tanh-FreeDet | $31.87\pm 0.00$ | $0.59\pm 0.00$ | $0.00\pm 0.00$

As seen in Table 3, for $d=6$ the seed has less influence in all cases but
(1+1)-ES-10. In this latter, we can see huge variance on the results. In the
other cases, there is barely randomness in the accuracies of the models
trained with the same procedure and different seeds. Even though the
approaches deriving from CMA-ES-12 are promising in theory, the results are
worst than those achieved with $C=I$ as in (1+1)-ES-10. This is probably due
to the non-convergence of the covariance matrix. Empirically, we observed that
the spectral norm of the difference of consecutive covariance matrices is not
decreasing as it would be expected. Thus, the approaches of CMA-ES described
in Section 3.2.2 are not a good alternative for efficiently training robust
models. This alternative would probably improve if updating the covariance
matrix for more than just 2 times with fixed parameters. Nevertheless, at the
cost of having more forward passes, one could use Worst-of-$k$ which would
approach the grid search performance as $k$ increases (see Section 3.1).

The best accuracies for affine adversarial attacks are obtained with the model
trained with W-10. In 2 out of the 3 models trained with (1+1)-ES the results
were very similar to those obtained with W-10. In particular, the accuracies
were in both cases 98.04, 71.41 and 33.16 for natural, W-10 and W-100
evaluation respectively. Nevertheless, the results with the third seed were
much worse, which shows the instability of this method with $d=6$. Thus in
this case W-10 would clearly be the preferred method for robust training.

##  5 Conclusions and future work

To sum up, the proposed methods allow to effectively train robust models
against affine transform attacks. The initialization, handling of
deterministic behavior and update strategy gives room for the improvement in
the use of Algorithm 1 for finding better adversarial examples.

Other more complex adversarial attacks such as JPEG compression with different
quality levels such as in (Hendrycks and Dietterich, 2019) could be
constrained within a ball defined using the Structural Similarity Index (SSIM)
described in (Zhou Wang et al., 2004) as the metric to determine which inputs
are valid adversarial examples. This metric was actually conceived in the
image processing literature taking into account the human visual system and
aiming to measure the image quality based on an initial uncompressed or
distortion-free image as reference. Hence, it is more meaningful to define
which adversarial examples do not change the semantics of the image instead of
the definitions using $\ell_{p}$ norms or handcrafted limits as in (Engstrom
et al., 2017).

In general, the possibility of implicitly give $\mathcal{S}$ by a membership
oracle for (1+1)-ES allows using more complex metrics and looking for affine-
invariant alternatives that give meaningful transforms for a ball of fixed
radius. Using an implicit definition of $\mathcal{S}$ may also be possible for
CMA-ES but the presented formulation failed to terminate in practice. In order
to tackle with the infinite resampling, we introduced two variants that don’t
need to resample and thus can compute adversarial examples at a speed
comparable to the Worst-of-$k$ method with $k$ being the same number of
forward passes. Nevertheless this approaches are not compatible with
$\mathcal{S}$ given by a membership oracle, but need a projection operator in
the projected version, or a constrained set described by an hyper-cube in the
variant with hyperbolic tangent.

##  6 Acknowledgments

I would like to express my sincere gratitude to my advisor Dr. Sebastian Stich
for introducing me to the marvelous field of adversarial robustness and for
its dedication and continuous support.

I am also particularly grateful for the assistance given by Maksym
Andriushchenko. He brought us back to earth by providing a paper implementing
what we thought was a great novel idea (Engstrom et al., 2017), and he also
provided his broad knowledge in adversarial robustness.

Last by not the least, I would like to thank Prof. Martin Jaggi for his
valuable feedback.

## References

  * Abadi et al. (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015.  Software available from tensorflow.org.
  * Akimoto (2019) Youhei Akimoto.  Cma-es with diagonal decoding.  https://gist.github.com/youheiakimoto/1180b67b5a0b1265c204cba991fa8518, 2019\.
  * Akimoto and Hansen (2018) Youhei Akimoto and Nikolaus Hansen.  Cma-es and advanced adaptation mechanisms.  In _GECCO_ , pages 720–744, 07 2018.  doi: 10.1145/3205651.3207854.
  * Auger (2009) Anne Auger.  Benchmarking the (1+1) evolution strategy with one-fifth success rule on the bbob-2009 function testbed.  _GECCO_ , 2009.  doi: 10.1145/1570256.1570342.
  * Bowyer et al. (2011) Kevin W. Bowyer, Nitesh V. Chawla, Lawrence O. Hall, and W. Philip Kegelmeyer.  SMOTE: synthetic minority over-sampling technique.  _CoRR_ , abs/1106.1813, 2011.
  * Cohen et al. (2019) Jeremy M Cohen, Elan Rosenfeld, and J. Zico Kolter.  Certified adversarial robustness via randomized smoothing, 2019.
  * Engstrom et al. (2017) Logan Engstrom, Brandon Tran, Dimitris Tsipras, Ludwig Schmidt, and Aleksander Madry.  Exploring the Landscape of Spatial Robustness.  _arXiv e-prints_ , art. arXiv:1712.02779, Dec 2017.
  * Fawzi and Frossard (2015) Alhussein Fawzi and Pascal Frossard.  Manitest: Are classifiers really invariant? _CoRR_ , abs/1507.06535, 2015.
  * Goodfellow et al. (2014) Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and Harnessing Adversarial Examples.  _arXiv e-prints_ , art. arXiv:1412.6572, Dec 2014.
  * Hendrycks and Dietterich (2019) Dan Hendrycks and Thomas Dietterich.  Benchmarking neural network robustness to common corruptions and perturbations.  _Proceedings of the International Conference on Learning Representations_ , 2019.
  * Madry et al. (2017) Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks, 2017.
  * Meunier et al. (2019) Laurent Meunier, Jamal Atif, and Olivier Teytaud.  Yet another but more efficient black-box adversarial attack: tiling and evolution strategies, 2019.
  * Shafahi et al. (2019) Ali Shafahi, Mahyar Najibi, Amin Ghiasi, Zheng Xu, John P. Dickerson, Christoph Studer, Larry S. Davis, Gavin Taylor, and Tom Goldstein.  Adversarial training for free! _CoRR_ , abs/1904.12843, 2019.
  * Sharif et al. (2018) Mahmood Sharif, Lujo Bauer, and Michael K. Reiter.  On the suitability of l${}_{\mbox{p}}$-norms for creating and preventing adversarial examples.  _CoRR_ , abs/1802.09653, 2018.
  * Simard et al. (2003) P. Y. Simard, D. Steinkraus, and J. C. Platt.  Best practices for convolutional neural networks applied to visual document analysis.  In _Seventh International Conference on Document Analysis and Recognition, 2003. Proceedings._ , pages 958–963, Aug 2003.  doi: 10.1109/ICDAR.2003.1227801.
  * Zhou Wang et al. (2004) Zhou Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ , 13(4):600–612, April 2004.  doi: 10.1109/TIP.2003.819861.