summaries: This paper tackles the important challenge of making sense of temporal measurements made in biological systems. Among other, those have the peculiarity that they are not independent but with a dependency structure, which can be encoded as a graph or a network. The authors claim that their approach, DyMoN is adapted to the many challenges of biological high-throughput data sets: noise, sparsity and lack of temporal resolution. The paper presents three very different use of the method in complex biological systems in Section 3: (i) Calcium imaging of visual cortex neurons, (ii) T--cell development in the thymus, and (iii) Human embryonic stem cell differentiation. Section 4 assesses the performance of the method on simulated data sets as well as on a face-recognition data set. Moreover, the authors demonstrate how the features of the NN can be interrogated to shed new insight about the process under scrutiny. They also show the gain in running time a comapred to other approaches.

summaries: > Even though the paper details the underlying Markovian setup in Section 2, it is unclear to the reader how this knits with the FFNN architecture, for example what are the Markovian functions at hidden layer and output layer. Are they all conditional probabilities? How do you prove that this is what occurs within each node?

abstractText: Biological data often contains measurements of dynamic entities such as cells or organisms in various states of progression. However, biological systems are notoriously difficult to describe analytically due to their many interacting components, and in many cases, the technical challenge of taking longitudinal measurements. This leads to difficulties in studying the features of the dynamics, for examples the drivers of the transition. To address this problem, we present a deep neural network framework we call Dynamics Modeling Network or DyMoN. DyMoN is a neural network framework trained as a deep generative Markov model whose next state is a probability distribution based on the current state. DyMoN is well-suited to the idiosyncrasies of biological data, including noise, sparsity, and the lack of longitudinal measurements in many types of systems. Thus, DyMoN can be trained using probability distributions derived from the data in any way, such as trajectories derived via dimensionality reduction methods, and does not require longitudinal measurements. We show the advantage of learning deep models over shallow models such as Kalman filters and hidden Markov models that do not learn representations of the data, both in terms of learning embeddings of the data and also in terms training efficiency, accuracy and ability to multitask. We perform three case studies of applying DyMoN to different types of biological systems and extracting features of the dynamics in each case by examining the learned model.

authors:
Dario Amodei
Sundaram Ananthanarayanan
Rishita Anubhai
Jingliang Bai
Eric Battenberg
Carl Case
Jared Casper
Bryan Catanzaro
Qiang Cheng
Guoliang Chen

title: Deep speech 2: Endto-end speech recognition in english and mandarin

authors:
Mikhail Belkin
Partha Niyogi

title: Laplacian eigenmaps and spectral techniques for embedding and clustering

venue: In Advances in neural information processing systems,

authors:
Sean C Bendall
Kara L Davis
El-ad David Amir
Michelle D Tadmor
Erin F Simonds
Tiffany J Chen
Daniel K Shenfeld
Garry P Nolan
Dana Peer

title: Single-cell trajectory detection uncovers progression and regulatory coordination in human b cell development

authors:
Miriam Bibel
Jens Richter
Emmanuel Lacroix
Yves-Alain Barde

title: Generation of a defined and uniform population of cns progenitors and neurons from mouse embryonic stem cells

venue: Nature protocols,

authors:
Paolo Bonifazi
Miri Goldin
Michel A Picardo
Isabel Jorquera
A Cattani
Gregory Bianconi
Alfonso Represa
Yehezkel Ben-Ari
Rosa Cossart

title: Gabaergic hub neurons orchestrate synchrony in developing hippocampal

authors:
Konstantinos Bousmalis
George Trigeorgis
Nathan Silberman
Dilip Krishnan
Dumitru Erhan

title: Domain separation networks

authors:
Samuel R Bowman
Luke Vilnis
Oriol Vinyals
Andrew Dai
Rafal Jozefowicz
Samy Bengio

title: Generating sentences from a continuous space

venue: In Proceedings of The 20th SIGNLL Conference on Computational Natural Language Learning,

authors:
Tsai-Wen Chen
Trevor J Wardill
Yi Sun
Stefan R Pulver
Sabine L Renninger
Amy Baohan
Eric R Schreiter
Rex A Kerr
Michael B Orger
Vivek Jayaraman

title: Ultrasensitive fluorescent proteins for imaging neuronal activity

authors:
Junyoung Chung
Kyle Kastner
Laurent Dinh
Kratarth Goel
Aaron C Courville
Yoshua Bengio

title: A recurrent latent variable model for sequential data

authors:
Ronald R. Coifman
Matthew J. Hirn

title: Diffusion maps for changing data

venue: Applied and Computational Harmonic Analysis,

authors:
Alexander Dubbs
James Guevara
Rafael Yuste

title: Moco: fast motion correction for calcium imaging

venue: Frontiers in neuroinformatics,

authors:
G.K. Dziugaite
D.M. Roy
Z. Ghahramani

title: Training generative neural networks via Maximum Mean Discrepancy optimization

authors:
Santiago Fernández
Alex Graves
Jürgen Schmidhuber

title: An application of recurrent neural networks to discriminative keyword spotting

venue: In Proceedings of the 17th International Conference on Artificial Neural Networks,

authors:
Niels Geijsen
Melissa Horoschak
Kitai Kim
Joost Gribnau
Kevin Eggan
George Q Daley

title: Derivation of embryonic germ cells and male gametes from embryonic stem cells

authors:
Anirudh Goyal
Nan Ke
Surya Ganguli
Yoshua Bengio

title: Variational walkback: Learning a transition operator as a stochastic recurrent net

authors:
Arthur Gretton
Karsten M Borgwardt
Malte J Rasch
Bernhard Schölkopf
Alexander Smola

title: A kernel two-sample test

venue: Journal of Machine Learning Research,

authors:
Laleh Haghverdi
Maren Buettner
F Alexander Wolf
Florian Buettner
Fabian J Theis

title: Diffusion pseudotime robustly reconstructs lineage branching

venue: Nature Methods,

authors:
Rafal Józefowicz
Oriol Vinyals
Mike Schuster
Noam Shazeer
Yonghui Wu

title: Exploring the limits of language modeling

authors:
Gabriela Kania
Przemyslaw Blyszczuk
Andrea Jochheim
Michael Ott
Anna M Wobus

title: Generation of glycogen-and albumin-producing hepatocyte-like cells from embryonic stem cells

venue: Biological chemistry,

authors:
Diederik P. Kingma
Jimmy Ba

title: Adam: A method for stochastic optimization

venue: In International Conference for Learning Representations,

authors:
Diederik P Kingma
Max Welling

title: Auto-encoding variational bayes

venue: arXiv preprint arXiv:1312.6114,

authors:
Kevin R Moon
David van Dijk
Zheng Wang
Daniel Burkhardt
William Chen
Antonia van den Elzen
Matthew J Hirn
Ronald R Coifman
Natalia B Ivanova
Guy Wolf
Smita Krishnaswamy

title: Visualizing transitions and structure for high dimensional data

venue: exploration. bioRxiv,

authors:
Toru Nakano
Hiroaki Kodama
Tasuku Honjo

title: In vitro development of primitive and definitive erythrocytes from different precursors

authors:
Sam T Roweis
Lawrence K Saul

title: Nonlinear dimensionality reduction by locally linear embedding

authors:
Manu Setty
Michelle D Tadmor
Shlomit Reich-Zeliger
Omer Angel
Tomer Meir Salame
Pooja Kathail
Kristy Choi
Sean Bendall
Nir Friedman
Dana Pe’er

title: Wishbone identifies bifurcating developmental trajectories from single-cell data

venue: Nature Biotechnology,

authors:
Uri Shaham
Kelly P Stanton
Jun Zhao
Huamin Li
Khadir Raddassi
Ruth Montgomery
Yuval Kluger

title: Removal of batch effects using distribution-matching residual

venue: networks. Bioinformatics,

authors:
Ronen Talmon
Ronald R. Coifman

title: Intrinsic modeling of stochastic dynamical systems using empirical geometry

authors:
Joanna Urban-Ciecko
Alison L Barth

title: Somatostatin-expressing neurons in cortical networks

venue: Nature Reviews Neuroscience,

authors:
D. van Dijk
R. Sharma
J. Nainys
K. Yim
P. Kathail
A. Carr
C. Burdsiak
K.R. Moon
C. Chaffer
D. Pattabiraman
B. Bierie
L. Mazutis
G. Wolf
S. Krishnaswamy
D. Pe’er

title: Recovering gene interactions from single-cell data using data diffusion

authors:
A. Waibel
T. Hanazawa
G. Hinton
K. Shikano
K. Lang

title: Phoneme recognition using time-delay neural networks

venue: IEEE Transactions on Acoustics, Speech and Signal Processing,

authors:
Kilian Q Weinberger
Fei Sha
Lawrence K Saul

title: Learning a kernel matrix for nonlinear dimensionality reduction

venue: In Proceedings of the twenty-first international conference on Machine learning,

authors:
Guy Wolf
Amir Averbuch

title: Linear-projection diffusion on smooth euclidean submanifolds

authors:
Guy Wolf
Aviv Rotbart
Gil David
Amir Averbuch

title: Coarse-grained localized diffusion

title: Somatostatin-expressing inhibitory interneurons in cortical circuits

venue: In Shaham et al

heading: 1 INTRODUCTION

text:
Many biological systems can be considered conceptually as stochastic dynamic processes. Examples include differentiation from stem cells to mature cellular lineages, responses to stimulation or drug perturbations, and responses to external stimulus or signaling. Such processes involve the slow and stochastic change in gene expression programs that alter over time, brain circuitry and neuronal communication where the firing of one set of neurons produces the firings of another set, or cancer progressions in which one subclone of cells out competes other clones. However, such biological systems are complex and in many cases, it is impossible to derive differential equations or analytical models for such processes in order to study them. In this paper, we propose to learn a generative neural network model of biological dynamic systems that we call a Dynamics Modeling Network, or DyMoN.
First, DyMoN provides a representational advantage by serving as an embodiment of the dynamics in lieu of a predetermined model, such as stochastic differential equations. This representation is deep and factored, in the sense that the “logic” of the dynamic transition is broken down into many increasingly abstract steps, each of which can be visualized or examined. In the cell trajectory context, a neural network can learn complex features that correspond to modules of genes that operate together to create a change in the expression program for the next cellular state. Further, the saliency of these features can be inferred from the Jacobian matrices of the neural network. Secondly, DyMoN provides a generative advantage as it generates new trajectories that have not been seen previously in the system. These can augment data and potentially reveal new paths (for instance novel differentiation paths in cellular development) that lead to desired outcomes. Third, it provides a denoising advantage, thanks to the the dimensionality reduction offered in neural networks, which can naturally denoise the high degree of noise in biological data, and in particular in biological trajectories. Fourth, the utilization of a deep network model offers a multitasking advantage: while previous models, such as HMMs, can simulate trajectories, and one may use PCA (or similar methods) to visualize trajectory
information, most previous methods are not designed or equipped to simultaneously learn multi-scale features (i.e., in many levels of abstraction), visualize the intrinsic underlying dynamics, and utilize them to generate new trajectories. Finally, the natural parallelizability of neural networks (e.g., with GPU-based implementations) offers computational advantages to deep learning approaches, as they can be used to process large volumes of noisy data as they become available for training. Furthermore, once trained, DyMoN can quickly generate new trajectories faster than most existing methods.
We have specifically designed DyMoN to learn a model of dynamics from the types of data available in biological systems: snapshot high dimensional data. By snapshot, we mean that biological data (for instance single cell RNA-sequencing data) is often collected only at one or a handful of time-points, thus the dynamics are either inferred from pseudotime (an axis of time progression based on different states of cell differentiation or other cellular programs) or from interpolation between discrete time points. Second, these types of trajectories are not a matter of sequence completion as is often done in a recurrent neural network, but rather the desired goal is to predict the near-future state of a cell given its current state. Thus, we propose a neural network framework which learns fixed-memory stochastic dynamics (e.g., memoryless or n-th order Markov process) in an observed system, and train this network by providing a current state (or n states) along with a probability distribution of next states. We penalize the network using maximal mean discrepancy (MMD) as a probabilistic distance between the desired and generated outputs. Within this framework, we show three variations of architectures that can realize this framework within three vastly different biological systems.
We train the DyMoN on three biological systems: dynamic calcium imaging data of neurons from a mouse visual cortex, mass cytometry of developing T cells in the mouse thymus, and single-cell RNA sequencing of human embryonic stem cells developing in an embryoid body. We interrogate the weights of the model as a ”transparent box” to understand the dynamics of the system learned by the network in order to both verify known gene-gene interactions in T cells and to uncover novel network interactions in neurons. We further show that the DyMoN is able to generate trajectories in biological data originally measured in discrete snapshots, and use this generative capacity to develop novel hypotheses for cellular programming in embryonic stem cell differentiation.

heading: 2 DYNAMICS MODELING NETWORK (DYMON)

text:
Let X ⊆ Rd be a finite dataset of d-dimensional states, and let T ⊆ X ×X be a (finite) collection of transitions between these states. Such transitions can either be observed by sampling a dynamical system, or constructed by geometry-revealing diffusion methods such as diffusion maps (Coifman & Lafon, 2006), Laplacian eigenmaps (Belkin & Niyogi, 2002), diffusion geometry methods (e.g., Wolf et al., 2012; Wolf & Averbuch, 2013; Coifman & Hirn, 2014; Talmon & Coifman, 2015), and PHATE (Moon et al., 2017). We propose a deep learning approach for learning the dynamics in T and the geometry represented by them as a stochastic velocity vector field of a Markov process. To do this, we use a feed-forward neural network we call Dynamics Modeling Network (DyMoN).
DyMoN is formed by a cascade of linear operations and nonlinear activations. These are controlled by optimized network weights, which we collectively denote by θ ∈ Θ, where Θ represents the space of possible weight values determined by the fixed network architecture. To capture the dynamics in T , DyMoN learns a velocity vector ∆θ(x) ∈ Rd and uses it to define a transition function T : Rd ×Θ→ Rd that generates a Markov process
xt = T (xt−1, θ) = xt−1 + ∆θ(xt−1) t = 1, 2, 3, . . . ,
given an initial state x0 ∈ X and optimized weights θ ∈ Θ. Further, to introduce stochasticity, we treat the output ∆θ(xt−1) as a random vector whose probability distribution determines the conditional probabilities P (xt|xt−1) of the Markov process. We provide three possible implementations of this framework in Figure 10.
Given a fixed architecture, we now describe how DyMoN is trained using the data points in X . Consider a source state x ∈ X with multiple transitions from it leading to target states Yx = {y : (x, y) ∈ T }. Let Px(y) = P (y | x), whose support is Yx, be the conditional probability of transitioning from x to y. For ease of notation, we allow repetitions in the set notation of Yx, and assume that data points in Yx are indeed distributed according to Px(y) and are sufficient for estimating Px(y). Further, we assume that these transitions are smooth in the sense that if x is similar to x′ then P (y | x) is similar to P (y | x′). Therefore in practice we may replace Yx in the following
Figure 1: (A) Schema of DyMoN architecture where x represents the neural network input vector, ŷ represents the predicted network output in (1), and is a random Gaussian noise vector. P (y|x) represents the distribution of outputs from many iterations of the stochastic dynamical process given x, and P (ŷ|x) represents the distribution of outputs of DyMoN given x and many different noise vectors. (B) Example input states. (C) Learned transition vectors (arrows) and output states from DyMoN.
training procedure with ∪x′≈xYx′ , which includes target points from neighbors of x to increase the robustness of DyMoN.
The stochastic output of DyMoN is enabled by a random input vector ε ∈ Rn sampled from a simple (e.g., normal or uniform) distribution F with zero mean and unit variance. This input can be explicitly written into DyMoN function as
Tε(x, θ) = x+ fθ(x, ε), (1)
where fθ : Rd ×Rn → Rd represents suitable feed forward layers for combining the training input x with the random input ε to provide an instantiation of the transition velocity vector ∆θ(x).
Given x ∈ X and θ ∈ Θ, one can consider the distribution P̂(θ)x of the random variable Tε(x, θ), ε ∼ F and estimate it by m i.i.d. instantiations = {εj ∼ F}mj=1 passed through the network to form Ŷ ( )(x,θ) = {Tε1(x, θ), . . . , Tεm(x, θ)}. To conform with the training data, DyMoN is optimized so that this distribution approximates the distribution Px, as captured from the training data by Yx. This optimization is given by arg minθ E [{Hx(θ) : x ∈ X ] with Hx(θ) = MMD(P̂ (θ) x ,Px) (see supplemental material for details on MMD). The MMD is computed using Ŷ ( )(x,θ) and Yx. Given trained weights θ ∈ Θ and an initial state x0 ∈ X , a random walk is generated by xt = Tεt(xt−1, θ), where t = 1, 2, . . ., and ε1, ε2, . . .
i.i.d.∼ F . We note that DyMoN can be used to learn deterministic time series by considering Dirac conditional probabilities in the Markov process, essentially ignoring the random input. Additionally, DyMoN can be extended to encode high-order Markov processes, which have fixed-size memory rather than being memoryless. To this end, an n-th order DyMoN takes n input states and have the distribution of ∆θ(xt−1, xt−2, . . . , xt−n) determine the conditional probability P (xt|xt−1, xt−2, . . . , xt−n) of an n-th order Markov process. The training and application of both a deterministic and a high-order DyMoN extend naturally from the stochastic memoryless (i.e., first-order) DyMoN.

heading: 3 BIOLOGICAL CASE STUDIES

heading: 3.1 BIOLOGICAL CASE STUDY 1: CALCIUM IMAGING OF VISUAL CORTEX NEURONS

text:
We apply DyMoN to dynamic calcium imaging data from 33 Somatostatin-expressing GABAergic (SOM) interneurons from a mouse visual cortex (Urban-Ciecko & Barth, 2016). SOM interneurons are a unique set of GABAergic cells and represent about 25% of all cortical interneurons. and are characterized by high basal firing activity. They are also thought to be connected in dense networks, though the architecture of these networks is not well-known. The data consists of a total of 129,500 frames at 30Hz, lasting 72 minutes, each containing a recording of the change of fluorescence (dF/F) indicating calcium activity over time . We train DyMoN to predict the the next state (which we define as 10 frames later) given the current state of the 33 neurons.
Surprisingly, though SOM cells are not connected to each other by chemical or electrical synapses and do not receive strong feedforward input, we see co-activation relationships in the data. We see in Figure 2A that neurons 9 and 18 are strongly affected by two large clusters of neurons in opposite directions. There is a group of nine neurons which activate 9 and 18, with neuron 19 being the strongest activator. A second group of 9 neurons which inhibit neurons 9 and 18, with neuron 4 being
the strongest inhibitor. We speculate that these interactions are the result of either 1. co-activation by lateral excitatory inputs from superficial pyramidal neurons, 2. co-modulation by neuromodulators, or 3. top-down feedback input regulating the activity of a specific subset of SOM cells (Yavorska & Wehr, 2016).
A subset of the SOM cells appear to function as hubs in an excitatory pattern (see Figure 2B). Early in postnatal development, a subset of inhibitory neurons function as hub neurons and regulate the activity of surrounding neurons. Hub neurons have mainly been characterized in the hippocampus and entorhinal cortex, and a hub role for sensory cortical interneurons would be highly novel (Bonifazi et al., 2009). Our future work involves examining whether these results hold in additional mouse samples.

heading: 3.2 BIOLOGICAL CASE STUDY 2: T CELL DEVELOPMENT IN THE THYMUS

text:
Next, we use DyMoN to learn transitions in single-cell data. Single-cell data has recently gained popularity in biology as a way of dissecting cellular heterogeneity. Recent works such as Diffusion Pseudotime (Haghverdi et al., 2016), Wanderlust (Bendall et al., 2014), and Wishbone (Setty et al., 2016) take the view that cells from a single experimental sample are at different stages of development and derive a pseudo-temporal developmental ordering to the cells. The idea of pseudo-time is that while we cannot follow a single cell through time with current high-dimensional single cell technologies, the entire population of cells can be used to derive potential cell trajectories within a sample.
Here, we use a mass cytometry dataset measuring actively developing T cells, which are adaptive immune cells, obtained from a mouse thymus (Setty et al., 2016). We preprocess this data using the MAGIC algorithm (van Dijk et al., 2018). After preprocessing, we are left with 17,000 cells, measured in 33 protein dimensions (Figure 3). We then train DyMoN on a trajectory obtained by using a diffusion-based Markov affinity matrix (as described in Moon et al. (2017)) and the ordering obtained by Wishbone.
After training we sample trajectories by initializing DyMoN with undifferentiated cells. We obtain two types of trajectories. To investigate the difference between the two trajectories we examine their marker expression as a function of trajectory progression (Figure 3ii, iii). Both trajectories start out with low CD4 and CD8, which is expected form naive T cells. Both trajectories then increase in CD4 and CD8, signaling the change to the double-positive T cells. Then, at this point, we find that the two trajectories diverge with the main difference between the two trajectories that in one CD8 goes down whereas in the other CD4 goes down. This signifies the transitions into T helper cells (CD4+/CD8-) or cytotoxic T cells (CD4-/CD8+) respectively.
To investigate the internal representation that is learned by DyMoN, we compute the Jacobian of the transition vector with respect to the inputs, at two different points on the trajectory: at the branch point (Figure 3Di) and at the CD8+/CD4- branch (Figure 3Dii). DyMoN is capable of learning
different associations between genes at different points in the ambient space. We find a negative association between CD4 and CD8 in the Jacobian obtained at the branch point, reflecting the decision between downregulating either CD4 or CD8. In addition, we find that in the Jacobian of the CD8+ branch Gata3, a developmental marker, goes down with CD4 and CD8 confirming the developmental characteristic of the trajectory. Finally, in the same CD8+ branch Jacobian, Foxp3 and CD25 are positively associated, which is consistent with their role as regulatory T cell markers. The concordance of our analysis with the established literature on T cell development provides validation that DyMoN can both learn meaningful trajectories, and can be further interrogated to provide insights into the network’s dynamical model of the system.

heading: 3.3 BIOLOGICAL CASE STUDY 3: HUMAN EMBRYONIC STEM CELL DIFFERENTIATION

text:
In order to apply DyMoN’s generative capabilities on a biologically novel system, we generated single-cell RNA sequencing data of human embryonic stem cells (hESCs) differentiating in an embryoid body (EB) system. EB differentiation has been shown to mirror early human development in numerous cell types (e.g. neuronal, hematopoietic, muscle, pancreatic cells, etc), and is used as a model to develop cellular reprogramming protocols (Bibel et al., 2007; Nakano et al., 1996; Kania et al., 2004; Geijsen et al., 2004).
To investigate potentially novel reprogramming strategies in the EB system with DyMoN, we measured approximately 31,000 cells, using the 10x Chromium platform, equally distributed over a 27-day differentiation time course, and preprocessed the cells with MAGIC (van Dijk et al., 2018). Because this dataset has numerous lineages, it is not amenable to training by diffusion pseudotime. Instead, we train the DyMoN on Markovian transitions based on the diffusion geometry of the dataset. We sample neighbors y of x over weighted affinities defined by a Gaussian kernel, retaining only those neighbors for which y is defined to be ”later” than x, where time is defined here as a smoothed estimate of the discrete time variable denoting development within the experimental time course.
In Figure 4, we show two representative trajectories generated through the data by DyMoN, one of which samples neural progenitor development (A), and the other bone progenitor development (B), both starting out as embryonic stem cells. We show in a heatmap of genes with high mutual information with the path position that each of these paths initially express known stem cell markers
NANOG and POU5F1, and follow a common differentiation transition into the ectoderm (LHX2). Following this, the two paths diverge, with path A taking the neural progenitor branch distinguished by SOX1, and path B taking the bone progenitor branch distinguished by ALPL. We observe a novel set of transcription factors distinguishing each stage of differentiation in these generated trajectories (e.g. RBM17, FNDC3B, HAND1, and HES5, see Figure 4iii), and propose these transcription factors for a potentially novel reprogramming protocol to obtain each respective mature cell type. In this system, DyMoN provides a novel form of hypothesis generation enabled by the deep abstract model of the differentiation process, enabling researchers to perform more informed validation experiments based on hypotheses generated in an unsupervised and data-driven manner.

heading: 4 EMPIRICAL VALIDATION

text:
In this section, we demonstrate the performance and accuracy of DyMoN on a variety of datasets including simulated datasets and the Frey faces dataset (Roweis & Saul, 2000). We show that DyMoN is able to learn processes that are harmonic or chaotic, ergodic or null-recurrent, and compare its performance to a range of other methods used to perform prediction in dynamical systems. We further show that DyMoN’s output and hidden layers can be interrogated to learn about the process being modeled. Specifically, we show that we can reliably use DyMoN to 1. learn the transition probabilities at each state; 2. generate stochastic trajectories within a system; 3. visualize the data; and 4. learn the relevant features that drive the transitions. For training details, see the supplementary material.
Learning transition probabilities: We first demonstrate DyMoN on two test cases: a circle, which is traversed directionally and non-directionally; and a tree, which is traversed directionally starting from the root. The tree is generated by creating points sequentially along branches at random angles with one another, and then adding Gaussian noise. In both cases, the inputs to the network are two-dimensional Euclidean coordinates. Training examples are drawn from the tree where transitions crossing branch points select a branch with probability 0.5.
Figure 5 shows vector fields for these two examples indicating the distribution of transitions at each state and at unobserved states as predicted by DyMoN. DyMoN is able to randomly select a direction of travel at each branch point, allowing it to walk along a data manifold in the correct direction. Empirical testing on these test systems also motivates our choice to include a skip connection from
the inputs to the outputs, allowing the network to learn the velocity of the transition function rather than the position transition itself (shown in the supplemental material).
Trajectory generation: We now demonstrate the ability of an nth order DyMoN to generate paths on a single and double pendulum. Since these processes are deterministic, we do not add noise to the input of DyMoN and use the mean squared error as the loss function. For the single pendulum, we use a second-order DyMoN with the current and previous angle of the pendulum as inputs. For the double pendulum, we use a third-order DyMoN with the current and last two angles of both pendulums as inputs.
Figure 6 shows the Euclidean coordinates of DyMoN-generated paths of both pendulums over time, with only the second pendulum shown in the case of the double pendulum. Both predicted pendulums show smooth trajectories, with the single pendulum showing periodic behavior and the double pendulum showing chaotic behavior.
We also trained DyMoN on the Frey faces dataset (Roweis & Saul, 2000) in order to demonstrate DyMoN’s capacity to learn an empirical model of a stochastic system for which no generating distribution exists. Figure 7 shows 1000 samples generated by DyMoN on the Frey faces dataset visualized using PCA. The first 1000 samples are discarded. Ten generated samples with uniform spacing in time are also shown. DyMoN samples a large range of states and generates a realistic new trajectory.
Stationary distribution comparison to other methods: We train the DyMoN to sample a Gaussian mixture model in order to show that DyMoN reliably learns to generate samples of the system which match the expected data distribution, thus avoiding the ”mode collapse“ problem which affects Generative Adversarial Networks. We generate training data from a 1-dimensional Gaussian mixture model using Metropolis-Hastings sampling, and compare the performance of Recurrent Neural Networks (RNNs), Hidden Markov Models (HMMs) and Kalman Filters (KFs) on this simple stochastic system. Figure 8 shows 50,000 points generated in a chain from each method, and we compare the Earth Mover’s Distance (EMD) of these chains to a distribution of 200,000 points
generated by MCMC in Table 5. DyMoN provides performance on par with the HMM and MCMC, both in terms of accuracy of sampling the distribution (EMD) and inference time. RNNs fail to capture the stochasticity of the system, falling into an infinite loop reproducing the same data point. Kalman filters are inherently difficult to train on a system without predefined states and transitions, and as such under-perform in an unsupervised setting. Additionally, although the marginal distribution is similar, the dynamics of the HMM fail to capture the full variability of the MCMC sampling.
Pointwise sample comparison to other methods: We train the a DyMoN with 2 convolutional layers with 5x5 filters and 2x2 max pooling, a single fully connected layer of 3 hidden nodes, and 2 deconvolutional layers on a video of a rotating teapot (Weinberger et al., 2004). We compare the performance of RNNs, HMMs and KFs to DyMoN on this higher-dimensional example. We produce a single transition from each frame in the time series (giving a chain of prior states to the HMM, RNN and KF in order to facilitate their estimates of the current state) and measure the mean squared error between the produced output and the true frame, ten time points after the current state. Due to the relatively little training data available (400 samples), both classical statistical learning methods (HMM, KF) have difficulty generating accurate samples from this system given a single snapshot. RNNs perform well on the large majority of frames when given a sequence as input, but fail entirely when given only a single frame as input on inference, making them inherently unsuitable to inference in snapshot biological systems where sequential input is not available. Additionally, DyMoN is able to generate images faster than all other methods, and is only slower in training than the HMM, which produces significantly lower quality images. Examples of images produced by each method are shown in Figure 11.
Visualization: A DyMoN with a low-dimensional latent layer can also be used to produce a visualizable embedding of the data. Figure 9 shows the embedding layer of the DyMoN trained on the teapot data produces a single trajectory homeomorphic to a circle, while the PCA embedding of the same dataset produces spurious intersections and branches in the data. We see that the RNN trained
with a similar architecture (a 64 node LSTM layer is added to both sides of the 3 node layer) does not contain an interpretable embedding, as the RNN does not simply model the transition and has memory that may confound the embedding.

heading: 5 CONCLUSION

text: Here we presented a framework, which we call DyMoN, for designing neural networks that can model stochastic dynamics from biological data. The DyMoN framework is well-suited for creating a deep model of n-th order Markovian stochastic dynamics from the types of data that occur in biological settings, especially in systems for which the generative process cannot be described by differential equations. In addition to actual longitudinal samples of dynamics, this can include noisy or sparse data, snapshot data where neighborhoods represent probabilistic next states, pseudotime trajectories created over data, or even trained Markov models. The flexibility of this framework is enabled by several aspects of DyMoN. First, since the networks encode Markovian dynamics, they do not require long histories that are difficult to obtain in biological data. Secondly, the MMD penalty used to train the next-state generation, can enforced via samples or a known probability distribution, making the networks trainable on top of other more shallow models such as diffusion operators and pseudotime inferences. We show that creating a deep model of biological stochastic dynamics enables several different analyses of such systems, including the generation of new trajectories to augment the data, extraction of feature dependencies of the dynamics, and visualization of the dynamic process, as well as denoising trajectories from naturally noisy data. We believe that DyMoN will enable inference of driving forces (genes, transcription factors, mutations, etc) of progression dynamics from the large amounts of biomedical data now being generated in many settings.

heading: A BACKGROUND

text:
Markov Processes: Markov processes are stochastic memoryless models that describe dynamic systems, i.e., each state depends only on the last state. However, higher order processes can be modeled by making the state dependent on several previous states. A deterministic Markov chain maps each state to only one next state whereas a stochastic Markov chain models a probabilistic transition such that at each state multiple transitions are possible. DyMoN is a neural network trained to replicate a stochastic (with input noise) or deterministic (no input noise) Markov process.
Maximum Mean Discrepancy: Divergences, such as KL divergence, are used to measure distances between probability distributions. However, they require density estimation and as such divergences are hard to compute for high dimensional systems. Maximum Mean Discrepancy (MMD) (Gretton et al., 2012) offers a solution by using the kernel trick to circumvent the curse of high dimensionality. Instead of comparing empirical distributions, MMD is defined on the inter- and intra-sample pairwise affinities, which are computed using a kernel. With MMD we can therefore compute the distributional distance (i.e. divergence) between two samples, without the need for density estimation. This distance is defined as
MMD2(µ, ν) = ∫∫ k(x, x′)dµdµ+ ∫∫ k(y, y′)dνdν − 2 ∫∫ k(x, y)dµdν
where k(·, ·) is a kernel function and x and y are sampled from the two distributions represented by the probability measures µ and ν. In practice, this distance is estimated using summation over two finite sets of samples.
MMD has been used in Dziugaite et al. (2015) to translate samples from one distribution into another distribution using a deep neural network architecture. Samples were generated by sampling random values from some simple distribution, e.g. Gaussian, and running them through a deep neural network. The network was trained by minimizing the MMD between generated samples and real samples. In Shaham et al. (2017) this architecture was extended by using a residual network with the application of removing undesirable batch effects that are associated with measurements. We use MMD to learn a conditional distribution that represents the possible next states in the Markov chain.
Calcium Imaging Data (Urban-Ciecko & Barth, 2016) Details: To express GCaMP in SOM+ cells, the mouse was anesthetized with 1- 2% isoflurane mixed with pure oxygen, and we injected the two 100nL injections of adenoassociated virus AAVdj-EF1a-fDIO-GCaMP6m in primary visual cortex (V1) at a rate of 75nL/min. For implantation of the imaging window, the mouse was anesthetized using a mixture of ketamine (100 mg/kg) and xylazine (10 mg/kg), and a 3 mm diameter craniotomy was opened over V1, where we inserted and fixed a small rectangular glass piece attached to a 5mm circular cover glass into the craniotomy. A custom titanium head post was secured to the skull with Metabond. Two weeks after implantation, mice were placed on the wheel and head-fixed under the microscope objective. Imaging was performed using a resonant scanner-based two-photon microscope (MOM, Sutter Instruments) coupled to a Ti:Sapphire laser (MaiTai DeepSee, Spectra Physics) tuned to 920 nm for GCaMP6. Images were acquired using ScanImage at 30 Hz, 512x512 pixels (580x580 um). Imaging of layer 2/3 was performed at 150-300 um depth relative to the brain surface. Mice ran freely on the wheel, and after 10 mins of recording, were presented with visual stimuli (drifting gratings with vary contrasts) for 50 mins (40 presentations of each contrast, contrasts ranging from 0%, i.e. gray screen, to 100%, steps of 10%). After the stimulation session was done, we recorded for another 10 mins of spontaneous activity. Analysis of imaging data was performed using ImageJ and custom routines in MATLAB. Motion artifacts and drifts in the Ca 2+ signal were corrected with the moco plug-in in ImageJ (Dubbs et al., 2016), and regions of interest (ROIs) were selected as described in (Chen et al., 2013). All pixels in a given ROI were averaged as a measure of fluorescence, and the neuropil signal was subtracted.

heading: B RELATED WORK

text:
Recurrent Networks (RNNs): RNNs are trained to predict the next state of a sequence and have been used for text analysis (Fernández et al., 2007), speech recognition (Amodei et al., 2016), and language modeling (Józefowicz et al., 2016), as well as other tasks that operate on time series data. In contrast to DyMoN, RNNs typically require a sequence or history to predict the next state.
Stochastic Generative Networks: There have been other networks that use stochasticity to learn conditional distributions. For example, given an input X and a stochastic input to a middle layer, a variational autoencoder (VAE) (Kingma & Welling, 2014) learns to transform these inputs into a Gaussian centered at a maximum likelihood point from which X is derived. Thus the variational network can denoise samples and also generate samples “like” given samples. However, VAEs focus on generating data points instead of transitions. Several stochastic RNNs have also been proposed (Bowman et al., 2016; Chung et al., 2015). For example, in Goyal et al. (2017), a stochastic generative RNN is trained to learn a process that converges to the full data distribution within a small number of steps given a simple initialization. In contrast, DyMoN learns the transitional probabilities from state to state of a Markov process instead of its stationary distribution.
Time Delay Neural Networks: DyMoN also has some conceptual connection to Time Delay Neural Networks (TDNNs) which include a contextual window samples as input. However, the main focus of TDNNs are to classify patterns with shift-invariance, such as recognizing phenomes in speech (Waibel et al., 1989). TDNNs generally achieve this by taking a time convolution through windows of time to train a classifier. In contrast, we focus on learning Markov processes with the goal of generating plausible next-states and analyzing the dynamics that drive these transitions.

heading: C DYMON DETAILS

text:
DyMoN architecture: One could envision many neural network architectures which fit the schema of the DyMoN that we propose. In this paper, we use three alternative models; these are shown in Figure 10. The first architecture generates transitions on the ambient space, most suitable to data of low dimensionality. The other two architectures generate transitions on the latent space of an autoencoder, which learns a latent space which, due to gradient descent propagating through the encoder and decoder, learns a latent space more suited to generating transitions than either the ambient space or latent spaces learned by other dimensionality reduction methods. Additionally, the generation of points in a sampled trajectory can be passed through the encoder and decoder in order to prevent accumulation of error over the course of many samples.
Training the DyMoN: We train DyMoN using leaky ReLU activations on the hidden nodes and linear activations on the residual output nodes. Stochastic DyMoNs are trained with Gaussian noise inputs, and all DyMoNs are trained with Gaussian corruption noise. To compute the MMD loss, we
use a multi-scale Gaussian kernel (Bousmalis et al., 2016) with 19 bandwidths ranging from 1e−6 to 1e6, evenly spaced on a log scale. The kernel is computed separately for each bandwidth and then MMD is computed on the sum of the kernels. The Adam optimizer (Kingma & Ba, 2015) is used for stochastic gradient descent in all cases.
Empirical validation of the skip connection: We trained a DyMoN on the circle and pendulum datasets both with and without a skip connection. The resulting training losses are shown in Table 3. Including a skip connection reduces the training loss.
Training time: Table 4 shows DyMoN training times for the datasets considered.

heading: D METHODS COMPARISON

text:
Competing methods were trained on the same dataset as for DyMoN, and with similar architecture where possible. We used Hidden Markov Models provided by the hmmlearn Python package and Kalman Filters provided by the pykalman package, in both cases learning all parameters by Expectation-Maximization (EM). The EM algorithm was run with default parameters. Recurrent Neural Networks were trained with the same number of hidden layers and convolutions as the DyMoN, except in the teapot example where we found performance improved by adding an additional two hidden layers of size 64. Performance of the competing methods using Earth Mover’s Distance
(EMD) on the GMM example is shown in Table 5 and examples of generated frames from the teapot data are shown in Figure 11.


summaries: This paper performs unsupervised classification where the number of classes is unknown. The main idea is to use the CycleGAN framework such that one can reconstruct the original image by first moving to latent space that represents another class (via the connector network), then moving back to the original latent space and going back into image space using a generator. Experiments are conducted on MNIST and CIFAR.

summaries: This paper develops an unsupervised classification algorithm using the idea of CycleGAN. Specifically, it constructs a piece-wise linear mapping (the Connector Network) between the discriminator network and the generator network. The learning objective is based on the cycle-consistency loss. Experiments show that it can achieve reasonable loss. This paper addresses an important problem, namely, unsupervised image classification, and may present interesting ideas. However, the paper is not in a good shape for publication in its current form.

abstractText: We propose a novel unsupervised classification method based on graph Laplacian. Unlike the widely used classification method, this architecture does not require the labels of data and the number of classes. Our key idea is to introduce a approximate linear map and a spectral clustering theory on the dimension reduced spaces into generative adversarial networks. Inspired by the human visual recognition system, the proposed framework can classify and also generate images as the human brains do. We build an approximate linear connector network C analogous to the cerebral cortex, between the discriminator D and the generator G. The connector network allows us to estimate the unknown number of classes. Estimating the number of classes is one of the challenging researches in the unsupervised learning, especially in spectral clustering. The proposed method can also classify the images by using the estimated number of classes. Therefore, we define our method as an unsupervised classification method.

title: How auto-encoders could provide credit assignment in deep networks via target propagation

venue: arXiv preprint arXiv:1407.7906,

authors:
Yoshua Bengio
Yann LeCun

title: Scaling learning algorithms towards ai

venue: Large-scale kernel machines,

authors:
Yoshua Bengio
Grégoire Mesnil
Yann Dauphin
Salah Rifai

title: Better mixing via deep representations

venue: In International Conference on Machine Learning,

authors:
David Berthelot
Thomas Schumm
Luke Metz

title: Began: boundary equilibrium generative adversarial networks

venue: arXiv preprint arXiv:1703.10717,

authors:
Deng Cai
Xiaofei He
Jiawei Han

title: Locally consistent concept factorization for document clustering

venue: IEEE Transactions on Knowledge and Data Engineering,

title: Algorithms for manifold learning

venue: Univ. of California at San Diego Tech. Rep,

authors:
Xi Chen
Yan Duan
Rein Houthooft
John Schulman
Ilya Sutskever
Pieter Abbeel

title: Infogan: Interpretable representation learning by information maximizing generative adversarial nets

venue: In Advances in Neural Information Processing Systems,

authors:
Ehsan Elhamifar
René Vidal

title: Sparse subspace clustering

venue: In Computer Vision and Pattern Recognition,

authors:
Robert Geirhos
David HJ Janssen
Heiko H Schütt
Jonas Rauber
Matthias Bethge
Felix A Wichmann

title: Comparing deep neural networks against humans: object recognition when the signal gets weaker

venue: arXiv preprint arXiv:1706.06969,

authors:
Xavier Glorot
Antoine Bordes
Yoshua Bengio

title: Domain adaptation for large-scale sentiment classification: A deep learning approach

venue: In Proceedings of the 28th international conference on machine learning

authors:
Ian Goodfellow
Honglak Lee
Quoc V Le
Andrew Saxe
Andrew Y Ng

title: Measuring invariances in deep networks. In Advances in neural information processing

authors:
Ian Goodfellow
Jean Pouget-Abadie
Mehdi Mirza
Bing Xu
David Warde-Farley
Sherjil Ozair
Aaron Courville
Yoshua Bengio

title: Generative adversarial nets

venue: In Advances in neural information processing systems,

authors:
Kaiming He
Xiangyu Zhang
Shaoqing Ren
Jian Sun

title: Deep residual learning for image recognition

venue: In Proceedings of the IEEE conference on computer vision and pattern recognition,

authors:
Pan Ji
Tong Zhang
Hongdong Li
Mathieu Salzmann
Ian Reid

title: Deep subspace clustering networks

title: Auto-encoding variational bayes

authors:
Marc T Law
Raquel Urtasun
Richard S Zemel

title: Deep spectral clustering learning

authors:
Mehdi Mirza
Simon Osindero

title: Conditional generative adversarial nets

venue: arXiv preprint arXiv:1411.1784,

authors:
Andrew Y Ng
Michael I Jordan
Yair Weiss

title: On spectral clustering: Analysis and an algorithm

authors:
Alec Radford
Luke Metz
Soumith Chintala

title: Unsupervised representation learning with deep convolutional generative adversarial networks

venue: arXiv preprint arXiv:1511.06434,

authors:
Uri Shaham
Kelly Stanton
Henry Li
Boaz Nadler
Ronen Basri
Yuval Kluger

title: Spectralnet: Spectral clustering using deep neural networks

venue: arXiv preprint arXiv:1801.01587,

authors:
Jianbo Shi
Jitendra Malik

title: Normalized cuts and image segmentation

venue: IEEE Transactions on pattern analysis and machine intelligence,

title: A tutorial on spectral clustering

venue: Statistics and computing,

authors:
Lihi Zelnik-Manor
Pietro Perona

title: Self-tuning spectral clustering

authors:
Junbo Zhao
Michael Mathieu
Yann LeCun

title: Energy-based generative adversarial network

venue: arXiv preprint arXiv:1609.03126,

authors:
Jun-Yan Zhu
Taesung Park
Phillip Isola
Alexei A Efros

title: Unpaired image-to-image translation using cycle-consistent adversarial networks

venue: In Computer Vision (ICCV),

heading: 1 INTRODUCTION

text:
Image classification has been one of the most important research in machine learning and artificial intelligence. For the predefined number k of classes, the classification models based on deep learning have performed well that the test accuracy is better than humans (Geirhos et al., 2017; He et al., 2016; Russakovsky et al., 2015; Szegedy et al., 2015). This result is, however, true only when the number of classes, k is predefined and we use the supervised learning with entire ground truth labels. Humans learn visual information by memorizing the compressed information (features) in the brain, reconstruct the object, and recognize the difference between objects. This capability allows humans to classify data where the number of classes is unknown into the reasonable numbers of classes. In this paper, we propose a method to mimic such human capability.
The main contribution of the proposed framework consists of two parts. We first show that there exists a approximate linear map between two separated neural networks of Generative Adversarial Networks (GANs). Secondly, we prove that the number of classes in a latent space is the same as the number of classes in the original space in terms of the multiplicity of eigenvalue 0 of a Laplacian matrix in spectral clustering theory. By combining the results, We derive a novel unsupervised classification framework, which estimates the number of classes k and then extract k vectors representing each class. With the framework, the approximate linear map between two networks can induce the information of connected components in a data graph.
The latent variable model is commonly used in the generative models. Variational auto-encoders (VAEs) (Kingma & Welling, 2013) and GANs (Goodfellow et al., 2014) are a few of the most famous cases that learn the data distribution successfully. GANs randomly samples the latent variables that are used to train the data distribution so we cannot figure out which latent variable generates a result that we are interested in. To avoid the drawback by random sampling, auxiliary conditions on latent variables were introduced (Chen et al., 2016; Mirza & Osindero, 2014). Introducing a label information into the condition of the estimated data probability is highly effective in classification where a data point x came from, however this method is not a unsupervised learning.
Spectral clustering is one of the most popular clustering algorithm that outperforms other wellknown clustering algorithms (Ng et al., 2002; Shi & Malik, 2000). Despite high performance, spectral clustering is difficult to use because it requires the full dataset that have very largeN data points. In this case, the dimension of the Laplacian matrix is N ×N and it takes at least O(N3). Therefore,
we need a novel dimension reduction technique to use spectral clustering in practical applications (Elhamifar & Vidal, 2009; Ji et al., 2017; Law et al., 2017; Shaham et al., 2018).
Recent image classification methods are based on supervised learning, in which a learning model processes a dataset consisting of pairs of an image and its ground-truth label. This training dataset usually consists of more than 10 K data point, some part of which need to be used for crossvalidation. It also takes a lot of costs to produce a reliable training dataset. In addition, estimating the number of classes is a key aspect of unsupervised learning. There have been a few heuristic method (Zelnik-Manor & Perona, 2005; Von Luxburg, 2007). We propose a method to find the number of classes analytically.
The main concepts of spectral clustering, GANs and the self-expressiveness property will be introduced Section 2. We prove the preliminary propositions that allow us to classify data with unsupervised learning methods in Section 3. Finally, we propose a unsupervised classification algorithm and test our algorithm in Section 4.

heading: 2 PRELIMINARIES

text: We first review cycle-consistent adversarial networks and sparse subspace clustering that we will use when extracting indicator vectors. We also overview the basic probabilistic property between a classification model and a generative model, which can be used as tools to generate indicator vectors.

heading: 2.1 CYCLE-CONSISTENT ADVERSARIAL NETWORKS

text:
Generative Adversarial Networks (GANs) is a generative model that a generator network learns a data distribution implicitly (Goodfellow et al., 2014). GANs is a two player minimax game with two neural networks: a discriminator D and a generator G. The discriminator D(x) is a probability measure that an input x follows the data distribution pdata(x) rather than the generator’s distribution pg . The generator G is a mapping from the latent space z ∼ pz to the data space. This minimax game can be represented as a stochastic optimization problem with value function given by
min G max D V (D,G) = Ex∼pdata(x) [logD(x)] + Ez∼pz(z) [log (1−D (G (z)))] . (1)
The cycle-consistent adversarial networks (cycleGANs) is a variant of GANs, which learns how to translate an representation from a source domain S to a target domain T without a paired dataset (Zhu et al., 2017). The cycle-consistency states that two mappings G : S → T and F : T → S encode the right permutation of images in S and T , so that, for an arbitrary data point s ∈ S, the mapping from S to T and back again, F (G(s)), should be the same as s. This property encourages adversarial networks G and D to have a connection between generated images and input data. We adopt two cycle consistency loss for the forward map and the backward map to learn a connection between two high-level representation domain in G and D.

heading: 2.2 SPARSE SUBSPACE CLUSTERING

text:
Sparse subspace clustering is a clustering algorithm that is based on spectral clustering. Spectral clustering has many variants depending on the similarity measure, Gaussian kernel similarity function (Ng et al., 2002), k-nearest neighbor graphs, -neighborhood graph. A similarity measure defined on the similarity graph G = (V,E) is a pairwise metric between two data points (nodes) in a dataset and the similarity wij of two data points xi, xj is equal to or greater than 0. The similarity wij is also same as wji because the above similarity graph is defined as an undirected graph. The similarity matrix W is composed of the elements which correspond to the similarity measure wij respectively and then W is symmetric. See more details about spectral clustering in (Von Luxburg, 2007).
Sparse subspace clustering (Elhamifar & Vidal, 2009) uses a novel similarity matrix that is based on the self-expressiveness property. When we assume that the data space is a union of linear subspaces S = ∪kl=1Sl, we can express an arbitrary data point in the dataset as a linear combination of other data points in the dataset, where the coefficients of the data points in the same subspace are non-zero and others should be 0 (Elhamifar & Vidal, 2009; Ji et al., 2017). Formally, each data point yi ∈ S can be reconstructed as yi = ∑N j=1 cijyj = Yci, where cii = 0, N is the size of the dataset,
ci := [ci1 ci2 · · · ciN ]T is the coefficient vector of yi and Y = [y1 y2 · · · yN ] is the dictionary matrix. The sparse solution of ci indicates that the node i is in the same subspace as the node j if the element cij is non-zero. In this paper, we expand the notion of the self-expressiveness property to reduce the dimensionality of L in Section 3.3.

heading: 3 GRAPH LAPLACIAN OF THE CONNECTOR NETWORK

text: We first introduce our proposed architecture of adversarial networks. As an autoencoder point of view, we can consider two adversarial networks as autoencoder. We build a connector network between two networks to improve cycle consistency, and extract the information of connected components (as shown in Figure 1-b). By training with the cycle-consistency loss, we can make the tuples (Z,X,H) for each data point X in a given dataset. Then, we show the relation between the spectral clustering and unsupervised classification. In order to show the relationship, we prove the dimension reduction of the mutual expressiveness property in Section 3.3.

heading: 3.1 DISENTANGLING UNFOLDS AND EXPANDING THE VOLUME OF MANIFOLD

text:
Deep neural networks learn a high-level abstract representation of a data manifold and a well trained abstraction can disentangle the underlying factors of variation (Bengio et al., 2007). These properties allow interpolation between data samples more naturally and smoothly. (Bengio et al., 2013) studied the disentangling effect and proposed three hypotheses on the shape of manifolds in different representation domains. These hypotheses are verified by several empirical results (Goodfellow et al., 2009; Glorot et al., 2011). The three hypotheses are summerized as follows.
• The deeper representations can better disentangle the underlying factors of variation. • Disentangled representations unfold the manifolds near which raw data concentrates ir-
regularly, and expand the relative volume occupied by high-probability points near these manifolds. Then these representations are with greater convexity.
• The underlying class factors of variations will be better disentangled than other factors, so the deeper layer has better discriminant capability of classes.
In our proposed architecture, there are two high-level (deeper) representation spaces which are called as feature spaces. These feature spaces fill the space more uniformly than the pixel space X . H and
Z also have smoother density. Therefore, H and Z have greater convexity than X and it permits linear combination on these two manifolds in the feature spaces H and Z, are the two manifolds are more convex than the manifold in the pixel space X . By manifold hypothesis, different class manifolds are well-separated by regions of very low density (Cayton, 2005). The feature spaces are high-level representation spaces so that the density is unfolded and a linear combination of features generates amore natural result in the pixel space Bengio (2014).
Now, we can define the effective basis of the subspace for each class. In the feature spaces, each class manifold is assumed to have the enough convexity. The class manifolds are well-separated, and show good discriminant capability. Using both properties, we can build the effective basis for each class subspaces and generate the original space by a direct sum of these subspaces, where the subspaces are orthogonal. We empirically tested an effective basis as shown in figure 3. Randomly sampled features in a subspace of the data space can be an effective basis, and are possible to generate all data points in the same subspace. The features in a different subspace cannot generate the given data samples. We can estimate the effective dimension with the number of the effective basis. We then can generate orthogonal subspaces of the entire classes if the effective dimension is small enough than the dimension of feature spaces.

heading: 3.2 CONNECTOR NETWORK AND CYCLE CONSISTENCY

text:
We will study the relationship among 3 different spaces, Z, X and H in the proposed architecture. As shown in figure 1-(b), we can consider two adversarial networks as a decoder and an encoder of an auto-encoder respectively (Berthelot et al., 2017; Zhao et al., 2016). The discriminator network D encodes the pixel space X to the feature space H , and the generator network G decodes the feature space Z to X . We define the connector network C which maps from H to Z. The map C between two feature spaces will be trained by the cycle consistency loss to obtain the tuples (Z,X,H) with the correct permutation, where all the elements are in the same-class manifold and shares the same learned features. Cycle consistency loss tends to minimize the difference between an arbitrary point and the returned point after traveling other representation domains. This cycleconsistency objective permits inverse maps between all representation spaces. When we apply the forward cycle-consistency loss ‖z − C(D(G(z)))‖1 and the backward cycle consistency loss ‖x− G(C(D(x)))‖1, the maps between all pairs of neighboring spaces become unique. Therefore, the feature-invariant one-to-one maps among G,C,D can be induced. All elements of an arbitrary tuple (Z,X,H) share same class features as shown In figure 1-(a). We now have the following assumption.
Assumption 3.1. All elements of the cycle consistent tuple (Z,X,H) shares the identical factors, such as class factor and other factors (e.g. rotation and brightness).
As shown in fig 3-(c, d), we can conjecture that each class subset is closed and the linear combinations in the certain class subset are also in the class because each class manifold is unfolded and expanded in Z and H . In addition, we showed there exists the effective basis for each class. By
the hypotheses (Bengio et al., 2013), we know that each class manifold has the compactness and the convexity. Now we can assume the following assumption. Assumption 3.2. Each class subspace Hi, where Hi ⊂ H are orthogonal, so as Zi where Zi ⊂ Z. We define sets for each class SHi and SZi in Hi, Zi respectively to be the the class manifolds which lies on the corresponding class subspaces. The sets SHi and SZi are compact and convex.
This assumption is more tight and realistic than the union of linear subspaces in Elhamifar & Vidal (2009). We can also conjecture that the feature space H and Z are descried as the direct sum of Hi and Zi respectively. Now, we can show the existence of an approximate linear map C from H to Z. Lemma 3.1. There exists a compact and convex linear map Ci from SHi to SZi , If ‖Ci‖ is closed.
Proof.
1. By the assumption 3.2, the class subset SHi and SZi are convex. Let h1, h2 ∈ SHi and z1 = Cih1, z2 = Cih2. We know that Ci is linear, then
tz1 + (1− t)z2 = tCih1 + (1− t)Cih2 = Ci(th1 + (1− t)h2). (2)
th1 + (1− t)h2 is also in SHi for t ∈ [0, 1]. Also, tz1 + (1− t)z2 is in SZi because SZi is convex. Therefore Ci is a convex linear map for all i.
2. If SHi is compact, then for h ∈ SHi , ‖h‖ is closed. Then z ∈ SZi is closed, because
‖z‖ = ‖Cih‖ ≤ ‖Ci‖‖h‖. (3)
We know that Ci and h are finite. Therefore SZi is also compact. Theorem 3.1. By the assumption that the data space are the direct sum of orthogonal class subspaces, the map from H to Z can be written as a block matrix, of which block component consists of Ci. C can maps from all class subset SHi to SZi .
Proof. By Lemma. 3.1., Ci is a compact and convex linear map from Hi to Zi. H and Z are the direct sum of the orthogonal subspaces Hi and Zi respectively. Therefore C can be represented as a block diagonal matrix for some ordered basis.

heading: 3.3 MUTUAL EXPRESSIVENESS PROPERTY AND DIMENSION REDUCTION

text:
As we discussed in Section 2, the self-expressiveness property shows that we can express an arbitrary data point as a linear combination of other N data points. By modifying this property slightly, we can write each data point x as a linear combination ofN sampled data points α = {x1,x2, · · · ,xN} from the data space X . This fact allows us to generate a set of feature variables, which satisfies the self-expressiveness property. The self-expressiveness property implies that the generated data points are in the class manifolds. Now, we relax the self-expressiveness property to generalize the problem. We divide the latent space H into H1, · · · , Hl where the self-expressiveness only satisfies in each
Hi not in H − Hi. We also know that it is possible to discard redundant coefficients in the selfexpressiveness expression by optimizing the matrix of coefficients (Elhamifar & Vidal, 2009), so there exists at least one sparse solution of x where most coefficients are zero. By observing this property, it is expected that the number of non zero coefficients approaches the dimension of latent space, m. That is because N is much higher than dimX , and X is a m-dimensional manifold. Image clustering methods, for instance, treat a dataset where the number of data points is much larger than the dimension of data points so we can optimize a solution to have sufficiently many zero coefficients in the self-expressive linear combination. Theorem 3.2. From now on, Z ′i = C ′iH ′i where Z ′i = (z1, · · · , zmi)T and H ′i = (h1, · · · , hmi)T .
Proof. By an observation on GANs, we have Ni feature variables βi = {hij}1≤j≤Ni for each 1 ≤ i ≤ l where D(xi) = hi and N = N1 + · · ·+Nl. By the self-expressiveness property in Hi, we can express each hi ∈ Hi by the linear combination of β′i = {h1, h2 · · · , hmi , hmi+1, · · ·hNi} where β′i is a reordered set of βi and the coefficents of the linear combination {cmi+1, cmi+2, · · · , cNi} are zero. In other words, the ordered set γi = {h1, · · · , hmi} is enough to represent all points in the class manifold Xi and Zi, which correspond to (Zi, Xi, Hi), i.e. span(γi) = Hi. Hence, we can choose C ′i by an mi ×mi real matrix. This theorem is empirically shown by the our effective basis conjecture. Now, we can conclude that the dimZi = dimHi = mi Ni. We define the connector network C as a ( ∑ mi) × ( ∑ mi) matrix. The similarity matrix is denoted by W = |C|+ |CT |. Theorem 3.3. Let G be an undirected graph with non-negative weights. Then, the number of connected components A1, · · · , Ak in G is the same as the multiplicity k of the eigenvalue 0 of L. If W is a k block diagonal matrix, then the multiplicity of the eigenvalue 0 of its Laplacian is also k.
Proof. the proof in (Von Luxburg, 2007). By the definition of W ,the number of sub-block matrices in W is the same as the number of connected components.
We conclude that the numbers of classes are identical between spectral clustering and our method.

heading: 3.4 UNSUPERVISED CLASSIFICATION

text:
Now, we consider a unsupervised classification problem where there is no ground-truth label. For example, when the cross-entropy loss is applied in the supervised classification, we need pairs of input data and ground-truth probability mass p = [0, , 1, · · · , 0], which contains a single 1 at the ith element, i.e. the class index c = i, to estimate the difference between the prediction and the ground-truth. We can formally write the probabilistic interpretation as
P (c = i|x) = P (c = i)P (x|c = i) P (x)
(4)
, where x is a data point. If we assume that the numbers of data points in the whole classes are same, then P (c) is uniform. Then, for the given point x, P (x) is fixed so the we can measure P (c|x) by the estimating the likelihood P (x|c). In Section 3.3, we showed that the indicator vectors of the connected components are the eigenvectors of the Laplacian matrix L, which are corresponding to the eigenvalue 0. The components of the indicator vectors span the eigenspaces of the matched connected components respectively. The orders of components for the Laplacian matrix LC and the connector network C are identical by the definition in Section 3.3. Therefore, for all connected component Ai, the indicator vectors 1Ai = h ∈ Hi spans the class subspace Hi. The indicator vector itself, however, might not be in the class manifold because it can exist outside the boundary of the manifold. We find out that the changes of norms between H and Z are different by the class. The connector network can be considered as a singular value decomposition for each class. This implies that each class data points move differently by the class. We define the mean vector of a class in H as h′ and the Gaussian noise n.
‖z‖2 = ‖Ch‖2 = hTCTCh (5) = (h′ + n)TCTC(h′ + n) (6)
The expectation of ‖z‖2 is, E[‖z‖2] = h′TCTCh′ + E[nTCTCn] (7)
= h′TCTCh′ + E[nTZnZ ] (8)
, where nZ is the scaled noise in Z. We now assume that the distribution of ‖z‖2 follows a Gaussian with a variance σnZ . Then, we can calculate the probability that ‖z‖2 is in the certain class. We also know that B = CTC is a normal operator, then B satisfies the following theorem. Theorem 3.4. |xTBx| ≤ max{|λ| : λ is an eigenvalue of B}‖x‖22
Proof. See in Chapter 5.6 of (Horn et al., 1990).
By ordering the eigenvalue of B, we can find out the interval that ‖z‖2 lies in. This implies that the maximum eigenvalue of the sub-block matrix in B for z is in the above interval. This eigenvalue works as the indicator of classifying class. Therefore we can build of Gaussian distributions for ‖z‖2 with the mean λ‖h‖2 and the standard deviation σnZ .
4 UNSUPERVISED CLASSIFICATION INTO UNKNOWN k CLASSES
In this section, we describe the loss function and the training algorithm of unsupervised classification into unknown k classes. It is known that estimating the number of classes k is hard and heuristic (Von Luxburg, 2007). We find the empirical observation on the hint of k. Using the observation, we suggest a conjecture that helps to estimate the number of classes k in a certain case.
Algorithm 1 Unsupervised Classification into Unknown k Classes (UCUC) 1: for the number of iterations do 2: Sample minibatch of m latent vectors z(0), · · · , z(m) from prior pg(z) 3: Sample minibatch of m data points x(0), · · · ,x(m) from dataset 4: Update the discriminator by ascending :
∇θd 1
m m∑ i=1 [ logD ( x(i) ) + log ( 1−D ( G ( z(i) )))] +λD ∥∥∥x(i) −G(C (D (x(i))))∥∥∥ 1
(9) 5: Update the generator and connect by descending their stochastic gradients:
∇θg 1
m m∑ i=1 [ log ( 1−D ( G ( z(i) )))] + λG ∥∥∥z(i) − C (D (G(z(i))))∥∥∥ 1
(10)
6: Update the connector network C with two cycle-consistency loss ∇C [∥∥∥z(i) − C (D (G(z(i))))∥∥∥ 1 + ∥∥∥x(i) −G(C (D (x(i))))∥∥∥ 1 ] (11)
7: end for 8: Compute the moving average of each ith eigenvalue of normalized symmetric Laplacian matrix Lsym = D − 12LD− 1 2 and find the estimated number k̂ by computing the second local minimum
9: Calculate the eigenvalues of CTC and find a matched class for an input
4.1 ESTIMATING THE NUMBER OF CLASSES k
We found an interesting observation on the number of classes k. The eigenvalue λ of the graph Laplacian matrix equals to Ncut of the given graph (Shi & Malik, 2000). In the ideal cases, we can assume that L is a block-diagonal matrix and each of blocks Li is a proper sub-graph Laplacian. If there is no connection between subgraphs, then L is a perfect block-diagonal matrix and the multiplicity of the eigenvalue 0 is the same as the number of blocks. In the real world dataset, there exists quite large off-diagonal noise in the similarity matrix, as shown in the previous result. These non-zero off-diagonal elements make the eigenvalue zero only at the first eigenvalue and increase the Ncut of the graph gradually. We can consider the ith eigenvalue as the Ncut by the block diagonal matrix with i blocks, which is applied on the perturbed block diagonal matrix of which the underlying connected component is k. We now have the following conjecture on the pattern of eigenvalues. We assume that the off-diagonal elements have the same value and the dimensions of
all eigenspaces are the same. Now, we estimate the number of classes k by starting from m k. We can consider the difference of Ncut between k-block diagonal matrix and m-block diagonal matrix as the the area of Sm − Sm ∩ Sk, where Sm = n 2 m and Sk = n2
k when the matrix size is n. As m increases, the difference of Ncut increases, but at k = m the value have the very small local minimum becuase two matrices are perfectly overlapped. Using this property we can derive the function of the difference of the Ncut over m as
Sm ∩ Sk = n2( 1
k2 (k −m+ 1) + m−1∑ i=1 [( i m − b ik m c1 k )2 + ((1 + b ik m c) 1 k − i m )2], (12)
where m ≤ k. if m > k, and switch m and k. This overlapped area is n 2
k at m = k trivially. The figure of the theoretical conjecture is in Figure-4a. This result can predict the empirical result of the pattern of the eigenvalue. Therefore, we can conjecture that the underlying number of classes is at the second deep local minimum.

heading: 5 EVALUATION

text:
We introduce the setting of experiments and the architecture of our adversarial networks (Radford et al., 2015). For the smaller size of dataset, such as MNIST and CIFAR10, we reduce the size of kernels properly. We have the result the expectation of the estimated number of k is 10.2 for MNIST and 10.1 for CIFAR10. Now we use the unsupervised clustering accuracy (ACC) (Cai et al., 2011; Shaham et al., 2018).
ACC(l, c) = 1
n max π∈Π n∑ i=1 1{li = π(ci)}, (13)
where li is true label, ci is the ground-truth label, and π is the collection of permutations.

heading: 6 CONCLUSION

text: We have proposed an unsupervised classification method that can automatically estimate the unknown number of classes k. Our method provides the learned class-features. The proposed method saves the tremendous cost of producing ground-truth labels in a large dataset.

heading: A SUPPLEMENTARY MATERIALS

text:
A.A MODEL ARCHITECTURE
A.B ADDITIONAL EXPERIMENTAL RESULTS
In the latent space H,Z, the each class subset has a Gaussian distribution, which implies the subset is convex. We observed that the connector networks Ci are not identical among the whole classes. Each connector network Ci is decomposed into an unique singular value decomposition matrices. Therefore, we can extract the difference for each class by analyzing the eigenvalues in Section 3.4.
We can observe thatW is a near block diagonal matrix, of which off-diagonal elements are not zero. This leads to the summation of some residue weights which should be zeros theoretically. When the number of connected components are matched, those residues are no longer summed up. Therefore, the deep local minimum in the eigengap exists.