text:
\begin{document}

\title{On Field Size and Success Probability in Network Coding}

\author{Olav Geil
\thanks{Department of Mathematical Sciences,
Aalborg University, Denmark, \protect\url{olav@math.aau.dk}},

Ryutaroh Matsumoto
\thanks{Department of Communications and Integrated Systems,
Tokyo Institute of Technology, Japan
\protect\url{ryutaroh@rmatsumoto.org}},

Casper Thomsen
\thanks{Department of Mathematical Sciences,
Aalborg University, Denmark,
\protect\url{caspert@math.aau.dk}}}

\maketitle

\begin{abstract}
    Using tools from algebraic geometry and Gr\"obner basis theory we solve two
    problems in network coding. First we present a method to determine the
    smallest field size for which linear network coding is feasible. Second we
    derive improved estimates on the success probability of random linear
    network coding. These estimates take into account which monomials occur in
    the support of the determinant of the product of Edmonds matrices. Therefore
    we finally investigate which monomials can occur in the determinant of the
    Edmonds matrix.\\

    \textbf{Keywords.} Distributed networking, linear network coding, multicast,
    network coding, random network coding.
\end{abstract}

\section{Introduction}\label{secintro}

In a traditional data network, an intermediate node only forwards data and never
modifies them.  Ahlswede et~al.~\cite{ahlswede00} showed that if we allow
intermediate nodes to process their incoming data and output modified versions
of them then maximum throughput can increase, and they also showed that the
maximum throughput is given by the minimum of maxflows between the source node
and a sink node for single source multicast on an acyclic directional network.
Such processing is called network coding.  Li et~al.~\cite{li03} showed that
computation of linear combinations over a finite field by intermediate nodes is
enough for achieving the maximum throughput.  Network coding only involving
linear combinations is called linear network coding.  The acyclic assumption
was later removed by Koetter and M\'edard \cite{koetter03}.

In this paper we shall concentrate on the error-free, delay-free multisource
multicast network connection problem where the sources are uncorrelated.
However, the proposed methods described can be generalized to deal with delays
as in \cite{ho06}. The only exception is the description in
Section~\ref{secpath}.

Considering multicast, it is important to decide whether or not all receivers
(called sinks) can recover all the transmitted information from the senders
(called sources).  It is also important to decide the minimum size $q$ of the
finite field $\F_q$ required for linear network coding.

Before using linear network coding we have to decide coefficients in linear
combinations computed by intermediate nodes. When the size $q$ of a finite field
is large, it is shown that random choice of coefficients allows all sinks to
recover the original transmitted information with high probability \cite{ho06}.
Such a method is called random linear network coding and the probability is
called success probability. As to random linear network coding the estimation or
determination of the success probability is very important.  Ho
et~al.~\cite{ho06} gave a lower bound on the success probability.

In their paper \cite{koetter03}, Koetter and M\'edard introduced an algebraic
geometric point view on network coding.  As explained in \cite{bn:cox},
computational problems in algebraic geometry can often be solved by Gr\"obner
bases. In this paper, we shall show that the exact computation of the minimum
$q$ can be made by applying the division algorithm for multivariate polynomials,
and we will show that improved estimates for the success probability can be
found by applying the footprint bound from Gr\"obner basis theory. These results
introduce a new approach to network coding study. As the improved estimates take
into account which monomials occur in the support of the determinant of a
certain matrix \cite{ho06} we study this matrix in details at the end of the
paper.

\section{Preliminary}\label{secpre}

We can determine whether or not all sinks can recover all the transmitted
information by the determinant of some matrix \cite{ho06}.  We shall review the
definition of such determinant.  Let $G=(V,E)$ be an directed acyclic graph with
possible parallel edges that represents the network topology.  The set of source
and sink nodes is denoted by $S$ and $T$ respectively. Assume that the source
nodes $S$ together get $h$ symbols in $\F_q$ per unit time and try to send them.

Identify the edges in $E$ with the integers $1, \ldots, |E|$. For an edge
$j=(u,v)$ we write $\text{head}(j)=v$ and $\text{tail}(j)=u$. We define the $|E|
\times |E|$ matrix $F=( f_{i,j})$ where  $f_{i,j}$ is a variable if
$\text{head}(i)=\text{tail}(j)$ and $f_{i,j}=0$ otherwise. The variable
$f_{i,j}$ is the coding coefficient from $i$ to $j$.

Index $h$ symbols in $\F_q$ sent by $S$ by $1, \ldots, h$.  We also define an $h
\times |E|$ matrix $A=(a_{i,j})$ where $a_{i,j}$ is a variable if the edge $j$
is an outgoing edge from the source $s \in S$ sending the $i$-th symbol and
$a_{i,j}=0$ otherwise. Variables $a_{i,j}$ represent how the source nodes send
information to their outgoing edges.

Let $X(l)$ denote the $l$-th symbol generated by the sources $S$, and let $Y(j)$
denote the information sent along edge $j$.  The model is described by the
following relation
\begin{equation*}
    Y(j)
    = \sum_{i=1}^h a_{i,j} X(i)+\sum_{i: \text{head}(i)=\text{tail}(j)}
        f_{i,j}Y(i)
    .
\end{equation*}

For each sink $t \in T$ define an $h \times |E|$ matrix $B_t$ whose $(i,j)$
entry $b_{t,i,j}$ is a variable if $\text{head}(j)=t$ and equals $0$ otherwise.
The index $i$ refers to the $i$-th symbol sent by one of the sources.  Thereby
variables $b_{t,i,j}$ represent how the sink $t$ process the received data from
its incoming edges.

The sink $t$ records the vector
\begin{equation*}
    \vec{b^{(t)}} = \big( b_1^{(t)}, \ldots , b_h^{(t)} \big)
\end{equation*}
where
\begin{equation*}
    b_i^{(t)} = \sum_{j: \text{head}(j)=t} b_{t,i,j}Y(j)
    .
\end{equation*}
We now recall from~\cite{ho06} under which conditions all informations sent by
the sources can always be recovered at all sinks. As in \cite{ho06} we define
the Edmonds matrix $M_t$ for $t\in T$ by
\begin{equation}
    \label{eqMt}
    M_t =
    \begin{pmatrix}
        A & 0 \\
        I-F & B_t^T
    \end{pmatrix}
    .
\end{equation}
Define the polynomial $P$ by
\begin{equation}
    \label{eqP}
    P = \prod_{t\in T} |M_t|
    .
\end{equation}
$P$ is a multivariate polynomial in variables $f_{i,j}$, $a_{i,j}$ and
$b_{t,i,j}$. Assigning a value in $\F_q$ to each variable corresponds to
choosing a coding scheme. Plugging the assigned values into $P$ gives an element
$k \in \F_q$. The following theorem from~\cite{ho06} tells us when the coding
scheme can be used to always recover the information generated at the sources
$S$ at all sinks in $T$.

\begin{theorem}
    Let the notation and the network coding model be as above. Assume a coding
    scheme has been chosen by assigning values to the variables $f_{i,j}$,
    $a_{i,j}$ and $b_{t,i,j}$. Let $k$ be the value found by plugging the
    assigned values into $P$. Every sink $t \in T$ can recover from
    $\vec{b^{(t)}}$ the informations $X(1), \ldots , X(h)$ no matter what they
    are, if and only if $k \neq 0$ holds.
\end{theorem}

\begin{proof}
    See~\cite{ho06}.
\end{proof}

\section{Computation of the Minimum Field Size}\label{secfieldsize}

We shall study computation of the minimum symbol size $q$. For this purpose we
will need the division algorithm for multivariate polynomials
\cite[Sec.~2.3]{bn:cox} to produce the remainder of a polynomial $F(X_1, \ldots
, X_n)$ modulo $(X_1^q-X_1, \ldots, X_n^q-X_n)$ (this remainder is independent
of the choice of monomial ordering). We adapt the standard notation for the
above remainder which is
\begin{equation*}
    F(X_1, \ldots, X_n) \rem (X_1^q-X_1, \ldots, X_n^q-X_n)
    .
\end{equation*}
The reader unfamiliar with the division algorithm can think of the above
remainder of $F(X_1, \ldots, X_n)$ as the polynomial produced by the following
procedure. As long as we can find an $X_i$ such that $X_i^q$ divides some term
in the polynomial under consideration we replace the factor $X_i^q$ with $X_i$
wherever it occurs. The process continues until the $X_i$-degree is less than
$q$ for all $i=1, \ldots, n$. It is clear that the above procedure can be
efficiently implemented.

\begin{proposition}\label{pro1}
    Let $F(X_1$, \ldots, $X_n)$ be an $n$-variate polynomial over $\F_q$.  There
    exists an $n$-tuple $(x_1, \ldots, x_n) \in \F_q^n$ such that $F(x_1,
    \ldots, x_n) \neq 0$ if and only if
    \begin{equation*}
        F(X_1, \ldots, X_n) \rem (X_1^q-X_1, \ldots, X_n^q-X_n)
        \neq 0
        .
    \end{equation*}
\end{proposition}

\begin{proof}
    As $a^q=a$ for all $a \in \F_q$ it holds that  $F(X_1, \ldots, X_n)$
    evaluates to the same as $R(X_1, \ldots, X_n):=F(X_1, \ldots, X_n) \rem
    (X_1^q-X_1, \ldots, X_n^q-X_n)$ in every $(x_1, \ldots, x_n) \in \F_q^n$.
    If  $R(X_1, \ldots, X_n)=0$ therefore $F(X_1, \ldots, X_n)$ evaluates to
    zero for every choice of $(x_1, \ldots, x_n) \in \F_q^n$. If $R(X_1, \ldots,
    X_n)$ is nonzero we consider it first as a polynomial in $\F_q(X_1, \ldots,
    X_{n-1})[X_n]$ (that is, a polynomial in one variable over the quotient
    field $\F_q(X_1, \ldots, X_{n-1})$). But the $X_n$-degree is at most $q-1$
    and therefore it has at most $q-1$ zeros. We conclude that there exists an
    $x_n\in \F_q$ such that $R(X_1, \ldots, X_{n-1},x_n) \in \F_q[X_1, \ldots,
    X_{n-1}]$ is nonzero. Continuing this way we find $(x_1, \ldots, x_n)$ such
    that $R(x_1, \ldots, x_n)$ and therefore also $F(x_1, \ldots, x_n)$ is
    nonzero.
\end{proof}

From~\cite[Th.\ 2]{ho06} we know that for all prime powers $q$ greater than
$|T|$ linear network coding is possible. It is now straightforward to describe
an algorithm that finds the smallest field $\F_q$ of prescribed characteristic
$p$ for which linear network coding is feasible. We first reduce the polynomial
$P$ from~\eqref{eqP} modulo the prime $p$.  We observe that although $P$ is a
polynomial in all the variables $a_{i,j}$, $b_{t,i,j}$, $f_{i,j}$ the variable
$b_{t,i,j}$ appears at most in powers of $1$. This is so as it appears at most
in a single entry in $M_t$ and does not appear elsewhere.  Therefore $\F_q$ can
be used for network coding if $P \rem p$ does not reduce to zero modulo the
polynomials $a_{i,j}^q-a_{i,j}$, $f_{i,j}^q-f_{i,j}$. To decide the smallest
field $\F_q$ of characteristic $p$ for which network coding is feasible we try
first $\F_q=\F_p$. If this does not work we then try $\F_{p^2}$ and so on. To
find an $\F_q$ that works we need at most to try $\lfloor \log_p(|T|)\rfloor$
different fields as we know that linear network coding is possible whenever $q >
|T|$.

Note that once a field $\F_q$ is found such that the network connection problem
is feasible the last part of the proof of Proposition 1 describes a simple way
of deciding coefficients $(x_1, \ldots, x_n) \in \F_q^n$ that can be used for
network coding.

From~\cite[Sec.~7.1.3]{frasol} we know that it is an NP-hard problem to find the
minimum field size for linear network coding. Our findings imply that it is
NP-hard to find the polynomial $P$ in~\eqref{eqP}.

\section{Computation of the Success Probability of Random Linear Network
Coding}\label{secrandom}

In random linear network coding we from the beginning fix for a collection
\begin{equation*}
    K \subseteq \{1, \ldots, h\} \times \{1, \ldots, | E | \}
\end{equation*}
the $a_{i,j}$'s with $(i,j) \in K$ and also we fix  for a collection
\begin{equation*}
    J \subseteq \{1, \ldots, | E| \} \times \{1, \ldots,|E| \}
\end{equation*}
the $f_{i,j}$'s with $(i,j) \in J$. This is done in a way such that there exists
a solution to the network connection problem with the same values for these
fixed coefficients. A priori of course we let $a_{i,j}=0$ if the edge $j$ is not
emerging from the source sending information $i$, and also a priori we of course
let $f_{i,j}=0$ if $j$ is not an adjacent downstream edge of $i$.  Besides these
a priori fixed values there may be good reasons for also fixing other
coefficients $a_{i,j}$ and $f_{i,j}$ \cite{ho06}. If for example there is only
one upstream edge $i$ adjacent to $j$ we may assume $f_{i,j}=1$.  All the
$a_{i,j}$'s and $f_{i,j}$'s which have not been fixed at this point are then
chosen randomly and independently. All coefficients are to be elements in
$\F_q$. If a solution to the network connection problem exists with the
$a_{i,j}$'s and the $f_{i,j}$'s specified,  it is possible to determine values
of $b_{t,i,j}$ at the sinks such that a solution to the network connection
problem is given. Let $\mu$ be the number of variables $a_{i,j}$ and $f_{i,j}$
chosen randomly. Call these variables $X_1, \ldots, X_\mu$. Consider the
polynomial $P$ in~\eqref{eqP} and let $\tilde{P}$ be the polynomial made from
$P$  by plugging in the fixed values of the $a_{i,j}$'s and the fixed values of
the $f_{i,j}$'s (calculations taking place in $\F_q$). Then $\tilde{P}$ is a
polynomial in $X_1, \ldots, X_\mu$. The coefficients of $\tilde{P}$ are
polynomials in the $b_{t,i,j}$'s over $\F_q$. Finally, define
\begin{equation*}
    \widehat{P} := \tilde{P} \rem (X_1^q-X_1, \ldots, X_\mu^q-X_\mu)
    .
\end{equation*}
The success probability of random linear network coding is the probability that
the random choice of coefficients will lead to a solution of the network
connection problem\footnote{This corresponds to saying that each sink can
recover the data at the maximum rate promised by network coding.} as in
Section~\ref{secpre}. That is, the probability is the number
\begin{align}
    \nonumber
    &\frac{| \{ (x_1, \ldots, x_\mu) \in \F_q^\mu | \tilde{P}(x_1, \ldots,
        x_\mu) \neq 0 \} |}{q^\mu}
    \\
    \label{eqprob2}
    &=\frac{| \{ (x_1, \ldots, x_\mu) \in \F_q^\mu | \widehat{P}(x_1, \ldots,
        x_\mu) \neq 0 \} |}{q^\mu}
    .
\end{align}

To see the first result observe that for fixed $(x_1, \ldots, x_\mu) \in
\F_q^\mu$, $\tilde{P}(x_1, \ldots, x_\mu)$ can be viewed as a polynomial in the
variables $b_{t,i,j}$'s with coefficients in $\F_q$ and recall that the
$b_{t,i,j}$'s occur in powers of at most $1$. Therefore, if $\tilde{P}(x_1,
\ldots, x_\mu) \neq 0$, then by Proposition~\ref{pro1} it is possible to choose
the $b_{t,i,j}$'s such that if we plug them into $\tilde{P}(x_1, \ldots, x_\mu)$
then we get nonzero.  The last result follows from the fact that $\tilde{P}(x_1,
\ldots, x_\mu) = \widehat{P}(x_1, \ldots, x_\mu)$ for all $(x_1, \ldots, x_\mu)
\in \F_q^\mu$. In this section we shall present a method to estimate the success
probability using Gr\"obner basis theoretical methods.

We briefly review some basic definitions and results of Gr\"obner bases.  See
\cite{bn:cox} for a more detailed exposition.  Let $\mathcal{M}(X_1, \ldots,
X_n)$ be the set of monomials in the variables $X_1, \ldots, X_n$.  A monomial
ordering $\prec$ is a total ordering on $\mathcal{M}(X_1, \ldots, X_n)$ such
that
\begin{equation*}
    L \prec M \Longrightarrow LN \prec MN
\end{equation*}
holds for all monomials $L$, $M$, $N \in \mathcal{M}(X_1, \ldots, X_n)$ and such
that every nonempty subset of $\mathcal{M}(X_1, \ldots, X_n)$ has a unique
smallest element with respect to $\prec$.  The leading monomial of a polynomial
$F$ with respect to $\prec$, denoted by $\textsc{lm}(F)$, is the largest
monomial in the support of $F$. Given a polynomial ideal $I$ and a monomial
ordering the footprint $\Delta_{\prec}(I)$ is the set of monomials that cannot
be found as leading monomials of any polynomial in $I$. The following
proposition explains our interest in the footprint (for a proof of the
proposition see~\cite[Pro.~8.32]{beckerweis}).

\begin{proposition}\label{profootgen}
    Let $\F$ be a field and consider the polynomials $F_1, \ldots, F_s \in
    \F[X_1, \ldots, X_n]$. Let $I=\langle F_1, \ldots, F_s \rangle \subseteq
    \F[X_1, \ldots, X_n]$ be the ideal generated by $F_1, \ldots, F_s$. If
    $\Delta_{\prec}(I)$ is finite then the number of common zeros of $F_1,
    \ldots, F_s$ in the algebraic closure of $\F$ is at most equal to $|
    \Delta_{\prec}(I)|$.
\end{proposition}

Proposition~\ref{profootgen} is known as the footprint bound. It
has the following corollary.

\begin{corollary}\label{corweakfoot}
    Let $F \in \F[X_1, \ldots,X_n]$ where $\F$ is a field containing $\F_q$.
    Fix a monomial ordering and let
    \begin{equation*}
        X_1^{j_1}\cdots X_n^{j_n}
        = \textsc{lm}\big( F \rem (X_1^q-X_1, \ldots, X_n^q-X_n) \big)
        .
    \end{equation*}
    The number of zeros of $F$ over $\F_q$ is at most equal to
    \begin{equation}
        \label{eqcirc1}
        q^n-\prod_{v=1}^n(q-j_v)
        .
    \end{equation}
\end{corollary}

\begin{proof}
    We have
    \begin{multline*}
        \Delta_{\prec}(\langle F, X_1^q-X_1, \ldots, X_n^q-X_n \rangle )
        \\
        \subseteq \Delta_{\prec}(\langle \textsc{lm}(F \rem (X_1^q-X_1,
        \ldots,X_n^q-X_n)),X_1^q, \ldots, X_n^q \rangle )
    \end{multline*}
    and the size of the latter set equals~\eqref{eqcirc1}. The result now
    follows immediately from Proposition~\ref{profootgen}.
\end{proof}

\begin{theorem}\label{theestimate}
    Let as above $\tilde{P}$ be found by plugging into $P$ some fixed values for
    the variables $a_{i,j}$, $(i,j) \in K$, and by plugging into $P$ some fixed
    values for the variables $f_{i,j}$, $(i,j) \in J$, and by leaving the
    remaining $\mu$ variables flexible. Assume as above that there exists a
    solution to the network connection problem with the same values for these
    fixed coefficients. Denote by $X_1, \ldots, X_\mu$ the variables to be
    chosen by random and define $\widehat{P}:= \tilde{P} \rem (X_1^q-X_1,
    \ldots, X_\mu^q-X_\mu)$. (Note that if $q > |T|$ then
    $\widehat{P}=\tilde{P}$). Consider $\widehat{P}$ as a polynomial in the
    variables $X_1, \ldots, X_\mu$  and let $\prec$ be any fixed monomial
    ordering. Writing $X_1^{j_1}\cdots X_\mu^{j_\mu} = \textsc{lm}(\widehat{P})$
    the success probability is at least
    \begin{equation}
        \label{eqtrio}
        q^{- \mu} \prod_{v=1}^\mu (q-j_v)
        .
    \end{equation}
    As a consequence the success probability is in particular at least
    \begin{equation}
        \label{eqtri}
        q^{- \mu} \min \bigg\{
        \prod_{i=1}^\mu (q-s_i)
        \bigg|
        X_1^{s_1}\cdots X_\mu^{s_\mu} \text{ is a monomial in the
        support of } \widehat{P}
        \bigg\}
        .
    \end{equation}
\end{theorem}

\begin{proof}
    Let $\F$ be the quotient field $\F_q(X_1, \ldots, X_\mu)$. The result
    in~\eqref{eqtrio} now follows by applying Corollary~\ref{corweakfoot}
    and~\eqref{eqprob2}. As the leading monomial of $\tilde{P}$ is of course a
    monomial in the support of $\tilde{P}$ \eqref{eqtri} is smaller or equal
    to~\eqref{eqtrio}.
\end{proof}

\begin{remark}
    The condition in Theorem~\ref{theestimate} that there exists a solution to
    the network connection problem with the coefficients corresponding to $K$
    and $J$ being as specified is equivalent to the condition that $\widehat{P}
    \neq 0$.
\end{remark}

We conclude this section by mentioning without a proof that Gr\"obner basis
theory tells us that the true success probability can be calculated as
\begin{equation*}
    q^{-\mu} \big( q^\mu-| \Delta_{\prec}(\langle \tilde{P}, X_1^q-X_1, \ldots,
    X_\mu^q-X_\mu \rangle )| \big)
    .
\end{equation*}
This observation is however of little value as it seems very difficult to
compute the footprint
\begin{equation*}
    \Delta_{\prec}(\langle \tilde{P}, X_1^q-X_1, \ldots, X_\mu^q-X_\mu \rangle)
\end{equation*}
due to the fact that $\mu$ is typically a very high number.

\section{The Bound by Ho et al.}\label{secho}

In~\cite{ho06} Ho et al.\ gave a lower bound on the success probability in terms
of the number of edges $j$ with associated random coefficients\footnote{We state
Ho et al.'s bound only in the case of delay-free acyclic networks.} $\{a_{i,j},
f_{l,j}\}$. Letting $\eta$ be the number of such edges \cite[Th.\ 2]{ho06} tells
us that if $q>|T|$ and if there exists a solution to the network connection
problem with the same values for the fixed coefficients, then the success
probability is at least
\begin{equation}
    \label{eqhosnabel1}
    p_\text{Ho} = \left( \frac{q-|T|}{q} \right)^\eta
    .
\end{equation}
The proof in~\cite{ho06} of~\eqref{eqhosnabel1} relies on two lemmas of which we
only state the first one.

\begin{lemma}\label{lemho3}
    Let $\eta$ be defined as above. The determinant polynomial of $M_t$ has
    maximum degree $\eta$ in the random variables $\{a_{i,j}, f_{l,j} \}$ and is
    linear in each of these variables.
\end{lemma}

\begin{proof}
    See~\cite[Lem.\ 3]{ho06}. Alternatively the proof can be derived as a
    consequence of Theorem~\ref{thepath} in Section~\ref{secpath}.
\end{proof}

Recall, that the polynomial $P$ in~\eqref{eqP} is the product of the
determinants $|M_t|$, $t\in T$. Lemma~\ref{lemho3} therefore implies that the
polynomial $\tilde{P}$ has at most total degree equal to $|T| \eta$ and that no
variable appears in powers of more than $|T|$. The assumption $q > |T|$ implies
$\widehat{P}=\tilde{P}$ which makes it particular easy to see that the same of
course holds for $\widehat{P}$.  Combining this observation with the following
lemma shows that the numbers in~\eqref{eqtrio} and \eqref{eqtri} are both at
least as large as the number~\eqref{eqhosnabel1}.

\begin{lemma}
    Let $\eta, |T|, q \in \mathbf{N}$, $|T| < q$ be some fixed numbers. Let
    $\mu, x_1, \ldots, x_\mu \in \mathbf{N}_0$ satisfy
    \begin{equation*}
        0 \leq x_1 \leq |T|, \ldots, 0 \leq x_\mu \leq |T|
    \end{equation*}
    and $x_1 + \cdots + x_\mu \leq |T|\eta$. The minimal value of
    \begin{equation*}
        \prod_{i=1}^\mu \left( \frac{q-x_i}{q} \right)
    \end{equation*}
    (taken over all possible values of $\mu, x_1, \ldots, x_\mu$) is
    \begin{equation*}
        \left( \frac{q-|T|}{q} \right)^\eta
        .
    \end{equation*}
\end{lemma}

\begin{proof}
    Assume $\mu$ and $x_1, \ldots, x_\mu$ are chosen such that the expression
    attains its minimal value. Without loss of generality we may assume that
    \begin{equation*}
        x_1 \geq x_2 \geq \cdots \geq x_\mu
    \end{equation*}
    holds. Clearly, $x_1+\cdots + x_\mu=|T|\eta$ must hold. If $x_i<|T|$ and
    $x_{i+1}>0$ then
    \begin{equation*}
        (q-x_i)(q-x_{i+1} )> (q-(x_i+1))(q-(x_{i+1}-1))
    \end{equation*}
    which cannot be the case. So $x_1=\cdots =x_{\eta}=|T|$. The remaining
    $x_j$'s if any all equal zero.
\end{proof}

\section{Examples}

In this section we apply the methods from the previous sections to two concrete
networks. We will see that the estimate on the success probability of random
linear network coding that was described in Theorem~\ref{theestimate} can be
considerably better than the estimate described in~\cite[Th.\ 2]{ho06}. Also we
will apply the method from Section~\ref{secfieldsize} to determine the smallest
field of characteristic two for which network coding can be successful.

As random linear network coding is assumed to take place at the nodes in a
decentralized manner, one natural choice is to set $f_{i,j}=1$ whenever the
indegree of the end node of edge $i$ is one and $j$ is the downstream edge
adjacent to $i$. Clearly, if $j$ is not a downstream edge adjacent to $i$ we set
$f_{i,j}=0$. Whenever none of the above is the case we may choose $f_{i,j}$
randomly. Also if there is only one source and the outdegree of the source is
equal to the number of symbols to be send we may enumerate the edges from the
source by the numbers $1, \ldots, h$ and set $a_{i,j}=1$ if $1 \leq i=j \leq h$
and set $a_{i,j}=0$ otherwise. This strategy can be generalized also to deal
with the case of more sources. In the following two examples we will choose the
variables in the manner just described. The network in the first example is
taken from~\cite[Ex.~3.1]{frasol} whereas the network in the second example is
new.

\begin{example}\label{exsidstenye}
    Consider the delay-free and acyclic network in Figure~\ref{figny}.
    \begin{figure}
        \begin{center}
            \scalebox{0.909}{
            $
            \psmatrix[colsep=0.8cm,rowsep=.4cm,mnode=circle]
            &v_1\\
            v_2&&v_3\\
            &v_4\\
            &v_5\\
            v_6&&v_7\\
            v_8&&v_9\\
            &v_{10}\\
            &v_{11}\\
            v_{12}&&v_{13}
            \psset{arrowscale=2}
            \ncline{->}{1,2}{2,1}^{1}
            \ncline{->}{1,2}{2,3}^{2}
            \ncline{->}{2,1}{3,2}^{3}
            \ncline{->}{2,1}{5,1}>{4}
            \ncline{->}{2,3}{3,2}^{5}
            \ncline{->}{2,3}{5,3}<{6}
            \ncline{->}{3,2}{4,2}<{7}
            \ncline{->}{4,2}{5,1}^{8}
            \ncline{->}{4,2}{5,3}^{9}
            \ncline{->}{5,1}{6,1}>{10}
            \ncline{->}{5,3}{6,3}<{11}
            \ncline{->}{6,1}{7,2}^{12}
            \ncline{->}{6,1}{9,1}>{13}
            \ncline{->}{6,3}{7,2}^{14}
            \ncline{->}{6,3}{9,3}<{15}
            \ncline{->}{7,2}{8,2}<{16}
            \ncline{->}{8,2}{9,1}^{17}
            \ncline{->}{8,2}{9,3}^{18}
            \endpsmatrix
            $
            }
        \end{center}
        \caption{The network from Example~\ref{exsidstenye}}
        \label{figny}
    \end{figure}
    There is one sender $v_1$ and two receivers $v_{12}$ and $v_{13}$. The
    min-cut max-flow number is two for both receivers so we assume that two
    independent random processes emerge from sender $v_1$.  We consider in this
    example only fields of characteristic $2$. Following the description
    preceding the example we set $a_{1,1}=a_{2,2}=1$ and $a_{i,j}=0$ in all
    other cases. Also we let $f_{i,j}=1$ except
    \begin{equation*}
        f_{3,7}, f_{5,7}, f_{4,10}, f_{8,10}, f_{9,11}, f_{6,11},
        f_{12,16}, f_{14,16}
    \end{equation*}
    which we choose by random. As in the previous sections we consider
    $b_{t,i,j}$ as fixed but unknown to us. The determinant polynomial becomes
    \begin{equation*}
        \tilde{P}=(b^2c^2e^2gh+c^2f^2gh+a^2d^2f^2gh)Q
        ,
    \end{equation*}
    where
    \begin{align*}
        a &=f_{3,7}  & b &=f_{5,7}  & c &=f_{4,10}  & d &=f_{8,10} \\
        e &=f_{9,11} & f &=f_{6,11} & g &=f_{12,16} & h &=f_{14,16}
    \end{align*}
    and $Q=|B_{v_{12}}'| \, |B_{v_{13}}'|$. Here, $B_{v_{12}}'$ respectively
    $B_{v_{13}}'$ is the matrix consisting of the nonzero columns of
    $B_{v_{12}}$ respectively the nonzero columns of $B_{v_{14}}$.  Restricting
    to fields $\F_q$ of size at least $4$ we have $\widehat{P}=\tilde{P}$ and we
    can therefore immediately apply the bounds in Theorem~\ref{theestimate}.
    Applying~\eqref{eqtri} we get the following lower bound on the success
    probability
    \begin{equation*}
        P_\text{new 2}(q) = \frac{(q-2)^3(q-1)^2}{q^5}
        .
    \end{equation*}
    Choosing as monomial ordering the lexicographic ordering $\plex$ with
    \begin{equation*}
        a \plex b \plex d \plex e \plex g \plex h \plex f \plex c
    \end{equation*}
    the leading monomial of $\tilde{P}$ becomes $c^2f^2gh$ and therefore
    from~\eqref{eqtrio} we get the following lower bound on the success
    probability
    \begin{equation*}
        P_\text{new 1}(q) = \frac{(q-2)^2(q-1)^2}{q^4}
        .
    \end{equation*}
    For comparison the bound~\eqref{eqhosnabel1} from~\cite{ho06} states that
    the success probability is at least
    \begin{equation*}
        P_{\text{Ho}}(q)=\frac{(q-2)^4}{q^4}
        .
    \end{equation*}
    We see that $P_\text{new 1}$ exceeds $P_\text{Ho}$ with a factor
    $(q-1)^2/(q-2)^2$, which is larger than 1. Also $P_\text{new 2}$ exceeds
    $P_\text{Ho}$. In Table~\ref{tabeins} we list values of $P_\text{new 1}(q)$,
    $P_\text{new 2}(q)$ and $P_\text{Ho}(q)$ for various choices of $q$.

    \begin{table}
        \centering
        \caption{From Example~\ref{exsidstenye}: Estimates on the success
        probability}
        \label{tabeins}
        \begin{tabular}{r ccccc}
            \toprule
            $q$ & 4 & 8 & 16 & 32 & 64
            \\
            $P_\text{new 1}(q)$ & 0.140 & 0.430 & 0.672 & 0.893 & 0.909
            \\
            $P_\text{new 2}(q)$ & $0.703 \times 10^{-1}$ & 0.322 & 0.588 &
            0.773 & 0.880
            \\
            $P_\text{Ho}(q)$  & $0.625 \times 10^{-1}$ & 0.316 & 0.586 & 0.772
            & 0.880
            \\
            \bottomrule
        \end{tabular}
    \end{table}
    We next consider the field $\F_2$. We reduce $\tilde{P}$ modulo $(a^2-a,
    \ldots, h^2-h)$ to get
    \begin{equation*}
        \widehat{P} = (bcegh + cfgh + adfgh) Q
        .
    \end{equation*}
    From~\eqref{eqtri} we see that the success probability of random network
    coding is at least $2^{-5}$. Choosing as monomial ordering the lexicographic
    ordering described above~\eqref{eqtrio} tells us that the success
    probability is at least $2^{-4}$. For comparison the
    bound~\eqref{eqhosnabel1} does not apply as we do not have $q>|T|$. It
    should be mentioned that for delay-free acyclic networks the network coding
    problem is solvable for all choices of $q\geq |T|$ \cite{jaggi03} and
    \cite{sanders03}. From this fact one can only conclude that the success
    probability is at least $2^{-8}$ (8 being the number of coefficients to be
    chosen by random).
\end{example}

\begin{example}\label{exone}
    Consider the network in Figure~\ref{figone}.
    \begin{figure}
        \begin{center}
            \scalebox{0.909}{
            $
            \psmatrix[colsep=0.8cm,rowsep=1cm,mnode=circle]
            &&v_1\\
            v_2&&v_3&&v_4\\
            &v_5&v_6&v_7\\
            &v_8&v_9&v_{10}\\
            v_{11}&&v_{12}&&v_{13}
            \psset{arrowscale=2}
            \ncline{->}{1,3}{2,1}<{1}
            \ncline{->}{1,3}{2,3}<{2}
            \ncline{->}{1,3}{2,5}>{3}
            \ncline{->}{2,1}{3,2}<{4}
            \ncline{->}{2,1}{3,3}^{5}
            \ncline{->}{2,1}{5,1}<{6}
            \ncline{->}{2,3}{3,2}^{7}
            \ncline{->}{2,3}{3,3}<{8}
            \ncline{->}{2,3}{3,4}^{9}
            \ncline{->}{2,5}{3,3}^{10}
            \ncline{->}{2,5}{3,4}>{11}
            \ncline{->}{2,5}{5,5}>{12}
            \ncline{->}{3,2}{4,2}<{13}
            \ncline{->}{3,3}{4,3}<{14}
            \ncline{->}{3,4}{4,4}>{15}
            \ncline{->}{4,2}{5,1}^{16}
            \ncline{->}{4,2}{5,3}>{17}
            \ncline{->}{4,3}{5,1}_{18}
            \ncline{->}{4,3}{5,3}>{19}
            \ncline{->}{4,3}{5,5}_{20}
            \ncline{->}{4,4}{5,3}_{21}
            \ncline{->}{4,4}{5,5}^{22}
            \endpsmatrix
            $
            }
        \end{center}
        \caption{The network from Example~\ref{exone}}
        \label{figone}
    \end{figure}
    The sender $v_1$ generates 3 independent random processes. The vertices
    $v_{11}$, $v_{12}$ and $v_{13}$ are the receivers. We will apply network
    coding over various fields of characteristic two. We start by considering
    random linear network coding over fields of size at least $4$. As $4>|T|=3$
    we know that this can be done successfully.

    We set $a_{1,1}=a_{2,2}=a_{3,3}=1$ and $a_{i,j}=0$ in all other cases.  We
    let $f_{i,j}=1$ except $f_{4,13}, f_{7,13}, f_{5,14}, f_{8,14}, f_{10,14},
    f_{9,15}, f_{11,15}$, which we choose by random. As in the last section we
    consider $b_{t,i,j}$ as fixed but unknown to us. Therefore
    $\tilde{P}=\widehat{P}$ is a polynomial in the seven variables  $f_{4,13},
    f_{7,13}, f_{5,14},f_{8,14},f_{10,14},f_{9,15},f_{11,15}$.  The determinant
    polynomial becomes
    \begin{equation*}
        \widehat{P}=(abcdefg+abce^2f^2+b^2c^2efg)Q
        ,
    \end{equation*}
    where
    \begin{align*}
        a &= f_{4,13} & b &= f_{5,14}  & c &= f_{7,13}  & d &= f_{8,14}\\
        e &= f_{9,15} & f &= f_{10,14} & g &= f_{11,15}
    \end{align*}
    and $Q=|B_{v_{11}}'| \, |B_{v_{12}}'| \, |B_{v_{13}}'|$. Here, $B_{v_{11}}'$
    respectively $B_{v_{12}}'$ respectively $B_{v_{13}}'$ is the matrix
    consisting of the nonzero columns of $B_{v_{11}}$ respectively the nonzero
    columns of $B_{v_{12}}$ respectively the nonzero columns of $B_{v_{13}}$.
    Choosing a lexicographic ordering with $d$ being larger than the other
    variables and applying~\eqref{eqtrio} we get that the success probability is
    at least
    \begin{equation*}
        P_\text{new 1}(q)=\frac{(q-1)^7}{q^7}
        .
    \end{equation*}
    Applying~\eqref{eqtri} we see that the success probability is at least
    \begin{equation*}
        P_\text{new 2}(q)=\frac{(q-1)^3(q-2)^2}{q^5}
        .
    \end{equation*}
    For comparison~\eqref{eqhosnabel1} tells us that success probability is at
    least
    \begin{equation*}
        P_{\text{Ho}}(q)=\frac{(q-3)^3}{q^3}
        .
    \end{equation*}
    Both bound~\eqref{eqtrio} and bound~\eqref{eqtri}
    exceed~\eqref{eqhosnabel1} for all values of $q \geq 4$.  In
    Table~\ref{tab1} we list $P_\text{new 1}(q)$, $P_\text{new 2}(q)$ and
    $P_\text{Ho}(q)$ for various values of $q$.

    \begin{table}
        \centering
        \caption{From Example~\ref{exone}: Estimates on the success
        probability}
        \label{tab1}
        \begin{tabular}{r ccccc}
            \toprule
            $q$ & 4 & 8 & 16 & 32 & 64
            \\
            $P_\text{new 1}(q)$ & 0.133 & 0.392 & 0.636 & 0.800 & 0.895
            \\
            $P_\text{new 2}(q)$ & 0.105 & 0.376 & 0.630 & 0.799 & 0.895
            \\
            $P_\text{Ho}(q)$  & $0.156 \times 10^{-1}$ & 0.244 & 0.536 & 0.744
            & 0.865
            \\
            \bottomrule
        \end{tabular}
    \end{table}

    We next consider the field $\F_2$. We reduce $\tilde{P}$ modulo $(a^2-a,
    \ldots, g^2-g)$ to get
    \begin{equation*}
        \widehat{P} = (abcdefg + abcef + bcefg) Q
        .
    \end{equation*}
    From~\eqref{eqtri} we see that the success probability of random network
    coding is at least $2^{-7}$. Choosing a proper monomial ordering we get
    from~\eqref{eqtrio}  that the success probability is at least $2^{-5}$. For
    comparison neither \cite{ho06}, \cite{jaggi03}, nor \cite{sanders03} tells
    us that linear network coding is possible.
\end{example}

\section{The Topological Meaning of $|M_t|$}\label{secpath}

Recall from Section~\ref{secho} that Ho et al.'s bound~\eqref{eqhosnabel1}
relies on the rather rough Lemma~\ref{lemho3}. The following theorem gives a
much more precise description of which monomials can occur in the support of $P$
and $\tilde{P}$ by explaining exactly which monomials can occur in $|M_t|$.
Thereby the theorem gives some insight into when the bounds \eqref{eqtrio} and
\eqref{eqtri} are much better than the bound~\eqref{eqhosnabel1}. The theorem
states that if $K$ is a monomial in the support of $|M_t|$ then it is the
product of $a_{i,j}$'s, $f_{i,j}$'s and $b_{t,i,j}$'s related to $h$ edge
disjoint paths $P_1, \ldots, P_h$ that originate in the senders and end in
receiver $t$.

\begin{theorem}\label{thepath}
    Consider a delay-free acyclic network.  If $K$ is a monomial in the support
    of the determinant of $M_t$ then it is of the form $K_1\cdots K_h$ where
    \begin{equation*}
        K_u = a_{u,l^{(u)}_1}f_{l_1^{(u)},l_2^{(u)}}  f_{l_2^{(u)},l_3^{(u)}}
        \cdots f_{l^{(u)}_{s_u-1},l^{(u)}_{s_u}}b_{t,v_u,l^{(u)}_{s_u}}
    \end{equation*}
    for $u=1, \ldots, h$. Here, $\{v_1, \ldots, v_h\} = \{1, \ldots, h\}$ holds
    and $l_1^{(1)}, \ldots,l_h^{(h)}$ respectively $l_{s_1}^{(1)}, \ldots,
    l_{s_h}^{(h)}$ are pairwise different. Further
    \begin{equation*}
        f_{l^{(u_1)}_i,l^{(u_1)}_{i+1}} \neq f_{l^{(u_2)}_j,l^{(u_2)}_{j+1}}
    \end{equation*}
    unless $u_1=u_2$ and $i=j$ hold. In other words $K$ corresponds to a product
    of $h$ edge disjoint paths.
\end{theorem}

\begin{proof}
    A proof can be found in the appendix.
\end{proof}

We illustrate the theorem with an example.

\begin{example}\label{exbutter}
    Consider the butterfly network in Figure~\ref{figurbutter}.
    \begin{figure}
        \begin{center}
            \scalebox{0.909}{
            $
            \psmatrix[colsep=0.8cm,rowsep=0.4cm,mnode=circle]
            &s\\
            v_1&&v_2\\
            &v_3\\
            &v_4\\
            t_1&&t_2
            \psset{arrowscale=2}
            \ncline{->}{1,2}{2,1}^{1}
            \ncline{->}{2,1}{3,2}^{4}
            \ncline{->}{1,2}{2,3}^{2}
            \ncline{->}{2,3}{3,2}^{5}
            \ncline{->}{2,1}{5,1}>{3}
            \ncline{->}{2,3}{5,3}<{6}
            \ncline{->}{3,2}{4,2}<{7}
            \ncline{->}{4,2}{5,1}^{8}
            \ncline{->}{4,2}{5,3}^{9}
            \endpsmatrix
            $
            }
        \end{center}
        \caption{The butterfly network}
        \label{figurbutter}
    \end{figure}
    A monomial $K$ is in the support of $|M_{t_1}|$ if and only if it is in the
    support of the determinant of
    \begin{multline*}
        N_{t_1}
        = (n_{i,j})
        = \begin{bmatrix}
            I+F&B_{t_1}^T\\
            A &0
        \end{bmatrix}
        \\
        = \begin{bmatrix}
            1&0&f_{1,3}&f_{1,4}&0&0&0&0&0&0&0\\
            0&1&0&0&f_{2,5}&f_{2,6}&0&0&0&0&0\\
            0&0&1&0&0&0&0&0&0&b_{t_1,1,3}&b_{t_1,2,3}\\
            0&0&0&1&0&0&f_{4,7}&0&0&0&0\\
            0&0&0&0&1&0&f_{5,7}&0&0&0&0\\
            0&0&0&0&0&1&0&0&0&0&0\\
            0&0&0&0&0&0&1&f_{7,8}&f_{7,9}&0&0\\
            0&0&0&0&0&0&0&1&0&b_{t_1,1,8}&b_{t_1,2,8}\\
            0&0&0&0&0&0&0&0&1&0&0\\
            a_{1,1}&a_{1,2}&0&0&0&0&0&0&0&0&0\\
            a_{2,1}&a_{2,2}&0&0&0&0&0&0&0&0&0
        \end{bmatrix}
    \end{multline*}
    By inspection we see that the monomial
    \begin{equation*}
        K = a_{1,1} a_{2,2} b_{t_1,1,3} b_{t_1,2,8} f_{7,8} f_{5,7} f_{2,5}
        f_{1,3}
    \end{equation*}
    is in the support of $|N_{t_1}|$.  We can write $K = K_1 K_2$ where
    \begin{equation*}
        K_1=a_{1,1}f_{1,3}b_{t_1,1,3}
        \quad \text{and} \quad
        K_2=a_{2,2}f_{2,5}f_{5,7}f_{7,8}b_{t_1,2,8}
        .
    \end{equation*}
    This is the description guaranteed by Theorem~\ref{thepath}. To make it
    easier for the reader to follow the proof of Theorem~\ref{thepath} in the
    appendix we now introduce some of the notations to be used there. By
    inspection the monomial $K$ can be written
    \begin{equation*}
        K = \prod_{i=1}^{11} n_{i,p(i)}
    \end{equation*}
    where the permutation $p$ is given by
    \begin{align*}
        p(1) &= 3 & p(2) &= 5  & p(3) &= 10 & p(4)  &= 4 & p(5)  &= 7 & p(6)=6\\
        p(7) &= 8 & p(8) &= 11 & p(9) &= 9  & p(10) &= 1 & p(11) &= 2
    \end{align*}
    Therefore if we index the elements in $\{1, \ldots, 11\}$ by
    \begin{align*}
        i_1 &= 10 & i_2 &=1 & i_3 &= 3 & i_4    &= 11 & i_5    &= 2 & i_6 &= 5\\
        i_7 &= 7  & i_8 &=8 & i_9 &= 4 & i_{10} &= 6  & i_{11} &= 9
    \end{align*}
    then we can write
    \begin{align*}
        K_1 &= n_{i_1,p(i_1)} n_{i_2,p(i_2)} n_{i_3,p(i_3)}
        \\
        K_2 &= n_{i_4,p(i_4)} n_{i_5,p(i_5)} n_{i_6,p(i_6)} n_{i_7,p(i_7)}
        n_{i_8,p(i_8)}
    \end{align*}
    and we have
    \begin{equation*}
        n_{i_9,p(i_9)}=n_{i_{10},p(i_{10})}=n_{i_{11},p(i_{11})} = 1
    \end{equation*}
    corresponding to the fact $p(i_9)=i_9$, $p(i_{10})=i_{10}$ and
    $p(i_{11})=i_{11}$.
\end{example}

\begin{remark}
    The procedures described in the proof of Theorem~\ref{thepath} can be
    reversed. This implies that there is a bijective map between the set of edge
    disjoint paths $P_1, \ldots, P_h$ in Theorem~\ref{thepath} and the set of
    monomials in $|M_t|$.
\end{remark}

Theorem~\ref{thepath} immediately applies to the situation of random network
coding if we plug into the $a_{i,j}$'s and into the $f_{t,i,j}$'s on the paths
$P_1, \ldots, P_h$ the fixed values wherever such are given. Let as in
Lemma~\ref{lemho3} $\eta$ be the number of edges for which some coefficients
$a_{i,j}, f_{i,j}$ are to be chosen by random. Considering the determinant as a
polynomial in the variables to be chosen by random with coefficients in the
field of rational expressions in the $b_{t,i,j}$'s we see that no monomial can
contain more than $\eta$ variables and that no variable occurs more than once.
This is because the paths $P_1, \ldots, P_h$ are edge disjoint. Hence,
Lemma~\ref{lemho3} is a consequence of Theorem~\ref{thepath}.

\section{Acknowledgments}

The authors would like to thank the anonymous referees for their helpful
suggestions.

\appendix

\section{Proof of Theorem~\ref{thepath}}

The proof of Theorem~\ref{thepath} calls for the following technical lemma.
\begin{lemma}\label{lemnocycle}
    Consider a delay-free acyclic network with corresponding matrix $F$ as in
    Section~\ref{secpre}. Let $I$ be the $|E| \times |E|$ identity matrix and
    define
    \begin{equation*}
        \Gamma=(\gamma_{i,j})=I+F
        .
    \end{equation*}
    Given a permutation $p$ on $\{1, \ldots, |E|\}$ write
    \begin{equation*}
        p^{(i)}(\lambda) = \overbrace{p(p(\cdots(\lambda )\cdots ))}^\text{i
        times}
    \end{equation*}
    If for some $\lambda \in \{1, \ldots, |E|\}$ the following hold
    \begin{enumerate}[label=(\arabic*)]
        \item \label{enu:lemma:a}
            $\lambda,p(\lambda), \ldots,p^{(x)}(\lambda)$ are pairwise different

        \item \label{enu:lemma:b}
            $p^{(x+1)}(\lambda) \in \{ \lambda, p(\lambda), \ldots,
            p^{(x)}(\lambda) \}$
        \item \label{enu:lemma:c}
            $\gamma_{\lambda,p(\lambda)}, \gamma_{p(\lambda),p(p(\lambda))},
            \ldots, \gamma_{p^{(x)}(\lambda),p^{(x+1)}(\lambda)}$ are all
            nonzero
    \end{enumerate}
    then $x=0$.
\end{lemma}

\begin{proof}
    Let $p$ be a permutation and let $x$ and $\lambda$ be numbers such that
    \ref{enu:lemma:a}, \ref{enu:lemma:b} and \ref{enu:lemma:c} hold. As $p$ is a
    permutation then \ref{enu:lemma:a} and \ref{enu:lemma:b} implies that
    $p(p^{(x)}(\lambda))=\lambda$. Aiming for a contradiction assume $x>0$. As
    $p(\eta)=\eta$ does not hold for any $\eta \in \{\lambda, p(\lambda),
    \ldots, p^{(x)}(\lambda)\}$,
    \begin{equation*}
        \gamma_{\lambda,p(\lambda)}, \gamma_{p(\lambda),p^{(2)}(\lambda)},
        \ldots, \gamma_{p^{(x)}(\lambda),p^{(x+1)}(\lambda)}
    \end{equation*}
    are all non-diagonal elements in $I+F$. By \ref{enu:lemma:c} we therefore
    have constructed a cycle in a cycle-free graph and the assumption $x>0$
    cannot be true.
\end{proof}

\begin{proof}[of Theorem~\ref{thepath}]
    A monomial is in the support of the determinant of $M_t$ if and only if it
    is in the support of the determinant of
    \begin{equation*}
        N_t =
        \begin{pmatrix}
            I+F&B_t^T\\
            A&0
        \end{pmatrix}
        = (n_{i,j})
        .
    \end{equation*}
    To ease the notation in the present proof we consider the latter matrix.
    Let $p$ be a permutation on $\{1, \ldots, |E|+h\}$ such that
    \begin{equation}
        \label{diamond1}
        \prod_{s=1}^{|E|+h}n_{s,p(s)} \neq 0
        .
    \end{equation}
    Below we order the elements in $\{1, \ldots, |E|+h\}$ in a particular way by
    indexing them $i_1, \ldots, i_{|E|+h}$ according to the following set of
    procedures.

    Let $i_1=|E|+1$  and define recursively
    \begin{equation*}
        i_s=p(i_{s-1})
    \end{equation*}
    until $|E| <p(i_s) \leq |E|+h$. Note that this must eventually happen due to
    Lemma~\ref{lemnocycle}.  Let $s_1$ be the (smallest) number such that $|E| <
    p(i_{s_1}) \leq |E|+h$ holds. This corresponds to saying that $n_{i_1,
    p(i_1)}$ is an entry in $A$, that $n_{i_2, p(i_2)},$ $\ldots, n_{i_{s_1-1},
    p(i_{s_1-1})}$ are entries in $I+F$, and that $n_{i_{s_1},p(i_{s_1})}$ is an
    entry in $B_t^T$.  Observe, that $p(i_r)=i_r$ cannot happen for $2 \leq r
    \leq s_1$ as already $p(i_{r-1})=i_r$ holds. As $n_{i_r,p(i_r)}$ is non-zero
    by~\eqref{diamond1} we therefore must have
    \begin{equation*}
        n_{i_r,p(i_r)}=f_{i_{r},p(i_r)}=f_{i_{r},i_{r+1}}
    \end{equation*}
    for $2\leq r < s_1$.  Hence,
    \begin{equation*}
        (n_{i_1,p(i_1)}, \ldots, n_{i_{s_1},p(i_{s_1})})
        = (a_{1,i_2},f_{i_2,i_3}, \ldots, f_{i_{s_1-1}, i_{s_1}},
        b_{t,v_1,i_{s_1}})
    \end{equation*}
    for some $v_1$. Denote this sequence by $P_1$. Clearly, $P_1$ corresponds to
    the polynomial $K_1$ in the theorem.

    We next apply the same procedure as above starting with $i_{s_1+1}=|E|+2$ to
    get a sequence $P_2$ of length $s_2$. Then we do the same with
    $i_{s_1+s_2+1}=|E|+3, \ldots, i_{s_1+\cdots s_{h-1}+1}=|E|+h$ to get the
    sequences $P_3, \ldots, P_h$. For $u=2, \ldots, h$ we have
    \begin{multline*}
        P_u
        = \left(
            n_{i_{s_1+\cdots +s_{u-1}+1}, p(i_{s_1+\cdots +s_{u-1}+1})}, \ldots
            , n_{i_{s_1+\cdots +s_{u}}, p(i_{s_1+\cdots +s_{u}})}
        \right)
        \\
        = \left(
            a_{u,i_{s_1+ \cdots +s_{u-1}+2}},
            f_{i_{s_1+\cdots +s_{u-1}+2}, i_{s_1+\cdots +s_{u-1}+3}}, \ldots,
            \right.
        \\
            \left.
            \phantom{a_{i_{s_{u-1}}}}
            f_{i_{s_1+\cdots +s_{u}-1}, i_{s_1+\cdots +s_{u}}},
            b_{t,v_u,i_{s_1+\cdots +s_{u}}}
        \right)
        .
    \end{multline*}
    Clearly, $P_u$ corresponds to $K_u$ in the theorem. Note that the sequences
    $P_1, \ldots, P_h$ by the very definition of a permutation are edge disjoint
    in the sense that
    \begin{enumerate}[label=(\arabic*)]
        \item $n_{i,j}$ occurs at most once in $P_1, \ldots, P_h$,
        \item if $n_{j,l_1}, n_{j,l_2}$ occur in $P_1, \ldots, P_h$ then
            $l_1=l_2$,
        \item if $n_{j_1,l}, n_{j_2,l}$ occur in $P_1, \ldots, P_h$ then
            $j_1=j_2$.
    \end{enumerate}
    Having indexed $s_1+\cdots + s_h$ of the integers in $\{1, \ldots, |E|+h\}$
    we consider what is left, namely
    \begin{equation*}
        \Lambda=\{ 1, \ldots, |E|+h\} \setminus \{i_1, \ldots, i_{s_1+\ldots
        +s_h}\}
        .
    \end{equation*}
    By construction we have $i_1=|E|+1$, $\ldots, i_{s_1+\cdots
    +s_{h-1}+1}=|E|+h$ and therefore $\Lambda \subseteq \{1, \ldots, |E|\}$.
    Also by construction for every
    \begin{equation*}
        \delta \in \{ 1, \ldots, |E| \} \cap \{ i_1, \ldots, i_{s_1+\cdots
        +s_h} \}
    \end{equation*}
    we have $\delta=p(\epsilon)$ for some $\epsilon \in \{i_1, \ldots, i_{s_1+
    \cdots +s_h} \}$. Therefore $p(\lambda) \in \Lambda$ for all $\lambda \in
    \Lambda$ holds. In particular $p^{(x)}(\lambda) \in \{1, \ldots, |E|\}$ for
    all $x$. From Lemma~\ref{lemnocycle} we conclude that $p(\lambda)=\lambda$
    for all $\lambda \in \Lambda$.
\end{proof}


TITLE: Decomposing factorized entropy
QUESTION [0 upvotes]: I am trying to figure out how the equation for factorized entropy below is derived. The equation for entropy is $H(Q) = -\sum_x Q(x)\log Q(x)$ where $Q$ is a probability distribution. Let $Q(x) = \prod_i Q_i(x_i)$ (a factorized distribution). Then, this is the derivation(see page 4, slide 1) I have seen for factorized entropy:

$$ H(Q)= - \sum_x Q(x)\log Q(x)$$
  $$=-\sum_x(\prod_{i\in x} Q(x_i))\log(\prod_{x_i \in x} Q(x_i))$$
  $$=-\sum_x(\prod_{i} Q(x_i))\sum_i \log Q(x_i))$$ $$=-\sum_i\sum_{x_i}Q(x_i)\log Q(x_i)$$

I can't figure out how the last line was derived(replacement of $\prod_iQ(x_i)$ by $\sum_{x_i}Q(x_i)$) . What confuses me is that $ab (\log a + \log b) \equiv \log (ab)^{ab}$ is not same as $(a+b)(\log a + \log b) \equiv \log (ab)^{a+b}$. I will appreciate pointers on how the last line is derived.

REPLY [0 votes]: $\text{Line 3 }=-\sum_x(\Pi_jQ_j(x_j))\sum_i\log Q_i(x_i)=-\sum_xQ(x)\sum_i\log Q_i(x_i)=-\sum_i\sum_xQ(x)\log Q_i(x_i)=-\sum_i\sum_{x_i,x_{-i}}Q_{-i}(x_{-i})Q_i(x_i)\log Q_i(x_i)=-\sum_i(\sum_{x_{-i}}Q_{-i}(x_{-i}))(\sum_{x_i}Q_i(x_i)\log Q_i(x_i))=-\sum_i\sum_{x_i}Q_i(x_i)\log Q_i(x_i)=\text{last line}$
$x_{-i}$ means non-$x_i$ terms.
$\sum_{x_i,x_j,x_k}Q_i(x_i)Q_j(x_j)Q_k(x_k)=\sum_{x_i}Q_i(x_i)\sum_{x_j}Q_j(x_j)\sum_{x_k}Q_k(x_k)$, as long as $x_i,x_j,x_k$ are pairwise/mutually independent, which is necessary for $Q$ to be factorized.


text:
\begin{document}

\maketitle

\begin{abstract}
In a recent work (arXiv-DOI: 1804.08072v1) we introduced the Modified Augmented Lagrangian Method (MALM) for the efficient minimization of objective functions with large quadratic penalty terms.
From MALM there results an optimality equation system that is related to that of the original objective function. But, it is numerically better behaved, as the large penalty factor is replaced by a milder factor.

In our original work, we formulated MALM with an inner iteration that applies a Quasi-Newton iteration to compute the root of a multi-variate function. In this note we show that this formulation of the scheme with a Newton iteration can be replaced conveniently by formulating a well-scaled unconstrained minimization problem.

In this note, we briefly review the Augmented Lagrangian Method (ALM) for minimizing equality-constrained problems. Then we motivate and derive the new proposed formulation of MALM for minimizing unconstrained problems with large quadratic penalties. Eventually, we discuss relations between MALM and ALM.
\end{abstract}

\noindent
\subsubsection*{Brief summary}
We propose the iteration in Section~\ref{sec:4.1}\,. The iteration is used to solve problem \eqref{eqn:UPNP}. The iteration is also meaningful and well-behaved when $\omega_E \rightarrow +0$.

\section{Introduction}
ALM is an iterative method that is made for solving equality-constrained problems \cite[Chapter~17]{NumOpt}. The benefit is that ALM solves the original problem by minimizing a sequence of unconstrained functions that have mild quadratic penalties.

In \cite{MALM} we introduced an approach to using ALM for the solution of unconstrained functions with large quadratic penalties. This scheme we called MALM. MALM provides the same benefit as ALM. Namely, MALM solves the original problem by minimizing a sequence of unconstrained functions that have mild quadratic penalties. Hence, the large quadratic penalties are replaced by smaller ones. This is beneficial for convergence.

In our original presentation, MALM was formulated with a Quasi-Newton system: A new iterate was constructed by computing the root of a multi-variate function. In this work, through careful analysis of this multi-variate function, we are able to characterize the root of this function as a critical point for an unconstrained minimization problem. This allows us to formulate the iteration of MALM in a way that is more similar to the iteration of ALM. Hence, we are able to compare both ALM and MALM in a more elegant and meaningful way. Also, the new formulation makes the iteration of MALM simpler, which is desirable when attempting to incorporate it into other existing optimization methods.

\section{Background}
Consider the \textit{equality-constrained nonlinear program}:
\begin{equation}
	\tag{ECNP}\label{eqn:ECNP}
	\begin{aligned}
		&\operatornamewithlimits{min}_{\bx \in \R^n} 	&\quad 	f(\bx)& 	\\
		&\text{subject to} 						&		c(\bx)&=\bO \in \R^m
	\end{aligned}
\end{equation}
with the Lagrangian function $\cL(\bx,\blambda) := f(\bx)- \blambda\t \cdot c(\bx)$.

ALM is an iterative method for solving \eqref{eqn:ECNP}. The method is given a sufficiently small penalty parameter $\omega >0$, an initial guesses $\bx_0 \in \R^n$ for the primal solution, and an initial guess $\blambda_0 \in \R^m$ for the dual solution. The iteration of ALM then generates sequences $\lbrace\bx_k\rbrace_{k \in \N_0}$, $\lbrace\blambda_k\rbrace_{k\in\N_0}$, where for $k \in \N$ it holds:
\begin{enumerate}
	\item $\bx_k$ minimizes
	\begin{align*}
		f_k(\bx) := \cL(\bx,\blambda_{k-1}) + \frac{1}{2 \cdot \omega} \cdot \|c(\bx)\|_2^2
	\end{align*}
	and is computed with a numerical method for unconstrained minimization, using the initial guess $\bx_{k-1}$\,.
	\item After $\bx_k$ has been computed, the dual vector $\blambda_k$ is computed as
	\begin{align*}
		\blambda_k := \blambda_{k-1} - \frac{1}{\omega} \cdot c(\bx_k)\,.
	\end{align*}
\end{enumerate}
Under mild conditions on the problem \eqref{eqn:ECNP} and a sufficiently small choice for $\omega$, that however is finite and does not need to converge to zero, the iteration of ALM converges to a stationary point $\bx_\infty,\blambda_\infty$ of the following optimality equations due to Karush, Kuhn and Tucker (KKT) \cite{NumOpt,Lancelot}:
\begin{align*}
	\nabla_\bx \cL(\bx,\blambda) &= \bO\in\R^n\,,\\
	\nabla_\blambda \cL(\bx,\blambda) &= \bO\in\R^m\,.
\end{align*}

\section{Motivation of MALM}
Sometimes, instead of \eqref{eqn:ECNP}, we may rather want to minimize the following problem, that we call \textit{unconstrained penalty nonlinear program}:
\begin{equation}
	\tag{UPNP}\label{eqn:UPNP}
	\begin{aligned}
		&\operatornamewithlimits{min}_{\bx \in \R^n} 	&\quad 	f(\bx) + \frac{1}{2 \cdot \omega_E} \cdot \|c(\bx)\|_2^2
	\end{aligned}
\end{equation}
where typically $\omega_E>0$ is a very small user-defined parameter.

There are various reasons why, for given functions $f,c$, one may rather prefer solving \eqref{eqn:UPNP} than \eqref{eqn:ECNP}.
\begin{enumerate}
	\item Problem \eqref{eqn:UPNP} is always feasible. Thus, we do not run into any numerical issues related to infeasibility. This is beneficial when the user simply asks for a solution under any circumstances.
	\item It is impossible to verify numerically whether \eqref{eqn:ECNP} is feasible. This is so already because of the reason that in all but very few cases the equations $c(\bx)=\bO$ do not have a digitally representable solution for $\bx$.
	\item Unless assuming a constraint qualification (e.g. LICQ), the optimality equations of \eqref{eqn:ECNP} are ill-posed, and hence do not admit a meaningful numerical treatment. In contrast, \eqref{eqn:UPNP} is always well-posed, with a condition number that can be controlled by the parameter $\omega_E>0$.
	\item There are applications where $m>n$. This immediately causes violation of all constraint qualifications. However, \eqref{eqn:ECNP} may still possess feasible points in such cases. In these situations, it constitutes a sound approach to generate numerical solutions of \eqref{eqn:ECNP} by solving \eqref{eqn:UPNP} for very small values of $\omega_E>0$.
	\item In general, the numerical objective of satisfying the constraints $c$ is in contradiction to the numerical objective of minimizing $f$. Hence, from a numerical point of view, it is impossible to think of any method that does not actually minimize a bias (merit function) between $f$ and a norm of $c$. In \eqref{eqn:UPNP}, at least we know this merit function and can analyze its meaning with respect to an application.
	\item In logical consequence of the former bullet, many numerical methods for treating \eqref{eqn:ECNP} are actually based on minimizing \eqref{eqn:UPNP} for decreasing values of $\omega_E>0$. The methods described in the following literature serve as examples: \cite{Amudson,ChenGoldfarb,ForsgrenGill}.
	\item Finally and not yet mentioned, there simply exist problems where people literally wish to minimize a function of the form \eqref{eqn:UPNP}. Examples arise from applications where $n,m$ are dimensions of discretized infinite spaces for $\bx,\blambda$. This is the case, e.g., in constrained optimization with partial differential equations. It is for instance in these situations, where it can make sense to consider discretizations where $m>n$, and treat the resulting problem as \eqref{eqn:UPNP}. To give an example, in \cite{myOCP} we present a discretization for optimal control problems that results in nonlinear programs of the form \eqref{eqn:UPNP}.
\end{enumerate}

If we were to minimize \eqref{eqn:UPNP} with a numerical method for unconstrained minimization, then there would be an issue. This issue would be, that the large quadratic penalty terms in the objective impede the rate of convergence of the numerical method. Hence, it seems desirable to apply a technique such as ALM, that manages to solve a minimization problem with only mild quadratic penalties; i.e., it only uses quadratic penalties with values $\omega \gg \omega_E$. Usually, these milder choices of $\omega$ form well-scaled objectives and hence promotive fast convergence.

However, it is unclear how ALM could be applied to treat \eqref{eqn:UPNP}. This is because after all, \eqref{eqn:UPNP} and \eqref{eqn:ECNP} are two entirely different problem statements.

\largeparbreak

In \cite{MALM}, we made an approach to reformulating \eqref{eqn:UPNP} as an instance of \eqref{eqn:ECNP}, by introducing an auxiliary variable $\bxi \in \R^m$:
\begin{subequations}
	\label{eqn:Subst}
	\begin{align}
		&\operatornamewithlimits{min}_{\bx \in \R^n\,,\bxi \in \R^m} 	&\quad 	f(\bx) + \frac{\omega_E}{2} \cdot \|\bxi\|_2^2& \label{eqn:SubstObjective}	\\
		&\text{subject to} 						&		c(\bx) + \omega_E \cdot \bxi &=\bO \in \R^m \label{eqn:SubstConstraints}
	\end{align}
\end{subequations}
This is advantageous, because the new objective \eqref{eqn:SubstObjective} is well-scaled in $\bx,\bxi$ even for the smallest values of $\omega_E$, and the constraints \eqref{eqn:SubstConstraints} are well-scaled in $\bx,\bxi$ as well and feasible for all values of $\omega_E>0$.

\largeparbreak

In \cite{MALM}, after formulating \eqref{eqn:Subst}, we then applied ALM to \eqref{eqn:Subst}. We arrived at an optimality system that was linear in $\bxi$. We eliminated $\bxi$ from that system. We were then able to establish an iterative scheme, that we called MALM. This scheme is effectively the iteration of ALM applied to \eqref{eqn:Subst}.

The iteration goes as follows. Provided are initial guesses $\bx_0,\blambda_0$, and a sufficiently small value $\omega \gg \omega_E$. The iteration of MALM then generates sequences $\lbrace\bx_k\rbrace_{k \in \N_0}$, $\lbrace\blambda_k\rbrace_{k\in\N_0}$, where for $k \in \N$ it holds:
\begin{enumerate}
	\item $\bx_k$, together with an auxiliary variable $\tblambda \in \R^m$, solves the following equation system
	\begin{align*}
		F_k(\bx,\tblambda) := \begin{pmatrix}
			\nabla_\bx \cL(\bx,\blambda_{k-1} + \tblambda)\\
			c(\bx) + \omega_E \cdot \blambda_{k-1} + (\omega + \omega_E) \cdot \tblambda
		\end{pmatrix}\,.
	\end{align*}
	In our original work \cite{MALM}, we proposed that the root $(\bx_k,\tblambda)$ of $F_k$ be computed using a particular Quasi-Newton line-search method. To this end, we presented both a suitable inertia-correction scheme for the step-direction obtained from the Newton system and a suitable merit-function for which this step-direction is a descent-direction. We referred to literature results that proved global and fast local convergence of this scheme.
	\item After $\bx_k,\tblambda$ has been computed, the dual vector $\blambda_k$ is computed as
	\begin{align*}
	\blambda_k := \blambda_{k-1} + \tblambda
	\end{align*}
\end{enumerate}

\section{Associated objective for MALM}
The iteration of MALM is a superior approach over minimizing \eqref{eqn:UPNP} with an iteration for unconstrained minimization. But, MALM has one detail, that makes it un-elegant and difficult to handle. We point out this detail in the following:

The iteration of ALM defines $\bx_k$ as the local unconstrained minimizer of $f_k$. In contrast to that, the iteration of MALM defines $\bx_k$ as the root of a particular optimality system. It would be nice if instead we could define $\bx_k$ in MALM as the unconstrained minimizer of a function like $f_k$, as well.

\largeparbreak

The proposal of this note is that we found a suitable function for $f_k$. Namely, we can interpret $F_k$ as the optimality system for the following unconstrained minimization problem:
\begin{equation}
	\tag{UMP$_k$}\label{eqn:UMP}
	\begin{aligned}
		&\operatornamewithlimits{min}_{\bx \in \R^n} 	&\quad 	f_k(\bx) := \cL(\bx,\blambda_k) + \frac{1}{2 \cdot (\omega+\omega_E)} \cdot \|c(\bx)+\omega_E \cdot \blambda_{k-1}\|_2^2
	\end{aligned}
\end{equation}
This minimality problem we found as follows: We eliminate the second row in $F_k$ and obtain for $\tblambda$ the formula
$$ 	\tblambda = \frac{-1}{\omega + \omega_E} \cdot \big(\,c(\bx_k) + \omega_E \cdot \blambda_{k-1}\,\big) 	$$
in terms of $\bx_k,\blambda_{k-1}$. We insert it into the first row of $F_k$. The first row of $F_k$ with $\tblambda$ inserted is
$$ 	\nabla_\bx \cL(\bx_k,\blambda_{k-1}) - \nabla c(\bx_k) \cdot \frac{-1}{\omega + \omega_E} \cdot \big(\,c(\bx_k) + \omega_E \cdot \blambda_{k-1}\,\big) \,.	$$
It turns out that this expression coincides with the gradient of our proposed choice for $f_k$. Hence, minimizers of $f_k$ are roots of $F_k$.

\subsection{Proposed new form for MALM}\label{sec:4.1}
The formulation of $f_k$ allows us to simplify our initial proposal for MALM into the following new iterative scheme, that is described below.

The iteration goes as follows. Provided are initial guesses $\bx_0,\blambda_0$, and a sufficiently small value $\omega \gg \omega_E$. The iteration of MALM then generates sequences $\lbrace\bx_k\rbrace_{k \in \N_0}$, $\lbrace\blambda_k\rbrace_{k\in\N_0}$, where for $k \in \N$ it holds:
\begin{enumerate}
	\item $\bx_k$ minimizes
	\begin{align*}
	f_k(\bx) := \cL(\bx,\blambda_{k-1}) + \frac{1}{2 \cdot (\omega+\omega_E)} \cdot \|c(\bx)+\omega_E \cdot \blambda_{k-1}\|_2^2
	\end{align*}
	and is computed with a numerical method for unconstrained minimization, using the initial guess $\bx_{k-1}$\,.
	\item After $\bx_k$ has been computed, the dual vector $\blambda_k$ is computed as
	\begin{align*}
	\blambda_k := \blambda_{k-1} - \frac{1}{\omega+\omega_E} \cdot \big(\,c(\bx_k)+\omega_E \cdot \blambda_{k-1}\,\big)\,.
	\end{align*}
\end{enumerate}
The second step remains unchanged from the originally proposed MALM. We just inserted the above expression for $\tblambda$ in order to vanish it and clean up the iteration.

\subsection{Discussion and benefits of the result}
We find that our proposed new MALM is a direct generalization of ALM. This is because for $\omega_E=0$ the two schemes coincide.

The benefit of our method is that it is both capable of solving \eqref{eqn:ECNP} and \eqref{eqn:UPNP}, while ALM can only treat the latter. If we choose $\omega_E>0$, then MALM solves \eqref{eqn:UPNP} for that respective value of $\omega_E$. If instead we choose $\omega_E=0$, then MALM coincides with ALM and solves \eqref{eqn:ECNP}. Hence, for every $\omega_E \in [0,\infty)$ the iteration is well-defined. Therefore, MALM offers a convenient treatment for solving problems of both classes, and for solving problems \eqref{eqn:UPNP} with literally arbitrarily small values of $\omega_E$.

From a numerical perspective, the proposed new MALM iteration appears as suitable for computations as does ALM. We justify this: In ALM, for the update of the Augmented Lagrangian in step 2, there is a division through a small parameter $\omega$. In the new proposed MALM, this divisor is augmented with $\omega_E$. Certainly, this augmentation does not introduce any numerical difficulty, because it only increase of the divisor. Hence, we believe that MALM resembles the desirable practical convergence behavior of ALM, while being applicable to a more general setting.

\section{Summary}
We have modified the Augmented Lagrangian Method in a small detail. Namely, we added a parameter $\omega_E \in [0,\infty)$ into the iteration formulas. This modification allows us to treat problems of either form \eqref{eqn:ECNP} or \eqref{eqn:UPNP}.

Our Modified Augmented Lagrangian Method is numerically as well-behaved as the original Augmented Lagrangian Method. This relates to numerical well-posedness and to the rate of convergence. Our modified method is insensitive to small choices of $\omega_E>0$ and coincides with the original Augmented Lagrangian in the limit $\omega_E = 0$.

\largeparbreak

Further research may be related to the following subject. We understand that --\,whenever existent\,-- a solution of \eqref{eqn:ECNP} can be approached by solving a sequence of \eqref{eqn:UPNP} for $\omega_E \searrow +0$\,. Now, with our Modified Augmented Lagrangian Method, we even have a convenient tool at hand that can solve \eqref{eqn:UPNP} in the limit $\omega_E = 0$\,. This motivates research investigations in the following direction: What solution does the sequence \eqref{eqn:UPNP} for $\omega_E \searrow +0$ converge to when \eqref{eqn:ECNP} does not have any solutions. And, related to this situation, to which values does the (Modified) Augmented Lagrangian method with $\omega_E \searrow +0$ converge to? Or, does the limit exist in these cases?

Certainly, a progress in this regard is provided in our analysis of a method for convex quadratic programming. We show that for $f$ convex quadratic and $c$ linear the solution of \eqref{eqn:UPNP} for $\omega_E \searrow +0$ converges to a solution of the bi-objective minimization problem
\begin{equation*}
	\begin{aligned}
		&\operatornamewithlimits{min}_{\bx \in \R^n} 	&\quad 	f(\bx)&\,,\\
		&\text{subject to} 								&\quad \|c(\bx)\|_2& \text{ is locally minimal\,.}
	\end{aligned}
\end{equation*}
We conjecture that a similar result can also be verified for the more general case of non-convex nonlinear programming.

\largeparbreak

Finally, justified by our proposed method, it would be desirable if people could please provide optimization software that is capable of treating problems \eqref{eqn:UPNP}, and that does not return on error when $m>n$\,.

\FloatBarrier

\bibliography{NOTE_MALM_bib}
\bibliographystyle{plain}

\end{document}

config: arxiv


text:
\begin{document}

\title{Distribution free goodness of fit tests for regularly varying tail distributions

}

\author{Thuong\ T.\ M.\ Nguyen}

\institute{Thuong T. M. Nguyen \at
              School of Mathematics and Statistics,
              Victoria University of Wellington\\ Wellington 6140, New Zealand \\
              \email{thuong.nguyen@vuw.ac.nz}

           }

\maketitle

\begin{abstract}
We discuss in this paper a possibility of constructing a whole class of asymptotic distribution-free tests for testing regularly varying tail distributions. The idea is that we treat the tails of distributions as members of a parametric family and using MLE to estimate the exponent. No matter what the exponent's estimator is, we are able to transform the whole class into a specific distribution with a prefix exponent so that we are free from choosing any functional of the tail empirical process as a distribution-free test statistic. The asymptotic behaviour of some new tests, as examples from the whole class of new tests, are demonstrated as well.

\keywords{Regularly varying \and distribution free \and goodness of fit test\and unitary transformation.}

\end{abstract}

\section{Introduction}

Suppose that we have a random sample $X_1,X_2,\cdots, X_n$ of a random variable $X$ following some unknown distribution $F$. No matter how $F$
behaves, we are
mainly interested in the right tail behavior of $F$, i.e., in $F(x)$ when $F(x) \to 1$. Among the class of heavy-tail distributions, regularly varying
tail distributions has been attracted most attention due to its various applications. Briefly, a distribution $F$ is said to be {\it regularly varying
in tail}
if
\begin{align}
\lim_{x\to\infty} \frac{1-F(tx)}{1-F(x)} = t^{-\theta}
\end{align}
for all $ t >0$ where $\theta >0$ is called the {\it index} or {\it exponent} of regular variation. By this
definition, regularly varying tail distributions are often indicated as Pareto-type distributions.

Regular variation of the tail of a distribution often appears as a natural condition in various theoretical results of probability theory.
The most typical example is that it is the condition for the distribution of the partial maxima
to belong to the domain of attraction of extreme value distributions. This could be found in various references, we refer to books of de Haan and Ferreira
\cite{deHaanbook} and Resnick \cite{Resnick} among others.
Moreover, regular variation naturally arises as a common phenomenon in numerous practical case studies such as finance, insurance, physics, geology,
hydrology and
engineering, etc. Therefore, detecting and testing regular variation plays an important role in probabilistic areas as well as in practical applications.

\medskip

The problem of estimating the exponent $\theta$ has a very rich literature. Among various studies, the Hill estimator is the most
common
one, see Hill
\cite{Hill}. Recall that the Hill estimator of the exponent $\theta$ essentially takes on the following form
\begin{align}
 \hat{\theta}^H = \Big(\frac{1}{k} \sum_{i=1}^k \ln X_{i,n} - \ln X_{k,n}\Big)^{-1}
\end{align}
where $k = k(n)$ is a certain sample fraction, $k(n) \to \infty$ in an appropriate way and $X_{1,n} \le \cdots \le X_{n,n}$ denotes the upper order
statistics
of the sample $X_1,\cdots, X_n$. Many other estimators was proposed, most of them is also based on the upper order statistics and is not too difficult to
compute. These include, to name a few, popular estimators introduced by Dekkers et. al. \cite{Dekkersetal}, de Haan and Resnick \cite{DeHaan},
Pickands \cite{Pickands}, Teugels \cite{Teugels} among others. The rate of the convergence of an estimator was also discussed, for example, in Hall and
Welsh \cite{Hall84}.

\medskip

In contrast to the numerous number of approaches to estimate the exponent $\theta$, using goodness of fit (gof) for
testing regular variation has been addressed in a modest number of studies.  Not long ago, Beirlant et al. \cite{Beirlantetall} modified
the Jackson statistic - which was originally proposed as a gof for testing exponentiality - for testing Pareto-type data. Koning and Peng \cite{Koning}
examined the Kolmogorov-Smirnov, Berk-Jones and the estimated score tests and compare them in terms of Bahadur efficiency. In that paper, the Berk-Jones
and estimated score tests, which are not based on the empirical process, were shown to perform better than the Kolmogorov-Smirnov test.

\medskip

Note that regularly varying tail distributions are members of the family of the generalized Pareto distribution (GPD), which is
of the form
\begin{align}\label{gpd}
 F(x;k,\sigma) =\begin{cases} 1 - (1-\frac{kx}{\sigma})^{1/k}, \quad k\ne 0, \sigma >0\\
 1 - e^{-x/\sigma}, \qquad\quad\,\, k=0, \sigma >0
 \end{cases}
\end{align}
where $k$ and $\sigma$ are shape and scale parameters and $x >0$. To test the fit of data to a GDP, there has been several studies such as Davison and
Smith \cite{Davison}, Choulakian and Stephens \cite{Choulakian}. In these papers, the critical values for Cramer-von Mises
statistic
and
Anderson-Darling statistic, which is a weighted Cramer-von Mises statistic, were given.

\medskip

In this paper, we will introduce a wide class of asymptotically
distribution free gof tests for testing regularly varying tail distributions. Our approach follows a
new method introduced in Khmaladze \cite{Estate2}. In that paper, the author proposed a unitary transformation which enables us to create a class of
gof tests for both simple and parametric hypothesis testing problems. The method for the latter will be adopted for our problem, which we will present in
Section \ref{method}. Briefly speaking, we will use exactly the same transformation in \cite{Estate2} to derive a modification of the empirical process.
Only this time, the assumption that $F$ belongs to a parametric family of distributions is no longer available and the fact is that the right tail just
partially describes the distribution $F$. Hence, we can only consider the tail empirical process and transform it. The transformed tail empirical process,
under the hypothesis of interest, possesses a limit in distribution asymptotically free from any underlying distribution $F$ as well as the unknown
index/exponent of the regular variation. Therefore, any appropriate functionals of the transformed process could be used as asymptotically distribution
free test statistics. Section \ref{review} will be devoted for the main literature of the empirical process based on the tail.

\medskip

Some simulation results will be documented in Section \ref{simulation}. Namely, we will take the Kolmogorov-Smirnov (KS), Cramer-von Mises ($\Omega^2$)
and Anderson-Darling ($A^2$)
tests as
examples from the new class of asymptotically distribution free gof tests for demonstration. The asymptotically distribution free property of these test
statistics will be illustrated for different choices of the original distribution $F$.

\section{Main results}\label{method}

\subsection{Overall review}\label{review}

Suppose that from a random sample $X_1, \cdots, X_n$ we are only considering excesses over a certain threshold $x_0$ which
is considerably large. Let us denote the subsample of all observations
exceeding $x_0$ by $\Xtil_1,\cdots,\Xtil_m$. Generally, the choice of $x_0$ as well as the sample fraction $m/n$ may require some educated guessing. Hill
\cite{Hill} suggested $m$ be chosen as an adaptive, data-analytic basis.
Several methods for choosing $m$ or $x_0$ based on survey data could be found in Drees and Kaufmann \cite{Drees}, Danielsson et al. \cite{Danielsson}
and
Guillou and Hall \cite{Guillou}. Nevertheless, there has existed various research assuming that $m$ is known
such that $m\to \infty$ as $n\to \infty$ but $m=o(n)$. With this assumption and some prior knowledge about the underlying distribution function
$F$, Haeusler and Teugels \cite{Haeusler} derived a general condition which can be used to determine the optimal value $m$ explicitly. Hall \cite{Hall82}
also considered $m$ having a deterministic value and a quite common estimate of the exponent was introduced. Throughout this
paper, our approach will be based on the same such assumption, that also means $x_0$ fixed.

\medskip

Our main aim is to create a class of gof test for testing the hypothesis $H_0:$ ``$F$ is a regularly varying tail distribution" against the alternative
$H_1:$ ``$F$ is not a regularly varying tail distribution". It is sensible that the exponent $\theta$ needs to be estimated based on the right tail only.
Specifically, we only consider observed values $\set{\Xtil_1,\cdots,\Xtil_m}$, which now should be looked at as a sample of a different random
variable, let say, $\Xtil$.

\medskip

Denote $\Ttil = \frac{\Xtil}{x_0}$ and $T =\Ttil -1$. Under the hypothesis of interest and by the definition of the regular variation, the survival
distribution of $\Xtil$ conditional on the specified large value $x_0$ is
\begin{align}
\frac{\Pds\set{\Xtil\ge x_0t}}{\Pds\set{\Xtil\ge x_0}} = \frac{1-F(x_0t)}{1-F(x_0)} = t^{-\theta} + o(1), \quad t\ge 1,\, \theta\ge 0.
\end{align}
Then, the
distribution of the positive continuous random variable $T$ under the null hypothesis is
\begin{align}
 H_{\theta}(t) = 1 - (1+t)^{-\theta}, \qquad t \ge 0.
\end{align}
Clearly, the density function is $h(t) = \theta(1+t)^{-(\theta+1)}.$ Viewing $T$ as a new random variable, we can define the tail empirical process in a
similar way as of the standard empirical process, i.e., we have
\begin{align}
H_m(t) = \frac{1}{m} \sum_{i=1}^m 1_{\set{T_i \le t}}.
\end{align}

\noindent The tail parametric
empirical process $\vhat_{mH}$ is
\begin{align}\label{empro}
 \vhat_{mH}(t) = \sqrt{m}[H_m(t) - H_{\thetahat_m}(t)].
\end{align}
where $\thetahat_m$ is an estimator of $\theta$ calculated from the subsample $\Xtil_1,\cdots,\Xtil_m$ (or $T_1,\cdots, T_m$). Assume that the true
unknown exponent under the hypothesis $H_0$ is $\theta_0$. Let $\thetahat_m$ be the maximum likelihood estimator (MLE), then it is
the solution of the equation
\begin{align}
\sum_{i=1}^m \frac{\partial \log h(T_i, \theta)}{\partial \theta} =0.
\end{align}
That yields
\begin{align}
\thetahat_m = \frac{m}{\sum_{i=1}^m \log(T_i +1)},
\end{align}
which actually coincides with the Hill estimator $\thetahat^H$. This estimator was proved to be consistent in the sense that
\begin{align}
 \thetahat^H = \thetahat_m \stackrel{\Pds}{\to} \theta_0
\end{align}
under the condition that $m\to \infty$ as $n\to\infty$ such that $m/n \to 0$. The proof of this
convergence
could be found in Mason \cite{Mason}. Regarding the asymptotic normality of the estimator, we refer to Haeusler and Teugels \cite{Haeusler}
, Geluk et al. \cite{Geluketal.}, de Haan and Resnick \cite{deHaanRes} among various others.

We will spend few more lines here to review the property of the limit in distribution of the process $\vhat_{mH}$ because it is
essential for the method we present below. For all the concepts and terminologies, we follow and keep the same as in \cite{Estate2}. For the sake of
lucidity, we extract and represent the method only for our particular problem.

Consider the space $\Lcal_2(H)$. Recall that a function $\phi$ is integrable with respect to $H$ if $\int_0^{\infty} \phi(s) H(ds) <\infty$ and square
integrable if
$\int_0^{\infty} \phi^2(s) H(ds) <\infty$. The space $\Lcal_2(H)$ consists of all square integrable functions with respect to $H$. The inner product and
norm in $\Lcal_2(H)$ are defined as usual. That is, for any $\phi,\phitil \in
\Lcal_2(H)$  we have
\begin{align*}
 \norm{\phi}_H^2 &= \int_0^{\infty} \phi^2(s) H(ds),\\
 \langle \phi, \tilde{\phi}\rangle_H &= \int_{0}^{\infty} \phi(s) \tilde{\phi}(s)H(ds).
\end{align*}

Denote by $w_H(\phi)$ a function-parametric $H$-Brownian motion where $\phi$ is a square integrable function in
$\Lcal_2(H)$. That
means,
$w_H(\phi)$ for each $\phi$ is a Gaussian random variable with expected value $0$ and variance $\Esf w^2_H(\phi) =
\norm{\phi}_H^2.$ This also implies that the covariance between $w_H(\phi)$ and $w_H(\tilde{\phi})$ is

\begin{align*}
 \Esf w_H(\phi)w_H(\tilde{\phi}) = \langle \phi, \tilde{\phi}\rangle_H.
\end{align*}
If $\phi_t(t^*) = 1_{\set{t^*\le t}}$ then $w_H(\phi_t):= w_H(t)$ is simply the Brownian motion in time $H(t).$ As usual, a linear transformation of
$w_H(t)$ which is of the form $v_H(t)=w_H(t) - H(t) w_H(\infty)$ is the Brownian bridge in time $H(t).$

Assuming that under the null hypothesis $H_0,$ the true unknown exponent of the regular variation is $\theta_0$. Denote by $v_{mH}$ the usual tail
empirical process, that is,
\begin{align}
 v_{mH}(t) = \sqrt{m}[H_m(t) - H_{\theta_0}(t)].
\end{align}

\noindent Let us stress that, since we assumed $x_0$ is known, the limit in distribution of
$v_{mH}(t)$, as a {\it conditional} tail empirical process, is
similar to that of the standard empirical process, which is the Brownian bridge $v_H(t)$. For a full description on the property of general tail empirical
process, we refer to Einmahl \cite{Einmahl90,Einmahl92}. It was proved in these papers that if such process is {\it unconditional} on $x_0$, the limit in
distribution of
$v_{mH}(t)$ is a Brownian motion in time
$H(t)$. On a recent approach on distribution free gof test for testing tail copula by Can et al. \cite{Canetal}, the tail empirical process - constructed
on tail copula, was mapped to a standard Brownian motion. Their construction
based on the innovative martingale method in Khmaladze \cite{Estate1} which is known as the Khmaladze transformation. Our approach here use a
different Khmaladze transformation in \cite{Estate2}, and so let us call it Khmaladze-2.

\medskip

Considers the function-parametric version of the tail empirical process
\begin{align}\label{func.em.pr}
v_{mH}(\phi) = \int_0^{\infty} \phi(t) v_{mH}(dt) = \frac{1}{\sqrt{m}} \sum_{i=1}^m [\phi(T_i) - \Esf \phi(T_i)]
\end{align}
where $\phi$ is a function in $\Lcal_2(H)$. In a similar way to achieve the limit of $v_{mH}(t)$, the limit in distribution of $v_{mH}(\phi)$ with some
proper restriction on functions $\phi$ is
called a
function-parametric $H$-Brownian bridge. Specifically,
\begin{align}\label{limitvH}
 v_H(\phi) = w_H(\phi) - \langle \phi, \mathbb{1} \rangle_H w_H(\mathbb{1}),
\end{align}
where $\mathbb{1}$ stands for the function identically equals to $1$.

\medskip

We can also say more about the asymptotic behavior of the tail parametric empirical process where $\thetahat_m$ is
the MLE which is also the Hill's estimator. It is well-known since long ago that under some usual and mild constraints, MLE
$\thetahat_m$ possesses the asymptotic property
\begin{align}\label{mle-es}
 \sqrt{m}(\thetahat_m -\theta_0) = \Gamma_{H}^{-1} \int_0^{\infty} \frac{\dot{h}_{\theta_0}(t)}{h_{\theta_0}(t)} v_{mH}(dt,\theta_0) + o_P(1),
\qquad m \to \infty
\end{align}
where $h_{\theta_0}$ denotes the hypothetical density and $\dot{h}_{\theta_0}$ its derivatives in $\theta$. Denote by
\begin{align}
 \Gamma_{H} = \int_0^{\infty} \frac{\dot{h}_{\theta_0}^2(t)}{h_{\theta_0}^2(t)} H(dt)
\end{align}
the Fisher information. Then it is easy to check that $\Gamma_{H} = 1/\theta_0^2.$ As a consequence of \eqref{mle-es}, we can expand the process
$\vhat_{mH}(t)$ as
\begin{align}\label{em.pr-expand}
 \vhat_{mH}(t) &= v_{mH}(t) - \sqrt{m}[H_{\thetahat_m}(t)-H_{\theta_0}(t)]\nonumber\\
 &= v_{mH}(t) - \int_0^t \frac{\dot{h}_{\theta_0}(t)}{h_{\theta_0}(t)} H_{\theta_0}(dt) \Gamma_{H}^{-1}  \int_0^{\infty}
\frac{\dot{h}_{\theta_0}(t)}{h_{\theta_0}(t)} v_{mH}(dt,\theta_0) + o_P(1)\nonumber\\
&= v_{mH}(t) - \int_0^t \beta_H^T(t) H_{\theta_0}(dt) \int_0^{\infty} \beta_H(t) v_{mH}(dt,\theta_0)+ o_P(1)
\end{align}
where
\begin{align}
 \beta_H(t) = \Gamma_{H}^{-1/2} \frac{\dot{h}_{\theta_0}(t)}{h_{\theta_0}(t)}
\end{align}
denotes the normalized score function. This expression represents the limit in distribution $\vhat_H$ of the process $\vhat_{mH}$, that is, an
orthogonal
projection of the process $v_{H}$ parallel to the normalized score function $\beta_H$. Note that $\beta_H$ by its definition is of unit norm in the space
$\Lcal_2(H)$. Moreover, functions $\beta_H$ and $\mathbb{1}$ are orthogonal. Therefore, from \eqref{em.pr-expand} we have the limit in
distribution of the process $\vhat_{mH}(\phi)$ is
\begin{align}\label{limit-of-em.pr}
 \vhat_H (&\phi) = v_H(\phi) - \langle \phi, \beta_H \rangle_H v_H(\beta_H)\nonumber\\
 &\stackrel{\tiny{\mbox{from} (\ref{limitvH})}}{=} w_H(\phi) - \langle \mathbb{1}, \phi \rangle_H w_H(\mathbb{1}) - \langle \beta_H, \phi \rangle_H
w_H(\beta_H).
\end{align}
This implies that $\vhat_H$ is an orthogonal projection of the function-parametric $H$-Brownian motion $w_H$ parallel to the subspace generated by 2
functions $\set{\mathbb{1}, \beta_H}$. The process $\vhat_H$ depends not only on the true unknown exponent $\theta_0$ but also the score function
$\beta_H$. In the terminology of \cite{Estate2}, the process $\vhat_H$ is called a $\beta_H$-projected H-Brownian motion.

\subsection{Method}\label{method}

Our approach follows the method in Khmaladze \cite{Estate2} for parametric family of distributions. The main idea can be briefly explained as follows:
Under the null hypothesis, the empirical process $\vhat_{mH}(\phi)$ (see \eqref{empro} and \eqref{func.em.pr}) constructed on the tail of an unknown
distribution $F$, from a fixed threshold $x_0$, possesses an ``unspecified'' limit in distribution (see \eqref{limit-of-em.pr}). ``Unspecified'' here
means that the limit depends on some unknown parameter. We will map the
process $\vhat_{mH}$ into another process $\vhat_{mG}$ (see \eqref{trans-em-pr-Gm}) whose limit in distribution is specified. The mapping procedure is
one-to-one to guarantee that no statistical information is lost.

\medskip

We chose $G$ to be the exponential distribution $G(t)= 1-e^{-t}$ simply by preference with some reasons. That is, both $G$ and $H$ are members of
the GPD family (see \eqref{gpd}). The distribution $G$ is the limiting distribution of the GPD$(k,\sigma)$ as $k \to
0$ and scaled by $\sigma=1$. A research by Davison and Smith \cite{Davison} employed Kolmogorov-Smirnov and Anderson-Darling statistics to
test the fit of the GPD to data using the critical value derived from an exponential distribution for the purpose. They discussed that the approach
may be suspect since the exponential distribution is just a member of the whole family though.

\medskip

It is obvious that distributions $G$
and $H$ are equivalent or in other words, mutually absolutely continuous. It is also easy to check that the Fisher
information of $G$ is $\Gamma_G =1$. Put
\[
\ell(t) = \sqrt{\frac{dG}{dH}(t)} = \theta^{-\frac{1}{2}} (1+t)^{\frac{\theta+1}{2}} e^{-\frac{t}{2}},
\]
then this function belongs to $\Lcal_2(H)$. In addition, if $\phi\in\Lcal_2(G)$ then $\ell\phi\in\Lcal_2(H)$ and $\norm{\phi}_G = \norm{\ell\phi}_H.$ Note
that when
$H$ and $G$ are equivalent, it is known that a transformation from a $H$-Brownian motion into a $G$-Brownian motion is straightforward by multiplication.
Namely, $w_H(\ell \phi) = w_G(\phi)$ is a G-Brownian motion in $\Lcal_2(G)$. However, mapping a Brownian
bridge such as $v_H$ or $\vhat_H$ to another Brownian bridge is not that straightforward any more. The fact is that $v_H(\ell\phi)$ depends on both $H$
and $G$ and so does $\vhat_H(\ell\phi)$.

\medskip

Consider a subspace $\Lcalhat$ of $\Lcal_2(H)$ generated by four functions $\set{\mathbb{1}, \beta_H, \ell,\ell\beta_G}$. These functions are of unit norm
in
$\Lcal_2(H)$. More explicitly, the score functions $\beta_H$ and $\beta_G$ are
\begin{align}\label{scorefuns}
\beta_H(t) = 1- \theta \log(1+t),\qquad \qquad \beta_G(t) = 1-t.
\end{align}

\noindent For any function $f$ and $g$ in $\Lcal_2(H)$, define the unitary operator $K_{f,g}$ as
\begin{align*}
K_{f,g} = I - \frac{1}{1-\langle f,g\rangle_H}(g-f) \langle g-f, \cdot \rangle_H
\end{align*}
where $I$ is the identity function. This unitary operator will turn the function $f$ into $g$ and reversely $g$ to $f$ and any function which is orthogonal
to $f$ and $g$ in the subspace $\Lcal_2(H)$ into itself.

\medskip

Step into the method, we first consider the unitary operator
\begin{align}
K_{\mathbb{1},\ell} = I - \frac{1}{1-\langle \mathbb{1}, \ell\rangle_H} (\ell-\mathbb{1})\langle \ell-\mathbb{1}, \cdot \rangle_H.
\end{align}
This operator will map $\ell$ to $\mathbb{1}$ and $\mathbb{1}$ to $\ell.$ Next, consider the image of the function $\ell \beta_G$ via
$K_{\mathbb{1},\ell}$, which is
\begin{align*}
\widetilde{\ell\beta}_G &= \ell \beta_G - \frac{1}{1-\langle \mathbb{1}, \ell\rangle_H}(\ell -\mathbb{1}) \langle \ell-\mathbb{1},\ell\beta_G\rangle_H\\
&= \ell \beta_G - \frac{1}{1-\int_0^{\infty} \ell(s)h(s)ds} (\ell -\mathbb{1})\int_0^{\infty} (\ell(s) -1) \ell(s) \beta_G(s) h(s) ds.
\end{align*}
Then consider the operator $K_{\beta_H,\widetilde{\ell\beta}_G}$ defined as
\begin{align*}
K_{\beta_H,\widetilde{\ell\beta}_G} = I - \frac{1}{1-\langle \beta_H,\widetilde{\ell\beta}_G\rangle_H} (\widetilde{\ell\beta}_G -\beta_H) \langle
\widetilde{\ell\beta}_G - \beta_H,\cdot\rangle_H.
\end{align*}

Set $\Khat = K_{\beta_H,\widetilde{\ell\beta}_G}K_{\mathbb{1},\ell}$, then this unitary operator will map $\ell$ to $\mathbb{1}$ and $\ell \beta_G$ to
$\beta_H$. The non-uniqueness of such unitary operator like $\Khat$ was discussed thoroughly in Khmaladze \cite{Estate2}, Section
3.4. Nevertheless, we believe that this operator $\Khat$ is simple enough for practical purpose, especially with only one parameter.

The main result for testing composite hypothesis was stated as Theorem 7 in \cite{Estate2}, it is essential so we restate it here accordingly to our
notations.

\begin{theorem}\label{theorem} (Restatement of Theorem 7 in \cite{Estate2})

 If $\vhat_H$ is a $\beta_H$-projected $H$-Brownian motion and $G$ is absolutely continuous with respect to $H$, then
 \begin{align}
  \vhat_G(\phi) = \vhat_H(\Khat(\ell\phi)) = \Khat(\vhat_H(\ell\phi))
 \end{align}
is a $\beta_G$-projected G-Brownian motion.
\end{theorem}

\noindent As a consequence, transform the function-parametric tail empirical process $\vhat_{mH}(\phi)$ by $\Khat$, we obtain another process
\begin{align}\label{trans-em-pr-Gm}
 \vhat_{mG}(\phi) = \vhat_{mH}(\Khat(\ell\phi)) = \Khat(\vhat_{mH}(\ell\phi)),
\end{align}
which has $\vhat_G(\phi)$ as a limit in distribution.

\noindent Let
\begin{align*}
\phi_x(t) = 1_{\set{t\le x}}
\end{align*}
be series of indicator functions defined on $\Lcal_2(H)$ depending on $x$. Let $x$ runs from $0$ to $\infty$. From Theorem \ref{theorem}, we know that
the
limit in distribution $\vhat_H(\Khat(\ell\phi_x))$ of the process $\vhat_{mH}(\Khat(\ell\phi_x))$ is a $\beta_G$-projected G-Brownian motion in
$\phi_x$. Hence, any statistic as an appropriate function based on $\vhat_{mH}(\Khat(\ell\phi_x))$ will be asymptotically distribution free. Denote by
$\phitil_x$ the image of $\ell\phi_x$ via the operator $\Khat$, that is,
\begin{align*}
 \phitil_x = \Khat(\ell \phi_x) = K_{\beta_H,\widetilde{\ell\beta}_G}K_{\mathbb{1},\ell} (\ell\phi_x).
\end{align*}
For practical purpose, the form of $\phitil_x$ is
\begin{align}
 \phitil_x &= \ell \phi_x - \frac{1}{1-\langle \mathbb{1},\ell\rangle_H} (\ell-1) \langle \ell-1,\ell\phi_x\rangle_H - \frac{1}{1-\langle
\beta_H,\widetilde{\ell\beta_G}\rangle_H}(\widetilde{\ell\beta_G}-\beta_H)\\\nonumber
 &\times \Big[\langle
\widetilde{\ell\beta_G} -\beta_H,\ell\phi_x\rangle_H - \frac{\langle \ell-\mathbb{1},\ell\phi_x\rangle_H}{1-\langle
\mathbb{1},\ell\rangle_H}\langle\widetilde{\ell\beta_G}-\beta_H,\ell-\mathbb{1}\rangle_H\Big].
\end{align}

\noindent Applying the unitary operator $\Khat$ on the process $\vhat_{mH}(\ell\phi_x)$ we have
\begin{align}\label{trans-em-pr-G}
 \vhat_{mG}(\phi_x) = \Khat(\vhat_{mH}(\ell\phi_x)) &= \vhat_{mH}(\phitil_x) = \int_0^{\infty} \phitil_x(t) \vhat_{mH}(dt)\nonumber\\
 &= \frac{1}{\sqrt{m}} \sum_{i=1}^m \big[\phitil_x(\Ttil_i) - \Esf_{\thetahat_m} \phitil_x(\Ttil_i)\big].
\end{align}
Here $\Esf_{\thetahat_m}$ denotes the expected value with respect to the distribution $H_{\thetahat_m}$. That means,
$\Esf_{\thetahat_m} f(\Ttil) = \int_0^{\infty} f(s) h_{\thetahat_m}(s) ds$ for any integrable function $f$.
As a result of the Theorem \ref{theorem}, the limit in distribution of process $ \vhat_{mH}(\phitil_x)=\Khat(\vhat_{mH}(\ell\phi_x))$ is $\vhat_G(\phi_x)$
- a projected $G$-Brownian motion. We demonstrate in the next Section that any functionals of $\vhat_{mH}(\phitil_x)$ is asymptotically distribution
free.

\section{Simulation results}\label{simulation}

The main purpose of this section is to show the asymptotically distribution free property of the new test statistics based on the transformed empirical
process $\vhat_{mH}(\phitil_x)$. To do
so, we will take some prevalent gof tests as examples, namely the Kolmogorov-Smirnov (KS), the Cramer-von Mises and the Anderson-Darling tests. These
test statistics are some specific functionals of the transformed process $\vhat_{mH}(\phitil_x)$. Namely, an analogue version of the
Kolmogorov-Smirnov test statistic is
\begin{align}
 KS = \max_{x}  \abs{\vhat_{mH}(\phitil_x)}.
\end{align}
The Cramer-von Mises statistics should be of the form
\begin{align}
 \Omega^2 = \int_0^{\infty} \vhat^2_{mH}(\phitil_x) dG(x),
\end{align}
and its weighted version called Anderson-Darling statistic is
\begin{align}
 A^2 = \int_0^{\infty} \frac{\vhat^2_{mH}(\phitil_x)}{G(x)(1-G(x))} dG(x).
\end{align}

\begin{figure}
	 \begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{KS-x0-3-1000ss}
		\subcaption{\it Threshold $x_0=3,$ sample size $n=1000$}
		\label{fig:x0-3-ss-1000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{KS-x0-5-5000ss}
		\subcaption{\it Threshold $x_0=5,$ sample size $n=5000$}
		\label{fig:x0-5-ss-5000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{KS-x0-10-8000ss}
		\subcaption{\it Threshold $x_0=10$, sample size $n=8000$}
		\label{fig:x0-10-ss-8000}
	\end{subfigure}
	\caption{\it Distribution of the KS test statistics. Solid line: Pareto distribution with $\theta_0=3$; Dashed line: Cauchy distribution}
	\label{fig:KS-test}
\end{figure}
\begin{figure}
	 \begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-x0-3-1000ss}
		\subcaption{\it Threshold $x_0=3,$ sample size $n=1000$}
		\label{fig:omega-x0-3-ss-1000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-x0-5-5000ss}
		\subcaption{\it Threshold $x_0=5,$ sample size $n=5000$}
		\label{fig:omega-x0-5-ss-5000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-x0-10-8000ss}
		\subcaption{\it Threshold $x_0=10$, sample size $n=5000$}
		\label{fig:omega-x0-10-ss-5000}
	\end{subfigure}
	\caption{\it Distribution of the $\Omega^2$ test statistics. Solid line: Pareto distribution with $\theta_0=2$; Dashed line: Cauchy distribution}
	\label{fig:omega-test}
\end{figure}
\begin{figure}
	 \begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-weighted-x0-3-1000ss}
		\subcaption{\it Threshold $x_0=3,$ sample size $n=1000$}
		\label{fig:omega-weighted-x0-3-ss-1000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-weighted-x0-5-5000ss}
		\subcaption{\it Threshold $x_0=5,$ sample size $n=5000$}
		\label{fig:omega-weighted-x0-5-ss-5000}
	\end{subfigure}
	\hspace{.03\linewidth}
	\begin{subfigure}{.3\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{Omega-weighted-x0-10-8000ss}
		\subcaption{\it Threshold $x_0=10$, sample size $n=5000$}
		\label{fig:omega-weighted-x0-10-ss-5000}
	\end{subfigure}
	\caption{\it Distribution of the $A^2$ test statistics. Solid line: Pareto distribution with $\theta_0=0.5$; Dashed line: Cauchy
distribution}
	\label{fig:anderson-darling-test}
\end{figure}

\newpage
As shown in Section \ref{method}, the limit in distribution of the process $\vhat_{mH}(\phitil_x)$ is a projected Brownian motion in time
$G(x)$. Hence, theoretically these test statistics are asymptotically distribution free. Especially, it makes sense that we integrate with respect to
nothing else but $G(x)$ in order to get the $\Omega^2$ and $A^2$ tests. In principle, $x$ runs from $0$ to infinity. However, we choose $x\in
\set{0.1,0.2, \cdots, 7.9,8}$ since $G(8) = 0.9997 \approx 1$ hence we practically just need $x_{\max}=8$ as the maximum value for $x$. The $\Omega^2$ and
$A^2$ tests can be approximated as follows:
\begin{align}
 \Omega^2 &\approx 0.1 \times \sum_{x=0.1}^{8}  \vhat_{mH}(\phitil_x) e^{-x},\nonumber\\
 A^2 &\approx 0.1 \times  \sum_{x=0.1}^{8}  \frac{\vhat_{mH}(\phitil_x)}{1-e^{-x}}.\nonumber
\end{align}

To create the curves of the cumulative distributions of the new tests, we choose Pareto and Cauchy distributions as the underlying distributions. It is
known that the tail of the Cauchy distribution is regularly varying with the exponent $\theta_0 =1$. For the Pareto distributions, the exponent is
arbitrarily selected and positive. We did choose some $\theta_0$ ranging from $0.5$ to $10$. Sample size $n$ for simulation needs to be large to guarantee
that the tail is sufficiently big, namely, $m$ is not less than $40$. Hence, we chose $n$ at least equals $1000$. We also chose different threshold $x_0$,
namely $x_0 =3,5,10$. Usually, each curve is produced by $5000$ iterations of simulation.

Figures \ref{fig:KS-test}, \ref{fig:omega-test} and \ref{fig:anderson-darling-test} show the plots of the cumulative functions for $KS, \Omega^2, A^2$
test statistics respectively with different choices of the threshold $x_0$, different sample sizes and different underlying distributions. As we can see
in these figures, the two curves of the cumulative
distribution functions in each plot are not distinguishable, which practically demonstrate the asymptotically distribution free property of the new test
statistics.
Moreover, as we notice, for different $x_0$, the difference between these curves is also very minor.

Regarding the time of the procedure, it took approximately $1$ hour to create the cumulative distribution functions with $x_0 =5$ and sample size $n=5000$
by $5000$ iterations for two different original distributions $F$. Therefore, we believe that the method is easy to implement and also quick.

\section{Discussion}

It is possible to extend this approach to testing multidimensional regularly varying tail distributions. As long as we are able to find a transparent expression for the class of regularly varying tail distributions with an explicit parametric family form, we are able to transform the whole family into a specific distribution, from that we can build up a whole class of asymptotic distribution-free test statistics.

\bigskip

\noindent {\bf Acknowledgements}

The author is greatly indebted to Professor Estate Khmaladze for his invaluable suggestion for approaching this
problem
and many useful discussions during this research. This research was conducted during the author's doctorate of philosophy degree, so the author would also like to thank Victoria University of Wellington for providing the doctoral scholarship to conduct this research.


TITLE: One periodic binomial model
QUESTION [1 upvotes]: I need to look into a one-period Binomial model $(B_t, S_t)$ with interest rate $r = 0.1$ , $S_0 = 100$ and
$$
S_t= 120 \, \text{with probability}\, 0.5
$$
$$
S_t= 60\, \text{with probability}\, 0.5
$$
a) I need to show this market has no arbitrage.
b) Find the replicating portfolio for the contingent claim $X = S^2_T$.
c) Find the price of $X$ that leads to no arbitrage.
d) Find the risk-neutral measure for this market.
Now I guess that in a) I need to show that the state-price vector is positive, but how do I do that? For the rest I have no idea and would welcome your help!

REPLY [3 votes]: To get you started, let $\omega_1$ and $\omega_2$ denote the two market states at terminal time $T$.
There are two securities, a bond with initial price $B_0 = 1$ and terminal prices $B_T(\omega_1) = B_T(\omega_2) = 1.1$ and,  a stock with initial price $S_0 = 100$ and terminal prices $S_T(\omega_1) = 120$ and $S_T(\omega_2) = 60$.
For (a), the state-price vector $(q_1,q_2)'$ is the solution of
$$\pmatrix{B_T(\omega_1) & B_T(\omega_2) \\S_T(\omega_1) & S_T(\omega_2)}\pmatrix{q_1\\q_2 }=\pmatrix{B_0\\S_0 },$$
which reduces to
$$1.1 q_1 + 1.1 q_2 = 1\\ 120 q_1 + 60 q_2 = 100$$
The existence of the solution $q_1 \approx 0.7576, \, q_2 \approx 0.1515$ implies no arbitrage.
For (b), let $\alpha$ and $\beta$ be the holdings of stock and bond in the portfolio that replicates the contingent claim $X = S_T^2$. The value of the portfolio must be identical to the value of the claim in all terminal states, implying
$$\alpha S_T(\omega_1) + \beta B_T(\omega_1) = S_T(\omega_1)^2\\ \alpha S_T(\omega_2) + \beta B_T(\omega_2) = S_T(\omega_2)^2$$
This reduces to
$$120 \alpha +1.1 \beta = (120)^2\\ 60 \alpha + 1.1 \beta = (60)^2$$
with the solution $\alpha = 180$ and $\beta \approx -6545.45$.
From here it should be obvious how to find the fair value of the claim for part (c) as the initial value of the replicating portfolio. I'll let you think about how to find the risk-neutral measure.