\section{Introduction}

Unnormalized probability models have drawn significant attention over the years. These models arise, for example, in energy based models, where the normalization constant is intractable to compute, and are thus relevant to numerous settings. Particularly, they have been extensively used in the context of restricted Boltzmann machines \citep{smolensky1986information,hinton2002training}, deep belief networks \citep{hinton2006fast,salakhutdinov2009deep}, Markov random fields \citep{carreira2005contrastive,hinton2006reducing}, and recently also with deep neural networks \citep{xie2016theory,song2019generative,du2019implicit,grathwohl2019your,nijkamp2019learning}.

Fitting an unnormalized density model to a dataset is challenging due to the missing normalization constant of the distribution. A naive approach is to employ approximate maximum likelihood estimation (MLE). This approach relies on the fact that the likelihood's gradient can be approximated using samples from the model, generated using Markov Chain Monte Carlo (MCMC) techniques. However, a good approximation requires using very long chains and is thus impractical. This difficulty motivated the development of a plethora of more practical approaches, like score matching \citep{hyvarinen2005estimation}, noise contrastive estimation (NCE) \citep{gutmann2010noise}, and conditional NCE (CNCE) \citep{ceylan2018conditional}, which replace the log-likelihood loss with objectives that do not require the computation of the normalization constant or its gradient.

Perhaps the most popular method for learning unnormalized models is contrastive divergence (CD) \citep{hinton2002training}. CD's advantage over MLE stems from its use of short Markov chains initialized at the data samples. CD has been successfully used in a wide range of domains, including modeling images \citep{hinton2006fast}, speech \citep{mohamed2010phone}, documents \citep{hinton2009replicated}, and movie ratings \citep{salakhutdinov2007restricted}, and is continuing to attract significant research attention \citep{liu2017learning,gao2018learning,qiu2019unbiased}.

Despite CD's popularity and empirical success, there still remain open questions regarding its theoretical properties. The primary source of difficulty is an unjustified approximation used to derive its objective's gradient, which biases its update steps \citep{carreira2005contrastive,bengio2009justifying}. The difficulty is exacerbated by the fact that CD's update steps cannot be expressed as the gradients of any fixed objective \citep{tieleman2007some,sutskever2010convergence}.

In this paper, we present an alternative derivation of CD, which relies on completely different principles and requires no approximations. Specifically, we show that CD's update steps are the gradients of an adversarial game in which a discriminator attempts to classify whether a Markov chain generated from the model is presented to it in its original or a time-reversed order (see Fig.~\ref{fig:adverserial_cd}). Thus, our derivation sheds new light on CD's success: Similarly to modern generative adversarial methods \citep{goodfellow2014generative}, CD's discrimination task becomes more challenging as the model approaches the true distribution. This keeps the update steps effective throughout the entire training process and prevents early saturation as often happens in non-adaptive methods like NCE and CNCE. In fact, we derive CD as a natural extension of the CNCE method, replacing the fixed distribution of the contrastive examples with an adversarial adaptive distribution.

\begin{figure}[t]
    \centering
    \includegraphics[width=\textwidth]{adverserial_cd.pdf}
    \caption{
        Contrastive divergence as an adversarial process. In the first step, the distribution model is used to generate an MCMC process which is used to generate a chain of samples. In the second step the distribution model is updated using a gradient descent step, using the MCMC transition rule.
    }
    \label{fig:adverserial_cd}
\end{figure}

CD requires that the underlying MCMC be exact, which is not the case for popular methods like Langevin dynamics. This commonly requires using Metropolis-Hastings (MH) rejection, which ignores some of the generated samples. Interestingly, our derivation reveals an alternative correction method for inexact chains, which does not require rejection.

\section{Background}

\subsection{The Classical Derivation of CD}

Assume we have an unnormalized distribution model $p_\theta$. Given a dataset of samples $\{x_i\}$ independently drawn from some unknown distribution $p$, CD attempts to determine the parameters $\theta$ with which $p_\theta$ best explains the dataset. Rather than using the log-likelihood loss, CD's objective involves distributions of samples along finite Markov chains initialized at~$\{x_i\}$. When based on chains of length $k$, the algorithm is usually referred to as CD-$k$.

Concretely, let $q_\theta(x'|x)$ denote the transition rule of a Markov chain with stationary distribution $p_\theta$, and let $r^m_\theta$ denote the distribution of samples after $m$ steps of the chain. As the Markov chain is initialized from the dataset distribution and converges to $p_\theta$, we have that $r_\theta^0=p$ and $r_\theta^\infty=p_\theta$. The CD algorithm then attempts to minimize the loss
\begin{align}
\label{eq:cd_loss}
\ell_{\text{CD-}k}
&=D_{\mathrm{KL}}(r^0_\theta||r^{\infty}_\theta)-D_{\mathrm{KL}}(r^k_\theta||r^{\infty}_\theta)\nonumber\\
&=D_{\mathrm{KL}}(p||p_{\theta})-D_{\mathrm{KL}}(r^{k}_\theta||p_{\theta}),
\end{align}
where $D_{\mathrm{KL}}$ is the Kullback-Leibler divergence. Under mild conditions on $q_\theta$ \citep{cover1994processes} this loss is guaranteed to be positive, and it vanishes when $p_{\theta}=p$ (in which case $r^{k}_\theta=p_\theta$).

To allow the minimization of (\ref{eq:cd_loss}) using gradient-based methods, one can write
\begin{align}
\label{eq:cd_grad}
\nabla_{\theta}\ell_{\text{CD-}k}=
&\mathbb{E}_{\tilde{X}\sim r^k_{\theta}} [\nabla_{\theta}\log p_{\theta}(\tilde{X})]
-\mathbb{E}_{X\sim p} [\nabla_{\theta}\log p_{\theta}(X)] +\frac{dD_{\mathrm{KL}}(r^k_{\theta}||p_{\theta})}{d r^k_{\theta}}\nabla_{\theta}r^k_{\theta}.
\end{align}
Here, the first two terms can be approximated using two batches of samples, one drawn from $p$ and one from $r^k_\theta$. The third term is the derivative of the loss with respect only to the $\theta$ that appears in $r_{\theta}^k$, ignoring the dependence of $p_{\theta}$ on $\theta$. This is the original notation from \citep{hinton2002training}; an alternative way to write this term would be $\nabla_{\tilde{\theta}}D_{\mathrm{KL}}(r^{k}_{\tilde{\theta}}||p_{\theta})$. This term turns out to be intractable and in the original derivation, it is argued to be small and thus neglected, leading to the approximation
\begin{equation}
\label{eq:cd_grad_approx}
\nabla_{\theta}\ell_{\text{CD-}k}\approx
\frac{1}{n}\sum_i\left(
\nabla_{\theta}\log p_{\theta}(\tilde{x}_i)
-\nabla_{\theta}\log p_{\theta}(x_i)
\right)
\end{equation}
Here $\{x_i\}$ is a batch of $n$ samples from the dataset and~$\{\tilde{x}_i\}$ are $n$ samples generated by applying $k$ MCMC steps to each of the samples in that batch. The intuition behind the resulting algorithm (summarized in App.~\ref{appendix:algo}) is therefore simple. In each gradient step $\theta\leftarrow\theta-\eta\nabla_{\theta}\ell_{\text{CD-}k}$, the log-likelihood of samples from the dataset is increased on the expense of the log-likelihood of the contrastive samples $\{\tilde{x}_i\}$, which are closer to the current learned distribution~$p_\theta$.

Despite the simple intuition, it has been shown that without the third term, CD's update rule (\ref{eq:cd_grad}) generally cannot be the gradient of any fixed objective \citep{tieleman2007some,sutskever2010convergence} except for some very specific cases. For example, \cite{hyvarinen2007connections} has shown that when the Markov chain is based on Langevin dynamics with a step size that approaches zero, the update rule of CD-1 coincides with that of score-matching \cite{hyvarinen2005estimation}. Similarly, the probability flow method of \cite{sohl2011minimum} has been shown to be equivalent to CD with a very unique Markov chain. Here, we show that regardless of the selection of the Markov chain, the update rule is in fact the exact gradient of a particular adversarial objective, which adapts to the current learned model in each step.

\subsection{Conditional Noise Contrastive Estimation}
\label{sec:cnce}

Our derivation views CD as an extension of the CNCE method, which itself is an extension of NCE. We therefore start by briefly reviewing those two methods.

In NCE, the unsupervised density learning problem is transformed into a supervised one. This is done by training a discriminator $D_\theta(x)$ to distinguish between samples drawn from $p$ and samples drawn from some preselected contrastive distribution $p_\text{ref}$. Specifically, let the random variable $Y$ denote the label of the class from which the variable $X$ has been drawn, so that $X|(Y=1)\sim p$ and $X|(Y=0)\sim p_\text{ref}$. Then it is well known that the discriminator minimizing the binary cross-entropy (BCE) loss is given by
\begin{equation}
D_\text{opt}(x)=\mathbb{P}(Y=1|X=x)=\frac{p(x)}{p(x)+p_{\text{ref}}(x)}.
\end{equation}
Therefore, letting our parametric discriminator have the form
\begin{equation}
D_\theta(x)=\frac{p_{\theta}(x)}{p_{\theta}(x)+p_{\text{ref}}(x)},
\end{equation}
and training it with the BCE loss, should in theory lead to $D_\theta(x)=D_\text{opt}(x)$ and thus to $p_\theta(x)=p(x)$. In practice, however, the convergence of NCE highly depends on the selection of $p_{\text{ref}}$. If it significantly deviates from $p$, then the two distributions can be easily discriminated even when the learned distribution $p_{\theta}$ is still very far from~$p$. At this point, the optimization essentially stops updating the model, which can result in a very inaccurate estimate for $p$. In the next section we provide a precise mathematical explanation for this behavior.

The CNCE method attempts to alleviate this problem by drawing the contrastive samples based on the dataset samples. Specifically, each dataset sample $x$ is paired with a contrastive sample~$\tilde{x}$ that is drawn conditioned on $x$ from some predetermined conditional distribution $q(\tilde{x}|x)$ (\emph{e.g.} $\mathcal{N}(x,\sigma^2I)$). The pair is then concatenated in a random order, and a discriminator is trained to predict the correct order. This is illustrated in Fig.~\ref{fig:cnce}. Specifically, here the two classes are of pairs $(A,B)$ corresponding to $(A,B)=(X,\tilde{X})$ for $Y=1$, and $(A,B)=(\tilde{X},X)$ for $Y=0$, and the discriminator minimizing the BCE loss is given by
\begin{align}
D_\text{opt}(a,b)&=\mathbb{P}(Y=1|A=a,B=b) =\frac{q(b|a)p(a)}{q(b|a)p(a)+q(a|b)p(b)}.
\end{align}
Therefore, constructing a parametric discriminator of the form
\begin{equation}\label{eq:cnce_discriminator}
D_{\theta}(a,b)
=\frac{q(b|a)p_\theta(a)}{q(b|a)p_\theta(a)+q(a|b)p_\theta(b)}
=\left(1+\frac {q(a|b)p_{\theta}(b)}{q(b|a)p_{\theta}(a)}\right)^{-1},
\end{equation}
and training it with the BCE loss, should lead to $p_\theta\propto p$. Note that here $D_{\theta}$ is indifferent to a scaling of $p_{\theta}$, which is thus determined only up to an arbitrary multiplicative constant.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.5\textwidth}\centering
        \includegraphics[width=0.95\textwidth]{cnce.pdf}
        \caption{CNCE}\label{fig:cnce}
    \end{subfigure
    \begin{subfigure}[b]{0.5\textwidth}\centering
        \includegraphics[width=0.95\textwidth]{cnce_mcmc.pdf}
        \caption{CD-$1$}\label{fig:cnce_mcmc}
    \end{subfigure}
    \caption{\textbf{From CNCE to CD-$1$.}
        (a) In CNCE, each contrastive sample is generated using a fixed conditional distribution $q(\cdot|\cdot)$ (which usually corresponds to additive noise). The real and fake samples are then concatenated and presented to a discriminator in a random order, which is trained to predict the correct order. (b) CD-$1$ can be viewed as CNCE with a $q(\cdot|\cdot)$ that corresponds to the transition rule of a Markov chain with stationary distribution~$p_\theta$. Since~$q$ depends on $p_\theta$ (hence the subscript $\theta$), during training the distribution of contrastive samples becomes more similar to that of the real samples, making the discrimination task harder.
    }
    \label{fig:cnce_CD1}
\end{figure}

CNCE improves upon NCE, as it allows working with contrastive samples whose distribution is closer to $p$. However, it does not completely eliminate the problem, especially when $p$ exhibits different scales of variation in different directions. This is the case, for example, with natural images, which are known to lie close to a low-dimensional manifold. Indeed if the conditional distribution $q(\cdot|\cdot)$ is chosen to have a small variance, then CNCE fails to capture the global structure of $p$. And if $q(\cdot|\cdot)$ is taken to have a large variance, then CNCE fails to capture the intricate features of $p$ (see Fig.~\ref{fig:toy_experiment}). The latter case can be easily understood in the context of images (see Fig.~\ref{fig:cnce}). Here, the discriminator can easily distinguish which of its pair of input images is the noisy one, without having learned an accurate model for the distribution of natural images (\emph{e.g.,}  simply by comparing their smoothness). When this point is reached, the optimization essentially stops.

In the next section we show that CD is in fact an adaptive version of CNCE, in which the contrastive distribution is constantly updated in order to keep the discrimination task hard. This explains why CD is less prone to early saturation than NCE and CNCE.

\begin{figure}[tb!]
    \centering
    \begin{subfigure}[t]{0.19\textwidth}
        \centering
        \includegraphics[scale=0.5, trim=0cm -2.8inch 0cm 0cm]{media/spiral_3d.pdf}
        \caption{The toy model}
        \label{fig:toy_model}
    \end{subfigure}

    \begin{subfigure}[t]{0.76\textwidth}
        \centering
        \includegraphics[scale=0.5]{media/spiral_cnce_vs_acd.pdf}
        \caption{Comparing CNCE with CD}
        \label{fig:cnce_vs_acd}
    \end{subfigure}
    \caption{
        A toy example illustrating the importance of the adversarial nature of CD. Here, the data lies close to a 2D spiral embedded in a 10-dimensional space. (a)~The training samples in the first~3 dimensions. (b)~Three different approaches for learning the distribution: CNCE with large contrastive variance (top), CNCE with small contrastive variance (middle), and CD based on Langevin dynamics MCMC with the weight adjustment described in Sec.~\ref{sec:ACD} (bottom). As can be seen in the first two columns, CD adapts the contrastive samples according to the data distribution, whereas CNCE does not. Therefore, CNCE with large variance fails to learn the distribution because the vast majority of its contrastive samples are far from the manifold and quickly become irrelevant (as indicated by the weights $\alpha_\theta$ in the third column). And CNCE with small variance fails to learn the global structure of the distribution because its contrastive samples are extremely close to the dataset samples. CD, on the other hand, adjusts the contrastive distribution during training, so as to generate samples that are close to the manifold yet traverse large distances along it.
    }\label{fig:toy_experiment}
\end{figure}

\section{An Alternative Derivation of CD}

We now present our alternative derivation of CD. In Sec.~\ref{sec:decom} we identify a decomposition of the CNCE loss, which reveals the term that is responsible for early saturation. In Sec.~\ref{sec:CD1}, we then present a method for adapting the contrastive distribution in a way that provably keeps this term bounded away from zero. Surprisingly, the resulting update step turns out to precisely match that of CD-$1$, thus providing a new perspective on CD learning. In Sec.~\ref{sec:CDk}, we extend our derivation to include CD-$k$ (with $k\geq1$).

\subsection{Reinterpreting CNCE}
\label{sec:decom}

Let us denote
\begin{equation}\label{eq:W1}
w_{\theta}(a,b)\triangleq\frac{q(a|b)p_{\theta}(b)}{q(b|a)p_{\theta}(a)},
\end{equation}
so that we can write CNCE's discriminator (\ref{eq:cnce_discriminator}) as
\begin{equation}\label{eq:CNCEdiscrimnator}
D_{\theta}(a,b)=\left(1+w_{\theta}(a,b)\right)^{-1}.
\end{equation}
Then we have the following observation (see proof in App.~\ref{appendix:cnce_grad}).

\begin{observation}
\label{obs:CNCEgrad}
The gradient of the CNCE loss can be expressed as
\begin{align}
\label{eq:cnce_grad}
\nabla_{\theta}\ell_{\text{CNCE}}=\mathbb{E}\!_{\substack{X\sim p \\ \tilde{X}|X\sim q}}\!\left[\alpha_\theta(X,\tilde{X})\left(\nabla_{\theta}\log p_{\theta}(\tilde{X})-\nabla_{\theta}\log p_{\theta}(X)\right)\right],
\end{align}
where
\begin{equation}\label{eq:weight}
\alpha_\theta(x,\tilde{x})\triangleq (1+w_\theta(x,\tilde{x})^{-1})^{-1}.
\end{equation}
\end{observation}

Note that (\ref{eq:cnce_grad}) is similar in nature to the (approximate) gradient of the CD loss~(\ref{eq:cd_grad_approx}). Particularly, as in CD, the term $\nabla_{\theta}\log p_{\theta}(\tilde{X})-\nabla_{\theta}\log p_{\theta}(X)$ causes each gradient step to increase  the log-likelihood of samples from the dataset on the expense of the log-likelihood of the contrastive samples. However, as opposed to CD, here we also have the coefficient $\alpha_{\theta}(x,\tilde{x})$, which assigns a weight between $0$ and $1$ to each pair of samples $(x,\tilde{x})$. To understand its effect, observe that
\begin{equation}
\alpha_\theta(x,\tilde{x})=1-D_\theta(x,\tilde{x})=D_{\theta}(\tilde{x},x).
\end{equation}
Namely, this coefficient is precisely the probability that the discriminator assigns to the incorrect order of the pair. Therefore, this term gives a low weight to ``easy'' pairs (\emph{i.e.,}  for which $D_\theta(x,\tilde{x})$ is close to $1$) and a high weight to ``hard'' ones.

This weighting coefficient is of course essential for ensuring convergence to $p$. For example, it prevents  $\log p_{\theta}$ from diverging to $\pm\infty$ when the discriminator is presented with the same samples over and over again. The problem is that a discriminator can often correctly discriminate all training pairs, even with a $p_\theta$ that is still far from $p$. In such cases, $\alpha_\theta$ becomes practically zero for all pairs and the model stops updating. This shows that a good contrastive distribution is one which keeps the discrimination task hard throughout the training. As we show next, there is a particular choice which provably prevents $\alpha_\theta$ from converging to zero, and that choice results in the CD method.

\subsection{From CNCE to CD-\(1\)}
\label{sec:CD1}

To bound $\alpha_\theta$ away from $0$, and thus avoid the early stopping of the training process, we now extend the original CNCE algorithm by allowing the conditional distribution $q$ to depend on $p_\theta$ (and thus to change from one step to the next). Our next key observation is that in this setting there exists a particular choice that keeps $\alpha_\theta$ constant.

\begin{observation}
If $q$ is chosen to be the transition probability of a reversible Markov chain with stationary distribution $p_\theta$, then
\begin{equation}
\alpha_{\theta}(x,\tilde{x})=\frac{1}{2}, \quad\forall x,\tilde{x}.
\end{equation}
\end{observation}
\begin{proof}
A reversible chain with transition $q$ and stationary distribution $p_\theta$, satisfies the detailed balance property
\begin{equation}\label{eq:detailed_balance}
    q(\tilde{x}|x)p_{\theta}(x)=q(x|\tilde{x})p_{\theta}(\tilde{x}),\quad \forall x,\tilde{x}.
\end{equation}
Substituting (\ref{eq:detailed_balance}) into (\ref{eq:W1}) leads to $w_{\theta}(x,\tilde{x})=1$, which from (\ref{eq:weight}) implies $\alpha_\theta(x,\tilde{x})=\tfrac{1}{2}$.
\end{proof}

This observation directly links CNCE to CD. First, the suggested method for generating the contrastive samples is precisely the one used in CD-$1$. Second, as this choice of $q$ leads to $\alpha_\theta(x,\tilde{x})=\tfrac{1}{2}$, it causes the gradient of the CNCE loss (\ref{eq:cnce_grad}) to become
\begin{equation}
\label{eq:cnce_mcmc_grad_exact}
\nabla_{\theta}\ell_{\text{CNCE}}=\tfrac{1}{2}\mathbb{E}_{\substack{X\sim p \\ \tilde{X}|X\sim q}}
\left[\nabla_{\theta}\log p_{\theta}(\tilde{X})-\nabla_{\theta}\log p_{\theta}(X)\right],
\end{equation}
which is exactly proportional to the CD-$1$ update (\ref{eq:cd_grad_approx}). We have thus obtained an alternative derivation of CD-$1$. Namely, rather than viewing CD-$1$ learning as an \emph{approximate} gradient descent process for the loss (\ref{eq:cd_loss}), we can view each step as the \emph{exact} gradient of the CNCE discrimination loss, where the reference distribution $q$ is adapted to the current learned model $p_\theta$. This is illustrated in Fig.~\ref{fig:cnce_mcmc}.

Since $q$ is chosen based on $p_\theta$, the overall process is in fact an adversarial game. Namely, the optimization alternates between updating $q$, which acts as a generator, and updating $p_{\theta}$, which defines the discriminator. As $p_{\theta}$ approaches $p$, the distribution of samples generated from the MCMC also becomes closer to $p$, which makes the discriminator's task harder and thus prevents early saturation.

It should be noted that formally, since $q$ depends on $p_\theta$, it also indirectly depends on $\theta$, so that a more appropriate notation would be $q_\theta$. However, during the update of $p_{\theta}$ we fix $q_{\theta}$ (and vise versa), so that the gradient in the discriminator update does not consider the dependence of $q_\theta$ on $\theta$. This is why (\ref{eq:cnce_mcmc_grad_exact}) does not involve the gradient of $\tilde{X}$ which depends on $q_{\theta}$.

The reason for fixing $q_{\theta}$ comes from the adversarial nature of the learning process. Being part of the chain generation process, the goal of the transition rule $q_{\theta}$ is to generate chains that appear to be time-reversible, while the goal of the classifier, which is based on the model $p_{\theta}$, is to correctly classify whether the chains were reversed. Therefore, we do not want the optimization of the classifier to affect~$q_{\theta}$. This is just like in GANs, where the generator and discriminator have different objectives, and so when updating the discriminator the generator is kept fixed.

\begin{comment}
\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{cnce_mcmc.pdf}
    \caption{
        CD-$1$ viewed as CNCE with a single MCMC step. Each contrastive sample is generated by taking a single step in an MCMC process starting with a real sample from dataset. The two samples are then concatenated in a random order, and a discriminator is trained to predict the correct order. As the training advances, the distribution of contrastive samples becomes more similar to that of real samples, and thus the discrimination task becomes harder.
    }
    \label{fig:cnce_mcmc}
\end{figure}
\end{comment}

\subsection{From CD-\(1\) to CD-\(k\)}
\label{sec:CDk}

To extend our derivation to CD-$k$ with an arbitrary $k\geq1$, let us now view the discrimination problem of the previous section as a special case of a more general setting. Specifically, the pairs of samples presented to the discriminator in Sec.~\ref{sec:CD1}, can be viewed as Markov chains of length two (comprising the initial sample from the dataset and one extra generated sample). It is therefore natural to consider also Markov chains of arbitrary lengths. That is, assume we initialize the MCMC at a sample $x_i$ from the dataset and run it for $k$ steps to obtain a sequence $(x^{(0)},x^{(1)},\ldots,x^{(k)})$, where  $x^{(0)}=x_i$. We can then present this sequence to a discriminator either in its original order, or time-reversed, and train the discriminator to classify the correct order. We coin this a \emph{time-reversal classification task}. Interestingly, in this setting, we have the following.

\begin{observation}\label{obv:cdk}
When using a reversible Markov chain of length $k+1$ with stationary distribution $p_\theta$, the gradient of the BCE loss of the time-reversal classification task is given by
\begin{equation}\label{eq:CDkgrad}
\nabla_{\theta}\ell_{\text{CNCE}}=\tfrac{1}{2}\mathbb{E}\left[\nabla_{\theta}\log p_{\theta}(X^{(k)})-\nabla_{\theta}\log p_{\theta}(X^{(0)})\right],
\end{equation}
which is exactly identical to the CD-$k$ update (\ref{eq:cd_grad_approx}) up to a multiplicative factor of $\tfrac{1}{2}$.
\end{observation}

This constitutes an alternative interpretation of CD-$k$. That is, CD-$k$ can be viewed as a time-reversal adversarial game, where in each step, the model $p_\theta$ is updated so as to allow the discriminator to better distinguish MCMC chains from their time-reversed counterparts.

Two remarks are in order. First, it is interesting to note that although the discriminator's task is to classify the order of the whole chain, its optimal strategy is to examine only the endpoints of the chain, $x^{(0)}$ and $x^{(k)}$. Second, it is insightful to recall that the original motivation behind the CD-$k$ loss (\ref{eq:cd_loss}) was that when $p_\theta$ equals~$p$, the marginal probability of each individual step in the chain is also $p$. Our derivation, however, requires more than that. To make the chain indistinguishable from its time-reversed version, the joint probability of all samples in the chain must be invariant to a flip of the order. When $p_\theta=p$, this is indeed the case, due to the detailed balance property (\ref{eq:detailed_balance}).

\begin{proof}[Proof of Observation~\ref{obv:cdk}]
We provide the outline of the proof (see full derivation in App. \ref{appendix:cnce_chain_grad}). Let $(A^{(0)},A^{(1)},\ldots,A^{(k)})$ denote the input to the discriminator and let $Y$ indicate the order of the chain, with $Y=1$ corresponding to $(A^{(0)},A^{(1)},\ldots,A^{(k)})=(X^{(0)},X^{(1)},\ldots,X^{(k)})$ and $Y=0$ to $(A^{(0)},A^{(1)},\ldots,A^{(k)})=(X^{(k)},X^{(k-1)},\ldots,X^{(0)})$. The discriminator that minimizes the BCE loss is now given by
\begin{align}\label{eq:discriminator_long_chain}
D(a_0,a_1,\ldots,a_k)&=\mathbb{P}(Y=1|A^{(0)}=a_0,A^{(1)}=a_1,\ldots,A^{(k)}=a_k) \nonumber\\
&=\left(1+\frac
{q(a_0|a_1)\cdots q(a_{k-1}|a_k)p(a_k)}
{q(a_k|a_{k-1})\cdots q(a_1|a_0)p(a_0)}
\right)^{-1} \nonumber\\
&=\left(1+\prod_{i=1}^k w_{\theta}(a_{i-1},a_i)\right)^{-1}.
\end{align}
The CNCE paradigm thus defines a discriminator $D_\theta$ having the form of (\ref{eq:discriminator_long_chain}) but with $p$ replaced by~$p_\theta$. Recall that despite the dependence of the transition probability~$q$ on the current learned model~$p_\theta$, it is regarded as fixed within each discriminator update step. We therefore omit the subscript $\theta$ from $q$ here. Similarly to the derivation of (\ref{eq:cnce_grad}), explicitly writing the gradient of the BCE loss of our discriminatrion task, gives
\begin{align}
\label{eq:cnce_chain_grad}
\nabla_{\theta}\ell_{\text{chain}}
&=\mathbb{E}
\left[\left(1+\prod_{i=1}^k w_{\theta}(X^{(i-1)},X^{(i)})^{-1}\right)^{-1}\left(\nabla_{\theta}\log p_{\theta}(X^{(k)})-\nabla_{\theta}\log p_{\theta}(X^{(0)})\right)
\right]\nonumber \\
&=\mathbb{E}
\left[\alpha_\theta(X^{(0)},\ldots,X^{(k)})\left(\nabla_{\theta}\log p_{\theta}(X^{(k)})-\nabla_{\theta}\log p_{\theta}(X^{(0)})\right)
\right].
\end{align}
where we now defined
\begin{equation}\label{eq:weight_chain}
\alpha_\theta(a_0,\ldots,a_k)\triangleq \left(1+\prod_{i=1}^k w_{\theta}(a_{i-1},a_i)^{-1}\right)^{-1}.
\end{equation}
Note that (\ref{eq:cnce_grad}) is a special case of (\ref{eq:cnce_chain_grad}) corresponding to $k=1$, where $X$ and $\tilde{X}$ in (\ref{eq:cnce_grad}) are $X^{(0)}$ and $X^{(1)}$ in (\ref{eq:cnce_chain_grad}). As before, when $q$ satisfies the detailed balance property (\ref{eq:detailed_balance}), we obtain $w_{\theta}=1$ and consequently the weighting term $\alpha_\theta$ again equals $\tfrac{1}{2}$. Thus, the gradient (\ref{eq:cnce_chain_grad}) reduces to (\ref{eq:CDkgrad}), which is exactly proportional to the CD-$k$ update (\ref{eq:cd_grad_approx}).
\end{proof}

\subsection{MCMC Processes That do not Have Detailed Balance}
\label{sec:ACD}

In our derivation, we assumed that the MCMC process is reversible, and thus exactly satisfies the detailed balance property (\ref{eq:detailed_balance}). This assumption ensured that $w_{\theta}=1$ and thus $\alpha_\theta=\tfrac{1}{2}$. In practice, however, commonly used MCMC methods satisfy this property only approximately. For example, the popular discrete Langevin dynamics process obeys detailed balance only in the limit where the step size approaches zero. The common approach to overcome this is through Metropolis-Hastings (MH) rejection \citep{hastings1970monte}, which guarantees detailed balance by accepting only a portion of the proposed MCMC transitions. In this approach, the probability of accepting a transition from $x$ to $\tilde{x}$ is closely related to the weighing term $w_{\theta}$, and is given by
\begin{equation}
A(x,\tilde{x})=\min\left(1,w_{\theta}(x,\tilde{x})\right).
\end{equation}
Interestingly, our derivation reveals an alternative method for accounting for lack of detailed balance.

Concretely, we saw that the general expression  for the gradient of the BCE loss (before assuming detailed balance) is given by (\ref{eq:cnce_chain_grad}). This expression differs from the original update step of CD-$k$ only in the weighting term $\alpha_\theta(x^{(0)},\ldots,x^{(k)})$. Therefore, all that is required for maintaining correctness in the absence of detailed balance, is to weigh each chain by its ``hardness''  $\alpha_\theta(x^{(0)},\ldots,x^{(k)})$ (see Alg.~\ref{algo:acd} in App.~\ref{appendix:algo}). Note that in this case, the update depends not only on the end-points of the chains, but rather also on their intermediate steps. As can be seen in Fig. \ref{fig:cd_flavors}, this method performs just as well as MH, and significantly better than vanilla CD without correction.

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.5]{media/spiral_cd_variants_2d.pdf}
    \caption{
        Here, we use different CD configurations for learning the model of Fig.~\ref{fig:toy_experiment}. All configurations use Langevin dynamics as their MCMC process, but with different ways of compensating for the lack of detailed balance. From left to right we have the ground-truth density, CD w/o any correction, CD with Metropolis-Hastings rejection, and CD with our proposed adjustment.
        }
    \label{fig:cd_flavors}
\end{figure}

\section{Illustration Through a Toy Example}
To illustrate our observations, we now conclude with a simple toy example (see Fig.~\ref{fig:toy_experiment}). Our goal here is not to draw general conclusions regarding the performance of CNCE and CD, but rather merely to highlight the adversarial nature of CD and its importance when the data density exhibits different scales of variation along different directions.

We take data concentrated around a 2-dimensional manifold embedded in 10-dimensional space. Specifically, let $e^{(1)},\dots,e^{(10)}$ denote the standard basis in $\mathbb{R}^{10}$. Then each data sample is generated by adding Gaussian noise to a random point along a 2D spiral lying in the $e^{(1)}$-$e^{(2)}$ plane. The STD of the noise in the $e^{(1)}$ and $e^{(2)}$ directions is 5 times larger than that in the other 8 axes. Figure \ref{fig:toy_model} shows the projections of the the data samples onto the the first 3 dimensions. Here, we use a multi-layer perceptron (MLP) as our parametric model, $\log p_\theta$, and train it using several different learning configurations (for the full details see App.~\ref{appendix:experiment_details}).

Figure~\ref{fig:cnce_vs_acd} visualizes the training as well as the final result achieved by each configuration. The first two rows show CNCE with Gaussian contrastive distributions of two different STDs. The third row shows the adjusted CD described in Sec.~\ref{sec:ACD} with Langevin Dynamics as its MCMC process. As can be seen, for CNCE with a large STD, the contrastive samples are able to explore large areas around the original samples, but this causes their majority to lie relatively far from the manifold (see their projections onto the $e^{(1)}$-$e^{(3)}$ plane). In this case, $\alpha_\theta$ decreases quickly, causing the learning process to ignore most samples at a very early stage of the training. When using CNCE with a small STD, the samples remain relevant throughout the training, but this comes at the price of inability to capture the global structure of the distribution. CD, on the other hand, is able to enjoy the best of both worlds as it adapts the contrastive distribution over time. Indeed, as the learning progresses, the contrastive samples move closer to the manifold to maintain their relevance. Note that since we use the adjusted version of CD, the weights in this configuration are not precisely $1$. We chose the step size of the Langevin Dynamics so that the median of the weights is approximately $10^{-2}$.

Figure~\ref{fig:cd_flavors} shows the results achieved by different variants of CD. As can be seen, without correcting for the lack of detailed balance, CD fails to estimate the density correctly. When using MH rejection to correct the MCMC, or our adaptive CD (ADC) to correct the update steps, the estimate is significantly improved.

\section{Conclusion}
The classical CD method has seen many uses and theoretical analyses over the years. The original derivation presented the algorithm as an approximate gradient descent process for a certain loss. However, the accuracy of the approximation has been a matter of much dispute, leaving it unclear what objective the algorithm minimizes in practice. Here, we presented an alternative derivation of CD's update steps, which involves no approximations. Our analysis shows that CD is in essence an adversarial learning procedure, where a discriminator is trained to distinguish whether a Markov chain generated from the learned model has been time-flipped or not. Therefore, although predating GANs by more than a decade, CD in fact belongs to the same family of techniques. This provides a possible explanation for its empirical success.

\paragraph{Acknowledgement} This research was supported by the Technion
Ollendorff Minerva Center.


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

Let us consider an unlabeled training set ${\cal T} = \{{\bf x}_1, \ldots, {\bf x}_n\}$ of $n$ independent random vectors (i.e., patterns) in a $d$-dimensional feature space, say $\mathbb{R}^d$. The patterns are assumed to be identically distributed according to an unknown probability density function (pdf) $p({\bf x})$. Density estimation consists in finding a model for $p({\bf x})$ based on ${\cal T}$. This requires exploiting the statistical knowledge on $p({\bf x})$ implicitly underlying the very data sample ${\cal T}$. To this aim, a suitable algorithm is applied, either rooted in statistics or in machine learning. The algorithm is expected to come up with a model of $p({\bf x})$ that fits the data well, in compliance with some well-defined statistical criterion  (the maximum likelihood criterion \cite{DudHar73} being possibly the most popular). The present paper proposes and investigates mixtures of multivariate component densities realized via deep neural networks (DNNs) for density estimation. It extends to multivariate mixture densities our previous work on univariate pdfs, presented as a workshop communication in \cite{Trentin-NMM-ANNPR18}. The reader is referred to that paper for a representative list of up-to-date applications where (implicitly or explicitly) density estimation is fundamental and still an open issue.

\subsection{Intrinsic difficulties of density estimation via machine learning}

As pointed out by Vapnik \cite{Vapnik95}, density estimation is an intrinsically difficult problem, and it is still open nowadays. This latter fact is mostly due to the shortcomings of established statistical approaches, either parametric or non-parametric (the reader is referred to \cite{TrentinLC18} for a list of the major drawbacks of the statistical techniques), and by the technical difficulties that arise from attempting to use artificial neural networks (ANNs) or machine learning for pdf estimation. Such difficulties stem from: (1) the unsupervised nature of the learning task, (2) the numerical instability problems entailed by pdfs, whose codomains may span the interval  $[0, +\infty)$, and (3) the requirement for the resulting model to respect the axioms of probability. Furthermore, the use of maximum-likelihood (ML) training in ANNs tends to result in the ``divergence problem'', observed first in the realm of hybrid ANN/hidden Markov models \cite{Trentin15_ML}. It consists in the progressive divergence of the value of the ANN connection weights as ML training proceeds, resulting in an unbounded growth of the integral of the pseudo-pdf computed by the ANN. The problem does not affect radial basis functions (RBF) networks whose hidden-to-output weights were constrained to be positive and to sum to one, as in the RBF/echo state machine for sequences proposed in \cite{TrentinSS15_ESN}, or in the RBF/graph neural network presented in \cite{BonRigTre2018} for the estimation of generalized random graphs. Unfortunately, the use of RBFs in the latter contexts is justified by its allowing for a proper algorithmic hybridization with models devised specifically for sequence/structure processing, but using RBFs as a stand-alone paradigm for density estimation is of neglectable practical interest, since they end up realizing plain Gaussian mixture models (GMM) estimated via ML.

\subsection{Earlier approaches to ANN-based density estimation}

In spite of these difficulties, several approaches to pdf estimation via ANNs are found in the literature \cite{TrentinFreno09}. A critical survey of the literature is found in \cite{TrentinLC18}. Most of these approaches suffer from some limitations. First, a ML technique is presented in \cite{ModhaFainman94} where the ``integral equals 1'' requirement is satisfied numerically by dividing the output of a multilayer Perceptron (MLP)  by the numerical integral of the function the MLP computes. No algorithms for computing the numerical integral over high-dimensional spaces are handed out in \cite{ModhaFainman94}. Nonetheless, this approach is related to the technique presented in this paper, insofar that ML will be exploited herein. Differently from \cite{ModhaFainman94}, a multi-dimensional ad-hoc numeric integration method will be used in the following, jointly with hard constraints, over a mixture of DNNs.

Instead of estimating the pdf directly, techniques have been proposed that focus on the (theoretically equivalent) estimation of the corresponding cumulative distribution function (cdf) \cite{Vapnik00supportvector,IsmailAtaiya2002den}. Supervised training is used , e.g. the backpropagation algorithm (BP), relying on the empirical cdf of the data for generating the target outputs required in supervised training.  After training the MLP model $\phi(\cdot)$ of the cdf, the pdf can be retrieved by taking derivatives of $\phi(\cdot)$. Unfortunately, there are drawbacks to these approaches, as well (see \cite{TrentinLC18}). In particular, a good approximation of the cdf does not entail a good estimate of its derivative. Negative values of $\frac{\partial \phi(x)}{\partial x}$ may even occur, since a linear combination of logistics is not necessarily monotonically increasing. Finally, these algorithms apply naturally to univariate cases, but their application to  multivariate scenarios is way trickier (to say the least).

The idea of exploiting synthetically-generated target outputs was applied to DNN-based pdf estimation in \cite{TrentinLC18} and \cite{Trentin_Soft_NEPL}.  The former relies on a unbiased variant of the Parzen-window (PW) density estimation technique \cite{DudHar73} for labeling the training set for a DNN, known as Parzen neural network. Like in the traditional $k_n$-Nearest Neighbor ($k_n$-NN) statistical technique \cite{DudHar73}, the resulting model does not satisfy the axioms of probability. To the contrary, the algorithm for density estimation via feed-forward DNNs presented in \cite{Trentin_Soft_NEPL} uses a modified loss function, optimized via gradient descent. Such a criterion function combines two terms: the loss between the network output and a synthetically-generated non-parametric estimate of the pdf at hand evaluated over the corresponding input pattern, and a loss between the integral of the function realized by the MLP and the unity value (that is the training ``target'' for the integral of the DNN). Efficient integration methods based on Markov chain Monte Carlo with importance sampling are used to compute the integral of the DNN and its derivatives w.r.t. the DNN parameters within the gradient-descent scheme used for the optimization. The asymptotic convergence of the resulting technique to the correct solution was formally proven in \cite{TrentinMathematics2020}. The ideas behind such integration methods are exploited in this paper, as well.

\subsection{Mixture densities}

A generalization of plain pdf models stems from the adoption of mixture densities, where the unknown pdf is rather modeled in terms of a combination of any number of component densities \cite{DudHar73}. GMMs are the most popular instance of mixture densities \cite{Bishop:2006:PRM:1162264}. Mixture densities were originally intended as real-life extensions of the single-pdf parametric model (e.g., one Gaussian may not be capable to explain the data distribution but $K$ Gaussian pdfs might as well be). Yet, there is much more than this behind the notion of mixture density. In fact, in mixture densities the different component pdfs are specialized to explain distinct latent phenomena (e.g., stochastic processes) that underlie the overall data generation process, each such phenomenon having different likelihood of occurrence w.r.t. others at diverse regions of the feature space. This suites particularly those situations where the statistical population under analysis is composed of several sub-populations, each having different distribution. Relevant examples are the following:

\begin{enumerate}

\item the {\em statistical study of heterogeneity} in meta-analysis \cite{metanalysis}, where samples drawn from disjoint populations (e.g., adults and children, male and female subjects, etc.) are collectively collected and have to be analyzed as a whole;

\item the modeling of {\em unsupervised} or {\em partially-supervised} \cite{SchwenkerT14a_Survey} data samples in statistical pattern recognition \cite{DudHar73}, where each sub-population corresponds to a class or category;

\item the {\em distribution of financial returns} on the stock market depending on latent phenomena such as a political crisis or a war \cite{finance-crisis};

\item the {\em assessment of projectile accuracy} in the military science of ballistics when shots at the same target come from multiple locations and/or from different munition types \cite{CIS-106239}.

\end{enumerate}

As it happens, the sub-populations in a mixture are unlikely to be distributed individually according to simple (e.g., Gaussian) pdfs. Consequently, parametric models  (like GMMs) are generally not a very good fit. In fact, let $\xi_1, \ldots, \xi_K$ be $K$ disjoint states of nature (the outcomes of a discrete, latent random variable $\Xi$, each outcome corresponding to a specific sub-population), and let $p({\bf x}|\xi_i)$ be the pdf that explains the distribution of the random observations ${\bf x}$ given the $i$-th state of the latent variable, for $i=1, \ldots, K$. At the whole population level the data will then be distributed according to the mixture $p({\bf x}) = \sum_{i=1}^K P(\xi_i)p({\bf x}|\xi_i)$. Attempts to apply a GMM to model $p({\bf x})$ will not necessarily result in a one-to-one relationship between the Gaussian components in the GMM and the state-specific generative models $p({\bf x}|\xi_i)$. In general, at the very least, more than one Gaussian component will be needed to model $p({\bf x}|\xi_i)$. Although mixtures of mixture models offer increased modeling capabilities over plain mixture models to this end, they turned out to be unpopular due to the difficulties of estimation of their parameters \cite{ASteBFT15} and their excessive sensitivity to local maxima of their criterion function (generally the likelihood of their parameters given the data).

\subsection{The proposed DNN-based model}

Aiming at overcoming the aforementioned difficulties encountered with the established approaches, and due to the unexploited potential of using mixture models, the paper proposes a plausible solution in the form of a mixture model built on deep neural networks. The model, called deep neural mixture model (DNMM), relies on a convex combination of component pdfs estimated by component-specific DNNs. The DNMM belongs to the broad family of  non-parametric pdf estimation techniques, insofar that (differently form, say, GMMs) no assumptions on the form of the individual component pdfs are made \cite{DudHar73}. Due to the learning and generalization capabilities of DNNs, the DNMM can actually learn a general form for the mixture at hand, overcoming the drawbacks of the traditional non-parametric techniques, as well. A ML training algorithm is devised, satisfying (at least numerically) a combination of hard and soft constraints required in order to guarantee a proper probabilistic interpretation of the estimated model. The resulting machine can also be seen as a novel, special case of mixture of experts \cite{Yuksel2012TwentyYO} having a specific task, a ML-based unsupervised training algorithm, and a particular probabilistic strategy for assigning credit to its individual experts.

This paper is the extended, journal version of a previous workshop communication that we presented in \cite{Trentin-NMM-ANNPR18}. Over that communication, the present article extends the treatment to DNNs, defines formally the modeling capabilities of DNMM, presents the cross-validated likelihood algorithm for DNMM model selection, and reports on the results of experiments conducted over complex, multi-dimensional pdfs (in \cite{Trentin-NMM-ANNPR18} only experiments on univariate setups are reported and analyzed).

\subsection{Overview of the paper}

The paper is organized as follows. Section \ref{sec:algorithm} is devoted to the ML-based learning algorithm for DNMMs. The family of pdfs that can actually be estimated to any degree of precision via DNMM is defined formally in Section \ref{sec:theory}. An automatic model selection algorithm that exploits the probabilistic nature of the DNMM along with the cross-validated likelihood criterion is presented in Section \ref{sec:model-selection}.  The  experimental evaluation is reported in Section \ref{sec:exp}, where the DNMM compares favorably with respect to established statistical techniques (parametric as well as non-parametric) in the task of estimating non-trivial mixtures densities having different number of components and diverse dimensionality of their definition domain. Finally, Section \ref{sec:conclusions} draws the conclusions and outlines the current research directions.

\section{DNMM: formal definition and training algorithm}
\label{sec:algorithm}
\sloppy

Hereafter we rely on the notation introduced at the beginning of Section \ref{sec:intro}, such that our goal is to estimate $p({\bf x})$ from the information encapsulated within the training sample ${\cal T}$. To this aim, we introduce a deep neural mixture model ${\tilde p}({\bf x}|W)$ having the following form:

\begin{equation}
{\tilde p}({\bf x}|W) = \sum_{i=1}^{K} c_i  {\tilde p}_i({\bf x}|W_i)
  \label{eq:NMM}
\end{equation}

\noindent where $W$ represents the whole set of parameters in the DNMM (that is $c_1, \ldots, c_K, W_1, \ldots, W_k$). The mixing parameters $c_i$ shall satisfy $c_i \in [0,1]$ for $i=1, \ldots, K$ and $\sum_{i=1}^{K} c_i =1$. The $i$-th component pdf ${\tilde p}_i({\bf x}|W_i)$ is defined, for $i=1, \ldots, K$, as

\begin{equation}
{\tilde p}_i({\bf x}|W_i) = \frac{\varphi_i({\bf x}, W_i)}{\int \varphi_i({\bf x}, W_i) d{\bf x}}
  \label{eq:component}
\end{equation}

\noindent where $\varphi_i({\bf x}, W_i)$ is the function computed by a component-specific DNN whose set of learnable parameters is $W_i$. We say that this DNN realizes the $i$-th deep neural component of the DNMM. Henceforth, plain feed-forward DNNs are assumed (with arbitrary activation functions), but the following calculations can be adapted to a variety of alternate families of DNNs in a (mostly) straightforward manner. In order for the DNMM to satisfy Kolmogorov's axioms of probability, a constraint on $\int \varphi_i({\bf x}, W_i) d{\bf x}$ shall be imposed shortly. It goes without saying that each  DNN in the DNMM has $d$ input neurons (matching the dimensionality of the feature space) and a single output neuron (yielding the value of the estimated pdf), and it is expected to have as many hidden layers as needed (the automatic model selection procedure presented in Section \ref{sec:model-selection} is suitable for architecture selection, as well). Without loss of generality for all the present intents and purposes, we assume\footnote{It is seen that, in practical applications, any data normalization approach may be applied to the training sample in order to ensure respect of this assumption.}  that the random vectors of interest lie within a compact $S \subset \mathbb{R}^d$, such that $S$ can be regarded as the definition domain of $\varphi_i({\bf x}, W_i)$ for all $i=1, \ldots, K$. In so doing, numerical integration algorithms are viable to compute $\int \varphi_i({\bf x}, W_i) d{\bf x}$, as well as the other integrals required shortly. Equation (\ref{eq:component}) reduces to ${\tilde p}_i({\bf x}|W_i) = \frac{\varphi_i({\bf x}, W_i)}{\int_S \varphi_i({\bf x}, W_i) d{\bf x}}$.

The choice of the form of the activation function $f_i(.)$ used in the output layer of the $i$-th DNN requires some precautions. Due to the very nature of pdfs, $f_i(.)$ is expected to have (at least in principle) a codomain defined as $[0, +\infty)$. This may be granted in several different ways. In this paper we use a logistic sigmoid with component-specific adaptive amplitude $\lambda_i \in \mathbb{R}^+$, that is $f_i(a_i) = \lambda_i/(1+\exp(-a_i))$ as described in \cite{Trentin2001}, where $a_i$ represents the current activation value for the output neuron of the $i$-th DNN in the mixture. Therefore, each neural component in the DNMM can stretch its output over any (component-specific) range $[0, \lambda_i)$, which, instead of being predefined by the user,  is learned from the data (just like any other parameter in $W_i$) in compliance with the nature of the corresponding component pdf\footnote{Other advantages entailed by the use of adaptive amplitudes are pointed out in \cite{Trentin2001}.}.

The DNMM training algorithm revolves around a learning rule for the mixture parameters $W$ given ${\cal T}$, such that at the end of the training process the quantity ${\tilde p}({\bf x}| W)$ results in a robust estimate of $p({\bf x})$. This is achieved by pursuing two purposes: (1) exploiting the information encapsulated in ${\cal T}$ to approximate the unknown pdf; (2) preventing the DNNs in the DNMM from developing spurious solutions, by enforcing the constraints $\int_S \varphi_i({\bf x}, W_i) d{\bf x} = 1$ for all $i=1, \ldots, K$. To this aim, a constrained algorithm is devised that builds on the stochastic gradient-ascent maximization of the point-wise likelihood ${\tilde p}({\bf x}_j | W)$ of the DNMM given the current training pattern ${\bf x}_j$. The stochastic optimization step has to be applied iteratively for $j=1, \ldots, n$. This is achieved by means of an on-line, differentiable loss function $C(.)$ defined as

\begin{equation}
  C(W,{\bf x}_j) = {\tilde p}({\bf x}_j| W)  - \rho \sum_{i=1}^{K} \frac{1}{2} \left ( 1 - \int_S \varphi_i({\bf x}, W_i) d{\bf x} \right )^2
\label{eq:criterion}
\end{equation}

\noindent to be maximized with respect to the DNMM parameters $W$ under the (hard) constraints that $c_i \in [0,1]$ for $i=1, \ldots, K$ and $\sum_{i=1}^{K} c_i =1$. The second term in the loss function is rather a ``soft'' constraint that enforces a unit
integral of ${\tilde p}_i({\bf x}, W_i)$ over $S$ for all $i=1, \ldots, K$, as sought, such that $\int_S {\tilde p}({\bf x} | W) d {\bf x} \simeq 1$. The hyper-parameter
$\rho \in \mathbb{R}^+$ balances the importance of the constraints, and
may be used in real-life applications in order to tackle potential numerical issues. The learning rule $\Delta w$ for a generic parameter $w$ in the DNMM is written as
$\Delta w = \eta \frac{\partial C(.)}{\partial w}$, where $\eta \in \mathbb{R}^+$ is the learning rate.  Different calculations are needed, depending on $w$ being either (i) a mixing coefficient, i.e. $w = c_k$, or (ii) a parameter\footnote{A connection weight, bias, adaptive amplitude, or any other trainable parameter.} within any of the DNN-based component densities. In case (i), we introduce $K$ unconstrained latent  variables $\gamma_1, \dots, \gamma_K$, and we let

\begin{equation}
c_k = \frac{\varsigma(\gamma_k)}{\sum_{i=1}^{K}\varsigma(\gamma_i)}
\label{eq:mixing-par}
\end{equation}

\noindent for $k=1,\ldots,K$, where $\varsigma(x) = 1/(1+e^{-x})$. In so doing, hereafter each $\gamma_k$ can be treated as the unknown parameter to be estimated, in place of the corresponding $c_k$. Consequently, higher-likelihood mixing parameters $c_k$ that satisfy the required constraints are implicitly yielded by the application of the learning rule. The latter can be written as:

\begin{eqnarray}
  \Delta \gamma_k & = & \eta \frac{\partial C(.)}{\partial \gamma_k}\\  \nonumber
  & = & \eta \frac{\partial {\tilde p}({\bf x}_j| W)}{\partial \gamma_k}\\  \nonumber
  & = & \eta \frac{\partial \sum_{i=1}^{K} c_i  {\tilde p}_i({\bf x}_j|W_i)}{\partial \gamma_k}\\  \nonumber
& = & \eta \sum_{i=1}^{K} {\tilde p}_i({\bf x}_j| W_i)  \frac{\partial}{\partial \gamma_k}\left(\frac{\varsigma(\gamma_i)}{\sum_{\ell=1}^{K}\varsigma(\gamma_\ell)}\right)  \\ \nonumber
& = &  \eta \left\{ {\tilde p}_k({\bf x}_j| W_k) \frac{\varsigma^\prime(\gamma_k)}{\sum_{\ell=1}^{K} \varsigma(\gamma_\ell)} - \sum_{i=1}^{K} {\tilde p}_i({\bf x}_j| W_i) \frac{\varsigma(\gamma_i)\varsigma^\prime(\gamma_k)}{[ \sum_\ell \varsigma(\gamma_\ell)]^2} \right\}  \\ \nonumber
& = & \eta \frac{\varsigma^\prime(\gamma_k)}{\sum_\ell \varsigma(\gamma_\ell)}
\left\{ {\tilde p}_k({\bf x}_j| W_k) -  {\tilde p}({\bf x}_j| W) \right\}
\label{eq:delta-mixing-par}
\end{eqnarray}

Next, let us focus on scenario (ii), i.e. let $w$ be a parameter within one of the DNNs. Taking the partial derivative of $C(W,{{\bf x}_j})$ with respect to $w$ requires calculating the derivatives of the first and the second terms in the right-hand side of Equation (\ref{eq:criterion}), respectively. For notational convenience, hereafter we assume that $w$ belongs to the generic $k$-th DNN. The partial derivative of the first term can be written as

\small
\begin{eqnarray}
  \label{eq:partialC}
  \frac{\partial {\tilde p}({\bf x_j}| W)}{\partial w}  & = & \frac{\partial}{\partial w} \sum_{i=1}^{K} c_i {\tilde p}_i({\bf x_j} | W_i)   \\ \nonumber
  & = & \frac{\partial}{\partial w} \{c_k {\tilde p}_k({\bf x_j} | W_k)\}   \\ \nonumber
  & = &  c_k \frac{\partial}{\partial w} \left\{ \frac{\varphi_k({\bf x}_j, W_k)}{\int_S \varphi_k({\bf x}, W_k)d{\bf x}} \right\} \\ \nonumber
  & = & c_k \left\{ \frac{1}{\int_S \varphi_k({\bf x}, W_k)d{\bf x}} \frac{\partial \varphi_k({\bf x}_j, W_k)}{\partial w} - \frac{{\tilde p}_k({\bf x}_j, W_k)}{\int_S \varphi_k({\bf x}, W_k)d{\bf x}} \frac{\partial}{\partial w} \int_S \varphi_k({\bf x}, W_k)d{\bf x} \right\}  \\ \nonumber
  & = & \frac{c_k}{\int_S \varphi_k({\bf x}, W_k)d{\bf x}} \left\{  \frac{\partial \varphi_k({\bf x}_j, W_k)}{\partial w} - \frac{ \varphi_k({\bf x}_j, W_k)}{\int_S \varphi_k({\bf x}, W_k)d{\bf x}} \int_S \frac{\partial \varphi_k({\bf x}, W_k)}{\partial w} d{\bf x} \right\}
\end{eqnarray}
\normalsize

\noindent where we used Leibniz rule in the last step of the
calculations. In passing, it is worth noting that Equation (\ref{eq:partialC}) is the mathematical statement of the very rationale behind the specific impact that the current training pattern ${\bf x}_j$ has on the learning process for distinct neural components of the DNMM. In fact, the amount of parameter change $\Delta w$ is proportional to the probabilistic ``credit'' $c_k$ of the neural component at hand.  Moreover, the quantities within brackets in Equation (\ref{eq:partialC}) depend on the value of the $k$-th DNN output over ${\bf x}_j$, as well as on its derivative. If, at any given time during the training process, $\varphi_k(.)$ does not change significantly in a neighborhood of ${\bf x}_j$ (e.g. if ${\bf x}_j$ lies in a high-likelihood plateau or, vice versa, in a close-to-zero plateau of $\varphi_k(.)$) then the contribution of the first quantity within brackets is neglectable. Moreover, if $\varphi_k({\bf x}_j) \simeq 0$ then the second term within brackets turns out to be neglectable, as well.  To the contrary, the contribution of ${\bf x}_j$ to the parameter adaptation of $k$-th component DNN turns out to be paramount if $\varphi_k(.)$ returns a high likelihood over ${\bf x}_j$ and significant variations in its surroundings.

Next, Leibniz rule is used again in the calculation of the derivative
of the second term in the right-hand side of Equation (\ref{eq:criterion}), which can be written as

\begin{eqnarray}
\label{eq:constraint}
\frac{\partial}{\partial w} \left\{  \rho \sum_{i=1}^{K} \frac{1}{2} \left ( 1 - \int_S \varphi_i({\bf x}, W_i) d{\bf x} \right )^2 \right\} =  \\ \nonumber
= \frac{\partial}{\partial w}  \left\{ \frac{\rho}{2} \left ( 1 - \int_S \varphi_k({\bf x}, W_k) d{\bf x} \right )^2 \right\} =  \\ \nonumber
= - \rho \left ( 1 - \int_S \varphi_k({\bf x}, W_k) d{\bf x} \right ) \frac{\partial}{\partial w} \int_S \varphi_k({\bf x}, W_k) d{\bf x}   \\ \nonumber
= - \rho \left ( 1 - \int_S \varphi_k({\bf x}, W_k) d{\bf x} \right )  \int_S  \frac{\partial \varphi_k({\bf x}, W_k)} {\partial w} d{\bf x}.
\end{eqnarray}

\noindent In order to compute the right-hand side of Equations (\ref{eq:partialC}) and (\ref{eq:constraint}) we need suitable algorithms for the computation of $ \frac{\partial \varphi_k({\bf x}_j, W_k)}{\partial w}$, $\int_S \varphi_k({\bf x}, W_k)d{\bf x}$, and $ \int_S  \frac{\partial}{\partial w} \varphi_k({\bf x}, W_k)d{\bf x}$. As regards the computation of $ \frac{\partial \varphi_k({\bf x}_j, W_k)}{\partial w}$ we proceed as in traditional BP, or as in \cite{Trentin2001} in case $w = \lambda_k$. As for the integrals, deterministic numerical quadrature integration techniques (e.g., Simpson's method, trapezoidal rule, etc.) are viable only if $d =1$. In fact, in terms of computational burden, they cannot realistically scale up to higher dimensions ($d \geq 2$). This is all the more critical in the light of the fact that  $ \int_S  \frac{\partial}{\partial w} \varphi_k({\bf x}, W_k)d{\bf x}$ shall be iteratively computed for each parameter of each DNN in the DNMM. Moreover, deterministic integration methods do not exploit at all the nature of the function under integration. Roughly speaking, herein the integrand is expected to be an ``approximation'' of  the pdf (say, $p_k({\bf x})$) that explains the distribution of that specific sub-sample of ${\cal T}$ that is drawn from the $k$-th component of the mixture. To the contrary, accounting for the pdf of the data should drive the integration algorithm towards integration points that cover ``interesting'' regions\footnote{That is, regions having high component-specific likelihood.} of the domain of the integrand. For these reasons, we apply a component-oriented version of the approach we presented in \cite{Trentin_Soft_NEPL}. The resulting approach can be seen as an instance of Markov chain Monte Carlo \cite{MCMC-03}. It is a non-deterministic, multi-dimensional integration technique which accounts for the component-specific probability distribution of the data. Let $\phi_k({\bf x})$ denote the integrand at hand (i.e., $\varphi_k({\bf x}, W_k)$ or $\frac{\partial \varphi_k({\bf x}, W_k)}{\partial w}$). Monte Carlo with importance sampling \cite{ImportanceSampling} yields an approximation of the integral of  $\phi_k({\bf x})$ over $S$ in the form $\int_S \phi_k({\bf x}) d{\bf x} \simeq \frac{V(S)}{m} \sum_{\ell=1}^m  \phi_k(\dot{{\bf x}}_\ell)$  where $m$ properly sampled integration points $\dot{{\bf x}}_1, \ldots, \dot{{\bf x}}_m$ are used. Sampling of the $\ell$-th integration point $\dot{{\bf x}}_\ell$ (for $\ell =1, \ldots, m$) is obtained by drawing it at random from the mixture pdf $p^{(k)}_u({\bf x})$ defined as

\begin{equation}
  p_u^{(k)}({\bf x}) = \alpha(t) u({\bf x}) + (1 - \alpha(t)){\tilde p}_k({\bf x}|W_k)
  \label{eq:mixture}
\end{equation}

\noindent where $u({\bf x})$ is the uniform distribution over $S$, and $\alpha:\mathbb{N} \rightarrow (0,1)$ is a decaying function of the number $t$ of the DNMM training epochs\footnote{In the present context, a training epoch is a completed re-iteration of Equations (\ref{eq:partialC}) and (\ref{eq:constraint}) for all the parameters of the DNMM over all the observations in ${\cal T}$.} for $t=1, \ldots, T$, such that $\alpha(1) \simeq 1.0$ and $\alpha(T) \simeq 0.0$. As in \cite{Trentin_Soft_NEPL} we let  $\alpha(t) = 1/ (1 + e^{\frac{t/T - 1/2}{\theta} })$, where $\theta$ is a hyperparameter. Equation (\ref{eq:mixture})  is such that the importance sampling mechanism it implies  accounts for the (estimated) component density ${\tilde p}_k({\bf x} | W_k)$ of the $k$-th latent subsample of the data, therefore respecting the natural distribution of such sub-population and focusing integration on the relevant integration points (i.e., the points having high component-specific likelihood).  At the same time, since the estimates of this component pdfs are unfit during the early stage of the DNMM training, Equation (\ref{eq:mixture}) prescribes a (safer) sampling from a uniform distribution at the beginning (practically behaving like a plain Monte Carlo algorithm). As long as the robustness of the DNMM estimate increases, i.e. as $t$ increases, sampling from  ${\tilde p}_k({\bf x} | W_k)$ replaces progressively the sampling from $u({\bf x})$, ending up in purely non-uniform importance sampling. The form of $\alpha(t)$ is defined accordingly. Since $\varphi_k({\bf x}, W_k)$ is intrinsically non-negative, for $t \rightarrow T$ the sampling occurs substantially from $\left\vert{\varphi_k({\bf x}, W_k)}\right\vert /\int_S \left\vert{\varphi_k({\bf x}, W_k)}\right\vert d{\bf x}$, that is a sufficient condition for granting that the variance of the estimated integral vanishes and the corresponding error goes to zero \cite{Ohl:1998jn}.

Sampling  from $p^{(k)}_u({\bf x})$ requires a viable technique for sampling from the $k$-th DNN in the DNMM. A variant of Markov chain Monte Carlo, namely the Metropolis--Hastings (M-H) algorithm \cite{Newman1999book}, is exploited in this paper. M-H is robust to the fact that, during training, $\varphi_k({\bf x}, W_k)$ may not be properly normalized but it is proportional by construction to the corresponding pdf estimate (which is normalized properly, instead, by definition) \cite{Newman1999book}. Due to its efficiency and ease of
sampling, a multivariate logistic pdf with radial basis, having location ${\bf x}$ and scale $\sigma$, is used as the {\em proposal} pdf $q({\bf x}^\prime | {\bf x})$ required by M-H to generate a new candidate sample ${\bf x}^\prime = (x^\prime_1, \ldots, x^\prime_d)$ from the current sample ${\bf x} = (x_1, \ldots, x_d)$. Formally, such a proposal pdf is defined as $q({\bf x}^\prime | {\bf x}, \sigma) = \prod_{i=1}^d \frac{1}{\sigma} e^{(x^\prime_i - x_i)/\sigma} (1 + e^{(x^\prime_i -x_i)/\sigma})^{-2}$ which can be easily sampled via the inverse transform sampling technique. The hyperparameters needed (i.e., the scale $\sigma$ of the proposal pdf and the burn-in period for M-H) are fixed empirically as part of the overall model selection process (see Section \ref{sec:model-selection}).

\section{Class of pdfs that can be modeled accurately via DNMMs}
\label{sec:theory}

It is seen that not all pdfs can be estimated in an accurate manner via DNMMs. A couple of simple, univariate examples should be more than enough to convince the reader of this: the Dirac's Delta (which is not continuous and not bounded) and the standard exponential pdf (defined as $p_e(x) =0$ if $x<0$, $p_e(x) = \exp(-x)$ otherwise). In \cite{TrentinLC18} we introduced the class of {\em nonpaltry pdfs}, that basically embrace all the ``interesting'' pdfs that are of practical interest (and, that can  be estimated to any degree of precision by means of Parzen Neural Networks). The arguments handed out by \cite{TrentinLC18} can be extended to DNMMs, as well. The formal definition of this class of function goes as follows. Let $S$ be a compact subset of $I_d$, where $I_d = [0,1]^d$ (the symbol $S$ has the same practical meaning it had in the previous section). A continuous pdf $p:\mathbb{R}^d \rightarrow \mathbb{R}$ is said to be {\em nonpaltry} over $S$ if  $\int_{S} p({\bf x}) d {\bf x} = 1$. Accordingly, we formally define the class ${\cal P}(S)$ of nonpaltry pdfs over $S$ as the set of all the density functions that are nonpaltry over $S$. In passing, note that mixture densities built on nonpaltry component densities are, in turn, nonpaltry. Theorem 3.1 in \cite{TrentinLC18} states that for any nonpaltry pdf $p(\cdot)$ over $S$ at least one DNN exists that computes $\phi^{(\epsilon)}(\cdot)$ s.t. $\| \phi^{(\epsilon)}(\cdot) - p(\cdot) \|_{\infty,S} < \epsilon$ for any $\epsilon \in \mathbb{R}^+$. This (trivially) implies that at least one DNMM (with $c$ =1) exists that computes $\phi_{DNMM}^{(\epsilon)}(\cdot)$ s.t. $\| \phi_{DNMM}^{(\epsilon)}(\cdot) - p(\cdot) \|_{\infty,S} < \epsilon$ for any $\epsilon \in \mathbb{R}^+$. Roughly speaking, we can estimate any nonpaltry pdf to any degree of precision by means of a DNMM.

\section{Automatic likelihood-driven model selection strategy}
\label{sec:model-selection}

The former section presented a modeling (or, approximation) capability of the family of pdfs realized by DNMMs, regardless of the actual convergence (during training) of a specific DNMM to the unknown pdf $p(\cdot)$ underlying the data ${\cal T}$ at hand. In practice, roughly speaking, the larger the cardinality of the training sample ${\cal T}$, the closer to $p(\cdot)$ it tends to become (in probability) the estimate realized by the DNMM. This means that, for sufficiently large values of $n$, the likelihood of the DNMM given a separate validation set ${\cal V}$ (independently drawn from $p(\cdot)$, as well) tends to ``increase\footnote{The term ``increase'' is used hereby with its qualitative meaning, that is to indicate a trend and not a strictly monotonic behavior: in fact, the likelihood of the DNMM does not necessarily increase monotonically as the training process proceeds.}'' with $n$. This informal reasoning suggests the adoption of the maximum-likelihood criterion as an evaluation function that can be used throughout the DNMM model selection process. This is made viable by the very probabilistic nature of the DNMM, and it is not available to generic, non-probabilistic DNNs or mixtures of experts.

In short, an iterative model/hyperparameter search strategy may be developed as follows.To fix ideas, let us first consider the issue of fixing the number of neurons in a certain hidden layer of a DNN within the DNMM\footnote{Note that there is no assumption in the definition of DNMM that forces the practitioner to use the same number of layers and/or neurons in the different DNNs realizing the components of the mixture, although this may be a computationally realistic choice.}. We start with an initial number of hidden neurons, and we increase it by $u$ units (for a certain integer $u$, e.g. $u = 1$) at each of the following steps of the selection procedure. At each iteration (starting from the initialization) the DNMM is trained on ${\cal T}$, and the likelihood of the resulting model given ${\cal V}$ is evaluated. The procedure is iterated for as long as this likelihood increases. We stop searching when, for $\tau$ consecutive iterations, the relative gain in likelihood between consecutive steps is lower than a fixed percentage $\nu \%$. The minimum description length principle is then applied \cite{MDL-J1978465}, selecting the simplest model amongst those (explored throughout the last $\tau$ iterations) which yielded comparable values of the corresponding validation likelihood given ${\cal V}$ (where ``comparable'' means that their difference is below a user-defined threshold). It is worth noting that this approach relies on an instance of the so-called cross-validated likelihood criterion \cite{Rust1985}. The same incremental strategy is viable for selecting the number of layers within the DNNs in the mixture (e.g., starting from a single hidden layer and deepening the architecture by adding one more layer at each step).

The present cross-validated likelihood search strategy can be applied, in a straightforward manner, to the problem of selecting the other hyperparameters controlling the behavior of the DNMM learning process. In order to limit the computational burden of such an automatic model selection process, random search of the values of the hyperparameters (each choice characterized by the corresponding likelihood given ${\cal V}$) is recommended instead of going through the step-by-step incremental approach we proposed for selecting the umber of neurons.

\section{Experiments}
\label{sec:exp}
\sloppy

Experiments are conducted on random samples drawn from multimodal mixtures having known form. Section \ref{sec:univariate} reports on the results obtained on univariate data using mixtures of 3-layer DNNs, while Section \ref{sec:multivariate} presents the results of experiments involving multivariate data (of different dimensionalities) using deeper DNMMs.

\subsection{Experiments with univariate data}
\label{sec:univariate}

In the present, univariate case, the random samples were drawn from mixtures $p(x)$ of $c$ Fisher-Tippett pdfs, i.e. $p(x) = \sum_{i=1}^{c}\frac{P_i}{\beta_i} \exp \left(-\frac{x-\mu_i}{\beta_i}\right) \exp \left\{ -\exp \left(-\frac{x-\mu_i}{\beta_i}\right) \right\}$. The mixing coefficients $P_1, \ldots, P_c$ were drawn at random from the uniform distribution over $[0.1]$ and normalized in such a way that $\sum_{i=1}^c  P_i = 1$.  The component densities of the Fisher-Tippett mixture are identified by their scales $\beta_i$ and locations $\mu_i$, for $i=1, \ldots, c$. The scales were drawn at random from the uniform distribution over $(0.01, 0.9)$, and the locations were randomly and uniformly distributed over $(0,10)$. Estimation tasks were created using $1200$ random patterns each, drawn from $p(x)$, and a variable number $c$ of component densities (namely $c=5$, $10$, $15$ and $20$). Each $c$-specific sample was split into a training set (with $n = 800$ patterns) and a validation set (the remaining $400$ patterns).  The integrated squared error (ISE) between $p(x)$ and its estimate ${\tilde p}(x)$, i.e. $\int \left(p(x) - {\tilde p}(x)\right)^2 d x$, is adopted in order to assess the performance of the different estimators. Simpson's method was applied to compute numerically the ISE.

In the present, illustrative setup we used DNMMs involving DNNs with a 3-layer architecture (input layer, hidden layer of  9 neurons, output layer). Sigmoid activation functions were used in the hidden and the output layers, having layer-wise  \cite{Trentin2001} adaptive $\lambda$. All the DNMM parameters underwent random initialization over zero-centered intervals, except for the values of $\lambda$ (that were initialized to $1$) and the mixing parameters that were initialized as $c_i = 1/K$ for $i=1, \ldots, K$. As in \cite{Trentin_Soft_NEPL} we used a function $\alpha(t)$ having $\theta = 0.07$, and we relied on $m=400$ integration point. The latter ones were sampled at the beginning of each training epoch using a scale $\sigma = 9$ for the logistic proposal pdf in M-H. The burn-in period of the Markov chain in M-H was stretched over the first $500$ states. The other hyperparameters of the DNMM were fixed via random-search based on the cross-validated likelihood criterion presented in Section \ref{sec:model-selection}. No normalization was applied to the values of the input data. Table \ref{tab:comp} reports the results. DNMMs having different values of $K$ (that is $K=4$, $8$ and $12$) were tested and compared with 8-GMM, 16-GMM, and 32-GMM (initialized via $k$-means and refined via iterative ML \cite{DudHar73}), $k_n$-NN with unbiased $k_n = 1 \sqrt{n}$ \cite{DudHar73}, and Parzen Window (PW) with standard width $h_n = 1/\sqrt{n}$ of its Gaussian kernels \cite{DudHar73}.

\begin{table}[htp]
\caption{Estimation of the Fisher-Tippett  mixture $p(x)$ (with $n=800$) in terms of integrated squared error as a function of the number $c$ of the Fisher-Tippett  component densities. Best results are shown in boldface. ({\em Legend}: $k$-GMM and $k$-DNMM denote the GMM and the DNMM with $k$ components, respectively).}
\label{tab:comp}
\vspace{0.2cm}
\begin{small}
\begin{tabular}{l l @{\hskip 0.4cm} l @{\hskip 0.4cm} l @{\hskip 0.4cm} l|l}
\hline
$c$  & 5 & 10 & 15 & 20 & {\em Avg. $\pm$ std. dev.}\\
\hline
8-GMM & $9.60\textrm{e}-3$ & $1.12\textrm{e}-2$ & $4.57\textrm{e}-2$ & $7.99\textrm{e}-2$ &  $(3.66 \pm 2.89) \textrm{e}-2$ \\
16-GMM & $6.33\textrm{e}-3$ & $9.29\textrm{e}-3$ & $3.78\textrm{e}-2$ & $4.24\textrm{e}-2$ & $(2.40 \pm 1.63) \textrm{e}-2$  \\
32-GMM & $7.15\textrm{e}-3$ & $9.82\textrm{e}-3$ & $2.41\textrm{e}-2$ & $3.03\textrm{e}-2$ & $(1.78 \pm 0.97) \textrm{e}-2$  \\
$k_n$-NN & $6.54\textrm{e}-3$ & $8.70\textrm{e}-3$ & $2.03\textrm{e}-2$ & $2.36\textrm{e}-2$ & $(1.48 \pm 0.73)\textrm{e}-2$  \\
PW & $6.02\textrm{e}-3$ & $8.94\textrm{e}-3$ & $2.14\textrm{e}-2$ & 1$.98\textrm{e}-2$ & $(1.40 \pm 0.67)\textrm{e}-2$  \\
4-DNMM & $6.41\textrm{e}-3$ & $7.06\textrm{e}-3$ & $1.09\textrm{e}-2$ & $1.40\textrm{e}-2$ & $(9.59 \pm 3.07) \textrm{e}-3$  \\
8-DNMM & ${\bf 5.89\textrm{e}-3}$ & ${\bf 6.02\textrm{e}-3}$ & $8.11\textrm{e}-3$ & $1.01\textrm{e}-2$ & $({\bf 7.53} \pm 1.73) \textrm{e}-3$ \\
12-DNMM & $6.38\textrm{e}-3$ & $6.27\textrm{e}-3$ & ${\bf 8.05\textrm{e}-3}$ & ${\bf 9.64\textrm{e}-3}$ & $(7.59  \pm {\bf 1.38}) \textrm{e}-3$
\end{tabular}
\end{small}
\end{table}

The Table shows that, regardless of the model used, the corresponding ISE increase as a function of $c$, as expected. The DNMMs improve systematically (and, significantly) over the statistical techniques, whatever the value of of $c$. On average, both the 8-DNMM and the 16-DNMM yield a $46\%$ relative ISE reduction over the PW (the latter turns out to be the most robust non-neural estimation algorithm). Welch's t-test (used in order to account for the different variances of the models) \cite{Japkowicz:2011} returns a level of confidence $> 90\%$ on the statistical significance of the gap between the 8- (or, 12-) DNMM and the PW. Furthermore, the DNMMs  turn out to be the stablest models overall, as proved by the values of the corresponding standard deviations (last column of the table). This is evidence of the fact that the estimation accuracy offered by the DNMMs is less sensitive to the complexity of the underlying Fisher-Tippett mixture (i.e., to $c$) than the accuracies yielded by the traditional statistical techniques. Finally, differences in terms of ISE are observed among the DNMMs depending on $k$. Nonetheless, differences between the 8-DNMM and the 12-DNMM turn out to be mild, and they depend on the complexity of the underlying pdf to be estimated (at least to some extent), as expected.

\subsection{Experiments with multivariate data}
\label{sec:multivariate}

Hereafter we use multivariate data drawn from mixtures of generalized extreme value distributions (GEV) \cite{castillo2004extreme} with null slope and having the following parametric form:

\begin{equation*}
{\text{\em m-GEV}}({\bf x};{\bf \Theta} ) = \sum_{k=1}^{c_{T}} \frac{1}{c_{T}} \prod_{i=1}^{d} \frac{1}{\beta_{ki}} \exp
\left(-\frac{x_i -\mu_{ki}}{\beta_{ki}}\right) \exp \left\{ -\exp
\left(-\frac{x_i - \mu_{ki}}{\beta_{ki}}\right) \right\}
\label{eq:GEV}
\end{equation*}

\noindent where the GEV parameters ${\bf \Theta} = (d, c_{T}, {\bm \mu}_1, \ldots, {\bm \mu}_{c_T}, {\bm \beta}_1, \ldots, {\bm \beta}_{c_T})$ are defined as follows:

\begin{enumerate}

\item $d$ is the dimensionality of the feature space as usual. Therefore, the generic real valued random vector ${\bf x} \in \mathbb{R}^d$ can be written as ${\bf x} = (x_1, \ldots, x_d)$;

\item $c_T$ is the number of component GEVs in the mixture, and $c_T = c^d$ where $c$ denotes the (fixed) number of diverse modes of ${\text{\em m-GEV}}({\bf x};{\bf \Theta} )$ along each one of the dimensions in the definition domain of the pdf;

\item ${\bm \mu}_k = (\mu_{k1}, \ldots, \mu_{kd})$ and ${\bm \beta}_k = (\beta_{k1}, \ldots, \beta_{kd})$ are the location and the scale vectors of $k$-th component GEV, respectively. They are defined in such a manner that the set $\{({\bm \mu}_k, {\bm \beta}_k) | k=1, \ldots, c^d\}$ embraces all possible combinations of $c$ dimension-specific parameters $({\bar \mu}_{i1}, {\bar \beta}_{i1}), \ldots, ({\bar \mu}_{ic}, {\bar \beta}_{ic})$  for $i=1, \ldots, d$ (for a total of $c_T = c^d$ combinations).

\end{enumerate}

\noindent  For each $i=1,\ldots,d$ and $j=1,\ldots,c$, the values ${\bar \mu}_{ij}$ are random quantities, independently and uniformly drawn from the interval $(0.1,0.9)$, and the iid random values ${\bar \beta}_{ij}$ are drawn from the interval $(0.03,0.05)$. In the following, numerical integration is computed over the interval $S = [0, 1.1]^d$.

\begin{table}[htp]
\caption{Estimation of the multivariate mixture of generalized extreme value distributions  ${\text {\em m-GEV}}({\bf x};{\bf \Theta})$ (with $n=800$) in terms of relative ($\%$) ISE reduction over the baseline, as a function of the number $C_T$ of the component densities of the ${\text {\em m-GEV}}({\bf x};{\bf \Theta})$ and of the dimensionality $d$ of the feature space.}
\label{tab:mult}
\vspace{0.2cm}
\begin{tabular}{l l @{\hskip 0.4cm} l @{\hskip 0.4cm} l @{\hskip 0.4cm} l|l}
\hline
 & $C_T = 4$ & $C_T = 9$ & $C_T = 16$ & $C_T = 25$ & {\em Avg.}\\
\hline
$d = 2$ & 11.31 & 10.52 & 7.38 & 9.02 & $9.56 \pm 1.50$  \\
$d = 4$ & 8.44 & -0.07 & 5.75 & 7.01 & $5.28 \pm 3.23$  \\
$d = 6$ & 4.98 & -1.64 & 5.80 & 8.13 &  $4.32 \pm 3.63$ \\
$d = 8$ & 6.20 & 7.34 & 8.63 & 8.00 & $7.54 \pm 0.90$ \\
\hline
{\em Avg.} & $7.73 \pm 2.41$ & $4.04 \pm 5.05$  & $6.89 \pm 1.20$ & $8.04 \pm 0.71$ & ${\bf 6.67 \pm 3.29}$
\end{tabular}
\end{table}

We generated different estimation tasks involving {\em m-GEVs} with an increasing number of components ($C_T = 4, 9, 16$ and $25$, respectively) and an increasing dimensionality of the feature space ($d = 2, 4, 6$ and $8$, respectively). For each combination of such values for $C_T$ and $d$, a sample of as many as $1200$ patterns was randomly drawn from ${\text {\em m-GEV}}({\bf x};{\bf \Theta})$. Each such data sample was split into a training set ($800$ random vectors) and a validation set (the remaining $400$ patterns), as we did before. The validation set was used to select the architectures and hyperparameters for the algorithms via ML-based random search (applying the model selection procedure presented in Section \ref{sec:model-selection}) on a $d$- and $c_T$-specific basis. This cross-validated likelihood criterion was exploited also in order to fix a statistical baseline for assessing the performance of the DNMM, as follows. First, three established and popular statistical estimation techniques suitable for multivariate, multimodal pdf estimation were evaluated, namely the PW, the $k_n$-NN, and the GMM (the latter was evaluated several times, for an increasing number of component Gaussians ranging from 4 to 32). In PW we used Gaussian kernels with bandwidth $h_n = h_1/\sqrt{n}$, where $h_1$ was selected over the $[10^{-1},10^1]$ range.  For the $k_n$-NN we let  $k_n = k_1 \sqrt{n}$, as usual, with $k_1$ selected via cross-validated ML. As for the GMMs, the $k$-means algorithm was used for the initialization of the parameters, and the latter were refined via iterative ML re-estimation \cite{DudHar73}. After ML selection of suitable hyperparameters, the best amongst these statistical models (in terms of likelihood on the validation set) was used for fixing the baseline (for the specific combination of $d$ and $c_T$ at hand, then the whole procedure was repeated for the other combinations). The baseline was quantified as the ISE with respect to the true {\em m-GEV} underlying the data sample for those specific values of $d$ and $c_T$. Next, a DNMM was selected (again, based on the procedure presented in Section \ref{sec:model-selection}). Selection involved fixing the architecture of the DNNs (number of layers between 3 and 6, number of neurons per layer, type of activation functions) and the hyperparameters used for training the DNMM. The same architectures and hyperparameters were applied for all the DNNs in the DNMM at any step of the random search process. DNMMs were initialized as in the previous Section, and the same parameters used therein for the numeric integration (burn-in period, form of $alpha(t)$, $m$, etc.) were used herein. Eventually, the ISE offered by the best after-training DNMM selected this way was evaluated w.r.t. the true {\em m-GEV}, and compared with the baseline. In so doing, the results (reported in table \ref{tab:mult}) can be expressed in terms of the percentage of relative ISE reduction ($\%$) yielded by the best DNMM over the baseline\footnote{To fix ideas, a $0\%$ relative ISE reduction means that the DNMM equals the baseline, a negative relative ISE reduction means that the DNMM performance worsened the baseline, while positive reductions represent gains over the baseline (the larger, the better: ideally, in the limit case, a $100\%$ relative ISE reduction would be observed if the DNMM could result in a perfect model of the true {\em m-GEV} in presence of an imperfect estimation baseline).}. Expressing the results of the experiments in this way has two significant advantages: (1) keeping the amount of values to be read and compared, and the size of the table, to a human-readable scale, and (2) conveying the sense of the gap between the DNMM and any other established statistical estimation technique.

It is seen that, in a quasi-systematic manner, the DNMM improves significantly over the statistical approaches. Only in two cases out of 16 the ISE yielded by the DNMM is just in line with (or, slightly worse than) the baseline. Overall, the DNMM offers an average $6.67 \pm 3.29$ relative ISE reduction over its closest competitor (that is, the baseline). The statistical significance of the improvement offered by the DNMM over the baseline (and, in turn, over any statistical technique) computed via Welch's t-test is high, namely $\geq 90\%$. Moreover, the values of the std. deviation (last column and last row of the table) show that the DNMM results in a substantially stable behavior overall.

We cannot help but pointing out that the gap between the DNMM and the best statistical estimates is not as wide as in the univariate setup. A possible rationale behind the reduced gap is likely to be found in the model selection procedure we applied: the random-search strategy applied within the cross-validated ML model selection algorithm over the multivariate data (much needed, in place of a more accurate grid search, in order to keep the overall computational burden to a realistic scale) turns out to be capable of finding suitable, yet quite sub-optimal architectures and hyperparameters. This affects significantly the model selection in the DNMM (insofar that the latter has many hyperparameters to be fixed, each choice resulting in profoundly different models and, therefore, performances), while the issue is much less severe in the case of the statistical estimators (which basically involve selecting just a few hyperparameters, even one only in the PW and the $k_n$-NN cases). Another reason for the reduced gap between the DNMM and the baseline lies in our inheriting from the univariate setup the hyperparameters for M-H: it is seen that the effectiveness of (say) a certain number $m$ of integration points (which is suitable to univariate data) reduces significantly as the dimensionality of the feature space increases (in fact, many more integration points would be needed in order to obtain a representative coverage of the multidimensional space that could grant an accurate numerical computation of the integral of the multivariate pdf). Nevertheless, also in the present setup the results are evidence of the improved estimation and modeling capabilities of the DNMM even in complex multivariate density estimation tasks.

\section{Conclusions}
\label{sec:conclusions}

The paper presented the DNMM as a sound DNN-based model for multivariate pdf estimation that satisfies Kolmogorov's axioms. The DNMM overcomes all the usual shortcomings of traditional statistical techniques. In particular, its nonparametric modeling capabilities allow for realizing flexible component densities that are not constrained to have a predefined form, as it happens in parametric statistics. Moreover, the DNMM does not suffer from the limitations of nonparametric statistical estimators, as well, insofar that it is a learning machine (instead of a memory-based algorithm), which entails generalization, smooth and compact solutions, as well as reduced time- and space-complexity. Furthermore, the probabilistic nature of the DNMM allows for the exploitation of the cross-validated likelihood criterion in order to carry out the model selection relying on the ML of the hyperparameters given the validation data.

The experimental results show significant improvements yielded by the DNMM over the baseline of the statistical techniques on both the univariate and the multivariate data drawn from several complex pdfs. Future work  revolves around the initialization procedure: non-uniform initialization of the mixing coefficients may turn out to be helpful in breaking potential ties, and initializing the individual DNN-based components via supervised learning of a subset of the components of a pre-es`<timated reference mixture model (i.e., a GMM) may improve the random initialization of the DNNs parameters.

\bibliographystyle{spmpsci}


\section{Introduction}

Topological Data Analysis (TDA) has recently gained a great amount of popularity in the machine learning community. Apart from Persistent Homology, one of the more prominent methods has been the Mapper, originally proposed by \cite{singh2007topological} due to its ability to provide a visualization of high-dimensional data. In short, the Mapper algorithm, as explained in the original paper (see \cite{singh2007topological}) consists of a few steps. Given a point cloud $X \in \mathbb{R}^d$, one uses a filtration function $f : \mathbb{R}^d \rightarrow \mathbb{R}^m$, where $m < d$. This function, as stated by \cite{singh2007topological}, \textit{"can be a function which reflects geometric properties of the data"}. Once reduced using the filtrations, the points are clustered using a clustering method of choice and mapped into overlapping bins forming vertices connected if they have common points, producing a graph that serves as an approximation of the Reeb graph (see \cite{munch2015convergence}).

One of the difficulties with the Mapper is however the choice of parameters such as deciding on the overlap between the bins and design of the filtering function $f$. These choices determine the construction of the resulting graph, and which filtering function that best captures the topology is not always entirely clear. While the choice of function allows for certain freedom, it might also pose as a barrier of entry for the regular practitioner and a difficulty in generating useful representations when using the Mapper algorithm. The optimal is of course to reduce the dimensionality in a way such that the local and the global topology of the data is preserved, while the summarization of the data lose the minimum amount of information. However, as \cite{kumari2020shapevis} points out, different filtering functions capture different things. \cite{carriere2018statistical} mention that the two most common choices are the eccentricity, i.e., the maximum distance to any other point in the point cloud, and the eigenfunctions of the covariance matrix as used in PCA. Other methods such as SVD (\cite{lum_singh_lehman_ishkanov_vejdemo-johansson_alagappan_carlsson_carlsson_2013}) as well as the kernel density estimation (\cite{singh2007topological}) have been extensively used. These do however all come with their drawbacks. SVD and PCA do not capture non-linear relations very well and might miss important non-linear properties of the data, although they preserve distance. Non-linear techniques such as t-SNE (\cite{maaten2008visualizing}) and IsoMAP (\cite{tenenbaum2000global}) are better at preserving local distances but lose the global perspective as pointed out by \cite{kumari2020shapevis}. The kernel density estimator relies on a bandwidth $\epsilon$, set by the user, introducing another parameter to tune. The choice of filtering function is therefore a choice for the user where trade-offs have to be made. Finding a filtering function without the aforementioned drawbacks and which could simultaneously capture both the global and local topologies would therefore be highly beneficial, since it would produce Mapper graphs representing the shape of the data more accurately.

\section{Topological Autoencoders as a filtering function}

We propose the idea that a candidate for a filtering function preserving global and local topology is the topological autoencoder (TAE), introduced by \cite{moor2019topological}. These are autoencoders (AEs) regularized with a Persistent Homology loss to preserve the topology. As seen in their paper (see \cite{moor2019topological}), the TAEs manage to preserve the shape and form of the data in lower dimensions to a certain extent, and manages to capture the global as well as local topology apart from other methods. As they also note, the TAE can learn multiple manifolds, which is a difficult task for many methods. When this topology is preserved for a distance metric of choice, we can derive a compact representation in lower dimensions that better represents the data's shape. TAEs do of course come with a longer training time, a risk of overfitting, neural architecture choices and a need for larger datasets. It does however provide a bridge between topology and deep learning.

\section{Initial experiments}

To evaluate whether TAEs result in more representative graphs than other filtering functions, we are looking at how the Mapper maps data belonging to different manifolds into vertices with different filtering functions. As a start, we use the synthetic data constructed by \cite{moor2019topological}, consisting of ten 100-dimensional hyperspheres inside one larger 100-dimensional hypersphere in 101 dimensions. We scale the data and use DBSCAN (see \cite{ram2010density}) as clustering method and set its $\varepsilon$ to 4 for all comparisons and the euclidean distance as the distance metric. For the TAE, we reproduce the experiments from \cite{moor2019topological}, with the exact same parameter choices as in their best runs. We use the same training data to train the representations and construct the Mapper graph, to then map the test data into the Mapper graph.

Our initial experiments\footnote{Code available at \url{https://github.com/Filco306/Topo-AE-Mapper}.} aim for a high separation of the spheres into different vertices as well as separate simplicial complexes for them, working under the presumption that vertices should contain points generated from only one sphere. We use the metric $m_{k, o, i}$ as seen in equation \ref{eq:eval-result}, where $V_{k, o, i}$ is the set of vertices in the mapper graph $\mathcal{G}_{k, o, i} := (V_{k, o, i}, E_{k, o, i})$ constructed using filtering function $k$, overlap $o$ and number of intervals $i$. $\mathcal{C} = \{c_1, ..., c_n\}$ is the set of different manifolds in the dataset to which each point can belong. In our initial experiments, $k \in \{\mbox{T-SNE}, \mbox{TAE}, \mbox{Kernel Density}, \mbox{SVD}, \mbox{PCA}, \mbox{Eccentricity}\}$ and $\mathcal{C} = \{c_1,..,c_{11}\}$, one for each hypersphere. For the trainable filtering functions, we train our filtering functions on a training set of data and visualize and build the graph on a test set. However, for TSNE, which is not made for predicting mapping points to its latent space, we use as many points (9000) to train, but only produce a graph on the test set points for a fair comparison.  We use a covering with overlappings rangings from 0.025 to 0.4 with a 0.025 step and an interval cardinality ranging between 5 and 45 with step of 5 and report the best and worst results in table \ref{table:res}. In other words, $o \in \mathcal{O} := \{0.025, 0.05, ..., 0.4\}$ and $i \in \mathcal{I} := \{5, 10, ..., 45\}$. Running over many choices can bring insights into a filtering function's robustness to choices of the aforementioned parameters.

\begin{equation} \label{eq:eval-result}
m_{k, o, i} = \frac{1}{|V_{k, o, i}|}\sum_{v \in V_{k, o, i}} \sum_{c \in \mathcal{C}} \mathds{1}(c \in v), \forall o \in \mathcal{O}, i \in \mathcal{I}
\end{equation}

The intuition behind our current metric $m_k$ is that fewer points belonging to different spheres in a vertex in the graph indicates a preservation of local as well as global topology. With a perfect separation of manifolds, every node will only contain points belonging to one (in our case, one sphere).
\section{Results}

Table \ref{table:res} display our initial results, indicating that TAEs are better at allocating the different spheres into separate parts than the other methods, capturing more of the topological properties of the data. Although t-SNE is close, it does not completely separate the spheres like the TAE. In figure \ref{fig:mappers} we provide visualizations of the best results from table \ref{table:res}, where we can see that the TAE forms a graph in which somewhat of a circle can be distinguished for the larger sphere, and simultaneously separate the smaller spheres to become separate. The kernel density and the eccentricity also achieves optimal results according to the metric, but for these optimal results, they do not capture the inherent shape as well as the TAE (bottom right).

\begin{figure}[h]
\centerline{\includegraphics[width=60mm]{best-res.png}}
\caption{Examples of Mapper graphs for each filtration with the best results from table \ref{table:res}. Node color indicate average label value of which sphere points in the node belong to - yellow is the large sphere, others are the smaller. }
\label{fig:mappers}
\end{figure}

\begin{table}[]
\centering
\caption{Best average number of spheres found in each node. Lower is better, 1 is optimal. }
\label{table:res}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\textbf{}        & \textbf{TAE}   & \textbf{PCA} & \textbf{T-SNE} & \textbf{SVD} & \textbf{Kernel Density} & \textbf{Eccentricity} \\ \hline
\textbf{Best}    & \textbf{1.000} & 1.1          & 1.026          & 1.119        & \textbf{1.0}            & \textbf{1.0}          \\ \hline
\textbf{Average} & \textbf{1.021} & 1.206        & 1.130          & 1.205        & 1.621                   & 1.137                 \\ \hline
\textbf{Worst}   & \textbf{1.115} & 1.311        & 1.259          & 1.311        & 2.25                    & 1.3                   \\ \hline
\end{tabular}
\end{table}

\begin{figure}[h]
\centerline{\includegraphics[width=100mm]{dists.png}}
\caption{Distribution of the performance over the grid search for the different filtering functions. As can be seen, the TAE is significantly more centered towards the optimal value of the metric used.  }
\label{fig:dists}
\end{figure}

By looking at the distribution of the metric $m_k$ in figure \ref{fig:dists} for the different filter functions, we can see that the TAE has a distribution significantly closer to the optimal result (being 1) compared to the other filtering functions. This shows that over a large range of parameters, the TAE outperforms the rest in generating representations separating the different spheres into different nodes.

\section{Discussion \& future work}

In this extended abstract, we have shown that topological autencoders are superior to other filtering functions in terms of representing high-dimensional data when combined with the Mapper algorithm on a set of manifolds. We plan to further extend the experiments to more simulated manifold datasets as well as the other datasets applied in \cite{moor2019topological} to confirm whether this holds generally. Furthermore, the metrics used to evaluate how a filtering function best captures the topology of the data as well as its robustness will be refined, which could elaborate the way to benchmark different filtering functions and their ability to better preserve the topology of higher-dimensional data. In essence, it will be beneficial to aim for a method preserving global and local topology while simultaneously gaining an empirical understanding of different filtering functions' strengths and shortcomings.

\section*{Acknowledgements}

This work was partially supported by Wallenberg Autonomous Systems Program (WASP).

\medskip