Online Estimation of Multiple Dynamic Graphs in Pattern Sequences

TL;DR

This paper introduces an online Bayesian method for estimating multiple overlapping graph structures in dynamic binary pattern sequences, improving over traditional models and capturing complex time-dependent interactions.

Contribution

It presents a sequential Bayes algorithm for real-time estimation of multiple graph weights and structures in an Ising-type network model.

Findings

Outperforms traditional orthogonal decomposition methods.

Successfully uncovers overlapping graphs in simulated data.

Effectively models neuronal activity dynamics.

Abstract

Sequences of correlated binary patterns can represent many time-series data including text, movies, and biological signals. These patterns may be described by weighted combinations of a few dominant structures that underpin specific interactions among the binary elements. To extract the dominant correlation structures and their contributions to generating data in a time-dependent manner, we model the dynamics of binary patterns using the state-space model of an Ising-type network that is composed of multiple undirected graphs. We provide a sequential Bayes algorithm to estimate the dynamics of weights on the graphs while gaining the graph structures online. This model can uncover overlapping graphs underlying the data better than a traditional orthogonal decomposition method, and outperforms an original time-dependent Ising model. We assess the performance of the method by simulated data, and demonstrate that spontaneous activity of cultured hippocampal neurons is represented by dynamics of multiple graphs.

Tables

Table 1. TABLE I: Algorithm for the online estimation of the graphs

1.= Set $t=1$ and initialize ${𝐉}_1$.

2.= Obtain the filter density of ${𝜽}_t$, using Eq. 23.

3.= Compute the gradient of the Q-function

===evaluated at ${𝐉}_t$, using Eq. 36.

4.= Obtain ${𝐉}_{t+1}$, using Eq. 33.

5.= Advance in time: $t\to t+1$.

6.= Repeat from step 2.

Equations

$$p(\mathbf{x}|\mathbf{j})=\exp \left[\sum _{i=1}^N{h}_i{x}_i+\sum _{i<j}{j}_{ij}{x}_i{x}_j-\psi (\mathbf{j})\right].$$

$$\mathbf{j}=(h_1,\dots ,h_N,j_{1,2},\dots ,j_{N-1,N}{)}^{\prime }.$$

$$p(\mathbf{x}|\mathbf{j})=\exp \left[{\mathbf{j}}^{\prime }\tilde{\mathbf{F}}(\mathbf{x})-\psi (\mathbf{j})\right],$$

$$\tilde{\mathbf{F}}(\mathbf{x})=(x_1,\dots ,x_N,x_1{x}_2,\dots ,x_{N-1}{x}_N{)}^{\prime }.$$

$$\psi (\mathbf{j})=\log \sum _{\mathbf{x}}\exp [{\mathbf{j}}^{\prime }\tilde{\mathbf{F}}(\mathbf{x})].$$

$$\tilde{\mathbf{\eta }}\equiv \frac{\partial \psi (\mathbf{j})}{\partial \mathbf{j}}=E_{\mathbf{X}|\mathbf{j}}\tilde{\mathbf{F}}(\mathbf{X}),$$

$$\tilde{\mathbf{G}}\equiv \frac{{\partial }^2\psi (\mathbf{j})}{\partial \mathbf{j}\partial {\mathbf{j}}^{\prime }}=E_{\mathbf{X}|\mathbf{j}}\tilde{\mathbf{F}}(\mathbf{X})\tilde{\mathbf{F}}(\mathbf{X}{)}^{\prime }-E_{\mathbf{X}|\mathbf{j}}\tilde{\mathbf{F}}(\mathbf{X})E_{\mathbf{X}|\mathbf{j}}\tilde{\mathbf{F}}(\mathbf{X}{)}^{\prime }.$$

$$p(\mathbf{x}|\mathbf{\theta },\mathbf{J})=\exp \left[\sum _{k=1}^D{\theta }^k\left[\sum _{i=1}^N{h}_i^k{x}_i+\sum _{j>i}{j}_{ij}^k{x}_i{x}_j\right]-\psi (\mathbf{\theta },\mathbf{J})\right],$$

$$p(\mathbf{x}|\mathbf{\theta },\mathbf{J})=\exp \left[{\mathbf{\theta }}^{\prime }\mathbf{F}(\mathbf{x},\mathbf{J})-\psi (\mathbf{\theta },\mathbf{J})\right].$$

$$\mathbf{F}(\mathbf{x},\mathbf{J})=\left[\begin{array}{c}\sum _{i=1}^N{h}_i^1{x}_i+\sum _{j>i}{j}_{ij}^1{x}_i{x}_j\\ \sum _{i=1}^N{h}_i^2{x}_i+\sum _{j>i}{j}_{ij}^2{x}_i{x}_j\\ ⋮\\ \sum _{i=1}^N{h}_i^D{x}_i+\sum _{j>i}{j}_{ij}^D{x}_i{x}_j\end{array}\right],$$

$$\mathbf{J}$$

$$\mathbf{F}(\mathbf{x},\mathbf{J})={\mathbf{J}}^{\prime }\tilde{\mathbf{F}}(\mathbf{x}).$$

$$\mathbf{\eta }$$

$$={\left[\begin{array}{cccc}{\mathbf{j}}_1^{\prime }\tilde{\mathbf{\eta }}& {\mathbf{j}}_2^{\prime }\tilde{\mathbf{\eta }}& \dots & {\mathbf{j}}_D^{\prime }\tilde{\mathbf{\eta }}\end{array}\right]}^{\prime },$$

$$\mathbf{G}$$

$$=E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}\left[{\mathbf{J}}^{\prime }\tilde{\mathbf{F}}(\mathbf{x})\tilde{\mathbf{F}}(\mathbf{x}{)}^{\prime }\mathbf{J}\right]-E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}{\mathbf{J}}^{\prime }\tilde{\mathbf{F}}(\mathbf{x})E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}\tilde{\mathbf{F}}(\mathbf{x}{)}^{\prime }\mathbf{J}$$

$$={\mathbf{J}}^{\prime }\left\{E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}[\tilde{\mathbf{F}}(\mathbf{x})\tilde{\mathbf{F}}(\mathbf{x}{)}^{\prime }]-E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}\tilde{\mathbf{F}}(\mathbf{x})E_{\mathbf{X}|\mathbf{\theta },\mathbf{J}}\tilde{\mathbf{F}}(\mathbf{x}{)}^{\prime }\right\}\mathbf{J}$$

$$={\mathbf{J}}^{\prime }\tilde{\mathbf{G}}\mathbf{J}.$$

$${\mathbf{\theta }}_t={\mathbf{\theta }}_{t-1}+{\mathbf{\xi }}_t(\lambda ).$$

$$p(\mathbf{x}|{\mathbf{\theta }}_t,{\mathbf{J}}_t)$$

$$=\exp \left[\sum _{k=1}^D{\theta }_t^k\left[\sum _{i=1}^N{h}_i^{k,t}{x}_i+\sum _{i<j}{j}_{ij}^{k,t}{x}_i{x}_j\right]-\psi ({\mathbf{\theta }}_t,{\mathbf{J}}_t)\right]$$

$$=\exp \left[{\mathbf{\theta }}_t^{\prime }\mathbf{F}(\mathbf{x},{\mathbf{J}}_t)-\psi ({\mathbf{\theta }}_t,{\mathbf{J}}_t)\right].$$

$$p({\mathbf{\theta }}_t|{\mathbf{X}}_{1:t},{\mathbf{J}}_{1:t})=\frac{p({\mathbf{X}}_t|{\mathbf{\theta }}_t,{\mathbf{J}}_t)p({\mathbf{\theta }}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t-1})}{p({\mathbf{X}}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t})}.$$

$$p({\mathbf{\theta }}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t})=p({\mathbf{\theta }}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t-1}).$$

$$p({\mathbf{\theta }}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t-1})$$

$$=\int p({\mathbf{\theta }}_t|{\mathbf{\theta }}_{t-1})p({\mathbf{\theta }}_{t-1}|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t-1})d{\mathbf{\theta }}_{t-1}.$$

$${\mathbf{\theta }}_{t|t-1}$$

$${\mathbf{W}}_{t|t-1}$$

$$f({\mathbf{\theta }}_t)$$

$$={\mathbf{\theta }}_t^{\prime }\mathbf{F}({\mathbf{X}}_t,{\mathbf{J}}_t)-\psi ({\mathbf{\theta }}_t)$$

$$-\frac{1}{2}({\mathbf{\theta }}_t-{\mathbf{\theta }}_{t|t-1}{)}^{\prime }{\mathbf{W}}_{t|t-1}^{-1}({\mathbf{\theta }}_t-{\mathbf{\theta }}_{t|t-1})+const.$$

$${\mathbf{\theta }}^{\mathrm{new}}={\mathbf{\theta }}^{\mathrm{old}}-[\nabla \nabla f({\mathbf{\theta }}^{\mathrm{old}}){]}^{-1}\nabla f({\mathbf{\theta }}^{\mathrm{old}}).$$

$$\nabla f(\mathbf{\theta })=\mathbf{F}({\mathbf{X}}_t,{\mathbf{J}}_t)-\mathbf{\eta }(\mathbf{\theta })-{\mathbf{W}}_{t|t-1}^{-1}(\mathbf{\theta }-{\mathbf{\theta }}_{t|t-1}),$$

$$\nabla \nabla f(\mathbf{\theta })=-\mathbf{G}(\mathbf{\theta })-{\mathbf{W}}_{t|t-1}^{-1},$$

$${\mathbf{W}}_{t|t}=-[{\nabla \nabla f(\mathbf{\theta })|}_{{\mathbf{\theta }}_{t|t}}{]}^{-1}=[\mathbf{G}({\mathbf{\theta }}_{t|t})+{\mathbf{W}}_{t|t-1}^{-1}{]}^{-1}.$$

$$l({\mathbf{J}}_t)\equiv \log p({\mathbf{X}}_t|{\mathbf{X}}_{1:t-1},{\mathbf{J}}_{1:t-1},{\mathbf{J}}_t).$$

Full text

I Introduction

Networks of interacting elements, or graphs, describe varieties of systems, such as communication systems, biological and ecological systems, and social media systems [1]. Examples include people’s communications via emails, infection of diseases, or communications between neurons in nervous systems. A great majority of these networks can be simplified to undirected graphs of binary units, where nodes take either the active or inactive state, and connections between two nodes do not have directionality. For example, spiking activity of neural populations has been investigated using models of binary patterns [2]. In addition, sequences of correlated binary patterns can represent many artificial time-series data, such as text, movies, and music [3].

Networks of pairwise connected elements are generally described as a Markov random field, or an undirected graphical model; with nodes and edges representing random variables and their relations respectively. This model can be written using the exponential family distribution. When the random variables are binary, their probability distribution follows the celebrated Ising model, or the Boltzmann machine– a stochastic extension of the Hopfield memory network. However, it is unlikely that observed binary patterns are sampled from such a stationary distribution except for controlled conditions. Hence, the model was extended so it can capture time-varying structures in the data [4, 5, 6, 7, 8, 9, 10, 11, 12].

Of these time-dependent Ising models, methods based on the state-space framework offer the sequential Bayes algorithm that provides estimates of the time-varying network with Bayesian credible intervals, and the expectation-maximization (EM) algorithm that optimizes parameters in the model [4, 9, 10, 12]. However, while estimating a time-dependent full Ising network is an important approach in elucidating binary activity dictated by all the pairwise interactions, a few dominant patterns may capture these time-series due to the innate correlations. Capturing essential structures in the time-varying data is an important tool to perform compression, prediction, and control of these patterns. In this study, we provide an online method to decompose the underlying structure, assuming different time-scales for hierarchical structures in the binary data generation.

Here we construct a low-dimensional representation of the Ising model that summarizes networks as a weighted combination of undirected graphs, assuming that the weights follow independent processes and that the underlying structure of the graphs changes slowly compared to the speed of weight changes. We provide a sequential Bayes algorithm to trace the dynamics of the weights, which we call network states, while updating the underlying graph structures online. The process of learning the underlying graphs is much slower than the network dynamics caused by weight changes, which allows separation of the two dynamics. The method is thus expected to extract temporal hierarchies in natural and artificial systems that generate binary patterns. This approach makes a contrast to previous offline unsupervised methods based on orthogonal or low-rank representation of the network dynamics to reduce dimensionality of time-dependent Ising models [13, 14].

This paper is constructed as follows. In Methods, we explain (A) the stationary Ising model of binary patterns and construct (B) a low-dimensional model composed of multiple weighted graphs. We then introduce (C) the state-space model in which the weights of the graphs change. We derive (D) the sequential Bayes algorithm for estimating the weights as well as (E) an online update algorithm of the graphs. In Results, (A) we corroborate the algorithm using simulated data, and demonstrate that the method extracts time-varying features better than a traditional orthogonal decomposition method. Next, (B) the model selection is performed to select the number of features in the data. We confirm the selected model exhibits improved goodness-of-fit over the original time-dependent Ising model. Finally, (C) we show the applicability of the method in analyzing neural data. In Discussion, we conclude the paper with possible extensions of the method.

II Methods

II-A The Ising model

The time-series model of binary data we propose in this study is based on the Ising/spin-glass model in statistical physics. In order to specify notations, here we describe the Ising model. This model will be extended to the multi-graph and time-dependent model in the next sections. For $N$ binary variables, the Ising model is given as

[TABLE]

Here the pattern $\mathbf{x}=(x_{1},\ldots,x_{N})^{\prime}$ is a $N$-tuple binary vector, where each entry indicates the active $x_{i}=1$ or inactive $x_{i}=0$ state of the $i^{th}$ node. Each node of the pattern can represent a pixel in image recognition or the activity of a neuron in neuroscience spiking data. $\psi(\mathbf{j})$ is a log normalization term. It ensures that the probabilities of all patterns sum up to 1. $h_{i}$ is the bias parameter for the $i^{th}$ node whereas $j_{ij}$ represents the interaction between the $i^{th}$ and the $j^{th}$ nodes. We express the parameters by a single column vector:

[TABLE]

We call this vector the graph of the Ising model. By adjusting the graph, we can embed correlated patterns such as binary images as the ones that frequently appear in the model. The Ising model is written in the canonical form as

[TABLE]

where $\tilde{\mathbf{F}}(\mathbf{x})$ is the feature vector computed from the binary variables as

[TABLE]

The log normalization is then computed as

[TABLE]

The computation of $\psi(\mathbf{j})$ requires considering all possible binary patterns. The first derivative of $\psi(\mathbf{j})$ with respect to $\mathbf{j}$ provides the expectation of the feature vector:

[TABLE]

where $\mathbf{X}$ is a sample of binary data, and $E_{\mathbf{X}|\mathbf{j}}$ represents the expectation of $\mathbf{X}$, using $p(\mathbf{x}|\mathbf{j})$. Furthermore, the second derivative yields the Fisher information matrix:

[TABLE]

We will utilize these relations frequently in the next sections.

II-B The multi-graph Ising model

Next, we construct an Ising model composed of multiple graphs, assuming that the data is sampled from a network that is a weighted combination of these graphs. Namely, we consider the following model:

[TABLE]

where $D$ is the number of undirected graphs. We call this model a multi-graph Ising model. $h_{i}^{k}$ and $j_{ij}^{k}$ correspond to the graph parameters of the $k^{th}$ graph. $\theta^{k}$ is the weight applied to the $k^{th}$ graph, and we define $\boldsymbol{\theta}=(\theta^{1},\theta^{2},\ldots,\theta^{D})^{\prime}$. If $h_{i}^{k}$ and $j_{ij}^{k}$ are known, the weights $\theta^{k}$ constitute the canonical parameters of the model. Given the graph parameters and $D<N+N(N-1)/2$, this model can represent the data with a smaller dimensionality. The multi-graph Ising model can be written in the canonical form as

[TABLE]

Here, the feature vector, $\mathbf{F}(\mathbf{x},\mathbf{J})$, is given by:

[TABLE]

and $\mathbf{J}$ is the matrix containing all graphs:

[TABLE]

Using $\mathbf{J}$, the feature vector can be further expressed as

[TABLE]

Note that $\tilde{\mathbf{F}}(\mathbf{x})$ is the feature vector of the original Ising model of Eq. 1. The log normalization term, $\psi(\boldsymbol{\theta},\mathbf{J})$, is then calculated in the same way as in the original Ising model: $\psi(\boldsymbol{\theta},\mathbf{J})=\log\sum_{\mathbf{x}}\exp[\boldsymbol{\theta}^{\prime}\mathbf{F}(\mathbf{x},\mathbf{J})]$.

Like in the original Ising model, the first and second-order derivatives of the log-normalization term respectively provide the expectation and the second-order cumulant of the new feature vector $\mathbf{F}(\mathbf{x},\mathbf{J})$. These can be expressed using the expectation parameters and the Fisher information of the original Ising model (Eqs. 6,7). The expectation of the feature vector is given as

[TABLE]

where $E_{\mathbf{X}|\boldsymbol{\theta},\mathbf{J}}$ represents the expectation of $\mathbf{X}$, using $p(\mathbf{x}|\boldsymbol{\theta},\mathbf{J})$. Similarly, the second-order cumulant of the feature vector, i.e., the Fisher information matrix $\mathbf{G}$, is given as

[TABLE]

II-C The state-space multi-graph Ising model

The model developed so far can only produce stationary data. In order to analyze the dynamics of binary patterns, it is essential to consider variations of model parameters in time. Figure 1 shows the state-space modeling framework used for this purpose. We assume the following state equation to link the states of two consecutive time bins:

[TABLE]

At each time step, we add a zero mean Gaussian noise, $\boldsymbol{\xi}_{t}$, that has a covariance given by $\lambda^{-1}\mathbf{I}$, where $\mathbf{I}$ is a $D\times D$ identity matrix. $\lambda$ is time-independent and will be fixed. This state equation assumes that the dynamics of the network states are independent Gaussian processes.

To use this state equation, it is necessary to assume the state of the model at the first time bin. We assume its probabilistic distribution to be Gaussian, with mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$: $\boldsymbol{\theta}_{1}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$. We will fix these hyper-parameters ($\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$). However, unlike the algorithm developed in [4, 9, 10, 12], here we assume that one of the hyper-parameters, namely $\mathbf{J}$, is also time-dependent, and denote it as $\mathbf{J}_{t}$ for $t=1,\ldots,T$.

At each time step, we assume that the data is sampled from the current multi-graph Ising model with network state $\boldsymbol{\theta}_{t}$ and graphs $\mathbf{J}_{t}$:

[TABLE]

This observation model is the main subject of this study. Eqs. 15 and II-C constitute the state-space model of the dynamic pattern sequences. With this augmentation, related parameters of the model ($\boldsymbol{\theta}_{t}$, $\mathbf{J}_{t}$, $\boldsymbol{\eta}_{t}$ and $\mathbf{G}_{t}$) now become time-dependent. The log normalization $\psi(\boldsymbol{\theta}_{t},\mathbf{J}_{t})$ is also time-dependent as it is a function of $\boldsymbol{\theta}_{t}$ and $\mathbf{J}_{t}$. We assume that the graph parameter $\mathbf{J}_{t}$ follows much slower dynamics than the changes of weights $\boldsymbol{\theta}_{t}$, which makes their respective effect on the binary patterns distinguishable. In the next sections, we aim to estimate $\boldsymbol{\theta}_{t}$ and $\mathbf{J}_{t}$ from time-varying data. The algorithm that will be developed in these sections is summarized in Table I.

II-D Sequential Bayes estimation

In this subsection, we provide the sequential Bayes algorithm to estimate the time-dependent weights $\boldsymbol{\theta}_{t}$ online. Here, we denote the graphs up to time $t$ as $\mathbf{J}_{1:t}=\left\{\mathbf{J}_{1},\mathbf{J}_{2},\ldots,\mathbf{J}_{t}\right\}$. We use the same notation for the data up to time $t$ ($\mathbf{X}_{1:t}$).

The online estimation of the parameters is done by the recurrent construction of the filter density. The filter density using Bayes’ theorem is given as

[TABLE]

To obtain the second term of the numerator, we used:

[TABLE]

This holds because $\mathbf{J}_{t}$ is independent of $\boldsymbol{\theta}_{t}$ given that $\mathbf{X}_{t}$, which connects $\boldsymbol{\theta}_{t}$ and $\mathbf{J}_{t}$, is marginalized. The estimate of the state at time $t$ ($\boldsymbol{\theta}_{t}$) given the data up to time $t-1$ is called one-step prediction. The one-step prediction density is derived from the Chapman-Kolmogorov equation [15]:

[TABLE]

Because the transition of the state is a Gaussian process, given that the filter density at $t-1$ is a Gaussian density with mean $\boldsymbol{\theta}_{t-1|t-1}$ and covariance $\mathbf{W}_{t-1|t-1}$, the one-step prediction density is also Gaussian with mean and covariance given by

[TABLE]

The filter density is a combination of the multi-graph Ising model and the one-step prediction density. Here, we approximate the posterior by a Gaussian using the Laplace method. As such, the mean of the approximated Gaussian density is identified as the mode of the posterior. We denote this mode as $\boldsymbol{\theta}_{t|t}$. It maximizes the log posterior given by

[TABLE]

Because of the log-concavity of both the observation and the state models, it is guaranteed that $f(\boldsymbol{\theta}_{t})$ contains a unique mode that can be found by convex optimization techniques. Here we use the Newton-Raphson method:

[TABLE]

In this equation, the first and second derivatives are given by

[TABLE]

where $\boldsymbol{\eta}$ is the expectation of the feature by $p(\mathbf{x}|\boldsymbol{\theta},\mathbf{J}_{t})$ computed as in Eq. 13. $\mathbf{G}$ is the Fisher information matrix calculated by Eq. 14, where the expectation is performed by $p(\mathbf{x}|\boldsymbol{\theta},\mathbf{J}_{t})$. Given the mode found by the above algorithm, the filtered covariance at each time step is approximated using the Hessian of the log posterior evaluated at $\boldsymbol{\theta}_{t|t}$:

[TABLE]

This approximate Gaussian filter is then used to make a prediction density at $t+1$, which completes the recursion.

II-E Identification of the underlying graphs

We combine the sequential Bayes estimation of the time-dependent weight parameters with the online estimation of the graphs $\mathbf{J}_{t}$. Our approach is based on the maximum likelihood estimation of $\mathbf{J}_{t}$ using the stochastic gradient method. We optimize $\mathbf{J}_{t}$ under the principle of maximizing the marginal likelihood at time step $t$ given the past observations up to $t-1$:

[TABLE]

This optimization will be performed by the stochastic gradient method:

[TABLE]

starting with a nominal initial estimate for $\mathbf{J}_{1}$. With this approach, we slightly update $\mathbf{J}_{t}$ in the direction toward the optimal parameters at time $t$ to obtain the graphs at time $t+1$, $\mathbf{J}_{t+1}$. This approach that adds the contribution of the data toward the optimal graph parameters at every time step allows us to estimate them in an online fashion.

To employ the stochastic gradient method, we use the fact that the derivative of the marginal log-likelihood function at time step $t$ can be evaluated by the derivative of an alternative lower bound as shown below. For this goal, we introduce the lower bound of the marginal log-likelihood obtained by the following Jensen’s inequality, similarly to the EM algorithm [15]. Here, in order to search an optimal graph at time $t$, we introduce a variable $\mathbf{J}^{\ast}_{t}$. Then,

[TABLE]

where $\left\langle\cdot\right\rangle_{\boldsymbol{\theta}_{t}|\mathbf{X}_{1:t},\mathbf{J}_{1:t}}$ is the expectation by the filter density at time $t$. The second term of the last equality is the entropy of the filter density, which is not a function of $\mathbf{J}^{\ast}_{t}$. The first term is the expected complete data log-likelihood, a.k.a. Q-function:

[TABLE]

In the second term of the last equality, we dropped $\mathbf{J}^{\ast}_{t}$, following the same argument made for Eq. 18. Hence, in the lower bound of the marginal log-likelihood, only the first term of the Q-function is a function of $\mathbf{J}^{\ast}_{t}$ :

[TABLE]

where $c$ is a function independent of $\mathbf{J}^{\ast}_{t}$.

We note that the derivative of the marginal log-likelihood in Eq. 28 can be replaced with the derivative of the Q-function since we expect

[TABLE]

if the approximate posterior is close to the exact posterior. The derivation of Eq. 32 for the exact posterior is given at the end of this section. This equation means that, at $\mathbf{J}_{t}$, the optimal direction to update the graphs for maximizing the Q-function is the same as the direction that optimally maximizes the marginal log-likelihood. According to this, the online update of the underlying graphs at each time step with a gradient ascent scheme is given as:

[TABLE]

Here we obtain the derivative of the Q-function, using

[TABLE]

[TABLE]

where $E_{\mathbf{X}|\boldsymbol{\theta}_{t},\mathbf{J}^{\ast}_{t}}$ represents the expectation by $p(\mathbf{x}|\boldsymbol{\theta}_{t},\mathbf{J}^{\ast}_{t})$. Using the above results, the derivative of the Q-function with respect to $\mathbf{J}^{\ast}_{t}$ is expressed as

[TABLE]

Updating $\mathbf{J}_{t}$ requires evaluating Eq. 36 at $\mathbf{J}^{\ast}_{t}=\mathbf{J}_{t}$, which requires the evaluation of the expectation with respect to the filter density at time $t$. This can be computed by sampling from the approximate Gaussian posterior distribution, $\boldsymbol{\theta}_{t}\sim\mathcal{N}(\boldsymbol{\theta}_{t|t},\mathbf{W}_{t|t})$. Finally, we sketch the equality for Eq.32 for an exact posterior density:

[TABLE]

This completes the stochastic gradient search of the underlying graphs under the maximum likelihood principle.

III Results

III-A Application to artificially generated data

Generation of artificial data To demonstrate the applicability of the method, we fitted the model to synthesized binary data representing dynamically changing images (Fig. 2), in which each node ($x_{i}$) corresponds to a pixel that can either be black ($x_{i}=1$) or white ($x_{i}=0$). At each step, a sample was generated by Eq. II-C, a mixture of multiple graphs $\mathbf{J}$ with time-dependent weights $\boldsymbol{\theta}_{t}$. Here we explain how we constructed $\mathbf{J}$ and $\boldsymbol{\theta}_{t}$ to generate data.

We generated the data from a mixture of two fixed graphs, each containing 9 nodes. We used a constant $\mathbf{J}$ to synthesize the data. If the graphs are used separately, the model most frequently generates an image similar to a ‘+’ or a ‘T’ respectively (Fig. 2A). Note that these images have overlapped components. For the time-dependent weights, we chose sinusoidal waves of period $T_{\rm ep}$=1500 time steps. The middle panel of Fig. 2B shows the dynamics of the weights of the underlying model (dashed lines) over one epoch of the data. We set their baseline and amplitude to 0.5 with a phase shift of $\pi$ so that each graph has in turn a dominant contribution while the other has a small contribution. We repeatedly generated the time-series according to the above procedure, and obtained 200 epochs of the binary data ($T$=$300\,000$ time steps).

On the top panel of Fig. 2B, we show images sampled from the underlying model at different time steps in a single epoch. In particular, because of the chosen dynamics of the weights, we expect to obtain with a high probability an image corresponding to a ’+’ and a ’T’ at $t=0.25T_{\rm ep}$ and $t=0.75T_{\rm ep}$ respectively. The other images are sampled from a mixture of the two graphs and do not necessarily resemble either one of the main images.

Estimating the state-space multi-graph Ising model On the middle panel of Fig. 2B, we also show the dynamics of the weights of the fitted model (solid lines) at the last epoch of the data. These are obtained by the sequential Bayes algorithm. We show their credible interval computed from the covariance matrix of the posterior density of the weights as $\boldsymbol{\theta}_{t|t}^{k}\pm 2\sqrt{\left[\mathbf{W}_{t|t}\right]_{k,k}}$, where $\left[\mathbf{W}_{t|t}\right]_{k,k}$ corresponds to the $k^{th}$ element on the diagonal of the covariance matrix. The dynamics of the fitted model seem to adequately concord with those of the underlying model, as the underlying weights are generally found in the credible intervals of the estimates. On the bottom panel of Fig. 2B, we show images sampled from the fitted model for the last epoch of the data. The sampled images from the fitted model became closer to those generated artificially as the online estimation of $\mathbf{J}$ progressed.

At a given time $t$, the marginal log-likelihood of the data can be approximated by (see [9]):

[TABLE]

The top panel of Fig. 2C shows the average of the marginal log-likelihoods obtained for each epoch. Namely, for epoch $r$, we computed:

[TABLE]

for $r=1,...,200$. The marginal log-likelihood of our fitted model increased with more epochs of the data.

For this simple example, we expected that the matrix $\mathbf{J}$ that generated the data would be reconstructed by the online algorithm. In order to verify the goodness of the estimated $\mathbf{J}$, we computed the correlation coefficient between the columns of the fitted and underlying $\mathbf{J}$ matrices:

[TABLE]

where $\mathbf{j}^{k}$ corresponds to the $k^{th}$ column of the $\mathbf{J}$ matrix and the dot denotes the inner product of two vectors. Both graphs had a correlation coefficient close to 0 at the beginning (the first guess of $\mathbf{J}$ was randomized as a zero-mean, unit-variance Gaussian) and it progressed toward 1 (Fig. 2C Bottom panel).

Finally, in order to show that the proposed method captured the features better than a traditional orthogonal decomposition method, we performed a principle component analysis (PCA) on the MAP estimate of the time-dependent full Ising model that includes only a single graph as in Eq. 1, but whose parameters all vary in time [4, 9]. We used the same binary data to fit the model. The data matrix considered to perform the PCA was constituted of the MAP estimates at every time step. Features extracted by PCA were less correlated with the generative graphs: the absolute values of the correlation coefficients between the columns of the generative $\mathbf{J}$ and PC1 and PC2 were below 0.6, which is lower than the result obtained by our method ($\sim$ 0.7).

III-B Model selection

In practical application of the method, the exact amount of graphs from which the data is sampled ($D$) is unknown. Hence, it is necessary to specify how many graphs we need to include in the model. To obtain the most predictive model that avoids overfitting to the data, we compare the models by their respective Akaike Information Criterion (AIC) given by

[TABLE]

as a proxy for selection criteria of the online method. Here, $l(\mathbf{J}_{T-T_{\rm ep}:T})=\sum_{t=T-T_{\rm ep}}^{T}l(\mathbf{J}_{t})$ and $m$ is the number of free parameters, which in our case is the number of elements in the $\mathbf{J}$ matrix. Here we compute the marginal log-likelihood of the models at the last epoch of the data. The model with the smallest AIC should be selected.

In Fig. 3, we show that the AIC properly identified the number of graphs in artificially generated data. We generated 200 epochs of data by constructing a model containing two 9-node, zero-mean, unit-variance random Gaussian graphs. Fig. 3A shows the dynamics of the weights used to generate the data over one epoch. We then fitted the model to this data using up to 4 graphs. The same values were used for all other hyper-parameters ($\boldsymbol{\mu}$, $\boldsymbol{\Sigma}$, $\lambda$). Fig. 3B shows the fitted weight dynamics. The AIC of each model is shown in the titles of the plots. The model containing 2 graphs has the lowest AIC, which means that using the AIC as a discriminating criterion can avoid over-fitted models.

Last, we compared the goodness-of-fit of our model to that of a time-dependent full Ising model [4, 9]. For the full model, the mean of the prior of $\mathbf{j}$ was the only free parameter, so $m$ is equal to the number of elements in $\mathbf{j}$. The AIC obtained ($27\,372$) is higher than the AIC of the multi-graph models. This confirms that reducing the dimensionality can lead to a better fit to data sampled from the multi-graph Ising model.

III-C Application to neural data

Next, we applied our model to neural data. We analyzed spontaneous activity of neurons recorded from rat hippocampal dissociated cultures using multi-electrode arrays [16, 17]. In these experiments, the authors recorded the cultures for about 5 weeks while neurons grew connections, and reported that the neural activity approached a critical state over time [16]. The recording length is approximately 1 hour per day. We analyzed one culture on day 34 that contained 85 neurons. We chose the 12 neurons that showed the highest firing rates. We used 10 ms bins to construct binary data (Fig. 4A), and fitted the multi-graph Ising models with the number of graphs ($D_{\rm{fit}}$) equal to 1, 2, and 3. Fig. 4B shows the estimated weights $\boldsymbol{\theta}_{t}$ for each model. Correlation coefficients between the learned $\mathbf{J}$ at the last step and the graphs obtained at each step confirms that the learning of $\mathbf{J}$ was completed by the first half of the observation period (Fig. 4C). The learning for models with more than 3 graphs was not completed in this period, hence these models were excluded from the analysis. Notably, in the latter half during which $\mathbf{J}$ is stable, the weight of the single-graph model is stationary whereas weights of multi-graph models vary in time. Since the learned graphs within the multi-graph models ($D_{\rm{fit}}=2,3$) are significantly correlated, the coordinated weight changes mostly capture the stationary activity of the population (not shown). Nevertheless differential dynamics of the multi-graphs may capture non-stationary activity of the population. To evaluate the predictive power of the models, we computed their AICs as its proxy, using the data in the latter half of the observation period. The model with 3 graphs was selected. This result suggests that the spontaneous activity of cultured neurons is not stationary, as reported in [18].

IV Discussion

We proposed an algorithm that can identify graphs underlying binary time-series data. Our model is composed of multiple time-dependent Ising graphs with weight parameters. The model is fitted by sequential Bayes estimation for the weights and stochastic gradient based on maximizing the likelihood for the graphs. The method was corroborated by artificially generated and neural time-series data.

It should be noted that our method extracts multiple graphs whose weights follow the assumption specified by the prior density (Eq. 15). In the current study, we assumed independent Markov processes for the weight dynamics. Given the observation model without the prior, the choice of $\mathbf{J}$ is undetermined because the following observation model with $D$ graphs is identical for any choice of invertible $D$ by $D$ matrix $\mathbf{Z}$:

[TABLE]

Because of the chosen prior distribution, the current method tends to extract independent weight processes with the same time-scale. Despite this assumption about the processes, the determination of the graphs still suffers ambiguity. For example, a single graph with a constant weight is captured by multiple graphs with significantly correlated weight dynamics if a linear combination of the graphs can reconstruct the original graph, as seen in the analysis of neural data. In future works, we need to impose different time-scales for weight processes and/or regularize $\mathbf{J}$ to disambiguate the graph estimation.

The current model can be extended in several ways toward predicting complex time-series data. First, the number of binary units that can be treated by the model needs to be increased. The exact computations of $\psi$, $\tilde{\boldsymbol{\eta}}$ and $\tilde{\mathbf{G}}$ using Eqs. 5, 6 and 7 are only practically realizable with a limited number of nodes (up to approximately 20). For this goal, we can use the pseudolikelihood and TAP/Bethe approximation methods for the state-space Ising model introduced in [12]. Second, we can use the particle filter to model a non-Gaussian posterior as long as the dimension of the state in this model is small. Finally, an important extension is to model non-linear transitions of the state and capture long-term relationships between time bins, similarly to the modern recurrent neural network models [3].

V Conclusion

In summary, the proposed model can estimate undirected graphs and their underlying dynamics in binary data in an online manner. This model is a dynamic probabilistic network that aims to predict sequences of binary patterns, assuming that they are sampled from a mixture of Ising graphs. While further development is necessary for practical use, the model potentially provides a probabilistic approach for pattern prediction that can be an alternative to classical recurrent networks.

Acknowledgment

The authors thank Prof. Shinsuke Koyama for helpful discussions.