\section{Introduction}
A randomly generated deterministic
finite
automaton (DFA)
on $n$ states
and $k$ symbols
is drawn as follows:
for each state and each
of the $k$
symbols
in the alphabet,
the transition arrow's destination is chosen uniformly at random among the $n$ states; the $n
k
$
random choices are independent\footnote{By symmetry, we may always take
the state $q=1$ to be the starting state.}.
Then each state is chosen to be accepting (or not)
independently with probability $1/2$.
This natural model for a ``typical'' DFA
goes back to \cite{Trakhtenbrot73} and
was considered in \cite{DBLP:conf/alt/AngluinEKR10,kont-nadler08}
in the context of learning theory. In particular,
in \cite{DBLP:conf/alt/AngluinEKR10} it is shown (perhaps surprisingly)
that random DFAs possess sufficient complexity to
embed nontrivial parity problems.

Let us define the {\em state complexity}
of a DFA $M$ as the number of states in the canonical (minimal) DFA
equivalent to $M$, and denote it by $\nrm{M}$.
In this paper, we study the
state complexity
of random DFAs in the model defined above.
\paragraph{Related work}
We are not aware of previous literature dealing with the specific problem we have posed.
The somewhat related problem of enumerating finite automata according to various criteria has
been extensively studied; see
\cite{MR1990452} and the references therein.
Some recent results include enumeration of minimal automata
\cite{DBLP:conf/stacs/BassinoDS12}, generation of random complete DFAs
\cite{MR2114149}, and enumeration and generation of accessible DFAs \cite{MR2347395}.
In a different line of enquiry, Pittel investigated the distributions induced by
transitive closures \cite{MR690140} and rumor spreading \cite{MR873245}.

\section{Background and notation}
We use standard automata-theoretic notation throughout; the reader is
referred to \cite{DBLP:books/daglib/0011126,ShengYu}
for background.
We put
$[n]=\set{0,\ldots,n-1}$.
Thus, $[k]$ is a $k$-ary alphabet and
$[k]^*$
is the
set of all finite
words (strings) over
this alphabet.
The notation
$\abs{
\cdot
}$
is used for both word length and set cardinality.
Standard order-of-magnitude notation $o(\cdot)$ and $O(\cdot)$ is used,
as well as their ``with high probability'' variants $o_P(\cdot)$
and $O_P(\cdot)$.
The $\tilde O(\cdot)$ notation ignores polylog factors.

An $n$-state $k$-ary Deterministic Finite-state Automaton
is
a
tuple
$M=(Q,q_0,A,\delta)$ where
\begin{itemize}
\item $Q=[n]$
is the set of states
\item $q_0=1$ is the starting state;
\item $A\subseteq [n]$ is the set of accepting states;
\item $\delta:[n]\times[k]\to [n]$ is the transition function.
\end{itemize}
The
transition function $\delta$ may be extended to
$[n]\times[k]^*$ via the recursion
\begin{eqnarray}
\label{eq:drecurs}
\delta(q,u_1u_2\cdots u_{n})
= \delta( \delta(q,u_1),u_2\cdots u_{n})
.
\end{eqnarray}

If the accepting states are unspecified,
the transition function $\delta$ induces
a directed multigraph on $n$ nodes with regular outdegree $k$,
called
a $k$-ary {\em semiautomaton}.

We recall the standard equivalence relation over the states of a DFA:
a word $x\in[k]^*$
{\em distinguishes} between
the states $p,q\in[n]$
if exactly one of
the states $\delta(p,x)$, $\delta(q,x)$ is accepting.
If no $x\in[k]^*$ distinguishes between $p$ and $q$,
these states are {\em equivalent}, denoted by $p\equiv q$.

A standard high-level algorithm\footnote{
Hopcraft's celebrated algorithm \cite{MR0403320}
for minimizing a DFA has runtime complexity $O(n\log n)$.}
 for minimizing a DFA
proceeds in two stages:
\begin{itemize}
\item {\tt REMOVE-UNREACHABLE}: Remove all states $q$ such that
there is no directed path from the starting state $q_0$ to $q$.
\item {\tt COLLAPSE-EQUIVALENT}: Collapse each set of mutually
equivalent states into a single state.
\end{itemize}

\section{Main results}
Our main result is the following estimate on the
state complexity of random DFAs:
\bethn
\label{thm:main}
Let $M_n\k$ be a random
DFA
on $n$ states and $k$ symbols
drawn uniformly from
$[n]^{[n]\times[k]}\times2^{[n]}$.
Then, for any fixed $k\ge2$ and sufficiently large $n$,
\begin{eqnarray}
\label{eq:Mnkdev}
\pr{\abs{\nrm{M_n\k}-
\alpha_k n}>\sqrt n\log n
}
&=&
\Theta\paren{\frac{
1
}{n^k}},
\end{eqnarray}
where $\alpha_k$
is unique positive root\footnote{
A closed-form expression for $\alpha_k$
is possible
via the Lambert $W$ function \cite{MR1809988}:
$\alpha_k = 1+W(-ke^{-k})/k$.
This constant also appears (as $\omega_k$) in \cite{DBLP:conf/stacs/BassinoDS12}
and seems to be
intimately related to
average reachability properties of semiautomata in
$[n]^{[n]\times[k]}$ under the uniform measure.
} of
$x=1-e^{-kx}$.
In particular,
\begin{eqnarray*}
\mathbb{E}\nrm{M_n\k} = \alpha_k n + O(\sqrt n\log n).
\end{eqnarray*}
\enthn
\begin{rem}
\label{rem:pittel}
Observe that
$
0.7968\approx
\alpha_2 < \alpha_3 <\ldots<\alpha_\infty=1$.
For $k=1$, the behavior
of $\nrm{M_n\k}$
is qualitatively different
than described in Theorem \ref{thm:main}.
The
equation $x=1-e^{-x}$
has no positive solution
and
$\mathbb{E}\nrm{M_n^{(1)}} = \Theta(\sqrt n)$,
which follows from the analysis in
\cite{MR690140}.
\end{rem}
\begin{rem}
\label{rem:lb}
The lower bound of $\Omega(1/n^k)$ in (\ref{eq:Mnkdev}) is trivial,
since
with probability
$1/n^k$,
all of state 1's arrows point back to itself and
$\nrm{M_n\k}=1$.
\end{rem}

Our proof of Theorem \ref{thm:main} proceeds in two principal stages.
First we show that in our model, a random semiautomaton has
roughly $\alpha_kn$ reachable states with high probability.
As in
\cite{Trakhtenbrot73},
we refer to the number of reachable states as the {\em accessibility spectrum} of the automaton.
\bethn
\label{thm:Rn}
Let
$R_n\k$
be
the
accessibility spectrum
of a
random
semiautomaton
on $n$ states and $k$ symbols
drawn uniformly from
$[n]^{[n]\times[k]}$.
Then, for every fixed $k\ge2$ and as $n\to\infty$,
\begin{eqnarray}
\label{eq:Rndev}
\pr{\abs{R_n\k-\alpha_k n}>\sqrt n\log n
}
&=&
O\paren{\frac{1}{n^k}}.
\end{eqnarray}
\enthn

Second, we show that with high probability, very few states are lost
when equivalent ones are merged. Define $E_n\k$ to be the number
of ``excess'' reachable states:
$$E_n\k = R_n\k - \nrm{M_n\k}.$$
Note that in principle we need only show that
the number of states lost due to merging is small {\em after} the unreachables
are removed, but we will actually show that this is true even without removing them.

\bethn
\label{thm:Enk}
For every fixed $k\ge2$,
\begin{eqnarray*}
\pr{E_n\k >
C_k
\frac{
\log n}{\log\log n}
} = O\paren{\oo{n^k}}
\end{eqnarray*}
for an appropriate constant $C_k$.
\enthn
\begin{rem}
Theorem \ref{thm:Enk}
continues to hold
when each state is accepting with probability $0<p<1$ instead of $1/2$;
only the constants $C_k$ and those implicit in $O(\cdot)$
will change.
\end{rem}

\section{Proofs}
\belen
\label{lem:F}
Define the function
\begin{eqnarray*}
F(t) =
n(1-(1-1/n)^t)-(t-1)/2,
\qquad 1\le t\le n
.
\end{eqnarray*}
Then, for sufficiently large $n$,
\begin{eqnarray*}
\frac{F^2(t)}{t} &\ge&
\left\{
\begin{array}{lll}
0.01t,&& t\le n/2
,\\
\Omega(
\log^2n
)
,&& t\in[n/2,
\alpha_2n-\sqrt n\log n
]
\cup
[\alpha_2n+\sqrt n\log n,\infty)
.
\end{array}\right.
\end{eqnarray*}
\enlen
\bepf
We have $F(0)=1/2$, and for $t\le n/2$
\begin{eqnarray*}
F'(t) &=& -\oo2 +n\log\paren{1+\oo{n-1}}\cdot\paren{1-\oo n}^t \\
&\ge& -\oo2 + n\paren{\oo{n-1}-\oo{2(n-1)^2}}\cdot\paren{1-\oo n}^{n/2} \\
&\ge& -\oo2 + \frac{n}{n-1}\cdot\frac{2n-3}{2n-2}\cdot\paren{\oo{\sqrt e}-o(1)}\\
&\ge& \oo{\sqrt e}-\oo2-o(1)>0.1.
\end{eqnarray*}
This proves the estimate on $F^2(t)/t$ in the range $[1,n/2]$.

Now consider $t\in[n/2,
\alpha_2n-\sqrt n\log n
]$, and observe that $F(t)=H(t)+O(1)$, where
\begin{eqnarray*}
H(t)=n-t/2-n\exp(-t/n).
\end{eqnarray*}
By the definition of $\alpha_2$,
we have
$H(\alpha_2n)=H(0)=0$.
Furthermore, $H''(t)=-\exp(-t)/n<0$, and so $H$ is concave
with $H'(n\log2)=0$, and therefore increasing on
$[0,n\log2]$ and decreasing on $[n\log2,\infty)$.
Hence, to lower-bound $H^2(t)/t$ in the given range, it
suffices to estimate $H$
at its right endpoint:
\begin{eqnarray*}
H(\alpha  n-\sqrt n\log n) &=& \oo2\sqrt n\log n-e^{-\alpha }\sqrt n\log n+O(\log^2 n)
=\Omega(\sqrt n\log n)
.
\end{eqnarray*}
Since $H'(t)<1/2-e^{-\alpha }$ for $t>\alpha n$, we have
$H(\alpha_2  n+x)=\Omega(x)$ for $x>0$,
which completes the proof.
\enpf

\bepf[Proof of Theorem \ref{thm:Rn}]
We will prove the theorem for $k=2$;
the general case is completely analogous --- only the
constants implicit in $O(1/n^k)$ will vary with $k$.
For readability, we will write $\alpha=\alpha_2
$
and
$R_n=R_n^{(2)}$.
It will be convenient to embed
$R_n
$
in a slightly more general random process.
Fix $n\ge1$, and
define the
sequence of
random variables $(\nu_t)_{t=1}^\infty$,
as follows:
\begin{eqnarray*}
\nu_1 &=& 1,\\
\nu_{t+1} &=&
\left\{
\begin{array}{ll}
\nu_t,& \text{with probability } \nu_t/n
,\\
\nu_t+1,& \text{with probability } 1-\nu_t/n.
\end{array}\right.
\end{eqnarray*}
Clearly,
$\nu_t$ is with probability 1 nondecreasing,
upper-bounded by $n$, and reaches $n$ after a finite
number of steps.
Let us also define
\begin{eqnarray}
\label{eq:ON}
\omega_t
=
2\nu_t
+
1-t,
\qquad t\ge1.
\end{eqnarray}

Now consider the
following process for generating random
directed multigraphs with regular outdegree 2.
For time steps $t=1,2,\ldots$,
we will maintain the set of nodes $N_t$,
reached from $q_0=1$ by time $t$,
and two sets of edges: {\em open edges} $O_t$
and {\em closed edges} $C_t$.
A closed edge $c$ is an ordinary directed arrow from a source node $p$ to a
destination node $q$ marked with a $\sigma\in[k]$ and
denoted by $c=(p\wed{\sigma}q)$.
An open edge $o$ has a specified source $p$ but an as yet unspecified destination;
such an edge will be denoted by $o=(p\wed{\sigma}\star)$.
We initialize $N_1=\set{1}$,
$C_1=\emptyset$
and $O_1=\set{(1\wed{0}\star),(1\wed{1}\star)}$.
At time $t+1$, some (arbitrarily chosen\footnote{It is easy to see
that the distribution of $N_t,C_t,O_t$ is unaffected by the order
in which the open edges are selected.})
open edge in $o
\in O_t$
(if one exists)
chooses a destination node $q$
as follows:
\begin{itemize}
\item[(i)]
$q\in N_t$
with probability $|N_t|/n$
(that is, $o$ will point to a previously reached node);
\item[(ii)]
$q\in[n]\setminus N_t$
with probability $1-|N_t|/n$.
\end{itemize}
In event (i), $O_{t+1}=O_t\setminus\set{o}$,
while
in event (ii), $O_{t+1}=(O_t\setminus\set{o})\cup
\set{(q\wed{0}\star),(q\wed{1}\star)}$;
in both cases, $N_{t+1}=N_t\cup\set{q}$
and $C_{t+1}=C_t\cup\set{o}$.

The random semiautomaton embeds into the process $(\nu_t,\omega_t)$
via the following natural
correspondence: $|N_t|=\nu_t$ and $|O_t|=\omega_t$ as long as
 the latter is nonnegative (in particular,
the correspondence breaks down
for $t>2n+1$,
since
$\omega_t$ becomes negative).
Let $\tau$ be
the smallest $t$ for which $\omega_t=0$
---
i.e., the first time
there are no longer any open edges to choose from.
Then the pair $(N_\tau,C_\tau)$ defines\footnote{
Since the quantity of interest is the
accessibility spectrum,
it is unnecessary to define transitions out of unreachable states.
}
a semiautomaton with
accessibility spectrum
$R_n=\nu_\tau$, drawn uniformly from
$[n]^{[n]\times\set{0,1}}$.
Hence, proving (\ref{eq:Rndev}) amounts to showing that
\begin{eqnarray}
\label{eq:taudev}
\pr{\abs{\tau-2\alpha  n}>
\sqrt n\log n
} &=& O\paren{\frac{
1
}{n^2}}.
\end{eqnarray}
Indeed, (\ref{eq:taudev}) implies that $\tau=(2+o_P(1))\alpha  n$,
and $\nu_\tau=(1+o_P(1))\nu_{\alpha n}$.
Since,
by definition, $\tau$ is the smallest $t$ for which
$\nu_t=(t-1)/2$,
we have
\begin{eqnarray}
\nonumber
\pr{\tau\in[a,b]}
&\le&
\pr{\exists t\in[a,b
]: \nu_t=(t-1)/2} \\
\label{eq:tau[a,b]}
&\le&
\pr{\exists t\in[a,b
]: \nu_t\le(t-1)/2}
.
\end{eqnarray}
We estimate the left tail of $\tau$ as follows:
\begin{eqnarray*}
\pr{\tau\le\alpha n-
\sqrt n\log n
} &\le&
P_0 + P_1 + P_2,
\end{eqnarray*}
where
\begin{eqnarray*}
P_0 &=& \pr{\tau\in[1,150\log n]},\\
P_1 &=& \pr{\tau\in[150\log n,n/2]},\\
P_2 &=& \pr{\tau\in[n/2,\alpha n-
\sqrt n\log n
]}.
\end{eqnarray*}
To bound $P_0$, we
first argue, by elementary combinatorics, that $\pr{\omega_3<3}=O(1/n^2)$.
Now we condition on the high-probability event that there are at least $3$ open arrows available
after $3$ steps.
If all of the open arrows have been exhausted between time $t=4$ and $t=T$, then certainly at least
three of these arrows must point back to the $O(T)$ previous states. Thus, for $T=150\log n$,
\begin{eqnarray*}
P_0 \in
O\paren{\frac{1}{n^2}
+
{\binom{T}{3}}\paren{\frac{T}{n}}^3}
\subset
O\paren{\frac{1}{n^2}}
.
\end{eqnarray*}

To bound $P_1$ and $P_2$, we observe that an alternate interpretation
is possible for $\nu_t$.
Namely,
when $t$ balls are thrown into $n$ bins
uniformly at random,
the number of non-empty bins
is distributed as $\nu_t$.
We also observe that
\begin{eqnarray}
\label{eq:Enut}
\mathbb{E}\nu_t = n(1-(1-1/n)^{t})
,
\qquad t\ge1
.
\end{eqnarray}
Now by the Chernoff bound for negatively associated random variables
\cite[Prop. 5, Thm. 13]{Dubhashi:1998:BBS:299633.299634},
\begin{eqnarray}
\label{eq:duvshani}
\pr{\nu_t-\mathbb{E}\nu_t \le
-\Delta
}
&\le& \exp(-2\Delta^2/t),
\qquad \Delta>0.
\end{eqnarray}
Hence,
\begin{eqnarray*}
P_1 &\le& \sum_{t=150\log n}^{n/2} \pr{\nu_t-\mathbb{E}\nu_t\le -F(t)
},\\
P_2 &\le& \sum_{t=n/2}^{\alpha n-
\sqrt n\log n
} \pr{\nu_t-\mathbb{E}\nu_t\le -F(t)},
\end{eqnarray*}
where $F(t)$
is as in Lemma \ref{lem:F}.
The estimates in the lemma and (\ref{eq:duvshani})
yield
\begin{eqnarray*}
P_1 &\le& \frac{n}{2}\exp(-3\log n) \in O\paren{\oo{n^2}}
\end{eqnarray*}
and
\begin{eqnarray*}
P_2 &\in& O(n)\exp(-\Omega(\log^2 n))
\subset O\paren{\oo{n^2}}.
\end{eqnarray*}
This proves the left-tail estimate in (\ref{eq:taudev}).
To prove the corresponding right-tail estimate,
we observe that, analogously to (\ref{eq:tau[a,b]}),
\begin{eqnarray*}
\pr{\tau\in[a,b]}
\le
\pr{\exists t\in[a,b
]: \nu_t\ge(t-1)/2}
.
\end{eqnarray*}
The deviation probability is bounded as in (\ref{eq:duvshani}):
\begin{eqnarray*}
\pr{\nu_t-\mathbb{E}\nu_t\ge
\Delta
}
&\le& \exp(-2\Delta^2/t),
\qquad \Delta>0.
\end{eqnarray*}
Hence
\begin{eqnarray*}
\pr{\tau>\alpha n+\sqrt n\log n
} &\le& \sum_{t=\alpha n+\sqrt n\log n}^n \pr{\nu_t-\mathbb{E}\nu_t\ge G(t)},
\end{eqnarray*}
where $G(t)=-F(t)$.
Invoking again Lemma \ref{lem:F}, we have
\begin{eqnarray*}
\pr{\tau>\alpha n+\sqrt n\log n
} &\in&
O(n)
\exp(-\Omega(\log^2 n))
\subset O\paren{\oo{n^2}}.
\end{eqnarray*}
\enpf

\bepf[Proof of Theorem \ref{thm:Enk}]
Again, for
ease of exposition, we only prove the claim for $k=2$.
We start by explaining the idea of the proof. We need to show that there are usually ``few'' pairs of equivalent states.
Let us start by describing two ``typical'' situations in which equivalent states emerge.
The first is where a state is mapped into itself by every member of $\set{0,1}$.
Two such states are equivalent if and only if both are accepting or both are rejecting, which happens with a probability
of $1/2$. More generally, if from each of the two states one can reach very few states, then there is a non-negligible probability
that the states are equivalent. Thus, we will show that there are few states
with small accessibility spectra.
In the preceding sentence,
``few'' means (with high probability) ``at most $2$'', while ``small'' means ``less than $4\log_2n$''.

The second principal reason for two states $q,q'$ to be equivalent is that
$\delta(q,0)=\delta(q',0)$
and
$\delta(q,1)=\delta(q',1)$.
Again, $q$ and $q'$ are equivalent in this case with probability $1/2$.
Thus, we will need to show that there are few pairs of states $q,q'$
for which there are few words in $\set{0,1}^*$ taking $q$ and $q'$ to distinct states.
Here, the first ``few'' means ``at most $C\log n/\log\log n$'' and the second means ``up to $4\log_2 n$''.

Let us now consider the above scenarios in more detail.
The (random) set of states reachable from $q$ is given by $\set{q,\delta(q,0),\delta(q,1),\delta(q,00),\delta(q,01),\ldots}$.
Thus, the states reachable from $q$ reside on a binary tree whose edges are marked by letters in $\set{0,1}$.
Each time the random DFA selects a state $p=\delta(q,w)$, if $p$ is already in the tree, the edge that would create a directed cycle is not drawn.
We refer to the resulting tree as the tree {\em growing} from $q$.
Its size is the {accessibility spectrum} of $q$, denoted by $S(q)$.

Let $C>0$ be a constant to be determined later.
A state's accessibility spectrum is said to be {\em small} if it is below $C\log_2n$.
As in the
proof of Theorem \ref{thm:Rn},
the probability of a given state having a small accessibility spectrum is $O(1/n^2)$.
A similar argument shows that the joint probability of any pair of states
$q,q'$
having small accessibility spectra is $\tilde O(1/n^4)$.
Indeed,
consider the event of $S(q')$ being small, conditioned on $S(q)$ being such.
Draw the states ${\delta(q',0),\delta(q',1),\delta(q',00),\ldots}$ similarly
to the proof of Theorem \ref{thm:Rn}. The event in question is contained in
the event whereby, in the course of the first $C\log_2n$ steps of the process
of ``closing'' the open edges, we encounter at least twice
either a state visited already or a state belonging to the tree growing from $q$.
The probability of the latter event is clearly $O(\log^4/n^2)$.
Hence,
\begin{eqnarray*}
\pr{S(q),S(q')\text{ are both small}}
=O\paren{\frac{\log^4n}{n^4}}.
\end{eqnarray*}
Carrying this line of reasoning over to triples,
we have that the probability of any
three states having small accessibility spectra is $\tilde O(1/n^6)$
--- and therefore,
\begin{eqnarray*}
\pr{\text{there are 3 distinct states with small accessibility spectra}}
&\in &
\tilde O\paren{\oo{n^6}
\binom{n}{3}}\\
&=&\tilde O\paren{\oo{n^3}}
\subset O\paren{\oo{n^2}}.
\end{eqnarray*}

In view of the discussion above, we may assume (after removing up to 2 states)
that all states have large accessibility spectra. Consider two states $q,q'$. Let $T$ be a tree
of size $m=C\log_2n$ growing from $q$ (this will typically be a subtree of a larger tree of size $O(\alpha n)$).
The nodes of $T$ are ${\delta(q,w_1),\delta(q,w_2),\ldots,\delta(q,w_m)}$ for certain words
$w_1,w_2,\ldots,w_m\in\set{0,1}^*$.
If $\delta(q,w_i)\neq\delta(q',w_i)$,
$i=1,2,\ldots,m$, then the probability of $q,q'$ being equivalent is at most $1/2^m=1/n^C$.
(Note that this holds even if the states $\delta(q',w_i)$, $1\le i\le m$, are not
mutually distinct, in fact even if they all coincide. Similarly,
it does not matter if some of the states $\delta(q',w_i)$ coincide with some of the
$\delta(q,w_j)$, as long as $i\neq j$.)
The probability that both $\delta(q,0)=\delta(q',0)$
and
$\delta(q,1)=\delta(q',1)$ is $1/n^2$.
Call a state pair satisfying these equalities a {\em dud}.
The union bound does not yield a non-trivial upper bound on the number of duds,
and a more refined analysis will be needed.
Clearly, the probability that $d$
specific
pairwise disjoint
pairs are duds is $1/n^{2d}$.
Now the probability that there exist $d$ disjoint duds is at most
\begin{eqnarray*}
\label{eq:duds}
\oo{n^{2d}}\binom{n}{2d}\cdot(2d-1)\cdot(2d-3)\cdot\ldots\cdot1
\le
\oo{n^{2d}}\cdot \frac{n^{2d}}{(2d)!}\cdot\frac{(2d)!}{2^dd!}=\oo{2^dd!}
.
\end{eqnarray*}
Choosing $d=3\log n/\log\log n$
and applying Stirling's formula, we see that
the probability
of there being $d$ disjoint duds is $O(1/n^2)$.

Other than duds --- pairs ``dying'' right away after 2 steps ---
we must consider pairs dying after $4, 6,\ldots,C\log_2n$ steps. However,
the probability that a pair will die after 4 steps is $O(1/n^3)$, after
6 steps --- $O(1/n^4)$, and so forth. Hence, the probability
that there will be 2 pairs for which the process dies after 4 steps is
\begin{eqnarray*}
\binom{n}{4}\cdot O\paren{\frac{1}{n^6}}= O\paren{\oo{n^2}},
\end{eqnarray*}
that there will be a pair that dies after 6 steps,
\begin{eqnarray*}
\binom{n}{2}\cdot O\paren{\frac{1}{n^4}}= O\paren{\oo{n^2}},
\end{eqnarray*}
and that some pair will die after
$t\in[8,4\log_2n]$ steps,
\begin{eqnarray*}
\binom{n}{2}\cdot\tilde O\paren{\frac{1}{n^5}}= \tilde O\paren{\oo{n^3}}\subset O\paren{\oo{n^2}}.
\end{eqnarray*}
Now for two states $q,q'$
reaching distinct states for
many words $w_i$,
the probability of being equivalent is at most $1/n^C$.
Thus, it suffices to take $C=4$
to bound the probability of any pair of states
growing large trees yet being equivalent by
\begin{eqnarray*}
\binom{n}{2}\cdot O\paren{\frac{1}{n^4}}= O\paren{\oo{n^2}}.
\end{eqnarray*}
\enpf

\bepf[Proof of Theorem \ref{thm:main}]
Follows immediately from Theorems \ref{thm:Rn} and \ref{thm:Enk}
since the former estimates the number of
states remaining after
{\tt REMOVE-UNREACHABLE}
and the latter bounds the number of states lost after
{\tt COLLAPSE-EQUIVALENT}.
\enpf

\bibliographystyle{plain}


\section{Introduction}

The process of lepton-nucleon deep-inelastic scattering (DIS)
is a clean source of basic information about the hadron substructure in terms
of the parton model. Moreover, the QCD corrections to the parton model
provide the connection of the DIS structure functions with the parameters
of the QCD Lagrangian, in particular to the strong coupling $\alpha_s$
and the heavy-quark masses.
The higher order QCD corrections are manifest in the scaling violations
of the structure functions w.r.t. the virtual photon momentum transfer $Q^2$.
This phenomenon was observed shortly after the discovery of the partonic structure
of the nucleon and provided one of the first constraints on $\alpha_s$.
With the dramatic improvement in the accuracy of the lepton-nucleon
DIS data and the progress in the theoretical calculations the value of $\alpha_s$
can in principle be determined with an accuracy of $O(1\%)$~\footnote{For a recent overview
on precision determinations of $\alpha_s(M_Z^2)$ see \cite{Bethke:2011tr}.}.
The scaling violations, however, are also sensitive to the parton
distribution functions (PDFs). Therefore the determination of $\alpha_s$
has to be performed simultaneously with the nucleon PDFs in global fits.
Furthermore, an elaborate theoretical description of the
QCD scaling violations is available for the leading-twist terms only.
In practice this requires a careful isolation of the higher-twist effects
and/or their independent phenomenological parameterization.
Another important aspect of DIS phenomenology is related to the
$c$- and $b$-quark contributions. The heavy-quark production cross
section is sensitive to the heavy-quark masses. Therefore the DIS
data provide a constraint on the $c$- and $b$-quark masses, $m_{c,b}$.
The structure functions of the semi-inclusive process with the heavy
quark in the final state are particularly useful for this purpose, although
the data on inclusive structure functions are competitive with the
semi-inclusive ones due to a much better accuracy. A major pitfall arising
in the analysis of heavy quark production is related to the account of
the higher-order QCD corrections. Due to the two scales appearing in the problem
the calculations are quite involved. Therefore the NNLO corrections
to the heavy-quark lepto-production are known in partial form
only~\cite{Kawamura:2012cr} at present. The problem of the high-order
corrections is bypassed in the so-called variable-flavor-number (VFN) scheme
assuming zero mass for the $c$- and $b$-quark. In this approximation
the available high-order massless DIS Wilson coefficients can be
employed for the calculation of the heavy-quark lepto-production rates.
As well-known, the VFN approximation is obviously inapplicable  at scales
$Q^2 \sim m_{c,b}^2$ and it is commonly supplemented by {\it modeling} of the Wilson
coefficient at low $Q^2$ in order to arrive in this way at a general-mass VFN (GMVFN)
scheme. In the present paper we essentially focus on the determination of
$\alpha_s$ and $m_c$ based on the fixed-flavor-number (FFN) scheme. Here
the mass effects are taken into account on field-theoretic grounds,
free of model ambiguities appearing in the so-called GMVFN schemes. Moreover, we employ
the massive Wilson coefficients derived using the running-mass definition, which provide
improved perturbative stability of the heavy-quark production rate~\cite{Alekhin:2010sv}.

The paper is organized as follows.
In Section~\ref{sec:basics} we outline the theoretical basis of the
analysis and describe the data used. Sections~\ref{sec:alphas} and \ref{sec:mass} contain
our results on the determination of $\alpha_s$ and the $c$-quark mass,
respectively. In Section~\ref{sec:vfn} we compare the FFN and
VFN schemes with particular emphasis on the uncertainties
in the determination of $\alpha_s$ and $m_c$.

\section{Theoretical and Experimental Ingredients of the Analysis}
\label{sec:basics}

Our determination of $\alpha_s$ and $m_c$ is based on the QCD analysis
of the DIS data obtained in fixed-target experiments and at the HERA
collider.
Only the proton- and deuteron-target samples are selected which
allows to minimize the impact of the nuclear corrections on the results~\footnote{For a
detailed description
of the data set used cf. Ref.~\cite{Alekhin:2012ig}.}.
 The DIS
data are combined with the ones on the fixed-target Drell-Yan process,
providing a supplementary constraint on the PDFs and
to facilitate the separation of the valence and sea quark distributions.
\begin{figure}[pt]
\centerline{\psfig{file=scans.eps,width=4.25in}}
\vspace*{8pt}
  \caption[]{
    The $\chi^2$-profile versus the value of $\alpha_s(M_Z^2)$,
    for the separate data subsets, all obtained in variants of
    the ABM11 analysis with the value of $\alpha_s$ fixed and
    all other parameters fitted (solid lines: NNLO fit, dashes: NLO fit);
    from Ref.~\protect\cite{Alekhin:2012ig}.}
  \protect\label{fig:scans}
\end{figure}
The main version of our analysis is performed at NNLO using
the three-loop anomalous dimensions in the PDF evolution
and corresponding Wilson coefficients
for the light-flavor DIS structure functions and the Drell-Yan process.
For the neutral-current (NC) heavy-quark contribution we employ the
approximate NNLO Wilson coefficients~\cite{Kawamura:2012cr}.
These terms were derived combining results obtained with the soft-gluon resummation
technique and the high-energy limit of the DIS structure
functions~\cite{Catani:1990eg}. These two
approaches provide a good approximation
at kinematics close to threshold of the heavy-quark production
and far beyond the threshold, respectively. Between these two
regimes the constraint coming from the available NNLO
massive operator-matrix-element (OME) Mellin moments~\cite{Bierenbaum:2009mv}
are employed.
A remaining uncertainty in the NNLO Wilson coefficient obtained in this way is
quantified by its margins, A and B. To find the best shape of the NNLO term
preferred by the data we use a linear interpolation between
these margins
\begin{equation}
  c_{2}^{\,(2)} \,=\,
  (1-d_N) c_{2}^{\,(2),A} + d_N c_{2}^{\,(2),B}
  \, .
\label{eq:inter}
\end{equation}
and fit the interpolation parameter $d_N$ to the data simultaneously with the
PDF parameters, $\alpha_s$ and $m_c$.
The charged-current (CC) heavy-quark production is calculated with account
of the NLO corrections~\cite{Gottschalk:1980rv,Gluck:1996ve,Blumlein:2011zu},
which are the highest-order ones presently available.
Both NC and CC massive Wilson coefficients used in our analysis are
used with the $\overline{\rm MS}$-definition of the heavy-quark mass.
If compared to the case using the pole-mass the present choice is perturbatively
more stable \cite{Alekhin:2010sv}. With the running-mass definition this contribution
basically vanishes since the mass can be defined at the typical renormalization/factorization
scale of the process considered.

The leading-twist terms provided by the QCD-improved parton model
are not sufficient at small $Q^2$ and/or final-state hadronic mass $W$, where
parton correlations cannot be neglected. To account for these
we add on the top of the leading-twist term the
twist-4 contribution to the DIS structure functions $F_{2,T}$,
parameterized in a  model-independent form using
spline interpolation~\footnote{The twist-6 terms were also checked in the
fit and found comparable to zero within errors
applying the cut of $Q^2 > 2.5~{\rm GeV}^2$ used in the present analysis.}.
The twist-4 spline
coefficients are fitted to the data together with other the parameters.
This is particularly important for the case of $\alpha_s$ in view of
its strong correlation with the higher twist terms.

\section{Strong Coupling Constant}
\label{sec:alphas}

\begin{figure}[pt]
\centerline{\psfig{file=comp.eps,width=3.75in}}
\vspace*{8pt}
  \caption{The values of $\alpha_s(M_Z^2)$ obtained in the ABM
   fit~\protect\cite{Alekhin:2012ig} in comparison to the ones obtained
   by JR~\protect\cite{JRnew}, CTEQ~\protect\cite{Gao:2013xoa},
   MSTW~\protect\cite{Martin:2009bu}, and
   NNPDF~\protect\cite{Ball:2011us}
   from the analysis of the DIS data (squares)
   and from the combination of the DIS and
   jet Tevatron data~\protect\cite{Abulencia:2007ez,Abazov:2008hua}
   (circles).
}
  \protect\label{fig:comp}
\end{figure}
The strong coupling constant $\alpha_s$ can be determined
by comparing the $Q^2$-dependence of
the DIS cross section measurements with the predictions based on the
QCD-improved paton model. In our analysis the value of $\alpha_s$ is
obtained simultaneously with the nucleon PDFs and the twist-4 terms.
This allows to take into account the correlation of $\alpha_s$ with
other parameters affecting the data $Q^2$-dependence.
The central value of $\alpha_s$ obtained in
this way depends on the perturbative order and reduces from NLO to NNLO.
In particular, the ABM11 fit~\cite{Alekhin:2012ig} yields
\begin{eqnarray}
  \label{eq:alphas-nlo}
  \alpha_s(M_Z^2) \,\,=&
  0.1180\, \pm 0.0012 (\text{exp})\,
  \hspace*{30mm}
  &{\rm NLO} \nonumber
  \, ,
  \\
  \label{eq:alphas}
  \alpha_s(M_Z^2) \,\,=&
  0.1134\, \pm 0.0011 (\text{exp})\,
  \hspace*{30mm}
  &{\rm NNLO}
  \, .
\end{eqnarray}
The values of $\alpha_s$ preferred by each particular DIS data set are
demonstrated in Fig.~\ref{fig:scans} by means of the $\chi^2$-profiles
obtained in the variants of the ABM11 fit with $\alpha_s$ fixed
at the values in the range of $0.104\div0.130$. The HERA and BCDMS
data sets have a similar $\chi^2$-shape with minima around the
values Eq.~(\ref{eq:alphas}), while the SLAC and NMC data pull
the value of $\alpha_s$ somewhat up and down, respectively.
Note that the two latter sets
are sensitive to higher-twist terms due to
substantial small-$Q^2$ contributions  in these samples. In contrast,
the HERA and BCDMS data are far less sensitive to higher twists terms
and it is worth noting that in the variant of
the ABM11 fit excluding the SLAC and NMC data and setting the higher twist
terms to zero we find $\alpha_s(M_Z^2)=0.1133\, \pm 0.0011 (\text{exp.})$ at NNLO.
The good agreement of this value with Eq.~(\ref{eq:alphas})
substantiates the consistency between different data sets in our analysis
once the higher twist terms are taken into account. Moreover, this
cross-check confirms that the combination of the BCDMS and HERA data
can be used for accurate determination of $\alpha_s(M_Z^2)$ since these two data
sets provide complimentary constraints on the PDFs~\cite{Adloff:2000qk},
see also~\cite{Alekhin:2013kla}.

The NNLO value of $\alpha_s$ Eq.~(\ref{eq:alphas}) is in a good agreement
with the results of the JR analysis~\cite{JRnew} and the recent
CTEQ determination~\cite{Gao:2013xoa},
while the MSTW~\cite{Martin:2009bu} and NNPDF~\cite{Ball:2011us}
groups report substantially bigger values,
cf. Fig.~\ref{fig:comp}.
The discrepancy with MSTW can be explained in part by impact
of the jet Tevatron data, which pull the value of $\alpha_s(M_Z^2)$ up
by $0.001\div 0.002$, depending on the fit details.
Furthermore, changing our fit ansatz in direction of the MSTW and NNPDF ones we
approach their value of $\alpha_s$. In particular,
dropping the higher twist terms simultaneously
with an additional cut of $W^2>12.5~{\rm GeV}^2$
imposed by MSTW and NNPDF we
obtain $\alpha_s(M_Z^2)=0.1191\, \pm 0.0006 (\text{exp.})$. In a similar way,
disregarding the error correlation in the HERA and NMC data, like
in the MSTW analysis, we obtain $\alpha_s(M_Z^2)$ shifted by
$+0.0026$, in the direction of the MSTW value.
Note that in this context the recently updated JR analysis~\cite{JRnew}
treats the  higher twist properly.
The CTEQ analysis~\cite{Gao:2013xoa} seems to be
less sensitive to the impact of the higher twist
contributions, compared to that by MSTW and NNPDF due to a more stringent cut on
$Q^2$.

\section{The Mass of the Charm Quark}
\label{sec:mass}
\begin{figure}[pt]
\centerline{\psfig{file=pull.eps,width=4.25in}}
\vspace*{8pt}
\caption{\protect\label{fig:herac}
        The combined HERA data on the reduced cross section for the open
  charm production~\protect\cite{Abramowicz:1900rp} versus $x$ at different
  values of $Q^2$ in comparison with the analysis of~\protect\cite{Alekhin:2012vu} at
  NLO (dashed line) and NNLO (solid line) together with
  a fit variant based on the option (A+B)/2 of the NNLO Wilson
  coefficients of Ref.~\protect\cite{Kawamura:2012cr},
  cf. Eq.~(\protect\ref{eq:inter}) (dotted line);
  from Ref.~\protect\cite{Alekhin:2012vu}.
}
\end{figure}
The sensitivity to the charm quark mass $m_c$ in our analysis appears
essentially due to the data on
the NC and CC inclusive DIS~\cite{Aaron:2009aa}, and
the CC semi-inclusive charm production in DIS~\cite{Bazarko:1994tt,Goncharov:2001qe}
with the most essential experimental constraint on $m_c$
coming from the semi-inclusive charm-production
HERA data~\cite{Abramowicz:1900rp}.
The latter sample comprises the statistics of the H1 and ZEUS experiments
obtained for different $c$-quark decay channels. The combination was
performed similarly to the case of the inclusive HERA data~\cite{Aaron:2009aa}
and allows to reduce the systematic errors of each experiment due to cross-calibration
of the experiments. The FFN scheme provides a good description of the
semi-inclusive HERA data up to the largest $Q^2$-values covered,
cf. Fig.~\ref{fig:herac}, with a value of $\chi^2/NDP=61/52$
at NNLO. Here  NDP denotes the
number of data points.
The $\overline{\rm MS}$-values of $m_c$ found in the analysis of~\cite{Alekhin:2012vu} are
\begin{eqnarray}
  \label{eq:mcres-nlo}
  m_c(m_c) \,\,=&
  1.15\, \pm 0.04 (\text{exp})\,^{+0.04}_{-0.00} (\text{scale})
  \hspace*{30mm}
  &{\rm NLO}
  \, ,
  \\
  \label{eq:mcres-nnlo}
  m_c(m_c) \,\,=&
  1.24\, \pm 0.03 (\text{exp})\,^{+0.03}_{-0.02} (\text{scale})\,^{+0.00}_{-0.07} (\text{th}),
  \hspace*{14mm}
  &{\rm NNLO_\text{approx}}
  \, ,
\end{eqnarray}
\begin{figure}[pt]
\centerline{\psfig{file=alpha.eps,width=2.75in}}
\vspace*{8pt}
\caption{\protect\label{fig:alpha}
  The values of $m_c(m_c)$ obtained in the NLO and NNLO variants of
  our analysis with the value of $\alpha_s(M_Z^2)$ fixed.
  The position of the star displays the result with the value of
  $\alpha_s(M_Z^2)$ fitted~\protect\cite{Alekhin:2012ig};
  from Ref.~\protect\cite{Alekhin:2012vu}.
}
\end{figure}
at NLO and NNLO, respectively, see also~\cite{Alekhin:2012un}.
The experimental accuracy of 30~MeV obtained at NNLO is quite competitive with other
determinations of $m_c$ based on the $e^+e^-$ data and the central value
Eq.~(\ref{eq:mcres-nnlo}) is in a good agreement with the world
average~\cite{Beringer:1900zz}.
The scale error in Eqs.~(\ref{eq:mcres-nlo}) and (\ref{eq:mcres-nnlo}) is obtained
varying the factorization scale by a factor of $1/2$ and $2$ around
the nominal value of $\sqrt{m_c^2+\kappa Q^2}$,
where $\kappa=4$ for NC and $\kappa=1$ for  CC
heavy-quark production, respectively.
In the NNLO case an additional error related to the uncertainty
in the massive Wilson coefficients contributes. The value of
Eq.~(\ref{eq:mcres-nnlo}) is obtained for the interpolation parameter
$d_N=-0.1$ being preferred by the fit, roughly corresponding to option A
of the Wilson coefficients of Ref.~\cite{Kawamura:2012cr}. Meanwhile,
option B is clearly excluded by the data with $\chi^2$/NDP=115/52.
Therefore the uncertainty due to the missing NNLO massive terms
is estimated as a variation between options A and (A+B)/2 which yields
the value of 70~MeV Eq.~(\ref{eq:mcres-nnlo}).
The value of $m_c(m_c)$ obtained in our analysis demonstrates remarkable
stability w.r.t. $\alpha_s(M_Z^2)$.
Performing variants of our analysis with the values of $\alpha_s$
fixed in the wide range around the best value preferred by the
data, we find a variation of $m_c(m_c)$ in the range of 10-20~MeV,
depending on the order, cf. Fig~\ref{fig:alpha}.

The NLO value of $m_c(m_c)$ obtained in our analysis is somewhat lower
than the one of $m_c(m_c)=1.26\pm0.05~(\text{exp})~{\rm GeV}$ from
the analysis based on the HERA data only~\cite{Abramowicz:1900rp}. To
understand the difference we checked the cases of a cut on
$Q^2=3.5~{\rm GeV}^2$, likewise in the HERA fit, and no
semi-inclusive Tevatron data~\cite{Bazarko:1994tt,Goncharov:2001qe}
included. As a result we obtain
shifts in $m_c(m_c)$ by $+30~\text{MeV}$ and $+40~\text{MeV}$, respectively.
The remaining discrepancy should be attributed to the particularities
of the HERA PDFs. The value of $m_c(m_c)$ was recently also determined by
the CTEQ collaboration~\cite{Gao:2013wwa}. In contrast to our case this
determination is based on the S-ACOT-$\chi$
prescription as GMVFN scheme. Furthermore, the $\overline{\rm MS}$
coefficients of Ref.~\cite{Gao:2013wwa} are obtained by straightforward
substitution of the pole- and running-mass matching relation
into the pole-mass coefficients. The expressions obtained in such a way
correspond to a mixed order in $\alpha_s$. Moreover, the advantage
of this approach is not evident in view of the poor perturbative convergence
of the mass matching relation. The central CTEQ value of
$m_c(m_c)=1.12 ^{+0.11}_{-0.17}~{\rm GeV}$ is lower than
the world average, while the CTEQ errors are much larger than those in
Eqs.~(\ref{eq:mcres-nlo}) and (\ref{eq:mcres-nnlo}) due to
the impact of the uncertainty in the GMVFN scheme modeling.
\begin{figure}[pt]
\centerline{\psfig{file=pdfder.eps,width=4.75in}}
\vspace*{8pt}
\caption{\label{fig:pdfder} The difference between the
$c$-quark PDFs derivatives
$\dot{c}(x,\mu^2)\equiv \frac{dc(x,\mu^2)}{d\ln\mu^2}$ calculated with
the FOPT matching condition and with the massless 4-flavor evolution
starting at the matching point $\mu_0=m_c=1.4~{\rm GeV}$ versus the
factorization scale $\mu^2$ at different values of $x$ in the LO,
NLO, and NNLO* approximations. The arrows display upper margin of the
HERA collider kinematics with the collision
c.m.s. energy squared $s=10^5~{\rm GeV}^2$ and the vertical lines correspond to
the matching point position $\mu_0$.}
\end{figure}
 %

\section{VFN Uncertainties}
\label{sec:vfn}
The choice of the factorization scheme plays an essential role in the analysis of
existing DIS data due to important constraints coming from the small-$x$ region,
where the heavy quark contribution is numerically large. While
our analysis is based on the FFN scheme, many other groups employ
different variants of the GMVFN scheme, which differ by {\it modeling} of the
low-$Q^2$ region. The spread between these variants
is rather substantial and thus implies a corresponding uncertainty in the basic parameters
determined in these GMVFN fits. In particular, the value of $m_c$ determined from a
combination of the inclusive and semi-inclusive HERA data with the
different versions of the
ACOT and RT prescriptions for the VFN scheme demonstrate a
spread of 400 MeV~\cite{Abramowicz:1900rp} .
There are also sources of the VFN scheme uncertainties, which are
common for all these prescriptions.
Firstly, the matching of heavy-quark PDFs is commonly performed
at the factorization scale $\mu$ equal to the
$\mu_0=m_{c}$, resp. $m_b$. Clearly, at these scales neither of the heavy flavors can
be dealt with as massless. The matching point $\mu_0$
is not fixed by theory
and in principle it can vary in a wide range being an artefact of the description
not contributing to the observables according to the renormalization group equations.
Additional
uncertainty emerge for the 4(5)-flavor PDFs
in the NNLO analysis. They are commonly matched with the
3(4)-flavor ones using {\it NLO matching condition},
since the NNLO OMEs are not yet known in the complete form~\footnote{
For progress in this field
cf.~\cite{Ablinger:2012ej,Ablinger:2012qm,Ablinger:2010ty,Blumlein:2012vq}.}.
However, to provide consistency with the NNLO Wilson coefficients
the evolution of these PDFs is performed in the NNLO approximation that
introduces an additional uncertainty due to missing higher-order
corrections into the analysis.
At the same time, the evolution of the 4(5)-flavor PDFs in the VFN scheme
leads to a resummation of the terms $\sim \ln(\mu^2)$
which in part reproduce the higher-order corrections, being known, however, not to be dominant.
Relations between the resummation effects and the VFN evolution
uncertainties are illustrated in Fig.~\ref{fig:pdfder} by comparison
of the $\mu$-derivatives for the $c$-quark distribution
being calculated in different ways. In one case the distributions
are matched at the scale of $\mu_0$ using the fixed-order-perturbation-theory
(FOPT) matching conditions and then evolved starting from $\mu_0$ with
the massless splitting functions. In another case they are calculated
with the FOPT matching conditions at all scales.
The difference between these two cases do not demonstrate a significant
rise with $\mu$.  The only exclusion is observed at $x\lesssim 0.0001$
and at scales outside of the kinematics being probed in experiment.
Therefore it cannot be attributed to the impact of the log-term resummation.
In contrast, there is a substantial difference in the derivatives
calculated in the NLO and NNLO*, i.e. the combination of the NLO
matching with the NNLO evolution. This difference yields an estimate of
the uncertainty in the VFN scheme due to the missing higher-orders, which
is obviously larger than the resummation effects.
Checking the impact of this uncertainty w.r.t. the value of $\alpha_s(M_Z^2)$
in combination with the variation of the matching
point for the 4-flavor PDFs in the range of $1.2\div 1.5~{\rm GeV}$
we find a value of $\pm 0.001 $. This is comparable to the experimental uncertainty
and makes the VFN schemes incompetitive with the FFN one in the
precision determination of $\alpha_s(M_Z^2)$.

\section*{Acknowledgments}

We thank P. Jimenez-Delgado and E. Reya for discussions. This work has been supported in part
by Helmholtz Gemeinschaft under contract VH-HA-101 ({\it Alliance Physics at the Terascale}),
DFG Sonderforschungsbereich/Transregio~9 and by the European Commission through contract
PITN-GA-2010-264564 ({\it LHCPhenoNet}).


\section*{\hfil #1\hfil}}
\bibliographystyle{abbrv}

\usepackage{appendix}

\newcommand{\Aightarrow}{\stackrel{a}{\Rightarrow}}

\newcommand{\comment}[1]{}

\usepackage{graphicx}

\usepackage{subfigure}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mdwlist}

\usepackage{algorithm}

\usepackage{hyperref}

\usepackage{xspace}
\usepackage{setspace}

\newcommand{\BSB}{{\tt Broadcast\_Binary}\xspace}
\newcommand{\DG}{{\tt Diag\_Graph}\xspace}
\newcommand{\TRUE}{{\tt TRUE}\xspace}
\newcommand{\FALSE}{{\tt FALSE}\xspace}

\newcommand{\teq}[1]{\stackrel{#1}{\equiv}}
\newcommand{\nteq}[1]{\stackrel{#1}{\nequiv}}

\newcommand{\sg}{{\mathcal G}}
\newcommand{\sv}{V}
\newcommand{\se}{{\mathcal E}}
\newcommand{\sw}{{\mathcal W}}
\newcommand{\ms}{{\mathcal S}}
\newcommand{\sn}{{\mathcal N}}

\newcommand{\nchoosek}[2]{{#1 \choose #2}}

\newcommand{\mybox}[1]{\vspace{5pt}\centerline{\framebox{\parbox[c]{\textwidth}{#1}}}\vspace{5pt}}

\newcommand{\BlackBox}{\rule{2.6mm}{2.6mm}}
\newenvironment{proof}{\noindent {\bf Proof:}~}{\hspace*{\fill}\(\Box\)}
\newenvironment{proofSketch}{\noindent{\bf Proof Sketch:}}{\hspace*{\fill}\(\Box\)}
\newtheorem{theorem}{Theorem}
\newtheorem{conjecture}{Conjecture}
\newtheorem{claim}{Claim}
\newtheorem{corollary}{Corollary}
\newtheorem{definition}{Definition}
\newtheorem{condition}{Condition}
\newtheorem{lemma}{Lemma}
\newtheorem{example}{Example}
\newcommand{\fillblackbox}{\hspace*{\fill}\(\BlackBox\)}
\newcommand{\fillbox}{\hspace*{\fill}\(\Box\)}
\newcommand{\fig}[1]{Figure~\ref{#1}}
\newcommand{\eqn}[1]{Equation~\ref{#1}}
\newcommand{\refsec}[1]{Section~\ref{#1}}
\newcommand{\num}[1]{(\romannumeral#1)}
\def\noflash#1{\setbox0=\hbox{#1}\hbox to 1\wd0{\hfill}}

\newcommand{\scriptf}{\mathcal{F}}
\newcommand{\scripte}{\mathcal{E}}
\newcommand{\scriptv}{\mathcal{V}}
\newcommand{\scriptl}{\mathcal{L}}
\newcommand{\scripts}{\mathcal{S}}

\newcommand{\matrixh}{\textbf{H}}
\newcommand{\graphh}{\textit{H}}
\newcommand{\bfF}{\textbf{F}}
\newcommand{\bfC}{\textbf{C}}
\newcommand{\bfR}{\textbf{R}}
\newcommand{\bfL}{\textbf{L}}

\newcommand{\shortdividerline}{\begin{center} \line(1,0){150} \end{center}}
\newcommand{\dividerline}{\begin{center}\hrule\end{center}}

\newcommand{\Zightarrow}{\rightarrow}

\newcommand{\propagate}[3]{#1\stackrel{#3}{\rightsquigarrow}{#2}}
\newcommand{\notpropagate}[3]{#1\stackrel{#3}{\not\rightsquigarrow}{#2}}

\newcommand{\StableVector}{{\tt StableVector}}
\newcommand{\ssss}{{\mbox{\it start}}}
\newcommand{\eeee}{{\mbox{\it end}}}

\newcommand{\vectorv}{{\normalfont\textbf{v}}}
\newcommand{\matrixm}{\textbf{M}}

\newcommand{\I}{{\mathcal I}}
\newcommand{\D}{{\textbf d_H}}
\newcommand{\M}{{\mathcal M}}
\newcommand{\R}{{\mbox{\it Verified\,}}}
\newcommand{\T}{{\mathcal T}}
\newcommand{\Q}{{\mathcal Q}}
\newcommand{\vmsg}{{\tt VerifyMsg}}
\newcommand{\HH}{{\mathcal H}}
\newcommand{\VV}{{\mathcal V}}
\newcommand{\Ha}{$H_a$}
\newcommand{\Hi}{$H_0$}
\newcommand{\Hl}{$H_L$}
\newcommand{\recv}{{\tt Recv}}
\newcommand{\send}{{\tt Send}}

\newcommand{\RBRecv}{{\tt RBRecv}}
\newcommand{\RBSend}{{\tt RBSend}}

\newcommand{\bfA}{{\bf A}}
\newcommand{\bfH}{{\bf H}}
\newcommand{\bfQ}{{\bf Q}}
\newcommand{\bfM}{{\bf M}}
\newcommand{\bfv}{v}

\begin{document}
\title{Byzantine Convex Consensus: Preliminary Version\footnote{\normalsize We present an optimal algorithm in our follow-up work \cite{Tseng_BCC_optimal}.}
\footnote{\normalsize This research is supported in part by National Science Foundation award CNS 1059540. Any opinions, findings, and conclusions or recommendations expressed here are those of the authors and do not necessarily reflect the views of the funding agencies or the U.S. government.}}

\author{Lewis Tseng$^{1,3}$, and Nitin Vaidya$^{2,3}$\\~\\
 \normalsize $^1$ Department of Computer Science,\\
 \normalsize $^2$ Department of Electrical and Computer Engineering,
 and\\ \normalsize $^3$ Coordinated Science Laboratory\\ \normalsize University of Illinois at Urbana-Champaign\\ \normalsize Email: \{ltseng3,
nhv\}@illinois.edu \\~\\ \normalsize Technical Report}

\date{July 3rd, 2013}
\maketitle

\begin{abstract}
{\normalsize

Much of the past work on asynchronous approximate Byzantine consensus has
assumed {\em scalar} inputs at the nodes \cite{AA_Dolev_1986, AA_nancy}.
Recent work has yielded approximate Byzantine consensus algorithms
for the case when the input at each node is a $d$-dimensional vector,
and the nodes must reach consensus on
a vector in the convex hull
of the input vectors at the fault-free nodes
\cite{herlihy_multi-dimension_AA,Vaidya_BVC}.
The $d$-dimensional vectors can be equivalently viewed as {\em points}
in the $d$-dimensional Euclidean space.
Thus, the algorithms in \cite{herlihy_multi-dimension_AA, Vaidya_BVC}
require the fault-free nodes to decide on a point
in the $d$-dimensional space.

~

In this paper, we generalize the problem to allow the decision to be
a {\em convex polytope} in the $d$-dimensional space, such that the decided
polytope is within the convex hull of the input vectors at the fault-free nodes.
We name this problem as {\em Byzantine convex consensus} (BCC), and present an
asynchronous approximate BCC algorithm with optimal fault tolerance.
Ideally, the goal here is to agree on a convex polytope that
is as large as possible.
While we do not claim that our algorithm satisfies this goal,
we show a bound on the output convex polytope chosen by our algorithm.

}
\end{abstract}

\thispagestyle{empty}
\newpage
\setcounter{page}{1}

\section{Introduction}
\label{s_intro}

Much of the past work on asynchronous approximate Byzantine consensus has
assumed {\em scalar} inputs at the nodes \cite{AA_Dolev_1986, AA_nancy}.
Recent work has yielded approximate Byzantine consensus algorithms
for the case when the input at each node is a $d$-dimensional vector,
and the nodes must reach consensus on
a vector in the convex hull
of the input vectors at the fault-free nodes
\cite{herlihy_multi-dimension_AA,Vaidya_BVC}.
The $d$-dimensional vectors can be equivalently viewed as {\em points}
in the $d$-dimensional Euclidean space.
Thus, the algorithms in \cite{herlihy_multi-dimension_AA, Vaidya_BVC}
require the fault-free nodes to decide on a point
in the $d$-dimensional space.
In this paper, we generalize the problem to allow the decision to be
a {\em convex polytope} in the $d$-dimensional space, such that the decided
polytope is within the convex hull of the input vectors at the fault-free nodes.
We name this problem as {\em Byzantine convex consensus} (BCC), and present an
asynchronous approximate BCC algorithm with optimal fault tolerance.

We consider Byzantine convex consensus (BCC) in an {\em asynchronous} system consisting of $n$ nodes,
of which at most $f$ may be Byzantine faulty. The Byzantine faulty nodes may behave
in an arbitrary fashion,
and may collude with each other. Each node $i$ has a $d$-dimensional vector of reals as its {\em input} $x_i$.
All nodes can communicate with each other directly on reliable and FIFO (first-in first-out)
channels.  Thus, the underlying communication graph can be modeled as a {\em complete graph},
with the set of nodes being $V=\{1,2,\cdots, n\}$.
The impossibility of {\em exact} consensus in asynchronous systems \cite{FLP_one_crash}
applies to BCC as well. Therefore, we consider
the {\em Approximate BCC} problem with the following requirements:

\begin{itemize}

\item \textbf{Validity}: The {\em output} (or {\em decision}) at each fault-free node must be
a convex polytope in the convex hull of the $d$-dimensional input vectors at the fault-free nodes.
(In a degenerate case, the output polytope may simply be a single {\em point}.)

\item \textbf{Approximate Agreement}: For any $\epsilon > 0$, the {\em Hausdorff distance} (defined below) between the output polytopes at any two fault-free nodes must be at most $\epsilon$.

\end{itemize}
Ideally, the fault-free nodes should reach consensus on the largest
possible convex polytope that satisfies the validity constraint. We present an optimal algorithm that agrees on a convex polytope that is as {\em large} as possible under adversarial conditions in our follow-up work \cite{Tseng_BCC_optimal}.
The motivation behind reaching consensus on a convex polytope is that
a solution to BCC is expected to also facilitate solutions to a large
range of consensus problems (e.g., Byzantine vector consensus \cite{herlihy_multi-dimension_AA,Vaidya_BVC}, or convex function optimization
over a convex hull of the inputs at fault-free nodes).
Future work will explore these potential applications.

To simplify the presentation,
we do not include a {\em termination} condition above.
We instead prove that approximate agreement condition is {\em eventually} satisfied
by our proposed algorithm (in addition to validity).
However, we can augment the proposed algorithm
 using techniques similar to those in \cite{abraham_04_3t+1_async,herlihy_multi-dimension_AA,Vaidya_BVC},
to terminate within a finite number of rounds.

\begin{definition}
\label{def:dist}
For two convex polytopes $h_1, h_2$, the {\em Hausdorff distance} is defined as \cite{Hausdorff}
\[
\D(h_1, h_2) ~~=~~ \max~~ \{~~ \max_{p_1 \in h_1}~\min_{p_2 \in h_2} d(p_1, p_2),~~~~ \max_{p_2 \in h_2}~\min_{p_1 \in h_1} d(p_1, p_2)~~\}
\]
\end{definition}
where $d(p,q)$ is the Euclidean distance between points $p$ and $q$.

\paragraph{Lower Bound on $n$:}
As noted above,
\cite{herlihy_multi-dimension_AA, Vaidya_BVC} consider the problem of
reaching approximate Byzantine consensus on a
vector (or a point) in the convex hull of the $d$-dimensional
input vectors at the fault-free nodes, and show that
$n \geq (d+2)f+1$ is necessary.
\cite{herlihy_colorless_async} generalizes
the same lower bound to colorless tasks.
The lower bound proof in
\cite{herlihy_multi-dimension_AA, Vaidya_BVC} also implies that
$n \geq (d+2)f+1$ is necessary to ensure that BCC is solvable.
We do not reproduce the lower bound proof here, but
in the rest of the paper, we assume that $n\geq (d+2)f+1$,
and also that $n\geq 2$ (because consensus is trivial when $n=1$).

\section{Preliminaries}
\label{s_ops}

Some notations introduced throughout the paper are summarized in Appendix \ref{app_s_notations}.
In this section,
we introduce operations $\HH$, $H_l$, $H$,
and a {\em reliable broadcast} primitive used later in the paper.
\begin{definition}
\label{def:hh}
Given a set of points $X$, $\HH(X)$ is defined
as the convex hull of the points in $X$.
\end{definition}
\begin{definition}
\label{def:linear_hull}
Suppose that $\nu$ convex polytopes $h_1, h_2,\cdots, h_\nu$, and $\nu$ constants $c_1, c_2, \cdots, c_\nu$
are given such that
(i) $0 \leq c_i \leq 1$ and $\sum_{i = 1}^\nu c_i = 1$,
and (ii) for $1\leq i\leq\nu$, if $c_i\neq 0$, then $h_i\neq\emptyset$.
Linear combination of these convex
polytopes, $H_l(h_1, h_2,\cdots, h_\nu;~c_1, c_2, \cdots, c_\nu)$,
 is defined as follows:
\begin{itemize}
\item Let $Q := \{ i ~|~ c_i\neq 0, ~ 1\leq i\leq \nu \}$.
\item
$p \in H_l(h_1, h_2,\cdots, h_\nu;~c_1, c_2,\cdots, c_\nu)$ if and only if
\begin{equation}
\label{eq:linear_hull}
\text{for each~} i\in Q,
\text{~there exists~} p_i \in h_i, ~~\text{such that}~~p = \sum_{i \in Q} c_i p_i
\end{equation}
\end{itemize}
\end{definition}
Note that a convex polytope may possibly consist of a single point.
Because $h_i$'s above are all convex, $H_l(h_1, h_2,\cdots, h_\nu;~c_1, c_2,\cdots, c_\nu)$ is also a
convex polytope (proof included in Appendix \ref{a:linear_convex} for completeness).
The parameters for $H_l$ consist of two lists, a list of polytopes
$h_1,\cdots,h_\nu$, and a list of weights $c_1,\dots,c_\nu$.  With an abuse of notation, we
will specify one or both of these lists as either a {\em row vector}
or a {\em multiset}, with the understanding that the {\em row vector} or
{\em multiset} here represents an ordered list of its elements.

Function $H$ below is called in our algorithm with parameters $(\VV,t)$ wherein
$t$ is a round index ($t\geq 0$) and $\VV$ is a set of tuples of the
form $(h,j,t-1)$, where $j$ is a node identifier; when $t=0$, $h$ is a point
in the $d$-dimensional Euclidean space, and when $t>0$, $h$
is a convex polytope.

\vspace*{2pt}

\noindent
\hrule
{\bf Function $H(\VV,t)$}
\begin{enumerate}
\item Define multiset $X:=\{h~|~ (h,j,t-1)\in \VV\}$.
In our use of function $H$, each $h\in X$ is always non-empty.

\item If $t=0$ then
${\tt temp}~:= ~ \cap_{\,C \subseteq X, |C| = |X| - f}~~\HH(C)$.\\
In our use of function $H$, when $t=0$, each $h\in X$ is simply a point.
The intersection above is over the convex hulls of the subsets of $X$ of size $|X| - f$.
\item If $t>0$ then ${\tt temp}~:=~ H_l(X; \frac{1}{|X|}, \cdots,\frac{1}{|X|})$.
~~Note that all the weights here are equal to $\frac{1}{|X|}$.
\item Return {\tt temp}.
\end{enumerate}
\hrule

\paragraph{{\em Reliable Broadcast}\, Primitive:}

We will use the {\em reliable broadcast} primitive
from \cite{abraham_04_3t+1_async}, which is also used in other related work \cite{herlihy_multi-dimension_AA, Vaidya_BVC}.
\begin{itemize}
\item
As seen later, our algorithm proceeds in asynchronous rounds, and the
nodes perform reliable broadcast of messages that each consist of
a 3-tuple of the form $(v,i,t)$: here $i$ denotes the sender node's identifier, $t$
is round index, and $v$ is message value (the value $v$ itself is often a tuple).
The operation $\RBSend(v,i,t)$ is used by node $i$ to perform {\em reliable broadcast} of $(v,i,t)$  in round $t$.
\item
When message $(v,i,t)$ is {\em reliably received} by some node $j$,
the event $\RBRecv(v,i,t)$ is said to have occurred at node $j$
 (note that $j$ may possibly be equal to $i$).
The second element in a reliably received message 3-tuple,
namely $i$ above, is always identical to the identifier
of the node that performed the corresponding reliable broadcast.
When we say that node $j$ reliably receives $(v,i,t)$ we mean that event $\RBRecv(v,i,t)$
occurs at node $j$.
\end{itemize}
Each fault-free node performs one reliable broadcast ($\RBSend$) in each round of our algorithm.
The reliable broadcast primitive has the properties listed below,
as proved previously \cite{abraham_04_3t+1_async,herlihy_colorless_async}.
\begin{itemize}
\item {\bf Fault-Free Integrity:} If a fault-free node $i$ {\em never} reliably broadcasts $(v, i, t)$,
then no fault-free node ever reliably receives $(v,i,t)$.

\item {\bf Fault-Free Liveness:} If a fault-free node $i$ performs reliable broadcast of $(v, i, t)$, then each fault-free node eventually reliably receives $(v, i, t)$.

\item {\bf Global Uniqueness:} If two fault-free nodes $i, j$ reliably receive $(v, k, t)$ and $(w, k, t)$, respectively, then $v = w$, even if node $k$ is faulty.

\item {\bf Global Liveness:} For any two fault-free nodes $i, j$, if $i$ reliably receives $(v, k, t)$, then $j$ will eventually reliably receive $(v, k, t)$, even if node $k$ is faulty.

\end{itemize}

\section{Proposed Algorithm: Verified Averaging}
\label{s_alg}

The proposed algorithm (named {\em Verified Averaging}) proceeds in asynchronous rounds.
The input at each node $i$ is a $d$-dimensional vector of reals, denoted as $x_i$.
 In each round $t~(t \geq 0)$, each node $i$ computes a state variable $h_i$, which represents a convex polytope in the $d$-dimensional Euclidean space. We will refer to the value of $h_i$ at the {\em end} of the $t$-th round performed by node $i$ as $h_i[t]$, $t\geq 0$. Thus, for $t\geq 1$, $h_i[t-1]$ is the value of $h_i$ at the {\em start} of the $t$-th round at node $i$.

Motivated by previous work that uses a mechanism to simulate omission failures in
presence of
Byzantine faults \cite{Welch_textbook}, our algorithm uses a similar technique,
named {\em verification}. Informally,
the {\em verification} mechanism ensures that
if a faulty node deviates from the algorithm specification
(except possibly choosing an invalid input vector),
then its incorrect messages will be ignored by
the fault-free nodes.
Thus, aside from choosing a bad input, a faulty node cannot cause any other damage
to the execution.
Before we present the algorithm, we introduce a convention for the brevity of presentation:
\begin{itemize}
\item When we say that $(*,i,t) \in \VV$, we mean that {\em there exists $z$ such that $(z,i,t)\in\VV$}.

\item When we say that $(*,i,t) \not\in \VV$, we mean that {\em $\forall z$, $(z,i,t)\notin\VV$}.
\end{itemize}
The proposed {\em Verified Averaging} algorithm for node $i\in V$ is presented below.
All references to line numbers in our discussion refer to numbers listed
on the right side of the algorithm pseudo-code.
Whenever a message is reliably received by any node, a handler is called
to process that message. Handlers for multiple reliably received messages
may execute {\em concurrently} at a given node. For correct behavior, lines 3-7 below
are executed {\bf atomically}, and similarly, lines 11-16 are executed
{\bf atomically}.

In the proposed algorithm, in round 0, each node $i$ uses $\RBSend$ to reliably broadcast
$(x_i,i,0)$ where $x_i$ is its
input (line 1).

Lines 2-7 specify the event handler for event $\RBRecv(x,j,0)$ at node $i$.
Whenever a new message of the form $(x,j,0)$ is reliably received by node $i$
(line 2),
the set $\R_i[0]$ is updated (line 3).
Note that the message received by node $i$ on line 2 may possibly have been reliably
broadcast by node $i$ itself.
When size of set
$\R_i[0]$ becomes at least $n-f$ for the first time (line 4),
node $i$ computes $h_i[0]$ (line 6); the $\R$ set used for
computing $h_i[0]$ is saved as $\R^c_i[0]$ (line 5).
Having computed $h_i[0]$, node $i$ can proceed to round 1 (line 7).
Note that new messages may still be added to $\R_i[0]$ afterwards
whenever event of the form $\RBRecv(x_j,j,0)$ occurs.
Thus, $\R_i[0]$ may continue to grow even after node $i$ has proceeded
to round 1; however, $\R^c_i[0]$ is not modified again.
Recall that lines 3-7 are performed atomically.

On entering
round $t$, $t\geq 1$, each node $i$ reliably broadcasts
$((h_i[t-1],\R^c_i[t-1]),i,t)$ (line 8).

Lines 9-16 specify the event handler for event $\RBRecv((h,\VV),j,t)$ at node $i$.
Whenever a message of the form $((h,\VV),j,t)$ is reliably received from node $j$ (line 10),
node $i$ first waits until its own set $\R_i[t-1]$ becomes
large enough to contain $\VV$. Note that $\R_i[t-1]$ is initially computed
in round $t-1$, but it may continue to grow even after node $i$ proceeds
to round $t$.
If the condition $\VV\subseteq\R_i[t-1]$ never becomes true, then this
message is not processed further.
The message $((h,\VV),j,t)$ is considered correct if all the following
conditions are true: (i) $\VV\subseteq \R_i[t-1]$,
(ii) $|\VV|\geq n-f$,
(iii) $(*,j, t-2)\in \VV$, and
(iv) $h=H(\VV,t-1)$.
Conditions (ii), (iii) and (iv) are tested in line 11,
and node $i$ does not reach line 11 until condition (i)
becomes true at line 10.
If the message $((h,\VV),j,t)$ is considered correct, then $(h,j,t-1)$ it is added to $\R_i[t]$ (line 12).

When both conditions $|\R_i[t]|\geq n-f$ and $(h_i[t-1],i,t-1)\in \R_i[t]$
are true for the first time (line 13),
node $i$ computes $h_i[t]$ (lines 14 and 15), and then proceeds to round $t+1$.
 Similar to round 0, the set used
in computing $h_i[t]$ is saved as $\R^c_i[t]$ (line 14).
While $\R^c[t]$ remains unchanged afterwards, $\R_i[t]$ may continue to grow, even
after node $i$ proceeds to round $t+1$, if new round $t$ messages are reliably received
later.
Recall that lines 11-16 are performed atomically.

\vspace*{8pt}
\hrule

\vspace*{2pt}

\noindent {\bf Verified Averaging Algorithm:} Steps at node $i$

\vspace*{4pt}

\hrule

\vspace*{8pt}

\noindent {\bf Initialization:} All sets used below are initialized to $\emptyset$.\\

\noindent {\bf Round 0:}

\begin{itemize}
\item $\RBSend(x_i,i,0)$ 		\hfill 1

\item
{\bf Event handler for event $\RBRecv(x, j, 0)$ at node $i$} : \hfill 2\\

\centerline{\underline{\em Lines 3-7 are performed atomically.}}

\begin{list}{}{}
\item[\hspace*{0.4in}$-$] $\R_i[0] := \R_i[0] \cup \{(x,j,-1)\}$ \hfill 3 \\

\indent {\em Comment}: The third element of the 3-tuple added to $\R$ above is set\\
 \indent as $-1$ to facilitate consistent treatment in round 0 and rounds $t>0$. \\

\item[\hspace*{0.4in}$-$] When $|\R_i[0]| \geq n-f$ \& $(x_i,i,-1) \in \R_i[0]$  \underline{both true for the first time},	\hfill 4

	\hspace*{1in} $\R^c_i[0]:= \R_i[0]$ \hfill 5

	\hspace*{1in} $h_i[0] := H(\R^c_i[0],0)$ \hfill 6

	\hspace*{1in} Proceed to round 1  \hfill 7
\end{list}

\end{itemize}

\noindent {\bf Round $t ~ (t \geq 1)$:}

\begin{itemize}
\item $\RBSend((h_i[t-1], \R^c_i[t-1]), i, t)$ \hfill 8

\item
{\bf Event handler for event $\RBRecv((h,\VV), j, t)$ at node $i$} : \hfill 9

\begin{itemize}

\item Wait until $\VV \subseteq \R_i[t-1]$ \hfill 10 \\

\centerline{\underline{\em Lines 11-16 are performed atomically.}}

\item
If $|\VV|\geq n-f$ and $(*,j,t-2)\in \VV$ and $h = H(\VV,t-1)$  \hfill 11 \\
then $\R_i[t] := \R_i[t] \cup \{(h,j,t-1)\}$ \hfill 12 \\

\item
When $|\R_i[t]| \geq n-f$ \& $(h_i[t-1],i,t-1) \in \R_i[t]$ \underline{both true for first time} \hfill 13

	\hspace*{1.6in} $\R^c_i[t]:= \R_i[t]$ \hfill 14

	\hspace*{1.6in} $h_i[t] := H(\R^c_i[t],t)$ \hfill 15

	\hspace*{1.6in} Proceed to round $t+1$  \hfill 16

\end{itemize}
\end{itemize}

\hrule

\vspace*{8pt}

~

\begin{definition}
\label{d_verified}
A node $k$'s execution of round $r$, $r\geq 0$, is said to be
\underline{\bf verified by a fault-free node $i$} if,
eventually node $i$ reliably receives message
of the form $((h,\VV),k,r+1)$ from node $k$, and subsequently adds $(h,k,r)$ to $\R_i[r+1]$.
Note that node $k$ may possibly be faulty.
Node $k$'s execution of round $r$ is said to be \underline{\bf verified} if it
is verified by at least one fault-free node.
\end{definition}
We now introduce some more notations (which are also summarized in Appendix \ref{app_s_notations}):
\begin{itemize}
\item For a {\em given} execution of the proposed {\em Verified Averaging} algorithm, let $F$ denote the {\em actual} set of faulty nodes in the execution.
Let $|F|=\phi$. Thus, $0\leq \phi\leq f$.
\item For $r\geq 0$, let $F_v[r]$ denote the set of faulty nodes whose round $r$ execution is verified by at least one fault-free node, as per Definition \ref{d_verified}.
Note that $F_v[r]\subseteq F$.
\item Define $\overline{F_v}[r]=F-F_v[r]$, for $r\geq 0$.
\end{itemize}
For each faulty node $k\in F_v[r]$, by
	Definition \ref{d_verified}, there must exist
	a fault-free node $i$ that eventually reliably receives
	a message of the form $((h,\VV),k,$r+1$)$ from
	node $k$, and adds $(h,k,r)$ to $\R_i[r+1]$.
	Given these $h$ and $\VV$, for future reference, let us define
	\begin{eqnarray}
	 h_k[r]&=&h \label{e_faulty_h} \\
	 \R^c_k[r]&=&\VV \label{e_faulty_V}
	\end{eqnarray}
	Node $i$ verifies node $k$'s round $r$ execution after node $i$
	has entered its round $r+1$.
	Since round $r$ execution of faulty node $k$ above is verified by fault-free
	node $i$, due to the check performed by node $i$ at line 11, the equality below holds
	for $h_k[r]$ and $\R^c_k[r]$ defined in (\ref{e_faulty_h}) and (\ref{e_faulty_V}).
	\begin{eqnarray}
	h_k[r]& = & H(\R^c_k[r],r) \label{e_faulty_H}
	\end{eqnarray}
(The proof of Claim \ref{c_ver} in Appendix \ref{a_claims} elaborates on the
above equality.)
	While the algorithm requires each node $k$ to maintain
	variables $h_k[r]$ and $\R^c_k[r]$, we cannot assume correct behavior on
	the part of the faulty nodes. However, from the perspective of
	each fault-free node that verifies the round $r$ execution of faulty node $k\in F_v[r]$,
	node $k$ behaves ``as if'' these local variable take the values specified
	in (\ref{e_faulty_h}) and (\ref{e_faulty_V}) that satisfy (\ref{e_faulty_H}).
	Note that if the round $r$ execution of a faulty node $k$ is verified by more than
	one fault-free node, due to the {\em Global Uniqueness} of reliable broadcast,
	all these fault-free nodes must have reliably received identical round $r+1$ messages
	from node $k$.

Proofs of Lemmas \ref{l_progress}, \ref{lemma:J_in_H0} and \ref{lemma:always_notFv_not_verified}
below
are presented in Appendices \ref{a_lemma_progress}, \ref{a_lemma_J}, and
\ref{a_always}, respectively.
These lemmas are used to prove the correctness of the {\em Verified Averaging} algorithm.

\begin{lemma}
\label{l_progress}
If all the fault-free nodes progress to the start of round $t$, $t \geq 0$, then all the
fault-free nodes will eventually progress to the start of round $t+1$.
\end{lemma}

\begin{lemma}
\label{lemma:J_in_H0}
For each node $i\in V-\overline{F_v}[0]$, the polytope $h_i[0]$ is non-empty.
\end{lemma}

\begin{lemma}
\label{lemma:always_notFv_not_verified}
For $r\geq 0$,
if $b \in \overline{F_v}[r]$, then for all $\tau\geq r$,
\begin{itemize}
\item $b \in \overline{F_v}[\tau]$, and
\item for all $i\in V-\overline{F_v}[\tau+1]$, $(*,b,\tau)\not\in\R^c_i[\tau+1]$.
\end{itemize}
\end{lemma}

\section{Correctness}
\label{s_correctness}

We first
introduce some terminology and definitions related to matrices. Then, we develop a
{\em transition matrix} representation of the proposed algorithm, and use
that to prove its correctness.

\subsection{Matrix Preliminaries}

We use boldface upper case letters to denote matrices, rows of matrices, and their elements. For instance, $\bfA$ denotes a matrix, $\bfA_i$ denotes the $i$-th row of matrix $\bfA$, and $\bfA_{ij}$ denotes the element at the intersection of the $i$-th row and the $j$-th column of matrix $\bfA$.

\begin{definition}
\label{d_stochastic}
A vector is said to be stochastic if all its elements
are non-negative, and the elements add up to 1.
A matrix is said to be row stochastic if each row of the matrix is a
stochastic vector.
\end{definition}
For matrix products, we adopt the ``backward'' product convention below, where $a \leq b$,
\begin{equation}
\label{backward}
\Pi_{\tau=a}^b \bfA[\tau] = \bfA[b]\bfA[b-1]\cdots\bfA[a]
\end{equation}
For a row stochastic matrix $\bfA$,
 coefficients of ergodicity $\delta(\bfA)$ and $\lambda(\bfA)$ are defined as
follows \cite{Wolfowitz}:
\begin{eqnarray*}
\delta(\bfA) & = &   \max_j ~ \max_{i_1,i_2}~ \| \bfA_{i_1\,j}-\bfA_{i_2\,j} \| \label{e_zelta} \\
\lambda(\bfA) & = & 1 - \min_{i_1,i_2} \sum_j \min(\bfA_{i_1\,j} ~, \bfA_{i_2\,j}) \label{e_lambda}
\end{eqnarray*}
\begin{claim}
\label{claim_zelta}
For any $p$ square row stochastic matrices $\bfA(1),\bfA(2),\dots, \bfA(p)$,
\begin{eqnarray*}
\delta(\Pi_{\tau=1}^p \bfA(\tau)) ~\leq ~
 \Pi_{\tau=1}^p ~ \lambda(\bfA(\tau)).
\end{eqnarray*}
\end{claim}
Claim \ref{claim_zelta} is proved in \cite{Hajnal58}.
Claim \ref{c_lambda_bound} below follows directly from the definition of $\lambda(\cdotp)$.
\begin{claim}
\label{c_lambda_bound}
If there exists a constant $\gamma$, where $0<\gamma\leq 1$, such
that, for any pair of rows $i,j$ of matrix $\bfA$, there exists a column
$g$ (that may depend on $i,j$) such that,
$\min (\bfA_{ig},\bfA_{jg}) \geq \gamma$,
then
 $\lambda(\bfA)\leq 1-\gamma<1$.
\end{claim}

Let $\vectorv$ be a column vector with $n$ elements, such that the $i$-th element of vector $\vectorv$, namely $\vectorv_i$, is a convex polytope in the $d$-dimensional Euclidean space. Let $\bfA$ be a $n\times n$ row stochastic square matrix. Then multiplication of
matrix $\bfA$ and vector $\vectorv$ is performed by multiplying each row of $\bfA$
with column vector $\vectorv$ of polytopes. Formally,
\begin{equation}
\label{eq:multiplication}
\bfA \vectorv = [H_l(\vectorv^T; \bfA_1)
~~~~~H_l(\vectorv^T; \bfA_2)
~~~~~...~~~~~
H_l(\vectorv^T; \bfA_n)]^T
\end{equation}
where $^T$ denotes the transpose operation (thus, $\vectorv^T$ is the transpose of $\vectorv$).
$H_l$ is defined in Definition \ref{def:linear_hull}. Thus, the result of the multiplication
$\bfA\vectorv$ is a column vector consisting of $n$ convex polytopes.
Similarly, product of row vector $\bfA_i$ and above vector $\vectorv$ is
obtained as follows, and it is a polytope.
\begin{eqnarray}
\label{e_r_c}
\bfA_i\vectorv &=& H_l(\vectorv^T\,;~\bfA_i)
\end{eqnarray}

\subsection{Transition Matrix Representation of {\em Verified Averaging}}

\noindent
Let $\vectorv[t]$, $t\geq 0$, denote a column vector of length $|V|=n$.
In the remaining discussion,
we will refer to $\vectorv[t]$ be the state of the system at the end of round $t$.
In particular,
$\vectorv_i[t]$ for $i\in V$ is viewed as
the state of node $i$ at the end of round $t$.
We define $\vectorv[0]$ as follows:
\begin{itemize}
\item[(I1)] For each fault-free node $i\in V-F$, $\vectorv_i[0]:=h_i[0]$.
\item[(I2)] For each faulty node $k\in F_v[0]$,
	$\vectorv_k[0]:=h_k[0]$, where
	$h_k[0]$ is defined in (\ref{e_faulty_h}).
\item[(I3)] For each faulty node $k\in \overline{F_v}[0]$,
$\vectorv_k[0]$ is {\em arbitrarily} defined as the origin,
or the all-0 vector.
We will justify this arbitrary choice later.
\end{itemize}

We will show that the state evolution
can be represented in a matrix form as in (\ref{matrix:alg1}), for a suitably
chosen $n\times n$ matrix  $\matrixm[t]$.
$\matrixm[t]$ is said to be the {\em transition matrix} for round $t$.
\begin{equation}
\label{matrix:alg1}
\vectorv[t] = \matrixm[t]~\vectorv[t-1], ~~~~~ t\geq 1
\end{equation}
For all $t\geq 0$,
Theorem \ref{t_M} below proves that,
for each $i\in V-\overline{F_v}[t]$, $h_i[t]=\vectorv_i[t]$.

Given a particular execution of the algorithm, we construct the transition matrix
$\bfM[t]$
for round $t\geq 1$ using
the following procedure.

\vspace*{10pt}
\hrule

\vspace*{2pt}

\noindent {\bf Construction of the Transition Matrix for Round $t~(t \geq 1)$}

\vspace*{4pt}

\hrule

\vspace*{8pt}

\begin{itemize}
\item For each node $i \in V-\overline{F_v}[t]$, and each $k\in V$:

\begin{list}{}{}
\item{} If $(h_k[t-1],k,t-1) \in \R^c_i[t]$, then
\begin{equation}
\label{eq:matrix_i}
\matrixm_{ik}[t] = \frac{1}{|\R^c_i[t]|}
\end{equation}

\item{} Otherwise,
\begin{equation}
\label{eq:matrix_i-2}
\matrixm_{ik}[t] = 0
\end{equation}

\end{list}

{\em Comment:}
For a faulty node $i\in F_v[t]$, $h_i[t]$ and $\R^c_i[t]$
are defined in (\ref{e_faulty_h}) and (\ref{e_faulty_V}).\\

\item For each node $j \in \overline{F_v}[t]$,
and each $k \in V$,
\begin{eqnarray}
\matrixm_{jk}[t] &=& \frac{1}{n}
\label{e_fv}
\end{eqnarray}

\end{itemize}

\hrule

\vspace*{8pt}

~

\begin{theorem}
\label{t_M}
For $r\geq 0$,
with state evolution
specified as $\vectorv[r+1]=\bfM[r+1]\vectorv[r]$ using $\bfM[r+1]$ constructed above,
for all $i\in V-\overline{F_v}[r]$,
(i) $h_i[r]$ is non-empty, and
 (ii) $h_i[r]=\vectorv_i[r]$.
\end{theorem}
\begin{proof}

The proof of the theorem is by induction. The theorem holds for $r=0$ due to
Lemma \ref{lemma:J_in_H0}, and the choice
of the elements of $\vectorv[0]$, as specified in (I1), (I2) and (I3) above.

Now, suppose that the theorem holds for $r=t-1$ where $t-1\geq 0$,
and prove it for $r=t$.
Thus, by induction hypothesis, for
all $i\in V-\overline{F_v}[t-1]$, $h_i[t-1]=\vectorv_i[t-1]\neq \emptyset$.
Now, $\vectorv[t] = \bfM[t]\vectorv[t-1]$.
\begin{itemize}
\item In round $t\geq 1$, each fault-free node $i\in V-F$ computes its new state
$h_i[t]$
at line 15 using function $H(\R^c_i[t],t)$. The function $H(\R^c_i[t],t)$ for $t\geq 1$ then computes a
linear combination of $|\R^c_i[t]|$ convex hulls, with all the weights being equal to
$\frac{1}{|\R^c_i[t]|}$.
Also,
by Definition \ref{d_verified} and the definition of $\overline{F_v}[t-1]$, if $(h,j,t-1)\in \R^c_i[t]$, then $j\not\in \overline{F_v}[t-1]$ (i.e., $j\in V-\overline{F_v}[t-1]$).
Therefore, if $(h,j,t-1)\in \R^c_i[t]$, then either $j$ is fault-free, or it is faulty and its
round $t-1$ execution is verified: thus, $h=h_j[t-1]$.

Also, by induction hypothesis, $h=h_j[t-1]\neq\emptyset$.
This implies that $h_i[t]=H(\R^c_i[t],t)$ is non-empty.

Then observe that, by defining $\bfM_{ik}[t]$ elements as in (\ref{eq:matrix_i})
and (\ref{eq:matrix_i-2}), we ensure that
$\bfM_i[t]\vectorv[t-1]$ equals $H(\R^c_i[t],t)$, and hence equals $h_i[t]$.

\item
For $i\in F_v[t]$ as well, as shown in (\ref{e_faulty_H}), $h_i[t]=H(\R^c_i[t],t)$, where
$h_i[t]$ and $\R^c_i[t]$ are as defined in (\ref{e_faulty_h})
and (\ref{e_faulty_V}).
The function $H(\R^c_i[t],t)$ for $t\geq 1$ then computes a
linear combination of $|\R^c_i[t]|$ convex hulls, with all the weights being equal to
$\frac{1}{|\R^c_i[t]|}$.
Consider an element $(h,j,t-1)$ in $\R^c_i[t]$. We argue that $j \in V - \overline{F_v}[t-1]$. Suppose this is not true, i.e., $j \in \overline{F_v}[t-1]$.
By Definition \ref{d_verified}, node $i$'s round $t$ execution is verified by some fault-free node $k$, which implies that eventually, $\R^c_i[t] \subseteq \R_k[t]$. However, since $k$ is fault-free, and $(h,j,t-1) \not\in \R_k[t]$, a contradiction.
Hence, if $(h,j,t-1)\in \R^c_i[t]$, then $j\in V-\overline{F_v}[t-1]$.
That is, if $(h,j,t-1)\in \R^c_i[t]$, then either $j$ is fault-free, or it is faulty and its
round $t-1$ execution is verified: thus, $h=h_j[t-1]$.

Also, by induction hypothesis, $h=h_j[t-1]\neq\emptyset$.
This implies that $h_i[t]=H(\R^c_i[t],t)$ is non-empty.

Then observe that, by defining $\bfM_{ik}[t]$ elements as in (\ref{eq:matrix_i})
and (\ref{eq:matrix_i-2}), we ensure that
$\bfM_i[t]\vectorv[t-1]$ equals $H(\R^c_i[t],t)$, and hence equals $h_i[t]$.

\end{itemize}

\comment{+++++++++ old incorrect proof+++++++++
Now, suppose that the theorem holds for $r=t-1$ where $t-1\geq 0$,
and prove it for $r=t$.
Thus, by induction hypothesis, for
all $i\in V-\overline{F_v}[t-1]$, $h_i[t-1]=\vectorv_i[t-1]\neq \emptyset$.
Now, $\vectorv[t] = \bfM[t]\vectorv[t-1]$.
\begin{itemize}
\item In round $t\geq 1$, each fault-free node $i\in V-F$ computes its new state
$h_i[t]$
at line 15 using function $H(\R^c_i[t],t)$. The function $H(\R^c_i[t],t)$ then computes a
linear combination of $|\R^c_i[t]|$ convex hulls, with all the weights being equal to
$|\R^c_i[t]|$.
Also, by Lemma \ref{lemma:always_notFv_not_verified},
if $(h,j,t-1)\in \R^c_i[t]$, then $j\not\in \overline{F_v}[t]$ (i.e., $j\in V-\overline{F_v}[t]$).
Therefore, if $(h,j,t-1)\in \R^c_i[t]$, then $h=h_j[t-1]$.
Also, by Lemma \ref{lemma:always_notFv_not_verified}, since $j\not\in \overline{F_v}[t]$,
we have $j\not\in\overline{F_v}[t-1]$.
Therefore, by induction hypothesis, $h=h_j[t-1]\neq\emptyset$.
This implies that $h_i[t]=H(\R^c_i[t],t)$ is non-empty.

Then observe that, by defining $\bfM_{ik}[t]$ elements as in (\ref{eq:matrix_i})
and (\ref{eq:matrix_i-2}), we ensure that
$\bfM_i[t]\vectorv[t-1]$ equals $H(\R^c_i[t],t)$, and hence equals $h_i[t]$.

\item
For $i\in F_v[t]$ as well, as shown (\ref{e_faulty_H}),
$h_i[t]=H(\R^c_i[t])$, where
$h_i[t]$ and $\R^c_i[t]$ are as defined in (\ref{e_faulty_h})
and (\ref{e_faulty_V}).
The rest of the argument in this case is similar to that
above for $i\in V-F$, and $\bfM_i[t]\vectorv[t-1]=h_i[t]$ for $i\in F_v[t]$.
\end{itemize}
+++++++++++++}
\end{proof}

\comment{++++++++ old text++++++
From the discussion in the above proof,
and Lemma \ref{lemma:always_notFv_not_verified}, it follows that
the state $\vectorv_j[t]$ for each node $j\in\overline{F_v}[t]$ does not affect the state
of the nodes $V-\overline{F_v}[\tau]$, for $\tau\geq t+1$. Therefore,
the value of $\vectorv_j[t]$ for $j\in \overline{F_v}[t]$ is
of no consequence to the other nodes.
++++++++++}
Now, we argue that for $t \geq 0$, the state $\vectorv_j[t]$ for each node $j\in\overline{F_v}[t]$ does not affect the state
of the nodes $V-\overline{F_v}[\tau]$, for $\tau\geq t+1$.
From the discussion in the above proof, we see that for $j\in\overline{F_v}[t]$, $(*,j,t)\not\in \R^c_i[t+1]$ for $i \in V-\overline{F_v}[t+1]$. Thus, the sate $\vectorv_j[t]$ does not affect the state $h_i[t+1]$. Then, by Lemma \ref{lemma:always_notFv_not_verified}, if $j\in\overline{F_v}[t]$, then $j\in\overline{F_v}[\tau]$, for $\tau\geq t+1$. Thus, by the same argument, the sate $\vectorv_j[\tau]$ does not affect the state $h_i[t+1]$.
This justifies the somewhat arbitrary choice of
$\vectorv_j[0]$ for $j\in\overline{F_v}[0]$,
and $\matrixm_{jk}[t]$ in (\ref{e_fv}) for $j\in \overline{F_v}[t], ~t\geq 1$.
This choice does simplify the remaining proof somewhat.

The above discussion shows that, for $t\geq 1$, the evolution of $\vectorv[t]$ can be written as
in (\ref{matrix:alg1}), that is, $\vectorv[t]=\matrixm[t]\vectorv[t-1]$.
Given the matrix product definition in (\ref{eq:multiplication}),
it is easy to verify that
\[ \bfM[\tau+1] ~ \left(\bfM[\tau] \vectorv[\tau-1]\right) ~=~ \left(\bfM[\tau+1]\bfM[\tau]\right)~\vectorv[\tau-1]
\text{~ for~} \tau \geq 1.
\]
Therefore, by repeated application of (\ref{matrix:alg1}),
we obtain:
\begin{eqnarray}
\vectorv[t] & = & \left(\,\Pi_{\tau=1}^t \matrixm[\tau]\,\right)\, \vectorv[0],
 ~~~~ t\geq 1
\label{e_unroll}
\end{eqnarray}
Recall that we adopt the ``backward'' matrix product convention presented in
(\ref{backward}).

\begin{lemma}
\label{lemma:transition_matrix}
For $t\geq 1$,
transition matrix $\bfM[t]$ constructed using the above procedure satisfies the following conditions.
\begin{itemize}
\item For $i,j \in V$, there exists a fault-free node $g(i,j)$ such that $\bfM_{ig(i,j)}[t] \geq \frac{1}{n}$.

\item $\bfM[t]$ is a row stochastic matrix, and $\lambda(\bfM[t])\leq 1-\frac{1}{n}$.

\end{itemize}

\end{lemma}
The proof of Lemma \ref{lemma:transition_matrix} is presented in Appendix \ref{a_t_matrix}.

\subsection{Correctness of {\em Verified Averaging}}

\begin{definition}
\label{def:valid_hull}
A convex polytope $h$ is said to be {\em valid} if every point in $h$ is in the convex hull of the inputs at the fault-free nodes.
\end{definition}

Lemmas \ref{lemma:valid_initial_hull}
and \ref{lemma:linear_valid} below are proved in
Appendices \ref{app_s_lemma:valid_initial_hull} and \ref{app_s_lemma:linear_valid}, respectively.

\begin{lemma}
\label{lemma:valid_initial_hull}
$h_i[0]$ for each node $i \in \sv - \overline{F_v}[0]$ is valid.
\end{lemma}

\begin{lemma}
\label{lemma:linear_valid}
Suppose non-empty convex polytopes $h_1, h_2, \cdots, h_k$ are all valid. Consider $k$ constants $c_1, c_2, \cdots, c_k$ such that $0 \leq c_i \leq 1$ and $\sum_{i = 1}^k c_i = 1$.
Then the linear combination of these convex polytopes,
$H_l(h_1, h_2, \cdots, h_k; c_1, c_2, \cdots, c_k)$, is valid.
\end{lemma}

\begin{theorem}
\label{thm:correctness}
Verified Averaging satisfies the {\em validity} and
{\em approximate agreement} properties
after a large enough number of asynchronous rounds.
\end{theorem}

\begin{proof}
Repeated applications of
Lemma \ref{l_progress} ensures that the fault-free nodes will progress
from round 0 through round $r$, for any $r\geq 0$, allowing us to use (\ref{e_unroll}). Consider round $t \geq 1$.
Let
\begin{eqnarray}\bfM^* &=& \Pi_{\tau=1}^t \bfM[\tau].
\label{e_Mstar}
\end{eqnarray}
(To simplify the presentation, we do not include the round index $[t]$
in the notation $\bfM^*$ above.)
 Then $\vectorv[t]=\bfM^* \vectorv[0]$.
By Lemma \ref{lemma:transition_matrix},
each $\bfM[t]$ is a {\em row stochastic} matrix,
therefore, $\bfM^*$ is also row stochastic.
By Lemma \ref{lemma:valid_initial_hull},
$h_i[0]=\vectorv_i[0]$ for each $i\in \sv - \overline{F_v}[0]$ is valid. Therefore,
by Lemma \ref{lemma:linear_valid}, $\bfM_i^*\vectorv[0]$ for each $i\in V-F$
is valid. Also, by Theorem \ref{t_M} and (\ref{e_unroll}), $h_i[t]=\bfM_i^*\vectorv[0]$
for $i\in V-F$.
Thus, $h_i[t]$ is valid, and {\em Verified Averaging} satisfies the
validity condition for all $t\geq 0$.

Let us define $\alpha = 1-\frac{1}{n}$.
By Lemma \ref{lemma:transition_matrix}, $\lambda(\bfM[t]) \leq 1-\frac{1}{n}=\alpha$. Then by Claim \ref{claim_zelta},
\begin{eqnarray}
\label{e_alpha}
\delta(\bfM^*)=
\delta(\Pi_{\tau=1}^t \bfM[\tau])
~ \leq ~ \lim_{t\rightarrow\infty} \Pi_{\tau=1}^{t} \lambda(\bfM[\tau])
 ~\leq~ {\left(1-\frac{1}{n}\right)}^t ~=~\alpha^t
\end{eqnarray}

Consider any two fault-free nodes $i,j\in V-F$.
By (\ref{e_alpha}), $\delta(\bfM^*)\leq \alpha^t$.
Therefore, by the definition of $\delta(\cdot)$, for $1 \leq k \leq n$,
\begin{equation}
\label{eq:delta1}
\| \bfM^*_{ik} -  \bfM^*_{jk}\| \leq \alpha^t
\end{equation}
By Lemma \ref{lemma:always_notFv_not_verified}, and construction of the transition matrices,
it should be easy to see that $\bfM^*_{ib}=0$ for $b\in \overline{F_v}[0]$.
Then, for
any point $p_i^*$ in $h_i[t]=\bfM^*_i\vectorv[0]$,
there must exist, for all $k\in V-\overline{F_v}[0]$, $p_k \in h_k[0],$ such that
\begin{equation}
\label{pi}
p_i^* = \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{ik} p_k
= \left(\sum_{k\in V-\overline{F_v}[0]} \bfM^*_{ik} p_k(1),~~\sum_{k\in V-\overline{F_v}[0]} \bfM^*_{ik} p_k(2), \cdots, \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{ik} p_k(d)\right)
\end{equation}
where $p_k(l)$ denotes the value of $p_k$'s $l$-th coordinate.
Now choose point $p_j^*$ in $h_j[t]$ defined as follows.
\begin{equation}
\label{pj}
p_j^* = \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{jk} p_k
= \left(\sum_{k\in V-\overline{F_v}[0]} \bfM^*_{jk} p_k(1),~~ \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{jk} p_k(2), \cdots,  \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{jk} p_k(d)\right)
\end{equation}

Then the Euclidean distance between $p_i^*$ and $p_j^*$ is $d(p_i^*,p_j^*)$.
The following derivation is obtained by simple algebraic manipulation, using (\ref{eq:delta1}),
(\ref{pi}) and (\ref{pj}).
The omitted steps in the algebraic manipulation are shown in Appendix \ref{a_hausdorff}.
\begin{align}
d(p_i^*, p_j^*) &=  \sqrt{\sum_{l=1}^d (p_i^*(l) - p_j^*(l))^2} = \sqrt{\sum_{l=1}^d \left(\sum_{k\in V-\overline{F_v}[0]} \bfM^*_{ik} p_k(l) - \sum_{k\in V-\overline{F_v}[0]} \bfM^*_{jk} p_k(l)\right)^2}\nonumber\\
&\leq \alpha^t \sqrt{\sum_{l=1}^d \left(\sum_{k\in V-\overline{F_v}[0]} \| p_k(l) \|\right)^2}\label{eq:d_pipj1} \leq ~~ \alpha^t\Omega
\end{align}
where
$\Omega = \max_{p_k \in h_k[0],k\in V-\overline{F_v}[0]} \sqrt{\sum_{l=1}^d (\sum_{k\in V-\overline{F_v}[0]} \| p_k(l) \|)^2}$.
Because the $h_k[0]$'s in the definition of $\Omega$ are all valid (by
Lemma \ref{lemma:valid_initial_hull}),
$\Omega$ can itself be upper bounded by a function of the input vectors at the fault-free
nodes.
Since $\alpha=1-\frac{1}{n}<1$, for large enough $t$, $\alpha^t\Omega<\epsilon$, for any given
$\epsilon$. Therefore, for fault-free $i,j$, for large enough $t$, for each point $p_i^*\in h_i[t]$
there exists a point $p_j^*[t]\in h_j[t]$ such that $d(p^*_i,p^*_j)<\epsilon$ (and, similarly, vice-versa).
Thus, by Definition \ref{def:dist}, eventually
Hausdorff distance $\D(h_i[t],h_j[t]) <\epsilon$.
Since this holds true for any pair of fault-free nodes
$i,j$, approximate agreement property is eventually satisfied.
\end{proof}

\section{Convex Polytope Obtained by {\em Verified Averaging}}
\label{s_size}

Recall that $|F|=\phi\leq f$.
Let $G = \cup_{i \in \{\sv - F\}} x_i$ be the set of the inputs at all fault-free nodes.
Thus, $|G|=n-\phi\geq n-f$.
Define a convex polytope $I$ as follows.
\begin{equation}
\label{eq:I}
I = \cap_{D \subset G, |D| = n - 2f-\phi}\, \HH(D)
\end{equation}

\begin{lemma}
\label{lemma:opt}
For all $i\in V-F$ and $t\geq 0$,
 $I\subseteq h_i[t]$.
\end{lemma}

The lemma is proved in Appendix \ref{a_opt}.
The lemma establishes a ``lower bound'' on the convex polytope that the fault-free
nodes decide on.
Due to Theorem \ref{t_M}, $h_i[t]$ is always non-empty for $i\in V-F$.
However, $I$ may possibly be empty, depending on the inputs at the fault-free
nodes, and the total number of nodes.
We believe that it may be possible to improve the above ``lower bound'' by using
a somewhat more complex algorithm (using the {\em stable-vectors} primitive
from \cite{herlihy_colorless_async} in round 0).
This improvement is left as a topic for future work. In the follow-up work \cite{Tseng_BCC_optimal}, we present an optimal algorithm that agrees on a convex polytope that is as {\em large} as possible under adversarial conditions

\section{Summary}

This paper addresses {\em Byzantine Convex Consensus} (BCC), wheresin each
node has a $d$-dimensional vector as its {\em input}, and each fault-free node has decide on a
polytope that is in the {\em convex hull} of the input vectors at the fault-free nodes.
We present an {\em asynchronous approximate} BCC algorithm with optimal fault tolerance,
and establish a lower bound on the convex polytope agreed upon.