Deep Temporal Sigmoid Belief Networks \\ for Sequence Modeling
Zhe Gan
Chunyuan Li
Ricardo Henao
David Carlson and Lawrence Carin
Department of Electrical and Computer Engineering
Duke University
Durham
NC 27708
\{zhe.gan
chunyuan.li
r.henao
david.carlson
lcarin\
% Deep dynamic generative models are developed to learn sequential dependencies in time-series data. The multi-layered model is designed by constructing a hierarchy of temporal sigmoid belief networks (TSBNs), defined as a sequential stack of sigmoid belief networks (SBNs). Each SBN has a contextual hidden state, inherited from the previous SBNs in the sequence, and is used to regulate its hidden bias. Scalable learning and inference algorithms are derived by introducing a recognition model that yields fast sampling from the variational posterior. This recognition model is trained jointly with the generative model, by maximizing its variational lower bound on the log-likelihood. Experimental results on bouncing balls, polyphonic music, motion capture, and text streams show that the proposed approach achieves state-of-the-art predictive performance, and has the capacity to synthesize various sequences. %
DeepTSBN.tex
\documentclass{article} % For LaTeX2e
\usepackage{nips15submit_e,times}
\include{subtex/mlVecMat}

\usepackage{xr}
\externaldocument{DeepTSBNsupp}

\title{Deep Temporal Sigmoid Belief Networks \\ for Sequence Modeling}


\author{
Zhe Gan, Chunyuan Li, Ricardo Henao, David Carlson and Lawrence Carin \\
Department of Electrical and Computer Engineering\\
Duke University, Durham, NC 27708\\
\texttt{\{zhe.gan, chunyuan.li, r.henao, david.carlson, lcarin\}@duke.edu} \\
}

% The \author macro works with any number of authors. There are two commands
% used to separate the names and addresses of multiple authors: \And and \AND.
%
% Using \And between authors leaves it to \LaTeX{} to determine where to break
% the lines. Using \AND forces a linebreak at that point. So, if \LaTeX{}
% puts 3 of 4 authors names on the first line, and the last on the second
% line, try using \AND instead of \And before the third author name.

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\nipsfinalcopy % Uncomment for camera-ready version

\begin{document}


\maketitle

\input{subtex/abstract}
\input{subtex/introduction}
\input{subtex/model}
\input{subtex/inference}
\input{subtex/relatedwork}
\input{subtex/experiments}
\input{subtex/conclusion}
\input{subtex/acknowledgements}

\newpage
\bibliographystyle{unsrt}
\small{\bibliography{references}}

\input{subtex/SupplementaryMaterial}

\end{document}

SupplementaryMaterial.tex

\appendix
%
\section{Outline of the NVIL algorithm}
%
The outline of the NVIL algorithm for computing gradients are shown below (reproduced from \cite{mnih2014neural}). $C_{\lambdav} (\vv_t)$ represents the data-dependent baseline, and $\alpha=0.8$ throughout the experiments.
%
\vspace{-8mm}
\begin{center}
\begin{minipage}[t]{0.6\linewidth}
  \vspace{0pt}  
  \begin{algorithm}[H] %[th!]
  \caption{Compute gradient estimates for the model parameters and recognition parameters.}\label{alg:nvil}
  \small
  \begin{algorithmic}
  	\STATE $\Delta \thetav \leftarrow 0, \Delta \phiv \leftarrow 0, \Delta \lambdav \leftarrow 0$
  	\STATE $\Lcal \leftarrow 0$
  	\FOR{$t\leftarrow 1$ {\bfseries to} $T$}
  		\STATE $\hv_t \sim q_{\phiv} (\hv_t|\vv_t)$
  		\STATE $l_t \leftarrow \log p_{\thetav}(\vv_t,\hv_t) - \log q_{\phiv}(\hv_t|\vv_t)$
  		\STATE $\Lcal \leftarrow \Lcal + l_t$
  		\STATE $l_t \leftarrow l_t - C_{\lambdav} (\vv_t)$
  	\ENDFOR
  	\STATE $c_b \leftarrow \mbox{mean} (l_1,\ldots,l_T)$
  	\STATE $v_b \leftarrow \mbox{variance} (l_1,\ldots,l_T)$
  	\STATE $c \leftarrow \alpha c + (1-\alpha)c_b$
  	\STATE $v \leftarrow \alpha v + (1-\alpha)v_b$
  	\FOR{$t\leftarrow 1$ {\bfseries to} $T$}
  		\STATE $l_t \leftarrow \frac{l_t - c}{\max(1,\sqrt{v})}$
  		\STATE $\Delta \thetav \leftarrow \Delta \thetav + \nabla_{\thetav}\log p_{\thetav}(\vv_t,\hv_t)$
  		\STATE $\Delta \phiv \leftarrow \Delta \phiv + l_t\nabla_{\phiv}\log q_{\phiv}(\hv_t|\vv_t)$
  		\STATE $\Delta \lambdav \leftarrow \Delta \lambdav + l_t\nabla_{\lambdav} C_{\lambdav} (\vv_t)$
  	\ENDFOR
  \end{algorithmic}
  \end{algorithm}
\end{minipage}%
\end{center}

%
\section{Learning and Inference Details on TSBN}
%
For $t=1,\ldots,T$, consider $\vv_t \in \{0,1\}^M$, $\hv_t \in \{0,1\}^J$, the model parameters $\thetav$ are specified as $\Wmat_1 \in \R^{J \times J}$, $\Wmat_2 \in \R^{M \times J}$, $\Wmat_3 \in \R^{J \times M}$, $\Wmat_4 \in \R^{M \times M}$, $\bv \in \R^J$, and $\cv \in \R^M$. The generative model is expressed as
%
\begin{align}
p(h_{jt}=1|\hv_{t-1},\vv_{t-1}) &= \sigma(\wv_{1j}^\top \hv_{t-1} + \wv_{3j}^\top \vv_{t-1} + b_j)\,, \\
p(v_{mt} = 1 | \hv_t,\vv_{t-1}) &= \sigma(\wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1} + c_m)\,, \label{eqn:likelihood}
\end{align}
%
The recognition model is expressed as
%
\begin{align} \label{recogn}
q(h_{jt}=1|\hv_{t-1},\vv_t,\vv_{t-1}) = \sigma(\uv_{1j}^\top \hv_{t-1} + \uv_{2j}^\top \vv_t + \uv_{3j}^\top \vv_{t-1} + d_j )\,, 
\end{align}
%
where the recognition parameters are specified as $\Umat_1 \in \R^{J \times J}$, $\Umat_2 \in \R^{J \times M}$, $\Umat_3 \in \R^{J \times M}$, and $\dv \in \R^J$. $\hv_0$ and $\vv_0$, needed for $p(\hv_1)$, $p(\vv_1|\hv_1)$ and $q(\hv_1|\vv_1)$, are defined as zero vectors, for conciseness.

In order to implement the NVIL algorithm described in \cite{mnih2014neural}, we need to calculate the lower bound and also the gradients. Specifically, we have the variational lower bound $\Lcal = \sum_{t=1}^T \E_{q_{\phiv}(\hv|\vv)} [ l_t]$, where $l_t$ is expressed as 
%
\begin{align}
l_t &= \sum_{j=1}^J \left( \psi_{jt}^{(1)} h_{jt} - \log (1+\exp(\psi_{jt}^{(1)})) \right) + \sum_{m=1}^M \left( \psi_{mt}^{(2)} v_{mt} - \log (1+\exp(\psi_{mt}^{(2)})) \right) \\
&- \left[ \sum_{j=1}^J \left( \psi_{jt}^{(3)} h_{jt} - \log (1+\exp(\psi_{jt}^{(3)})) \right) \right] \,, \nonumber
\end{align}
%
and we have defined
\begin{align} 
\psi_{jt}^{(1)} &= \wv_{1j}^\top \hv_{t-1} + \wv_{3j}^\top \vv_{t-1} + b_j \,, \\
\psi_{mt}^{(2)} &= \wv_{2m}^\top \hv_{t} + \wv_{4m}^\top \vv_{t-1} + c_m \,, \\
\psi_{jt}^{(3)} &= \uv_{1j}^\top \hv_{t-1} + \uv_{2j}^\top \vv_{t} + \uv_{3j}^\top \vv_{t-1} + d_j \,. 
\end{align}
%
By further defining 
%
\begin{align}
\chi_{jt}^{(1)} = h_{jt} - \sigma(\psi_{jt}^{(1)}), \,\,\,\, \chi_{mt}^{(2)} = v_{mt} - \sigma(\psi_{mt}^{(2)}), \,\,\,\,
\chi_{jt}^{(3)} = h_{jt} - \sigma(\psi_{jt}^{(3)}),
\end{align}
%
The gradients for the model parameters $\thetav$ are expressed as
%
\begin{align}
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{1jj^{\prime}}} &=  \chi_{jt}^{(1)} h_{j^{\prime}t-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{3jm}} &=  \chi_{jt}^{(1)} v_{mt-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial b_j} &=  \chi_{jt}^{(1)} \,, \label{eqn:begin} \\ 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{2mj}} &=  \chi_{mt}^{(2)} h_{jt}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{4mm^{\prime}}} &=  \chi_{mt}^{(2)} v_{m^{\prime}t-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial c_m} &=  \chi_{mt}^{(2)} \,. 
\end{align}
%
The gradients for the recognition parameters $\phiv$ are expressed as
%
\begin{align}
\frac{\partial \log q_{\phiv}(\hv_t|\vv_t)}{\partial u_{1jj^{\prime}}} &=  \chi_{jt}^{(3)} h_{j^{\prime}t-1}, & 
\frac{\partial \log q_{\phiv}(\hv_t|\vv_t)}{\partial u_{2jm}} &=  \chi_{jt}^{(3)} v_{mt}\,, \\ 
\frac{\partial \log q_{\phiv}(\hv_t|\vv_t)}{\partial u_{3jm}} &=  \chi_{jt}^{(3)} v_{mt-1}, & 
\frac{\partial \log q_{\phiv}(\hv_t|\vv_t)}{\partial d_j} &=  \chi_{jt}^{(3)} \,. \label{eqn:end}
\end{align}
%
\subsection{Modeling Real-valued Data}
%
When modeling real-valued data, we substitute (\ref{eqn:likelihood}) with $p(\vv_t|\hv_t,\vv_{t-1}) = \Ncal (\muv_t,\mbox{diag}(\sigmav_t^2))$, where 
%
\begin{align}
\mu_{mt} = \wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1}+c_m, \,\,\,\
\log \sigma_{mt} = (\wv_{2m}^{\prime})^\top \hv_t + (\wv_{4m}^{\prime})^\top \vv_{t-1}+c_m^{\prime},
\end{align}
%
and we have $\Wmat_2^{\prime} \in \R^{M \times J}$ and $\Wmat_4^{\prime} \in \R^{M \times M}$. The recognition model remains the same as in (\ref{recogn}). Let $\tau_{mt} = \log \sigma_{mt}$, we obtain 
%
\begin{align}
l_t &= \sum_{j=1}^J \left( \psi_{jt}^{(1)} h_{jt} - \log (1+\exp(\psi_{jt}^{(1)})) \right) - \sum_{m=1}^M \left( \frac{1}{2} \log 2\pi + \tau_{mt} + \frac{(v_{mt}-\mu_{mt})^2}{2e^{2\tau_{mt}}} \right) \\
&- \left[ \sum_{j=1}^J \left( \psi_{jt}^{(3)} h_{jt} - \log (1+\exp(\psi_{jt}^{(3)})) \right) \right] \,. \nonumber
\end{align}
%
All the gradient calculation remains the same as (\ref{eqn:begin})-(\ref{eqn:end}), except the following.
%
\begin{align}
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{2mj}} &=  \chi_{mt}^{(4)} h_{jt}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{4mm^{\prime}}} &=  \chi_{mt}^{(4)} v_{m^{\prime}t-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial c_m} &=  \chi_{mt}^{(4)} \,, \\ 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{2mj}^{\prime}} &=  \chi_{mt}^{(5)} h_{jt}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{4mm^{\prime}}^{\prime}} &=  \chi_{mt}^{(5)} v_{m^{\prime}t-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial c_m^{\prime}} &=  \chi_{mt}^{(5)} \,,
\end{align}
%
where we have defined
%
\begin{align}
\chi_{mt}^{(4)} = \frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial \mu_{mt}} = \frac{v_{mt}-\mu_{mt}}{e^{2\tau_{mt}}}, \,\,\,\ \chi_{mt}^{(5)} = \frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial \tau_{mt}} = \frac{(v_{mt}-\mu_{mt})^2}{e^{2\tau_{mt}}} -1 \,.
\end{align}
%
\subsection{Modeling Count Data}
%
We also introduce an approach for modeling time-series data with count observations, by replacing (\ref{eqn:likelihood}) with $p(\vv_t|\hv_t,\vv_{t-1}) = \prod_{m=1}^M y_{mt}^{v_{mt}}$, where
%
\begin{align}
y_{mt} = \frac{\exp(\wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1}+c_m)}{\sum_{m^{\prime}=1}^M \exp(\wv_{2m^{\prime}}^\top \hv_t + \wv_{4m^{\prime}}^\top \vv_{t-1}+c_{m^{\prime}})} \,.
\end{align}
%
The recognition model still remains the same as in (\ref{recogn}). The $l_t$ now is expressed as
%
\begin{align}
l_t &= \sum_{j=1}^J \left( \psi_{jt}^{(1)} h_{jt} - \log (1+\exp(\psi_{jt}^{(1)})) \right) + \sum_{m=1}^M \left( \psi_{mt}^{(2)} v_{mt} - v_{mt} \log \sum_{m^{\prime}=1}^M \exp(\psi_{mt}^{(2)}) \right) \\
&- \left[ \sum_{j=1}^J \left( \psi_{jt}^{(3)} h_{jt} - \log (1+\exp(\psi_{jt}^{(3)})) \right) \right] \,. \nonumber
\end{align}
%
All the gradient calculations remain the same as (\ref{eqn:begin})-(\ref{eqn:end}), except the following
%
\begin{align}
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{2mj}} &=  \chi_{mt}^{(6)} h_{jt}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial w_{4mm^{\prime}}} &=  \chi_{mt}^{(6)} v_{m^{\prime}t-1}, & 
\frac{\partial \log p_{\thetav}(\vv_t,\hv_t)}{\partial c_m} &=  \chi_{mt}^{(6)} \,.
\end{align}
%
where we have defined $\chi_{mt}^{(6)} = v_{mt} - y_{mt}\sum_{m^{\prime}=1}^M v_{m^{\prime}t}$.
%
\section{Learning and Inference Details on Deep TSBN}
%
For the ease of notation, we consider a two-hidden-layer deep TSBN here, which can be readily extended to a deep model with any depth. For $t=1,\ldots,T$, we consider the observation as $\vv_t \in \{0,1\}^M$. The top hidden layer is denoted as $\zv_t \in \{0,1\}^J$.
%
\input{subtex/gmodel_supp}
%
\subsection{Using stochastic hidden layer}
%
Denote the first stochastic hidden layer as $\hv_t \in \{0,1\}^K$. The generative model is expressed as
%
\begin{align}
p(z_{jt} = 1) &= \sigma(\wv_{1j}^\top \zv_{t-1} + \wv_{3j}^\top \hv_{t-1} + b_{1j}) \,, \\
p(h_{kt} = 1) &= \sigma(\wv_{2k}^\top \zv_t + \wv_{4k}^\top \hv_{t-1} + \wv_{6k}^\top \vv_{t-1}+b_{2k}) \,, \\
p(v_{mt} = 1) &= \sigma(\wv_{5m}^\top \hv_t + \wv_{7m}^\top \vv_{t-1}+b_{3m}) \,,
\end{align}
%
where we have defined $\Wmat_1 \in \R^{J \times J}$, $\Wmat_2 \in \R^{K \times J}$, $\Wmat_3 \in \R^{J \times K}$, $\Wmat_4 \in \R^{K \times K}$, $\Wmat_5 \in \R^{M \times K}$, $\Wmat_6 \in \R^{K \times M}$, and $\Wmat_7 \in \R^{M \times M}$. The bias terms are $\bv_1 \in \R^{J \times 1}$, $\bv_2 \in \R^{K \times 1}$ and $\bv_3 \in \R^{M \times 1}$.
%
The corresponding recognition model is expressed as
%
\begin{align}
q(h_{kt} = 1) &= \sigma (\uv_{5k}^\top \vv_t + \uv_{4k}^\top \hv_{t-1} + \uv_{6k}^\top \vv_{t-1} + c_{2k}) \\
q(\zv_{jt} = 1) &= \sigma (\uv_{2j}^\top \hv_t + \uv_{1j}^\top \zv_{t-1} + \uv_{3j}^\top \hv_{t-1} + c_{1j})
\end{align}
%
where the recognition parameters are specified as $\Umat_1 \in \R^{J \times J}$, $\Umat_2 \in \R^{J \times K}$, $\Umat_3 \in \R^{J \times K}$, $\Umat_4 \in \R^{K \times K}$, $\Umat_5 \in \R^{K \times M}$ and $\Umat_6 \in \R^{K \times M}$. The bias terms are $\cv_1 \in \R^{J \times 1}$ and $\cv_2 \in \R^{K \times 1}$. 
%
Now, $l_t$ is expressed as
%
\begin{align}
l_t &= \sum_{j=1}^J \left( \psi_{jt}^{(1)} z_{jt} - \log (1+\exp(\psi_{jt}^{(1)})) \right) + \sum_{k=1}^K \left( \psi_{kt}^{(2)} h_{kt} - \log (1+\exp(\psi_{kt}^{(2)})) \right) \nonumber \\
&+ \sum_{m=1}^M \left( \psi_{mt}^{(3)} v_{mt} - \log (1+\exp(\psi_{mt}^{(3)})) \right) \\
&- \left[\sum_{k=1}^K \left( \psi_{kt}^{(4)} h_{kt} - \log (1+\exp(\psi_{kt}^{(4)})) \right) +  \sum_{j=1}^J \left( \psi_{jt}^{(5)} z_{jt} - \log (1+\exp(\psi_{jt}^{(5)})) \right) \right] \,, \nonumber
\end{align}
%
and we have defined
\begin{align} 
\psi_{jt}^{(1)} &= \wv_{1j}^\top \zv_{t-1} + \wv_{3j}^\top \hv_{t-1} + b_{1j} \,, \\
\psi_{kt}^{(2)} &= \wv_{2k}^\top \zv_t + \wv_{4k}^\top \hv_{t-1} + \wv_{6k}^\top \vv_{t-1}+b_{2k} \,, \\
\psi_{mt}^{(3)} &= \wv_{5m}^\top \hv_t + \wv_{7m}^\top \vv_{t-1}+b_{3m} \,, \\
\psi_{kt}^{(4)} &= \uv_{5k}^\top \vv_t + \uv_{4k}^\top \hv_{t-1} + \uv_{6k}^\top \vv_{t-1} + c_{2k} \,, \\
\psi_{jt}^{(3)} &= \uv_{2j}^\top \hv_t + \uv_{1j}^\top \zv_{t-1} + \uv_{3j}^\top \hv_{t-1} + c_{1j} \,. 
\end{align}
%
All the gradients can be calculated readily as in (\ref{eqn:begin})-(\ref{eqn:end}).
%
\subsection{Using deterministic hidden layer}

\begin{table}[t!]
\begin{center}
\begin{minipage}{0.75\linewidth}
%\vskip -0.3in
\caption{\small{Average prediction error and the average negative log-likelihood per frame for the bouncing balls dataset. $(\diamond)$ taken from \cite{mittelman2014structured}.}}\label{Table:bouncing_ball}
\vskip -0.1in
\centering
%\raggedleft
\small
\begin{sc}
\begin{adjustbox}{scale=1.0,tabular=llcrc,center}
\hline
% \abovespace \belowspace
Model & Dim & Order & Pred. Err. & Neg. Log. Like. \\
\hline
DTSBN-s  & 100-100 & 2 & {\bf 2.79} $\pm$ 0.39 & {\bf 69.29} $\pm$ 1.52 \\
DTSBN-d  & 100-100 & 2 & 2.99 $\pm$ 0.42 & 70.47 $\pm$ 1.52\\
DTSBN-s  & 100-100 & 1 & 10.39 $\pm$ 0.38 & 78.63 $\pm$ 0.92 \\
\hline
TSBN  & 100 & 4 & 3.07 $\pm$ 0.40  & 70.41 $\pm$ 1.55 \\
TSBN  & 100 & 2 & 4.00 $\pm$ 0.45  & 73.32 $\pm$ 1.75 \\
TSBN  & 100 & 1 & 9.48 $\pm$ 0.38  & 77.71 $\pm$ 0.83 \\
HMSBN & 100 & 1 & 23.94 $\pm$ 0.41 & 86.27 $\pm$ 0.80 \\
\hline
AR    & 0  & 2 & 3.63  $\pm$ 0.42 & 73.80 $\pm$ 1.46 \\
AR    & 0  & 1 & 11.01 $\pm$ 0.24 & 93.61 $\pm$ 0.67 \\
\hline 
RTRBM$^\diamond$  & 3750 & 1 & 3.88 $\pm$ 0.33 & $-$\\
SRTRBM$^\diamond$  & 3750 & 1 & 3.31 $\pm$ 0.33 & $-$\\
\hline
\end{adjustbox}
\end{sc}
\end{minipage}
\end{center}
\vspace{-5mm}
\end{table}
%
For the generative model, denote the deterministic hidden layer as $\hv_t^g \in \R^K$. For the recognition model, denote the deterministic hidden layer as $\hv_t^r \in \R^K$. $\Wmat_3$ and $\Umat_3$ are set to be zero matrices for the ease of gradient calculation. The generative model is expressed as
%
\begin{align}
p(z_{jt} = 1) &= \sigma(\wv_{1j}^\top \zv_{t-1} + b_{1j}) \,, \\
h_{kt}^g &= f(\wv_{2k}^\top \zv_t + \wv_{4k}^\top \hv_{t-1}^g + \wv_{6k}^\top \vv_{t-1}+b_{2k}) \,, \\
p(v_{mt} = 1) &= \sigma(\wv_{5m}^\top \hv_t + \wv_{7m}^\top \vv_{t-1}+b_{3m}) \,,
\end{align}
%
The corresponding recognition model is expressed as
%
\begin{align}
h_{kt}^r &= f (\uv_{5k}^\top \vv_t + \uv_{4k}^\top \hv_{t-1} + \uv_{6k}^\top \vv_{t-1} + c_{2k}) \\
q(\zv_{jt} = 1) &= \sigma (\uv_{2j}^\top \hv_t + \uv_{1j}^\top \zv_{t-1} + c_{1j})
\end{align}
%
Hence, $\Lcal = \sum_{t=1}^T \E_{q_{\phiv}(\hv|\vv)} [ l_t]$, and $l_t$ is expressed as
%
\begin{align}
l_t &= \sum_{j=1}^J \left( \psi_{jt}^{(1)} z_{jt} - \log (1+\exp(\psi_{jt}^{(1)})) \right) 
+ \sum_{m=1}^M \left( \psi_{mt}^{(2)} v_{mt} - \log (1+\exp(\psi_{mt}^{(2)})) \right) \\
&- \left[\sum_{j=1}^J \left( \psi_{jt}^{(3)} z_{jt} - \log (1+\exp(\psi_{jt}^{(3)})) \right) \right] \,, \nonumber
\end{align}
%
and we have defined
\begin{align} 
\psi_{jt}^{(1)} &= \wv_{1j}^\top \zv_{t-1} + b_{1j} \,, \\
\psi_{mt}^{(2)} &= \wv_{5m}^\top \hv_t^g + \wv_{7m}^\top \vv_{t-1}+b_{3m} \,, \\
\psi_{jt}^{(3)} &= \uv_{2j}^\top \hv_t^r + \uv_{1j}^\top \zv_{t-1} + c_{1j} \,. 
\end{align}
%
The gradients \emph{w.r.t.} $\Wmat_1,\Wmat_5,\Wmat_7, \Umat_1$ and $\Umat_2$ can be calculated easily. In order to calculate the gradients \emph{w.r.t.} $\Wmat_2, \Wmat_4, \Wmat_6, \Umat_4, \Umat_5$ and $\Umat_6$, we need to obtain $\frac{\partial \Lcal}{\partial h_{kt}^g}$ and $\frac{\partial \Lcal}{\partial h_{kt}^r}$, which can be calculated recursively via the back-propagation through time algorithm. Specifically, $\frac{\partial \Lcal}{\partial h_{kt}^g} = \frac{\partial \Qcal_1}{\partial h_{kt}^g}$ and we have defined
%
\begin{align}
\Qcal_1 = \sum_{t=1}^T \sum_{m=1}^M \left( \psi_{mt}^{(2)} v_{mt} - \log (1+\exp(\psi_{mt}^{(2)})) \right) \,.
\end{align}
%
We observe that $\Qcal_1$ can be computed recursively using
%
\begin{align}
\Qcal_{t} &= \sum_{\tau=t}^T \sum_{m=1}^M \left( \psi_{m\tau}^{(2)} v_{m\tau} - \log (1+\exp(\psi_{m\tau}^{(2)}))  \right) \\
&= \Qcal_{t+1} + \sum_{m=1}^M \left( \psi_{mt}^{(2)} v_{mt} - \log (1+\exp(\psi_{mt}^{(2)}))  \right) \,,
\end{align}
%
where $\Qcal_{T+1} = 0$. Using the chain rule, we have
%
\begin{align}
\frac{\partial \Qcal_t}{\partial h_{kt}^g} &= \sum_{k^{\prime}} \frac{\partial \Qcal_{t+1}}{\partial h_{k^{\prime}t+1}^g} \cdot \frac{\partial h_{k^{\prime}t+1}^g}{\partial h_{kt}^g} + \sum_{m=1}^M w_{5mk} (v_{mt}-\sigma(\psi_{mt}^{(2)})) \\
&= \sum_{k^{\prime}} \frac{\partial \Qcal_{t+1}}{\partial h_{k^{\prime}t+1}^g} \cdot f^{\prime} (\psi_{k^{\prime}t+1}^{(4)}) w_{4k^{\prime}k} + \sum_{m=1}^M w_{5mk} (v_{mt}-\sigma(\psi_{mt}^{(2)})) \,,
\end{align}
%
where we have defined 
%
\begin{align} 
\psi_{kt}^{(4)} &= \wv_{2k}^\top \zv_t + \wv_{4k}^\top \hv_{t-1}^g + \wv_{6k}^\top \vv_{t-1}^g + b_{2k} \,, 
\end{align}
%
and 
\begin{align}
\frac{\partial \Qcal_T}{\partial h_{kt}^g} = \sum_{m=1}^M w_{5mk} (v_{mT}-\sigma(\psi_{mT}^{(2)})) \,.
\end{align}
$\frac{\partial \Lcal}{\partial h_{kt}^r}$ can be calculated similarly.

\section{Additional Results}
%
\subsection{Generated Data}
%
The generated, synthetic motion capture data, and polyphonic music data can be downloaded from \url{https://drive.google.com/drive/u/0/folders/0B1HR6m3IZSO_SWt0aS1oYmlneDQ}.
%
\subsection{Bouncing balls dataset}
%
Additional experimental results are shown in Table \ref{Table:bouncing_ball}. AR represents an auto-regressive Markov model without latent variables.
%
\subsection{MIT motion capture dataset}
%
We randomly select 10\% of the dataset as the test set. Quantitative results are shown in Table \ref{Table:mocap}.

\begin{table}[t!]
\begin{center}
%
\begin{minipage}{0.45\linewidth}
%\vskip -0.3in
\caption{\small{Average prediction error obtained for the MIT motion capture dataset.}}\label{Table:mocap}
\vskip -0.1in
\centering
%\raggedleft
\small
\begin{sc}
\begin{adjustbox}{scale=1.0,tabular=lr,center}
\hline
% \abovespace \belowspace
Model & Pred. Err. \\
\hline
DTSBN-s   & {\bf 3.71} $\pm$ 0.03  \\
DTSBN-d   & 4.19 $\pm$ 0.01  \\
TSBN   & 3.86 $\pm$ 0.02  \\
HMSBN  & 17.49 $\pm$ 0.20  \\
\hline
\end{adjustbox}
\end{sc}
\end{minipage}
\end{center}
\vspace{-5mm}
\end{table}


abstract.tex
%!TEX root = ../DeepTSBN.tex

\begin{abstract}
%
Deep dynamic generative models are developed to learn sequential dependencies in time-series data. The multi-layered model is designed by constructing a hierarchy of temporal sigmoid belief networks (TSBNs), defined as a sequential stack of sigmoid belief networks (SBNs). Each SBN has a contextual hidden state, inherited from the previous SBNs in the sequence, and is used to regulate its hidden bias. Scalable learning and inference algorithms are derived by introducing a recognition model that yields fast sampling from the variational posterior. This recognition model is trained jointly with the generative model, by maximizing its variational lower bound on the log-likelihood. Experimental results on bouncing balls, polyphonic music, motion capture, and text streams show that the proposed approach achieves state-of-the-art predictive performance, and has the capacity to synthesize various sequences.
%
\end{abstract}

acknowledgements.tex

%\section*{Acknowledgements}
\vspace{-0.1in}

\paragraph{Acknowledgements} This research was supported in part by ARO, DARPA, DOE, NGA and ONR. 


conclusion.tex
\vspace{-0.1in} 
\section{Conclusion} \label{sec:conclusion}
\vspace{-0.1in} 
%
We have presented the Deep Temporal Sigmoid Belief Networks, an extension of SBN, that models the temporal dependencies in high-dimensional sequences. To allow for scalable inference and learning, an efficient variational optimization algorithm is developed. Experimental results on several datasets show that the proposed approach obtains superior predictive performance, and synthesizes interesting sequences. 

In this work, we have investigated the modeling of different types of data individually. One interesting future work is to combine them into a unified framework for dynamic multi-modality learning. Furthermore, we can use high-order optimization methods to speed up inference \cite{fan2015fast}. 

experiments.tex
%!TEX root = ../DeepTSBN.tex

\section{Experiments} \label{sec:exp}
%
\begin{figure}[t!]
%\vskip -0.1in
\centering
\begin{minipage}{0.31\linewidth}
	\includegraphics[width=0.9\linewidth]{figure/bouncing_balls_dictionary}
\end{minipage}
%
\begin{minipage}{0.38\linewidth}
	\includegraphics[width=0.9\linewidth]{figure/music_generation}
\end{minipage}
%
\begin{minipage}{0.27\linewidth}
	\includegraphics[width=\linewidth]{figure/topic29} \\
	\includegraphics[width=\linewidth]{figure/topic30} \\
	\includegraphics[width=\linewidth]{figure/topic130}
\end{minipage}
%
%\vskip -0.1in
\caption{(Left) Dictionaries learned using the HMSBN for the videos of bouncing balls. (Middle) Samples generated from the HMSBN trained on the polyphonic music. Each column is a sample vector of notes. (Right) Time evolving from 1790 to 2014 for three selected topics learned from the STU dataset. Plotted values represent normalized probabilities that the topic appears in a given year. Best viewed electronically.}
\label{fig:performance}
\vskip -0.1in
\end{figure}
%
We present experimental results on four publicly available datasets: the bouncing balls \cite{sutskever2009recurrent}, polyphonic music \cite{boulanger2012modeling}, motion capture \cite{taylor2006modeling} and state-of-the-Union \cite{NIPS2014_RLike_Shao}. To assess the performance of the TSBN model, we show sequences generated from the model, and report the average log-probability that the model assigns to a test sequence, and the average squared one-step-ahead prediction error per frame. Code is available at \url{https://github.com/zhegan27/TSBN_code_NIPS2015}.

The TSBN model with $\Wmat_3 = \boldsymbol{0}$ and $\Wmat_4 = \boldsymbol{0}$ is denoted Hidden Markov SBN (HMSBN), the deep TSBN with stochastic hidden layer is denoted DTSBN-S, and the deep TSBN with deterministic hidden layer is denoted DTSBN-D.

Model parameters were initialized by sampling randomly from $\Ncal (\boldsymbol{0},0.001^2 \Imat)$, except for the bias parameters, that were initialized as 0. The TSBN model is trained using a variant of RMSprop \cite{graves2013generating}, with momentum of 0.9, and a constant learning rate of $10^{-4}$. The decay over the root mean squared gradients is set to 0.95. The maximum number of iterations we use is $10^5$. The gradient estimates were computed using a single sample from the recognition model. The only regularization we used was a weight decay of $10^{-4}$. The data-dependent baseline was implemented by using a neural network with a single hidden layer with 100 tanh units. 

For the prediction of $\vv_t$ given $\vv_{1:t-1}$, we (\emph{i}) first obtain a sample from $q_{\phiv}(\hv_{1:t-1}|\vv_{1:t-1})$; (\emph{ii}) calculate the conditional posterior $p_{\thetav}(\hv_t|\hv_{1:t-1},\vv_{1:t-1})$ of the current hidden state ; (\emph{iii}) make a prediction for $\vv_t$ using $p_{\thetav}(\vv_t|\hv_{1:t},\vv_{1:t-1})$. On the other hand, synthesizing samples is conceptually simper. Sequences can be readily generated from the model using ancestral sampling.
%
\subsection{Bouncing balls dataset}
%
We conducted the first experiment on synthetic videos of 3 bouncing balls, where pixels are binary valued. We followed the procedure in \cite{sutskever2009recurrent}, and generated 4000 videos for training, and another 200 videos for testing. Each video is of length 100 and of resolution $30 \times 30$. 

The dictionaries learned using the HMSBN are shown in Figure \ref{fig:performance} (Left). Compared with previous work \cite{sutskever2009recurrent,boulanger2012modeling}, our learned bases are more spatially localized. In Table \ref{Table:bouncing_ball},  we compare the average squared prediction error per frame over the 200 test videos, with recurrent temporal RBM (RTRBM) and structured RTRBM (SRTRBM). As can be seen, our approach achieves better performance compared with the baselines in the literature. Furthermore, we observe that a high-order TSBN reduces the prediction error significantly, compared with an order-one TSBN. This is due to the fact that by using a high-order TSBN, more information about the past is conveyed. We also examine the advantage of employing deep models. Using stochastic, or deterministic hidden layer improves performances. More results, including log-likelihoods, are provided in Supplementary Section D.
%
\begin{table} [t!]
\begin{center}
\begin{minipage}{0.45\linewidth}
%\vskip -0.3in
\caption{\small{Average prediction error for the bouncing balls dataset. $(\diamond)$ taken from \cite{mittelman2014structured}.}}\label{Table:bouncing_ball}
\vskip -0.1in
\centering
%\raggedleft
\small
\begin{sc}
\begin{adjustbox}{scale=0.9,tabular=llcr,center}
\hline
% \abovespace \belowspace
Model & Dim & Order & Pred. Err. \\
\hline
DTSBN-s  & 100-100 & 2 & {\bf 2.79} $\pm$ 0.39 \\
DTSBN-d  & 100-100 & 2 & 2.99 $\pm$ 0.42 \\
TSBN  & 100 & 4 & 3.07 $\pm$ 0.40 \\
TSBN  & 100 & 1 & 9.48 $\pm$ 0.38 \\
\hline 
RTRBM$^\diamond$  & 3750 & 1 & 3.88 $\pm$ 0.33 \\
SRTRBM$^\diamond$  & 3750 & 1 & 3.31 $\pm$ 0.33 \\
\hline
\end{adjustbox}
\end{sc}
\end{minipage}
%
\hspace{0.3cm}
%
\begin{minipage}{0.45\linewidth}
%\vskip -0.3in
\caption{\small{Average prediction error obtained for the motion capture dataset. $(\diamond)$ taken from \cite{mittelman2014structured}.}}\label{Table:mocap}
\vskip -0.1in
\centering
%\raggedleft
\small
\begin{sc}
\begin{adjustbox}{scale=0.9,tabular=llrr,center}
\hline
% \abovespace \belowspace
Model & Walking & Running \\
\hline
DTSBN-s   & {\bf 4.40} $\pm$ 0.28 & {\bf 2.56} $\pm$ 0.40 \\
DTSBN-d   & 4.62 $\pm$ 0.01 & 2.84 $\pm$ 0.01 \\
TSBN   & 5.12 $\pm$ 0.50 & 4.85 $\pm$ 1.26 \\
HMSBN  & 10.77 $\pm$ 1.15 & 7.39 $\pm$ 0.47 \\
\hline 
ss-SRTRBM$^\diamond$   & 8.13 $\pm$ 0.06 & 5.88 $\pm$ 0.05 \\
g-RTRBM$^\diamond$   & 14.41 $\pm$ 0.38 & 10.91 $\pm$ 0.27 \\
\hline
\end{adjustbox}
\end{sc}
\end{minipage}
\end{center}
\vspace{-5mm}
\end{table}
%
\vspace{-0.1in} 
\subsection{Motion capture dataset}
\vspace{-0.05in} 
%
In this experiment, we used the CMU motion capture dataset, that consists of measured joint angles for different motion types. We used the 33 running and walking sequences of subject 35 (23 walking sequences and 10 running sequences). We followed the preprocessing procedure of \cite{mittelman2014structured}, after which we were left with 58 joint angles. We partitioned the 33 sequences into training and testing set: the first of which had 31 sequences, and the second had 2 sequences (one walking and another running). We averaged the prediction error over 100 trials, as reported in Table \ref{Table:mocap}. The TSBN we implemented is of size 100 in each hidden layer and order 1. It can be seen that the TSBN-based models improves over the Gaussian (G-)RTRBM and the spike-slab (SS-)SRTRBM significantly. 
%
%\footnote{\tiny{\url{http://mocap.cs.cmu.edu/}}}
%
\begin{figure}[h]
%\vskip -0.1in
\centering
\includegraphics[width=0.32\linewidth]{figure/mocap_walk}
\includegraphics[width=0.32\linewidth]{figure/mocap_run_run_walk}
\includegraphics[width=0.32\linewidth]{figure/mocap_run_walk}
%\vskip -0.1in
\caption{Motion trajectories generated from the HMSBN trained on the motion capture dataset. (Left) Walking. (Middle) Running-running-walking. (Right) Running-walking. }
\label{fig:mocap}
\vskip -0.1in
\end{figure}

Another popular motion capture dataset is the MIT dataset\footnote{\small{Quantitative results on the MIT dataset are provided in Supplementary Section D.}}. To further demonstrate the directed, generative nature of our model, we give our trained HMSBN model different initializations, and show generated, synthetic data and the transitions between different motion styles in Figure \ref{fig:mocap}. These generated data are readily produced from the model and demonstrate realistic behavior. The smooth trajectories are walking movements, while the vibrating ones are running. Corresponding video files (AVI) are provided as mocap 1, 2 and 3 in the Supplementary Material.
%
%Depending on initialization, we can synthesize various motions, and learn the transitions between them. 
%
%\footnote{\tiny{\url{http://people.csail.mit.edu/ehsu/work/sig05stf/}}} 
%
\vspace{-0.1in} 
\subsection{Polyphonic music dataset}
\vspace{-0.05in} 
%
The third experiment is based on four different polyphonic music sequences of piano \cite{boulanger2012modeling}, \emph{i.e.}, Piano-midi.de (Piano), Nottingham (Nott), MuseData (Muse) and JSB chorales (JSB). Each of these datasets are represented as a collection of 88-dimensional binary sequences, that span the whole range of piano from A0 to C8. 

The samples generated from the trained HMSBN model are shown in Figure \ref{fig:performance} (Middle). As can be seen, different styles of polyphonic music are synthesized. The corresponding MIDI files are provided as music 1 and 2 in the Supplementary Material. Our model has the ability to learn basic harmony rules and local temporal coherence. However, long-term structure and musical melody remain elusive. The variational lower bound, along with the estimated log-likelihood in \cite{boulanger2012modeling}, are presented in Table \ref{Table:music}. The TSBN we implemented is of size $100$ and order $1$. Empirically, adding layers did not improve performance on this dataset, hence no such results are reported. The results of RNN-NADE and RTRBM \cite{boulanger2012modeling} were obtained by only $100$ runs of the annealed importance sampling, which has the potential to overestimate the true log-likelihood. Our variational lower bound provides a more conservative estimate. Though, our performance is still better than that of RNN.
%
\subsection{State of the Union dataset}
%\vspace{-0.05in}
%
\begin{table} [t!]
\begin{center}
\begin{minipage}{0.45\linewidth}
%\vskip -0.3in
\caption{\small{Test log-likelihood for the polyphonic music dataset. $(\diamond)$ taken from \cite{boulanger2012modeling}.}}\label{Table:music}
\vskip -0.1in
%\centering
\raggedright
\small
\begin{sc}
\begin{adjustbox}{scale=0.85,tabular=lcccc,center}
\hline
Model & Piano. & Nott. & Muse. & JSB. \\
\hline
TSBN  & -7.98 & -3.67 & -6.81 & -7.48 \\
\hline 
RNN-NADE$^\diamond$  & {\bf -7.05} & {\bf -2.31} & {\bf -5.60} & {\bf -5.56} \\
RTRBM$^\diamond$  & -7.36 & -2.62 & -6.35 & -6.35 \\
RNN$^\diamond$ & -8.37 & -4.46 & -8.13 & -8.71 \\
\hline
\end{adjustbox}
\end{sc}
\end{minipage}
%
\hspace{0.4cm}
%
\begin{minipage}{0.45\linewidth}
	%\vskip -0.3in
	\caption{\small{Average prediction precision for STU. $(\diamond)$ taken from \cite{GP_DFA_AISTATS2015}.}}\label{Table:STU}
	\vskip -0.1in
	%\centering
	\raggedleft
	\small
	\begin{sc}
		\begin{adjustbox}{scale=0.8,tabular=llcc,center}
			\hline
			% \abovespace \belowspace
			Model & Dim & MP & PP \\
			\hline
			HMSBN  & 25 & ${\bf 0.327}\pm0.002 $ & $0.353\pm0.070$  \\
			DHMSBN-s  & 25-25 &$0.299\pm0.001$ & ${\bf 0.378}\pm0.006$  \\
			\hline 
			GP-DPFA $^\diamond$  & 100 &   $0.223\pm0.001$     & $0.189\pm0.003$ \\
			DRFM$^\diamond$  & 25  &$0.217\pm0.003$  & $0.177\pm0.010$  \\
			\hline
		\end{adjustbox}
	\end{sc}
\end{minipage}
\end{center}
\vspace{-5mm}
\end{table}

The State of the Union (STU) dataset contains the transcripts of $T = 225$ US State of the Union addresses, from 1790 to 2014. Two tasks are considered, \emph{i.e.}, prediction and dynamic topic modeling. 
\vspace{-0.2in} 
%
\paragraph{Prediction}
The prediction task is concerned with estimating the held-out words. We employ the setup in \cite{GP_DFA_AISTATS2015}. After removing stop words and terms that occur fewer than 7 times in one document or less than 20 times overall, there are 2375 unique words. The entire data of the last year is held-out.  For the documents in the previous years, we randomly partition the words of each document into 80\%/20\% split. The model is trained on the 80\% portion, and the remaining 20\% held-out words  are used to test the prediction at each year. The words in both held-out sets are ranked according to the probability estimated from (\ref{softmax}).

To evaluate the prediction performance, we calculate the precision $@$top-$M$as in~\cite{GP_DFA_AISTATS2015}, which is given by the fraction of the top-$M$ words, predicted by the model, that matches the true ranking of the word counts.  $M = 50$ is used.  Two recent works are compared, GP-DPFA \cite{GP_DFA_AISTATS2015} and DRFM \cite{NIPS2014_RLike_Shao}. The results are summarized in Table~\ref{Table:STU}. Our model is of order 1. The column MP denotes the mean precision over all the years that appear in the training set. The column PP denotes the predictive precision for the final year. Our model achieves significant improvements in both scenarios.
\vspace{-0.1in}
% 
\paragraph{Dynamic Topic Modeling }
The setup described in \cite{NIPS2014_RLike_Shao} is employed, and the number of topics is $200$. To understand the temporal dynamic per topic, three topics are selected and the normalized probability that a topic appears at each year are shown in Figure~\ref{fig:performance} (Right). Their associated top 6 words per topic are shown in Table~\ref{table:stu_topics}. The learned trajectory exhibits different temporal patterns
across the topics. Clearly, we can identify jumps associated with some key historical events. For instance, for Topic 29, we observe a positive jump in 1986 related to military and paramilitary activities in and against Nicaragua brought by the U.S. Topic 30 is related with war, where the War of 1812, World War II and Iraq War all spike up in their corresponding years. In Topic 130, we observe consistent positive jumps from 1890 to 1920, when the American revolution was taking place. Three other interesting topics are also shown in Table~\ref{table:stu_topics}. Topic 64 appears to be related to education, Topic 70 is about Iraq, and Topic 74 is Axis and World War II. We note that the words for these topics are explicitly related to these matters. 
%
% After removing stop words and terms that occur fewer than 3 times in one document or less than 10 times overall, there are 7518 unique words. The number of hiddien units $J=200$, the bias terms are initialized as $-1$ to impose sparsity. 
\vskip -0.1in
\begin{table}[h]
\begin{center}
	\caption{\small{Top 6 most probable words associated with the STU topics. }}
	\vskip -0.1in
	\begin{adjustbox}{scale=0.9,tabular=cccccc,center}
		\hline
		Topic $ \# 29$ &  Topic $ \# 30$ & Topic $ \# 130$ & 
		Topic $ \# 64$ & Topic $ \# 70$ & Topic $ \# 74$  \\
		\hline
		family &   officer     & government & generations & Iraqi & Philippines   \\
		budget &  civilized  & country       &  generation & Qaida & islands   \\
		Nicaragua &  warfare  & public      & recognize & Iraq & axis   \\
		free        &  enemy  & law             & brave & Iraqis & Nazis   \\
		future     &  whilst & present         & crime & AI & Japanese   \\
		freedom & gained  & citizens       & race  & Saddam & Germans    \\
		%excellence & lake & united           & balanced & ballistic & mines   \\
		%drugs    &  safety & house          & streets  & terrorists & sailors   \\
		%families & American & foreign     & college & hussein & Nazi   \\
		%God       &  militia & gentlemen   &  school & failures & hats   \\
		\hline
	\end{adjustbox} 	\label{table:stu_topics}
\end{center}
\vspace{-5mm}
\end{table}  


gmodel_supp.tex
%\documentclass{article} % For LaTeX2e
%\usepackage{nips15submit_e,times}
%\include{subtex/mlVecMat}

%\begin{document}

\begin{figure}[h]
\begin{center}
\begin{adjustbox}{minipage=0.7\linewidth,scale=0.85}
\begin{tikzpicture}[ bend angle=45, >=latex, font=\normalsize ]
\tikzstyle{obs} = [ circle, thick, draw = black!100, fill = blue!20, minimum size = 6mm ]
\tikzstyle{lat} = [ circle, thick, draw = black!100, fill = red!0, minimum size = 6mm ]
\tikzstyle{lat2} = [ circle, thick, draw = black!5, fill = black!5, minimum size = 0mm ]
\tikzstyle{par} = [ circle, draw, fill = black!100, minimum width = 1pt, inner sep = 0pt]
\tikzstyle{every label} = [black!100]
\begin{scope}[ node distance = 1.25cm,rounded corners = 4pt ] %
%------------------- the first part ---------------------------
% \zv
\node [lat] (dgz0) [] {};
\node [lat] (dgz1) [ right of = dgz0, node distance = 1.6 cm ] {};
\node [lat] (dgz2) [ right of = dgz1, node distance = 1.6 cm ] {};
% \hv
\node [lat] (dgh0) [ below of = dgz0, label = 180:, node distance = 1.6 cm] {};
\node [lat] (dgh1) [ below of = dgz1, node distance = 1.6 cm] {};
\node [lat] (dgh2) [ below of = dgz2, node distance = 1.6 cm] {};
% \vv
\node [obs] (dgv0) [ below of = dgh0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (dgv1) [ below of = dgh1, node distance = 1.6 cm] {};
\node [obs] (dgv2) [ below of = dgh2, node distance = 1.6 cm] {};
% line from z0 -> z1
\draw [->,line width=1.0pt,] (dgz0) -- (dgz1); 
% line from z1 -> z2
\draw [->,line width=1.0pt,] (dgz1) -- (dgz2);
% line from h0 -> h1
\draw [->,line width=1.0pt,] (dgh0) -- (dgh1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (dgh1) -- (dgh2);
% line from v0 -> v1
\draw [->,line width=1.0pt,] (dgv0) -- (dgv1); 
% line from v1 -> v2
\draw [->,line width=1.0pt,] (dgv1) -- (dgv2); 
% line from z0 -> h0
\draw [->,line width=1.0pt,] (dgz0) -- (dgh0); 
% line from h0 -> v0
\draw [->,line width=1.0pt,] (dgh0) -- (dgv0); 
% line from z1 -> h1
\draw [->,line width=1.0pt,] (dgz1) -- (dgh1); 
% line from h1 -> v1
\draw [->,line width=1.0pt,] (dgh1) -- (dgv1); 
% line from z2 -> h2
\draw [->,line width=1.0pt,] (dgz2) -- (dgh2); 
% line from h2 -> v2
\draw [->,line width=1.0pt,] (dgh2) -- (dgv2);
% line from h0 -> z1
\draw [->,line width=1.0pt,] (dgh0) -- (dgz1); 
% line from h1 -> z2
\draw [->,line width=1.0pt,] (dgh1) -- (dgz2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (dgv0) -- (dgh1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (dgv1) -- (dgh2);
% words
\draw (0.8,0.2) node {{$\Wmat_1$}};
\draw (-0.3,-0.8) node {{$\Wmat_2$}};
\draw (0.5,-0.6) node {{$\Wmat_3$}};
\draw (0.8,-1.4) node {{$\Wmat_4$}};
\draw (-0.3,-2.4) node {{$\Wmat_5$}};
\draw (0.5,-2.2) node {{$\Wmat_6$}};
\draw (0.8,-3.0) node {{$\Wmat_7$}};
\draw (1.6,-4.0) node {(a) Generative model};
%
%
%------------------- the second part ---------------------------
% \zv
\node [lat] (drz0) [right of = dgz2, node distance = 1.6 cm] {};
\node [lat] (drz1) [ right of = drz0, node distance = 1.6 cm ] {};
\node [lat] (drz2) [ right of = drz1, node distance = 1.6 cm ] {};
% \hv
\node [lat] (drh0) [ below of = drz0, label = 180:, node distance = 1.6 cm] {};
\node [lat] (drh1) [ below of = drz1, node distance = 1.6 cm] {};
\node [lat] (drh2) [ below of = drz2, node distance = 1.6 cm] {};
% \vv
\node [obs] (drv0) [ below of = drh0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (drv1) [ below of = drh1, node distance = 1.6 cm] {};
\node [obs] (drv2) [ below of = drh2, node distance = 1.6 cm] {};
% line from z0 -> z1
\draw [->,line width=1.0pt,] (drz0) -- (drz1); 
% line from z1 -> z2
\draw [->,line width=1.0pt,] (drz1) -- (drz2);
% line from h0 -> h1
\draw [->,line width=1.0pt,] (drh0) -- (drh1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (drh1) -- (drh2);
% line from z0 -> h0
\draw [->,line width=1.0pt,] (drh0) -- (drz0); 
% line from h0 -> v0
\draw [->,line width=1.0pt,] (drv0) -- (drh0); 
% line from z1 -> h1
\draw [->,line width=1.0pt,] (drh1) -- (drz1); 
% line from h1 -> v1
\draw [->,line width=1.0pt,] (drv1) -- (drh1); 
% line from z2 -> h2
\draw [->,line width=1.0pt,] (drh2) -- (drz2); 
% line from h2 -> v2
\draw [->,line width=1.0pt,] (drv2) -- (drh2);
% line from h0 -> z1
\draw [->,line width=1.0pt,] (drh0) -- (drz1); 
% line from h1 -> z2
\draw [->,line width=1.0pt,] (drh1) -- (drz2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (drv0) -- (drh1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (drv1) -- (drh2);
% words
\draw (5.6,0.2) node {{$\Umat_1$}};
\draw (4.5,-0.8) node {{$\Umat_2$}};
\draw (5.3,-0.6) node {{$\Umat_3$}};
\draw (5.6,-1.4) node {{$\Umat_4$}};
\draw (4.5,-2.4) node {{$\Umat_5$}};
\draw (5.3,-2.2) node {{$\Umat_6$}};
\draw (6.4,-4.0) node {(b) Recognition model};
\end{scope}
\end{tikzpicture}
\end{adjustbox}
\end{center}
\vspace{-2mm}
\caption{\small{Generative and recognition model of a two-layer Deep TSBN.}} \label{fig:gmodel}
	%\vspace{-5mm}
\end{figure}

%\end{document}

graphicalmodel.tex
%\documentclass{article} % For LaTeX2e
%\usepackage{nips15submit_e,times}
%\include{subtex/mlVecMat}

%\begin{document}

\begin{figure}[t!]
\begin{center}
\begin{adjustbox}{minipage=1.08\linewidth,scale=0.85}
\begin{tikzpicture}[ bend angle=45, >=latex, font=\normalsize ]
\tikzstyle{obs} = [ circle, thick, draw = black!100, fill = blue!20, minimum size = 6mm ]
\tikzstyle{lat} = [ circle, thick, draw = black!100, fill = red!0, minimum size = 6mm ]
\tikzstyle{lat2} = [ circle, thick, draw = black!5, fill = black!5, minimum size = 0mm ]
\tikzstyle{par} = [ circle, draw, fill = black!100, minimum width = 1pt, inner sep = 0pt]
\tikzstyle{every label} = [black!100]
\begin{scope}[ node distance = 1.25cm,rounded corners = 4pt ] %
%------------------- the first part ---------------------------
% \hv
\node [lat] (h0) [] {};
\node [lat2] (hidden) [ below of = h0, node distance = 0.8 cm ] {};
\node [lat] (h1) [ right of = h0, node distance = 2.0 cm ] {};
\node [lat] (h2) [ right of = h1, node distance = 2.0 cm ] {};
% \vv
\node [obs] (v0) [ below of = h0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (v1) [ below of = h1, node distance = 1.6 cm] {};
\node [obs] (v2) [ below of = h2, node distance = 1.6 cm] {};
% line from h0 -> h1
\draw [->,line width=1.0pt,] (h0) -- (h1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (h1) -- (h2);
% line from v0 -> v1
\draw [->,line width=1.0pt,] (v0) -- (v1);
% line from v1 -> v2
\draw [->,line width=1.0pt,] (v1) -- (v2);
% line from h0 -> v0
\draw [->,line width=1.0pt,] (h0) -- (v0); 
% line from h1 -> v1
\draw [->,line width=1.0pt,] (h1) -- (v1); 
% line from h2 -> v2
\draw [->,line width=1.0pt,] (h2) -- (v2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (v0) -- (h1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (v1) -- (h2); 
% word from h0 -> h1
\draw (1,0.2) node {{$\Wmat_1$}};
% word from h0 -> v1
\draw (-0.3,-0.8) node {{$\Wmat_2$}};
% word from v0 -> h1
\draw (0.7,-0.6) node {{$\Wmat_3$}};
% word from v0 -> h1
\draw (1,-1.4) node {{$\Wmat_4$}};
% words
\draw (2,-2.4) node {(a) Generative model};
%
%
%------------------- the second part ---------------------------
% \hv
\node [lat] (rh0) [below of = v0, node distance = 1.6 cm] {};
\node [lat] (rh1) [ right of = rh0, node distance = 2.0 cm ] {};
\node [lat] (rh2) [ right of = rh1, node distance = 2.0 cm ] {};
% \vv
\node [obs] (rv0) [ below of = rh0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (rv1) [ below of = rh1, node distance = 1.6 cm] {};
\node [obs] (rv2) [ below of = rh2, node distance = 1.6 cm] {};
% line from h0 -> h1
\draw [->,line width=1.0pt,] (rh0) -- (rh1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (rh1) -- (rh2);
% line from v0 -> h0
\draw [->,line width=1.0pt,] (rv0) -- (rh0); 
% line from v1 -> h1
\draw [->,line width=1.0pt,] (rv1) -- (rh1); 
% line from v2 -> h2
\draw [->,line width=1.0pt,] (rv2) -- (rh2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (rv0) -- (rh1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (rv1) -- (rh2); 
% word from h0 -> h1
\draw (1,0.2-3.2) node {{$\Umat_1$}};
% word from h0 -> v1
\draw (-0.3,-0.8-3.2) node {{$\Umat_2$}};
% word from v0 -> h1
\draw (0.7,-0.6-3.2) node {{$\Umat_3$}};
% words
\draw (2,-2.4-3.2) node {(b) Recognition model};
%
%
%------------------- the third part ---------------------------
% \zv
\node [lat] (dgz0) [right of = hidden, node distance = 5.6 cm] {};
\node [lat] (dgz1) [ right of = dgz0, node distance = 1.6 cm ] {};
\node [lat] (dgz2) [ right of = dgz1, node distance = 1.6 cm ] {};
% \hv
\node [lat] (dgh0) [ below of = dgz0, label = 180:, node distance = 1.6 cm] {};
\node [lat] (dgh1) [ below of = dgz1, node distance = 1.6 cm] {};
\node [lat] (dgh2) [ below of = dgz2, node distance = 1.6 cm] {};
% \vv
\node [obs] (dgv0) [ below of = dgh0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (dgv1) [ below of = dgh1, node distance = 1.6 cm] {};
\node [obs] (dgv2) [ below of = dgh2, node distance = 1.6 cm] {};
% line from z0 -> z1
\draw [->,line width=1.0pt,] (dgz0) -- (dgz1); 
% line from z1 -> z2
\draw [->,line width=1.0pt,] (dgz1) -- (dgz2);
% line from h0 -> h1
\draw [->,line width=1.0pt,] (dgh0) -- (dgh1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (dgh1) -- (dgh2);
% line from v0 -> v1
\draw [->,line width=1.0pt,] (dgv0) -- (dgv1); 
% line from v1 -> v2
\draw [->,line width=1.0pt,] (dgv1) -- (dgv2); 
% line from z0 -> h0
\draw [->,line width=1.0pt,] (dgz0) -- (dgh0); 
% line from h0 -> v0
\draw [->,line width=1.0pt,] (dgh0) -- (dgv0); 
% line from z1 -> h1
\draw [->,line width=1.0pt,] (dgz1) -- (dgh1); 
% line from h1 -> v1
\draw [->,line width=1.0pt,] (dgh1) -- (dgv1); 
% line from z2 -> h2
\draw [->,line width=1.0pt,] (dgz2) -- (dgh2); 
% line from h2 -> v2
\draw [->,line width=1.0pt,] (dgh2) -- (dgv2);
% line from h0 -> z1
\draw [->,line width=1.0pt,] (dgh0) -- (dgz1); 
% line from h1 -> z2
\draw [->,line width=1.0pt,] (dgh1) -- (dgz2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (dgv0) -- (dgh1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (dgv1) -- (dgh2);
% words
\draw (7.2,-4.8) node {(c) Generative model};
% time
\draw [->,line width=1.0pt,] (5.6,-0.2) -- (8.8,-0.2);
\draw (7.2,0) node {Time};
%
%
%------------------- the fourth part ---------------------------
% \zv
\node [lat] (drz0) [right of = dgz2, node distance = 1.6 cm] {};
\node [lat] (drz1) [ right of = drz0, node distance = 1.6 cm ] {};
\node [lat] (drz2) [ right of = drz1, node distance = 1.6 cm ] {};
% \hv
\node [lat] (drh0) [ below of = drz0, label = 180:, node distance = 1.6 cm] {};
\node [lat] (drh1) [ below of = drz1, node distance = 1.6 cm] {};
\node [lat] (drh2) [ below of = drz2, node distance = 1.6 cm] {};
% \vv
\node [obs] (drv0) [ below of = drh0, label = 180:, node distance = 1.6 cm] {};
\node [obs] (drv1) [ below of = drh1, node distance = 1.6 cm] {};
\node [obs] (drv2) [ below of = drh2, node distance = 1.6 cm] {};
% line from z0 -> z1
\draw [->,line width=1.0pt,] (drz0) -- (drz1); 
% line from z1 -> z2
\draw [->,line width=1.0pt,] (drz1) -- (drz2);
% line from h0 -> h1
\draw [->,line width=1.0pt,] (drh0) -- (drh1); 
% line from h1 -> h2
\draw [->,line width=1.0pt,] (drh1) -- (drh2);
% line from z0 -> h0
\draw [->,line width=1.0pt,] (drh0) -- (drz0); 
% line from h0 -> v0
\draw [->,line width=1.0pt,] (drv0) -- (drh0); 
% line from z1 -> h1
\draw [->,line width=1.0pt,] (drh1) -- (drz1); 
% line from h1 -> v1
\draw [->,line width=1.0pt,] (drv1) -- (drh1); 
% line from z2 -> h2
\draw [->,line width=1.0pt,] (drh2) -- (drz2); 
% line from h2 -> v2
\draw [->,line width=1.0pt,] (drv2) -- (drh2);
% line from h0 -> z1
\draw [->,line width=1.0pt,] (drh0) -- (drz1); 
% line from h1 -> z2
\draw [->,line width=1.0pt,] (drh1) -- (drz2); 
% line from v0 -> h1
\draw [->,line width=1.0pt,] (drv0) -- (drh1); 
% line from v1 -> h2
\draw [->,line width=1.0pt,] (drv1) -- (drh2);
% words
\draw (12.0,-4.8) node {(d) Recognition model};
% time
\draw [->,line width=1.0pt,] (10.4,-0.2) -- (13.6,-0.2);
\draw (12.0,0) node {Time};

%------------ draw the background -----------------------------
\begin{pgfonlayer}{background}
\filldraw [line width = 1pt, draw=black!50, fill=black!5]
(4.5cm,0.5cm)  rectangle (-0.7cm,-6.0cm); % right, top, left, bottom
\end{pgfonlayer}

\begin{pgfonlayer}{background}
\filldraw [line width = 1pt, draw=black!50, fill=black!5]
(14.4cm,0.5cm)  rectangle (5.0cm,-5.4cm); % right, top, left, bottom
\end{pgfonlayer}

%\draw (-3.1,-2.0) node {{$_{k=1,\ldots,K}$}};
%\begin{pgfonlayer}{background}
%\filldraw [line width = 1pt, draw=black!50, fill=black!5]
%(-2.4cm,0.5cm)  rectangle (-3.8cm,-2.2cm); % right, top, left, bottom
%\end{pgfonlayer}

\end{scope}
\end{tikzpicture}
\end{adjustbox}
\end{center}
\vspace{-2mm}
\caption{\small{Graphical model for the Deep Temporal Sigmoid Belief Network. (a,b) Generative and recognition model of the TSBN. (c,d) Generative and recognition model of a two-layer Deep TSBN.}} \label{fig:gmodel}
	\vspace{-5mm}
\end{figure}

%\end{document}

inference.tex
%!TEX root = ../DeepTSBN.tex

\section{Scalable Learning and Inference} \label{sec:inference} 
%
Computation of the exact posterior over the hidden variables in \eqref{eq:model} is intractable. Approximate Bayesian inference, such as Gibbs sampling or mean-field variational Bayes (VB) inference, can be implemented \cite{gan2015learning,neal1992connectionist}. However, Gibbs sampling is very inefficient, due to the fact that the conditional posterior distribution of the hidden variables does not factorize. The mean-field VB indeed provides a fully factored variational posterior, but this technique increases the gap between the bound being optimized and the true log-likelihood, potentially resulting in a poor fit to the data. To allow for tractable and scalable inference and parameter learning, without loss of the flexibility of the variational posterior, we apply the \emph{Neural Variational Inference and Learning} (NVIL) algorithm described in \cite{mnih2014neural}.
%
\subsection{Variational Lower Bound Objective}
%
We are interested in training the TSBN model, $p_{\thetav}(\Vmat,\Hmat)$, described in (\ref{eq:model}), with parameters $\thetav$. Given an observation $\Vmat$, we introduce a fixed-form distribution, $q_{\phiv}(\Hmat|\Vmat)$, with parameters $\phiv$, that approximates the true posterior distribution, $p(\Hmat|\Vmat)$. We then follow the variational principle to derive a lower bound on the marginal log-likelihood, expressed as\footnote{This lower bound is equivalent to the marginal log-likelihood if $q_{\phiv}(\Hmat|\Vmat)$ = $p(\Hmat|\Vmat)$.}
%
\begin{align} \label{lbound}
\Lcal(\Vmat,\thetav,\phiv) = \E_{q_{\phiv}(\Hmat|\Vmat)} [\log p_{\thetav}(\Vmat,\Hmat)-\log q_{\phiv}(\Hmat|\Vmat)] \,.
\end{align}
%
We construct the approximate posterior $q_{\phiv}(\Hmat|\Vmat)$ as a recognition model. By using this, we avoid the need to compute variational parameters per data point; instead we compute a set of parameters $\phiv$ used for all $\Vmat$. In order to achieve fast inference, the recognition model is expressed as
%
\begin{align}
q_{\phiv} (\Hmat|\Vmat) = q(\hv_1|\vv_1)\cdot \prod_{t=2}^T q(\hv_t|\hv_{t-1},\vv_t,\vv_{t-1}) \,,
\end{align}
%
and each conditional distribution is specified as
%
\begin{align} \label{recong}
 q(h_{jt}=1|\hv_{t-1},\vv_t,\vv_{t-1}) = \sigma(\uv_{1j}^\top \hv_{t-1} + \uv_{2j}^\top \vv_t + \uv_{3j}^\top \vv_{t-1}+d_j) \,,
\end{align}
%
where $\hv_0$ and $\vv_0$, for $q(\hv_1|\vv_1)$, are defined as zero vectors. The recognition parameters $\phiv$ are specified as $\Umat_1 \in \R^{J\times J}$, $\Umat_2 \in \R^{J\times M}$, $\Umat_3 \in \R^{J\times M}$. For $i=1,2,3$, $\uv_{ij}$ is the transpose of the $j$th row of $\Umat_i$, and $\dv=[d_1,\ldots,d_J]^\top$ is the bias term. The graphical model is shown in Figure \ref{fig:gmodel}(b).

The recognition model defined in (\ref{recong}) has the same form as in the approximate inference used for the TRBM \cite{sutskever2007learning}. Exact inference for our model consists of a forward and backward pass through the entire sequence, that requires the traversing of each possible hidden state. Our feedforward approximation allows the inference procedure to be fast and implemented in an online fashion.
%
%the computational complexity of which is $O(T\cdot 2^{2J})$
%
\subsection{Parameter Learning}
%
To optimize (\ref{lbound}), we utilize Monte Carlo methods to approximate expectations and stochastic gradient descent (SGD) for parameter optimization. The gradients can be expressed as
%
\begin{align}
\nabla_{\thetav} \Lcal(\Vmat) &= \E_{q_{\phiv}(\Hmat|\Vmat)} [\nabla_{\thetav} \log p_{\thetav}(\Vmat,\Hmat)], \label{model_rule} \\
\nabla_{\phiv} \Lcal(\Vmat) &= \E_{q_{\phiv}(\Hmat|\Vmat)} [(\log p_{\thetav}(\Vmat,\Hmat)-\log q_{\phiv}(\Hmat|\Vmat)) \times \nabla_{\phiv} \log q_{\phiv} (\Hmat|\Vmat) ]. \label{recogn_rule}
\end{align}
%
Specifically, in the TSBN model, if we define $\hat{v}_{mt} = \sigma(\wv_{2m}^\top \hv_t  + \wv_{4m}^\top \vv_{t-1} + c_m)$ and $\hat{h}_{jt} = \sigma(\uv_{1j}^\top \hv_{t-1} + \uv_{2j}^\top \vv_t + \uv_{3j}^\top \vv_{t-1} + d_j)$, the gradients for $\wv_{2m}$ and $\uv_{2j}$ can be calculated as
%
\begin{align} \label{gradient}
\frac{\partial \log p_{\thetav}(\Vmat,\Hmat) }{\partial w_{2mj}} = \sum_{t=1}^T ( v_{mt} - \hat{v}_{mt} ) \cdot h_{jt}, \qquad
\frac{\partial \log q_{\phiv} (\Hmat|\Vmat) }{\partial u_{2jm}} = \sum_{t=1}^T (h_{jt} - \hat{h}_{jt} ) \cdot v_{mt}.
\end{align}
%
Other update equations, along with the learning details for the TSBN variants in Section \ref{sec:tsbn_variants}, are provided in the Supplementary Section B. We observe that the gradients in (\ref{model_rule}) and (\ref{recogn_rule}) share many similarities with the \emph{wake-sleep} algorithm \cite{hinton1995wake}. Wake-sleep alternates between updating $\thetav$ in the wake phase and updating $\phiv$ in the sleep phase. The update of $\thetav$ is based on the samples generated from $q_{\phiv}(\Hmat|\Vmat)$, and is identical to (\ref{model_rule}). However, in contrast to (\ref{recogn_rule}), the recognition parameters $\phiv$ are estimated from samples generated by the model, \emph{i.e.}, $\nabla_{\phiv} \Lcal(\Vmat) = \E_{p_{\thetav}(\Vmat,\Hmat)} [\nabla_{\phiv} \log q_{\phiv} (\Hmat|\Vmat)]$. This update does not optimize the same objective as in (\ref{model_rule}), hence the wake-sleep algorithm is not guaranteed to converge \cite{mnih2014neural}.

Inspecting (\ref{recogn_rule}), we see that we are using $l_{\phiv}(\Vmat,\Hmat)= \log p_{\thetav}(\Vmat,\Hmat)-\log q_{\phiv}(\Hmat|\Vmat)$ as the \emph{learning signal} for the recognition parameters $\phiv$. The expectation of this learning signal is exactly the lower bound (\ref{lbound}), which is easy to evaluate. However, this tractability makes the estimated gradients of the recognition parameters very noisy. In order to make the algorithm practical, we employ the variance reduction techniques proposed in \cite{mnih2014neural}, namely: (\emph{i}) centering the learning signal, by subtracting the data-independent baseline and the data-dependent baseline; (\emph{ii}) variance normalization, by dividing the centered learning signal by a running estimate of its standard deviation. The data-dependent baseline is implemented using a neural network. Additionally, RMSprop \cite{tieleman2012lecture}, a form of SGD where the gradients are adaptively rescaled by a running average of their recent magnitude, were found in practice to be important for fast convergence; thus utilized throughout all the experiments. The outline of the NVIL algorithm is provided in the Supplementary Section A.
%
\subsection{Extension to deep models} \label{sec:deep_inference}
%
The recognition model corresponding to the deep TSBN is shown in Figure \ref{fig:gmodel}(d). Two kinds of deep architectures are discussed in Section \ref{sec:deep_model}. We illustrate the difference of their learning algorithms in two respects: (\emph{i}) the calculation of the lower bound; and (\emph{ii}) the calculation of the gradients. 

The top hidden layer is stochastic. If the middle hidden layers are also stochastic, the calculation of the lower bound is more involved, compared with the shallow model; however, the gradient evaluation remain simple as in (\ref{gradient}). On the other hand, if deterministic middle hidden layers (\emph{i.e.}, recurrent neural networks) are employed, the lower bound objective will stay the same as a shallow model, since the only stochasticity in the generative process lies in the top layer; however, the gradients have to be calculated recursively through the \emph{back-propagation through time} algorithm \cite{werbos1990backpropagation}. All details are provided in the Supplementary Section C.

introduction.tex
%!TEX root = ../DeepTSBN.tex

\section{Introduction} \label{sec:intro}
%
Considerable research has been devoted to developing probabilistic models for high-dimensional time-series data, such as video and music sequences, motion capture data, and text streams. Among them, Hidden Markov Models (HMMs) \cite{rabiner1986introduction} and Linear Dynamical Systems (LDS) \cite{kalman1963mathematical} have been widely studied, but they may be limited in the type of dynamical structures they can model. An HMM is a mixture model, which relies on a single multinomial variable to represent the history of a time-series.  To represent $N$ bits of information about the history, an HMM could require $2^N$ distinct states. On the other hand, real-world sequential data often contain complex non-linear temporal dependencies, while a LDS can only model simple linear dynamics.



Another class of time-series models, which are potentially better suited to model complex probability distributions over high-dimensional sequences, relies on the use of Recurrent Neural Networks (RNNs) \cite{hermans2013training,martens2011learning,pascanu2012difficulty,graves2013generating}, and variants of a well-known \emph{undirected} graphical model called the Restricted Boltzmann Machine (RBM) \cite{taylor2006modeling,sutskever2007learning,sutskever2009recurrent,boulanger2012modeling,mittelman2014structured}. One such variant is the Temporal Restricted Boltzmann Machine (TRBM) \cite{sutskever2007learning}, which consists of a sequence of RBMs, where the state of one or more previous RBMs determine the biases of the RBM in the current time step. Learning and inference in the TRBM is non-trivial. The approximate procedure used in \cite{sutskever2007learning} is heuristic and not derived from a principled statistical formalism.

Recently, deep \emph{directed} generative models \cite{kingma2013auto,mnih2014neural,rezende2014stochastic,gan2015learning} are becoming popular. A \emph{directed} graphical model that is closely related to the RBM is the Sigmoid Belief Network (SBN) \cite{neal1992connectionist}.  In the work presented here, we introduce the Temporal Sigmoid Belief Network (TSBN), which can be viewed as a temporal stack of SBNs, where each SBN has a contextual hidden state that is inherited from the previous SBNs and is used to adjust its hidden-units bias. Based on this, we further develop a deep dynamic generative model by constructing a hierarchy of TSBNs. This can be considered as a deep SBN \cite{gan2015learning} with temporal feedback loops on each layer. Both stochastic and deterministic hidden layers are considered.

Compared with previous work, our model: (\emph{i}) can be viewed as a generalization of an HMM with distributed hidden state representations, and with a deep architecture; (\emph{ii}) can be seen as a generalization of a LDS with complex non-linear dynamics; (\emph{iii}) can be considered as a probabilistic construction of the traditionally deterministic RNN; (\emph{iv}) is closely related to the TRBM, but it has a fully generative process, where data are readily generated from the model using ancestral sampling; (\emph{v}) can be utilized to model different kinds of data, \emph{e.g.}, binary, real-valued and counts.

The ``explaining away'' effect described in \cite{hinton2006fast} makes inference slow, if one uses traditional inference methods. Another important contribution we present here is to develop fast and scalable learning and inference algorithms, by introducing a recognition model \cite{kingma2013auto,mnih2014neural,rezende2014stochastic}, that learns an inverse mapping from observations to hidden variables, based on a loss function derived from a variational principle. By utilizing the recognition model and variance-reduction techniques from \cite{mnih2014neural}, we achieve fast inference both at training and testing time.

mlVecMat.tex
\usepackage{hyperref}
\usepackage{url}
%\documentstyle[nips14submit_09,times,art10]{article} % For LaTeX 2.09

\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{wrapfig}
\usepackage{bm}
\usepackage{subcaption}
\usepackage{psfrag}
\usepackage{tikz}
\usepackage{epstopdf}

\usepackage{color}
\usepackage{tikz}
\usetikzlibrary{arrows,shapes,snakes,automata,backgrounds,fit,petri}
\usepackage{adjustbox}

% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}


\newcommand{\beq}{\vspace{0mm}\begin{equation}}
\newcommand{\eeq}{\vspace{0mm}\end{equation}}
\newcommand{\beqs}{\vspace{0mm}\begin{eqnarray}}
\newcommand{\eeqs}{\vspace{0mm}\end{eqnarray}}
\newcommand{\barr}{\begin{array}}
\newcommand{\earr}{\end{array}}
\newcommand{\Amat}[0]{{{\bf A}}}
\newcommand{\Bmat}{{\bf B}}
\newcommand{\Cmat}{{\bf C}}
\newcommand{\Dmat}{{\bf D}}
\newcommand{\Emat}[0]{{{\bf E}}}
\newcommand{\Fmat}[0]{{{\bf F}}\xspace}
\newcommand{\Gmat}{{\bf G}}
\newcommand{\Hmat}{{\bf H}}
\newcommand{\Imat}{{\bf I}}
\newcommand{\Jmat}[0]{{{\bf J}}\xspace}
\newcommand{\Kmat}[0]{{{\bf K}}\xspace}
\newcommand{\Lmat}[0]{{{\bf L}}}
%\newcommand{\Mmat}[0]{{{\bf M}}\xspace}
\newcommand{\Mmat}{{\bf M}}
\newcommand{\Nmat}[0]{{{\bf N}}\xspace}
\newcommand{\Omat}[0]{{{\bf O}}}
\newcommand{\Pmat}{{\bf P}}
\newcommand{\Qmat}[0]{{{\bf Q}}\xspace}
\newcommand{\Rmat}[0]{{{\bf R}}}
\newcommand{\Smat}[0]{{{\bf S}}}
\newcommand{\Tmat}[0]{{{\bf T}}}
\newcommand{\Umat}[0]{{{\bf U}}}
\newcommand{\Vmat}[0]{{{\bf V}}}
\newcommand{\Wmat}[0]{{{\bf W}}}
\newcommand{\Xmat}[0]{{{\bf X}}}
\newcommand{\Ymat}{{\bf Y}}
%\newcommand{\Ymat}[0]{{{\bf Z}}}
\newcommand{\Zmat}{{\bf Z}}

\newcommand{\av}[0]{{\boldsymbol{a}}}
\newcommand{\bv}[0]{{\boldsymbol{b}}}
\newcommand{\cv}[0]{{\boldsymbol{c}}}
\newcommand{\dv}{\boldsymbol{d}}
\newcommand{\ev}[0]{{\boldsymbol{e}}\xspace}
\newcommand{\fv}[0]{{\boldsymbol{f}}\xspace}
\newcommand{\gv}[0]{{\boldsymbol{g}}\xspace}
\newcommand{\hv}[0]{{\boldsymbol{h}}}
\newcommand{\iv}[0]{{\boldsymbol{i}}\xspace}
\newcommand{\jv}[0]{{\boldsymbol{j}}\xspace}
\newcommand{\kv}[0]{{\boldsymbol{k}}\xspace}
\newcommand{\lv}[0]{{\boldsymbol{l}}}
\newcommand{\mv}[0]{{\boldsymbol{m}}}
\newcommand{\nv}[0]{{\boldsymbol{n}}\xspace}
\newcommand{\ov}[0]{{\boldsymbol{o}}\xspace}
\newcommand{\pv}[0]{{\boldsymbol{p}}}
\newcommand{\qv}[0]{{\boldsymbol{q}}\xspace}
\newcommand{\rv}{\boldsymbol{r}}
\newcommand{\sv}[0]{{\boldsymbol{s}}}
\newcommand{\tv}[0]{{\boldsymbol{t}}\xspace}
\newcommand{\uv}{\boldsymbol{u}}
\newcommand{\vv}{\boldsymbol{v}}
\newcommand{\wv}{\boldsymbol{w}}
\newcommand{\xv}{\boldsymbol{x}}
\newcommand{\yv}{\boldsymbol{y}}
\newcommand{\zv}{\boldsymbol{z}}
\newcommand{\cdotv}{\boldsymbol{\cdot}}

\newcommand{\Gammamat}[0]{{\boldsymbol{\Gamma}}\xspace}
\newcommand{\Deltamat}[0]{{\boldsymbol{\Delta}}\xspace}
\newcommand{\Thetamat}{\boldsymbol{\Theta}}
\newcommand{\Betamat}{\boldsymbol{\Beta}}
\newcommand{\Lambdamat}{\boldsymbol{\Lambda}}
\newcommand{\Ximat}[0]{{\boldsymbol{\Xi}}\xspace}
\newcommand{\Pimat}[0]{{\boldsymbol{\Pi}}\xspace}
\newcommand{\Sigmamat}[0]{{\boldsymbol{\Sigma}}}
\newcommand{\Upsilonmat}[0]{{\boldsymbol{\Upsilon}}\xspace}
\newcommand{\Phimat}{\boldsymbol{\Phi}}
\newcommand{\Psimat}{\boldsymbol{\Psi}}
\newcommand{\Omegamat}[0]{{\boldsymbol{\Omega}}}

\newcommand{\alphav}{\boldsymbol{\alpha}}
\newcommand{\betav}[0]{{\boldsymbol{\beta}}}
\newcommand{\gammav}[0]{{\boldsymbol{\gamma}}\xspace}
\newcommand{\deltav}[0]{{\boldsymbol{\delta}}\xspace}
\newcommand{\epsilonv}{\boldsymbol{\epsilon}}
\newcommand{\zetav}{\boldsymbol{\zeta}}
\newcommand{\etav}[0]{{\boldsymbol{\eta}}\xspace}
\newcommand{\ellv}[0]{{\boldsymbol{\ell}}}
\newcommand{\thetav}{\boldsymbol{\theta}}
\newcommand{\iotav}[0]{{\boldsymbol{\iota}}}
\newcommand{\kappav}[0]{{\boldsymbol{\kappa}}\xspace}
\newcommand{\lambdav}[0]{{\boldsymbol{\lambda}}}
\newcommand{\muv}[0]{{\boldsymbol{\mu}}}
\newcommand{\nuv}[0]{{\boldsymbol{\nu}}}
\newcommand{\xiv}[0]{{\boldsymbol{\xi}}\xspace}
\newcommand{\omicronv}[0]{{\boldsymbol{\omicron}}\xspace}
\newcommand{\piv}{\boldsymbol{\pi}}
\newcommand{\rhov}[0]{{\boldsymbol{\rho}}\xspace}
\newcommand{\sigmav}[0]{{\boldsymbol{\sigma}}}
\newcommand{\tauv}[0]{{\boldsymbol{\tau}}\xspace}
\newcommand{\upsilonv}[0]{{\boldsymbol{\upsilon}}\xspace}
\newcommand{\phiv}{\boldsymbol{\phi}}
\newcommand{\chiv}[0]{{\boldsymbol{\chi}}\xspace}
\newcommand{\psiv}{\boldsymbol{\psi}}
\newcommand{\varthetav}{\boldsymbol{\vartheta}}
\newcommand{\omegav}[0]{{\boldsymbol{\omega}}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\Z}{\mathbb{Z}}

\newcommand{\Xcal}{\mathcal{X}}
\newcommand{\Ycal}{\mathcal{Y}}
\newcommand{\NNcal}{\mathcal{N}}
\newcommand{\Lcal}{\mathcal{L}}
\newcommand{\Ncal}{\mathcal{N}}
\newcommand{\Bcal}{\mathcal{B}}
\newcommand{\Dcal}{\mathcal{D}}
\newcommand{\Tcal}{\mathcal{T}}
\newcommand{\Pcal}{\mathcal{P}}
\newcommand{\Ical}{\mathcal{I}}
\newcommand{\Gcal}{\mathcal{G}}
\newcommand{\Qcal}{\mathcal{Q}}
model.tex
%!TEX root = ../DeepTSBN.tex

\section{Model Formulation} \label{sec:model}
%
\input{subtex/graphicalmodel}
%
\subsection{Sigmoid Belief Networks}
%
Deep dynamic generative models are considered, based on the Sigmoid Belief Network (SBN) \cite{neal1992connectionist}. An SBN is a Bayesian network that models a binary visible vector $\vv \in \{0,1\}^M$, in terms of binary hidden variables $\hv \in \{0,1\}^J$ and weights $\Wmat \in \R^{M \times J}$ with
%
\begin{align}
p(v_m=1|\hv) = \sigma(\wv_m^\top \hv + c_m), \qquad p(h_j=1) = \sigma(b_j),
\end{align}
%
where $\vv = [v_1,\ldots,v_M]^\top$, $\hv = [h_1,\ldots,h_J]^\top$, $\Wmat = [\wv_1,\ldots,\wv_M]^\top$, $\cv = [c_1,\ldots,c_M]^\top$, $\bv = [b_1,\ldots,b_J]^\top$, and the logistic function, $\sigma(x) \triangleq 1 / (1 + e^{-x})$. The parameters $\Wmat$, $\bv$ and $\cv$ characterize all data, and the hidden variables, $\hv$, are specific to particular visible data, $\vv$.

The SBN is closely related to the RBM \cite{hinton2002training}, which is a Markov random field with the same bipartite structure as the SBN. The RBM defines a distribution over a binary vector that is proportional to the exponential of its \emph{energy}, defined as $-E(\vv,\hv) = \vv^\top \cv + \vv^\top \Wmat \hv + \hv^\top \bv$. The conditional distributions, $p(\vv|\hv)$ and $p(\hv|\vv)$, in the RBM are factorial, which makes inference fast, while parameter estimation usually relies on an approximation technique known as Contrastive Divergence (CD) \cite{hinton2002training}.

The energy function of an SBN may be written as $-E(\vv,\hv) = \vv^\top \cv + \vv^\top \Wmat \hv + \hv^\top \bv- \sum_m \log (1+\exp(\wv_m^\top\hv+c_m))$. SBNs explicitly manifest the generative process to obtain data, in which the hidden layer provides a directed ``explanation'' for patterns generated in the visible layer. However, the ``explaining away'' effect described in \cite{hinton2006fast} makes inference inefficient, the latter can be alleviated by exploiting recent advances in variational inference methods \cite{mnih2014neural}.

\subsection{Temporal Sigmoid Belief Networks}
%
The proposed Temporal Sigmoid Belief Network (TSBN) model is a sequence of SBNs arranged in such way that at any given time step, the SBN's biases depend on the state of the SBNs in the previous time steps. Specifically, assume we have a length-$T$ binary visible sequence, the $t$th time step of which is denoted $\vv_t \in \{0,1\}^M$. The TSBN describes the joint probability as
%
\begin{align}\label{eq:model}
p_{\thetav}(\Vmat,\Hmat) = p(\hv_1)p(\vv_1|\hv_1) \cdot \prod_{t=2}^T p(\hv_t|\hv_{t-1},\vv_{t-1}) \cdot p(\vv_t|\hv_t,\vv_{t-1}) ,
\end{align}
%
where $\Vmat = [\vv_1,\ldots,\vv_T]$, $\Hmat = [\hv_1,\ldots,\hv_T]$, and each $\hv_t \in \{0,1\}^J$ represents the hidden state corresponding to time step $t$. For $t=1,\ldots,T$, each conditional distribution in (\ref{eq:model}) is expressed as
%
\begin{align}
p(h_{jt}=1|\hv_{t-1},\vv_{t-1}) &= \sigma (\wv_{1j}^\top \hv_{t-1} + \wv_{3j}^\top \vv_{t-1}+b_j), \label{eq:ph} \\
p(v_{mt}=1|\hv_t,\vv_{t-1}) &= \sigma (\wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1}+c_m), \label{eqn:likelihood}
\end{align}
%
where $\hv_0$ and $\vv_0$, needed for the prior model $p(\hv_1)$ and $p(\vv_1|\hv_1)$, are defined as zero vectors, respectively, for conciseness. The model parameters, $\thetav$, are specified as $\Wmat_1 \in \R^{J\times J}$, $\Wmat_2 \in \R^{M\times J}$, $\Wmat_3 \in \R^{J\times M}$, $\Wmat_4 \in \R^{M\times M}$. For $i=1,2,3,4$, $\wv_{ij}$ is the transpose of the $j$th row of $\Wmat_i$, and $\cv=[c_1,\ldots,c_M]^\top$ and $\bv=[b_1,\ldots,b_J]^\top$ are bias terms. The graphical model for the TSBN is shown in Figure \ref{fig:gmodel}(a).

By setting $\Wmat_3$ and $\Wmat_4$ to be zero matrices, the TSBN can be viewed as a Hidden Markov Model \cite{rabiner1986introduction} with an exponentially large state space, that has a compact parameterization of the transition and the emission probabilities. Specifically, each hidden state in the HMM is represented as a \emph{one-hot} length-$J$ vector, while in the TSBN, the hidden states can be any length-$J$ binary vector. We note that the transition matrix is highly structured, since the number of parameters is only quadratic \emph{w.r.t.} $J$. Compared with the TRBM \cite{sutskever2007learning}, our TSBN is fully directed, which allows for \emph{fast sampling} of ``fantasy'' data from the inferred model.

\vspace{-0.05in}
\subsection{TSBN Variants} \label{sec:tsbn_variants} 
\vspace{-0.05in}
%
\paragraph{Modeling real-valued data} 
%
The model above can be readily extended to model real-valued sequence data, by substituting (\ref{eqn:likelihood}) with $p(\vv_t|\hv_t,\vv_{t-1}) = \Ncal (\muv_t,\mbox{diag}(\sigmav_t^2))$, where
%
\begin{align}
\mu_{mt} = \wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1}+c_m, \,\,\,\
\log \sigma_{mt}^2 = (\wv_{2m}^{\prime})^\top \hv_t + (\wv_{4m}^{\prime})^\top \vv_{t-1}+c_m^{\prime},
\end{align}
%
and $\mu_{mt}$ and $\sigma_{mt}^2$ are elements of $\muv_t$ and $\sigmav_t^2$, respectively. $\Wmat_{2}^{\prime}$ and $\Wmat_{4}^{\prime}$ are of the same size of $\Wmat_{2}$ and $\Wmat_{4}$, respectively. Compared with the Gaussian TRBM \cite{sutskever2009recurrent}, in which $\sigma_{mt}$ is fixed to 1, our formalism uses a diagonal matrix to parameterize the variance structure of $\vv_t$.
%
\paragraph{Modeling count data} 
%
We also introduce an approach for modeling time-series data with count observations, by replacing (\ref{eqn:likelihood}) with $p(\vv_t|\hv_t,\vv_{t-1}) = \prod_{m=1}^M y_{mt}^{v_{mt}}$, where
%
\begin{align} \label{softmax}
y_{mt} = \frac{\exp(\wv_{2m}^\top \hv_t + \wv_{4m}^\top \vv_{t-1}+c_m)}{\sum_{m^{\prime}=1}^M \exp(\wv_{2m^{\prime}}^\top \hv_t + \wv_{4m^{\prime}}^\top \vv_{t-1}+c_{m^{\prime}})} \,.
\end{align}
%
This formulation is related to the Replicated Softmax Model (RSM) described in \cite{hinton2009replicated}, however, our approach uses a \emph{directed} connection from the binary hidden variables to the visible counts, while also learning the dynamics in the count sequences.

Furthermore, rather than assuming that $\hv_t$ and $\vv_t$ only depend on $\hv_{t-1}$ and $\vv_{t-1}$, in the experiments, we also allow for connections from the past $n$ time steps of the hidden and visible states, to the current states, $\hv_t$ and $\vv_t$. A sliding window is then used to go through the sequence to obtain $n$ frames at each time. We refer to $n$ as the \emph{order} of the model.

\subsection{Deep Architecture for Sequence Modeling with TSBNs} \label{sec:deep_model}
%
Learning the sequential dependencies with the shallow model in (\ref{eq:model})-(\ref{eqn:likelihood}) may be restrictive. Therefore, we propose two deep architectures to improve its representational power: (\emph{i}) adding stochastic hidden layers; (\emph{ii}) adding deterministic hidden layers. The graphical model for the deep TSBN is shown in Figure \ref{fig:gmodel}(c). Specifically, we consider a deep TSBN with hidden layers $\hv_{t}^{(\ell)}$ for $t=1,\ldots,T$ and $\ell= 1,\ldots,L$. Assume layer $\ell$ contains $J^{(\ell)}$ hidden units, and denote the visible layer $\vv_t=\hv_{t}^{(0)}$ and let $\hv_t^{(L+1)}= \boldsymbol{0}$, for convenience. In order to obtain a proper generative model, the top hidden layer $\hv^{(L)}$ contains stochastic binary hidden variables. 

For the middle layers, $\ell = 1,\ldots,L-1$, if stochastic hidden layers are utilized, the generative process is expressed as
$ p(\hv_t^{(\ell)}) = \prod_{j=1}^{J^{(\ell)}} p(h_{jt}^{(\ell)}|\hv_t^{(\ell+1)},\hv_{t-1}^{(\ell)},\hv_{t-1}^{(\ell-1)})$,
where each conditional distribution is parameterized via a logistic function, as in (\ref{eqn:likelihood}). If deterministic hidden layers are employed, we obtain $\hv_{t}^{(\ell)} = f(\hv_t^{(\ell+1)},\hv_{t-1}^{(\ell)},\hv_{t-1}^{(\ell-1)})$, where $f(\cdot)$ is chosen to be a rectified linear function. Although the differences between these two approaches are minor, learning and inference algorithms can be quite different, as shown in Section \ref{sec:deep_inference}. 

% ; this configuration has been shown to
%work better than sigmoid non-linearities \cite{nair2010rectified}.

relatedwork.tex
%!TEX root = ../DeepTSBN.tex
\section{Related Work} \label{sec:relatedwork}
%
The RBM has been widely used as building block to learn the sequential dependencies in time-series data, \emph{e.g.}, the conditional-RBM-related models \cite{taylor2006modeling,taylor2009factored}, and the temporal RBM \cite{sutskever2007learning}. To make \emph{exact} inference possible, the recurrent temporal RBM was also proposed \cite{sutskever2009recurrent}, and further extended to learn the dependency structure within observations \cite{mittelman2014structured}.

In the work reported here, we focus on modeling sequences based on the SBN \cite{neal1992connectionist}, which recently has been shown to have the potential to build deep generative models \cite{mnih2014neural,gan2015learning,gan2015scalable}. Our work serves as another extension of the SBN that can be utilized to model time-series data. Similar ideas have also been considered in \cite{henderson2010incremental} and \cite{hinton1995helmholtz}. However, in \cite{henderson2010incremental}, the authors focus on grammar learning, and use a feed-forward approximation of the mean-field VB to carry out the inference; while in \cite{hinton1995helmholtz}, the wake-sleep algorithm was developed. We apply the model in a different scenario, and develop a fast and scalable inference algorithm, based on the idea of training a recognition model by leveraging the stochastic gradient of the variational bound.

There exist two main methods for the training of recognition models. The first one, termed Stochastic Gradient Variational Bayes (SGVB), is based on a reparameterization trick \cite{kingma2013auto,rezende2014stochastic}, which can be only employed in models with continuous latent variables, {\emph{e.g.}}, the variational auto-encoder \cite{kingma2013auto} and all the recent recurrent extensions of it \cite{bayer2014learning,fabius2014variational,chung2015recurrent}. The second one, called Neural Variational Inference and Learning (NVIL), is based on the log-derivative trick \cite{mnih2014neural}, which is more general and can also be applicable to models with discrete random variables. The NVIL algorithm has been previously applied to the training of SBN in \cite{mnih2014neural}. Our approach serves as a new application of this algorithm for a SBN-based time-series model.
%
%\textcolor{red}{We also make a connection with the variational recurrent auto-encoder} \cite{fabius2014variational}, where Gaussian-distributed latent variables are applied based on the SGVB algorithm, but their experimental results are more limited. In the work presented here, we apply our framework to the modeling of binary, real-valued and count data.

DeepTSBN.bbl
\begin{thebibliography}{10}

\bibitem{rabiner1986introduction}
L.~Rabiner and B.~Juang.
\newblock An introduction to hidden markov models.
\newblock In {\em ASSP Magazine, IEEE}, 1986.

\bibitem{kalman1963mathematical}
R.~Kalman.
\newblock Mathematical description of linear dynamical systems.
\newblock In {\em J. the Society for Industrial \& Applied Mathematics, Series
  A: Control}, 1963.

\bibitem{hermans2013training}
M.~Hermans and B.~Schrauwen.
\newblock Training and analysing deep recurrent neural networks.
\newblock In {\em NIPS}, 2013.

\bibitem{martens2011learning}
J.~Martens and I.~Sutskever.
\newblock Learning recurrent neural networks with hessian-free optimization.
\newblock In {\em ICML}, 2011.

\bibitem{pascanu2012difficulty}
R.~Pascanu, T.~Mikolov, and Y.~Bengio.
\newblock On the difficulty of training recurrent neural networks.
\newblock In {\em ICML}, 2013.

\bibitem{graves2013generating}
A.~Graves.
\newblock Generating sequences with recurrent neural networks.
\newblock In {\em arXiv:1308.0850}, 2013.

\bibitem{taylor2006modeling}
G.~Taylor, G.~Hinton, and S.~Roweis.
\newblock Modeling human motion using binary latent variables.
\newblock In {\em NIPS}, 2006.

\bibitem{sutskever2007learning}
I.~Sutskever and G.~Hinton.
\newblock Learning multilevel distributed representations for high-dimensional
  sequences.
\newblock In {\em AISTATS}, 2007.

\bibitem{sutskever2009recurrent}
I.~Sutskever, G.~Hinton, and G.~Taylor.
\newblock The recurrent temporal restricted boltzmann machine.
\newblock In {\em NIPS}, 2009.

\bibitem{boulanger2012modeling}
N.~Boulanger-Lewandowski, Y.~Bengio, and P.~Vincent.
\newblock Modeling temporal dependencies in high-dimensional sequences:
  Application to polyphonic music generation and transcription.
\newblock In {\em ICML}, 2012.

\bibitem{mittelman2014structured}
R.~Mittelman, B.~Kuipers, S.~Savarese, and H.~Lee.
\newblock Structured recurrent temporal restricted boltzmann machines.
\newblock In {\em ICML}, 2014.

\bibitem{kingma2013auto}
D.~P. Kingma and M.~Welling.
\newblock Auto-encoding variational {B}ayes.
\newblock In {\em ICLR}, 2014.

\bibitem{mnih2014neural}
A.~Mnih and K.~Gregor.
\newblock Neural variational inference and learning in belief networks.
\newblock In {\em ICML}, 2014.

\bibitem{rezende2014stochastic}
D.~Rezende, S.~Mohamed, and D.~Wierstra.
\newblock Stochastic backpropagation and approximate inference in deep
  generative models.
\newblock In {\em ICML}, 2014.

\bibitem{gan2015learning}
Z.~Gan, R.~Henao, D.~Carlson, and L.~Carin.
\newblock Learning deep sigmoid belief networks with data augmentation.
\newblock In {\em AISTATS}, 2015.

\bibitem{neal1992connectionist}
R.~Neal.
\newblock Connectionist learning of belief networks.
\newblock In {\em Artificial intelligence}, 1992.

\bibitem{hinton2006fast}
G.~Hinton, S.~Osindero, and Y.~Teh.
\newblock A fast learning algorithm for deep belief nets.
\newblock In {\em Neural computation}, 2006.

\bibitem{hinton2002training}
G.~Hinton.
\newblock Training products of experts by minimizing contrastive divergence.
\newblock In {\em Neural computation}, 2002.

\bibitem{hinton2009replicated}
G.~Hinton and R.~Salakhutdinov.
\newblock Replicated softmax: an undirected topic model.
\newblock In {\em NIPS}, 2009.

\bibitem{hinton1995wake}
G.~Hinton, P.~Dayan, B.~Frey, and R.~Neal.
\newblock The ``wake-sleep'' algorithm for unsupervised neural networks.
\newblock In {\em Science}, 1995.

\bibitem{tieleman2012lecture}
T.~Tieleman and G.~Hinton.
\newblock Lecture 6.5-rmsprop: Divide the gradient by a running average of its
  recent magnitude.
\newblock In {\em COURSERA: Neural Networks for Machine Learning}, 2012.

\bibitem{werbos1990backpropagation}
P.~Werbos.
\newblock Backpropagation through time: what it does and how to do it.
\newblock In {\em Proc. of the IEEE}, 1990.

\bibitem{taylor2009factored}
G.~Taylor and G.~Hinton.
\newblock Factored conditional restricted boltzmann machines for modeling
  motion style.
\newblock In {\em ICML}, 2009.

\bibitem{gan2015scalable}
Z.~Gan, C.~Chen, R.~Henao, D.~Carlson, and L.~Carin.
\newblock Scalable deep poisson factor analysis for topic modeling.
\newblock In {\em ICML}, 2015.

\bibitem{henderson2010incremental}
J.~Henderson and I.~Titov.
\newblock Incremental sigmoid belief networks for grammar learning.
\newblock In {\em JMLR}, 2010.

\bibitem{hinton1995helmholtz}
G.~Hinton, P.~Dayan, A.~To, and R.~Neal.
\newblock The helmholtz machine through time.
\newblock In {\em Proc. of the ICANN}, 1995.

\bibitem{bayer2014learning}
J.~Bayer and C.~Osendorfer.
\newblock Learning stochastic recurrent networks.
\newblock In {\em arXiv:1411.7610}, 2014.

\bibitem{fabius2014variational}
O.~Fabius, J.~R. van Amersfoort, and D.~P. Kingma.
\newblock Variational recurrent auto-encoders.
\newblock In {\em arXiv:1412.6581}, 2014.

\bibitem{chung2015recurrent}
J.~Chung, K.~Kastner, L.~Dinh, K.~Goel, A.~Courville, and Y.~Bengio.
\newblock A recurrent latent variable model for sequential data.
\newblock In {\em NIPS}, 2015.

\bibitem{NIPS2014_RLike_Shao}
S.~Han, L.~Du, E.~Salazar, and L.~Carin.
\newblock Dynamic rank factor model for text streams.
\newblock In {\em NIPS}, 2014.

\bibitem{GP_DFA_AISTATS2015}
A.~Acharya, J.~Ghosh, and M.~Zhou.
\newblock Nonparametric {B}ayesian factor analysis for dynamic count matrices.
\newblock In {\em AISTATS}, 2015.

\bibitem{fan2015fast}
K.~Fan, Z.~Wang, J.~Beck, J.~Kwok, and K.~Heller.
\newblock Fast second-order stochastic backpropagation for variational
  inference.
\newblock In {\em NIPS}, 2015.

\end{thebibliography}
