\section{Introduction}
Urn models have provided simple and useful models to study different phenomena since the early works \cite{Polya, PolyaEggenberger}. The model starts with given quantities of objects of finitely many colors. A color is chosen at random, and further objects of each color are added to the urn, where the amount of each color added will depend on the color chosen. The additional amount, which will be nonnegative real valued, of each color can be summarized by a square matrix, called the replacement matrix, where the row corresponding to the color chosen will provide the amount to be added. The replacement matrices can be random, such that they are ``close'' to an irreducible replacement matrix in an appropriate sense. A square matrix $\boldsymbol{A}$ with nonnegative entries is called \textit{irreducible} if, for any $i, j$, there exists a positive integer $N\equiv N_{ij}$, such that the $(i,j)$-th entry of $\boldsymbol{A}^N$ is positive. By Perron-Frobenius theory, the \textit{dominant eigenvalue} of an irreducible matrix, namely the eigenvalue having the largest real part, is simple and positive and has corresponding left and right eigenvectors with all coordinates strictly positive.

In \cite{AthreyaNey}, the replacement matrices formed an i.i.d.\ sequence, with $L \log_+ L$ moments for the entries and irreducible mean matrix. Using branching process techniques, it was shown in \cite{AthreyaNey}, cf.\ Chapter~V.9.3, the normalized color count vector converges almost surely to the (normalized to probability) left eigenvector of the mean matrix corresponding to its dominant eigenvalue. However, the i.i.d.\ assumption on the replacement matrices can be restrictive for applications. The replacement matrices may depend on the entire past. In \cite{BaiHu} and \cite{LaruellePages}, the replacement matrices at any step are independent of the color choice in that step given the entire history. Martingale methods were used in \cite{BaiHu}, while stochastic approximation was used in \cite{LaruellePages}. In both works, the conditional means of the replacement matrices were assumed to concentrate around an irreducible, balanced (i.e., the row sums being same) matrix and almost sure convergence and central limit behavior of the proportion vector was obtained under second or higher moments. In \cite{Zhang12}, the balanced condition was removed and the mode of concentration of the conditional mean replacement matrices was relaxed. It then proved the almost sure convergence under $L (\log_+ L)^p$ moments, for some $p>1$, in general; and, under $L \log_+ L$ moments only for i.i.d.\ replacement matrices. Note that the latter result is essentially same as the one obtained in \cite{AthreyaNey}. Further, convergence in probability was proved in \cite{GangoMaulik} under $L^1$ moments alone.

In this article, we significantly relax the assumptions to obtain almost sure convergence results applying the method of stochastic approximation. Our method obtains almost sure convergence under $L (\log_+ L)^p$ condition for some $p>1$, as in \cite{Zhang12}. However, we relax the moment condition further; and use only $L \log_+ L$ condition with appropriate majorization condition on the norms of the replacement matrices. The additional moments in \cite{Zhang12} or the majorization condition provide the required uniform integrability. We also establish almost sure convergence result for i.i.d.\ replacement matrices under $L^1$ conditions alone, instead of $L \log_+ L$ condition in \cite{AthreyaNey, Zhang12}. As in \cite{GangoMaulik, Renlund, Renlund2}, we use random step sizes instead of deterministic ones considered in \cite{Zhang}, leading to a Lotka-Volterra type first order quadratic differential equation.

In Section~\ref{sec: lv}, we state the model and the main results. In Corollary~\ref{cor: iid}, we extend the result for i.i.d.\ replacement matrices in \cite{AthreyaNey, Zhang12} and prove it under $L^1$ assumptions rather than $L \log_+ L$ conditions. Corollary~\ref{cor: L log L} studies adapted replacement matrices, when they are majorized by a $L \log_+ L$ random variable. In Section~\ref{sec:urn}, the main tool for the proof -- an appropriate result on stochastic approximation -- was provided and the main results were proven. Finally, in Section~\ref{sec:erw}, we apply the results for the urn model to study a delayed elephant random walk on $d$-dimensional nonnegative orthant with randomly reinforced memory.

\section{Urn Model and Main Result} \label{sec: lv}
In this Section, we first collect some notations. In two subsequent subsections, we describe the model; and state the main result. In the following, all the vectors are row vectors, while the transpose of a row vector $\x$ will be $\x^T$. The vector $\bds{1}$ has all components $1$, the dimension being clear from the context. Also, for an event $A$, we shall denote its indicator function by $\bone_A$. For a matrix $\boldsymbol{A}$, define the norm $\mynorm{\boldsymbol{A}} := \max_i\sum_{j}\left|A_{ij}\right|.$ We shall also use the notation $\sigma\left(\boldsymbol{A}\right) := \min_i\sum_{j}\left|A_{ij}\right|$.

\subsection{Urn Model}
Consider an urn model with $K$ colors indexed by $\{1, 2, \ldots, K\}$. The composition vector after $n$-th trial is $\C_n = (C_{n1}, \ldots, C_{nK})$. The total content of the urn after $n$-th trial is $S_n = \sum_{i=1}^K C_{ni}$. The replacement matrix for $n$-th trial will be a random $K\times K$ matrix $\R_n$, with nonnegative entries. For $n$-th trial, $\Chi_n$ is the indicator of the color chosen in that trial, taking value $\bds{e}_i$, the $i$-th coordinate vector in $\ensuremath{\mathbb{R}}^K$ if the color $i$ is chosen in $n$-th trial. In this case, $i$-th row of $\R_n$ is added to the composition vector. Therefore, for $n\geq 1$, the urn composition evolves as:
\begin{equation}
\C_n = \C_{n-1} + \Chi_n\R_n.
\label{EvolutionEqn}
\end{equation}
The count vector after $n$-th trial, $\N_n$, giving the number of times each color is drawn till $n$-th trial, is $\N_n = \sum_{k=1}^n \Chi_k$. For $n\geq 1$, $\FF_{n}$ is the $\sigma$-field containing the information till time $n$, namely, $\FF_n$ is the $\sigma$-field generated by $\C_0$, $\left(\R_m, \Chi_m\right)_{m=1}^{n}$. The conditional expectations of the replacement matrices, denoted as $\bH_{n-1} := \CExp{\R_n}{n-1}$ and their truncated version, denoted as $\widetilde{\bH}_{n-1} := \CExp{\R_n \indic{\mynorm{\R_n}\le n}}{n-1}$, called the \textit{generating matrices} and the \textit{truncated generating matrices} respectively.

We make the following assumptions on choice of colors and reinforcement matrices.
\begin{assumption} \label{Assmp: basic}
The adapted sequence $((\Chi_n, \R_n), \FF_n)$ has distribution, which satisfies:
\begin{enumerate}[label={\normalfont (\roman*)},ref={\normalfont \roman*}]
\item The initial configuration $\C_0$ is nonzero with nonnegative entries and finite mean. \label{Assmp: init}
\item For all $i$, $\Prob\left(\Chi_n=\boldsymbol{e}_i|\C_0,(\R_m)_{m=1}^{n-1}, (\Chi_m)_{m=1}^{n-1}\right)=C_{n-1,i}/S_{n-1}$.
\item Given the past information $\FF_{n-1}$, $\Chi_n$ and $\R_n$ are conditionally independent.
\end{enumerate}
\end{assumption}

We also make uniform integrability type conditions on the replacement matrices.
\begin{assumption}
The sequence $(\mynorm{\R_n})$ satisfy one of the following: \label{Assmp: unif tail}
\begin{enumerate}[label={\normalfont (\alph*)},ref={\normalfont \alph*}]
  \item \label{Assmp: maj} The unconditional distributions of $\mynorm{\R_n}$ are majorized: there exists $c\in (0,\infty)$ and a positive random variable $R$ with finite expectation such that, for all $x>0$, \[\Prob\left(\mynorm{\R_n}>x\right) \le c \Prob(R>x).\]
  \item \label{Assmp: alt tail} For some nonnegative function $\phi$ on $[0, \infty)$ which is eventually positive and nondecreasing, satisfying $\sum 1/(n\phi(n))< \infty$ and $x/\phi(x)$ eventually monotone nondecreasing, we have, the following bounded moment condition:
      \begin{equation}
      \sup_n \Exp(\mynorm{\R_n}\phi(\mynorm{\R_n}) < \infty. \label{eq: exp bd}
      \end{equation}
\end{enumerate}
\end{assumption}

\begin{remark}
Under Assumption~\ref{Assmp: unif tail}, for all $n, i, j$, $R_{nij}$ have finite mean. \label{Rem: fin mean}
\end{remark}

\begin{remark} \label{rem: identical}
Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj} holds if $\left(\mynorm{\R_n}\right)$ have same marginal distribution.
\end{remark}

\begin{remark} \label{Rem: alt tail} One of the useful choices for $\phi(x) = (\log_+ x)^p$ for some $p>1$, has been considered by~\cite{Zhang12}. Other choices can include $\phi(x) = \log_+x (\log_+ \log_+ x)^p$ for some $p>1$.
\end{remark}

Finally, we state the conditions for closeness of the replacement matrices to an irreducible matrix.
\begin{assumption}
There exists a (possibly random) almost surely irreducible matrix $\bH$ with nonnegative entries, such that the truncated generating matrices $\widetilde{\bH}_n$ converge to $\bH$ in Cesaro sense in the operator norm almost surely, i.e., \label{Assmp: ae}
\[\frac1n \sum_{k=0}^{n-1} \mynorm{\widetilde{\bH}_{k}-\bH}\to 0 \text{ almost surely.}\]
The dominant eigenvalue of $\bH$ will be denoted by $\pfval$. By Perron-Frobenius theory, there exists unique left eigenvector $\lpf$ of $\bH$ corresponding to $\pfval$ such that $\lpf$ has all coordinates strictly positive and is normalized to be a probability vector.
\end{assumption}

\begin{remark} \label{Rem: sigma}
Since $\bH$ is irreducible, no row can be the zero vector and hence $\sigma(\bH)>0$.
\end{remark}

\begin{remark} \label{Rem: L1}
Both Assumptions~\ref{Assmp: unif tail}\eqref{Assmp: maj} and~\ref{Assmp: unif tail}\eqref{Assmp: alt tail} imply uniform integrability of  $\left(\mynorm{\R_n}\right)$ and, hence, of $\left(\R_n\indic{\mynorm{\R_n}\le n}\right)$. Thus, under either condition, the convergence in Assumption~\ref{Assmp: ae} is also in $L^1$.
\end{remark}

\begin{remark}
Different sufficient notions of closeness between the generating matrix and the limiting matrix $\bH$ have been considered in \cite{BaiHu, GangoMaulik, LaruellePages, Zhang12, Zhang}. See \cite{GangoMaulik} for further discussion on such assumptions.
\end{remark}

\subsection{Main Results} \label{subsec: result}
In this subsection, we state the main result and some corollaries covering the most relevant situations.
\begin{theorem} \label{Main}
Under Assumptions~{\normalfont \ref{Assmp: basic}},~{\normalfont \ref{Assmp: unif tail}} and~{\normalfont \ref{Assmp: ae}}, we have, almost surely, as well as in $L^1$,
\begin{equation}\label{strong conv}
  \frac{\C_n}{S_n} \to \boldsymbol{\pi}_{\bH}; \quad \frac{S_n}n \to \lambda_{\bH}; \quad
  \frac{\C_n}n \to \lambda_{\bH} \boldsymbol{\pi}_{\bH}; \quad \text{and} \quad
  \frac{\N_n}n \to \boldsymbol{\pi}_{\bH}.
\end{equation}
\end{theorem}

\begin{remark}
It is interesting to compare the assumptions in this article with those required to establish convergence in probability and in $L^1$ in \cite{GangoMaulik}. Note that Assumption~\ref{Assmp: basic} is same as Assumption~3.1 of \cite{GangoMaulik}. As noted in Remark~\ref{Rem: fin mean}, the finite mean of replacement matrices, as in Assumption~3.1~(ii) of \cite{GangoMaulik} follows from Assumption~\ref{Assmp: unif tail}.

As noted in Remark~\ref{Rem: L1}, Assumption~\ref{Assmp: unif tail} is the analog of the uniform integrability of $(\mynorm{\R_n})$ in Assumption~3.4 of \cite{GangoMaulik}. However, in this article, we need somewhat stronger conditions than uniform integrability, namely, either majorization condition in Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj} or appropriate uniform moment bounds in Assumption~\ref{Assmp: unif tail}\eqref{Assmp: alt tail}, as has been considered by~\cite{Zhang12}.

Finally, Assumption~\ref{Assmp: ae} corresponds to Assumptions~3.2 and~3.3 of \cite{GangoMaulik}. While Assumption~3.2 of \cite{GangoMaulik} on the properties of the matrix $\bH$ remain unchanged, given the almost sure convergence in conclusion, we naturally strengthen the mode of convergence to almost sure from convergence in probability. We further require the truncated generating matrices $\widetilde{\bH}_n$ to converge to $\bH$ instead of the generating matrices, as considered in \cite{GangoMaulik}. We shall provide sufficient conditions in Proposition~\ref{prop: suff tail ce} for convergence of $\bH_n$, as in \cite{GangoMaulik}, to become equivalent to Assumption~\ref{Assmp: ae}.
\end{remark}

There are simpler sufficient conditions which imply Assumptions~\ref{Assmp: basic},~\ref{Assmp: unif tail} and~\ref{Assmp: ae} and improve existing results. Chapter~V.9.3 of \cite{AthreyaNey}, and Theorem~2.1(a) of \cite{Zhang12} considered the urn model with i.i.d.\ replacement matrices having $L \log_+ L$ moments. In the following corollary we show $L^1$ moments suffice.

\begin{corollary} \label{cor: iid}
Let the adapted sequence $\left( \left( \Chi_n, \R_n\right), \FF_n\right)$ satisfy:
\begin{enumerate}[label={\normalfont (\roman*)},ref={\normalfont \roman*}]
\item The initial configuration $\C_0$ is nonzero with nonnegative entries and finite mean.
\item The replacement matrices $\R_n$ are i.i.d.\ with $\Exp\left(\mynorm{\R_1}\right)<\infty$. \item The colors chosen form the sequence $\left(\Chi_n\right)$, obtained using an i.i.d.\ Uniform$(0,1)$ sequence $(U_n)$ independent of the sequence $(\R_n)$ as:
\[\Chi_n = \sum_{i=1}^d \boldsymbol{e}_i \indic{U_n \in (\sum_{j=1}^{i-1} C_{n-1,j}/S_{n-1}, \sum_{j=1}^i C_{n-1,j}/S_{n-1}]}.\]
\end{enumerate}
Then the convergence~\eqref{strong conv} in Theorem~{\normalfont \ref{Main}}, namely,
\[\frac{\C_n}{S_n} \to \boldsymbol{\pi}_{\bH}; \quad \frac{S_n}n \to \lambda_{\bH}; \quad \frac{\C_n}n \to \lambda_{\bH} \boldsymbol{\pi}_{\bH}; \quad \text{and} \quad \frac{\N_n}n \to \boldsymbol{\pi}_{\bH}\]
hold both almost surely and in $L^1$.
\end{corollary}

While only $L^1$ condition suffices for i.i.d.\ replacement matrices, $L \log_+ L$ condition allows us to study adaptive urns. It is also natural to consider the generating matrices in Assumption~\ref{Assmp: ae} instead of the truncated ones.
\begin{assmprime}{Assmp: ae} \label{Assmp: aeprime}
There exists a (possibly random) matrix $\bH$ with nonnegative entries, which is almost surely irreducible, such that the generating matrices $\bH_n$ converge to $\bH$ in Cesaro sense in the operator norm almost surely:
    \[\frac1n \sum_{k=0}^{n-1} \mynorm{\bH_{k}-\bH}\to 0 \text{ almost surely.}\]
\end{assmprime}

\begin{remark} \label{Rem: tail ce}
Assumption~\ref{Assmp: aeprime} becomes equivalent to Assumption~\ref{Assmp: ae}, when considered along with the following negligibility condition on the tail conditional expectation:
\begin{equation}
\frac1n \sum_{k=1}^{n} \CExp{\mynorm{\R_{k}}\indic{\mynorm{\R_{k}}>k}}{k-1}\to 0 \text{ almost surely.} \label{eq: tail ce neg}
\end{equation}
\end{remark}

The following corollary shows Assumption~\ref{Assmp: aeprime}, along with $L \log_+ L$ condition on the majorizing variable in Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj} is again sufficient.
\begin{corollary} \label{cor: L log L}
The adapted sequence $\left( \left( \Chi_n, \R_n\right), \FF_n\right)$ satisfy Assumptions~{\normalfont\ref{Assmp: basic}},~{\normalfont\ref{Assmp: unif tail}\eqref{Assmp: maj}} and~{\normalfont\ref{Assmp: aeprime}}, where the majorizing random variable $R$, considered in Assumption~{\normalfont\ref{Assmp: unif tail}\eqref{Assmp: maj}} satisfy $\Exp(R \log_+ R) < \infty$. Then the convergence~\eqref{strong conv} in Theorem~{\normalfont \ref{Main}}, namely, ${\C_n}/{S_n} \to \boldsymbol{\pi}_{\bH}$; ${S_n}/n \to \lambda_{\bH}$; ${\C_n}/n \to \lambda_{\bH} \boldsymbol{\pi}_{\bH}$ and ${\N_n}/n \to \boldsymbol{\pi}_{\bH}$ hold both almost surely and in $L^1$.
\end{corollary}
Corollary~\ref{cor: L log L} extends the result on i.i.d.\ sequence of replacement matrices given in \cite{AthreyaNey} to adapted sequence with similar $L \log_+ L$ moment for the majorizing random variable. This corollary is also similar to the result in~\cite{Zhang12} for adapted sequence with $L (\log_+ L)^p$-bounded moments, for some $p>1$. In this corollary $L (\log_+ L)^p$-boundedness conditions is replaced by $L \log_+ L$ moment of the majorizing random variable.

The following corollary generalizes $L (\log_+ L)^p$ condition, for some $p>1$, used in Assumption~2.2(a) of \cite{Zhang12} in terms of Assumption~\ref{Assmp: unif tail}\eqref{Assmp: alt tail}. See also Remark~\ref{Rem: alt tail}.
\begin{corollary} \label{cor: phi}
Assumptions~{\normalfont \ref{Assmp: basic}},~{\normalfont \ref{Assmp: unif tail}\eqref{Assmp: alt tail}} and~\ref{Assmp: aeprime} imply the convergence~\eqref{strong conv} in Theorem~{\normalfont \ref{Main}} both almost surely and in $L^1$.
\end{corollary}

\section{Proofs of the Main Results} \label{sec:urn}
In the first subsection, we prove the relevant result on stochastic approximation and a couple other technical results. In the next subsection, we check the conditions of stochastic approximation for the urn model. In the last two subsections, we prove the theorem and corollaries of Section~\ref{subsec: result}.

\subsection{Stochastic Approximation and Other Results}
The main tool for proving Theorem~\ref{Main} is Stochastic Approximation. We shall use the following theorem on Stochastic Approximation, which follows from a special case of Theorem~9.2.8 of \cite{Duflo}, which in turn was adapted from \cite{KushnerClark}.
\begin{theorem} \label{Thm: SA}
Let $\left(\x_n\right)_{n=0}^\infty$ be a sequence that takes values in a compact convex subset $S$ of $\ensuremath{\mathbb{R}}^K$. Let $\left( a_n\right)_{n=0}^\infty$ be a nonincreasing sequence of positive real numbers, called the step sizes, such that
\begin{equation} \label{eq: step order}
0<\liminf {n}{a_n} \le \limsup {n}{a_n} <\infty.
\end{equation}
Define for $n\geq 0$ and $T>0$
\begin{equation} \label{eq: tau}
\Tau{n}{T}:=\inf\left\{  k\geq n: a_n+\cdots+a_k\geq T\right\}.
\end{equation}
Let $\left(\bgamma_n\right)_{n=0}^\infty$ be a sequence, called the error sequence, in $\ensuremath{\mathbb{R}}^K$ such that $\left(\bgamma_n\right)$ goes to $0$ in Cesaro sense, i.e., $\frac1n \bds{\Gamma}_n \to \bds{0}$, where $\bds{\Gamma}_n := \sum_{i=1}^n \bgamma_i$.
Let $h:S\to\ensuremath{\mathbb{R}}^K$ be a continuous function. Suppose for all $n\geq 0$
\[
\x_{n+1}=\x_n+a_n h(\x_n)+a_n\bgamma_n.
\]
If the ODE $\dot{\x}=h(\x)$ has a unique solution $\z:\ensuremath{\mathbb{R}}\to S$, then $\lim_{n\to\infty}\x_n=\z(0)$.
\end{theorem}
\begin{proof}
It is immediate from~\eqref{eq: step order} that $\lim_{n\to\infty} a_n=0$ and $\sum_{n=0}^\infty a_n=\infty$. Then, using Theorem~9.2.8 of \cite{Duflo}, it is enough to show that, for all $T>0$,
\[\lim_{n\to\infty}\sup_{n\leq j\leq\Tau{n}{T}}\Lone{\sum_{k=n}^j a_k\bgamma_k}=0.\]

Fix $T>0$. First observe that
\[a_n\bgamma_n = \left( a_n\bds{\Gamma}_n - a_{n-1}\bds{\Gamma}_{n-1}\right) + a_n a_{n-1}\bds{\Gamma}_{n-1} \left( \frac1{a_n} - \frac1{a_{n-1}} \right).\]
Since $(a_n)$ is nonincreasing, for $n\leq j\leq\Tau{n}{T}$, we have,
\[\Lone{\sum_{k=n}^j a_k\bgamma_k} \le \Lone{a_j\bds{\Gamma}_j} + \Lone{a_{n-1}\bds{\Gamma}_{n-1}} + \sup_{k\ge n-1} \Lone{a_k \bds{\Gamma}_k} \frac{a_n}{a_j}.\]
Hence, we have
\[\sup_{n\leq j\leq\Tau{n}{T}}\Lone{\sum_{k=n}^j a_k\bgamma_k} \le \sup_{j\ge n-1} \Lone{a_j\bds{\Gamma}_j} \left( 2 + na_n \cdot \frac{\Tau{n}{T}}{n} \cdot \frac1{\Tau{n}{T}a_{\Tau{n}{T}}}\right).\]
Again, using~\eqref{eq: step order}, we have $n a_n$ bounded away from both $0$ and $\infty$, and, further, $a_n \bds{\Gamma}_n \to \bds{0}$. So, it is enough to show that $\Tau{n}{T}/n$ is bounded for any $T>0$. Since $\liminf na_n >0$, choose $a^*>0$ such that $n a_n > a^*$ for all large enough $n$. Further note that $a_n < T$ for all large enough $n$. Hence, from~\eqref{eq: tau}, we have $\Tau{n}{T}>n$ for all large enough $n$ and
\[a^*\left[\log (\Tau{n}{T}-1) - \log (n-1) -1\right] \le \frac{a^*}{n} + \cdots + \frac{a^*}{\Tau{n}{T}-1} \le a_n + \ldots + a_{\Tau{n}{T}-1} \le T.\]
Thus, $\Tau{n}{T}/n \le (\Tau{n}{T}-1)/(n-1) \le \exp(T/a^*+1)$, as required.
\end{proof}

Next we prove two technical results. The first result studies average of a dominated adapted sequence, where the terms are centered by a truncated conditional expectation.
\begin{proposition} \label{prop: trunc mg}
Let $(X_n, \FF_n)$ be a adapted sequence, such that $|X_n|\le c' \mynorm{\R_n}$ for some $c'\in(0,\infty)$. Then, under Assumptions~{\normalfont \ref{Assmp: basic}} and~{\normalfont \ref{Assmp: unif tail}}, we have, with probability 1,
$\frac1n \sum_{k=1}^n \left( X_k - \CExp{X_k \indic {\mynorm{\R_k}\le k}}{k-1} \right) \to 0$. \end{proposition}
\begin{proof}
The proof follows Theorem~2.19 of \cite{HH}. By Kronecker's lemma, it is enough to show
\[
\sum_{k=1}^n\frac1{k} \left( X_k - \CExp{X_k \indic {\mynorm{\R_k}\le k}}{k-1}\right)< \mbox{ converges almost surely.}
\]
We write
\begin{multline} \label{eq: trunc mg}
  \sum_{k=1}^n \frac1k \left( X_k - \CExp{X_k \indic{\mynorm{\R_k}\le k}}{k-1} \right)\\
  = \sum_{k=1}^n \frac1k \left[\MGD{X_k \indic{\mynorm{\R_k}\le k}}{k-1}\right] + \sum_{k=1}^n \frac1k X_k \indic{\mynorm{\R_k}> k}.
\end{multline}

The first sum is a martingale and, since $|X_n|\le c' \mynorm{\R_n}$, it is almost surely convergent if
\[
\sum_n \frac{1}{n^2}\CExp{\mynorm{\R_n}^2 \indic{\mynorm{\R_n}\le n}}{n-1}< \infty\mbox{ almost surely.}
\]
Under Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj}, we show its expectation to be finite. We have
\begin{align*}
  &\sum_{n\ge1} \frac1{n^2} \Exp\left(\mynorm{\R_n}^2 \indic{\mynorm{\R_n}\le n}\right) = \sum_n \frac1{n^2} \int_0^n 2r \Prob(\mynorm{\R_n}>r) dr\\
  &\qquad \qquad \le 2c \sum_{j\ge1} \sum_{n\ge j} \frac1{n^2} \int_{j-1}^j r \Prob(R>r) dr \le 8c \int_0^\infty \Prob(R>r) dr < \infty.
\end{align*}
Under Assumption~\ref{Assmp: unif tail}\eqref{Assmp: alt tail}, for all large enough $n$, we have
\[\CExp{\mynorm{\R_n}^2 \indic{\mynorm{\R_n}\le n}}{n-1} \le \frac{n}{\phi(n)} \CExp{\mynorm{\R_n} \phi(\mynorm{\R_n})}{n-1},\]
since $x/\phi(x)$ is eventually monotone nondecreasing. The result follows from summability of $1/(n\phi(n))$.

For the second sum on the right side of~\eqref{eq: trunc mg}, note that $\Prob(\mynorm{\R_n}>n)$ is summable under Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj} by finiteness of the moment of the majorizing distribution. Markov type inequality gives
\[\Prob(\mynorm{\R_n}>n) \le \frac1{n \phi(n)} \Exp{\mynorm{\R_n} \phi(\mynorm{\R_n})},\] which is again summable under Assumption~\ref{Assmp: unif tail}\eqref{Assmp: alt tail}. Then, by Borel-Cantelli lemma, the second sum on the right of~\eqref{eq: trunc mg} has at most finitely many terms with nonzero indicator functions and hence converges almost surely.
\end{proof}

The next result provides SLLN for martingale differences.
\begin{proposition} \label{prop: mg SLLN}
Let $\left( U_n, \FF_n \right)$ be a martingale difference sequence. Let $\left( V_n \right)$ be a sequence of random variables having common distribution with finite mean, such that $V_n$ is independent of $\FF_{n-1}$ and $|U_n|\le V_n$ for all $n$. Then $\frac1n \sum_{k=1}^{n} U_k \to 0$ almost surely.
\end{proposition}
\begin{proof}
Since the martingale difference sequence is dominated by a sequence of random variables with common distribution with finite mean, as in derivation of~(2.20) of~\cite{HH}, we have
$\frac1n \sum_{k=1}^{n} \left( U_k - \CExp{U_k \indic{[|U_k|\le k]}}{k-1} \right) \to 0$ almost surely.

Further, since $V_n$ dominates $U_n$, we have
$\CExp{U_n \indic{[|U_n| > n]}}{n-1} \le \Exp{V_n \indic{[V_n| > n]}} \to 0$, which takes care of the remainder term in the centering.
\end{proof}

\subsection{Analysis of the Urn Model} \label{subsec: analysis}
In this subsection, we check that the evolution equation~\eqref{EvolutionEqn} of the urn model forms stochastic approximation and check the conditions of Theorem~\ref{Thm: SA}.

Total amount added to the urn in each step is denoted by $Y_n$ for $n\ge 0$. Define $Y_0:=S_0$ and for $n\geq 1$, $Y_n:=S_n-S_{n-1}=\Chi_n\R_n\boldsymbol{1}^T$. It is immediate to observe that
\begin{equation}
  0\leq Y_n=\Chi_n\R_n\boldsymbol{1}^T\leq\mynorm{\R_n} \text{ almost surely.} \label{Y bd}
\end{equation}

Following~\cite{GangoMaulik}, we rewrite the evolution equation \eqref{EvolutionEqn} as
\begin{equation} \label{eq: Mod SA}
\csd{n} = \csd{n-1} + \frac{1}{S_n}h_{\bH}\left(\csd{n-1}\right) + \frac{\bds{\delta}_n}{S_n} + \frac{\bxi_n}{S_n},
\end{equation}
where $h_{\bH}$, $\bds{\delta}_n$, $\bxi_n$ are defined as follows. The drift $h_{\bH}$, indexed by $K\times K$ matrices, is defined as $h_{\bH}\left(\bX\right) := \bX\bH - \bX\left(\bX\bH\boldsymbol{1}^T\right)$. For each $n\geq 1$, define the martingale difference term
\begin{align*}
\bds{\delta}_n := &\left(\Chi_n\R_n-\csd{n-1}Y_n\right) - \CExp{\left(\Chi_n\R_n-\csd{n-1}Y_n\right) \indic{\mynorm{\R_n}\le n}}{n-1},
\intertext{and adjusted truncated conditional expectation term}
\bxi_n := &\csd{n-1}\left(\widetilde{\bH}_{n-1}-\bH\right) - \csd{n-1}\left(\csd{n-1}\left(\widetilde{\bH}_{n-1}-\bH\right)\boldsymbol{1}^T\right)
= h_{\widetilde{\bH}_{n-1}-\bH} \left(\csd{n-1}\right),
\end{align*}
We shall refer to~\eqref{eq: Mod SA} as the stochastic approximation equation. Note that the stochastic approximation for the sequence $\left(\csl{n}\right)$ takes values in the closed bounded convex set of probability simplex in $\ensuremath{\mathbb{R}}^K$. The corresponding ODE $\dot{\x}=h_{\bH}\left(\x\right)$ defined in terms of $h_{\bH}$ is a first order quadratic equation of Lotka-Volterra type. Since $h_{\bH}(\x)$ is polynomial in $\x$, it is continuous, as required in Theorem~\ref{Thm: SA}. Further, we have, from Proposition~3.3 of~\cite{GangoMaulik}, that the Lotka-Volterra equation has unique solution in the probability simplex given by $\x(t)=\lpf$ for all $t$.

Next we consider the step sizes. Unlike the usual step size $(1/n)$, we use, as in \cite{GangoMaulik} the random step size $1/S_{n+1}$ of the stochastic approximation equation~\eqref{eq: Mod SA}. The step size $1/S_{n+1}$ has also been used to analyze two color urn models with fixed replacement unbalanced matrices in \cite{Renlund, Renlund2}. The choice of this random step size makes $\C_n/S_n$ take values in the $k$-dimensional probability simplex, which is compact and convex, as required in Theorem~\ref{Thm: SA}. Since $S_n$'s are partial sum of nonnegative $Y_n$'s, the step sizes are nonincreasing, as required in Theorem~\ref{Thm: SA}. In the next proposition, we check the condition~\eqref{eq: step order} as required in Theorem~\ref{Thm: SA}.

\begin{proposition} Under Assumptions~{\normalfont \ref{Assmp: basic}},~{\normalfont \ref{Assmp: unif tail}} and~{\normalfont \ref{Assmp: ae}}, we have:
\[0 < \sigma\left(\bH\right) \leq \liminf_{n\to\infty} \frac{S_n}n \leq \limsup_{n\to\infty} \frac{S_n}n \leq \mynorm{\bH} < \infty.\]
\end{proposition}
\begin{proof} Define
\[
\eta_n := \frac{S_n}n - \frac1n \sum_{m=1}^n \CExp{Y_m \indic {\mynorm{\R_m}\le m}}{m-1} + \frac1n \sum_{m=1}^n \csd{m-1}\left(\widetilde{\bH}_{m-1}-\bH\right)\boldsymbol{1}^T.
\]
Since $0\le Y_n\le \mynorm{\R_n}$, using Proposition~\ref{prop: trunc mg}, sum of the first two terms of $\eta_n$ satisfies
\begin{equation}\label{eq: eta two}
  \frac{S_n}n - \frac1n \sum_{m=1}^n \CExp{Y_m \indic {\mynorm{\R_m}\le m}}{m-1} \to 0 \text{ almost surely.}
\end{equation}
The third term of $\eta_n$ is bounded by $\frac1n \sum_{m=1}^n \mynorm{\widetilde{\bH}_{m-1}-\bH}$, which is also negligible almost surely by Assumption~\ref{Assmp: ae}. Thus $\eta_n\to 0$ almost surely. Further, we have
\[
0 < \sigma\left(\bH\right) \leq \left|\frac{S_n}n - \eta_n\right| \leq \mynorm{\bH} < \infty,
\]
which, along with Remark~\ref{Rem: sigma}, proves the result.
\end{proof}

Finally, we show Cesaro negligibility of the error terms $(\bds{\delta}_n)$ and $(\bxi_n)$ of the stochastic approximation equation~\eqref{eq: Mod SA}, as required in Theorem~\ref{Thm: SA}. \begin{lemma} Under Assumptions~{\normalfont \ref{Assmp: basic}} and~{\normalfont \ref{Assmp: unif tail}}, we have, with probability $1$, $\frac1n\sum_{m=1}^n \bds{\delta}_m \to \bds{0}$ and $\frac1n\sum_{m=1}^n \bxi_m \to \bds{0}$.
\end{lemma}
\begin{proof}
Each coordinate of $(\Chi_n\R_n-\csl{n-1}Y_n)$ is bounded by $2\mynorm{\R_n}$. Hence, using Proposition~\ref{prop: trunc mg}, each coordinate of $\frac1n\sum_{m=1}^n \bds{\delta}_m$ is negligible. Further, recall that \[\bxi_n = \csd{n-1}\left(\widetilde{\bH}_{n-1}-\bH\right) - \csd{n-1}\left(\csd{n-1}\left(\widetilde{\bH}_{n-1}-\bH\right)\boldsymbol{1}^T\right).\] Then $\Lone{\bxi_n} \le 2 \mynorm{\widetilde{\bH}_{n-1}-\bH}$, and by Assumption~\ref{Assmp: ae}, the result follows.
\end{proof}

\subsection{Proof of the Main Theorem}
We prove Theorem~\ref{Main} in this Subsection.

\begin{proof}[Proof of Theorem~{\normalfont\ref{Main}}]
Since the relevant conditions of Theorem~\ref{Thm: SA} have been checked for the proportion vector $\left(\C_n/S_n\right)$ in Section~\ref{subsec: analysis}, we have the almost sure convergence of the proportion vector to $\boldsymbol{\pi}_{\bH}$.

Next we consider the almost sure convergence of $S_n/n$. Using~\eqref{eq: eta two}, it is enough to show that
\begin{equation} \label{eq: Cesaro}
\frac{1}n\sum_{m=1}^n\CExp{Y_m\indic {\mynorm{\R_m}\le m}}{m-1}\to\pfval
\end{equation}
almost surely. Recall that
\[\frac{1}{n}\sum_{m=1}^n\CExp{Y_m\indic {\mynorm{\R_m}\le m}}{m-1}\\ = \frac{1}{n}\sum_{m=1}^n\csd{m-1}\left(\widetilde{\bH}_{m-1}-\bH\right)\boldsymbol{1}^T + \frac{1}{n}\sum_{m=1}^n\csd{m-1}\bH\boldsymbol{1}^T.\]
The first term on the right side converges to 0 almost surely by Assumption~\ref{Assmp: ae}. Further, since $\csl{n}\to\lpf$ almost surely, we get
$\left(\C_n/S_n\right)\bH\boldsymbol{1}^T\to\lpf\bH\boldsymbol{1}^T=\pfval\lpf\boldsymbol{1}^T=\pfval$
almost surely and~\eqref{eq: Cesaro} follows using Cesaro. Then, the almost sure limit of $\C_n/n$ follows easily from those of $\C_n/S_n$ and $S_n/n$.

From the almost sure convergence of the proportion vector, $\frac{1}{n}\sum_{m=0}^{n-1}\C_m/S_m \to \lpf$ holds almost surely. Hence, the almost sure limit of $\N_n/n$ will follow from
\[
\frac{\N_n}{n} - \frac{1}{n}\sum_{m=0}^{n-1}\csd{m} = \frac{1}{n}\sum_{m=1}^n\left(\Chi_m-\csd{m-1}\right) \to \boldsymbol{0} \text{ almost surely.}
\]
By Kronecker's Lemma, it is enough to show that $\sum_n\frac{1}{n}\left(\Chi_n-\C_{n-1}/S_{n-1}\right)$ is convergent almost surely, which holds since $\Lone{\Chi_n/n}^2=1/n^2$ is
summable.

Next we consider $L^1$ convergence. Since $\csl{n}$ and $\N_n/n$ are bounded, their almost sure convergence will imply convergence in $L^1$ as well by Dominated Convergence Theorem. If the convergence of $S_n/n$ is shown to hold in $L^1$, then by Lemma~3.5 of \cite{GangoMaulik}, the convergence of $\C_n/n$ will be in $L^1$ as well. So we consider $L^1$ convergence of $S_n/n$.

Observe that, from~\eqref{Y bd}, $0\leq Y_n\leq\mynorm{\R_n}$ for $n\ge 1$ and $Y_0= \C_0 \bds{1}^T$ is integrable by Assumption~\ref{Assmp: basic}\eqref{Assmp: init}. Further, as noted in Remark~\ref{Rem: L1}, Assumption~\ref{Assmp: basic}\eqref{Assmp: maj} or bounded moment condition~\eqref{eq: exp bd} of Assumption~\ref{Assmp: unif tail}\eqref{Assmp: alt tail} implies uniform integrability of $(\mynorm{\R_n})$. Hence $\left( Y_n\right)_{n=0}^\infty$, and thus $S_n/n=\sum_{k=0}^n Y_k/n$ are also uniformly integrable. So the almost sure convergence of $S_n/n$ becomes $L^1$ as well.
\end{proof}

\subsection{Proofs of Corollaries~\ref{cor: iid},~\ref{cor: L log L} and~\ref{cor: phi}}
Since the conclusions of Corollaries~\ref{cor: iid},~\ref{cor: L log L} and~\ref{cor: phi} are same as those of Theorem~\ref{Main}, we check that the assumptions of the corollaries imply those of the theorem.

\begin{proof}[Proof of Corollary~{\normalfont \ref{cor: iid}}]
  Assumption~\ref{Assmp: basic} clearly holds. Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj} holds as in Remark~\ref{rem: identical}. Further, the i.i.d.\ sequences $(\Chi_n)$ and $(\R_n)$ give Assumption~\ref{Assmp: aeprime}. Also, finite mean of $\mynorm{\R_n}$ implies~\eqref{eq: tail ce neg}. Assumption~\ref{Assmp: ae} then follows from Remark~\ref{Rem: tail ce}.
\end{proof}

Next we consider the other two corollaries.
\begin{proof}[Proof of Corollaries~~{\normalfont\ref{cor: L log L}} and~{\normalfont\ref{cor: phi}}]
Note that Corollaries~\ref{cor: L log L} and~\ref{cor: phi} assume Assumption~\ref{Assmp: basic},~\ref{Assmp: unif tail} and~\ref{Assmp: aeprime}, together with $L \log_+ L$ condition on the majorizing random variable. Hence, using Remark~\ref{Rem: tail ce}, it is enough to check~\eqref{eq: tail ce neg}, which we do next.
\end{proof}

\begin{proposition} \label{prop: suff tail ce}
Assumption~{\normalfont \ref{Assmp: basic}} and either of the following is sufficient for~\eqref{eq: tail ce neg}:
\begin{enumerate}[label={\normalfont (\roman*)},ref={\normalfont \roman*}]
\item Assumption~{\normalfont \ref{Assmp: unif tail}\eqref{Assmp: maj}} and the majorizing variable $R$ satisfies $\Exp(R \log_+R)<\infty$.
\item Assumption~{\normalfont \ref{Assmp: unif tail}\eqref{Assmp: alt tail}}.
\label{Assmp: phi}
\end{enumerate}
\end{proposition}

\begin{proof}
By Kronecker's lemma, it is enough to show
$\sum_{n=1}^{\infty} \frac1n \CExp{\mynorm{\R_{n}}\indic{\mynorm{\R_{n}}>n}}{n-1} < \infty \text{ almost surely}$.
\begin{enumerate}[label={\normalfont (\roman*)},ref={\normalfont \roman*}]
\item We show that the random variable on the left side above has finite expectation. Now, using Assumption~\ref{Assmp: unif tail}\eqref{Assmp: maj}, we have
\begin{align*}
&\sum_{n=1}^{\infty} \frac1n \Exp\left(\mynorm{\R_{n}}\indic{\mynorm{\R_{n}}>n}\right) \leq c\sum_{n=1}^{\infty} \frac1n \Exp\left( R\indic{R>n}\right)\\
&\qquad \qquad \le c \sum_{n=1}^{\infty} \sum_{j=n}^\infty \frac1n \Exp\left( R\indic{R\in (j, j+1]}\right) \le c \sum_{j=1}^{\infty} \Exp\left( R\indic{R\in (j, j+1]}\right) \sum_{n=1}^j \frac1n\\
&\qquad \qquad \le c \sum_{j=1}^{\infty} \Exp\left( R\indic{R\in (j, j+1]} (1+ \log_+ R)\right) \le c \Exp\left( R (1+ \log_+ R)\right) < \infty.
\end{align*}
\item Since $\phi$ is eventually nondecreasing, we have, for some $m$,
\[
\sum_{n=m}^\infty \frac1n \Exp{\mynorm{\R_{n}}\indic{\mynorm{\R_{n}}>n}} \le \sum_{n=m}^\infty \frac1{n \phi(n)} \Exp{\mynorm{\R_{n}}\phi\left(\mynorm{\R_{n}}\right)}
\]
and the result follows from the summability of $1/(n\phi(n))$.
\end{enumerate}
\end{proof}

\section{Applications to Elephant Random Walk} \label{sec:erw}
In this section, we consider the delayed elephant random walk (ERW) on nonnegative integer lattice in dimension $d$ with randomly reinforced memory. The walk was introduced in \cite{ERWoriginal} and has been studied extensively in both statistical physics and literature. We analyze the model using the results obtained in this article utilizing an interesting connection between the urn model and ERW, which was discovered in \cite{Bertoin}.

The location of the walk at the epoch $n$ is denoted by $\bds{L}_n$, while the step taken at that epoch is denoted by $\bds{X}_n$, where $\bds{L}_n = \bds{L}_{n-1} + {\bds{X}_n}$. At every epoch, the step can be a no move or can move unit step along one of the $d$ coordinate axes. The corresponding values of the step are $\bds{e}_0, \bds{e}_1, \ldots, \bds{e}_d$. The walk starts at the origin, i.e., $\bds{L}_0 = \bds{e}_0$. The walk takes the first step $\bds{X}_1$ at random among all the possibilities $\bds{e}_0, \ldots, \bds{e}_d$.

Like an elephant, the walk has memory about all the past steps. At every epoch $n>1$, the walk chooses $\tau_n$, one of the past epochs based on the memories associated with them and the step $\bds{X}_{\tau_n}$, associated with the selected epoch, has its memory reinforced randomly based on its value. The memory at epoch $j$ for the step $\bds{X}_n$ is $M_n^{(j)}$. It starts growing only after the step is taken. Thus $M_n^{(j)}=0$ for $j = 0, 1, \ldots, n-1$. For any $u = 1, \ldots, n-1$, we have $P[ X_{\tau_n} = X_u | \mathcal F_{n-1}] = {M_u^{(n-1)}}/{\left( M_1^{(n-1)}+\cdots+M_{n-1}^{(n-1)}\right)}$, where $\mathcal F_n$ denotes the history up to and including time $n$. The reinforcement is summarized by an i.i.d.\ sequence of random vectors $\left((A_{0n}, A_{1n}, \ldots, A_{dn})\right)_{n\ge 1}$ with nonnegative coordinates and common finite first moments $a$. If the selected step $\bds{X}_{\tau_n} = \bds{e}_i$ for some $i = 0, 1, \ldots, d$, then the associated memory $M_{\tau_n}^{(n-1)}$ is reinforced by $A_{in}$, that is $M_{\tau_n}^{(n)} = M_{\tau_n}^{(n-1)} + A_{in}$. Other memories remain unchanged.

To choose the $n$-th step, $\bds{X}_n$ the model requires two other parameters $p, q \in (0, 1)$ other than the mean memory reinforcement $a$. If the selected step $\bds{X}_{\tau_n}$ is no move, the new step $\bds{X}_n$ is also $\bds{e}_0$ with probability $p$ or one of the directions along the axes with equal probability $(1-p)/d$. If the selected step is a move along a direction $\bds{e}_i$ for some $i = 1, \ldots, d$, then the next step $\bds{X}_n$ is the same step $\bds{e}_i$ with probability $q$ or it is no move, i.e., $\bds{e}_0$ with probability $1-q$. The step $\bds{X}_n$ is chosen independently of the entire past till epoch $(n-1)$, $\tau_n$ and $(A_{0n}, A_{1n}, \ldots, A_{dn})$. Finally, the new step $\bds{X}_n$ is given memory $1$ i.e., $M_n^{(n)} = 1$.

The asymptotic behavior of the ERW above is given by the theorem below.

\begin{theorem} Let us consider a delayed elephant random walk on the positive orthant in dimension $d$ with random reinforcement of memory, parametrized by $a$, $p$ and $q$ as described above. Then the location $\bds{L}_n$ at epoch $n$ grows linearly and
$$\frac1n \bds{L}_n \to \frac{1}{d(1-q+a(1-p))} (a(1-p), \ldots, a(1-p)) \quad \text{almost surely and in $L^1$.}$$
\end{theorem}
\begin{proof} We use the connection between ERW and urn model established in \cite{Bertoin}. We consider the memory vector of each type of step at epoch $n$, denoted by $\boldsymbol{W}^{(n)}:=(W_0^{(n)}, \ldots, W_d^{(n)})$, where, for $i = 0, 1, \ldots, d$, we have $W_i^{(n)}:=\sum_{k\leq n} M_k^{(n)}\indic{\bds{X}_k=\bds{e}_i}.$ Every time a step is taken, memory of the type corresponding to the selected past epoch is increased by the random reinforcement. Further, memory of the type corresponding to the new step is increased by 1. Thus $\left(\boldsymbol{W}^{(n)}\right)_{n\geq 0}$ behaves as an urn model of $(d+1)$ types (indexed by $0, 1, \ldots, d$) with i.i.d.\ replacement matrix sequence and mean replacement matrix
\[ \R =
\begin{pmatrix}
  1+pa & (1-p)\frac{a}{d} & \cdots & (1-p)\frac{a}{d}\\
  1-q & q+a & \cdots & q+a\\
  \vdots & \vdots & \ddots & \vdots\\
  1-q & q+a & \cdots & q+a
\end{pmatrix}.\]
Clearly, the dominant eigenvalue of $\R$ is $(1+a)$ and the corresponding left eigenvector normalized to probability is
\[\frac{1}{d(1-q+a(1-p))} \left( d(1-q), a(1-p), \ldots, a(1-p)\right).\]

Finally, note that, for $i = 1, \ldots, d$,
\[W_i^{(n)}=\sum_{k=1}^{n} \left( 1 + \sum_{l=k+1}^n A_{il} \indic{\tau_l = k} \right) \indic{\bds{X}_k=\bds{e}_i} = \bds{L}_{ni} + a \sum_{l=1}^n \indic{\bds{X}_{\tau_l}=\bds{e}_i} + \sum_{l=1}^n (A_{il}-a) \indic{\bds{X}_{\tau_l}=\bds{e}_i}.\] By Corollary~\ref{cor: iid}, we have, for the composition vector, $W_i^{(n)}/n \to (1+a) {a(1-p)}/(d(1-q+a(1-p)))$ and, for the count vector, $\sum_{l=1}^n \indic{\bds{X}_{\tau_l}=\bds{e}_i}/n \to {a(1-p)}/(d(1-q+a(1-p)))$. Finally, the last term on the right is a sum of martingale differences, which are dominated by $(A_{il}+a)$, and hence, by Proposition~\ref{prop: mg SLLN}, the last term scaled by $n$ converges to 0 almost surely. Combining the result follows.
\end{proof}

\section*{Acknowledgment}
The research of the second author was partly supported by MATRICS grant number MTR/2019/001448 from SERB, Govt.\ of India.

\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
}
\providecommand{\href}[2]{#2}


\section{Jeans models}

In our study, we employ the JASMINE code (\cite[Posacki et al., 2013]{pos13}) to construct axisymmetric stellar systems, composed of stars, a Dark Matter halo (DM) and a central Black Hole (BH), through the solution of the Jeans Equations (JEs) in cylindrical coordinates. The total stellar distribution is modeled as the superposition of $i=1,...\,,N$ different components: from a structural point of view, each stellar component is defined by an assigned axisymmetric density distribution $\rho_{*i}$, with total mass $M_{*i}$ and possible additional physical scales (such as characteristic length $r_{*i}$, truncation radius, etc.) and parameters (such as flattening $q_i$, in case of oblate ellipsoidal component). A DM halo of density distribution $\rho_\mathrm{h}$ and total mass $M_\mathrm{h}$, and a central BH of mass $M_\mathrm{BH}$, can be included. The total gravitational potential due to the stars and the DM halo is numerically computed through the solution of the Poisson equation.

Each stellar density component of the resulting models, with symmetry axis aligned with the $z$-axis, is assumed to be described by a two-integral distribution function $f_i(E,J_\mathrm{z})$. As well known, under these assumptions, $\overline{v_\mathrm{R}}_i=\overline{v_\mathrm{z}}_i=0$, and the only non-zero rotational velocity can be in the azimuthal direction, $\overline{v_\varphi}_i$. In addition, the velocity dispersion tensor is diagonal and aligned with the coordinate system, and it can have only azimuthal anisotropy (while $\sigma_{\mathrm{R}i}=\sigma_{\mathrm{z}i}\equiv \sigma_i$). It follows that the JEs for the stellar distribution $\rho_{*i}$, embedded in a total gravitational potential $\Phi_\mathrm{tot}=\sum_i\phi_{*i}+\phi_\mathrm{h}+\phi_\mathrm{BH}$ (sum of the contributions of stars, DM and BH), are given by:
\begin{equation} \label{eq:vertJEs}
 \displaystyle{\frac{\partial \big(\rho_{*i}\sigma_i^2\big)}{\partial z} = -\rho_{*i} \frac{\partial \Phi_\mathrm{tot}}{\partial z}} \,,
\end{equation}
\begin{equation} \label{eq:radJEs}
 \displaystyle{\frac{\partial \big(\rho_{*i}\sigma_i^2\big)}{\partial R} = \frac{\rho_{*i} \Big(\overline{v_\varphi^2}_i-\sigma_i^2\Big)}{R} -\rho_{*i} \frac{\partial \Phi_\mathrm{tot}}{\partial R}} \,,
\end{equation}
(see e.g. \cite[Binney \& Tremaine, 2008]{b&t08}).

In order to split the azimuthal kinetic energy ($\overline{v_\varphi^2}_i$) between streaming and random motions, we adopt the $k$-decomposition (\cite[Satoh, 1980]{satoh80}), with the introduction of a constant parameter $k_i$, in principle different for each stellar component, so that
\begin{equation} \label{eq:satoh}
 \overline{v_\varphi}_i^2 = k_i^2 \Big(\overline{v_\varphi^2}_i-\sigma_i^2\Big) \,, \qquad
 \sigma_{\varphi i}^2 \equiv \overline{v_\varphi^2}_i - \overline{v_\varphi}_i^2  = \sigma_i^2 + \big(1-k_i^2\big) \Big(\overline{v_\varphi^2}_i-\sigma_i^2\Big) \,.
\end{equation}
Usually $0\leq k_i \leq 1$, but more general choices are possible (see e.g. \cite[Ciotti \& Pellegrini, 1996]{ciot-pel96}). If $k_i=1$, the component into consideration is an isotropic rotator ($\sigma_{\varphi i}=\sigma_i$) and it is fully rotationally supported; if $k_i=0$, there is no net rotation ($\overline{v_\varphi}_i=0$), so the component is supported only by the azimuthal anisotropy of the velocity dispersion tensor. In order to introduce anisotropy between the radial and vertical velocity dispersions (implicitly relaxing the assumption of a two-integral distribution function), i.e. to have $\sigma_{\mathrm{R}i}\neq\sigma_{\mathrm{z}i}$, the JASMINE code can be easily modified to include, for example, the Cappellari (\cite[2008]{capp08}) anisotropy.

Once the kinematic fields are derived, we construct their mass and luminosity-weighted projections along a given line-of-sight direction, and we obtain the corresponding 2D maps. In turn, these maps can be compared with the analogous maps of real stellar systems, as obtained for instance from Integral Field Spectroscopy.

The linearity of the JEs, with respect to $\rho_{*i}$ and $\Phi_\mathrm{tot}$, allows for a significant amount of \textit{scaling} of the solutions. Indeed, since the general solutions are built by adding the solutions pertaining to each stellar component $\rho_{*i}$, in the gravitational field produced by $\Phi_\mathrm{tot}$, it is possible to considerably reduce the amount of numerical work required. In practice, for a given set of mass components, we numerically solve the JEs for a "basis" of solutions (given by the solutions for \textit{mass-normalised} $\rho_{*i}$, $\phi_{*i}$, $\phi_\mathrm{h}$, $\phi_\mathrm{BH}$) and we project them. Then we exploit the linearity previously discussed, in addition with the linearity of the Satoh decomposition (eq. \ref{eq:satoh}) and of the projection operators, in order to recombine them in the post-processing tuning, when only the weight (mass, luminosity, anisotropy) coefficients are modified. This allows to build \textit{families} of models, by running the JASMINE code only once.

\section{Applications}

The applications of the analytical and numerical framework briefly illustrated in the previous section are numerous. It is now possible to reproduce, for example, different stellar kinematic components, such as counter rotating disks and kinematically decoupled cores. Furthermore, we can obtain luminosity-weighted total properties, by assuming a constant mass-to-light ratio different for each component. From a stellar population point of view, indeed, each component can be considered as a Simple Stellar Population, characterised by a certain age, metallicity, IMF, mass-to-light ratio (in agreement with Evolutionary Population Synthesis models, e.g. \cite[Maraston, 2005]{mar05}). The combination of different populations allows, for example, to study the systematic effects of IMF variations and mass-to-light ratio gradients on luminosity-weighted properties. The natural field of application of our study is the modeling of Early-Type Galaxies (ETGs), however the same approach can be also applied to modeling dwarf Spheroidal galaxies and Globular Clusters.

For illustrative purposes, we present a model of an elliptical galaxy, composed of two stellar populations, a DM halo and a central supermassive BH, whose parameters are reported in Table \ref{tab:parameters}. The total and first stellar components are modeled by ellipsoidal Jaffe density profiles, with a different flattening, mass and scale length. The second component is defined as the \textit{difference} of the two chosen distributions, $\rho_{*2}(R,z) = \rho_*(R,z) - \rho_{*1}(R,z)$, so that $M_{*2} = M_* - M_{*1}$, and in general it is \textit{not} ellipsoidal. In order to have a non-negative $\rho_{*2}(R,z)$ over the whole space, a \textit{non-negativity condition} is derived, constraining $q$, $q_1$, $M_{*1}/M_*$ and $r_{*1}/r_*$. In Figure \ref{fig:rho}, the stellar density profiles on the equatorial plane are shown. The DM halo is modeled by a spherical NFW density profile (\cite[Navarro, Frenk \& White, 1996]{NFW}), with given mass, scale radius and concentration. In Figure \ref{fig:sigma}, the radial (and vertical) velocity dispersions are shown for the two stellar components and for the total stars, as an example of the dynamics of a multi-component stellar system.
\begin{table}
	\centering
	\begin{tabular}{cccccccccc}
		\hline
		\multicolumn{2}{c}{Flattenings} &  \multicolumn{4}{c}{Masses} &
		\multicolumn{2}{c}{Scale lengths} & \multicolumn{2}{c}{Satoh parameters} \\
		$q$ & $q_1$ & $M_*$ & $M_{*1}/M_*$ & $M_\mathrm{h}/M_*$ & $M_\mathrm{BH}/M_*$ & $r_*$ & $r_{*1}/r_*$ & $k_1$ & $k_2$ \\
		\hline
		0.9 & 0.8 & $7.7 \cdot \mathrm{10^{10}\,M_\odot}$ & 0.3 & $10^2$ & $10^{-3}$ & 4.3 $\mathrm{kpc}$ & 0.5 & 1.0 & 0.0 \\
		\hline
	\end{tabular}
	\caption{Physical parameters of an illustrative model of an elliptical galaxy. This model is composed of two stellar populations, a DM halo and a central supermassive BH.}
	\label{tab:parameters}
\end{table}

\begin{figure}[ht]
 \centering
 \includegraphics[scale=0.8]{t11_rho_r}
 \caption{Density profiles on the equatorial plane of the two stellar components and total stars. The stellar component $\rho_{*2}$ is given by the difference of the two ellipsoidal Jaffe stellar models $\rho_*$ and $\rho_{*1}$. See Table \ref{tab:parameters} for the physical parameters adopted for this galaxy model.}
 \label{fig:rho}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.55]{t11_sigma}
	\caption{Radial (and vertical) components of the velocity dispersion of the total stars $\rho_*$ (left), and the two stellar components $\rho_{*1}$ (center) and $\rho_{*2}$ (right). In each panel, the contours represent the iso-densities of the corresponding component. See Table \ref{tab:parameters} for the physical parameters adopted for this galaxy model.}
	\label{fig:sigma}
\end{figure}

\newpage


\section{Introduction}

Neural machine translation (NMT), the current prevalent approach to automatic translation, is known to require large amounts of parallel training sentences and an extensive amount of training time on dedicated hardware.
The total training time significantly increases, especially when training strong baselines, searching for best hyperparameters or training multiple models for various language pairs.

\perscite{schwartz2019green} analyzed 60 papers from top AI conferences and found out that 80\% of them target accuracy over efficiency, and only a small portion of papers argue for a new efficiency result. They also noted that the increasing financial cost of the computations could make it difficult for researchers to engage in deep learning research or limit training strong baselines. Furthermore, increased computational requirements have also an environmental cost.
\perscite{strubell2019energy} estimated that training a single Transformer ``big'' model produces 87 kg of $\text{CO}_2${} and that the massive Transformer architecture parameter search produced 298 tonnes of $\text{CO}_2${}.\footnote{The paper reports numbers based on the U.S. energy mix.}

However, a lot of research has been already invested into cutting down the long training time by the design of NMT model architectures, promoting self-attentive \parcite{vaswani:2017} or convolutional \parcite{gehring2017convolutional} over recurrent ones \parcite{bahdanau:2014:corr} or the implementation of heavily optimized toolkits \parcite{mariannmt}.

In this paper, we propose a novel view on re-using already trained ``parent'' models without the need to prepare a parent model in advance or modify its training hyper-parameters.
Furthermore, we propose a second method based on a vocabulary transformation
technique that makes even larger improvements, especially for languages using an
alphabet different from the re-used parent model. Our transfer learning approach
leads to better performance as well as faster convergence speed of the
``child'' model compared to training the model from scratch. We document that
our methods are not restricted only to low-resource languages, but they can be used even for high-resource ones.

Previous transfer learning techniques \parcite{neubig-hu:2018:EMNLP,kocmi-adaption} rely on a shared vocabulary between the parent and child models. As a result, these techniques separately train parent model for each different child language pair. In contrast, our approach can re-use one parent model for multiple various language pairs, thus further lowering the total training time needed.

In order to document that our approach is not restricted to parent models trained by us, we re-use parent model trained by different researchers: we use the winning model of WMT 2019 for Czech-English language pair \parcite{popel-EtAl:2019:WMT}.

The paper is organized as follows: \cref{sec:method_desc} describes the method of Direct Transfer learning, including our improvement of vocabulary transformation.
\cref{sec:experiments} presents the model, training data, and our experimental setup.
\cref{sec:results} describes the results of our methods followed by the analysis in \cref{sec:analysis}. Related work is summarized in \cref{sec:related} and we conclude the discussion in \cref{sec:conclusion}.

\section{Transfer Learning}
\label{sec:method_desc}

In this work, we present the use of transfer learning to reduce the training time and improve the performance in comparison to training from random initialization even for high-resource language pairs.

Transfer learning is an approach of using training data from a related task to improve the accuracy of the main task in question \parcite{tan2018survey}. One of the first transfer learning techniques in NMT was proposed by \perscite{zoph-EtAl:2016:EMNLP2016}. They used word-level NMT and froze several model parts, especially embeddings of words that are shared between parent and child model.

We build upon the work of \perscite{kocmi-adaption}, who simplified the transfer learning technique thanks to the use of subword units \parcite{wu2016google} in contrast to word-level NMT transfer learning \parcite{zoph-EtAl:2016:EMNLP2016} and extended the applicability to unrelated languages.

Their only requirement, and also the main disadvantage of the method, is that the vocabulary has to be shared and constructed for the given parent and child languages jointly, which makes the parent model usable only for the particular child language pair. This substantially increases the overall training time needed to obtain the desired NMT system for the child language pair.

The method of \perscite{kocmi-adaption} consists of three steps: (1) construct the vocabulary from both the parent and child corpora, (2) train the parent model with the shared vocabulary until convergence, and (3) continue training on the child training data.

\perscite{neubig-hu:2018:EMNLP} call such approaches warm-start, where we use the child language pair to influence the parent model. In our work, we focus on the so-called cold-start scenario, where the parent model is trained without a need to know the language pair in advance. Therefore we cannot make any modifications of the parent training to better handle the child language pair.
The cold-start transfer learning is expected to have slightly worse performance than the warm-start approach. However, it allows reusing one parent model for multiple child language pairs, which reduces the total training time in comparison to the use of warm-start transfer learning.

We present two approaches: Direct Transfer that ignores child-specific vocabulary altogether; and Transformed Vocabulary, which modifies vocabulary of the already trained parent. Thus, one parent model can be used for multiple child language pairs.

\subsection{Direct Transfer}

Direct Transfer can be seen as a simplification of \perscite{kocmi-adaption}. We ignore the specifics of the child vocabulary and train the child model using the parent vocabulary. We suppose that the subword vocabulary can handle the child language pair, although it is not optimized for it.

We take an already trained model and use it as initialization for a child model using a different language pair. We continue the training process without any change to the vocabulary or hyper-parameters. This applies even to the training parameters, such as the learning rate or moments.

This method of continued training on different data while preserving hyper-parameters is used under the name ``continued training'' or ``fine-tuning'' \parcite{hinton2006reducing,micelibarone2017EMNLP2017}, but it is mostly used as a domain adaptation within a given language pair.

Direct Transfer relies on the fact that the current NMT uses subword units instead of words.
The subwords are designed to handle unseen words or even characters, breaking the input into shorter units, possibly down to individual bytes as implemented, for example, by Tensor2Tensor \parcite{tensor2tensor}. This property ensures that the parent vocabulary can, in principle, serve for any child language pair, but it can be highly suboptimal, segmenting child words into too many subwords.

We present an example of a Russian phrase and its segmentation based on English-Czech or English-Russian vocabulary in \cref{tab:vocab_variant_usage}. When using child-specific vocabulary, the segmentation works as expected, splitting the phrase into three tokens. However, when we use a vocabulary that contains only the Cyrillic alphabet\footnote{This happened solely due to noise in the Czech-English parent training data.} and not many longer sequences of characters, the sentence is split into 13 tokens. We can notice that English-Czech wordpiece vocabulary is missing a character \foreignlanguage{russian}{``Л''}, thus it breaks it into the byte representation ``\verb|\1051;|''.

\begin{table}[t]
\centering
\begin{tabular}{l|cc|cc}
\multicolumn{1}{c}{} & \multicolumn{2}{c|}{Child-specific} & \multicolumn{2}{c}{EN-CS vocab.}  \\
\multicolumn{1}{r}{Avg. \# per:} & Sent. & Word & Sent. & Word \\
\hline
Odia     & 95.8 & 3.7 & 496.8 & 19.1\\
Estonian & 26.0 & 1.1 &  56.2 &  2.3\\
Finnish  & 22.9 & 1.1 &  55.9 &  2.6\\
German   & 27.4 & 1.3 &  55.4 &  2.5\\
Russian  & 33.3 & 1.3 & 134.9 &  5.3\\
French   & 42.0 & 1.6 &  65.7 &  2.5\\
\end{tabular}
\caption{Average number of tokens per sentence (column ``Sent.'') and average
number of tokens per word (column ``Word'') when the training corpus is segmented
by child-specific or parent-specific vocabulary. ``Child-specific'' represents
the effect of using vocabulary customized for examined language. ``EN-CS''
corresponds to the use of English-Czech vocabulary.}
\label{tab:tokenized_dataset}
\end{table}

\begin{figure}[t]
\centering
\begin{tabular}{l|ll}
 & Segmented sentence  \\
\hline
Original & \foreignlanguage{russian}{Сьерра-Леоне}  \\
EN-RU & \foreignlanguage{russian}{Сьерра\_█-\_█Леоне\_}  \\
EN-CS  & \foreignlanguage{russian}{С█ь█ер█ра█\_█-\_█\textbackslash█10█51█;█е█о█не\_} \\
\end{tabular}
\caption{Illustration of segmentation of Russian phrase (gloss: Sierra Leone) with English-Czech and English-Russian vocabulary from our experiments. The character █ represents splits.}
\label{tab:vocab_variant_usage}
\end{figure}

We examine the influence of parent-specific vocabulary on the training dataset of the child. \cref{tab:tokenized_dataset} documents the segmenting effect of different vocabularies. If we compare the child-specific and parent-specific (``EN-CS'') vocabulary, the average number of tokens per sentence or per word increases more than twice.
For example, German has twice as many tokens per word compared to its child-specific vocabulary, and Russian has four times more tokens due to Cyrillic. Odia is affected even more.

Thus, we see that ignoring the vocabulary mismatch introduces a problem for NMT models in the form of an increasing split ratio of tokens. As expected, this is most noticeable for languages using different scripts.

\begin{algorithm}[t]
\begin{center}
\framebox{
\parbox{.9\columnwidth}{
 \SetKwInput{KwData}{Input}
 \KwData{Parent vocabulary (an ordered list of parent subwords) and the training corpus for the child language pair.}
 Generate child-specific vocabulary with the maximum number of subwords equal to the parent vocabulary size\;
 \For{subword S in parent vocabulary}{
  \eIf{S in child vocabulary}{
  continue\;
  }{
  Replace position of S in the parent vocabulary with the first unused child subword not contained in the parent\;
  }
 }
 \KwResult{Transformed parent vocabulary}
}}
 \vspace{1em}
 \caption{Transforming parent vocabulary to contain child subwords and match positions for subwords common for both of language pairs.}
 \label{algorithm}
\end{center}
\end{algorithm}
\subsection{Vocabulary Transformation}
\label{sec:vocabulary_transformation}

Using parent vocabulary roughly doubles the number of subword tokens per word, as we showed in the previous section. This problem would not happen with child-specific vocabulary. However, we are using an already trained parent with its vocabulary. Therefore, we propose a vocabulary transformation method that replaces subwords in the parent wordpiece \parcite{wu2016google} vocabulary with subwords from the child-specific vocabulary.

NMT models associate each vocabulary item with its vector representation
(embedding). When transferring the model from the parent to the child, we decide
which subwords should preserve their embedding as trained in the parent model
and which embeddings should be remapped to new subwords from the child
vocabulary. The goal is to preserve embeddings of subwords that are contained in both parent and child vocabulary. In other words, we reuse embeddings of subwords common to both parent and child vocabularies and reuse the vocabulary entries of subwords not occurring in the child data for other, unrelated, subwords that the child data need. Obviously, the embeddings for these subwords will need to be retrained.

\begin{table*}[t]
\small
\centering
\begin{tabular}{l|r|rrr}
Language pair  & Pairs  & Training set    & Development set    & Test set \\
\hline
EN - Odia & 27k & \perscite{parida2018odiacorpus} & \perscite{parida2018odiacorpus} & \perscite{parida2018odiacorpus} \\
EN - Estonian & 0.8M & Europarl, Rapid & WMT dev 2018 & WMT 2018 \\
EN - Finnish & 2.8M & Europarl, Paracrawl, Rapid  & WMT 2015 & WMT 2018\\
EN - German & 3.5M & Europarl, News commentary, Rapid & WMT 2017 & WMT 2018\\
EN - Russian & 12.6M & News Commentary, Yandex, and UN Corpus & WMT 2012 & WMT 2018\\
EN - French & 34.3M & \begin{tabular}{@{}r@{}}Commoncrawl, Europarl, Giga FREN, \\ News commentary,  UN corpus\end{tabular} & WMT 2013 & WMT dis. 2015\\
\end{tabular}
\caption{Corpora used for each language pair. The names specify the corpora from WMT 2018 News Translation Task data. Column ``Pairs'' specify the total number of sentence pairs in training data.}
\label{tab:dataset_sources}
\end{table*}

Our Transformed Vocabulary method starts by constructing the child-specific vocabulary with the size equal to the parent vocabulary size (the parent model is trained, thus it has a fixed number of embeddings). Then, as presented in Algorithm~\ref{algorithm}, we generate an ordered list of child subwords, where subwords known to the parent vocabulary are on the same positions as in the parent vocabulary, and other subwords are assigned arbitrarily to places where parent-only subwords were stored.

We experimented with several possible mappings between the parent and child vocabulary. We tried to assign subwords based on frequency, by random assignment, or based on Levenshtein distance of parent and child subwords. However, all the approaches reached comparable performance; neither of them significantly outperformed the others. One exception is when assigning all subwords randomly, even those that are shared between parent and child. This method leads to worse performance, having several BLEU points lower than other approaches.
Another approach would be to use pretrained subword embeddings similarly as proposed \perscite{kim2019effective}. However, in this paper, we focus on showing, that transfer learning can be as simple as not using any modifications at all.

\section{Experiments}
\label{sec:experiments}

In this section, we first provide the details of the NMT model used in our experiments and the examined set of language pairs. We then discuss the convergence and a stopping criterion and finally present the results of our method for recycling the NMT model as well as improvements thanks to the vocabulary transformation.

\subsection{Parent Model and its Training Data}
\label{sec:model_setting}

In order to document that our method functions in general and is not restricted to our laboratory setting, we do not train the parent model ourselves. Instead, we recycle two systems trained by \perscite{popel-EtAl:2019:WMT}, namely the English-to-Czech and Czech-to-English winning models of WMT 2019 News Translation Task. It is important to note, that we use two parent models and for experiments we always use the parent model with English on the same side, e.g. English-to-Russian child has English-to-Czech as a parent. We leave experimenting with different parents or various combinations for future works, because the goal of this work is to make approach most simple.

We decided to use this model for several reasons. It is trained to translate into Czech, a high-resource language that is dissimilar from any of the languages used in this work.\footnote{The linguistically most similar language of our language selection is Russian, but we do not transliterate Cyrillic into Latin script. Therefore, the system cannot associate similar Russian and Czech words based on appearance.} At the same time, it is trained using the state-of-the-art Transformer architecture as implemented in the Tensor2Tensor framework.\footnote{\url{https://github.com/tensorflow/tensor2tensor}
} \parcite{tensor2tensor}. We use Tensor2Tensor in version 1.8.0.

The model is described in \perscite{popel:2018:WMT}. It is based on the ``Big GPU Transformer'' setup as defined by \perscite{vaswani:2017} with a few modifications. The model uses reverse square root learning rate decay with 8000 warm-up steps and a learning rate of 1. It uses the Adafactor optimizer, the batch size of 2900 subword units, disabled layer dropout.

Due to the memory constraints, we drop training sentences longer than 100 subwords.
We use child hyper-parameter setting equal to the parent model. However, some hyper-parameters like learning rate, dropouts, optimizer, and others could be modified for the training of the child model. We leave these experiments for future work.

We train models on single GPU GeForce 1080Ti with 11GB memory. In this setup, 10000 training steps take on average approximately one and a half hours. \perscite{popel-EtAl:2019:WMT} trained the model on 8 GPUs for 928k steps, which means that on the single GPU, the parent model would need at least 7424k steps, i.e. more than 45 days of training.

In our experiments, we train all child models up to 1M steps and then take the model with the best performance on the development set. Because some of the language pairs, especially the low-resource ones, converge within first 100k steps, we use a weak early stopping criterion that stops the training whenever there was no improvement larger than 0.5\% of maximal reached BLEU over the past 50\% of training evaluations (minimum of training steps is 100k). This stopping criterion makes sure that no model is stopped prematurely.

\def\rlap{ $\ddagger$*}{\rlap{ $\ddagger$*}}
\begin{table*}[t]
\centering
\begin{tabular}{c|rr|rr|r@{~~~~~~~}rrr}
Language pair & \multicolumn{2}{c}{Baseline} & \multicolumn{2}{|c|}{Direct Transfer} & \multicolumn{4}{c}{Transformed Vocab}\\
                   &  BLEU & Steps & BLEU  &Steps &  BLEU & Steps & $\Delta$ BLEU & Speed-up\\
\hline
English-to-Odia    &  3.54 &   45k &  0.26 &  47k & \bf 6.38\rlap{ $\ddagger$*}{} & \bf 38k & 2.84 & 16 \%  \\
English-to-Estonian& 16.03 &   95k & \bf 20.75\significant{} & \bf 75k & 20.27\significant{} & \bf 75k & 4.24 & 21 \%\\
English-to-Finnish & 14.42 &  420k & 16.12\significant{} & \bf 255k & \bf 16.73\rlap{ $\ddagger$*}{} & 270k & 2.31 & 36 \% \\
English-to-German  & 36.72 &  270k & 38.58\significant{} & 190k & \bf 39.28\rlap{ $\ddagger$*}{} & \bf 110k & 2.56 & 59 \%\\
English-to-Russian & 27.81 & 1090k & 27.04 & 630k & \bf 28.65\rlap{ $\ddagger$*}{} & \bf 450k & 0.84 & 59 \% \\
English-to-French  & 33.72 &  820k & 34.41\significant{} & \bf 660k & \bf 34.46\significant{} & 720k & 0.74 & 12 \%\\
\hline
Estonian-to-English & 21.07 & 70k & 24.36\significant{} & \bf 30k  & \bf 24.64\rlap{ $\ddagger$*}{}  & 60k & 3.57 & 14 \%\\
Russian-to-English  & 30.31 & 980k & 23.41 & \bf 420k & \bf 31.38\rlap{ $\ddagger$*}{}  & 700k & 1.07 & 29 \%\\
\end{tabular}
\caption{Translation quality and training time. ``Baseline'' is trained from
scratch with its own vocabulary and child corpus only. ``Direct Transfer'' is
initialized with parent model using the parent vocabulary and continues
training. ``Transformed Vocab'' has the same initialization but merges the
parent and child vocabulary as described in
\cref{sec:vocabulary_transformation}. Best score and lowest training time in
each row in bold. The statistical significance is computed against the baseline ($\ddagger${}) or against ``Direct Transfer'' (*).
Last two columns show improvements of Transformed Vocabulary in comparison to the baseline.
}
\label{tab:main_results}
\end{table*}

\subsection{Studied Language Pairs}
\label{sec:studied_languages}

We use several child language pairs to show that our approach is useful for various sizes of corpora, language pairs, and scripts. To cover this range of situations, we select languages in \cref{tab:dataset_sources}. Future works could focus also on languages outside from Indo-European family, such as Chinese.

Another decision behind selecting these language pairs is to include language pairs reaching various levels of translation quality. This is indicated by automatic scores of the baseline setups ranging from 3.54 BLEU (English-to-Odia) to 36 BLEU (English-to-German)\footnote{The systems submitted to WMT 2018 for English-to-German translation have better performance than our baseline due to the fact, that we decided not to use Commoncrawl, which artificially made English-German parallel data less resourceful.}, see \cref{tab:main_results}.

The sizes of corpora are in \cref{tab:dataset_sources}. The smallest language pair is English-Odia, which uses the Brahmic writing script and contains only 27 thousand training pairs. The largest is the high-resource English-French language pair.

For most of the language pairs, we use training data from WMT \parcite{wmt_2018}.\footnote{\url{http://www.statmt.org/wmt18/}} We use the training data without any preprocessing, not even tokenization.\footnote{While the recommended best practice in past WMT evaluations was to use Moses tokenizer. It is not recommended for Tensor2Tensor with its build-in tokenizer any more.} See \cref{tab:dataset_sources} for the list of used corpora for each language pair. For some languages, we have opted out from using all available corpora in order to experiment on languages containing various magnitudes of parallel sentences.

For high-resource English-French language pair, we perform a corpora cleaning using language detection Langid.py \cite{lui-baldwin:2012:Demo}.
We drop all sentences that are not recognized as the correct language. It removes 6.5M (15.9
\%) sentence pairs from the English-French training corpora.

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

All reported results are calculated on the test data and evaluated with SacreBLEU \parcite{sacrebleu}. The results are in \cref{tab:main_results}. We discuss separately the training time, automatically assessed translation quality using the parent and the Transformed Vocabulary, and comparison to \perscite{kocmi-adaption} in the following sections.

Baselines use the same architecture, and they are trained solely on the child training data with the use of child-specific vocabulary. We compute statistical significance with a paired bootstrap resampling \parcite{bootstrap-koehn:2004}. We use 1000 samples and a confidence level of 0.05. Statistically significant improvements are marked by $\ddagger${}.

\subsection{Direct Transfer Learning}

First, we compare the Direct Transfer learning in contrast to the baseline. We see that Direct Transfer learning is significantly better than the baseline in both translation directions in all cases except for Odia and Russian, which we will discuss later. We get improvements for various language types, as discussed in \cref{sec:studied_languages}. The largest improvement is of 4.72 BLEU for the low-resource language pair of Estonian-English, but we also get an improvement of 0.69 BLEU for the high-resource pair French-English.

The results are even more surprising when we take into account the fact that the model uses the parent vocabulary, and it is thus segmenting words into considerably more subwords. This suggests that the Transformer architecture generalizes very well to short subwords.
However, the worse performance of English-Odia and English-Russian can be attributed to the different writing script. The Odia script is not contained in the parent vocabulary at all, leading to segmenting of each word into individual bytes, the only common units with the parent vocabulary.
Therefore, to avoid problems with filtering, we increase the filtering limit of long sentences during training from 100 to 500 subwords for these two language pairs (see \cref{sec:model_setting}).

\subsection{Results with Transformed Vocabulary}

As the results in \cref{tab:main_results} confirm, Transformed Vocabulary successfully tackles the problem of the child language using a different writing script.
We see ``Transformed Vocab'' delivering the best performance for all language
pairs except for English-to-Estonian, significantly improving over baseline and even
over ``Direct Transfer'' in most cases.

\subsection{Training Time}

In the introduction, we discussed that recent development in NMT focuses mainly on the performance over efficiency \parcite{schwartz2019green}. Therefore, in this section, we discuss the amount of training time required for our method to converge. We are reporting the number of updates (i.e. steps) needed to get the model used for evaluation.\footnote{Another possibility would be to report wall-clock time. However, that is influenced by server load and other factors. The number of steps is better for the comparison as long as the batch size stays the same across experiments.}

We see in \cref{tab:main_results} that both our methods converged in a lower number of steps than the baseline. For the Transformed Vocabulary method, we get a speed-up of 12--59 \%. The reduction in the number of steps is most visible in English-to-German and English-to-Russian.
It is important to note that the number of steps to the convergence is not precisely comparable, and some tolerance must be taken into account. It is due to the fluctuation in the training process. However, in neither of our experiments, Transformed Vocabulary is slower than baseline. Thus we conclude that our Transformed Vocabulary method takes fewer training steps to finish training than training a model from scratch.

\begin{table}[t]
\centering
\small
\begin{tabular}{ll|rrr}
&Language &  & Transf. &     Warm\\
&pair     & Baseline & vocab & Start \\
\hline
\multirow{4}{*}{\rot{BLEU}}&To Estonian   & 16.03  & 20.27 & \bf 20.75 \\
&To Russian    & 27.81 & 28.65 &\bf 29.03\significant{}  \\
&From Estonian & 21.07 & 24.64 &\bf 26.00\significant{}  \\
&From Russian  & 30.31 & \bf 31.38  & 31.15   \\
\hline
\multirow{4}{*}{\rot{Steps}}&To Estonian   & 95k  & \bf 75k & 735k  \\
&To Russian    & 1090k  & \bf 450k& 1510k \\
&From Estonian & 70k & \bf 60k& 700k   \\
&From Russian  & 980k& \bf 700k & 1465k  \\
\end{tabular}
\caption{Comparison of our Transformed Vocabulary method with \perscite{kocmi-adaption} (abridged as ``Warm Start''). The top half of table compares results in BLEU, the bottom half the number of steps needed to convergence. Steps of \perscite{kocmi-adaption} method are reported as the sum of parent and child training, due to the  nature of the method.
}
\label{tab:transfer_compar}
\end{table}

\begin{figure*}
\begin{tikzpicture}
\begin{axis}[
    ybar,
    bar width=.2cm,
    width=8.5cm,
    height=5cm,
    enlargelimits=0.15,
    legend style={at={(1,0)},
      anchor=south east,font=\scriptsize,legend columns=2},
    ylabel={BLEU},

    symbolic x coords={Direct T., Transf. V., Direct T, Transf. V},
    xtick=data,

    title=English-to-Estonian,
    xlabel style={text width=8.5cm},
    xlabel=\hspace{33pt} Freeze only one \hspace{17pt} Freeze all but one
    ]
\addplot +[black, fill=yellow] coordinates {(Direct T.,20.94)(Transf. V.,19.18)(Direct T,13.99)(Transf. V,16.47)};
\addplot +[black, fill=cyan, postaction={pattern=north east lines}] coordinates {(Direct T.,21.09)(Transf. V.,20.01)(Direct T,16.13)(Transf. V,14.22)};
\addplot +[black, fill=green, postaction={pattern=north west lines}] coordinates {(Direct T.,17.94)(Transf. V.,18.33)(Direct T,20.62)(Transf. V,17.68)};
\addplot +[black, fill=red] coordinates {(Direct T.,20.83)(Transf. V.,20.40)(Direct T,19.70)(Transf. V,17.33)};
\addplot +[black, fill=white] coordinates {(Direct T.,20.75)(Transf. V.,20.27)(Direct T,20.75)(Transf. V,20.27)};
\node[] at (axis cs: Direct T.,21.5) {\circled{1}};
\node[] at (axis cs: Transf. V.,21.5) {\circled{2}};
\node[] at (axis cs: Direct T,21.5) {\circled{3}};
\node[] at (axis cs: Transf. V,21.5) {\circled{4}};
\end{axis}

\begin{scope}[shift={(7.7,0)}]
\begin{axis}[
    ybar,
    bar width=.2cm,
    width=8.5cm,
    height=5cm,
    enlargelimits=0.15,
    legend style={at={(0,0)},
      anchor=south west,font=\scriptsize,legend columns=2},

    symbolic x coords={Direct T., Transf. V., Direct T, Transf. V},
    xtick=data,

    title=Estonian-to-English,
    xlabel style={text width=8.5cm},
    xlabel=\hspace{33pt} Freeze only one \hspace{17pt} Freeze all but one
    ]
\legend{Embedings,Encoder,Decoder,Attention,Train all}
\addplot +[black, fill=yellow] coordinates {(Direct T,21.30)(Transf. V,22.92)(Direct T.,24.96)(Transf. V.,21.54)};
\addplot +[black, fill=cyan, postaction={pattern=north east lines}] coordinates {(Direct T,25.23)(Transf. V,22.24)(Direct T.,21.01)(Transf. V.,20.38)};
\addplot +[black, fill=green, postaction={pattern=north west lines}] coordinates {(Direct T,16.86)(Transf. V,12.11)(Direct T.,25.29)(Transf. V.,23.97)};
\addplot +[black, fill=red] coordinates {(Direct T,22.93)(Transf. V,20.17)(Direct T.,24.39)(Transf. V.,23.46)};
\addplot +[black, fill=white] coordinates {(Direct T,24.36)(Transf. V,24.64)(Direct T.,24.36)(Transf. V.,24.64)};
\node[] at (axis cs: Direct T.,26) {\circled{5}};
\node[] at (axis cs: Transf. V.,26) {\circled{6}};
\node[] at (axis cs: Direct T,26) {\circled{7}};
\node[] at (axis cs: Transf. V,26) {\circled{8}};
\end{axis}
\end{scope}
\end{tikzpicture}
\caption{Child BLEU scores when trained with some parameters frozen. The left plot shows English-to-Estonian and the right is Estonian-to-English. In both plots, the first two groups are experiments where one component is frozen and the second two are when all components but one are frozen.}
\label{tab:freezing_analysis}
\end{figure*}

\subsection{Comparison to \perscite{kocmi-adaption}}

We replicated the experiments of \perscite{kocmi-adaption} with the identical
framework and hyperparameter setting in order to compare their method to ours.
We experiment with Estonian-English and Russian-English language pair in both
translation directions. Their approach needs an individual parent for every
child model, so we train four models: two English-to-Czech and two
Czech-to-English on the same parent training data as \perscite{kocmi-adaption}.
All vocabularies contain 32k subwords.
We compare their method with our Transformed Vocabulary. Furthermore, the results of Direct Transfer in \cref{tab:main_results} are also comparable with this experiment.

In \cref{tab:transfer_compar}, we see that their method reaches a slightly better performance in three translation models, where English-to-Russian and Estonian-to-English are significantly ($\ddagger${}) better than Transformed Vocabulary technique; the other two are on par with our method, which is understandable. The Transformed Vocabulary cannot outperform the warm-start technique since the warm-start parent model has the advantage of being trained with the vocabulary prepared for the investigated child.

However, when we compare the total number of steps needed to reach the performance, both our approaches are significantly faster than \perscite{kocmi-adaption}. The most substantial improvements are roughly ten times faster for Estonian-to-English, and the smallest difference for English-to-Russian is two times faster.
This is mostly because their method first needs to train the parent model that is specialized for the child, while our method can directly re-use any already trained model. Moreover, we can see that their method is even slower than the baseline model.

\section{Analysis by Freezing Parameters}
\label{sec:analysis}

To discover which transferred parameters are the most helpful for the child
model and which need to be changed the most, we follow the analysis used by
\perscite{thompson-EtAl:2018:WMT}: When training the child, we freeze some of
the parameters.

Based on the internal layout of the Transformer model in Tensor2Tensor, we divide the model into four components. (i) Word embeddings (shared between encoder and decoder) map each subword unit to a dense vector representation.
(ii) The encoder component includes all the six feed-forward layers converting the input sequence to the deeper representation. (iii) The decoder component consists again of six feed-forward layers preparing the choice of the next output subword unit. (iv) The multi-head attention is used throughout encoder and decoder, as self-attention layers interweaved with the feed-forward layers.

We run two sets of experiments: either freeze only one out of the four components and leave the rest of the model updating or freeze everything but the examined component. We also test it on two translation directions: \langpair{English}{Estonian} in the left hand part of \cref{tab:freezing_analysis} and \langpair{Estonian}{English} in the right hand part. In both cases, English-Czech (in the corresponding direction, i.e. with English on the correct side) serves as the parent.
We discuss individual components separately, indexing the experiments
\circled{1} to \circled{8}.

Similarly to \perscite{thompson-EtAl:2018:WMT} in domain adaptation, we observe
that parent embeddings serve well in Direct Transfer, freezing them has a
minimal impact compared to the baseline in \circled{1} and \circled{5}.
The frozen embeddings in Transformed
Vocabulary (\circled{2}, \circled{6}) results in significant performance drops
which can be attributed to the arbitrary assignment of embeddings to new subwords.

The comparison of all but embeddings frozen in \circled{4} and
\circled{8} (Transformed Vocabulary) is interesting. In \circled{8},
the performance of the network can be recovered close to the
baseline by retraining either parent source embeddings or the encoder. These
two components can compensate for each other.
This differs from
the case with English reused in the source (\circled{4}) where updating
embeddings to the child language is insufficient: the decoder must be updated to produce
fluent output in the new target language and even with the decoder updated, the
loss compared to the baseline is quite substantial.

The most important component for transfer learning is generally the
component handling the new language:
decoder in English-to-Estonian
and encoder in
the reverse.
With this component fixed,
the performance drops the most with this component fixed (\circled{1},
\circled{2}, \circled{5}, \circled{6}) and among the least with this component free to
update (\circled{3}, \circled{4}, \circled{7}, \circled{8}).
This confirms that at least for examined language pair, the Transformer model
lends itself very well to encoder or decoder re-use.

Other results in \cref{tab:freezing_analysis} reveal that the architecture can compensate for some of the training deficiencies. Freezing the encoder \circled{1}, \circled{2} (resp. decoder for Estonian-to-English \circled{5}, \circled{6}) or attention is not that critical as the frozen decoder (resp. encoder). The bad result of the encoder \circled{3}, \circled{4} (resp. decoder \circled{7}, \circled{8}) being the only non-frozen component shows that model is not capable of providing all the needed capacity for the new language, unlike the self-attention where the loss is not that large. This behaviour correlates with our intuition that the model needs to update the most the component that handles the differing language with the parent model (in our case Czech).

All in all, these experiments illustrate the robustness of the Transformer model that it is able to train and reasonably well utilize pre-trained weights even if they are severely crippled.

\section{Related Work}
\label{sec:related}

This paper focuses on re-using an existing NMT model in order to improve the performance in terms of training time and translation quality without any need to modify the model or pre-trained weights.

\perscite{dynamic_vocabulary} presented two model modifications for multilingual MT and showed that transfer learning could be extended to transferring from the parent to the first child, followed by the second child and then the third one. They achieved improvements with dynamically updating embeddings for the vocabulary of a target language.

The use of other language pairs for improving results for the target language pair has been approached from various angles. One option is to build multilingual models \parcite{liu2020multilingual}, ideally so that they are capable of zero-shot, i.e. translating in a translation direction that was never part of the training data.
\perscite{TACL1081} and \perscite{lu-EtAl:2018:WMT1} achieve this with a unique language tag that specifies the desired target language.
The training data includes sentence pairs from multiple language pairs, and the model implicitly learns translation among many languages. In some cases, it achieves zero-shot and can translate between languages never seen together. \perscite{gu-EtAl:2018:N18-1} tackled the problem by creating universal embedding space across multiple languages and training many-to-one MT system. \perscite{firat-cho-bengio:etal:2016} propose multi-way multi-lingual systems. Their goal is to reduce the total number of parameters needed to train multiple source and target models. In all cases, the methods are dependent on a special training schedule.

The lack of parallel data in low-resource language pairs can also be tackled by unsupervised translation \parcite{artetxe2017unsupervised,lample2018phrase}. The general idea is to train monolingual autoencoders for both source and target languages separately, followed by mapping both embeddings to the same space and training simultaneously two models, each translating in a different direction. In an iterative training, this pair of NMT systems is further refined, each system providing training data for the other one by back-translating monolingual data \parcite{sennrich-haddow-birch:2016:monolingual}.

For very closely related language pairs, transliteration can be used to generate training data from a high-resourced pair to support the low-resourced one as described in \perscite{Karakanta2018}.

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

In this paper, we focus on a setting where existing models are re-used without any preparation for knowledge transfer of original model ahead of its training. This is a relevant and prevailing situation in academia due to computing restrictions, and industry, where updating existing models and scaling to more language pairs is essential.
We evaluate and propose methods of re-using Transformer NMT models for any ``child'' language pair regardless of the original ``parent'' training languages and especially showing, that no modification is better than training from scratch.

The techniques are simple, effective, and applicable to models trained by others which makes it more likely that our experimental results will be replicated in practice.
We showed that despite the random assignment of subwords, the Transformed Vocabulary improves the performance and shortens the training time of the child model compared to training from random initialization.

Furthermore, we showed that this approach is not restricted to low-resource languages, and we documented that the highest improvements are (expectably) due to the shared English knowledge.
Moreover, we confirmed the robustness of the Transformer and its ability to achieve good results in adverse conditions like very fragmented subword units or parts of the network frozen.

The warm-start approach by \perscite{kocmi-adaption} performs slightly better than our Transformed Vocabulary, but it needs to be trained for a significantly longer time.
This leaves room for approaches that also focus on the efficiency of the training process. We perceive our approach as a technique for increasing the performance of a model without an increase in training time. Thus, re-using older models in cold-start scenario of transfer learning can be used in standard NMT training pipelines without any performance or speed losses instead of random initialization as is the common practice currently.

\section*{Acknowledgements}

This study was supported in parts by the grants 18-24210S of the Czech Science Foundation and 825303 (Bergamot) of the European Union.
This work has been using language resources and tools stored and distributed by the LINDAT/CLARIN project of the Ministry of Education, Youth and Sports of the Czech Republic (LM2015071).