Data-driven inference of hidden nodes in networks

TL;DR

This paper introduces a novel data-driven method for inferring hidden nodes and their interactions in complex networks, improving accuracy over existing methods and applicable to biological, financial, and image data.

Contribution

The paper presents a new approach for identifying hidden variables and their interactions in complex systems from observed data, including estimating the number of hidden nodes.

Findings

Outperforms existing methods in simulated kinetic Ising models

Accurately infers hidden nodes in neuronal and stock market networks

Identifies approximately 60 clusters in MNIST digit classification

Abstract

The explosion of activity in finding interactions in complex systems is driven by availability of copious observations of complex natural systems. However, such systems, e.g. the human brain, are rarely completely observable. Interaction network inference must then contend with hidden variables affecting the behavior of the observed parts of the system. We present a novel data-driven approach for model inference with hidden variables. From configurations of observed variables, we identify the observed-to-observed, hidden-to-observed, observed-to-hidden, and hidden-to-hidden interactions, the configurations of hidden variables, and the number of hidden variables. We demonstrate the performance of our method by simulating a kinetic Ising model, and show that our method outperforms existing methods. Turning to real data, we infer the hidden nodes in a neuronal network in the salamander retina and a stock market network. We show that predictive modeling with hidden variables is significantly more accurate than that without hidden variables. Finally, an important hidden variable problem is to find the number of clusters in a dataset. We apply our method to classify MNIST handwritten digits. We find that there are about 60 clusters which are roughly equally distributed amongst the digits.

Equations

$$P({\sigma }_i(t+1)|\mathbf{\sigma }(t))=\frac{\exp ({\sigma }_i(t+1)H_i(t))}{\exp (H_i(t))+\exp (-H_i(t))},$$

$$W_{ij}=\sum _k⟨\delta H_i\delta {\sigma }_k⟩[C^{-1}{]}_{kj},$$

$$H_i(t)\leftarrow \frac{{\sigma }_i(t+1)}{⟨{\sigma }_i(t+1){⟩}_{\text{model}}}{H}_i(t)={\sigma }_i(t+1)\frac{H_i(t)}{\tanh H_i(t)}.$$

$$\mathcal{L}=\prod _{t=1}^{L-1}\prod _{i=1}^{N}P({\sigma }_i(t+1)|\mathbf{\sigma }(t)).$$

$$D_v\equiv \sum _{i=1}^{N_v}{D}_i(W)$$

$$D\equiv D_{v}\bigg{(}1+\frac{N_{h}}{N_{v}}\bigg{)}.$$

$$P({\sigma }_i^v(t)=\pm 1|{\mathbf{\sigma }}^h(t))=\frac{\exp (\pm H_i(t))}{\exp (H_i(t))+\exp (-H_i(t))}$$

Full text

I Introduction

To go from observations to predictive understanding is to go from stamp-collecting to science. Absent principled quantitative laws, biological and social systems can be generally described as networks of interacting nodes, with time-series data providing a window on the dynamics of the underlying system. In the present era of big data, the network reconstruction problem has attracted considerable interest in research areas ranging from neuroscience Schneidman et al. (2006); Dombeck et al. (2007); Nguyen et al. (2016); Bernal-Casas et al. (2017) and genomics Lezon et al. (2006); Hickman and Hodgman (2009); Bar-Joseph et al. (2012) to finance  Pincus and Kalman (2004); Tse et al. (2010); Tabak et al. (2010); Bury (2013). A fundamental caveat is that such reconstructions always rely on partial observation of these complex networks. For example, it is hopeless to follow the simultaneous spiking activity of every neuron in the brain, the transcription of every gene in the genome, and every fluctuating factor in a financial system.

The problem of accounting for the unobserved constituents of any system is ill-posed without further information, simply because the number, the interactions, and the configurations of these hidden nodes must all be identified from the observed data and, a priori, one can make the former two as large and as complicated, respectively, as one pleases. To render the problem well-defined, one can first choose a theoretical model structure and then account for the unobserved nodes within this structure. Given the importance of this problem, much work has been devoted to it.

A simple approach is to maximize the likelihood of observed configurations after marginalizing unobserved configurations Dunn and Roudi (2013). Another effective approach is the Expectation Maximization (EM) algorithm for hidden variables that contains two alternating steps, inferring all interactions of observed and hidden variables from configurations of observed variables, and reconstructing the configurations of hidden variables consistent with these inferred interactions Dempster et al. (1977). As one might expect, this algorithm is computationally impractical for even moderately large systems if the fraction of unobserved variables is significant. Furthermore, hidden variable configuration reconstruction accuracy is greatly dependent on interaction inference accuracy, a factor that becomes significant for limited datasets. Therefore, recent network reconstruction methods have considered alternative approaches, such as mean field approximations Dunn and Roudi (2013); Tyrcha and Hertz (2014) and replica methods Bachschmid-Romano and Opper (2014); Battistin et al. (2015). However, the mean field approximations work only for weak and dense interactions Tyrcha and Hertz (2014), whereas the replica methods allows to infer strong and sparse interactions, but impose the stringent assumption of the independence between hidden variables Battistin et al. (2015). In addition to non-interacting hidden variables, random interaction strengths and the thermodynamic limit are two prerequisites for the exact inference of the replica methods Bachschmid-Romano and Opper (2014).

We recently formulated a new approach Hoang et al. (2018a, b) to network reconstruction for observed variables that is significantly more accurate inference-wise in the limit of sparse sampling and orders of magnitude faster computation-wise than previous methods. Based on this foundation, we propose a new approach for network reconstruction including hidden variables, by replacing the inference step with our approach. This does not, by itself, address the crucial question of the number of unobserved variables, so we complete our proposal by formulating a simple quantitative test of model complexity to determine this number.

This paper is organized as follows: We briefly review our inference method and outline its extension to hidden variables, paying especial attention to the determination of the number and interactions of hidden variables, amongst themselves and with observed variables. We then validate our method with simulated data from kinetic Ising models, showing the accurate determination of the number and interactions of hidden variables for a range of observed fractions of systems, going up to $40\%$ hidden variables. Turning to real data, we apply our method to reconstruct a neural network from partially observed neuronal activities, and a stock-market network using data of opening and closing stock prices of 25 American companies. We validate our network reconstructions by reproducing observed neuronal activities by pinning just a few neuron configurations, and by exhibiting a profitable stock trading strategy based on our inferred network. Finally, we demonstrate that our approach is suited to unsupervised data clustering, as well, since cluster membership is a type of hidden variable. We estimate the number of hidden features that can explain the MNIST hand-written digit dataset. Complete source code with documentation is available Hoang et al. (2018c).

II Method

We explain our approach in the context of a concrete example for ease of understanding. Consider a stochastic dynamical system in which a vector of $N$ binary ($\pm 1$) variables $\bm{\sigma}=(\sigma_{1},\cdots,\sigma_{N})$ evolves stochastically according to the conditional probability:

[TABLE]

for $i=1,\cdots,N$. The local field $H_{i}(t)=\sum_{j}W_{ij}\sigma_{j}(t)$ represents the summed influence of the present state $\sigma_{j}(t)$ on the future state $\sigma_{i}(t+1)$ through the weight $W_{ij}.$ This kinetic Ising model has a model expectation, $\langle\sigma_{i}(t+1)\rangle_{\text{model}}=\tanh H_{i}(t)$. Generating $\bm{\sigma}(t)$ given $W_{ij}$ is easy, but inferring $W_{ij}$ given $\bm{\sigma}(t)$ is not trivial. Although numerous methods exist for the inverse problem Roudi and Hertz (2011); Mézard and Sakellariou (2011); Zeng et al. (2013), we recently proposed a new approach Hoang et al. (2018a, b). We give here a simplified intuitive account. The first step is the linear regression of $H_{i}=\sum_{j}W_{ij}\sigma_{j}$ between $H_{i}$ and $\sigma_{j}$. Suppose we know $H_{i}(t)$ and $\sigma_{j}(t).$ The coefficient $W_{ij}$ can then be obtained as usual:

[TABLE]

where $C_{jk}\equiv\langle\delta\sigma_{j}\delta\sigma_{k}\rangle$ is the covariance matrix for $\bm{\sigma}(t),$ with $\langle f\rangle\equiv L^{-1}\sum_{t=1}^{L}f(t)$ and $\delta f\equiv f-\langle f\rangle.$ The second step is the update of the observable,

[TABLE]

The multiplicative update of $H_{i}(t)$ corrects the magnitude and sign of $H_{i}(t)$ based on the ratio of observed $\sigma_{i}(t+1)$ and model expectation $\langle\sigma_{i}(t+1)\rangle_{\text{model}},$ which is always larger than unity in absolute magnitude. A critical aspect of Eq. (3) is that the limit $|H_{i}|\downarrow 0$ gives $H_{i}(t)\leftarrow\sigma_{i}(t+1),$ independent of $H_{i}.$ Therefore, the update in Eq. (3) avoids being entirely multiplicative for determining $W_{ij}.$ These two steps, $H_{i}(t)\rightarrow W_{ij}$ and $W_{ij}\rightarrow H_{i}(t),$ provide a powerful iterative method. We continue this iteration until the discrepancy between data and model expectation D_{i}(W)\equiv\sum_{t}\big{[}\sigma_{i}(t+1)-\langle\sigma_{i}(t+1)\rangle_{\text{model}}\big{]}^{2} is minimized. We derived the linear regression in Eq. (2) using the concept of free energy in statistical mechanics Hoang et al. (2018a) so we call this method Free Energy Minimization (FEM). Notice that the parameter update in Eqs. (2-3) is completely independent of the computation of $D_{i}.$ This crucial feature allows the small sample size inference to avoid overfitting because the minimization of $D_{i}$ is used only as a stopping criterion.

Now we propose to apply the FEM method to infer interactions from/to hidden variables. The system has $N_{v}$ observable (visible) and $N_{h}$ hidden variables ($N=N_{v}+N_{h}$). As a variant of the EM algorithm, we first assign random configurations for hidden variables. We then infer interaction weights $W_{ij}$ for observed-to-observed, hidden-to-observed, observed-to-hidden, and hidden-to-hidden variables with the FEM method. Given $W_{ij}$, we can update the configurations of hidden variables with a probability $\mathcal{L}_{2}/(\mathcal{L}_{1}+\mathcal{L}_{2})$ where $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ represent the likelihoods $\mathcal{L}$ of the system before and after flipping,

[TABLE]

Note that the independent terms in the update of hidden states (at each $t$) of $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ cancel in the update ratio, $\mathcal{L}_{2}/(\mathcal{L}_{1}+\mathcal{L}_{2}).$ Therefore, we just need to calculate the dependent terms. The iterations between the parameter optimization (M step) and the variable update (E step) provide accurate inference of the interaction weights, $W_{ij},$ and the unknown configurations of hidden variables.

We must now consider the problem of determining the number of hidden variables. A simple measure would be the same discrepancy between observation $\sigma^{v}_{i}(t+1)$ and model expectation $\langle\sigma^{v}_{i}(t+1)\rangle_{\text{model}},$

[TABLE]

but this is clearly not taking the hidden variables into account. On the other hand, extending the sum in Eq. (5) to include hidden variables is useless because the E step update is minimizing these additional terms already. Since the error in inference of hidden variable states cannot be set by a scale smaller than the model discrepancy in the observed part, we define the scaled discrepancy of the entire system based on the observed part as

[TABLE]

The first term in Eq. 6 represents the goodness of fit for observed variables and the second term represents model complexity, so our criterion balances the two. Because $D_{v}\varpropto-\log\mathcal{L}_{v}$ where $\mathcal{L}_{v}$ represents the likelihood of observed variables, Eq. 6 can be rewritten as $D\varpropto-\log\mathcal{L}_{v}(1+N_{h}/N_{v})$. Our criterion is thus similar in spirit to the Akaike information criterion Akaike (1974) and Bayesian information criterion Schwarz (1978) with log-likelihood of observation ($-\log\mathcal{L}_{v}\sim D_{v}$) and model degrees of freedom ($N_{v}+N_{h}$).

Finally, our method can be summarized as the following set of steps:

For a range of numbers of hidden variables, in parallel and independently,

(i) Assign configurations of hidden variables at random;

(ii) Infer interaction weights $W_{ij}$ including observed-to-observed, hidden-to-observed, observed-to-hidden, and hidden-to-hidden from the configurations of observed and hidden variables using FEM;

(iii) Flip the states of hidden variables with probability $\mathcal{L}_{2}/(\mathcal{L}_{1}+\mathcal{L}_{2})$ (see Eq. (4)).

(iv) Repeat steps (ii) and (iii) until the discrepancy of observed variables is minimized. The final values of $W_{ij}$ and hidden states are the inferred coupling weights and configurations of hidden spins, respectively.

Pick the number of hidden variables that minimizes Eq. (6).

III Results

III.1 Kinetic Ising model

To demonstrate the performance of our method, we synthesized binary time series of $N=100$ spins by using the Sherington-Kirkpatrick model Sherrington and Kirkpatrick (1975). The update of spin $\bm{\sigma}$ follows Eq. (1) with preset coupling strengths $W_{ij}$ (Fig. 1A). Our goal is to reconstruct all of $W_{ij}$ from observations of a fraction of $\bm{\sigma}(t).$ Suppose that we only observe the time series of 60 spins with 40 spins hidden (Fig. 1B). When we reconstructed the interactions $W_{ij}$ between observed node $i$ and observed node $j$ using FEM, the reconstructed $W_{ij}$, based on the partial observations, showed a large error (Fig. 1C). We introduced 40 hidden variables, and applied the EM algorithm outlined in the previous section. The reconstructed $W_{ij}$ was close to the true $W_{ij}$ (Fig. 1D). How well are the hidden variable configurations recovered? For the case of 40 hidden variables, the true configurations of the hidden variables were recovered with an accuracy of 96.6%. The reconstruction accuracy increased with fewer spins hidden (Fig. 2). For instance, when 90 spins were observed with 10 spins hidden, the accuracy was 97.6%.

The number of hidden variables is usually unknown in real-world problems. When we reconstructed $W_{ij}$ with different numbers of hidden variables, the mean square errors of observed-to-observed interaction strengths, MSE$=N_{v}^{-2}\sum_{i,j\in\text{obs}}(W_{ij}-W_{ij}^{\textrm{true}})^{2}$, were minimal at the right number of hidden variables (Fig. 3, upper panel). The MSE is also inaccessible in real-world problems, but the minimum of $D$ (Eq. (6)) captured the correct value of $N_{h}$ (Fig. 3, lower panel, red lines).

To reconstruct $W_{ij}$ from observed and hidden variables, we used FEM. For the M step, mean field methods such as naïve, Thouless-Anderson-Palmer, and exact mean field methods (nMF, TAP, and eMF), and maximum likelihood estimation (MLE) can also be used. A brief review of these methods can be found in Ref. Hoang et al. (2018a). Given partial observations, mean field approaches were not successful in reconstructing $W_{ij}$(Fig. 4). For a small percentage of hidden variables (90 observable and 10 hidden), FEM and MLE showed a similar performance in the reconstruction of observed-to-observed and hidden-to-observed interactions. However, FEM outperformed MLE in reconstructing observed-to-hidden and hidden-to-hidden interactions (Fig. 4A-D). For a large percentage of hidden variables (60 observable and 40 hidden variables), FEM showed significantly better performance even for observed-to-observed and hidden-to-observed interactions (Fig. 4E-H). We quantified the reconstruction performance by measuring MSE between $W_{ij}$ and $W_{ij}^{\text{true}}$. FEM showed more accurate reconstruction of $W_{ij}$ with lower MSE in every case. More importantly, in addition to better performance, FEM took approximately 100 times less computation time than MLE due to its multiplicative update (see Hoang et al. (2018a, b, c) for details).

III.2 Neuronal network

The analysis of real data brings out issues far more clearly than simulated validations. Therefore, we applied our method to infer hidden nodes and their contributions in a real neuronal network. We used the time series data of the 80 most active neurons from published multi-channel recordings of neuronal firing in the salamander retina Tkačik et al. (2014). Considering the existence of unobserved hidden neurons, we modeled the evolution of neuronal activities by defining a local field, $H_{i}(t)=H_{i}^{\text{ext}}+\sum_{j}W_{ij}\sigma_{j}(t)$, that determines the future activity of $\sigma_{i}(t+1)$. The external local field $H_{i}^{\text{ext}}$ represents the bias of the $i$th neuron that sets its threshold. For various numbers of hidden neurons $N_{h}$, we computed $H_{i}^{\text{ext}}$ and $W_{ij}$. The activities of observed neurons were explained better and $D_{v}$ kept decreasing with a larger $N_{h}$ of hidden neurons (Fig. 5B). However, once we considered the overall discrepancy $D,$ an optimal number of hidden neurons was $N_{h}^{*}=4.$ Thus, the inclusion of four hidden neurons best explained the activities of observed neurons, taking model complexity into account. Given these four hidden neurons, the connection weights $W_{ij}$ were reconstructed as shown in Fig. 5C.

Since $W_{ij}^{\text{true}}$ is unknown for the neuronal network, we validated our reconstruction in two different ways. First, we selected some neurons as input neurons, and then based on the activities of these input neurons, we generated the activities of remaining neurons by using the reconstructed $W_{ij}$. Here we selected the input neurons based on having the strongest influence to other neurons by gauging $\sum_{i}|W_{ij}|$. Given varying numbers of input neurons, we could successfully reconstruct the actual activities of the remaining neurons (Fig. 5D). As the number of input neurons increased, the reconstruction accuracy increased (Fig. 5F). Moreover, once the four hidden neurons were considered the reconstruction accuracy was significantly improved. Second, given $\bm{\sigma}(t)$, we predicted $\bm{\sigma}(t+1)$, and then calculated the covariance $C_{ij}=\langle\delta\sigma_{i}(t+1)\delta\sigma_{j}(t)\rangle$. The reconstructed covariance was comparable with the actual covariance from the observation (Fig. 5E).

III.3 Stock network

Our method has a wide range of practical applications. As a demonstration, we reconstruct a stock market network with possible hidden nodes. We used stock price time series of 25 major companies in the S&P 500 index in five different sectors: technology (AAPL, GOOGL, MSFT, INTC, IBM), finance (BRK.B, JPM, WFC, BAC, C), health care (JNJ, PFE, UNH, MRK, AMGN), consumer discretionary (AMZN, WMT, HD, DIS, EBAY), and energy and industrial (XOM, CVX, GE, BA, MMM) Fusion Media Limited (2018). We examined the price difference between daily opening and closing stock prices. Their fluctuations from January 2005 to July 2018 are shown in (Fig. 6A). First, instead of considering the continuous price fluctuations, we defined a discretized measure of price changes. If the daily price increased at time $t$ for the $i$th company (opening price $<$ closing price), we defined $\sigma_{i}(t)=+1$. However, if the price decreased (opening price $>$ closing price), then $\sigma_{i}(t)=-1$. Finally, if the price was unchanged (opening price $=$ closing price), we defined $\sigma_{i}(t)=\sigma_{i}(t-1)$.

We applied our method to this discretized data, and inferred external factors $H_{i}^{\text{ext}}$ and interacting factors $\sum_{j}W_{ij}\sigma_{j}(t)$ that stochastically determine $\sigma_{i}(t+1)$. Since FEM works well even for small sample sizes Hoang et al. (2018a), we divided the data into two-year periods to probe possible slower temporal changes in the interactions $W_{ij}$ between stock prices. In particular, we show results from more recent data for 2014 to 2016 (Fig. 6B-D) and 2016 to 2018 (Fig. 6E-F). The discrepancy $D_{v}$ between observed $\sigma^{v}_{i}(t+1)$ and model expectation $\langle\sigma^{v}_{i}(t+1)\rangle_{\text{model}}$ kept decreasing as expected when more hidden nodes were introduced (Fig. 6B and E). However, the entire discrepancy $D$, considering the model complexity with hidden variables, showed a minimum at $N^{*}_{h}\approx 4$-$5$ hidden nodes. The inferred stock market network including interactions between observed and hidden nodes is visualized in Fig. 6C and F. When we generated time series of stock prices using the reconstructed network, we found that the covariance of the generated sequences was consistent with the covariance of original sequences (Fig. 6D and G).

An accurate predictive network reconstruction should enable profitable trades. In particular, does our discrete reduction of the price data still contain enough information to be useful? First, we reconstructed the interactions between companies including the appropriate number of hidden nodes by using stock price data for the most recent $T$ days: $\bm{\sigma}^{v}(t-T+1),\bm{\sigma}^{v}(t-T+2),\cdots,\bm{\sigma}^{v}(t)$. Then, we predicted the price change direction $\bm{\sigma}^{v}(t+1)$ for the next day. Our strategy was to buy the stock $i$ that had the highest probability of increasing with a maximum $H_{i}(t)$, and to sell the stock $j$ that had the highest probability of decreasing with a minimum $H_{j}(t)$ at the beginning of the day. This trading strategy is expected to have a maximum profit bounded by (close price($i$) $-$ open price($i$)) $+$ (open price($j$) $-$ close price($j$)). The trading simulation from 2008 to 2018 with a moving time window $T=500$ days obtained 350% cumulative profit (Fig. 6H). This profit was significantly higher than the profit of 50% using random trades, which is due to the secular rise of the entire stock index. Furthermore, the reconstructed network including hidden nodes showed a larger profit than the 250% profit from the reconstructed network ignoring the hidden nodes. Next, we refined the trading strategy by buying/selling the stock that has the highest probability of increasing/decreasing but only if its price has decreased/increased on the previous day. In particular, this may result in only buying or only selling on any specific day. This new strategy produced the same cumulative profit in total, but it doubled the profit per transaction (Fig. 6I and J). Finally, we confirmed that the optimal time window for the highest profit was about $T=500$ days (Fig. 6J).

III.4 Classification of handwritten digits

Another potential application of our method, interpreting hidden states as labels, is for unsupervised classification. We demonstrate this idea with the MNIST data of handwritten digits Lecun et al. (1998). The data has 60,000 digit samples of 28$\times$28 pixel gray-scale (between 0 and 255) images obtained from 500 different individuals. Some of sample digits are shown in Fig. 7A. Our goal is to classify the 60,000 images into distinct clusters where each cluster represents different digits as well as different writing styles without using true labels. We formulated the classification problem as follows. Different digits and writing style combinations are encoded in hidden variable states $\bm{\sigma}^{h}$. Then, one realization of $\bm{\sigma}^{h}(t)$ generates a digit image $\bm{\sigma}^{v}(t)$, where $t$ is now being used to index the MNIST images. The feature has $N_{h}$ degrees of freedom with $\sigma^{h}_{J}(t),J=1,\cdots,N_{h}$. In particular, for simplicity, we adopted one-hot encoding by assigning only one nonzero element $\sigma^{h}_{J}(t)=1$ among $N_{h}$ elements of $\bm{\sigma}^{h}(t)$. Then, the generated image has binary values of $\sigma^{v}_{i}(t)=1$ (gray $>1$) or $\sigma^{v}_{i}(t)=-1$ (otherwise) for the $i$th pixel, which is determined by the conditional probability,

[TABLE]

where $H_{i}(t)\equiv\sum_{J}W_{iJ}\sigma_{J}^{h}(t)$ represents a local field acting on the $i$th pixel. Here, we ignored observed pixels $i$ if more than 95% samples had the same value. The threshold 95% showed similar results as a more restrictive threshold of 99%. Thus, for this setup, our reconstruction method considers only hidden-to-observed interactions. Briefly summarizing the inference procedure, we (i) assign a random binary vector $\bm{\sigma}^{h}(t)$, in which only one element has nonzero value ($\sigma_{J}^{h}(t)=1$); (ii) apply FEM to reconstruct the interaction strength $W_{iJ}$ from hidden label $J$ to observed pixel $i$; (iii) update the hidden states by assigning $\sigma_{J}^{h}(t)=1$ for the label $J$ that makes the likelihood of the observed pixels of sample $t$ the highest and $\sigma_{J}^{h}(t)=0$ for the other $N_{h}-1$ elements; (iv) repeat steps (ii) and (iii) until the discrepancy $D_{v}$ between $\sigma^{v}_{i}(t)$ and $\langle\sigma^{v}_{i}(t)\rangle_{\text{model}}$ saturates. Then, the one-hot hidden states $\bm{\sigma}^{h}$ represent distinct classes of MNIST images $\bm{\sigma}^{v}$.

We examined various possible numbers (10 to 100) of labels by controlling the number $N_{h}$ of hidden variables. As the hidden degrees of freedom $N_{h}$ increased, the model generated images of $\sigma^{v}_{i}(t)$ closer to the originals. In other words, the discrepancy $D_{v}$ kept decreasing as $N_{h}$ increased (Fig. 7B). However, once the model complexity was penalized with the overuse of the hidden degrees of freedom, an optimal degrees of freedom $N_{h}^{*}$ was determined with a minimum overall discrepancy $D.$ The estimate $N_{h}^{*}\approx 60$ means that the 60,000 MNIST images can be optimally clustered into about 60 classes of digits and writing styles. The mean images $1/N_{c}\sum_{t\in c}\sigma^{v}_{i}(t)$ corresponding to the 60 labels are shown in Fig. 7C. Here $N_{c}$ is the number of samples corresponding to label $c.$ It is of particular interest that each digit was divided into approximately six classes, suggesting that, in the MNIST dataset, about six different writing styles exist for every digit. To confirm the robustness of this result, we repeated the analysis with only 20,000 of the MNIST images, and obtained a similar conclusion (dashed lines in Fig. 7B).

IV Summary

Given partial observations of systems, complete network reconstruction is a longstanding problem in inference. In this paper, we propose a new iterative approach based on free energy minimization (FEM) and expectation maximization. We demonstrated on simulated systems that our method can accurately estimate the actual number of hidden variables from partial observations. Furthermore, network reconstruction was successful in recovering not only observed-to-observed interactions but also those involving hidden variables (hidden-to-observed, observed-to-hidden, and hidden-to-hidden). Hidden-to-hidden interactions are challenging to reconstruct with mean-field methods Dunn and Roudi (2013); Tyrcha and Hertz (2014). We applied this method to reconstruct a real neuronal network and a stock market network with the inclusion of possible hidden variables. The reconstructed networks were then validated by reproducing real neuronal activities and by a profitable trade simulation, respectively. Finally, as another potential application to unsupervised pattern classification, we found hidden labels in hand-written digit data.

FEM is more effective for network reconstruction than maximum likelihood estimation (MLE), because it separates the cost function evaluation from the independent multiplicative parameter update. This has two major benefits that are crucial for the application to hidden variable problems to succeed. The first is that the cost function can be used as a stopping criterion to avoid overfitting for small sample sizes, important when considering large numbers of possible hidden variables. The second is that the multiplicative update is computationally much more efficient (approximately 100 times faster than usual MLE-based network reconstruction methods), also critical for determining the configurations of hidden variables. Since the algorithm reconstructs interactions strengths $W_{ij}$ from the $j$th node to the $i$th node independently for the $i$th node, the network reconstruction can be easily parallelized, and therefore scaled to large system sizes.

Acknowledgment

This work was supported by Intramural Research Program of the National Institutes of Health, NIDDK (D.-T.H.,V.P.), and by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2016R1D1A1B03932264) (J.J.).