text:
\section{Introduction}\label{section-introduction}
Bayesian inference software such as \proglang{Stan}~\citep{stan-software:2014}, \proglang{OpenBUGS}~\citep{thomas2006making}, and \proglang{JAGS}~\citep{plummer2004jags} provide high-level, domain-specific languages (DSLs) to specify and sample from probabilistic Directed Acyclic Graphs (DAGs). In some Bayesian projects, the convenience of using such DSLs comes at the price of reduced flexibility in model specification, and suboptimality of the underlying sampling algorithms used by the compilers. Furthermore, for large projects the end-goal might be to implement all or part of the sampling algorithm in a high-performance - perhaps parallel - language. In such cases, researchers may choose to start their development work by `rolling their own' joint probability distributions from the DAG specification, followed by application of their choice of a sampling algorithm to the joint distribution.

Many Monte Carlo Markov Chain (MCMC) algorithms have been proposed over the years for sampling from complex posterior distributions. Perhaps the most widely-known algorithm is Metropolis~\citep{metropolis1953equation} and its generalization, Metropolis-Hastings (MH)~\citep{hastings1970monte}. These multivariate algorithms are very easy to implement, but they can be slow to converge without a carefully-selected proposal distribution. A particular flavor of MH is the Stochastic Newton Sampler~\citep{qi2002hessian}, where the proposal distribution is a multivariate Gaussian based on the second-order Taylor series expansion of the log-probability. This method has been implemented in the \proglang{R} package \pkg{sns}~\citep{mahani2014sns}. This algorithm can be quite effective for twice-differentiable, log-concave distributions such as those encountered in Generlized Linear Regression (GLM) problems. Hamiltonian Monte Carlo (HMC) algorithms~\citep{girolami2011riemann,neal2011mcmc} have also gained popularity due to development of techniques for automated tuning of their parameters~\citep{hoffman2014no}.

Univariate samplers tend to have few tuning parameters and thus well suited for black-box MCMC software. Two important examples are Adaptive Rejection Sampling~\citep{gilks1992adaptive} (or ARS) and Slice Sampling~\citep{neal2003slice}. ARS requires log-density to be concave, and needs the first derivative, while slice sampler is generic and derivative-free. To apply these univariate samplers to multivariate distributions, they must be applied one-coordinate-at-a-time according to the Gibbs sampling algorithm~\citep{geman1984stochastic}, where at the end of each univariate step the sampled value is used to update the conditional distribution for the next coordinate. \pkg{MfUSampler} encapsulates this logic into a library function, providing a fast and reliable path towards Bayesian model estimation for researchers working on novel DAG specifications.

When posterior distribution exhibits strong correlation structure, one-coordinate-at-a-time algorithms can become inefficient as they fail to capture important geometry of the space~\citep{girolami2011riemann}. This has been a key motivation for research on black-box multivariate samplers, such as adaptations of slice sampler~\citep{thompson2011slice}.

The rest of this article is organized as follows. In Section~\ref{section-implementation} we provide a brief overview of the extended Gibbs sampling framework used in \pkg{MfUSampler}. In Section~\ref{section-usage} we illustrate how to use the software with an example. Section~\ref{section-workflow} shows how \pkg{MfUSampler} can be used to expedite the prototyping stage in Bayesian modeling problems. Finally, Section~\ref{section-summary} provides a summary and concluding remarks.

\section{Implementation}\label{section-implementation}
\pkg{MfUSampler} relies on three components:
\begin{enumerate}
\item Univariate MCMC samplers: As of version 0.9.1 of \pkg{MfUSampler}, two such samplers are supported: Univariate Slice Sampler with Stepout and Shrinkage~\citep{neal2003slice} and Adaptive Rejection Sampler~\citep{gilks1992adaptive}. For slice sampler, we have imported - with small modifications - Radford Neal's \proglang{R} code, posted on his website\footnote{\url{http://www.cs.toronto.edu/~radford/ftp/slice-R-prog}}, while for ARS we use the \proglang{R} package \pkg{ars}~\citep{rodrigues2014ars}. For technical details on the univariate sampling algorithms, see aforementioned publications or statistical textbooks~\citep{robert1999monte}.
\item Main sampling routine, \code{MfU.Sample}: This function is essentially a \code{for} loop for applying the underlying univariate sampler to each coordinate of the multivariate distribution. We refer to this as `extended Gibbs sampling' (Section~\ref{subsection-extended-gibbs}). The function \code{MfU.Control} allows the user to set the tuning parameters of the univariate sampler. Both \code{MfU.Sample} and \code{MfU.Control} are public functions, and their detailed behavior can be examined by consulting the package documentation.
\item Wrapper functions \code{MfU.fEval}, \code{MfU.fgEval.f} and \code{MfU.fgEval.g} (all internal functions) that return the conditional distribution and its gradients for each coordinate, using the underlying joint distribution and its gradient vector (Section~\ref{subsection-conditional-joint}).
\end{enumerate}

\subsection{Extended Gibbs Sampling}\label{subsection-extended-gibbs}
The \code{for} loop inside \code{MfU.Sample} is a direct implementation of Gibbs sampling~\citep{bishop2006pattern}, with one conceptual extension: rather than requiring an independent sample from each coordinate's conditional distribution, we expect a Markov transition for which the conditional distribution is an invariant distribution. Among the current univariate samplers implemented in \pkg{MfUSampler}, Adaptive Rejection Sampler produces a standard Gibbs sampler while the Slice Sampler falls under the `extended' Gibbs sampler. The following lemma provides the proof of invariance. (For a discussion of ergodicity for slice sampler, see ~\cite{roberts1999convergence}).
\begin{lemma}
If a coordinate-wise Markov transition leaves the conditional distribution invariant, it will also leave the joint distribution invariant.
\end{lemma}

\begin{proof}
The premise can be mathematically expressed as
\begin{equation}\label{equation-premise}
p(x'_k | \mathbf{x}_{\setminus k}) = \int_{x_k} T(x'_k, x_k | \mathbf{x}_{\setminus k}) p(x_k | \mathbf{x}_{\setminus k}) \, \mathrm{d} x_k,
\end{equation}
while the conclusion can be expressed as
\begin{equation}\label{equation-conclusion}
p(x'_k , \mathbf{x}_{\setminus k}) = \int_{x_k} T(x'_k, x_k | \mathbf{x}_{\setminus k}) p(x_k , \mathbf{x}_{\setminus k}) \, \mathrm{d} x_k.
\end{equation}
In the above $\mathbf{x}_{\setminus k}$ denotes all coordinates except for $x_k$ and $T(x'_k, x_k | \mathbf{x}_{\setminus k})$ denotes the coordinate-wise Markov transition density from $x'_k$ to $x_k$. Employing the Product Rule of Probability, we have $p(x_k , \mathbf{x}_{\setminus k}) = p(x'_k | \mathbf{x}_{\setminus k}) \times p(\mathbf{x}_{\setminus k})$. Since the coordinate-wise Markov transition does not change $\mathbf{x}_{\setminus k}$, we can factor $p(\mathbf{x}_{\setminus k})$ out of the integral, thereby easily reducing Equation ~\ref{equation-conclusion} to Equation~\ref{equation-premise}.
\end{proof}
Note that standard Gibbs sampling is a special case of the above lemma where $T(x'_k, x_k | \mathbf{x}_{\setminus k}) = p(x'_k | \mathbf{x}_{\setminus k})$. (The reader can easily verify that this special transition density satifies the premise.) A full Gibbs cycle is simply a succession of coordinate-wise Markov transitions, and since each one leaves the target distribution invariant according to the above lemma, same is true of the resulting composite Markov transition density.

\subsection{Producing the Conditional Distributions}\label{subsection-conditional-joint}
The internal function \code{MfU.fEval} is responsible for producing coordinate-wise conditional distributions used by the slice sampler:
\begin{Schunk}
\begin{Sinput}
R> MfU.fEval <- function(xk, k, x, f, ...) {
+    x[k] <- xk
+    return (f(x, ...))
+  }
\end{Sinput}
\end{Schunk}
The implementation is deceptively simple, but warrants some explanation. The funcion accepts \code{xk} -  the value of the \code{k}'th coordinate - and inserts it into the $K$-dimensional vector \code{x}. It then evaluates and returns the joint distribution \code{f} at \code{x} (fixed arguments are passed via \code{...}). Returning the joint distribution in lieu of the conditional distribution is correct because, from the perspective of each coordinate, the two are proportional:
\begin{equation}
p(x_k | \mathbf{x}_{\setminus k}) = \frac{p(x_k , \mathbf{x}_{\setminus k})}{p(\mathbf{x}_{\setminus k})} \propto p(x_k , \mathbf{x}_{\setminus k})
\end{equation}
where the first step follows from the Product Rule and in the second step, we have taken advantage of the fact that $p(\mathbf{x}_{\setminus k})$ is constant in terms of $x_k$. A multiplicative constant for density translates into an additive constant for log-density, and can be safely ignored in most MCMC algorithms, including the slice sampler and ARS.

Functions \code{MfU.fgEval.f} and \code{MfU.fgEval.g} produce log-density and its gradient for each coordinate, to be consumed by ARS. Since \proglang{R} package \pkg{ars} expects vector forms of input and output for log-density and its gradient, the above two functions implement such vectorization.

\section[Using MfUSampler]{Using \pkg{MfUSampler}}\label{section-usage}
In this section, we illustrate the use of \pkg{MfUSampler} using both slice sampler and ARS as underlying univariate distributions.
\subsection{Example 1: Bayesian Logistic Regression}\label{subsection-usage-example}
We use a Bayesian logistic regression problem ($N$ observations, $K$ coefficients) as an example. The DAG corresponding to this simple problem is shown in Figure~\ref{figure-logistic-dag}. We assume a non-informative, Gaussian prior on each of the $K$ regression coefficients in $\boldsymbol\beta$ (assuming proper scaling of the covariate matrix $\mathbf{X}$), using a mean of $\mu=0.0$ and a standard deviation of $\sigma=1e+6$. The full joint distribution corresponding to this DAG is given in Equation~\ref{equation-logprob-logistic}.

\begin{figure}[t]
\centering \includegraphics[width=0.4\textwidth]{fig_logistic_dag.jpg}
\caption{
Directed Acyclic Graph representing a Bayesian logistic regression problem, with $N$ observations and $K$ coefficients. A non-informative Gaussian prior is imposed on each of the $K$ elements of $\boldsymbol\beta$ using mean $\mu=0.0$ and standard deviation $\sigma=$1e+6. Plate notation is used to show the links from $\boldsymbol\beta$ to each of the $N$ response values $y_n$.
}
\label{figure-logistic-dag}
\end{figure}

\begin{equation} \label{equation-logprob-logistic}
L(\boldsymbol\beta) = - \sum_{n=1}^N \left\{ (1-y_n) \mathbf{x}_n^t \boldsymbol\beta + \mathrm{log}(1 + \exp(-\mathbf{x}_n^t \boldsymbol\beta)) \right\} - \frac{1}{2 \sigma^2} \sum_{k=1}^K (\beta_k - \mu)^2 + C,
\end{equation}
where the first term corresponds to the log-likelihood, the second term represents the log-prior, and $C$ captures any terms that are independent of $\boldsymbol\beta$. Note that, while the log-prior term is separable in $\boldsymbol\beta$, same is not true for log-likelihood. This means the joint log-distribution $L(\boldsymbol\beta)$ cannot be written as $\sum_k L_k(\beta_k)$. Therefore, conditional distributions used in each of the $K$ steps in a Gibbs sampling cycle will depend on the latest draws from the other $K-1$ elements.

\subsection[]{Slice Sampling using \pkg{MfUSampler}}\label{subsection-usage-slice}
First, we must load the package into an \proglang{R} session. We also set the random seed for reporduciblity of results:
\begin{Schunk}
\begin{Sinput}
R> library("MfUSampler")
R> set.seed(0)
\end{Sinput}
\end{Schunk}
Applying slice sampler to a log-density in \pkg{MfUSampler} is quite straightforward. First we must implement a function that returns the log-density. For the log-density of Equation~\ref{equation-logprob-logistic}, we can do:
\begin{Schunk}
\begin{Sinput}
R> logit.f <- function(beta, X, y, mu=0.0, sigma=1e+6) {
+    Xbeta <- X
+    return (-sum((1-y) * Xbeta + log(1 + exp(-Xbeta)))
+            - sum((beta - mu)^2)/(2*sigma^2))
+  }
\end{Sinput}
\end{Schunk}
where we have ignored the constant term $C$. Note that the first argument to \code{logit.f} is the argument for which we want to draw samples, which is \code{beta} in this case. Also note that \code{X} is a $N$-by-$K$ matrix whose n'th row corresponds to $\mathbf{x}_n$ in Equation~\ref{equation-logprob-logistic}. To test \pkg{MfUSampler} on \code{logit.f}, we first generate some simulated data:
\begin{Schunk}
\begin{Sinput}
R> N <- 1000
R> K <- 5
R> X <- matrix(runif(N*K, -0.5, +0.5), ncol=K)
R> beta <- runif(K, -0.5, +0.5)
R> y <- 1*(runif(N) < 1/(1+exp(-X
\end{Sinput}
\end{Schunk}
We now initialize $\boldsymbol\beta$ with zeros, draw samples from \code{logit.f}, and save the sampled coefficients to \code{beta.smp}:
\begin{Schunk}
\begin{Sinput}
R> nsmp <- 100
R> beta.ini <- rep(0.0, K)
R> beta.smp <- array(NA, dim=c(nsmp,K))
R> for (i in 1:nsmp) {
+    beta.ini <- MfU.Sample(beta.ini, f=logit.f, uni.sampler="slice", X=X, y=y)
+    beta.smp[i,] <- beta.ini
+  }
\end{Sinput}
\end{Schunk}
We can compare the mean of the sampled coefficients (\code{beta.mcmc}) with the MLE estimate (\code{beta.glm}), throwing away the first half of the samples for burn-in:
\begin{Schunk}
\begin{Sinput}
R> beta.mcmc <- colMeans(beta.smp[(nsmp/2+1):nsmp,])
R> beta.glm <- glm(y~X-1, family="binomial")$coefficients
R> cbind(beta.glm,beta.mcmc)
\end{Sinput}
\begin{Soutput}
     beta.glm  beta.mcmc
X1 -0.4265100 -0.4013649
X2 -0.5416799 -0.5814078
X3 -0.1684548 -0.1717733
X4 -0.4056165 -0.4323484
X5  0.4914152  0.5592871
\end{Soutput}
\end{Schunk}
Increasing \code{nsmp} should bring \code{beta.mcmc} closer to \code{beta.glm}.

\subsection[]{Adaptive Rejection Sampling using \pkg{MfUSampler}}\label{subsection-usage-ars}
In order to perform ARS using \pkg{MfUSampler}, we must construct a function that outputs either the log-density or its gradient vector, depending on the value of the boolean argument \code{grad}:
\begin{Schunk}
\begin{Sinput}
R> logit.fg <- function(beta, X, y, mu=0.0, sigma=1e+6, grad) {
+    Xbeta <- X
+    if (grad) return (t(X)
+                      - (beta - mu)/(2*sigma^2))
+    return (logit.f(beta, X, y, mu, sigma))
+  }
\end{Sinput}
\end{Schunk}
We can now apply \code{MfU.Sample}, this time setting \code{uni.sampler} argument to \code{"ars"}:
\begin{Schunk}
\begin{Sinput}
R> beta.ini <- rep(0.0, K)
R> beta.smp <- array(NA, dim=c(nsmp,K))
R> for (i in 1:nsmp) {
+    beta.ini <- MfU.Sample(beta.ini, f=logit.fg, uni.sampler="ars", X=X, y=y)
+    beta.smp[i,] <- beta.ini
+  }
R> beta.mcmc <- colMeans(beta.smp[(nsmp/2+1):nsmp,])
R> cbind(beta.glm,beta.mcmc)
\end{Sinput}
\begin{Soutput}
     beta.glm  beta.mcmc
X1 -0.4265100 -0.4118904
X2 -0.5416799 -0.4960775
X3 -0.1684548 -0.1691510
X4 -0.4056165 -0.4288610
X5  0.4914152  0.4751731
\end{Soutput}
\end{Schunk}
Again, larger values of \code{nsmp} would lead to closer match between \code{beta.mcmc} and \code{beta.glm}.

So far, we ignored the control parameters of Slice sampler and ARS, thus falling back to the default values provided for each algorithm. Sometimes it is necessary to override the default values, which is possible by calling \code{MfU.Control}. To see a description of the control parameters, type \code{?MfU.Control} in your \proglang{R} session. More details on each algorithm and its tuning parameters can be found by consulting the paper~\citep{neal2003slice} for the Slice sampler, and the software documentation~\citep{rodrigues2014ars} as well as the paper~\citep{gilks1992adaptive} for ARS.

The Bayesian logistic regression example is somewhat contrived. For example, the log-posterior of Equation~\ref{equation-logprob-logistic} can be shown to have a negative-definite Hessian matrix and therefore log-concave, thus being eligible for the multivariate Stochastic Newton Sampling~\citep{mahani2014sns}. The real power of \pkg{MfUSampler}, however, comes from its application to novel Bayesian models where developers seek to rapidly prototype flexible DAGs to validate their model specification, while having a path towards high-peformance applications. Section~\ref{section-workflow} contains one such example.

\section[]{\pkg{MfUSampler} for Bayesian Prototyping}\label{section-workflow}
With \pkg{MfUSampler}, prototyping a novel Bayesian model can be quite fast. Encapsulation of Gibbs cycles inside the \code{MfU.Sample} function circumvents subtle bugs during implementation of Gibbs sampling, and allows the researcher to focus his/her mental power on other, more important complexities of the problem such as model specification, MCMC diagnostics, etc.

\subsection{Example 2: Heteroscedastic Linear Regression}\label{subsection-example-heteroscedastic}
Consider a linear regression problem where we suspect heteroscedasticity of regression residuals, and that the non-uniform residual variance is itself dependent on the same covariates that are used to explain the mean response. To model such behavior, we assume the $N$ response measurements, $y_n$,  are independently drawn from this normal distribution:
\begin{equation} \label{equation-regression-heteroscedastic}
y_n \sim \mathcal{N}(\mathbf{x}_n^t \boldsymbol\beta, \sigma_{\mathrm{max}}^2/(1+\exp(-\mathbf{x}_n^t \boldsymbol\gamma)))
\end{equation}
This specification digresses from ordinary linear regression by making the variance parameter data-dependent. The nonlinear transformation ensures that the variance for all points is bounded between $0$ and $\sigma_{\mathrm{max}}^2$. Given the matrix of covariates $\mathbf{X}$ and the response vector $\mathbf{y}$, our goal is to estimate $\boldsymbol\beta$, $\boldsymbol\gamma$ and $\sigma_{\mathrm{max}}$. There can be several reasons for our interest in a Bayesian treatment of this problem, and in drawing samples from the joint posterior distribution on $(\boldsymbol\beta, \boldsymbol\gamma, \sigma_{\mathrm{max}})$: 1) future extension of the model in a probabilistic framework, e.g. using a Hierrarchical Bayesian approach to pool data across heterogneous observation units, 2) our suspicion that the posterior distribution is not globally convex, has a complex shape with multiple possible maxima, and thus sampling the its entire landscape can be safer than looking for a single, local optimum, and 3) fully probabilistic treatment of parameter estimation as well as response prediction.

To avoid clutter, we assume non-informative priors on all parameters and ignore them to focus only on the log-likelihood function corresponding to Equation~\ref{equation-regression-heteroscedastic}, which can be implemented using this \proglang{R} function:
\begin{Schunk}
\begin{Sinput}
R> loglike <- function(beta, gamma, sigmamax, X, y) {
+    mean.vec <-
+    sd.vec <- sigmamax/sqrt(1+exp(-
+    return (sum(dnorm(y, mean.vec, sd.vec, log=TRUE)))
+  }
\end{Sinput}
\end{Schunk}
To feed the log-likelihood function into \code{MfU.Sample}, we write a thin wrapper around it:
\begin{Schunk}
\begin{Sinput}
R> loglike.wrapper <- function(coeff, X, y) {
+    K <- ncol(X)
+    beta <- coeff[1:K]
+    gamma <- coeff[K+1:K]
+    sigmamax <- coeff[2*K+1]
+    return (loglike(beta, gamma, sigmamax, X, y))
+  }
\end{Sinput}
\end{Schunk}
We can now generate some data and draw samples from the resulting log-likelihood. Note the use of \code{MfU.Control} to set a lower bound of \code{1e-3} on \code{sigmax}:
\begin{Schunk}
\begin{Sinput}
R> # generate simulated data
R> K <- 5
R> N <- 1000
R> X <- matrix(runif(N*K, -0.5, +0.5), ncol=K)
R> beta <- runif(K, -0.5, +0.5)
R> gamma <- runif(K, -0.5, +0.5)
R> sigmamax <- 0.75
R> mu <-
R> var <- sigmamax^2/(1+exp(-
R> y <- rnorm(N, mu, sqrt(var))
R> # initialize and sample
R> coeff <- c(rep(0.0, 2*K), 0.5)
R> mycontrol <- MfU.Control(n = 2*K+1, slice.lower = c(rep(-Inf,2*K), 0.001))
R> coeff.smp <- array(NA, dim=c(nsmp, 2*K+1))
R> t <- proc.time()[3]
R> for (i in 1:nsmp) {
+    coeff <- MfU.Sample(coeff, f=loglike.wrapper, X=X, y=y, control = mycontrol)
+    coeff.smp[i,] <- coeff
+  }
R> t <- proc.time()[3]-t
R> cat("time:", t, "\n")
\end{Sinput}
\begin{Soutput}
time: 1.718
\end{Soutput}
\begin{Sinput}
R> beta.est <- colMeans(coeff.smp[(nsmp/2+1):nsmp, 1:K])
R> gamma.est <- colMeans(coeff.smp[(nsmp/2+1):nsmp, K+1:K])
R> sigmamax.est <- mean(coeff.smp[(nsmp/2+1):nsmp, 2*K+1])
R> cbind(beta, beta.est, gamma, gamma.est)
\end{Sinput}
\begin{Soutput}
            beta    beta.est     gamma  gamma.est
[1,] -0.08169418  0.01082307 0.3845259  0.3442652
[2,] -0.36962596 -0.33175088 0.4858229  0.4563613
[3,]  0.46721080  0.42836681 0.2954812  0.1705552
[4,] -0.17247123 -0.25379361 0.3079317  0.3556921
[5,]  0.33722306  0.23397001 0.1617440 -0.3678595
\end{Soutput}
\begin{Sinput}
R> c(sigmamax, sigmamax.est)
\end{Sinput}
\begin{Soutput}
[1] 0.750000 0.753041
\end{Soutput}
\end{Schunk}
Note that we applied \code{MfU.Sample} to \code{sigmamax} even though it is a scalar rather than a vector. In principle, we could directly call the slice sampler, but using the same higher level function \code{MfU.Sample} call keeps the code more organized and easier to track. Despite \code{nsmp} being small, we see general agreement between actual and estimated parameters, especially for \code{sigsq} and \code{beta}. Of course, proper MCMC diagnostics including trace plot examination, histogram examination, and effective size calculation must be done to adjust sampling parameters and determine the next step in the modeling process, including model re-specification.

\subsection[]{Performance Improvement Techniques}\label{subsection-workflow-efficiency}
Feeding the full joint distribution for a DAG into \code{MfU.Sample} is often a good, first step but it can be computationally brute-force. There are several opportunities for performance improvement once the foundation is laid and the model structure is somewhat validated. The root-cause of inefficiency in our brute-force approach is that we are evaluating the full, joint density during univariate sampling of each coordinate of the state space. This may not be necessary for several reasons:
\begin{enumerate}
\item For some variables or variable groups, exact sampling techniques may be possible. This often arises when the likelihood and prior terms in a hierarchical model are conjugate, making the posterior distribution have the same functional form as the likelihood. When conjugacy allows for exact sampling, it is often the preferred route compared to MCMC sampling~\citep{robert1999monte}.
\item If some variables (or groups of variables) are conditionally-independent, their joint distribution, conditioned on the remaining variables, is separable~\citep{wilkinson2006parallel}. Therefore, while sampling each variable, we only need to evaluate a subset of additive terms comprising log-likelihood. This situation can arise, for example, in Hierarchical Bayesian regression models~\citep{rossi2003bayesian,gelman2006data}. In addition to permitting lighter conditional posterior evaluations, conditional independence can also be taken advantage of in parallel Gibbs sampling~\citep{wilkinson2006parallel}.
\item Even when conditional independence does not exist, we may still find opportunities to drop some of the additive terms in the log-posterior for all or a subset of the variables, thereby reducing the time needed to evaluate conditional posteriors during Gibbs cycles.
\end{enumerate}

In the above cases, the general strategy is to split the full state space into subspaces and apply the more efficient techniques within each subspace. Bayesian compilers are effective to varying degrees at identifying and taking advantage of such optimization opportunities. A detailed discussion of these topics is beyond the scope of this paper. Here we illustrate the last item in the above list, continuing the heteroscedastic linear regression example of Section~\ref{subsection-example-heteroscedastic}.

The \code{loglike} function of Section~\ref{subsection-example-heteroscedastic} can be expanded as below (ignoring constant terms):
\begin{eqnarray}\label{equation-heteroscedastic-terms}\nonumber
L(\boldsymbol\beta, \boldsymbol\gamma, \sigma_{\mathrm{max}}) &=& -N \mathrm{log} \sigma_{\mathrm{max}} + \frac{1}{2} \sum_{n=1}^N \mathrm{log} \left( 1 + \exp(- \mathbf{x}_n^t \boldsymbol\gamma) \right) \\
&& - \frac{1}{2 \sigma_{\mathrm{max}}^2} \sum_{n=1}^N \left\{ (y_n - \mathbf{x}_n^t \boldsymbol\beta)^2 \left( 1 + \exp (- \mathbf{x}_n^t \boldsymbol\gamma) \right) \right\}
\end{eqnarray}
We see that, of the above three additive terms, only the last term depends on all three variable blocks $\boldsymbol\beta$, $\boldsymbol\gamma$ and $\sigma_{\mathrm{max}}$, while the first two depend only on $\sigma_{\mathrm{max}}$ and $\boldsymbol\gamma$, respectively. Taking advantage of this, we can write simplified conditional log-likelihood functions as below:
\begin{align}
\left\{
\begin{array}{lll}
L(\boldsymbol\beta \,|\, \boldsymbol\gamma, \sigma_{\mathrm{max}}) &= - \frac{1}{2 \sigma_{\mathrm{max}}^2} \sum_{n=1}^N \left\{ (y_n - \mathbf{x}_n^t \boldsymbol\beta)^2 \left( 1 + \exp (- \mathbf{x}_n^t \boldsymbol\gamma) \right) \right\} \\ \\
L(\boldsymbol\gamma \,|\, \boldsymbol\beta, \sigma_{\mathrm{max}}) &=  \frac{1}{2} \sum_{n=1}^N \mathrm{log} \left( 1 + \exp(- \mathbf{x}_n^t \boldsymbol\gamma) \right) - \frac{1}{2 \sigma_{\mathrm{max}}^2} \sum_{n=1}^N \left\{ (y_n - \mathbf{x}_n^t \boldsymbol\beta)^2 \left( 1 + \exp (- \mathbf{x}_n^t \boldsymbol\gamma) \right) \right\} \\ \\
L(\sigma_{\mathrm{max}} \,|\, \boldsymbol\beta, \boldsymbol\gamma) &= -N \mathrm{log} \sigma_{\mathrm{max}} - \frac{1}{2 \sigma_{\mathrm{max}}^2} \sum_{n=1}^N \left\{ (y_n - \mathbf{x}_n^t \boldsymbol\beta)^2 \left( 1 + \exp (- \mathbf{x}_n^t \boldsymbol\gamma) \right) \right\}
\end{array}
\right.
\end{align}
Referring to the three terms on the right-hand side of Equation~\ref{equation-heteroscedastic-terms} as components 1-3, the \proglang{R} implementation of the above conditional distributions will be:
\begin{Schunk}
\begin{Sinput}
R> loglike.component1 <- function(sigmamax, N) -N*log(sigmamax)
R> loglike.component2 <- function(gamma, X) 0.5*sum(log(1 + exp(-
R> loglike.component3 <- function(beta, gamma, sigmamax, X, y) {
+    -sum((y-
+  }
R> loglike.beta <- function(beta, gamma, sigmamax, X, y) {
+    loglike.component3(beta, gamma, sigmamax, X, y)
+  }
R> loglike.gamma <- function(gamma, beta, sigmamax, X, y) {
+    loglike.component2(gamma, X) +
+      loglike.component3(beta, gamma, sigmamax, X, y)
+  }
R> loglike.sigmamax <- function(sigmamax, beta, gamma, sigma, X, y) {
+    loglike.component1(sigmamax, nrow(X)) +
+      loglike.component3(beta, gamma, sigmamax, X, y)
+  }
\end{Sinput}
\end{Schunk}
Each Gibbs cycle will be broken into 3 steps, corresponding to $\boldsymbol\beta$, $\boldsymbol\gamma$ and $\sigma_{\mathrm{max}}$:
\begin{Schunk}
\begin{Sinput}
R> beta.ini <- rep(0.0, K)
R> gamma.ini <- rep(0.0, K)
R> sigmamax.ini <- 0.5
R> mycontrol.sigmamax <- MfU.Control(n = 1, slice.lower = 0.001)
R> coeff.smp <- array(NA, dim=c(nsmp, 2*K+1))
R> t <- proc.time()[3]
R> for (i in 1:nsmp) {
+    beta.ini <- MfU.Sample(beta.ini, loglike.beta, gamma=gamma.ini
+                   , sigmamax=sigmamax.ini, X=X, y=y)
+    gamma.ini <- MfU.Sample(gamma, loglike.gamma, beta=beta.ini
+                   , sigmamax=sigmamax.ini, X=X, y=y)
+    sigmamax.ini <- MfU.Sample(sigmamax, loglike.sigmamax
+                   , beta=beta.ini, gamma=gamma.ini
+                   , X=X, y=y, control = mycontrol.sigmamax)
+    coeff.smp[i,] <- c(beta.ini, gamma.ini, sigmamax.ini)
+  }
R> t <- proc.time()[3]-t
R> cat("time:", t, "\n")
\end{Sinput}
\begin{Soutput}
time: 1.589
\end{Soutput}
\begin{Sinput}
R> beta.est <- colMeans(coeff.smp[(nsmp/2+1):nsmp, 1:K])
R> gamma.est <- colMeans(coeff.smp[(nsmp/2+1):nsmp, K+1:K])
R> sigmamax.est <- mean(coeff.smp[(nsmp/2+1):nsmp, 2*K+1])
R> cbind(beta, beta.est, gamma, gamma.est)
\end{Sinput}
\begin{Soutput}
            beta      beta.est     gamma  gamma.est
[1,] -0.08169418  0.0009252645 0.3845259  0.3583743
[2,] -0.36962596 -0.3330523615 0.4858229  0.4546368
[3,]  0.46721080  0.4348825919 0.2954812  0.0726346
[4,] -0.17247123 -0.2702465415 0.3079317  0.2755692
[5,]  0.33722306  0.2422290123 0.1617440 -0.3499572
\end{Soutput}
\begin{Sinput}
R> c(sigmamax, sigmamax.est)
\end{Sinput}
\begin{Soutput}
[1] 0.7500000 0.7532132
\end{Soutput}
\end{Schunk}
In this case, the time savings from our improvements is modest, but for other problems the impact can be more pronounced. For example, for HB regression problems the speedup from breaking down the conditional posterior across regression groups will roughly be proportional to the number of groups, even before applying any parallelization (which could theoretically offer another multiplicative speedup equal to number of groups, given sufficient number of processing cores available).

For statistical problems of moderate to large size (i.e. number of observations and/or covariates) and in the absence of conjugacy for coefficients that are directly involved in explaining the response, the majority of time is often spent in log-density evaluations, rather than other activities such as random number generation or the sampling algorithm itself. Therefore, the next most rewarding optimization step is likely to be porting of log-density functions to high-performance languages such as \proglang{C}, \proglang{C++}, \proglang{FORTRAN}. For large problems, even parallel hardware such as Graphic Processing Units (GPUs) can be utilized by writing log-density functions in languages such as \proglang{CUDA}, while continuing to take advantage of \pkg{MfUSampler} for sampler control logic. Minimizing data movement between processor and co-processor is a key performance factor in such cases. Finally, should further performance improvement necessitate a rewrite of the \pkg{MfUSampler} logic in a high-peformance language, the package source code can be used as a blue-print for efficient development.

\section{Summary}\label{section-summary}
The \proglang{R} package \pkg{MfUSampler} enables MCMC sampling of multivariate distributions using univariate algorithms. It relies on an extension of Gibbs sampling from univariate independent samplig to univariate Markov transitions, and proportionality of conditional and joint distributions. By encapsulating these two concepts in a library, it reduces the possibility of subtle mistakes by researchers while re-implementing the Gibbs sampler and thus allows them to focus on other, more innovative aspects of their Bayesian modeling. Brute-force application of \pkg{MfUSampler} allows researchers to get their project off the ground, maintain full control over model specification, and utilize robust univariate samplers. This can be followed by an incremental optimization approach by taking advantage of DAG properties such as conjugacy, conditional independence and by porting log-density functions to high-peformance languages and hardware.


text:
\section{Languages of South Africa}
\label{sec:african}

The majority of languages in Africa are considered under-resourced~\cite{orife+etal_ICLR20}. %
This includes the eleven
official languages of South Africa.
As show in Figure~\ref{fig:hierarchy}, nine of these languages %
belong to the larger Southern Bantu family: isiZulu (Zul), isiXhosa (Xho), Sepedi (Nso), Setswana (Tsn), Sesotho (Sot), Xitsonga (Tso), siSwati (Ssw), Tshivenda (Ven) %
and isiNdebele (Nbl).
Many of these languages are also spoken %
in countries neighbouring South Africa.
From the Southern Bantu family there exist two principal, %
Nguni-Tsonga and Sotho-Makua-Venda. %
In the Sotho-Makua-Venda subfamily, Ven is somewhat of a standalone.
The other two languages, Afrikaans (Afr) and English (Eng), are Germanic languages from the %
Indo-European %
family.
All of these languages are considered under-resourced, except for Eng~\cite{barnard+etal_sltu14}.
To give an intuitive idea of how related these languages are, most of the Bantu languages  that are grouped together at the lowest level of the hierarchy %
would (to an %
extent) be intelligible to a native speaker of another languages in the same group.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.99\linewidth]{figures/language_hierarchy.pdf}
    \captionsep
    \caption[Caption]{A family tree for the official South African languages.\footnotemark
    }

    \label{fig:hierarchy}
    \vspace{-5mm}
\end{figure}
\footnotetext{\scriptsize{\url{https://southafrica-info.com/arts-culture/11-languages-south-africa/}}}

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

We investigated the effect of training language choice when applying a multilingual acoustic word embedding model to a zero-resource language.
Using word discrimination and query-by-example search tasks on languages spoken in Southern Africa, we
showed that training a multilingual model on languages related to the target is beneficial.
We observed gains in absolute scores, but also in data efficiency: you can achieve similar performance with much less data when training on multiple languages from the same family as the target.
We showed that even including just one related language
already gives a large
gain.
From a practical perspective, these results indicate that one should prioritise collecting data from related languages (even in modest quantities) rather than collecting more extensive datasets from diverse unrelated families, when building multilingual acoustic word embedding models for a zero-resource language.

\vspace{2pt}
{\eightpt
\noindent \textbf{Acknowledgements.}
This work is supported by the South African NRF (120409), a Google Faculty Award, and support from the Stellenbosch University School of Data Science and Computational Thinking.
We thank Ewald van der Westhuizen for the NCHLT forced alignments.}

\section{Experimental setup}
\label{sec:experiments}

\textbf{Data.} We perform all our experiments on the NCHLT speech corpus~\cite{barnard+etal_sltu14}, which provides wide-band speech from each of the eleven official South African languages (\S\ref{sec:african}).
We use a version of the corpus
where all repeated utterances were removed, leaving
roughly 56 hours of speech from around 200 speakers in each language.
We use the default training, validation and test sets.
We treat six of the languages as well-resourced:
Xho$^*$, Ssw$^*$, Nbl$^*$, Nso$^\dagger$, Tsn$^\dagger$ and Eng$^\ddagger$.
We use labelled data from the well-resourced languages to train a single supervised multilingual model and then apply the model to the target zero-resource languages.
More specifically,
we extract 100k true positive word pairs using forced alignments for each training language.
We select Zul$^*$, Sot$^\dagger$, Afr$^\ddagger$, and Tso as our zero-resource languages and use another language, Ven, for validation of each model.\footnote{We use superscripts to indicate the different language families: $^*$Nguni, $^\dagger$Sotho-Tswana, $^\ddagger$Germanic.}
Training and evaluation languages are
carefully selected such that for each evaluation language at least one language from the same
family is part of the training languages, except for Tso which is in a group of its own.

\label{subsec:awe_models}

\textbf{Models.}
All speech audio is parametrised as
13-dimensional
static Mel-frequency cepstral coefficients.
The encoder unit of the \tablesystem{ContrastiveRNN} models (\S\ref{sec:model}) consists of three unidirectional RNNs with 400-dimensional hidden vectors, with an embedding size of $M = 130$
dimensions.
Models are optimised using Adam optimisation~\cite{kingma+ba_iclr15} with a learning rate of $0.001$.
The temperature parameter $\tau$ in~\eqref{eqn:contrastive_loss} is set to 0.1.

\textbf{Evaluation.}
We consider two tasks for evaluating the performance of the AWE models. First, we
use the \textit{same-different}
word discrimination task~\cite{carlin+etal_icassp11} to measure the intrinsic quality of the
AWEs.
To evaluate a particular AWE model, a set of isolated test word segments is embedded.
We use roughly 7k isolated word instances per language from the test data.
For every word pair in this set, the cosine distance between their embeddings is calculated.
Two words can then be classified as being of the same or different type based on some distance threshold, and a precision-recall curve is obtained by varying the threshold.
The area under this curve is used as final evaluation metric, referred to as the average precision (AP).
We are particularly interested in obtaining embeddings that are speaker invariant.
As in~\cite{hermann+etal_csl21},
we therefore calculate AP by only taking the recall over instances of the same word spoken by different speakers.

The second task is QbE (\S\ref{sec:qbe}).
For each evaluation language we use approximately two hours of test
utterances as the search
collection.
Sub-segments for the utterances in the speech collection are obtained by embedding
windows stretching from 20 to 60 frames with a 3-frame overlap.
For each evaluation language we randomly
draw instances of 15 spoken query word types from a disjoint speech set (the development set---which we never use for any validation experiments) where we only consider query words with
at least 5 characters
for Afr and Zul and 3 for Sot.
There are
between 6 and 51 occurrences of each query word.
For each QbE test, we ensure that the relevant
multilingual AWE model
has not seen any of the search or query data
during training or validation.
We report precision at ten
($P@10$),
which is the fraction of the ten top-scoring retrieved utterances from the search collection that
 contains the given query.

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

Developing robust speech systems for zero-resource languages---where no transcribed speech resources are available for model training---remains a challenge.
Although full speech recognition is not possible in most zero-resource settings, researchers have proposed methods for applications such as speech search~\cite{levin+etal_icassp15,huang+etal_arxiv18,yuan+etal_interspeech18},
word discovery~\cite{park+glass_taslp08,jansen+vandurme_asru11,ondel+etal_interspeech19,rasanen+blandon_arxiv20}, and segmentation and clustering~\cite{kamper+etal_asru17,seshadri+rasanen_spl19,kreuk+etal_interspeech20}, making use of only
unlabelled speech audio.
Many of these applications require speech segments of different lengths to be compared.
This is conventionally done using alignment (e.g.\ with dynamic time warping).
But this can be slow and inaccurate.

Acoustic word embedding (AWE) models map a
variable duration speech segment
to a fixed dimensional vector~\cite{levin+etal_asru13}.
The goal is to map instances of the same word type to similar vectors.
Segments can then be efficiently compared by calculating the distance
in the embedding space.
Given the advantages AWEs have over
alignment methods, several AWE models have
been proposed~\cite{bengio+heigold_interspeech14,he+etal_iclr17,audhkhasi+etal_stsp17,wang+etal_icassp18,chen+etal_slt18,holzenberger+etal_interspeech18,chung+glass_interspeech18,haque+etal_icassp19,shi+etal_slt21,palaskar+etal_icassp19,settle+etal_icassp19,jung+etal_asru19}. Many of these are for the supervised setting, using labelled data to train a discriminative model.

For the zero-resource setting, a number of unsupervised AWE approaches have also been explored, many relying on autoencoder-based neural models trained on unlabelled data in the target language~\cite{chung+etal_interspeech16,kamper+etal_icassp16,holzenberger+etal_interspeech18,kamper_icassp19}.
However, there still exists a large performance gap between these unsupervised models and
their supervised counterparts~\cite{levin+etal_asru13,kamper_icassp19}.
A recent alternative for obtaining AWEs on a zero-resource language is to use multilingual transfer learning~\cite{ma+etal_arxiv20,kamper+etal_icassp20,kamper+etal_taslp21,hu+etal_slt21, hu+etal_interspeech20}.
The goal is
to have the benefits of supervised learning by training a model on labelled data from multiple well-resourced languages, but to then
 apply
the model to an unseen target zero-resource language without fine-tuning it---a form of \textit{transductive} transfer learning~\cite{ruder_phd19}.
This multilingual transfer approach has been shown to outperform unsupervised monolingual AWE models~\cite{kamper+etal_taslp21,jacobs+etal_slt21}.

Although
there is clear benefit in applying multilingual AWE models
to an unseen zero-resource language, it is still unclear how the particular choice of training languages affects subsequent performance.
Preliminary experiments~\cite{kamper+etal_taslp21} show improved scores
when training a monolingual model on one language and applying it to another from the same family.
But this has not been investigated systematically and there are still several unanswered questions:
Does the benefit of training on related languages diminish as we train on more languages (which might or might not come from the same family as the target zero-resource language)?
When training exclusively on related languages, does performance suffer when adding an unrelated language?
Should we prioritise data set size or language diversity when collecting data for multilingual AWE transfer?

We try to answer these questions using
a corpus of under-resourced languages spoken in Southern Africa.
These languages can be grouped into different families
based on their linguistic links.
We specifically
want to see whether it is beneficial to use closely related languages when training
a multilingual model catered to a specific zero-resource language, similar to~\cite{yi+etal_taslp19}.
We divide the corpus into training and test languages, with (some of) the test languages coming from families that also occur in
training.
We conduct several experiments where we add data from different language families, and also
control for the amount of data per language.
AWEs are evaluated
in an isolated word discrimination task and in query-by-example (QbE) speech search on full utterances.
To our knowledge, only one other
study~\cite{hu+etal_slt21} has done
AWE-based QbE
using multilingual transfer.

Our main findings are as follows.
(i)~Training a multilingual model using languages that are closely related to the target
improves
performance.
This is true, not just
because of the increase in data, but because of
language diversity.
(ii)~When
we systematically add training languages, the largest improvements are
gained from adding a single related language. Adding more related languages
gives small gains.
Adding unrelated languages generally gives small or no
gains, but also doesn't hurt.
(iii)~For a
target language, the
performance of a multilingual model trained on unrelated languages can be matched by a model trained with much less data from multiple related~languages.

\section{Acoustic word embedding model}
\label{sec:model}

We use the \tablesystem{ContrastiveRNN} AWE model of~\cite{jacobs+etal_slt21}.\footnote{We extend the code available at \scriptsize{\url{https://github.com/christiaanjacobs/globalphone_awe_pytorch}}.}
It performed the best of the model variants considered for multilingual transfer in~\cite{jacobs+etal_slt21}. %
The model consists of an encoder recurrent neural network (RNN) %
that produces fixed-dimensional representations from variable-length speech segments.
It is trained %
to minimise %
the distance between embeddings from %
speech segments %
of the same word type while %
maximising the distance between embeddings from multiple words of a different type.
Let's use $X = \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_T$ to denote a sequence of speech features.
Formally, given %
speech segments
$X_a$ and $X_p$ containing instances of the same word type and multiple negative examples $X_{n_{1}}, \ldots, X_{n_{K}}$, the \tablesystem{ContrastiveRNN} produces embeddings $\mathbf{z}_a, \mathbf{z}_p, \mathbf{z}_{n_{1}}, \ldots, \mathbf{z}_{n_{K}}$
(subscripts indicate anchor, positive and negative, respectively).
Each embedding is a fixed dimensional vector $\mathbf{z} \in \mathbb{R}^M$.
The model is illustrated in Figure~\ref{fig:contrastiveRNN}.
Let $\text{sim}(\mathbf{u}, \mathbf{v}) = \mathbf{u}^{\top}\mathbf{v}/\norm{\mathbf{u}}\norm{\mathbf{v}}$ denote the cosine similarity between vectors $\mathbf{u}$ and $\mathbf{v}$.
The loss given a positive pair $(X_a, X_p)$ and the set of negative examples is then defined as~\cite{chen+etal_icml20}:
\vspace*{-2.5pt}
\begin{equation}
J = -\text{log}\frac{\text{exp}\big\{\text{sim}(\mathbf{z}_a, \mathbf{z}_p)/\tau\big\}}{\sum_{j \in \{p, n_1, \hdots, n_K\}}^{}\text{exp}\big\{\text{sim}(\mathbf{z}_a, \mathbf{z}_j)/\tau\big\}}\,\text{,}
\label{eqn:contrastive_loss}
\end{equation}
where $\tau$ is a temperature parameter, tuned on development data.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.99\linewidth]{figures/contrastive_rnn.pdf}
    \captionsep
	\caption{The \tablesystem{ContrastiveRNN}.
	The model is trained to minimise the distance between the anchor and positive item $d(\mathbf{z}_a, \mathbf{z}_p)$ while maximising the distance between the anchor and multiple negatives $(\mathbf{z}_a, \mathbf{z}_{n_{k}})$.}
	\label{fig:contrastiveRNN}
	\vspace{-5mm}
\end{figure}

In the zero-resource setting we don't have labelled data in the target {language} %
to construct the positive and negative %
word pairs required for training.
We therefore follow the approach of~\cite{kamper+etal_taslp21}, and train a multilingual model on ground truth word pairs (extracted from forced alignments) from a number of %
languages for which we have labelled data.
Subsequently, at test time, we apply the encoder RNN from the multilingual model to extract AWEs for speech from
the target zero-resource language.

\section{Query-by-example speech search}
\label{sec:qbe}

For evaluating the different AWE {models}
we use
an isolated word discrimination task~\cite{carlin+etal_icassp11}.
Recent findings~\cite{algayres+etal_interspeech20} suggest, however, that this evaluation is not always indicative of downstream system performance.
We therefore also perform query-by-example (QbE) speech search, which in contrast to
{the word discrimination task}, does not assume a test set of isolated words, but instead operates on full unsegmented utterances.

Concretely, QbE speech search is the task of identifying the utterances in a speech collection that contain instances of a given spoken query.
A number of approaches have been put forward for AWE-based QbE~\cite{levin+etal_icassp15, chen+etal_icassp15, settle+etal_interspeech17}.
Here we use the simplified approach from~\cite{kamper+etal_icassp19}.
Using an AWE model, we first embed the query segment.
If we knew the word boundaries in the search collection, we could embed each of the words in an utterance and simply look up the closest embeddings to the query.
Instead, because we do not have word boundaries, each utterance is split into overlapping segments from some minimum to some maximum duration.
Each segment from each utterance is then embedded separately using the AWE model.
Finally, to do the QbE task, the query embedding is compared to each of the utterance sub-segment embeddings (using cosine distance), and the minimum distance over the utterance is then taken as the score for whether the utterance contains the given query.

\section{Experimental results}
\label{sec:results}

\begin{figure}[!t]
	\centering
	\includegraphics[width=0.99\linewidth]{figures/cross_lingual_heatmap.pdf}
    \captionsep
	\caption{AP (\%)
    when training a monolingual supervised \tablesystem{ContrastiveRNN} on each language (rows) and then evaluating it on each of the other languages (columns). Heatmap colours are normalised for each evaluation language (i.e. per column).}
	\label{fig:heatmap}
	\vspace{-4mm}
\end{figure}

\subsection{Cross-lingual evaluation}

Before looking at multilingual modelling, we first consider a cross-lingual evaluation where we treat each language as a training language, train a supervised monolingual AWE model, and then apply it to every other language.
This allows us to see the effect of training on related languages in a pairwise fashion.
The results are shown in Figure~\ref{fig:heatmap}.
For each evaluation language excluding Tso and Ven, which are in family groups of their own,  %
the best results are achieved from models trained on a language from the same family.
E.g.\ on Zul, Xho is the best training language giving an AP of 58.5\%.
Eng is the only exception where the model trained on Sot performs better than using Afr, the other Germanic language.
Although Ven is in its own group at the lowest layer of the family tree in Figure~\ref{fig:hierarchy}, some of the best results when evaluating on Ven are obtained using models trained on Sotho-Tswana languages (Sot, Tsn, Nso), which are in the same family at a higher level.
We also see that for all nine Bantu evaluation languages, worst performance is obtained from the two Germanic models (Afr, Eng).

\begin{table}[!t]
	\mytable
	\caption{AP (\%) on test data for multilingual models trained on different combinations of well-resourced languages. Models are applied to two zero-resource languages from different language families, Nguni and Sotho-Tswana. For each training language 100k word pairs were extracted.}
    \captionsep
	\vspace{-1mm}
    \eightpt
	\begin{tabularx}{1\linewidth}{Lcc}
		\toprule
		Multilingual model & Zul$^*$ & Sot$^\dagger$ \\
		\midrule
		\underline{\textit{Nguni:}} & &\\[2pt]
		Xho$^*$ + Ssw$^*$ + Nbl$^*$ &\textbf{68.6} &---  \\
		Xho$^*$ + Ssw$^*$ + Eng$^\ddagger$ &60.9 &---  \\
		Xho$^*$ + Nso$^\dagger$ + Eng$^\ddagger$ &55.7 &---  \\
		Tsn$^*$ + Nso$^\dagger$ + Eng$^\ddagger$ &37.5 &---  \\
		Xho$^*$ + Ssw$^*$ + Nbl$^*$ {(\scriptsize{$10$\% subset})} &58.6 &---  \\[2pt]
		\underline{\textit{Sotho-Tswana:}} \\[2pt]
		Nso$^\dagger$ + Tsn$^\dagger $&--- &\textbf{76.7} \\
		Nso$^\dagger$ + Eng$^\ddagger$ &--- &64.8 \\
		Xho$^*$ + Ssw$^*$ &--- &51.9 \\
		Xho$^*$ + Eng$^\ddagger$ &--- &52.5 \\
		Nso$^\dagger$ + Tsn$^\dagger$ {(\scriptsize{$10$\% subset})} &--- &58.4 \\
		\bottomrule
	\end{tabularx}
	\label{tbl:multilingual}
	\vspace{-4mm}
\end{table}

\begin{figure*}[t]
	\centering
	\includegraphics[scale=0.31]{figures/multilingual_increment.pdf}
    \captionsep
	\caption{Same-different results %
    from two sequences of multilingual models, trained by adding one language at a time. For each training language, 100k positive word pairs are used, which is indicated on the $x$-axis.} %
	\label{fig:multi_increment}

    \vspace{-3.5mm}

\end{figure*}

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.31]{figures/qbe_zul.pdf}
    \captionsep
    \caption{QbE results %
    on %
    Zul %
    using the same sequences of multilingual models as in Figure~\ref{fig:multi_increment}.}
    \label{fig:qbe_increment}
    \vspace{-4mm}
\end{figure}

\subsection{Multilingual evaluation}

The cross-lingual experiment above was in large part our inspiration for the subsequent analysis of multilingual models.
Concretely, we hypothesise that even better performance can be achieved by training on multiple languages from the same family as the target zero-resource language, and that this would be superior to multilingual models trained on unrelated languages.
Focusing on two evaluation languages, Zul and Sot from distinct language families, we
investigate this hypothesis
by training multilingual models with different language combinations, as shown in Table~\ref{tbl:multilingual}.
Firstly, we see the best result on a %
language when all the training languages comes from the same family as the target.
Secondly, we see how the performance %
gradually decrease as the number of training languages related to the target drop.
Furthermore, notice the performance boost from %
including even just one training language related to the target compared to not including any.
E.g.\ on Sot we see a increase of more than 12\% absolute when adding just one related language (from 52.5\% and 51.9\% to 64.8\%). %

To further demonstrate the benefit of using training languages from the same family, %
we train a multilingual model for each evaluation %
language on all its related languages using a 10\% subset of the original data.
For both Zul and Sot, %
the subset models outperform the models where no related languages are used. %
E.g.\ on Zul, the Xho+Ssw+Nbl subset model outperforms the full Tsn+Nso+Eng model (no related languages) by more than 20\% in AP.
Moreover, this subset model (58.6\% in AP) even outperforms the Xho+Nso+Eng model (55.7\%) where all the training data from one related language are included and almost matches the AP when using two related languages (Xho+Ssw+Eng, 60.9\%).

These comparisons do more than just show the benefit of training on related languages: they also show that it is beneficial to train on a diverse set of related languages. %

\subsection{Adding more languages} %

In the above experiments, we controlled for the amount of data per language and saw that training on languages from the same family improves multilingual transfer.
But this raises a question:
will adding additional unrelated languages harm performance? %
To answer this, we systematically train two sequences of multilingual models on all six well-resourced languages, %
evaluating each target language as a new training language is added.

Same-different results for all five evaluation languages are shown in Figure~\ref{fig:multi_increment} and QbE results for Zul are shown in Figure~\ref{fig:qbe_increment}. As for Zul, the trends in the same-different and QbE results track each other closely for the other evaluation languages (so these are not shown here).
In the fist sequence of multilingual models (green), we start by adding the three Nguni languages (Xho, Ssw, Nbl), followed by the two Sotho-Tswana languages (Nso, Tsn), and lastly the Germanic training language (Eng).
The second sequence (orange) does not follow a systematic procedure.

On Zul, the green sequence, which starts with a related language (Xho), initially achieves a higher score compared to the orange sequence in both Figures~\ref{fig:multi_increment} and~\ref{fig:qbe_increment}.
Then, the score gradually increase by adding more related languages (Ssw, Nbl). %
Thereafter, adding additional unrelated languages (Nso, Tsn, Eng) show no performance increase.
In fact, AP decreases slightly after adding the two Sotho-Tswana languages (Nso, Tsn), but not significantly.
The orange sequence starts low on Zul until the first related language (Ssw) is added, causing a sudden  increase.
Adding the Germanic language (Eng) has little effect.
Adding the last two related languages (Nbl, Xho) again causes the score to increase.
A similar trend follows for Sot and Afr, where adding related languages causes a noticeable performance increase, especially when adding the first related language; after this, performance seem to plateau %
when adding more unrelated languages.
(Afr is the one exception, with a drop when adding the last language in the orange sequence).
On Tso, which does not have any languages from the same family in the training set, AP gradually increases in both sequences without any sudden jumps.
Although Ven isn't in the same family as Nso and Tsn at the lowest level of the tree in Figure~\ref{fig:hierarchy}, it belongs to the same family {(Sotho-Tswana)} at a higher level. %
This explains why it
closely tracks the Sot results.

Summarising these results, we see that adding unrelated languages generally does not decrease scores, but also does not provide a big benefit (except if it is one of the earlier languages in the training sequence, where data is still limited). In contrast, it seems that training on languages from the same family is again beneficial; this is especially the case for the first related language, irrespective of where it is added in the sequence.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Solar-type stars in the mass range of $\sim$ 1.0 to 1.5~\msol, which are
mostly of the spectral type F, are  uniquely suited for examining issues
of stellar evolution, evolution of stellar activity, and the history of star
formation in the solar neighborhood. Their particular suitability is
due to the large spread in their evolutionary status and the availability
of large amounts of varied, pertinent data, as summarized below:
\begin{itemize}
\item
F stars are the most numerous type in the {\it Hipparcos} catalog, greatly
outnumbering all other types of stars with  measured
trigonometric parallaxes.
\item
For well over ten thousand F~stars, {\it uvby} photometric parameters are
available, which yield star's surface gravity along with temperature,
luminosity, and metal abundance.
\item
Thousands of these stars have measured radial velocities, which in
combination with proper motion and parallax from {\it Hipparcos} yield
spatial velocities.
\item
Over two thousand of the same population have {\it Rosat} data, which
translate to X-ray luminosities when combined with {\it Hipparcos} parallaxes.
\item
Single F stars are well differentiated from binary stars.
\item
F stars on and above the zero-age main sequence (ZAMS) straddle the entire age
 range of stars in the Galactic disk. With known distances and photometry,
this means the availability of their age---from standard stellar evolution
modeling---in the range of 0 to $\sim\!10$~Gyr.
\end{itemize}

The abundance of data from {\hipparcos} (ESA 1997), {\rosat} (Voges et al.
1999, 2000) and  \uvby photometry (Hauck \& Mermillod  1998), with
derived parameters such as temperature and metallicity and computed parameters
such as age, provides a tremendously rich background to explore the evolution
of stars and stellar  activity of solar-type stars over a period stretching
from the present epoch back to the time when the first stars
were formed in the nearby Galactic disk. Earlier, we used the indicated data
to address some aspects of these issues
(Suchkov \& McMaster 1999, Suchkov 2000, 2001, Griffin \& Suchkov 2003,
and Suchkov,  Makarov, \& Voges 2003). More recently Holmberg, Nordstrom,
\& Andersen (2009) published
refined \hipparcos parallaxes along with new data and derived parameters for
many \hipparcos stars. These newly measured parameters were radial velocity and
 multiplicity (differentiation between single, binary, etc. stars). The derived
parameters included temperature, metallicity, spatial (3D) velocity,
and age. In this paper, we take advantage of that body of data to revisit
and expand on our previous work.

\section{Data}

The base sample used in this paper is built from a set of 11900
F~stars described in Suchkov \& McMaster (1999). The astrometry and stellar
parameters are taken from Holmberg et al. (2009).
Along with the parameter constraints used in Suchkov \& McMaster (1999),
the current sample is additionally constrained specifically with respect
to the parameters from Holmberg et al. (2009):
\begin{itemize}
\item
Parallax error: $\sigma_{\pi}/\pi < 0.15$.
\item
Distance: $d < 150$~pc.
\item
Spatial velocity: $v < 200$~\kmsec.
\item
Metallicity: $-1.0 < \ferrh < 0.5$.
\item
Multiplicity: no binary or triple stars; also
no open cluster members.
\end{itemize}
The overluminosity parameter discussed below is constrained by
$-0.15 \leq \dmczero \leq 1.0$.
The subsample  of X-ray emitters has X-ray luminosity from Suchkov,
Makarov, \& Voges (2003) and is limited to $26 < \loglx  < 33$.

\section{Overluminosity}

Suchkov \& McMaster (1999) introduced an {\it overluminosity parameter}
defined as the difference between absolute magnitude $M_{c_0}$ derived
from the \uvby photometry and absolute magnitude $M_V$ based on
\hipparcos parallax:
\begeq
\dmczero = M_{c_0} - M_V.
\eneq
For F stars, the \uvby parameter $c_0$ is, in fact, sensitive to surface
gravity rather than luminosity.
It is a usable measure of luminosity as long as there is a one-to-one
relationship
between surface gravity and luminosity. The reasonably tight relationship
between $c_0$ and $M_V$ found for well-studied nearby single F~stars
implies that for many stars this is indeed so. But in general, such
a relationship cannot be presumed, and $\Delta M_{c_0}$ differentiates
single  F~stars according to their luminosity at a given surface gravity,
effective temperature, and chemical composition.

By definition, single stars of normal
luminosity will have   $\Delta M_{c_0} =0$. Given the uncertainty in the
measured values of the overluminosity parameter, the actual
 operational criterion of normal stars was set to $\Delta M_{c_0} =0 \pm 0.15$,
 which will also be adopted in this paper.
{\it Overluminous stars} (the term introduced
by Griffin \& Suchkov 2003), i.e., stars whose luminosity
is higher than that of normal stars of the same  temperature, gravity,
and chemical composition, are then defined as  those with
$\Delta M_{c_0} > 0.15$. At a given luminosity and temperature an overluminous
star has a larger surface gravity; hence, of the two stars located at the
same point in the $\logtef$--$\logl$ diagram an overluminous star
is more massive than a normal one.

Overluminosity distribution is continuous, peaking
at $\dmczero \approx 0$ and extending up to $\dmczero \sim 1.5$
(see Figure~1 in Suchkov 2001). The splitting of
F~stars into two groups, normal and overluminous, is just a convenient
way to present the issues associated with the overluminosity phenomenon rather
than a division into two distinct, separate populations.
As argued below, overluminosity reflects peculiarities in stellar evolution
of  stars intrinsically different from normally evolving stars.
In the course  of stellar evolution, the \dmczero of stars
that become overluminous progressively increases at a rate that depends
on the original difference between these stars and normally evolving stars,
so at any given value of \dmczero there is a mixture of stars of all ages.
The way we split our sample into two groups ensures that a substantial
fraction of the overluminous group consists of far evolved stars.
This explains why the overluminous group looks so different from
the normal group in the $\log T_e - M_V$ diagram (Figure~\ref{f1}).

\begin{figure}[ht]
\epsscale{1.0}
\plotone{f1.eps}
\caption{$\log T_e - M_V$ diagram for normal and overluminous F stars.
The number of stars is shown in the upper right corner.
A set of isochrones (labeled by age in Gyr), based on Chaboyer at al. 1999
for composition $Y=0.27$ and $Z=2\times10^{-2}$, and a line one magnitude
above the ZAMS are shown for reference.
Overluminous stars tend to avoid the area near the ZAMS and
are on average brighter than normal stars at the same temperature.
}
\label{f1}
\end{figure}

Originally, overluminous stars that were listed as single in the \hipparcos
catalog were believed to be mostly unidentified binary
stars with comparably bright normal components (Suchkov \& McMaster 1999).
An unresolved binary with equally bright normal components
would look overluminous by $\dmczero = 0.75$.
The lingering problem  was that there were
stars with \dmczero significantly larger than 0.75, which definitely did
not fit to the binary star hypothesis. Another puzzle was that
overluminous stars as a group turned out to be  much older than the
known binaries in the \hipparcos (Suchkov 2000). So the binaries were not
the complete answer, and something else must be at work.
Indeed, later we found that overluminous stars were comprised of at least
three different subgroups (Griffin \& Suchkov 2003). Of those three one
subgroup were binaries not identified in the \hipparcos and the other subgroup
were young pre-main-sequence stars. The third group was comprised of truly
single main-sequence stars. The stars in this last category are the main focus
 of this paper.

\section{Overluminosity: Age--velocity relation }

\begin{figure}[ht]
\epsscale{1.0}
\plotone{f2.eps}
\caption{Age--velocity relation for F~stars using radial,
tangential or spatial velocity. The velocity statistics
is the mean of absolute velocity value. The sample is
limited to stars with $\Delta M_{c_0} > -0.15$.
}
\label{f2}
\end{figure}

Older stars are long known to have on average higher velocities, which is
 demonstrated by the age--velocity relation (AVR) of F~stars
shown in Figure~\ref{f2}.  The effect is usually
attributed to some kind of a dynamic process that gradually increases
stellar velocities as time goes on. Thanks to this clear-cut relationship
between age and stellar velocities, the kinematics parameters
such as velocity dispersion and the mean of absolute velocities are
commonly used as age markers of stellar populations. If calibrated against
actual (absolute) physical age (for instance, age based on stellar evolution
models), the resulting ''kinematics'' age
can be used as a substitute for the actual age of physical groups of stars,
which proves especially useful when direct age determination (more accurately,
age from stellar evolution modeling) is difficult or impossible.

\begin{figure}[ht]
\epsscale{1.0}
\plotone{f3.eps}
\caption{Whatever the data source of velocity---Doppler shift (radial velocity),
  parallax plus proper motion (tangential velocity), or the combination of both
 (spatial velocity)---there is a remarkable discrepancy in age-velocity
relations of normal (circles)  and overluminous (diamonds) F stars, with
age derived from isochrone fitting.
}
\label{f3}
\end{figure}

The age parameter in Figure~\ref{f2} is taken from
Holmberg, Nordstrom, \& Andersen 2009. It is based on  ``standard'' stellar
evolution modeling  that requires initial conditions for only
three parameters: stellar mass, $M$, helium content, $Y$, and heavy element
abundance, $Z$.  Standard models provide
star's luminosity and temperature as  functions of age, thus allowing
one to compute age from observationally derived luminosity, temperature,
and chemical composition.
This ''{\it standard}`` age is commonly used in isochrone fitting technique
to determine  age of   star clusters or individual stars from effective
temperature, \tef, luminosity, $L$, and chemical composition parameters $Z$
and $Y$, so we will also be referring to it as {\it isochrone} age.

 Depending on how adequate the underlying model is and how  accurate the
derived stellar parameters are, standard age may differ from the actual age of
a star. Therefore, if a body of observational data reveals a discrepancy
between the standard age and actual age, that would  mean that basic premises
of the standard stellar evolution models need to be reconsidered.

\section{Age and kinematics}

With kinematics depending only on actual (absolute) age, the AVR
of any group of stars must  be a unique function of {\it standard}
age if the latter adequately represents the actual age.
But what if a derived AVR turns out to be a non-unique function?
What if groups of stars differentiated according to some property have
very different AVR?
The most obvious answer would be that stellar evolution modeling uses
invalid assumptions, missing  some underlying physics  and/or  some
parameter(s) that are as important as stellar mass and chemical composition.
Stars differing in such unaccounted for parameter(s)
would produce  different isochrone-based AVRs as their standard age does
not represent  actual age.

A careful examination of available data demonstrates that F stars do indeed
reveal a property suggesting that there is some physical factor  missing
from the current standard evolution models, which
 affects stellar evolution in a major way.
The problem can be seen in Figure~\ref{f3} that shows
age--velocity relation separately for
normal  and overluminous stars. Instead of being identical, and duplicating
the AVRs in Figure~\ref{f2}, the AVRs for these two subsets
turn out to be quite different. The overluminous stars have
much higher velocities at the same isochrone age, the discrepancy
increasing with age. This means that the actual age of overluminous
stars is significantly underestimated with respect to normal stars.
For instance, the kinematics in Figure~\ref{f3} suggests that, at
isochrone age of 5~Gyr, overluminous stars are in fact $\sim 2.5$~Gyr
 older than normal stars, i.e., older by about 50\%. This inference agrees
 with our results obtained earlier from our original data set
 (Suchkov \& McMaster 1999). It is to be noted that the lower panel
 in Figure~\ref{f3} is virtually identical to Figure~7 in Suchkov (2001)
 even though, unlike  Figure~\ref{f3}, the data in the latter Figure
 are based on the original \hipparcos parallaxes.

\section{Overluminosity: X-ray -- velocity relation}

Stellar X-ray emission of late-type stars is a manifestation of stellar
coronal activity that is associated with the star's outer convection zone.
The convection zone is thinner
in stars of larger mass, disappearing in stars earlier than
spectral type F. This accounts for a stronger coronal activity
of low-mass stars in comparison with more massive stars of the same age.
Outer convection apparently disappears altogether in stars earlier than
$\sim$F5, which is argued to be the reason for the drop of coronal
activity in this spectral range (e.g., Stauffer et al. 1994).

Stellar X-ray emission is commonly believed to decline with age as coronal
activity decays.
The decrease of X-ray luminosity of main-sequence stars with age was
inferred back in the 1980's from the {\it Einstein} surveys of open clusters
and field stars (e.g., Micela et al. 1988). For main-sequence stars, coronal
activity is associated with  magnetic fields generated in the star's
outer convection zone. Presumably, the field generation mechanism
(magnetic dynamo) is supported by the combined effect of convection and
rotation. As the  star ages, rotation slows down and the dynamo
gets less efficient, which reduces the star's coronal activity and
drives its X-ray emission down.

\begin{figure}[ht]
\epsscale{1.0}
\plotone{f4.eps}
\caption{X-ray luminosity of normal F stars (diamonds) slightly
 declines with increasing velocities, hence with age. In sharp contrast,
 for overluminous stars (circles), X-ray luminosity trends up
rather than down, contrary to the conventional picture of
the stellar activity evolution.
}
\label{f-f4}
\end{figure}

Figure~\ref{f-f4} shows the relationship between spatial velocity and
X-ray luminosity of normal and overluminous F stars.
Intriguingly, while  the behavior of normal stars in this diagram is
consistent with the concept of decaying coronal activity, the behavior of
the overluminous stars does not. Overluminous stars are consistently and
significantly  more luminous at higher velocities than at low velocities,
which implies  that coronal activity of these stars increases rather than
decreases with age.  The same conclusion was reached in our earlier study
(Suchkov, Makarov, \& Voges 2003).

\section{Discussion}

What makes a solar-mass star live longer than the current standard models of
stellar evolution predict? What can make it get
more active as it ages? What do the standard stellar evolution and stellar
activity modeling miss?
It appears that the only important initial characteristic not included
in standard modeling, a characteristic that differentiates stars of
the same mass and chemical composition at the start of their evolution,
is differential  rotation of
the stellar interior. Differential rotation of the inner of a star
can cause radial mixing of the stellar material, transporting fresh
hydrogen fuel into the hydrogen-burning core and thus extending the
star's main-sequence lifetime. Stars with initially different rotational
velocities and velocity radial profiles would experience different amount
of hydrogen enrichment, so their main-sequence lifetimes would be different.

Over the last decade, the probing of internal rotation of stars has
become a major activity in the rapidly developing field of asteroseismology.
As a result, a significant radial variation of rotation in the stellar
interior, from the star's fast rotating core out to the much slower
rotating upper layers,
 appears now to be a rather well established feature that is typical
for many stars.
The respective studies are driven to a great degree by the efforts to
solve the remaining problems in stellar evolution.
For instance, Pamyatnykh, Handler, and Dziembowski (2004) addressed
the problem of the angular momentum evolution and element mixing in
convectively stable layers.  Using asteroseismology data for
$\nu$~Eri, they inferred that rotation in the star's
 $\mu$-gradient zone is three times faster than in the envelope.
 With this kind of velocity gradient one may expect a massive amount
 of rotationally induced element mixing in the inner layers
 of this star.

Stellar evolution modeling that
includes internal rotation as a new factor additional
to mass and chemical composition is becoming a popular field in stellar
evolution theory. Eggenberger,
Maeder, and Meynet (2010) studied the effects of rotational
mixing on the properties of solar-type stars. They concluded that
"Rotational mixing has indeed a large impact on the properties of
the central layers by bringing fresh hydrogen fuel to the central
stellar core" and  "Due to rotational mixing,
the main-sequence lifetime is then larger for stellar
models including rotation" (see also Eggenberger et al. 2010).
 Ekstr{\"o}m et al. (2012) found a considerable increase
in the main-sequence lifetimes for stellar masses between about 1~\msol and
2~\msol, by about 25\%, once rotation is included. The new stellar
evolution modeling, therefore, leaves little doubt that internal rotation
does indeed prolong main-sequence lifetimes. It seems safe to assume
that further progress in development of models featuring
rotational mixing will fully explain large  differences between the lifetimes
of normal and overluminous stars discussed in this paper.

Stars' differential rotation is quite likely responsible
 not only for longer lifetimes but
also for stronger X-ray emission of overluminous stars.
One may speculate that angular momentum transfer from
the rapidly rotating core and inner layers of a star upwards would
counteract the secular decay of rotation in the upper layers,
which would ensure
that the rotationally induced magnetic field remains strong and the associated
chromospheric activity, hence X-ray luminosity, remains high.
The stronger X-ray emission
of old overluminous stars implies, in fact, that any such mechanism
must be capable of not just maintaining but enhancing that emission
as the star ages.

\section{Conclusions}

The observational data that became available during the last two decades
 reveal that the population of solar-type stars in the mass range
$M \sim 1$~to~$1.5 \msol$,
represented mostly by F stars, is non-uniform
 in terms of stars' stellar evolution rate
and also in terms of the direction of the evolution of coronal activity.
First, F stars that have been found to be too luminous for their temperature
and surface gravity -- the {\it overluminous} stars -- are much
older than estimated from the
standard stellar evolution models. Second, the X-ray emission of old
{\it overluminous} stars
is on average stronger than that of young stars, which is
opposite to the expectations from the commonly accepted evolutionary path
of coronal activity for stars in the solar-mass range.

These two findings  suggest that standard models of
stellar evolution and evolution
of stellar activity overlook some crucial physical processes.
Recent developments in stellar evolution modeling that include inner
rotation and angular momentum transfer as a new factor additional to
mass and chemical composition seem to provide the necessary parameters
and mechanisms to explain the existence of overluminous stars and their
properties. This is good news. What is not so good is that
computing stellar age from comparison of model calculations and
observational data is now becoming a task much more challenging
than before.