\section{Introduction}

One of the central debates in language acquisition is whether the syntactic structures of natural language are genetically specified or learned through experience. A key tool in this debate has been the use of formal mathematical analysis in order to determine what learners could or could not logically conclude from the type of data they observe. Perhaps the most famous is \emph{Gold's Theorem} \cite{gold1967language}, which holds that observers of positive examples of even a \emph{regular language} \cite{hopcroft1979introduction} could not identify the set of allowable (``grammatical'') strings with certainty. In linguistics and cognitive science, Gold's proof  was taken to mean that innate constraints must be present for learning to succeed. The theorem gave rise to rich formal theories of what could be learned under similar assumptions and correspondingly what must be innate \cite<e.g.>{wexler1983formal}, although its relevance for cognitive science is not universally accepted \cite{johnson2004gold}. Indeed, more recent work has undermined the view that learners must come with substantial built-in knowledge about the specific structures of language. \citeA{chater2007ideal} describe a theoretical learning framework without tight constraints that could identify the right language still from only positive data. \citeA{chater2007ideal}'s make two key differences from Gold-style analyes: first, \citeA{chater2007ideal} assume that sentences are \emph{sampled} from the parent's generative model, not a worst-case analysis. Second, the model considers \emph{all possible computations} as hypotheses, showing that the space need not be constrained innately for an idealized learner. The core underlying idea of this work is that learners attempt to find simple descriptions of the data they observe, much like work in Minimum Description Length \cite{grunwald2007minimum} and artificial intelligence more generally \cite{hutter2005universal}. This perspective has been further developed to correctly predict the difficulty of constructions \cite{hsu2010logical,hsu2011probabilistic}.

The present paper applies this perspective to learning the structures supporting natural language that was first presented in Chomsky's ``Three models for the description of language'' \cite{chomsky1956three, chomsky1957syntactic, chomsky1959certain}. Chomsky noted that many dependencies in natural language could be captured abstractly with fundamentally different kinds of computational devices. Some devices require a finite amount of memory (e.g. finite-state machines), some require a stack (e.g. phrase structure rules), and some require even more powerful computational systems (e.g. context-sensitive grammars). These relationships hold true for even very simple versions of natural language structures, including sets of strings---called \emph{formal languages}---over simple alphabets. For example, the simple recursion/iteration allowed by English adjectives (``The beautiful well-dressed young lady'') could be captured by the formal language $\lbrace a^n : n=1,2,3,\ldots \rbrace$ indicating that any number of adjectives ($a$s) could be put together into a valid (sub)string. The language $(ab)^n$ could be seen in, for instance, embedding in lists of determiner-noun pairs (``Bring me two cars, three accordions, and six pickles.'') or in subject-verb pairs of sentences in disourse (``John laughed. Sally cried. Mindy marveled.''). The formal language $\lbrace a^n b^n : n =1,2,\ldots\rbrace = \lbrace ab, aabb, aaabb, \ldots \rbrace$ might characterize the key dependencies in english \emph{if-then} relationships where every ``if'' (an ``a'') must be followed by a ``then'' (a ``b''), as in ``\textbf{If} \textbf{If} Mary cried \textbf{then} John was sad \textbf{then} John cares about Mary.''

In this study, we consider a learning model that, like the one of \citeA{chater2007ideal}, builds in a space of \emph{all computations}. In contrast to Chomsky's \emph{Universal Grammar} (UG), specifying such space, one only need to build in a few operations, a minimal Turing complete set. Results from computer science, complexity theory, and logic have shown that it takes very little formal mechanism to allow Turing-completeness. By building in only such a minimal mechanism---here, lambda calculus \cite{church1936unsolvable, hindley1986introduction}---we are able to construct a theory whose innate assumptions are explicit and extremely minimal. Out of a large, implicitly defined hypothesis space, learners are able to construct \emph{new computations} which correspond to different types of formal language systems. Our model shows that given a single such system, learners could construct or genuinely discover cognitive representations corresponding to finite, regular, context-free, even context-sensitive languages and beyond.

Previous study conducted by \citeA{perfors2011learnability} follows a similar paradigm, where they present a learner who does Bayesian inference over different kinds of grammatical systems, including finite, regular, and context-free grammars. They show that with very little positive evidence in the form of part-of-speech sequences from CHILDES \cite{macwhinney2000childes} learners could infer that natural language requires a context free grammar. However, the key difference between the present work and \cite{perfors2011learnability} is that we do not build in just a handful of possible alternatives; instead, we define an implicit space of hypotheses that is more concise or parsimonious to build in a larger, unconstrained space of hypotheses, allowing the description length of what must be assumed by the model to be very short. As the former has often been criticized to have ``built in'' \emph{more} than alternatives such as UG.

In the next section we describe a computational model that can acquire these kinds of fundamentally different computational structures. We then show that it can learn a variety of formal languages from positive evidence only, and that in doing so it is unproblematic to acquire languages with infinite cardinality, even given finite data. We then show that the model can combine several different computational sub-structures and acquire a grammar of simplified English, including non-context-free language structures described in prior literature. In the last section, we show that the model is able to captures several experimental findings from artificial language learning.

\section{A computational model}

The starting assumption for our model is~\citeA{fodor}'s \emph{Language of Thought} (LOT), the hypothesis that the structures of mental operations are structured and compositional, or language-like,. Following now a growing body of work in Bayesian LOT models \cite{goodman2008rational, piantadosi2012bootstrapping, probmods, piantadosi2011learning, piantadosi2016logical, piantadosi2016four, yildirim2013transfer}, we define a set of operations that compose and have learners consider all compositions of these operations as possible hypotheses. These meanings can be interpreted analogously to short computer programs which ``compute'' strings of part of speech sequences, although the input and output could be different in other domains. Importantly, the specification of these programs is \emph{domain general} meaning that essentially the same setup can be used to explain acquisition in other domains like counting, kinship, cross-modal transfer, and concept learning. In all cases, the task of the learner is to determine which compositions of primitives are likely to be correct, given the observed data.

\subsection{Hypothesis space}

\begin{table}[t!]
\centering
\caption{Primitive functions used in LOT.}
\label{prim}
\begin{tabular}{lllll}
\hline
\multicolumn{2}{l}{\bf Functions on strings} \\
 $pair(X, Y)$  &   concatenate lists $X$ and $Y$ &   \\
 $first(X)$  &  return the first element of list $X$  &    \\
 $rest(X)$  &  return the list $X$ without the first element& \\
   &  & \\
\multicolumn{2}{l}{\bf Logical functions} \\
 $flip(p)$  &  randomly return true (probability $p$) or false (probability $1-p$) & \\
 $empty(X)$  &  return true if string $X$ is empty, otherwise false & \\
 $if(B,X,Y)$ & return $X$ if $B$ else return $Y$ \\
   &  & \\
\multicolumn{2}{l}{\bf Function calls} \\
 $F_i(s)$  &  call the expression $F_i$ with argument $s$& \\ \hline
\end{tabular}
\end{table}

The primitives that we assume are compoased are motivated by minimalist, functional programming languages like \emph{scheme} \cite{sicp}, which also forms the basis of other LOT work \cite<e.g.>{probmods}. There are 3 classes of primitive functions, shown in Table.\ref{prim}. The first 3 are list operators: \emph{pair}, \emph{first} and \emph{rest}, which respectively join two elements together (somewhat like \citeA{chomsky1995minimalist}'s \emph{merge}), take the first element of a pair, or take the second element. The second class of operators are \emph{flip} and \emph{empty}, which return boolean values. \emph{flip} is notable in that it provides a stochastic element to the grammar, allowing nondeterminism in the generation of structures. \emph{if} allows for conditional expressions. Finally we allowa function to call another function $g_i$---potentially itself---via $function(list)$. We note the similarity between this grammar and others used in the conceptual modeling with a LOT in entirely different domains \cite{piantadosi2012bootstrapping}. While these operations look very simple, they actually permit a surprising range of computations to be expressed.

\begin{table}[tb!]
\centering
\caption{The grammar used to specify valid compositions of primitives.}
\label{rule}
\begin{tabular}{lll}
\hline
$START \to \lambda x. list$              & $atom \to a,b,c,\ldots$  & \\
$list \to pairs(list, list)$  & $list \to function(list)$ & \\
$list \to rest(list)$         & $list \to if(bool, list, list)$   & \\
$list \to first(list)$        & $bool \to flip(probability)$ & \\
$list \to \emptyset$          & $probability  \to 0.5, 0.6, \ldots 0.9$ & \\
$list \to x$                  & $bool \to empty(list)$ &\\
$list \to atom$ & & \\
\hline
\end{tabular}
\end{table}

Compositions of these primitives and grammar define a set of functions which act as hypotheses for learners. In our analogy to the Library of Babel, these primitives are the characters in the books, and a book (a composition of primitives) is a single hypothesis. The goal of the learner is to walk in the library (hypothesis space) and find the best book (hypothesis) that fits the observed data. For instance, a possible expression built out of these primitives is,
\begin{center}
\begin{tabular}{c}
\begin{lstlisting}[mathescape]
F1(x) = pair( a, if(flip(0.7), $\emptyset$, F1($\emptyset$)))
\end{lstlisting}
\end{tabular}
\end{center}
Here we have named this function ``F1'' and have allowed it to recursively call itself in its own definition via \emph{call}. When this expression evaluates, it concatenates (\emph{pair}s) an ``a'' with either an emptry string ($\emptyset$) or the outcome of calling itself with an empty string as an argument, \lstinline[mathescape]{call(F1, $\emptyset$)}. In terms of a set of strings, this function represents the formal language $\lbrace a, aa, aaa, aaaa, \ldots \rbrace = \lbrace a^n : n = 1,2,3,\ldots \rbrace$. However, this function also gives a \emph{probability distribution} over strings since some (e.g. $a$) are much more likely to be generated than others (e.g. $aaaa$) according to the probability of each \emph{flip}. In particular, this probability distribution this assigns is geometric, $P(a^n) = 0.7^{n-1}\cdot 0.3$. Note, however, that these operations permit a vast array of different types of languages and computations to be specified.

\subsection{Factorized hypothesis}

One important feature of nearly all working computational systems is that they have a mechanisms to handle complexity, in large part through re-use of subcomputations (For work on learning principled subcomputations to re-use, see \cite{o2015productivity}). One simple way to manage the complexity of hypotheses and efficiently structure search is to allow learners to incrementally build on hypotheses, extending partially working theories in order to cover new data. Here, we introduce a \emph{factorized} hypothesis consisting of a sequence of functions, each of which may be defined in terms of prior functions. This is best illustrated with an example:
\begin{center}
\begin{tabular}{lll}
F1(x) & = &
\begin{lstlisting}[mathescape]
pair(a, if(flip(0.7), $\emptyset$, F1($\emptyset$)))}}
\end{lstlisting} \\
F2(x) & = &
\begin{lstlisting}[mathescape]
if(empty(x), $\emptyset$, pair(first(x), pair(F2(rest(x)), b)))}
\end{lstlisting} \\
F3(x) & = &
\begin{lstlisting}[mathescape]
F2(F1($\emptyset$))
\end{lstlisting} \\
\end{tabular}
\end{center}
The first fucntion, F1, is simply the $a^n$ language shown above. The second function, F2, is more interesting: it is a function that takes an argument $x$, and recursively concatenates the first element of $x$ with a recursive call to itself, followed by a 'b'.  Because this is recursive on a shortening list (\emph{rest(x)}), one 'b' will be added for each element of $x$. For instance, if we called F2 on the string 'xyz', we would be left with the string 'xyzbbb'. Finally, F3 puts F1 and F2 together, calling F2 with F1 as an argument. Thus, the strings $a^n$ generated by F1 will be passed to F2, which will attach a single 'b' to for each 'a', giving rise to strings of the form $a^n b^n$.\footnote{Note this is just an example, as there are other, simpler ways to generate $a^nb^n$.} In this way, a complex computation can be broken down into pieces that are more manable, each of which might be more easily learned on their own. In the domain of language or grammar learning, such a division of labor is especially elegant since many valid strings have valid substrings, so a model that captures part of the data (the substrings) in F1 can be extended to a model that captures the longer strings as well in F2. Note that in permitting passing arguments to functions, this formamlism permits important types of variable abstraction.

In the implementation, we allow up to 10 different functions to be called in this way. Note that each function can only call itself and prior functions (e.g. F3 cannot call F4, but F4 can call F3 and F4). Also, it is assumed by convention that the output of the function is the final hypothesis called with the empty string as input.

\subsection{The probabilistic model}

Next, we turn to the specification of a probabilistic model over this space that supports tractable learning and inference. Following prior work using the Bayesian LOT, we convert the grammar of operations to a Probabilistic Context-Free Grammar (PCFG) \cite{manning1999foundations} that specifies the prior probability of any hypothesis. The use of a PCFG automatically penalizes complex hypotheses, fitting the idea that \emph{simplicity} is a core inductive pressure \cite{chater2003simplicity,feldman2000minimization}.

The goal of a learner is to observe some strings $D$ and infer a posterior on hypotheses, $P(H|D)$. Following Bayes rule, $P(H|D) \propto P(H)P(D|H)$, where $P(H)$ is given by the PCFG times a $2^{-K}$ penalty on having $K$ factorized components (F1, F2, $\ldots$, FK), and $P(D|H)$ is a multinomial likelihood specifying how likely the observed strings are to be generated by $H$. This likelihood include an outlier likelihood penalty of log value $-1000$ for strings in the data that are assigned zero probability by $H$. To manage the uncomputability inherent in this space, hypotheses that do not ``halt'' after a set number of recursive calls---e.g. those that are likely to be in an infinite loop---are assigned a zero prior.

In inference, most of the posterior probability mass will be held by just a few concise and well-fitting hypotheses. In our implementation (the code can be found on GitHub \cite{piantadosi2014lotlib}), we run 12 Monte-Carlo chains for 50000 steps of the tree-regeneration MCMC method presented in \citeA{goodman2008rational} for varying amounts of data. his algorithm essentially takes a random walk around hypotheses $H$ that is biased towards high probability regions. Over time, this will converge on hypotheses with high $P(H|D)$ and indeed can be shown to correctly be sampling from this posterior distribution. We collect 100 of the highest-posterior expressions found in each chain. These samples are taken as a fixed, finite hypothesis space for the purposes of computing summary statistics and plots. Code is available online in LOTlib \cite{piantadosi2014lotlib}. In each model run, we add to the grammar the required terminal symbols specific to each data set.

\section{Learning simple formal languages}

We first study the model's ability to learn several simple formal languages that are motivated by the structures present in natural language. Since Chomsky (1956), numerous authors have attempted to pinpoint what type of computational system supports natural language, with much of the debate focusing on what constructions evidence computational requirements more than context free languages. This has motivated attempts to map natural language structures onto very simple formal languages whose computational properties can be fully characterized.

Specifically, we examine the languages $a^n$, $a^nb^n$, $a^nb^{2n}$, $(ab)^n$, $a^nb^nc^n$, and the Dyck language (valid sequences of parentheses). Here, $a^n$ and $(ab)^n$ are regular languages and can be accepted by finite state automation;  $a^nb^n$, $a^nb^{2n}$ and the $Dyck$ language are context free languages, and are recognized by pushdown automation; $a^nb^nc^n$ is a context sensitive language that requires more than a context free grammar. Critically, many of these languages have analogs in natural language synta; and the Dyck language specifies the valid ways in which a sentence could be parsed.

\subsection{Results}

\begin{figure}
	\centering

	\begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/an.pdf}
    \caption{}
    \label{g_an_wf}
	\end{subfigure} \hfill
    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/abn.pdf}
    \caption{}
    \label{g_abn_wf}
	\end{subfigure}\hfill
    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/anbn.pdf}
    \caption{}
    \label{g_anbn_wf}
	\end{subfigure}

    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/anb2n.pdf}
    \caption{}
    \label{g_anb2n_wf}
	\end{subfigure} \hfill
    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/dyck.pdf}
    \caption{}
    \label{g_dyck_wf}
	\end{subfigure}\hfill
    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/anbncn.pdf}
    \caption{}
    \label{g_anbncn_wf}
	\end{subfigure}

    \caption{The graph shows the weighted F-score against the amount of observed data: (a) $a^n$; (b) $(ab)^n$; (c) $a^nb^n$; (d) $a^nb^{2n}$; (e) $Dyck$; (f) $a^nb^nc^n$. The F-score curve for Maximum a Posteriori (MAP) hypothesis is also shown. }
    \label{g_wf}
\end{figure}

Fig.\ref{g_wf} shows the posterior-weighted precision, recall and F-scores of the hypotheses given different sizes of data from these formal languages. Here, the data consist of observations of positive examples of strings (e.g. from $a^nb^n$ we might observe $\lbrace ab, ab, aabb, ab, aaabbb \rbrace$) and is provided to \emph{the same} learning model. The model is able to rapidly acquire the correct kind of computational system for each, often taking only 5-10 samples of strings before the model accurately learns the right system. The implication of this result is that ideal learner with domain-general constraints is able to learn the unobserved grammatical structures and computations that underlie a target set of strings rapidly, and from very little data. To illustrate, \ref{Tab:example_hyp} shows example $perfect$ hypotheses for each type of language.

\begin{table}[htb]
\centering
\caption{Example $perfect$ hypotheses for 6 classes for formal languages.}
\label{Tab:example_hyp}
\begin{tabular}{lllll}
\hline
$a^n$ & $a^nb^n$\\
\begin{lstlisting}[mathescape]
def $F_1$():
    return $pair$($\mb{a}$,$if$($flip()$,$recurse$(),$\emptyset$))
\end{lstlisting}&
\begin{lstlisting}[mathescape]
def $F_1$():
    def $F_x$():
        return $if$($flip()$,$\emptyset$,$F_1$())
    return $pair$($\mb{a}$,$pair$($F_x$(),$\mb{b}$))
\end{lstlisting} \\

\hline
$(ab)^n$ & $a^nb^{2n}$\\
\begin{lstlisting}[mathescape]
def $F_1$():
    def $F_x$():
        return $if$($flip()$,$\emptyset$,$F_1$())
    return $pair$($pair$($\mb{a}$,$\mb{b}$),$F_x$())
\end{lstlisting}&
\begin{lstlisting}[mathescape]
def $F_1$():
    def $F_x$():
        return $if$($flip()$,$\emptyset$,$F_1$())
    return $pair$($\mb{a}$,$pair$($F_x$(),$pair$($\mb{b}$,$\mb{b}$))
\end{lstlisting}\\

\hline
$Dyck$ & $a^nb^nc^n$\\
\begin{lstlisting}[mathescape]
def $F_1$():
    def $F_x$():
        return $if$($flip()$,$\emptyset$,$F_1$())
    return $pair$($\mb{a}$,$pair$($F_x$(),$pair$($\mb{b}$,$F_x$())))
\end{lstlisting}&
\begin{lstlisting}[mathescape]
def $F$():
    def $F_1$():
        return $pair$($\mb{b}$,$\mb{c}$)
    def $F_2$($g_1$):
        def $F_x$():
            return $pair$($pair$($\mb{b}$,$g_1$()),$\mb{c}$)
        return $pair$($\mb{a}$,$if$($flip$(),$g_1$(),$recurse$($F_x$)))
    return $F_2$($F_1$).
\end{lstlisting}\\
\hline
\end{tabular}
\end{table}

\subsection{Learning that a language is infinite}

One of the most striking properties of natural language is that there are infinitely many sentences of valid, grammatical English \cite<though see>{pullum2010recursion}. This fact might be considered to be a core aspect of our innate linguistic endowment, perhaps a consequence of a generative capacity that forms UG. On the other hand, once we consider learners who are computationally sophisticated, it may be possible to \emph{learn} that the generative system supporting language is infinite. After all, there are aspects of language which are \emph{not} infinite (e.g. the set of function words) and there have even been argued to be languages with finite cardinality \cite{futrellcorpus}. Interestingly, from a learning point of view, it is not entirely obvious what finite data might convince learners that their language permits infinitely many strings.

\begin{figure}
    \centering
	\includegraphics[width=0.5\textwidth]{figs/finite_vs_infinite.pdf}
    \caption{The sum of posterior probabilities of infinite hypotheses given finite and infinite target data.}
    \label{g_f_an}
\end{figure}

To investigate this, we generated data from two target hypothesis languages: (a) $a^n$ with $n \leq 3$ or (b) $a^n$ with $n=1,2,3,4,\ldots$. In both cases, we assume $P(a^n) \propto (\frac{1}{2})^n$. Fig.\ref{g_f_an} shows under each type of data the posterior probability of a hypothesis that generates an infinite language. The type of hypothsis the model learns in the infinite case looks like the example above for $a^n$; in the finite case, the model could learn a hypothesis like
\begin{center}
\begin{tabular}{lll}
\begin{lstlisting}[mathescape]
F1(x) = pair(a, pair(if(flip(0.5), $\emptyset$, a), if(flip(0.5), $\emptyset$, a)))
\end{lstlisting}
\end{tabular}
\end{center}
What should be clear is that as finite languages become more complex or have greater cardinality, they become more complex. In the absence of overwhelming data, learners will prefer a more concise hypothesis and this will tend to bias them towards hypotheses that generate infinite sets. It takes additional data to convince a simplicity-driven model that long strings ($n\geq 4$) are not permitted. This shows that languages with infinite cardinality may result from core infernetial properties of a computationally sophisticated learner \cite<see also>{piantadosi2017infinitely}.

\section{Learning a simplified English phrase structure grammar}

Of course, the problem faced in natural syntax is not just a series of unrelated formal languages, but data that combines many forms of computation together. We next test the learning model with a grammar of a simple subset of English:
\begin{center}
\begin{tabular}{lcl}
$S$ & $\overset{4}{\to}$ & $NP\; VP$ \\
$NP$ &$\overset{2}{\to}$ & $n$ \\
$NP$ &$\overset{1}{\to}$ & $d\; n$ \\
$NP$ &$\overset{1}{\to}$ & $d\; AP\; n$ \\
$AP$ &$\overset{3}{\to}$ & $a$ \\
$AP$ &$\overset{1}{\to}$ & $a\; AP$ \\
$VP$ &$\overset{2}{\to}$ & $v$ \\
$VP$ &$\overset{1}{\to}$ & $v\; NP$ \\
$VP$ &$\overset{1}{\to}$ & $v\; that\; S$ \\
\end{tabular}
\end{center}
Here, the number over each arrow gives the unnormalized probability of each expansion (so $AP$ is three times more likely to expand to ``a'' than to ``a AP''). This target grammar includes many interesting syntactic structures of English, including multiple expansions of a nonterminal type (e.g. NP), tail recursion (in AP), transitive and intransitive verbs (in VP), and sentential embedding (in VP).

Importantly, although we have written this as a phrase structure grammar, as above, the learning model is \emph{not} told that the strings it observes come from this grammar, or indeed from any kind of phrase structure at all. The strings that are provided have no parsing information or structure: for instance, one string the model might observe is 'if d a n v then n v that d n v d n', which corresponds to the part-of-speech sequence for a sentence like ``If an angry man smiles then Jane believes that the instrument is an accordion.'' The challenge for the model is to take these kinds of part-of-speech sequences and infer the computational process out of all compositions of primitives.

\begin{figure}[htb]
    \centering
    \hspace{-1cm}
    \begin{subfigure}[]{0.55\columnwidth}
            \includegraphics[width=\textwidth]{figs/simple3.pdf}
            \caption{}
            \label{simple_curve}
    \end{subfigure}
    \hspace{-1.5cm}
    \begin{subfigure}[]{0.5\columnwidth}
            \begin{tabular}[b]{cc}
        	\begin{lstlisting}[mathescape,basicstyle=\ttfamily\fontsize{7}{6}\selectfont]
        def $F$():
            def $F_1$():
        	if $flip$():
    		    return $pair$($\mb{a}$,$recurse$())
        	else:
    	 	    return $pair$($\mb{a}$,$\mb{n}$)

            def $F_2$($g_1$):
                if $flip$():
                    return $\emptyset$
                else:
            	    return $pair$($\mb{d}$,$rest$($g_1$()))

            def $F_3$($g_1$,$g_2$):
            	return $rest$($g_1$())

            def $F_4$($g_1$,$g_2$,$g_3$):
            	$x_0$ = $pair$($pair$($\mb{d}$,$g_3$($g_1$,$g_2$)),$\mb{v}$)
            	$x_1$ = $pair$($x_0$,$\mb{t}$)

                if $empty$($g_2$($g_1$)):
            	    return $pair$($x_0$,$g_2$($g_1$))
                else:
            	    return $pair$($x_1$,$recurse$($g_1$,$g_2$,$g_3$))

            return $F_4$($F_1$,$F_2$,$F_3$)
        	\end{lstlisting}
        	\end{tabular}
        \vspace{.1cm}
    	\caption{}
        \label{simple_example}
    \end{subfigure}
    \caption{(a)The graph shows the weighted F-score of hypotheses learning simple English case given different sizes of data; (b) Example hypothesis for simple English.}
    \label{simple}
\end{figure}

The learning curves are shown in Figure.~\ref{simple_curve} which gives the posterior probability (y-axis) of various types of grammars as a function of amount of data (x-axis). The model eventually settles on a hypotheses Fig.\ref{simple_example}. Importantly, the components of this grammar--recurive, tail-recursive, random-choices of expansions, etc. are not ``built in'' as necessary mechanisms for learners. Instead, the parts of this grammar which mimic a PCFG are constructed by composing function calls with the assumed primitives.

\section{Connections to human learning}

This computational model is interesting not only as a theoretical demonstration of the power of statistical learning, but also because its underlying principles can be directly connected to human behavior. In this section, we consider two prior experiments on human artificial grammar learning (AGL) and show that the model exhibits similar patterns as humans.

\subsection{Structure of the input}

A wealth research has sought to understand human of learning center-embedded recursion structure~\cite{fitch2004computational, lai1, lai2}. This structure, in our context, can be characterized as an ``$a^nb^n$'' language, and an example in natural language can be ``the car that I drove yesterday broke''. It is of particular interest because it is shown to be sufficiently complex and is a crucial feature in human language faculty~\cite{hauser}. \citeauthor{lai2} studied whether an unequal distribution of the training exemplars can influences the process of learning a center-embedded grammar. They argued that a biased distribution, where short exemplars are more likely to be presented to the learners than longer ones, can facilitate acquisition. In Lai et al.'s experiments, participants were asked to learn an artificial language generated with center-embedded recursive grammar which resembles the $a^nb^n$ language. They were exposed to 12 blocks of training data set and completed a test at the end of each block. Participants were divided into two groups. Training data for the first group was randomly sampled from the grammar (Random Case), where short and long exemplars has the same amount; the one for the other group has a biased distribution over exemplars. The amount of exemplars rises exponentially as they get shorter, e.g., the shortest exemplars are 2 times more than the second shortest and 4 times than the third shortest, and so on. These settings are illustrated in Fig.\ref{skewed_g} and Fig.\ref{random_g}. Their results showed an improvement in learning as measured by the rate of correctly recognizing the target language in test phrase: when given the data with biased distribution, participants performance increased by 30\%.

\begin{figure}[tb]
	\centering

	\begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/sf_data.pdf}
    \caption{}
    \label{skewed_g}
	\end{subfigure} \hfill
    \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/si_data.pdf}
    \caption{}
    \label{staged_g}
	\end{subfigure} \hfill
     \begin{subfigure}[]{0.3\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/rc_data.pdf}
    \caption{}
    \label{random_g}
    \end{subfigure}

    \caption{The compositions of data blocks in experimental conditions from Lai et al (2011, 2013): (a) skewed frequency, (b) staged input, and (c) random.}
    \label{fig:data_dist}
\end{figure}

\citeA{lai1} studied another factor that can facilitate the learning, i.e. staged input. They argued that AGL can be facilitated with an explicit control of exemplar difficulties in several learning stages. The experiments setup was similar to the previous one, except that learners in the test group now experience 3 stages of training, each containing 4 blocks. For each stage, the maximum length (or complexity) of exemplars was limited, and this limit increased as participants proceed to higher stage. This setting is illustrated in Fig.\ref{staged_g}. Experimental results showed improvement in performance by 40\% when learners were exposed to staged training.

We constructed three datasets of the same size (144 items) and language ($a^n b^n$) as that used in these experiments. We matched the input data composition to these studies. The only difference was the actual appearance of exemplars to be presented, which for simplicity we mapped to ``a'' and ``b'' instead of syllables.

\begin{figure}[htb]
	\centering
	\includegraphics[width=.4\textwidth]{figs/si_sf_un.pdf}
    \caption{The average F-score of MAP hypothesis given different types of data (skewed frequency, staged input and uniform input) as a function of the total MCMC steps taken.}
    \label{fig:si_sf_un}
\end{figure}

Figure Fig.\ref{fig:si_sf_un} shows the learning curves for each of these kinds of input. These results show that models with skewed frequency and staged input performs better than the random case---reaching a higher F-score faster---agreeing qualitatively with the human results. The computational explanation for this is also intuitive: data with a uniform distribution requires an equal chance to produce strings with every different level of complexity, which is hard to realize in our model: the main device for producing complex strings that follows a certain type of grammar is stochastic recursion, and that the chance for generating complex strings decreases exponentially is almost an innate property of this mechanism.

\subsection{Nonadjacent dependency}

The third case we consider is that of nonadjacent dependencies, one of the most characteristic features of human languages. Examples can be found in the dependencies between auxiliaries and inflectional morphemes (e.g. ``\ul{is} runn\ul{ing}, \ul{has} eat\ul{en}''), as well as dependencies involving number agreement (``The rock\ul{s} on the bluff \ul{are} jagged''); nested nonadjacent dependencies can be found in, for instance, the if-then constructions (above) \cite{chomsky1957syntactic}.

In the language learning literature, dependencies have been studied in very simple constructions without recursive embedding. \citeA{gomez} had participants learn an $aXb$ like language, meaning that the first ($a$) and last ($b$) word were fixed and $X$ was free to vary. They found that it was \emph{easier} to for people to discover the fixed dependency between $a$ and $b$ if the $X$ is \emph{more variable}. In their experiment, $X$ was chosen from a fixed pool of strings; as the cardinality of the pool increased, learners performed better. The result is intuitive: when there are more possible $X$es it is less likely to be predicted since there are more possibilities may free resources for discovering the association between $a$ and $b$.

\begin{figure}[htb]
    \centering

    \begin{subfigure}[]{0.4\textwidth}
    \centering
	\includegraphics[width=\textwidth]{figs/nonadjacent.pdf}
    \caption{}
    \label{non_g}
	\end{subfigure}
    \hspace{-1.6cm}
	\begin{subfigure}[]{0.3\textwidth}
    \centering
	\begin{tabular}[b]{cc}
        \begin{lstlisting}[mathescape,basicstyle=\ttfamily\fontsize{7}{6}\selectfont]
        Small pool size

        def $F$():
            def $F_1$():
                if $flip$():
                    return $pair$($\mb{a}$,$pair$($\mb{h}$,$pair$($\mb{h}$,$pair$($\mb{h}$,$\mb{b}$))))
                else $flip$():
                    return return $pair$($\mb{a}$,$pair$($\mb{h}$,$pair$($\mb{h}$,$pair$($\mb{i}$,$\mb{b}$))))

            def $F_2$($g_1$):
                return $g_1$()

            return $F_2$($F_1$)
        \end{lstlisting}
        \end{tabular}
    \caption{}
    \label{non_example}
	\end{subfigure}
    \begin{subfigure}[]{0.3\textwidth}
    	\centering
        \vspace{-1.4cm}
		\begin{tabular}[b]{cc}
        \begin{lstlisting}[mathescape,basicstyle=\ttfamily\fontsize{7}{6}\selectfont]
        Large pool size

        def $F$():
            def $F_1$():
                x = $\mb{h}$
                for _ in xrange(3):
                    x = $pair$(x,$\mb{h}$) if $flip$() else $pair$(x,$\mb{i}$)
                return x

            def $F_2$($g_1$):
                return $pair$($\mb{a}$,$pair$($g_1$(),$\mb{b}$))

            return $F_2$($F_1$)
        \end{lstlisting}
        \end{tabular}
	\end{subfigure}
    \caption{(a) The graph shows the weighted F-scores of adjacent and nonadjacent dependencies, and their differences with different pool sizes; (b) example hypotheses with small and large pool sizes (only one pair of dependency is shown for simplicity).}
    \label{non}
\end{figure}

We adopt the same paradigms in computational test. There are 3 pairs of nonadjacent dependencies. They are $aXb$, $cXb$ and $cXd$, where $X$ is a string with its length equals 3, and the characters in $X$ are $h$, $i$, $j$, $k$. We construct the pool that $X$ is drawn from by numerating all possible 3-characters strings' permutations (the total amount is, therefore, $3^4=81$). When generating the strings, we draw $X$ uniformly from the pool.

Gomez et al. obtained their results by showing learners strings with both trained and untrained dependencies (ones that do not appear in training strings), and then ask them to choose whether string is valid or not. Then the rates of learners accepting trained and untrained exemplars are computed respectively, and the difference of these two values is taken as the final measure of performance. In our computational test, we can keep track of how much attention our model puts to the nonadjacent dependencies by looking at the posterior probability of those hypotheses that generate such structures. The idea that links these two measures are straightforward: if learner assigns more posterior weight on correct hypotheses, then the probability of generating the untrained data given the posterior distribution state would be smaller, thus less chance to accept them as valid strings.

To get this measure, we first evaluate every hypothesis in our finite hypothesis set for 2048 runs, and check if all generated strings contain valid pairings, such as $aXb$ and $cXd$. In this case, those candidates with all their strings valid are taken as the hypotheses that contain the nonadjacent dependencies. Then we compute the posterior of those hypotheses given different pool size ranging from 2 to 25. The posterior probability curve is shown in Fig.\ref{non}.

This result reveals a similar learning preference as we observed in finite vs. infinite language case. When target language has a small cardinality (i.e. small pool size), the model tends to memorize all occurred strings (examples shown in Fig.\ref{non_example}). Since by this time, the cost of memorizing it still remains low---lower than that of exploring the underlying generative model of the target language. But as the pool size grows, memorizing observed data starts to increase the model complexity. So model tends to explore the potential rules for generating the data, leading to a increase in the posterior probability of those hypotheses containing nonadjacent dependencies.

\section{General discussion \& Conclusion}

Statistical language learning has often been criticized by formal linguists as failing to do justice to the structure of language. Indeed, the earliest statistical learning paradigms \cite{saffran1996statistical, aslin1998computation, saffran1999statistical} focused on domains with little structure such as word segmentation. But this was for good reason: the computational system behind syntax and structure is complex, and a considerable amount is learned about the general capacity of human learners---what general class of theories should science explore---by demonstrating their statistical acquisition in any domain. Indeed, these early statistical learning experiments motivated the development of richer computational models \cite<e.g.>{redington1998distributional}, as well as \citeA{perfors2011learnability}'s work that captured learners as doing model comparison. Perhaps the most interesting critique of Perfors et al. has argued that it ``builds in'' \emph{more} than alternatives such as Chomsky's \emph{Universal Grammar} (UG) because in Perfors et al.'s work, learners consider seven different grammars. Surely it is a more parsimonious theory of human nature that learners come equipped with one single grammar than seven. The arguement also applies to Chater et al.: they allow learners to consider an infinite number of grammars. How could ``building in'' an infinite number of grammars be a more plausible or parsimonious theory than building in just one?

It turns out that this thinking is wrongheaded---but for a very interesting reason. The argument can best be understood by considering Jorge Luis Borges' short story, \emph{The Library of Babel}. Borges imagines an infinite library full of every possible book---every possible sequence of characters. The curious fact about the library which contains \emph{every} possible book is that it actually contains \emph{essentially no information at all.} Its contents can be completely specified with just a few words (``all possible books''), or compressed into an extremely short generating computer program, yielding a negligible complexity or description length. Certainly, it is much easier to describe the entire library of Babel than to describe one of its typical books. In this way, it is often more concise or parsimonious to build in a larger, unconstrained space of hypotheses because a large space can easily have a more concise description (e.g. ``all computations'') compared to a constrained hypothesis space.

The model described in this paper shows in principle how learners can build in very little in terms of representation---just the specific operations in Table 1. This model is the first to show that statistical learning can \emph{construct} representations of provably different computational power, including regular, context-free, and context-sensitive hypotheses and in doing so, it requires surprisingly little data\footnote{There is an interesting sense in which models such as these challenge the very notion of innateness. We could ask: are the infinitely many hypotheses considered here ``innate''? Or is only the \emph{capacity for generating} hypotheses innate? We tend toward the latter sense because the amount that is ``built in'' is very small. When we download a C compiler, it would be very strange to think that all of the programs that could be written are ``innately'' included. Instead, the compiler has a generative capacity to evaluate many such programs, and it is left up to the programmer to construct the specific program that should be run.}.

The overarching view, then, is not one of comparing a few different built-in grammars, but applying likely the same remarkable capacity for discovering algorithms that is apparent in other cognitive domains. As the recent innovations in the Bayesian LOT and inductive computation \cite{chater2007ideal, hsu2010logical, hsu2011probabilistic} make clear, such statistical models are just getting started. However, is likely that ongoing research will push towards more detailed linguistic facts and correspondingly sophisticated learning theories.\footnote{Recently, for example, \citeA{mitchener2010computational} studied \citeA{berwick2011poverty}'s example ``*The child seems sleeping'' and showed that a hierarchical Bayesian model motivated by Perfors et al. and other models showed promise in acquiring the verb's syntax.} The work on the simplified English grammar has shown that the same model can begin to scale to more complex types of computations and syntactic processes.

\section{Conclusion}

We have shown that a statistical learner who operates over computational processes can discover many of the key features of syntactic structures, including regular, context-free, and context-sensitive langauges, using positive evidence only and a simple, unitary learning mechanism. This work opens the door to more sophisticated statistical learning theories that incorporate other areas of language like semantics and pragmatics into the generating process, and unifies theories of syntactic learning with state of the art machine learning models used in other areas of cognition.

\bibliographystyle{apacite}


\section{Introduction}\label{sec:intro}

The nature of dark matter (DM) in the Universe constitutes one of the most pressing questions in contemporary
fundamental physics. It has motivated an impressive joint effort in the field of particle- and astrophysics.
Among the various strategies, indirect detection searches are of particular importance as they
probe the self-annihilating nature of DM, potentially establishing a thermal freeze-out scenario as
an explanation of the measured DM relic density. With the recently published data of the AMS-02
experiment on-board the international space station the cosmic-ray (CR) precision era has begun,
allowing for a significant improvement in the sensitivity of DM searches. However, DM constraints
are largely affected by uncertainties in the description of CR propagation in the Galaxy.
This results in an order-of-magnitude uncertainty in the DM interpretation of CR fluxes
when derived for benchmark propagation models, like the MIN/MED/MAX
scenarios~\cite{Donato:2003xg}.

Here, we summarize our recent work~\cite{Cuoco:2016eej,Cuoco:2017rxb} where we exploit the
precise AMS-02 measurements of the fluxes of protons, helium and antiprotons. By preforming
a joint fit of the propagation parameters and DM signal we minimize the CR propagation
uncertainties. This treatment explores possible correlations and degeneracies between
CRs and DM, providing more robust and reliable DM constraints. While being able to put
strong constraints on the annihilation cross section at large masses, we find a statistically
significant (around $4.5\sigma$ local) indication of a signal from DM annihilation for lower
masses $m_\mathrm{DM}\lesssim100\,\mathrm{GeV}$ pointing to an annihilation cross section close to the
thermal one, $\langle\sigma v \rangle \simeq 3\times 10^{-26}\,\mathrm{cm}^3\mathrm{s}^{-1}$. We analyse
this indication for a variety of annihilation channels. Intriguingly, this possible signal from
DM is largely compatible with the Galactic center excess (GCE) seen in the Fermi-LAT
gamma-ray data.

\section{Cosmic-ray analysis}

The propagation of charged CRs can be described by a diffusion
equation~\cite{StrongMoskalenko_CR_rewview_2007} for the
particle density $\psi_i$ of species $i$ per volume and absolute
value of momentum $p$
\begin{equation}
\begin{split}
  \label{eqn::PropagationEquation}
  \frac{\partial \psi_i (\vec{x}, p, t)}{\partial t} = \;&
    q_i(\vec{x}, p) +
    \vec{\nabla} \cdot \left(  D_{xx} \vec{\nabla} \psi_i - \vec{V} \psi_i \right)
     +  \frac{\partial}{\partial p} p^2 D_{pp} \frac{\partial}{\partial p} \frac{1}{p^2} \psi_i   \\
   &- \frac{\partial}{\partial p} \left( \frac{\mathrm{d} p}{\mathrm{d} t} \psi_i  - \frac{p}{3} (\vec{\nabla \cdot V}) \psi_i \right)
    - \frac{1}{\tau_{f,i}} \psi_i - \frac{1}{\tau_{r,i}} \psi_i\,,
    \end{split}
\end{equation}
where $q_i(\bm{x}, p)$ is the source term of primary CRs and the terms
proportional to $D_{xx},  \bm{V}$, and $D_{pp}$ correspond to CR diffusion,
convection and reacceleration, respectively. The second line of
eq.~\eqref{eqn::PropagationEquation} includes the momentum gain or
loss rate~$\propto {\mathrm{d} p}/{\mathrm{d} t}$, adiabatic energy losses
$\propto \bm{\nabla \cdot V}$, and the loss by fragmentation and radioactive
decay $\propto 1/\tau_{f,i}$ and $1/\tau_{r,i}$, respectively.
We solve eq.~\eqref{eqn::PropagationEquation} numerically using
\textsc{Galprop}~\cite{Strong:1998fr} approximating the Galaxy by
a cylindrical volume assuming a steady state regime.
The source term for primary astrophysical CRs is modeled by a
double broken power law with smooth transitions. We allow for individual
spectral indices for proton and helium. Diffusion is modeled by a power law
in rigidity and convective winds are assumed to be orthogonal to the Galactic
plane. The respective parametrizations are detailed in~\cite{Korsmeier:2016kha}
and introduce a total of 11 background propagation parameters.

The source term for primary antiprotons from DM annihilation can be written as:
\begin{eqnarray}
  \label{eqn::DM_source_term}
  q_{\bar{p}}^{(\mathrm{DM})}(\bm{x}, E_\mathrm{kin}) =
  \frac{1}{2} \left( \frac{\rho(\bm{x})}{m_\mathrm{DM}}\right)^2
  \sum_f \left\langle \sigma v \right\rangle_f \frac{\mathrm{d} N^f_{\bar{p}}}{\mathrm{d} E_\mathrm{kin}} ,
\end{eqnarray}
where $m_\mathrm{DM}$ denotes the DM mass, $\rho(\bm{x})$ the DM density
profile, $\left\langle \sigma v \right\rangle_f$ the thermally averaged annihilation
cross section and $\mathrm{d} N^f_{\bar{p}}/\mathrm{d} E_\mathrm{kin}$ the corresponding
antiproton energy spectrum per DM annihilation, see~\cite{Cuoco:2016eej,Cuoco:2017rxb}
for further details. The factor $1/2$ in eq.~\eqref{eqn::DM_source_term}
corresponds to Majorana fermion DM, the corresponding factor for a Dirac fermion is $1/4$.
For a given annihilation channel $f$, ${\rm DM}\!+\!{\rm DM} \to f\!+\!\bar{f}$,
the source term introduces two additional parameters, $m_\mathrm{DM}$ and
$\left\langle \sigma v \right\rangle_f$.
The resulting 13 parameters are determined in a global fit of the AMS-02
proton and helium fluxes~\cite{Aguilar_AMS_Proton_2015, Aguilar_AMS_Helium_2015}
and the AMS-02 antiproton-to-proton ratio \cite{Aguilar:2016kjl}, supplemented by
proton and helium data from CREAM~\cite{Yoon_CREAM_CR_ProtonHelium_2011} and
VOYAGER~\cite{Stone_VOYAGER_CR_LIS_FLUX_2013}.
We use \textsc{MultiNest} \cite{Feroz_MultiNest_2008} for an efficient sampling.
The VOYAGER data, directly measuring the local interstellar space flux, is used to
constrain the solar modulation potential which we describe in the framework of
the force-field approximation.

\section{Cosmic-ray results}

We are able to put strong constraints on the annihilation cross section at large
DM masses, see figure~\ref{fig:CRfit} (left) for the $b\bar b$ annihilation
channel.\footnote{For CR antiproton constraints on a variety of other annihilation channels
see~\cite{Cuoco:2017iax}.} However, in the mass range around and below
100\,GeV the global fit of the CR antiproton data improves significantly when
adding a contribution of primary antiprotons from DM annihilation. This is due
to a sharp spectral feature in the antiproton flux at a rigidity of about 20\,GV, which is
well fitted by the spectrum from DM annihilation due to the kinematic cut-off set by
the DM mass. The feature cannot be described well by the secondary antiprotons
produced by the interactions of primary protons and helium nuclei on the interstellar
medium, since they only provide a smooth spectrum.

\begin{figure}[ht]
\centering
\setlength{\unitlength}{1\textwidth}
\begin{picture}(1,0.38)
  \put(-0.005,-0.015){\includegraphics[width=1.01\textwidth]{fig1.pdf}}
\end{picture}
\caption{Left:~Best fit region (1--3$\sigma$ contours) for a DM component of the
antiproton flux, and limits on the DM annihilation cross section into $b\bar{b}$
final states; from~\cite{Cuoco:2016eej}. The grey shaded band denotes the
systematic uncertainty. For comparison we show limits on the annihilation
cross section obtained from gamma-ray observations of dwarf
galaxies~\cite{Ackermann:2015zua,Fermi-LAT:2016uux}, and the thermal
value of the annihilation cross section,
$\left\langle \sigma v \right\rangle \simeq 3 \times 10^{-26}$~cm$^3$s$^{-1}$.
Right:~Comparison of the best fit regions (1--3$\sigma$ contours) for the
annihilation channels $gg$~(cyan), $WW^*$~(green), $b\bar b$~(red),
$ZZ^*$~(blue), $hh$~(pink) and $t\bar t$~(orange); from~\cite{Cuoco:2017rxb}.
}
\label{fig:CRfit}
\end{figure}

Figure~\ref{fig:CRfit} (right) shows the best fit regions in the DM mass and
annihilation cross section for the six considered annihilation channels.
The regions are obtained by minimizing the $\chi^2$ with respect to the
remaining 11 parameters in the fit. Hence, they include the uncertainties
in the CR source spectra and CR propagation.
All channels provide an improvement compared to a fit without DM\@.
For the latter we obtain a $\chi^2$/(number of degrees of freedom) of 71/165.
This is reduced to 50/163\,($gg$),  50/163\,($WW^*$), 46/163\,($ZZ^*$)
46/163\,($b \bar b$), 48/163\,($hh$) and 59/163\,($t \bar t$), respectively,
when adding a corresponding contribution from DM annihilation.
Formally, $\Delta \chi^2 \simeq 25$ (obtained for $b \bar b$ and $ZZ^*$)
corresponds to a statistical significance of $4.5\sigma$.\footnote{%
Note that this estimate does not account for possible systematic errors.}
Figure~\ref{fig:CRfit} (right) also shows that different annihilation channels imply
different preferred DM masses, ranging from $m_\mathrm{DM} \simeq 35$\,GeV for gluons
to $m_\mathrm{DM}$ near the Higgs or top mass for annihilation into $hh$ or $t\bar t$,
respectively. Note that for the channels $WW^*$ and $ZZ^*$ the best fit masses
lie below the corresponding on-shell production thresholds. We therefore take
into account the respective 3-body final state processes via one off-shell vector
boson, denoted by the asterisk. For all channels the fit points to an annihilation
cross section close to the value required for a thermal relic,
$\left\langle \sigma v \right\rangle \simeq 3 \times 10^{-26}\,\mathrm{cm}^3\mathrm{s}^{-1}$.

A deeper scrutiny of this potential signal would require a more accurate
determination of the antiproton production cross section, an accurate
modeling of solar modulation at low rigidities of less than about~5\,GV
as well as a better knowledge of correlations in the uncertainties of the
AMS-02 data.

\section{Joint fit with Fermi-LAT gamma-ray data}

Annihilation of DM would also result in a continuum flux of gamma rays,
predominantly from the decay of pions produced in
the fragmentation of primary annihilation products.
Accordingly, in this section we study whether the indication for a
DM annihilation signal in CRs found in the last section is compatible
with recent Fermi-LAT gamma-ray data. In particular, we consider data
from dwarf spheroidal galaxies~\cite{Fermi-LAT:2016uux} and the
so-called Galactic center excess (GCE) recently confirmed
by the Fermi-LAT collaboration~\cite{TheFermi-LAT:2015kwa}.

We consider a total of eleven dwarfs: the seven brightest confirmed dwarfs
analyzed in~\cite{Fermi-LAT:2016uux} as well as Willman~1, Tucana~III,
Tucana~IV and Indus~II. We use the public likelihood provided by
Fermi-LAT~\cite{Fermi-LAT:2016uux} and marginalize over the $J$-factors
according to their uncertainties. For the GCE we use the gamma-ray energy
spectrum and error covariance matrix obtained in~\cite{Calore:2014xka}.
We then perform a joint global fit using the gamma-ray and CR likelihoods
(marginalized over CR propagation parameters).
The gamma-ray fit follows the methodology described in~\cite{Cuoco:2016jqt}.
Besides the two model parameters, $\langle \sigma v \rangle$ and $m_\text{DM}$,
the fit introduces two nuisance parameters -- the $J$-factor for the Galactic center
and the local DM density affecting the DM antiproton flux, which are varied around
$\log (J_\mathrm{GC}/(\text{GeV}^{2}\text{cm}^{-5}))=53.54\pm 0.43$~\cite{Cuoco:2016jqt}
and $\rho_\odot=0.43\pm0.15$ GeV/cm$^3$~\cite{Salucci:2010qr}, respectively.
These two parameters take into account the uncertainty in the DM density profile
and can approximately be considered to be independent~\cite{Cuoco:2017rxb}.

Figure~\ref{fig:genfit} (left) shows the preferred range of DM masses and
annihilation cross sections, where we have marginalized over $J_\mathrm{GC}$ and
$\rho_\odot$. For most annihilation channels, we observe very good
agreement between the best fit regions of the CR signal and the DM
interpretation of the GCE while both are in accordance with limits from
dwarf spheroidal galaxies.
However, annihilation into gluons or top quarks is disfavored as an
explanation of both the CR antiproton flux and the GCE. While both
measurements individually prefer different regions of the DM mass for
annihilation into gluons, annihilation into  $t\bar t$ does not provide a
particularly good fit to both the GCE and the antiprotons.

\begin{figure}[t]
\centering
\setlength{\unitlength}{1\textwidth}
\begin{picture}(1,0.375)
  \put(-0.005,-0.015){\includegraphics[width=1.01\textwidth]{fig2.pdf}}
\end{picture}
\caption{Left: Same as figure~\protect\ref{fig:CRfit} (right) but taking into account the likelihood of Fermi-LAT dwarfs and the GCE in addition. Right: Individual and joint fits within the Higgs portal DM model; from~\cite{Cuoco:2017rxb}.}
\label{fig:genfit}
\end{figure}

\section{Implications for dark matter models}

Realistic DM models often feature an admixture of various annihilation channels.
It is thus not clear \emph{a priori} if a certain combination of annihilation channels,
considered individually in the last sections, can provide a good fit of the observations.
Here we consider the Higgs portal model where a DM field $X$ is directly coupled to
the Higgs bilinear through the operator $X^2H^\dagger H$.
Independent of the nature of $X$, below the Higgs-pair production threshold,
the relative weights of the annihilation channels are solely determined by the DM mass.
Accordingly, we can provide a joint fit of the CR and gamma-ray data within the Higgs portal
model in terms of the two parameters $m_\mathrm{DM}$ and $\langle\sigma v \rangle$. The result is shown in
figure~\ref{fig:genfit} (right), where we have, again, marginalized over $J_\mathrm{GC}$ and $\rho_\odot$.
The preferred region of DM masses is around $m_{\rm DM} \simeq 60$\,GeV,
where the model predicts annihilation into $b\bar b$, $WW^*\!$ and $gg$
with a weight of approximately 70\%, 20\%, and 10\%, respectively.
We find a $\chi_\mathrm{CR}^2$/(number of degrees of freedom) of 47/163
(compared to 71/165 without DM). Intriguingly the preferred mass of
$m_{\rm DM} \simeq 60$\,GeV in the fit coincides with the narrow mass
region around half of the Higgs mass which allows to reconcile existing
constraints from direct detection, the relic density and collider searches,
see~\cite{Cuoco:2017rxb} for a detailed analysis and discussion.

\section{Conclusion}\label{sec:summary}

With the recently released AMS-02 data the cosmic-ray precision era has started. In order to
reduce uncertainties with respect to the MIN/MED/MAX propagation scenarios we performed
a joint fit of propagation and DM parameters. Intriguingly we found a feature in the antiproton
flux which could be a possible hint for DM annihilation in the Galaxy.
We explore a wide class of annihilation channels including $gg$, $b \bar b$, $WW^{*}$,
$ZZ^{*}$, $hh$ and $t\bar t$ which provide a statistical improvement of the fit at the
4$\sigma$ level (3$\sigma$ level for $t \bar t$) with masses ranging from 35\,GeV to
the Higgs- and top-production thresholds. A further investigation of these potential DM
signals would have to examine systematic effects resulting from solar modulation,
secondary antiproton production, and correlations in AMS-02 data.

Interestingly, for most channels the signal is largely compatible with a DM interpretation
of the GCE as well as the constraints from dwarf spheroidal galaxies.
As an example, we perform a joint fit for the specific case of the Higgs portal DM model
which points to a DM mass near half of the Higgs mass -- a region which is also favored
by totally independent constraints, namely the interplay of the relic density and direct
detection limits.


\section{Introduction}

One central and universal problem in combinatorial optimization and
algorithmic graph theory is to find a partition of the vertex set
satisfying some properties.  The most basic formulation is to find a
$2$-partition that separates a source vertex $s$ from a target vertex
$t$, while the crossing edges have minimum cardinality or weight sum,
called the \emph{size} \cite{ford-62-flow}.  Such a partition, called
a minimum ($s,t$)-cut, can be efficiently computed.  This formulation
smoothly extends to the separation of two disjoint subsets of
vertices, both conceptually and computationally.  However, when we
have three or more (sets of) vertices to separate in a pairwise way
with the same minimization objective, it becomes computationally
intractable.

Dahlhaus et al.~\cite{dahlhaus-94-complexity-of-multiway-cuts} and
Cunningham \cite{cunningham-91-multiway-cut} formulated the
\textsc{multiterminal cut} problem, or \textsc{$p$-terminal cut} when
one wants to emphasize the number $p$ of terminals, and initiated its
study.  The classical but still must-read paper
\cite{dahlhaus-94-complexity-of-multiway-cuts} contains several
results.  First and foremost, it proves its NP-hardness, in general as
well as under special restraints, e.g., $p$ being a constant as small
as $3$, or the input graph being planar.  Another result, also
negative, is that the problem, in its general form or restricted to
constant number of terminals, is MAX SNP-hard.  The most cited
positive result of \cite{dahlhaus-94-complexity-of-multiway-cuts} is
an elegant approximation algorithm with ratio $2(1-1/p)$.  The main
observation behind this algorithm is that an objective $p$-partition
contains for any terminal a $2$-partition that separates it from other
terminals; called an \emph{isolating cut} for this terminal.  On the
one hand, a minimum isolating cut for each terminal can be easily
computed and its size is a lower bound for that of any isolating cut
for this terminal; on the other hand, each crossing edge in the
objective $p$-partition is incident to precisely two parts, and thus
counted twice in the size of all isolating cuts.  Therefore, the size
of any objective cut cannot be smaller than half of the total size of
minimum isolating cuts for all terminals.  Interestingly, the
approximation algorithm from this observation coincides with the
linear program (LP) relaxation from the dual of an integer program for
the multicommodity flow \cite{garg-04-weighted-multiway-cut}, which is
known to have half integrality property.
After a sequence of improvements and with great efforts, Buchbinder et
al.~\cite{buchbinder-13-exponential-clocks-and-mwc} managed to improve
the approximation ratio to $1.323$, which breaks another barrier to
this problem, i.e., $4/3$ for the approximation ratio.

Recently, \textsc{multiterminal cut} and its variations receive a lot
of interest from the perspective of parameterized complexity.  Recall
that a problem, parameterized by $k$, is {\em fixed-parameter
  tractable (FPT)} if there is an algorithm with runtime $f(k)\cdot
n^{O(1)}$, where $f$ is a computable function depending only on $k$.
Using the new concept \emph{important separator}, Marx
\cite{marx-06-parameterized-graph-separation} proved this problem,
parameterized by the cut size $k$, is FPT.  In sequel Chen et
al.~\cite{chen-09-multiway-cut} proposed the first single exponential
time algorithm, with runtime $4^k\cdot n^{O(1)}$.  It was further
improved to $2^k\cdot n^{O(1)}$ by Xiao
\cite{xiao-10-multiterminal-cut}.  Unlike its predecessors, which work
indifferently on \textsc{multiterminal cut} as well as its vertex
variation that is commonly known as \textsc{node multiway cut}, Xiao's
approach, however, applies exclusively to \textsc{multiterminal cut}.
A matching algorithm for \textsc{node multiway cut} was later reported
by Cygan et al.~\cite{cygan-11-multiway-cut}, who used a very
interesting branching scheme based on a novel usage of the LP for
multicommodity flow.  With the benefit of hindsight, we are able to
point out that techniques and results of Xiao
\cite{xiao-10-multiterminal-cut} and Cygan et
al.~\cite{cygan-11-multiway-cut} are essentially the same: the former
algorithm can also be re-interpreted with the same branching scheme
using the same flow-based LP.

The cut size $k$ is not the only parameter for \textsc{multiterminal
  cut}, and a natural alternative is $p$, the number of terminals.
However, as aforementioned, in general graphs it is already NP-hard
when $p=3$, hence very unlikely to be FPT.  On the other hand, when
$p$ is a constant, it can be solved in polynomial time on planar
graphs \cite{dahlhaus-94-complexity-of-multiway-cuts}.  Inspired by
this, it is natural to ask for its fixed-parameter tractability
parameterized by $p$ on planar graphs.  Marx
\cite{marx-12-planar-multiway-cut} recently gave a negative answer to
this question.  This apparently was a part of a systematic study of
this problem, of which other outcomes include a $2^{O(p)}\cdot
n^{O(\sqrt{p})}$ time algorithm \cite{klein-12-planar-multiway-cut},
improving previous one in time $n^{O(p)}$
\cite{dahlhaus-94-complexity-of-multiway-cuts}, and a polynomial time
approximation schema \cite{bateni-12-ptas-planar-multiway-cut}.

Naturally, the hardness of this problem should be attributed to the
interference of isolating cuts for different terminals, i.e., an
isolating cut for one terminal might help separate another pair of
terminals.  For example, consider the following unit-weighted graph: a
path of four vertices whose ends are terminals $t_1$ and $t_2$
respectively, and another terminal $t_3$ adjacent to both internal
vertices of the path.  Clearly, any isolating cut for $t_3$ has size
at least $2$, and there are two cuts with the minimum size.  Of which
one is a proper subset of the other, and the larger one with three
vertices would be more desirable.  This example suggests \emph{a
  containment relation of minimum isolating cuts for a terminal} and
\emph{the preference for the largest one among them}.  They had been
formalized in the classical work of Ford and Fulkerson
\cite{ford-62-flow} and Dahlhaus et al.
\cite{dahlhaus-94-complexity-of-multiway-cuts} (see
Theorems~\ref{thm:max-vol-min-cut} and \ref{thm:min-cut-safe} in
Section~\ref{sec:pre}).\footnote{Interestingly, these theorems seem to
  be widely and unfairly ignored, and hence are re-proved again and
  again even decades after their first appearance.  A plausible
  explanation might be that both pieces of work contain other more
  famous results that eclipse them.}

As a matter of fact, the {important separator} is an immediate
generalization of them.  With this technique, several cut related
problems were shown to be FPT.  Some of them, e.g., \textsc{directed
  feedback vertex set} \cite{chen-08-dfvs} and \textsc{multicut}
\cite{marx-14-multicut,thomasse-10-multicut}, had been open for quite
a long time.  Also related is the concept \emph{extreme set}, which
is very well-known in the study of cuts enumerations
\cite{nagamochi-10-algorithmic-connectivity}.

Here we manage to break the $2^k \cdot n^{O(1)}$ barrier, which
has withstood several rounds of attacks.
\begin{theorem}\label{thm:main}
  The \textsc{multiterminal cut} problem can be solved in $1.84^k
  \cdot n^{O(1)}$ time.
\end{theorem}
For a fixed number of terminals, we can even do better, especially
when $p$ is a small constant.
\begin{theorem}\label{thm:p-terminal}
  The \textsc{$3$-terminal cut} problem can be solved in $1.36^k \cdot
  n^{O(1)}$ time.  In general, \textsc{$p$-terminal cut} can be solved
  in $1.84^{(p-2)k/(p-1 )} \cdot n^{O(1)}$ time.
\end{theorem}

The main thrust of our algorithm is a careful analysis of vertices
\emph{close} to a terminal in a sense.  Recall that all previous
parameterized algorithms have the same branching scheme, which
considers each undecided vertex and branches on where to place it.
This operation either increases the size of the minimum isolating cut
for a terminal, or identifies a crossing edge; both can be bounded by
functions of $k$.  We show that a set of vertices whose inclusion to a
terminal increases the minimum isolating cut for it by one can be
grouped together, and then disposed of by a more careful branching
rule.  This rule, together with a measure based analysis
\cite{fomin-11-exact-algorithms}, delivers the claimed algorithm.  We would
like to point out its resemblance with \emph{extreme sets}
\cite{nagamochi-10-algorithmic-connectivity}, which might be used in
the preprocessing of non-terminal vertices, and hence shed some light
on the kernelization of \textsc{multiterminal cut}.  Our approach,
however, does not seem to be generalizable to \textsc{node multiway
  cut} in an easy way.  So we leave it open for a parameterized
algorithm for \textsc{node multiway cut} of time $c^k\cdot n^{O(1)}$
for some constant $c < 2$.

A very interesting generalization of the \textsc{multiterminal cut}
problem was recently proposed by Chekuri and Ene
\cite{chekuri-11-submodular-multiway-partition}, who studied the
\textsc{submodular multiway partition} problem, where the cut size is
replaced by general submodular set functions.  The
\textsc{$p$-terminal cut} problem is also related to the
\textsc{$p$-way cut} problem, which asks for a $p$ partition with the
minimum number/weight sum of crossing edges.  Also NP-hard in general,
it can be solved in polynomial time when $p$ is fixed
\cite{goldschmidt-94-k-cut}, and Kawarabayashi and Thorup
\cite{thorup-11-k-cut} recently proved that it is FPT parameterized by
cut size $k$.

\section{Preliminaries}\label{sec:pre}

For notational convenience, we omit braces for singleton sets, i.e.,
we use $t$ to denote both the element itself and the singleton set
that contains only this element.  As the precise meaning is always
clear from context, this abuse will not introduce ambiguities.  We use
$+$ and $-$ for operations set union and set difference,
respectively.

The graph, given by $(V,E)$, is simple and undirected.  The
cardinality of $V$ is denoted by $n$ throughout the paper.  Each edge
$e\in E$ is weighted by a positive integer $w(e)$; the weight of a
subset $E'$ of $E$ is defined as $w(E') = \sum_{e\in E'} w(e)$.  For
two (possibly intersecting) subsets $X,Y \subset V$, we denote by
$E(X,Y)$ the set of edges with one end in $X$ and the other in $Y$,
and use $w(X,Y)$ as a short hand for $w(E(X,Y))$.  For
a
subset $X$ of $V$,
the ordered partition $(X, V - {X})$ is called a \emph{cut} of $G$,
and denoted by $X$, whose \emph{size} is defined by $ d(X) = w(X, V -
{X})$.  By definition, $ d(\emptyset) = d(V) = 0$.  A cut $X$ is an
\emph{($S,T$)-cut} if $S\subseteq X \subseteq V - {T}$.  An
{($S,T$)-cut} with the minimum size is called a \emph{minimum
  ($S,T$)-cut}.

The following fact was first observed by Ford and Fulkerson
\cite[Section~1.5]{ford-62-flow}, as a by-product of the max-flow
min-cut theorem, and later rediscovered several times by different
authors.

\begin{theorem}[\cite{ford-62-flow}]\label{thm:max-vol-min-cut}
  Let $S, T \subset V$ be two disjoint nonempty subsets of vertices.
  There is a minimum ($S, T$)-cut $X$ such that all other minimum ($S,
  T$)-cuts are subsets of $X$.  Moreover, this cut can be constructed
  in polynomial time.
\end{theorem}

Such an ($S, T$)-cut will be called the \emph{maximum volume minimum
  ($S,T$)-cut}, or \emph{max-vol min-cut for ($S,T$)} for short.  Note
that the definition is asymmetric and the pair ($S, T$) is ordered;
the sets $S$ and $T$ are commonly referred to as the source and target
terminals, respectively.  Given any pair of terminals, their max-vol
min-cut can be found in polynomial time.  Indeed, most known
algorithms for minimum cut return the max-vol min-cut for either
($S,T$) or ($T,S$).  See also
\cite{nagamochi-10-algorithmic-connectivity} for an updated and more
comprehensive treatment.

Generalized in a natural way, given a set ${ T} = \{t_1,\dots,t_p\}$
of $p$ terminals, where $p\ge 3$, a partition $\{V_1,\dots,V_p\}$ of
$V$ is called a \emph{multiterminal cut} for ${T}$ if $t_i\in V_i$ for
each $1\le i\le p$.\footnote{We define multiterminal cut as a
  partition, instead of a disconnecting set, because it is more
  compatible with the original definition of cuts \cite{ford-62-flow},
  (see Sec.~1.6 of \cite{ford-62-flow} for an explanation on why Ford
  and Fulkerson chose partitions over disconnecting sets in the first
  place,) as well as its extension to submodular multiway partitions
  \cite{chekuri-11-submodular-multiway-partition}.}  The size of this
multiterminal cut is defined to be $\sum_{u\in V_i,v\in V_j, 1\le i<
  j\le p}w(u,v)$, that is, the weight sum of \emph{crossing edges},
edges with ends in different parts.  Since each edge is counted in
precisely two parts, the size is equivalent to
$\frac{1}{2}\sum^p_{i=1} d(V_i)$.  We are now ready for a formal
definition of the problem.

\medskip

\fbox{\parbox{0.9\linewidth}{
  \textsc{multiterminal cut} ($G,w,T,k$)

\begin{tabularx}{\linewidth}{rX}
  \textit{Input:} & a graph $G = (V, E)$, a weight function $w : E \to
  \mathbb{N}$, a set $T$ of $p$ distinct terminals, and a nonnegative
  integer $k$.
  \\
  \textit{Task:} & Either find a multiterminal cut for $T$ of size no
  more than $k$, or report no such a multiterminal cut exists.
\end{tabularx}
}}
\medskip

Needless to say, for each $1\le i\le p$, set $V_i$ makes a ($t_i, T -
t_{i}$)-cut.  In other words, $V_i$ isolates $t_i$ from all other
terminals; hence called an \emph{isolating cut for $t_i$}
\cite{dahlhaus-94-complexity-of-multiway-cuts}.  A \emph{min-iso-cut
  for $t_i$} is a short hand for \emph{minimum isolating cut for
  $t_i$}.  Noting that Theorem~\ref{thm:max-vol-min-cut} applies to
isolating cuts, the \emph{max-vol min-iso-cut for $t_i$} can be
defined naturally.

\begin{theorem}[\cite{dahlhaus-94-complexity-of-multiway-cuts}]
  \label{thm:min-cut-safe}
  Let $(G, w, T, k)$ be an instance of \textsc{multiterminal cut}.
  For any min-iso-cut $X$ for $t_i\in T$, there exists a minimum
  multiterminal cut $\{V_1,\dots,V_p\}$ for $T$ such that $X \subseteq
  V_i$.
\end{theorem}

We say that a set function $f: 2^V \to \mathbb{Z}$ is
\emph{submodular} if every pair of subsets $X,Y\subseteq V$ satisfies
\begin{equation}
  f(X) + f(Y) \ge f(X\cap Y) + f(X\cup Y).
\end{equation}

The following property can be easily verified by definition (see,
e.g., \cite{nagamochi-10-algorithmic-connectivity}), and shall be the
swiss army knife in what follows.
\begin{theorem}\label{thm:submodular}
  The cut size function, $d: 2^V\to \mathbb{Z}$, is submodular.
\end{theorem}

\begin{lemma}\label{lem:increment-isolating-cuts}
  Let ${T} = \{t_1, \dots, t_p\}$ be the set of terminals where $t_i$
  is the max-vol min-cut for ($t_i, T - t_i$) for each $1\le i\le p$.
  If the min-iso-cut size for a terminal decreases with the removal of
  an edge $e=uv$ , then any min-iso-cut for this terminal in $G - e$
  contains either $u$ or $v$, and its size decreases by at most
  $w(e)$.  The upper bound is reached if{f}

  the terminal is $u$ or $v$.  Moreover, the sizes of min-iso-cuts for
  at most two terminals in $G - e$ can be smaller than that in $G$.
\end{lemma}
  \begin{proof}
    Let $C_i$ be any min-iso-cut for $t_i$ in $G - e$.  If neither $u$
    nor $v$ is in $C_i$, then we must have $C_i = t_i$, as otherwise
    $d(C_i) > d(t_i)$, which contradicts the definition of
    min-iso-cut; for the same reason, $u$ and $v$ cannot be both
    contained in $C_i$.  Hence we assume $|\{u,v\}\cap C_i| = 1$.  The
    size of min-iso-cuts for $t_i$ decreases by $d(t_i) - (d(C_i) -
    w(e) )= w(e) + ( d(t_i) - d(C_i) ) \le w(e)$.  To achieve the
    upper bound, we must have $C_i = t_i$, which implies $t_i\in
    \{u,v\}$.

    Suppose to the contrary of the last assertion, both $C_i$ and
    $C_j$, where $i\ne j$, contain $u$.  Then neither of them contains
    $v$.  Let us consider $Y = C_i \cap C_j$.  The minimality of $C_i$
    implies $w(Y, C_i - Y) = w(Y, V - C_i)$; likewise $w(Y, C_j - Y) =
    w(Y, V - C_j)$.  We can conclude that $w(Y, V - (C_i \cup C_j)) =
    0$ and then $d(C_i\setminus C_j) = d(C_i)$.  Therefore,
    $C_i\setminus C_j$ is also a min-iso-cut for $t_i$ in $G - e$; a
    contradiction to the first assertion.  We have seen that $u$ is
    contained in the min-iso-cut for at most one terminal, and the
    same argument applies to $v$.  We can thus conclude the last
    assertion and the lemma.
  \end{proof}
\section{Distances to a terminal}
\label{sec:distances}

For the problem under concern, it will make no difference if a
terminal is replaced by a subset of vertices; they can always be
\emph{merged} into a single vertex.  By merging a subset $U$ of
vertices into a single vertex $u$, we remove all vertices in $U$,
introduce a new vertex $u$, and for each vertex $v\in N(U)$, set the
weight of the new edge $uv$ to be $w(U,v)$.  In particular, if $U$
contains a terminal $t_i$, then the newly introduced vertex will be
also identified as $t_i$; note that we never merge a vertex set
containing more than one terminal.  We remark that this operation does
not introduce non-integral weights.  If a non-terminal vertex $v$ is
decided to be in $V_i$, we can merge $\{v,t_i\}$ into $t_i$.
Theorem~\ref{thm:min-cut-safe} might be informally interpreted as: a
vertex whose merge into a terminal does not increase the min-iso-cuts
size for it can always be safely merged.  This observation inspires us
to consider the increment of the size of min-iso-cuts for $t_i$ with
the merge of a non-terminal vertex $v$ to it.  Noting that for two
disjoint sets $S$ and $T$, and a proper subset $S'$ of $S$, the
max-vol min-cut for ($S',T$) is a subset of that for ($S,T$), thereby
having a size no larger than the latter.

\begin{definition}
  Let $v$ be a non-terminal vertex and $t_i$ be a terminal.  Let $C'$
  and $C$ denote the max-vol min-cuts for ($t_i + v, T - t_i$) and,
  respectively, ($t_i, T - t_i$).  The \emph{extension} of $v$ from
  $t_i$ is defined to be $X_i(v) = C' - C$, and the \emph{distance} of
  $v$ from $t_i$ is defined to be $ ext_i(v) = d(C') - d(C)$.
\end{definition}
\\*
Using Theorem~\ref{thm:max-vol-min-cut}, we can compute the extension
and distance for any pair of non-terminal vertex $v$ and terminal
$t_i$ in polynomial time.  By definition, an extension contains no
terminal, and a distance is always nonnegative.  In particular, the
extension is empty if and only if the distance is $0$, then $v\in C$.
Although the set $C'$ might not be connected (even when $C$ is
connected), the extension always is:

\begin{lemma}
  \label{lem:X-connected}
  For any non-terminal vertex $v$ and terminal $t_i$ , the subgraph
  induced by $X_i(v)$ is connected.
\end{lemma}
  \begin{proof}
    We may assume $X_i(v)$ is nonempty.  Suppose, for contradiction,
    that $X_i(v)$ has a partition $\{Y,Z\}$ such that there is no edge
    between $Y$ and $Z$, and without loss of generality, $v \in Y$.
    Let $C$ be the max-vol min-iso-cut for $t_i$.  Since $Y$ and $Z$
    are disconnected and disjoint from $C$, we have $ d(C + Z) >
    d(C)$, and
    \[
     d(C + X_i(v)) -  d(C) =  d(C + Y) -  d(C)
    +  d(C + Z) -  d(C) >  d(C + Y) -  d(C),
    \]
    which means $C + Y$ is a strictly smaller cut for ($t_i + v, T -
    t_i$) than $C + X_i(v) = C + Y + Z$; a contradiction.
  \end{proof}

  Let us now put those vertices in the immediate vicinity to terminal
  $t_i$, i.e., vertices whose distances to $t_i$ are exactly $1$,
  under a closer scrutiny.

\begin{lemma}
  \label{lem:min+1}
  Let $v$ be a vertex at distance $1$ to terminal $t_i$.  There exists
  a minimum {multiterminal cut} for $T$ that keeps $X_i(v)$ in the
  same part.
\end{lemma}

\begin{proof}
  We denote by $C$ and $C'$ the max-vol min-cuts for ($t_i, T - t_i$)
  and ($t_i + v, T - t_i$), respectively, and let $X = X_i(v)$.  We
  prove this lemma by constructing a claimed cut.  Let ${\cal P} =
  \{V_1, V_2, \dots, V_p\}$ be a minimum multiterminal cut for $T$
  satisfying $C\subseteq V_i$; its existence is ensured by
  Theorem~\ref{thm:min-cut-safe}.  If $X\subset V_j$ for some $1\le
  j\le p$, then we are already done; hence we assume otherwise.  We
  claim that ${\cal P}' = \{V'_1, V'_2, \dots, V'_t\}$, where $V'_i =
  V_i \cup X$ and $V'_j = V_j \setminus X$ for $j\ne i$, is also a
  minimum multiterminal cut for $T$.  The remaining discussion is
  separated based on $X \cap V_i = \emptyset$ or not.

  If $X \cap V_i \neq \emptyset$, then $C$ is proper subset of $V_i
  \cap C'$, which implies $d(C) < d(V_i \cap C')$.  Noting all values
  are integral, we have $d(C') = d(C) + 1 \le d(V_i \cap C')$.  By the
  submodularity of the cut function (Theorem~\ref{thm:submodular}), we
  have $d(V_i) + d(C') \geq d(V_i \cup C') + d(V_i \cap C')$.
  Combining them, we get
    \[
      d(V_i) \geq  d(V_i \cup C') = d(V_i \cup X) =  d(V'_i).
    \]
    Noting that $E(V'_j,V'_l)\subseteq E(V_j,V_l)$ for any pair of
    $j,l$ different from $i$, this case is proved.

    Otherwise, $X \cap V_i = \emptyset$.  The crossing edges of ${\cal
      P}'$ but not of ${\cal P}$ are a subset of $E(X, V - {C'})$;
    while the crossing edges of ${\cal P}$ but not of ${\cal P}'$
    include $E(V_i,X)$ and at least one edge $e$ in $G[X]$: noting
    that by Lemma~\ref{lem:X-connected}, $G[X]$ was originally
    connected. Therefore the difference of the cut sizes of ${\cal
      P}'$ and ${\cal P}$ is at most:
    \[
    w(X, V - {C'}) - w(V_i, X) - w(e) \leq ext_i(u) - 1 = 0,
    \]
    which justifies the optimality of ${\cal P}'$ and finishes the
    proof.

\end{proof}

Lemma~\ref{lem:min+1} permits us to merge each extension $X_i(v)$ with
$ext_i(v) = 1$ into a single non-terminal vertex.  After both the
max-vol min-iso-cut for $t_i$ and $X_i(v)$ have been merged into $t_i$
and $v$ respectively, either they are adjacent, or the new vertex $v$
is only incident to a unit-weight edge.  As vertices of degree $1$ can
always be removed safely, we may assume that all vertices at unit
distances to $t_i$ are its neighbors.

\section{The algorithm}

From Theorems \ref{thm:max-vol-min-cut} and \ref{thm:min-cut-safe},
one can easily derive a bounded search tree algorithm as follows.  For
$1\le i\le p$, we initialize $V_i$ to be the max-vol min-iso-cut for
$t_i$ and merge it into $t_i$.  Then we grow it by including its
neighbors one by one, until no vertex is left out.  On each neighbor
$v$ of $t_i$, we have two options: including it and merging
$\{v,t_i\}$ into $t_i$, or excluding it and then edge $(t_i,v)$ is a
crossing edge.
This branching rule was first presented by Chen et
al.~\cite{chen-09-multiway-cut} and later used by Xiao
\cite{xiao-10-multiterminal-cut}.  The original analysis in
\cite{chen-09-multiway-cut} uses $2k - d(t_i)$ as the measure, which
delivers a $4^k$ bound on the number of leaves the algorithm traverses
in the search tree.  Inspired by the observation on isolating cuts,
Xiao \cite{xiao-10-multiterminal-cut} used the new measure $2k - \sum
d(t_i)$.  As it always holds that $\sum d(t_i)>k$ in a nontrivial
instance, the bound can be tightened to $2^k$.  This analysis is then
tight: there are instances of \textsc{multiterminal cut} on which this
branching rule will make a search tree of $2^k$ leaves.  In other
words, the base $2$ cannot be further lowered by analysis.

Let $m$ be the measure $2k - \sum d(t_i)$.  The bound $2^k$ follows
from $m\le k$ and the branching vector ($1, 1$): it forks into at most
two branches, in each of which the measure decreases by at least
$1$.\footnote{We shall use standard technique to analyze the bounded
  search tree algorithm (see, e.g.,
  \cite[Theorem~2.1]{fomin-11-exact-algorithms}).  We say that a branching rule
  has \emph{branch vector} ($\tau_1, \tau_2,\dots, \tau_r$) if, given
  an input instance of measure $m$, it branches into $r$ instances,
  and the measures in them are at most $m-\tau_1, m-\tau_2, \dots,
  m-\tau_r$, respectively.  With such a branch vector, the linear
  recurrence for the maximum number of leaves is
\[
{\cal T}(m) \le {\cal T}(m - \tau_1) + {\cal T}(m - \tau_2) + \dots
+{\cal T}(m - \tau_r),
\]
and then ${\cal T}(m) \le c^m$, where $c$ is the unique positive real
root of
\[
z^m - z^{m - \tau_1} - z^{m - \tau_2} - \dots - z^{m - \tau_r} = 0.
\]} To break the $2^k$ barrier, a branching vector strictly better
than ($1, 1$) is in need.  Observing that the inclusion of a vertex
$v$ always increases $\sum d(t_i)$ by the distance of $v$ to the
terminal, vertices at distance $2$ or more will not concern us; hence
we only need to concentrate on those vertices at unit distance to a
terminal.  On these vertices we design new branching rules, worst of
which has branching vector ($1, 2, 3$); thereby achieving the claimed
bound of Theorem~\ref{thm:main}.

Let $v$ be a non-terminal vertex.  If it has a unique neighbor $u$,
then $v$ and $u$ will always be together; if it has two neighbors $u$
and $u'$, and without loss of generality, $w(v,u) \ge w(v,u')$, then
there always exists a minimum multiterminal cut for $T$ that keeps $u$
and $v$ together.  This holds regardless whether $u$ and/or $u'$ is a
terminal or not.
In both cases, we may merge $v$ to $u$.  Now we are ready to present
our algorithm in Figure~\ref{fig:alg}.  Every time when the terminal
$t_p$ has been separated from the rest of the graph, i.e., the part
$V_p$ has been identified, we treat the rest of the graph as a new
instance with one less terminal.  Nonessential particulars on the
bookkeeping and recovering of merge operations is omitted for the
simplicity.  When the value of $h = \sum d(t_i)$ changes, it is marked
at the end of the line that makes it happen; $c^+$ means at least $c$.

\begin{figure*}[h]
\setbox4=\vbox{\hsize28pc \noindent\strut
\begin{quote}
  \vspace*{-5mm} \footnotesize
  {\bf Algorithm mcut($G, w, {T},k$)}\\
  {\sc input}: a graph $G$, a weight function $w$, a set $T$ of $p$
  terminals, and an integer $k$.

  \\
  {\sc output}: a multiterminal cut $\{V_1,\dots,V_p\}$ for $T$ of
  size at most $k$, or ``NO.''
  \\
  \\
  {\bf preprocessing}
  \\
  0.0 \hspace*{1mm} {\bf if} $p = 2$ {\bf then} find a min
  ($t_1,t_2$)-cut $X$ and {\bf return} ($X,V-X$);
  \\
  0.1 \hspace*{1mm} {\bf for} $ i = 1, \dots, p$ {\bf do}
  \\
  \hspace*{8mm} compute the max-vol min-iso-cut for $t_i$, and merge
  it into $t_i$;
  \\
  0.2 \hspace*{1mm} {\bf for} each vertex $v$ at distance $1$ to some
  terminal $t_i$ {\bf do}
  \\
  \hspace*{8mm} merge $X_i(v)$ into a single vertex $v$;
  \\
  0.3 \hspace*{1mm} $h = \sum d(t_i)$;
  $\quad$ \parbox[t]{0.85\linewidth}{ {\bf if} $2k < h$ {\bf then
      return ``NO''};
    \\
    {\bf if} $k \ge h$ {\bf then return} $\{t_1, \dots, t_{p-1}, V -
    {T} + t_p \}$;}
  \\
  0.4 \hspace*{1mm} dispose of non-terminal vertices of degree $\le
  2$;
  \\

  1 \hspace*{3mm} {\bf for} $ i = 1, \dots, p-1$ {\bf do}
  \\
  \hspace*{9mm} remove $E(t_i, T - t_i)$; $\quad k' = w(t_i, T -
  t_i)$; $\quad k=k- k'$; \hfill $\setminus\!\!\!\setminus h\downarrow
  2k'$
  \\
  2 \hspace*{3mm} {\bf if} $t_p$ has no neighbor {\bf then return}
  mcut($G - t_p, w, {T - t_p}, k$) and $t_p$;
  \\
  \hspace*{5mm} {\bf else} pick a vertex $v\in N(t_p)$;
  \\
  3 \hspace*{3mm} {\bf if} $ext_p(v) = 1$ and $N(v)\subseteq T$ {\bf
    then}
  \\
  \hspace*{9mm} merge $v$ into $t_p$; $\quad$ {\bf return} mcut($G, w,
  {T}, k$); \hfill $\setminus\!\!\!\setminus h\uparrow 1$
  \\*
  {\bf branching} $\qquad\qquad
  \setminus\!\!\setminus$ $N(v)\setminus T\ne\emptyset$ when $ext_p(v)
  = 1$
  \\
  4 \hspace*{3mm} {\bf case 1:}
  \\
  \hspace*{9mm} remove $(t_p, v)$; $\quad$ $k' = w(t_p, v)$; \hfill
  $\setminus\!\!\!\setminus h\downarrow 2k'-(1^+)$
  \\
  \hspace*{9mm} {\bf return} mcut($G, w, {T}, k - k'$);
  \\
  5 \hspace*{3mm} {\bf case 2a [$ext_p(v) = 1$ \& $w(v, t_j) = w(v,
    t_p)$ \& $\{u\} = N(v) \setminus T$]:}
  \\
  \hspace*{5mm} $\setminus\!\!\setminus$ $t_j$ is another terminal.
  \\

  \hspace*{9mm} merge $u,v$ into $t_p$; $\quad$ {\bf return} mcut($G,
  w, {T}, k$); \hfill $\setminus\!\!\!\setminus h\uparrow 2^+$
  \\
  6 \hspace*{3mm} {\bf case 2b [$ d(v)=3$ \& $\{u_1,u_2\} = N(v)
    \setminus T$]:}
  \\
  \hspace*{5mm} $\setminus\!\!\setminus$ in this case all three edges
  incident to $v$ have unit-weight, and we must have $ext_p(v)=1$.
  \\
  \hspace*{9mm} {\bf case 2b.1:} remove edge $(u_1, v)$; merge $v,u_2$
  into $t_p$; \hfill $\setminus\!\!\!\setminus h\uparrow 1$
  \\
  \hspace*{15mm} {\bf return} mcut($G, w, {T}, k - 1$);
  \\
  \hspace*{9mm} {\bf case 2b.2:} merge $v,u_1$ into $t_p$; \hfill
  $\setminus\!\!\!\setminus h\uparrow 2$
  \\
  \hspace*{15mm} {\bf return} mcut($G, w, {T}, k$);
  \\
  7 \hspace*{3mm} {\bf case 2c [otherwise]:}
  \\
  \hspace*{9mm} merge $v$ into $t_p$; {\bf return} mcut($G, w, {T},
  k$); \hfill $\setminus\!\!\!\setminus h\uparrow ext_p(v)$

\end{quote} \vspace*{-6mm} \strut} $$\boxit{\box4}$$
\vspace*{-9mm}
\caption{Algorithm for \textsc{multiterminal cut}}
\label{fig:alg}
\end{figure*}

In one case of our analysis, we will need the following lemma.
\begin{lemma}\label{lem:case1}
  In case 1 (step 4) of algorithm {\bf mcut}, $h$ decreases by at
  most $2w(t_p, v) - 1$.  Moreover, when $w(t_p, v) > 1$ and $ext_p(v)
  = 1$, the value of $h$ decreases by $2 w(t_p, v) - 1$ if and only if
  there is another terminal $t_j$ such that $w(t_p, v) = w(t_j, v)$.
\end{lemma}
\begin{proof}
  As all edges between $t_p$ and other terminals have been removed in
  step 1, the first assertion follows immediately from
  Lemma~\ref{lem:increment-isolating-cuts}.  The ``if'' direction of
  the second assertion is clear, and we now prove the ``only if''
  direction.

  Assume $t_j$ is the other affected terminal, then the size of
  min-iso-cuts for $t_j$ decreases by $w(t_p, v) - ext_j(v)$; thus
  $ext_j(v) = 1$, and $v$ and $t_j$ are adjacent.  Noting that
  $ext_p(v)$ is defined to be
    \begin{align*}
      & w(v, V - t_p) - w(v, t_p)\\
      =& w(v, t_j) + w(v, V - t_p - t_j) - w(v, t_p)\\
      =& w(v, V - t_j) - 1 + w(v, V - t_p - t_j) - w(v, t_p)\\
      =& w(v, V - t_p - t_j) - 1 + w(v, V - t_p - t_j),
    \end{align*}
    we can conclude $w(v, V - t_p - t_j) = 1$, and $w(t_p, v) = w(t_j,
    v) = \frac{ d(v) - 1}{2}$.

\end{proof}

\begin{theorem}
  Algorithm {\bf mcut} solves \textsc{multiterminal cut} in time
  $1.84^k \cdot n^{O(1)}$.
\end{theorem}
\begin{proof}
  Let us first verify its correctness.  Step 0.0 ensures that the
  instance has at least three terminals; otherwise it solves the
  instance by finding a minimum ($t_1, t_2$)-cut.  The correctness of
  preprocessing steps~0.1 and 0.2 follows immediately from
  Theorem~\ref{thm:min-cut-safe} and Lemma~\ref{lem:min+1},
  respectively.  After that, on any given multiterminal cut $\{V_1,
  \dots, V_p\}$ for $T$, we have $d(V_i)\ge d(t_i)$ for each $1\le
  i\le p$; in other words, this cut has size $1/2\sum^p_{i=1} d(V_i)
  \ge 1/2\sum^p_{i=1} d(t_i)= h/2$.  Therefore, if $h > 2k$ we can
  safely report that there is no solution of size $k$ for the input
  instance.  On the other hand, the multiterminal cut $\{t_1, \dots,
  t_{p-1}, V - {T} + t_p \}$ has size no more than $h$, and hence can
  be returned as a solution if $h \leq k$.  This justifies step 0.3.
  Step 0.4 is straightforward.  A posteriori, after the preprocessing,
  the max-vol min-iso-cut for each terminal consists of itself; and
  any non-terminal vertex has at least three neighbors.

  Steps 1 and 2 are clear: every edge connecting two terminals is a
  crossing edge; an isolated terminal can be disregarded for further
  consideration.  Step 3 takes care of the case that $v$ has
  unit-distance to $t_p$ and is only adjacent to terminals.  Since we
  have disposed of non-terminal vertices of degree $1$ or $2$ in step
  0.4, and step 1 deletes only edges between terminals, there are at
  least two other terminals adjacent to $v$.  Noting $ext_p(v) = 1$
  and all weights are integral, it follows that $w(v,t_p)$ is the
  largest among all edges incident to $v$, and thus putting $v$ into
  $V_p$ will minimize crossing edges incident to $v$.  As $v$ is
  nonadjacent to any other non-terminals, its choice has no effect on
  other vertices.  This justifies step 3.

  Hereafter it takes the branching steps, and enters exclusively one
  case of them, which then calls the algorithm recursively.  For each
  neighbor $v$ of $t_p$, we have to either include $v$ to $V_p$ or
  count $(v,t_p)$ as a crossing edge.  Step 4 deals with the latter
  option, where we remove the edge and decrease the parameter
  accordingly.  The other option, where $v$ is put into $V_p$, is
  handled by steps 5-7.

  Step 5 handles the case where $v$ is balanced in two terminals $t_p$
  and $t_j$, then there must be another non-terminal vertex $u$ such
  that $w(v,u) = 1$.  As $w(v, t_p) = w(v, t_j)\ge 1 = w(v, u)$, we
  may assume $v$ is in either $V_p$ or $V_j$.  If $u$ is not in the
  same part as $v$, then it does not matter where to put $v$.  We may
  put it in $V_j$ and thus count $(v,t_p)$ as a crossing edge, which
  is already covered by case 1 (step 4).  Therefore, we may assume
  that $v$ is in $V_p$ only if $u$ is also in $V_p$.

  Step 6 handles the case where $v$ is adjacent to three neighbors
  each with a unit-weight edge.  If an optimal partition puts the
  three neighbors of $v$ into different parts, then $v$ must be in one
  of them and it does not matter which one.  This situation has
  already been covered in case 1.  Hence we assume either $N(v)$ is in
  the same part, which is covered by case 2b.2; or separated into two
  parts, which is covered by case 1, case 2b.1, and case 2b.2
  respectively.

  The correctness of step 7 is straightforward.  This completes the
  proof of the correctness of the algorithm.

  It is clear that every step takes polynomial time.  Hence we will
  concentrate on the branching steps, where it makes 2 or 3 recursive
  calls to itself.  In particular, it either executes case 1, or one
  of the cases 2a, 2b, and 2c.  To precisely determine the time
  complexity, we define $m = 2k - h$ as the measure, and consider its
  decrease in each branch.  By assumption that $t_p$ is the max-vol
  min-iso-cut for $t_p$, and $X_p(v) = v$ when $ext_p(v) = 1$, for
  every other non-terminal vertex $u$ (different from $v$), we always
  have
  \begin{equation}\label{eq:2}
    d(\{t_p, v, u\})\ge d(t_p) + 2.
  \end{equation}
  Also note that the measure will not increase as a side effect in
  other steps.

  In case 1, $k$ decreases by $w(t_p,v)$; by Lemma~\ref{lem:case1},
  $h$ decreases by at most $2w(t_p,v) - 1$.  In total, $m$ decreases
  by at least $1$.
  \begin{itemize}
  \item Case 2a.  According to \eqref{eq:2}, $h$ increases by at least
    $2$; no edge is deleted and $k$ is unchanged.  Thus, $m$ decreases
    by at least $2$.
  \item Case 2b is further separated into two sub-cases.
    \begin{itemize}
    \item Case 2b.1.  According to \eqref{eq:2}, the min-iso-cut size
      for $t_p$ increases by 1.  By
      Lemma~\ref{lem:increment-isolating-cuts} and noting that the
      only deleted edge $u_1 v$ has unit-weight and $u_1$ is not a
      terminal, the min-iso-cut for other terminals are unchanged.
      Thus, $h$ increases by at least $1$.  As $k$ decreases by $1$,
      the measure decreases by at least $3$.
    \item Case 2b.2.  According to \eqref{eq:2}, $h$ increases by at
      least $2$ while $k$ remains unchanged.  Thus, $m$ decreases by at
      least $2$.
    \end{itemize}
  \item Case 2c.  The value of $h$ increases by $ext_p(v)$ while $k$
    remains unchanged; hence $m$ decreases by $ext_p(v)$.  According
    to Lemma~\ref{lem:case1}, $h$ decreases by $1$ in case 1 only if
    the condition of case 2b is satisfied.  In other words, if $m$
    decreases $1$ in case 2c, then it decreases by at least $2$ in
    case 1.
  \end{itemize}

  In summary, the worst case is when the algorithm branches into case
  1 and 2b.  In this situation, this algorithm makes three recursive
  calls to itself, with new measures $m-1$, $m-2$, and $m-3$,
  respectively.  Recall that we only start the branching when $h > k$,
  which means the initial value of the measure is $m = 2k - h < k$.
  We terminate the branching by returning ``NO'' as long as $m$
  decreases to a negative value, which means that the measure
  decreases by at most $k$, and therefore the number of recursive
  calls is upper bounded by $1.84^k$.  This finishes the analysis of
  time complexity.
\end{proof}

This theorem implies Theorem~\ref{thm:main}.

Given a set of isolating cuts $(C_1,\dots,C_p)$, taking any $p - 1$ of
them while leaving all remained vertices in one part will make a
multiterminal cut.  If we omit the one with the largest size, then the
cut obtained as such will have size at most $(1-1/p)\sum^p_{i=1}
d(C_i)$.  This observation was first used to improve the approximation
ratio from $2$ to $2(1-1/p)$
\cite{dahlhaus-94-complexity-of-multiway-cuts}, and is also applicable
for our analysis.  The improvement is especially significant when $p$
is small.

\begin{proof}[Proof of Theorem~\ref{thm:p-terminal}]
  The exit condition $k \geq h$ in step 0.3 of algorithm \textbf{mcut}
  can be replaced by $k \geq (1-1/p) h$.  In other words, we only
  branch when $\frac{p}{p-1}k < h \leq 2k$, which means the measure $m
  = 2k - h$ can decrease by at most
  \[
  2k - \frac{p}{p-1}k = \frac{p-2}{p-1}k,
  \]
  from which and the branching vector
  ($1,2,3$) Theorem~\ref{thm:p-terminal} follows.

\end{proof}

\small