Optimized Realization of Bayesian Networks in Reduced Normal Form using Latent Variable Model

TL;DR

This paper presents algorithmic and structural solutions to reduce the computational and memory costs of Bayesian networks in reduced normal form, including an efficient online learning algorithm and a C++ library, with a focus on latent variable models.

Contribution

It introduces cost-reduction techniques and an online learning algorithm for Bayesian networks in FGrn, implemented in a C++ library, enhancing practical usability especially for latent variable models.

Findings

Significant reduction in computational costs.

Efficient online learning algorithm comparable to batch methods.

Improved performance in latent variable model structures.

Abstract

Bayesian networks in their Factor Graph Reduced Normal Form (FGrn) are a powerful paradigm for implementing inference graphs. Unfortunately, the computational and memory costs of these networks may be considerable, even for relatively small networks, and this is one of the main reasons why these structures have often been underused in practice. In this work, through a detailed algorithmic and structural analysis, various solutions for cost reduction are proposed. An online version of the classic batch learning algorithm is also analyzed, showing very similar results (in an unsupervised context); which is essential even if multilevel structures are to be built. The solutions proposed, together with the possible online learning algorithm, are included in a C++ library that is quite efficient, especially if compared to the direct use of the well-known sum-product and Maximum Likelihood (ML) algorithms. The results are discussed with particular reference to a Latent Variable Model (LVM) structure.

Equations

$$p_V(v)\propto f_V(v)b_V(v),$$

$${\mathbf{p}}_V\propto {\mathbf{f}}_V\odot {\mathbf{b}}_V,$$

$$\begin{array}[]{llll}b_{V^{(i)}}(v)&\propto&\displaystyle\prod_{\mathclap{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}}^{m}f_{V^{(j)}}(v)\quad\displaystyle\prod_{\mathclap{k=m+1}}^{M}b_{V^{(k)}}(v),&\;i=1:m\\ f_{V^{(i)}}(v)&\propto&\displaystyle\prod_{j=1}^{m}f_{V^{(j)}}(v)\quad\displaystyle\prod_{\mathclap{\begin{subarray}{c}k=m+1\\ k\neq i\end{subarray}}}^{M}b_{V^{(k)}}(v),&\;i=m+1:M\end{array},$$

$$\begin{array}[]{llll}{\bf b}_{V^{(i)}}&\propto&\displaystyle\bigodot_{\mathclap{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}}^{m}{\bf f}_{V^{(j)}}\quad\displaystyle\bigodot_{\mathclap{k=m+1}}^{M}{\bf b}_{V^{(k)}},&\quad i=1:m\\ {\bf f}_{V^{(i)}}&\propto&\displaystyle\bigodot_{j=1}^{m}{\bf f}_{V^{(j)}}\quad\bigodot_{\mathclap{\begin{subarray}{c}k=m+1\\ k\neq i\end{subarray}}}^{M}{\bf b}_{V^{(k)}},&\quad i=m+1:M\end{array}$$

$$P(Y|V)=[Pr\{Y=y_j|V=v_i\}{]}_{j=1:|\mathcal{Y}|}^{i=1:|\mathcal{V}|}=[{\theta }_{ij}{]}_{j=1:|\mathcal{Y}|}^{i=1:|\mathcal{V}|},$$

$$P(Y|V)=\left[\begin{array}{ccc}P(Y=y_1|V=v_1)& \cdots & P(Y=y_{|\mathcal{Y}|}|V=v_1)\\ P(Y=y_1|V=v_2)& \cdots & P(Y=y_{|\mathcal{Y}|}|V=v_2)\\ ⋮& \ddots & ⋮\\ P(Y=y_1|V=v_{|\mathcal{V}|})& \cdots & P(Y=y_{|\mathcal{Y}|}|V=v_{|\mathcal{V}|})\end{array}\right]=\left[\begin{array}{ccc}{\theta }_{11}& \dots & {\theta }_{1|\mathcal{Y}|}\\ {\theta }_{21}& \dots & {\theta }_{2|\mathcal{Y}|}\\ ⋮& \ddots & ⋮\\ {\theta }_{|\mathcal{V}|1}& \dots & {\theta }_{|\mathcal{V}||\mathcal{Y}|}\end{array}\right].$$

$$f_Y(y_i)\propto \sum _{j=1}^{|\mathcal{V}|}{\theta }_{ij}{f}_V(v_j)\text{and}{b}_V(v_j)\propto \sum _{i=1}^{|\mathcal{Y}|}{\theta }_{ij}{b}_Y(y_i),$$

$${\mathbf{f}}_Y\propto P(Y|V{)}^T{\mathbf{f}}_V\text{and}{\mathbf{b}}_V\propto P(Y|V){\mathbf{b}}_Y.$$

$$f_V(v_i)={\pi }_V(v_i),i=1:|\mathcal{V}|,$$

$${\mathbf{f}}_V={\mathbf{\pi }}_V.$$

$${\theta }_{lm}^{(1)}⟵\frac{{\theta }_{lm}^{(0)}}{{\sum }_{n=1}^N{f}_{V[n]}(l)}\sum _{n=1}^N\frac{f_{V[n]}(l)b_{Y[n]}(m)}{{\mathbf{f}}_{V[n]}^T{\theta }^{(0)}{\mathbf{b}}_{Y[n]}}.$$

$$P(Y_1,\dots ,Y_N,S_1,\dots ,S_H)=P(Y_1|S_1,\dots ,S_H)\cdots P(Y_N|S_1,\dots ,S_H)P(S_1)\cdots P(S_H).$$

$$P((S_1{S}_2\dots S_H{)}^{(1)}|S_1)$$

$$⋮$$

$$P((S_1{S}_2\dots S_H{)}^{(H)}|S_H)$$

$$P(Y_1,...,Y_N,L)=P(Y_1|L)\cdots P(Y_N|L)P(L).$$

$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\\ & ⋮& \end{array}\right]⟶\left[\begin{array}{c}0\\ 1\\ 2\\ 0\\ ⋮\end{array}\right].$$

$${\left[\begin{array}{c}{\upsilon }_1\\ {\upsilon }_2\\ ⋮\\ {\upsilon }_d\end{array}\right]|}_{norm}=\left[\begin{array}{c}{\zeta }_1\\ {\zeta }_2\\ ⋮\\ {\zeta }_d\end{array}\right]⟶{\zeta }_i=\frac{{\upsilon }_i}{{\sum }_{j=1}^d{\upsilon }_j}.$$

$$\left[\begin{array}{ccc}{\xi }_1& \dots & {\xi }_l\end{array}\right]\left[\begin{array}{ccc}{p}_{11}& \cdots & p_{1d}\\ ⋮& \ddots & ⋮\\ p_{l1}& \cdots & p_{ld}\end{array}\right]$$

$$\left[\begin{array}{ccc}{p}_{11}& \cdots & p_{1d}\\ ⋮& \ddots & ⋮\\ p_{l1}& \cdots & p_{ld}\end{array}\right]\left[\begin{array}{c}{\xi }_1\\ ⋮\\ {\xi }_d\end{array}\right]$$

$${\mathbf{b}}_{V^{(0)}}$$

$${\mathbf{f}}_{V^{(1)}}$$

$${\mathbf{f}}_{V^{(2)}}$$

$${\mathbf{f}}_{V^{(3)}}$$

$${\bar{\mathbf{b}}}_{Y_i}=\left[\begin{array}{c}\frac{1}{|{\mathcal{Y}}_i|}\\ ⋮\\ \frac{1}{|{\mathcal{Y}}_i|}\end{array}\right].$$

$${\mathbf{b}}_{S^{(i)}}=P(Y_i|S){\bar{\mathbf{b}}}_{Y_i}$$

$$=\left[\begin{array}{c}\sum _{j=1}^{|{\mathcal{Y}}_i|}{\theta }_{1j}{\bar{b}}_{Y_i}({\xi }_j)\\ ⋮\\ \sum _{j=1}^{|{\mathcal{Y}}_i|}{\theta }_{|\mathcal{S}|j}{\bar{b}}_{Y_i}({\xi }_j)\end{array}\right]=\left[\begin{array}{c}\frac{1}{|{\mathcal{Y}}_i|}\sum _{j=1}^{|{\mathcal{Y}}_i|}{\theta }_{1j}\\ ⋮\\ \frac{1}{|{\mathcal{Y}}_i|}\sum _{j=1}^{|{\mathcal{Y}}_i|}{\theta }_{|\mathcal{S}|j}\end{array}\right]=\left[\begin{array}{c}\frac{1}{|{\mathcal{Y}}_i|}\\ ⋮\\ \frac{1}{|{\mathcal{Y}}_i|}\end{array}\right].$$

$${\theta }_{lm}^{(1)}⟵{\theta }_{lm}^{(0)}\sum _{n=1}^N\frac{f_{V[n]}(l)b_{Y[n]}(m)}{{\mathbf{f}}_{V[n]}^T{\theta }^{(0)}{\mathbf{b}}_{Y[n]}},$$

$${\theta }^{(1)}⟵{\theta }^{(0)}\odot \sum _{n=1}^N\frac{{\mathbf{f}}_{V[n]}{\mathbf{b}}_{Y[n]}^T}{{\mathbf{f}}_{V[n]}^T{\theta }^{(0)}{\mathbf{b}}_{Y[n]}}.$$

$${\theta }^0\odot {\mathbf{f}}_{V[n]}{\mathbf{b}}_{Y[n]}^T$$

$$=\left[\begin{array}{ccc}{\theta }_{11}^{(0)}& \cdots & {\theta }_{1|\mathcal{Y}|}^{(0)}\\ ⋮& \ddots & ⋮\\ {\theta }_{|\mathcal{V}|1}^{(0)}& \cdots & {\theta }_{|\mathcal{V}||\mathcal{Y}|}^{(0)}\end{array}\right]\odot \left[\begin{array}{ccc}{\varphi }_1{\beta }_1& \cdots & {\varphi }_1{\beta }_{|\mathcal{Y}|}\\ ⋮& \ddots & ⋮\\ {\varphi }_{|\mathcal{V}|}{\beta }_1& \cdots & {\varphi }_{|\mathcal{V}|}{\beta }_{|\mathcal{Y}|}\end{array}\right]$$

$$=\left[\begin{array}{ccc}{\theta }_{11}^{(0)}{\varphi }_1{\beta }_1& \cdots & {\theta }_{1|\mathcal{Y}|}^{(0)}{\varphi }_1{\beta }_{|\mathcal{Y}|}\\ ⋮& \ddots & ⋮\\ {\theta }_{|\mathcal{V}|1}^{(0)}{\varphi }_{|\mathcal{V}|}{\beta }_1& \cdots & {\theta }_{|\mathcal{V}||\mathcal{Y}|}^{(0)}{\varphi }_{|\mathcal{V}|}{\beta }_{|\mathcal{Y}|}\end{array}\right],$$

$$\sum _{l=1}^{|\mathcal{V}|}\sum _{m=1}^{|\mathcal{Y}|}{\theta }^{(0)}\odot {\mathbf{f}}_{V[n]}{\mathbf{b}}_{Y[n]}^T=\sum _{l=1}^{|\mathcal{V}|}{\varphi }_l\sum _{m=1}^{|\mathcal{Y}|}{\theta }_{lm}^{(0)}{\beta }_m,$$

$$\begin{array}{rl}{\mathbf{f}}_{V[n]}^T{\theta }^{(0)}{\mathbf{b}}_{Y[n]}& =\left[\begin{array}{ccc}{\varphi }_1& \cdots & {\varphi }_{|\mathcal{V}|}\end{array}\right]\left[\begin{array}{ccc}{\theta }_{11}^{(0)}& \cdots & {\theta }_{1|\mathcal{Y}|}^{(0)}\\ ⋮& \ddots & ⋮\\ {\theta }_{|\mathcal{V}|1}^{(0)}& \cdots & {\theta }_{|\mathcal{V}||\mathcal{Y}|}^{(0)}\end{array}\right]\left[\begin{array}{c}{\beta }_1\\ ⋮\\ {\beta }_{|\mathcal{Y}|}\end{array}\right]\\ & =\left[\begin{array}{ccc}{\varphi }_1& \cdots & {\varphi }_{|\mathcal{V}|}\end{array}\right]\left[\begin{array}{c}{\sum }_{m=1}^{|\mathcal{Y}|}{\theta }_{1m}^{(0)}{\beta }_m\\ ⋮\\ {\sum }_{m=1}^{|\mathcal{Y}|}{\theta }_{|\mathcal{V}|m}^{(0)}{\beta }_m\end{array}\right]=\sum _{l=1}^{|\mathcal{V}|}{\varphi }_l\sum _{m=1}^{|\mathcal{Y}|}{\theta }_{lm}^{(0)}{\beta }_m.\end{array}$$

$${\theta }^{(1)}⟵\sum _{n=1}^N\frac{{\theta }^{(0)}\odot {\mathbf{f}}_{V[n]}{\mathbf{b}}_{Y[n]}^T}{{\mathbf{f}}_{V[n]}^T{\theta }^{(0)}{\mathbf{b}}_{Y[n]}}$$

Full text

Keywords: Bayesian Networks, Belief Propagation, Factor Graphs, Latent Variable, Optimization

1 Introduction

The Factor Graph (FG) representation, and in particular the so-called Normal Form (FGn) (Forney, 2001; Loeliger, 2004), is a very appealing formulation to visualize and manipulate Bayesian graphs; representing their relative joint probability by assigning variables to arcs and functions (or factors) to nodes. Furthermore, in the Factor Graph in Reduced Normal Form (FGrn), through the use of replicator units (or equal constraints), the graph is reduced to an architecture in which each variable is connected to two factors at most (Palmieri, 2016); with belief messages that flow bidirectionally into the network. This paradigm has demonstrated its extensive modularity and flexibility in managing variables of different types and cardinalities (Palmieri and Buonanno, 2015), and can also be used to build multi-layer network (Palmieri and Buonanno, 2014; Buonanno and Palmieri, 2015b, c).

In a previous work, a Simulink library for the rapid prototyping of an FGrn network was already implemented (Buonanno and Palmieri, 2015a), but because of the limitations imposed by Simulink, it is not particularly suitable for treating large amounts of data, and/or complex architectures with many variables.

Truthfully, despite the large presence of such arguments in the literature (Koller and Friedman, 2009; Barber, 2012; Murphy, 2012), this type of structure always suffers from high computational and memory costs, due to the lack of attention given to the specific algorithmic implementation. Therefore, this work aims to improve complexity efficiency in both inference and learning (overcoming software limitations), and all the solutions obtained have been included in a C++ library (https://github.com/mlunicampania/FGrnLib). This manuscript represents an extended version of the work presented at IEEE International Workshop on Machine Learning for Signal Processing 2018 (Di Gennaro et al., 2018).

The various problems, related to the propagation and learning of probabilities within the FGrn paradigm, are addressed by focusing on the implementation of the Latent Variable Model (LVM) (Bishop, 1999; Murphy, 2012), also called Autoclass (Cheeseman and Stutz, 1996). LVMs can be used in a large number of applications and can be seen as a basic building block for more complex architectures.

After a brief introduction of the FGrn paradigm in section 2 and the LVM in section 3, necessary to provide the fundamental elements for subsequent discussions, the C++ library project is described in section 4. In section 5 a detailed analysis of the computational complexity of the various algorithmic elements is presented for each bulding block. In section 6 some simulation results that verify how the proposeed algorithms produce indisputable advantages are presented. Finally, in section 7, an incremental learning algorithm is introduced by modifying the ML recursions, that not only presents a significant decrease in terms of memory costs (allowing learning even in the presence of huge datasets), but shows in some cases better performance in avoiding local minima traps.

2 Factor Graphs in Reduced Normal Form

A FGrn requires only the combination of the elements shown in Figure 1:

1. (a)

The Variable $V$, which can take a single discrete value $v\in{\cal V}=\{v_{1},...,v_{|{\cal V}|}\}$, is represented as an oriented edge with two discrete messages, proportional ($\propto$) to distributions, that travel in both directions.

Depending on the direction assigned to the variable, the two messages are respectively called forward $f_{V}(v)$ and backward $b_{V}(v)$, and can be also represented as $|{\cal V}|$-dimensional column vectors: ${\bf f}_{V}$ and ${\bf b}_{V}$. Note that the Marginal Distribution $p_{V}(v)$, which is proportional to the posterior given the observations anywhere else in the network, is proportional to the product

[TABLE]

or in vector form

[TABLE]

where $\odot$ denotes the Hadamard (element-by-element) product. 2. (b)

The Replicator Block (Diverter) represents the equality constraint of the connected variables. The constraint of equality between the variables is obtained by making sure that each replica can carry different forward and backward messages. A replicator acts like a bus with messages combined and diverted towards the connected branches. The combination rule (product rule) is such that outgoing messages are the product of all the incoming ones, except for the one that belongs to the same variable

[TABLE]

or in vector form

[TABLE] 3. (c)

The SISO block, that is the core of the FGrn paradigm, represents the conditional probability matrix $P(Y|V)$ of $Y$ given $V$. Assuming that the output variable $Y$ takes values in the alphabet ${\cal Y}=\{y_{1},...,y_{|{\cal Y}|}\}$, this probability matrix is a $|{\cal V}|\times|{\cal Y}|$ row-stochastic matrix

[TABLE]

or more explicitly

[TABLE]

Outgoing messages are

[TABLE]

or in vector form

[TABLE] 4. (d)

The Source block defines an independent $|{\cal V}|$-dimensional source variable $V$ with its prior distribution $\pi_{V}$. Therefore, the outgoing message is

[TABLE]

or in vector form

[TABLE]

It should be emphasized that the rules presented are a rigorous translation of the total probability theorem and Bayes’ rule.

Note that the only parameters that need to be learned during the training phase are the matrices inside the SISO blocks and the priors inside the Sources. Although different variations are possible (Koller and Friedman, 2009; Barber, 2012; Palmieri, 2016), training algorithms are derived mainly as maximization of the likelihood on the observed variables that can be anywhere in the network. Furthermore, within the FGrn paradigm, learning takes place locally; that is, the parameters inside the SISO blocks and the Sources can be learned using only the backward and forward messages to that particular element. This also means that parameter learning in this representation can be addressed in a unified way because we can use a single rule to train any SISO or Source block in the system, simultaneously and independently of its position. Therefore, learning is done iterating over three simple steps:

present the observations to the network in the form of distributions; they can be anywhere in the network as backward or forward messages; 2. 2.

propagate the messages in the whole network, in accordance with the mathematical rules just described; 3. 3.

perform the update of SISO blocks and Sources using incoming messages.

The prior of a source is learned simply by calculating the new marginal probability (using Equation 1), due to the changes of the backward message to the Source. On the other hand, learning the matrices inside the SISO blocks is more complex. According to our experience (Palmieri, 2016), the best algorithm for learning the conditional probability matrices inside the SISO blocks is the Maximum Likelihood (ML) algorithm, which uses the following update

[TABLE]

Equation 4 represents the heart of the ML algorithm and usually requires multiple cycles in order to achieve convergence. However, changing the matrices also changes the propagated messages. For this reason, the whole learning process (starting from the presentation of the evidence) also needs to be performed several times; namely for a fixed number of epochs (an hyperparameter).

3 Latent Variable Models

Factor Graphs and in particular FGrn can be used to represent a joint probability distribution by appropriately using the dependencies/independence between the variables. In many cases, however, the probabilistic description of an event can be further simplified using a set of unobserved variables, which represent the unknown “causes” generating the event. In these models the variables are typically divided into Visible and Hidden: the first relating to the inputs observed during the inference, the latter belonging to the internal representation. Obviously, the main advantage of this type of model lies in the fact that it has fewer parameters than the complete model, producing a sort of compressed representation of the observed data (which are therefore easier to manage and analyse).

The simplest and most complete model that includes both visible and hidden variables is the bipartite graph of Figure 2, named Latent Variable Model (LVM) (Murphy, 2012; Bishop, 1999) or Autoclass (Cheeseman and Stutz, 1996); any other hidden variables model can ultimately be reduced to such a model. Within the figure, it is possible to distinguish $H$ latent variables, $S_{1},\ldots,S_{H}$, and $N$ observed variables, $Y_{1},\ldots,Y_{N}$; where typically $N\gg H$. It should be noted that although the given nomenclature seems to subdivide the variables according to their position, where the variables below are known while those above need to be estimated, the bidirectional structure of the network remains quite general, including cases in which some of the above variables may be known and some of the bottom variables need to be estimated.

The FGrn of a Bayesian network represented by the general LVM structure of Figure 2 is shown in Figure 3. The system is intrinsically represented as a generative model; that is, choosing to direct the variables downwards and positioning the marginally independent sources $S_{1},...,S_{H}$ at the top.

Moreover, it should be noted that the supposed independence of all known variables given the hidden variables (provided by the latent variable model) allows to greatly simplify the analysis of the total joint probability, which can now be represented by the factorization

[TABLE]

As said, the structure allows to represent totally heterogeneous variables, but it is good to clarify that (being a representation of a Bayesian network) the variables must be presented as probability vectors; that is, through vectors containing the “degree of similarity” of the single variable to each particular element of its discrete alphabet.

The source variables, which have prior distributions $\pi_{S_{1}}$, … $\pi_{S_{H}}$, are mapped to the product space ${\cal P}$, of dimensions $|{\cal P}|=|{\cal S}_{1}|\times\cdots\times|{\cal S}_{H}|$, via the fixed row-stochastic matrices (shaded blocks in Figure 3)

[TABLE]

where $\otimes$ denotes the Kronecker product, $1_{K}$ is a $K$-dimensional column vector with all ones, and $I_{K}$ is the $K\times K$ identity matrix (Palmieri, 2016). The conditional probability matrix is such that each variable contributes to the product space with its value, and it is uniform on the components that compete to the other source variables. This is the FGrn counterpart of the Junction Tree reduction procedure because it is equivalent to “marry the parents” in Bayesian Graphs (Koller and Friedman, 2009), but here there are explicit branches for the product space variable.

For this reason, although the messages traveling bi-directionally and the initial Bayesian network present many loops (Figure 2), the FGrn architecture will not show any convergence problem because the LVM has been reduced to a tree.

Finally, the $j$-th SISO block at the bottom of Figure 3, with $j=1,...,N$, represents the conditional probability matrices $P(Y_{j}|S_{1}S_{2}...S_{H})$, which than will have dimensions $|{\cal P}|\times|{\cal Y}_{j}|$.

3.1 The LVM with One Hidden Variable

When $H>1$ we have a Many-To-Many LVM model, which was already discussed in a previous work (Palmieri and Buonanno, 2015); calling it Discrete Independent Component Analysis (DICA) because it uses the same generative model of the Independent Component Analysis (ICA), but on discrete variables.

Vice versa, when $H=1$, we obtain a One-To-Many LVM model, in which there is just one latent factor (parent) that conditions $Y_{1},\ldots,Y_{N}$ (children) and where obviously fixed matrices (previously represented by the shaded blocks) are no longer necessary. Although the general paradigm has the great advantage of allowing the presence of Sources of different types, for simplicity in the tests performed we have preferred to focus exclusively on this more manageable architecture (Figure 4). The figure shows the most general case (used in the final tests), in which a Class Variable $L$ is also added to the bottom variables. This configuration is used in the case of supervised learning, allowing (after learning) to perform various tasks, including:

1. (a)

Pattern classification, achievable by injecting the observations as delta distributions on the backwards of $Y_{1},...,Y_{N}$ and leaving the backward on L uniform. The classification will be returned to $L$ through its forward message $f_{L}$. 2. (b)

Pattern completion, where only some of the observed variables $Y_{1},...,Y_{N}$ are available in inference (and then injected into the network through delta distributions on their backwards) while the others are unknown (and therefore represented by uniform backward distributions).

Also $L$ may be unknown (uniform distribution), partially known (generic distribution), or known (delta distribution). The forward distributions at the unknown variables complete the pattern, and at $L$ provides the best inference in case it is not perfectly known. The posterior on $L$ is obtained by multiplying the two messages on it (Equation 1); avoidable step if $L$ is not known at the beginning (because in this case the backward is uniform). 3. (c)

Prototype inspection, obtainable by injecting only a delta distribution at $L$ on the $j$th label. The forward distributions $f_{Y_{1}},...,f_{Y_{N}}$ will represent the separable prototype distribution of that class.

Another way to use the previous network is to make the class variable coincide with the hidden variable $S=L$, forcing the corresponding SISO block matrix to be diagonal. This constraint will create a so-called “Naive Bayes Classifier” (Barber, 2012), further simplifying the factorization into

[TABLE]

In this case, usually all the variables are observed during training, and the typical use in inference is to obtain $L$ from observed $Y_{1},...,Y_{N}$ .

Note that the case related to unsupervised learning can be obtained from the general model presented simply by eliminating the variable $L$. In this case, after learning, the elements of the alphabet ${\cal S}=\{s_{1},...,s_{|{\cal S}|}\}$ of the hidden variable $S$ represent “Bayesian clusters” of the data, which follows the prior distribution $\pi_{S}$ (learned blindly). The network can be used in inference both for the pattern completion, in the case where (as seen previously) only some of the underlying variables are known and we try to estimate the others through the corresponding forward messages, and to create a so-called embedding representation, in which the backward message becomes a different representation of the underlying known variables. In the latter case, in order to understand the representation that the network has created, we can look at the j-th centroid of the bayesian clusters injecting as ${\bf f}_{S}$ a delta distribution ${\bf\delta}_{j}=[0\ldots 1\ldots 0]^{T}$, where the $1$ is at the j-th position. The set of forward distributions $f_{Y_{1}},...,f_{Y_{N}}$ generated by the network will represent the marginal distributions around the centroid of the $j$th cluster.

4 Design of FGrnLib

There are several software packages, known in the literature, that can be used to design and/or simulate Bayesian networks (an updated list can be found in (Murphy, 2014)). Unfortunately, many of them are in closed packages and/or run only on private servers, preventing proper performance analysis. Others either have limitations on the number of variables and the size of the network, or do not use the FG architecture. Therefore, the main purpose of this work was to design an optimized library, called FGrnLib, for the realization of a Bayesian network through the use of the FGrn model; that is open and contains an efficient implementation of the elements in Figure 1.

The FGrnLib library has been written in C++, following the classic object-oriented paradigm, and it has been adapted for parallel computing (on multiprocessor systems with shared memory) through the use of the OpenMP application interface.

The various algorithmic operations have been implemented to limit as much as possible the computational complexity without significantly affecting memory requirements.

4.1 Data Structures

Before starting the analysis of individual operations, it is necessary to focus on the structure of the main classes (Figure 5) and their individual roles. These classes correspond to the main elements presented in Figure 1 (and mathematically described in section 2), but also have subtle design choices that need to be clarified (especially in reference to what was previously done in (Buonanno and Palmieri, 2015a)). The main classes are:

• •

The Link class, which represents a single discrete variable of the model. This class contains the two forward and backward messages for the variable, and is designed to ensure that each message at each time step is represented by a single vector. This means that in every instant of time every single link takes on only two messages (in the two different directions); providing better control of the information traveling on the network but preventing the possibility of carrying out stages of learning through the simultaneous presentation of all the evidence.

• •

The Diverter class, which imposes the constraint of equality among the variables. This class has been created to be as general as possible, in the sense that it can automatically adapt the parameters of the net to the number of variables. For this reason, it includes not only the process of replication of variables but also the creation and control of product space matrices (Equation 5). For space complexity, these (sparse and row-stochastic) matrices are stored in column vectors, whose elements represent the index of the active column in that particular row of the matrix

[TABLE]

.

• •

The SISOBlock class, which represents the probability of the output variable $Y$ given an input variable $V$, contains the row-stochastic conditional probability matrix $P(Y|V)$. Due to the fact that the Link class allows only two vector at a time, the SISO blocks are realized to permit the storage and retrieval of all messages that reached the blocks during the batch learning phase; avoiding transmitting the evidences all at the same time but reusing the single Link memory units for each step.

• •

The Source class, which represents the independent variables inside the model, with their prior probabilities.

4.2 Data Flow Management

Since the FGrn paradigm allows to work only through local information, the flow of messages within the network can be parallelized. However, parallelizing the message flow in the network imposes essential changes to the Diverter class. In fact, the multiplication of the messages internally to the Diverter can take place only after all the messages of the variables relating thereto have been updated. Being the only responsible for the combination of messages coming from different directions, the Diverter must necessarily act as a “barrier” (with the meaning that this term assumes in parallel programming); solving the synchronization problems simply by not updating the output values until it receives the activation signal from a supervisor. Although within the FGrnLib a Supervisor class has been defined to more easily manage some predefined network, the meaning of a supervisor is quite general here because it refers to any class that possesses all the references to the single elements of the network realized.

In Figure 6 it is shown how the supervisor handles the message scheduling, relatively to both the inference mode and the batch learning mode. Note that, in order for messages that travel in parallel to be propagated anywhere on the network, a number of steps equal to the diameter of the graph is required (Pearl, 1988).

Recalling therefore that, in graph theory, the diameter of the network is the longest possible path (in terms of number of crossed arcs) among the smallest existing paths able to connect each vertex inside the graph, it is easy to understand that a simple one-layer LVM network (like in Figure 3 or in Figure 4.a) will only need three propagation steps.

4.2.1 Inference

The inference phase, depicted in Figure 6.a, is relatively simple: the various messages proceed in parallel until they reach the Diverter. All that is necessary is to block the messages at the Diverter, to prevent the start of the multiplication process before all the messages have become available. The supervisor must therefore perform only two phases: parallelize the input variables and then activate the Diverter. It should be noted that typically all messages are initialized to uniform vectors, and updated according to the rules described in section 2.

For what has been said, in the first step of Figure 6.a the supervisor performs the parallelization of the input variables $Y_{1},\dots,Y_{N}$, placing in the backward messages appropriate distributions for the known (or partially known) variables and uniform distributions for the unknown ones. Hence, in the second step, the messages start to be propagated in the network (using the second in Equation 3). At this point, when all the variables have performed the second step, the supervisor must activate the Diverter to allow the propagation of the messages to continue. The last message propagation step allows to obtain the desired output values (using the first in Equation 3).

4.2.2 Batch Learning

Every single epoch of the batch learning phase (Figure 6.b) is not so different from what we have just seen for the inference phase, since the supervisor basically performs only two more operations. The first one, which only applies at the beginning of the whole learning phase (that is, it is not reiterated at each epoch), consists in activating the batch learning mode inside the SISO blocks and the Sources. This step enables storage within the SISO blocks, which consequently memorize the incoming vectors (from both sides) as soon as they are updated, and the Sources, which will begin to add together all incoming backward messages within a temporary vector. It should be noted that, in the case of multi-layer structures, the supervisor will also have to worry about preventing the SISO blocks of the next level from propagating incoming messages, thus transforming the latter into “barriers” as well. In fact, in order to get the classic layer-by-layer approach, the information should not propagate above the Diverter until the underlying batch learning phase is complete, but only if the top layer does not provide Sources. So as not to complicate the Diverter too much, forcing it to know the connected elements, the SISO blocks also provide a pause mode, which if enabled prevents forward propagation of messages.

After activation of the batch learning mode the messages propagate within the network through the same previous steps, being blocked again to the Diverter. However, as we can see in Figure 6.b, the descending phase will not include the production of the final messages, which will only be stored in the SISO blocks in order to avoid performing the related mathematical operations.

The second different operation, performed by the supervisor at the end of each epoch, consists in the activation of the actual learning procedure, which will execute the ML algorithm using the stored vectors (up to a sufficient convergence or at most for a fixed number of times) and change the prior of the Source. After this last phase, the procedure can be repeated for another epoch, presenting the evidence to the network again. When it will be necessary to end the learning, the last message of the supervisor will also modify the functioning modalities of the SISO blocks and the Sources, thus making them ready for the subsequent phases of inference.

5 Complexity and Efficient Algorithms

Having to pay particular attention to the computational and memory costs, in the creation of the library we worked on the details of each individual element. This, together with the probabilistic nature of the Bayesian networks, has led to the preliminary definition of particular basic data structures that it is perhaps necessary to analyse quickly. In fact, the library also defines classes that represent probability vectors and row-stochastic matrices, to facilitate the interpretation and definition of the variables and to easily manage all the algebraic operations.

5.1 Probability Vector

A probability vector is a vector with elements in $[0,1]$ that sum to one. Although the network does not use only probability vectors, the execution of every operation that uses them must then provide normalization. Every normalization consists of $d-1$ sums and $d$ divisions, so the computational cost will be $\mathcal{O}(d)$

[TABLE]

5.2 Row-stochastic Matrix Multiplication

In a row-stochastic matrix $P$ each row is a probability vector. It is important to observe that the premultiplication or postmultiplication of a row-stochastic matrix for a vector (with the appropriate dimensions)

[TABLE]

will consist of $ld$ multiplications and $(l-1)d$ or $l(d-1)$ sums respectively, producing the same computational cost $\mathcal{O}(ld)$.

5.3 Diverter

From a computational point of view, the most critical structure in the implementation of a Bayesian network using FGrn is represented by the Diverter, where, obviously, the greatest criticality is in the efficient implementation of the internal multiplication process (Equation 2).

In fact, the possible presence of zeros, in different positions of the single vectors, obliges to perform the calculation of the outgoing vectors individually. In the general case of Figure 3, indicating the total value of all the variables connected to the diverter with $M=H+N$ and assuming that the input vectors have dimensions equal to $d$, bare application of the multiplication rule would require $M(M-2)d$ multiplications. Regardless of the size of the input vectors, the computational cost is therefore polynomial and equal to $\mathcal{O}(M^{2})$ vectorial multiplications. Considering as an example the Diverter of Figure 7 (with $M=4$), the simple application of the product rule in the production of the outgoing messages $({\bf b}_{V^{(0)}},{\bf f}_{V^{(1)}},{\bf f}_{V^{(2)}},{\bf f}_{V^{(3)}})$ would requires eight vectorial multiplications:

[TABLE]

This process can be performed more efficiently by defining an order among the variables connected to the Diverter and by performing a double cascade process in which each variable is responsible only for passing the correct value to the neighboring variable. In this way, the variables at the ends of the chain will perform no multiplication while each variable inside the chain will perform only three multiplications, relative to the passage of the two temporary vectors along the chain and to the output of the outgoing message.

With reference to the example of Figure 7 we have the data flow represented in the Figure 8. In other words, the proposed solution exploits the presence of the same multiplication groups through a round-trip process. This reduces the computational complexity from quadratic to linear, $\mathcal{O}(M)$, finally requiring only $3(M-2)$ vector multiplications. Although there is obviously an increase in the memory required, due to temporary vectors along the chain, it remains linear with $M$, and has been further optimized (requiring only $M-1$ temporary vectors altogether) by choosing to reuse the same vectors ($a_{i}=a_{i}^{{}^{\prime}}$) when changing direction.

5.4 Unknown variables

A very attractive property of using a probabilistic paradigm, which makes it preferable in certain contexts, is represented by its ability to manage unknown inputs.

Even in the case of maximum uncertainty, that is when nothing is known about the particular variable in that particular observation set, a process of inference or learning can still be performed by making the corresponding message of that variable a uniformly distributed probability vector

[TABLE]

In this particular circumstance, the message propagation process can be optimized by avoiding the multiplication of backward vectors with the matrices inside the SISO blocks, noting that

[TABLE]

By not propagating the unknown variable (setting ${\bf b}_{S^{(i)}}$ as an expanded/reduced version of $\bar{\bf b}_{Y_{i}}$), for every single unknown variable present in input during the inferential and the learning process we can save $|{\cal S}|$ vector multiplications, improving overall network performance.

5.5 Efficient ML Implementation

Regarding the learning phase, particular attention has been given to the realization of an efficient implementation of the ML algorithm. First of all, since the matrix inside the SISO blocks is set to be row-stochastic by construction, it has been noted that the first divisor in the Equation 4 becomes unnecessary. For this reason, the equation can be rewritten as follows

[TABLE]

or in vector form

[TABLE]

Furthermore, it can be observed that the value obtained through the vector multiplications ${\bf f}_{V[n]}^{T}{\bf\theta}^{(0)}{\bf b}_{Y[n]}$ is actually equal to the sum of all elements of the matrix ${\bf\theta}^{0}\odot{\bf f}_{V[n]}{\bf b}_{Y[n]}^{T}$. This assertion is provable by observing that

[TABLE]

and realizing that the sum of all the elements of the matrix can be written in the form

[TABLE]

which precisely is equal to

[TABLE]

This suggests that moving the Hadamard product inside the summation

[TABLE]

and calculating the sum of all the elements of the matrix ${\bf\theta}^{0}\odot{\bf f}_{V[n]}{\bf b}_{Y[n]}^{T}$ at the same time of their generation, computational complexity of the algorithm can overall be reduced from $N(3|{\cal Y}|+1)|{\cal V}|$ to $2N|{\cal V}||{\cal Y}|$ multiplications, which is rather significant if considering that this reduction is relative to each algorithm recall.

6 Performance on LVM

To evaluate the computational advantages obtained by the proposed improvements let us consider the more general situation of Figure 3, in which the LVM (single-layer) model foresees $H$ sources (hidden) and $N$ output (observed) variables $Y_{1},...,Y_{N}$. Since the variables $Y_{1},...,Y_{N}$ can have different dimensions (never less than $2$), before starting it is important to underline that we will assume that all the observed variables have the same dimension $|{\cal Y}|$. This statement does not compromise in any way the goodness of the results obtained, for we can certainly decide to choose $|{\cal Y}|$ as the highest value possible being interested only in an upper-bound description of computational complexity (in terms of big-O notation). For the same reason, we will avoid considering the case (more advantageous) in which some variables are not know. Note that, the previous problem does not arise in the case of input variables to the Diverter because they already have the same dimension $|{\cal P}|$ by construction.

6.1 Cost of the inference phase

In the inference mode, at the beginning of the process, the backward messages of the $N$ output variables will be post-multiplied by the probabilistic matrix inside the SISO blocks; thus producing $\mathcal{O}(N|{\cal P}||{\cal Y}|)$ operations before being sent to the Diverter. As seen previously, $\mathcal{O}((H+N)|{\cal P}|)$ operations will be performed inside the Diverter, related to the multiplication of all the $H+N$ incoming messages to it. Finally, the messages must be propagated again to the SISO blocks, which will be pre-multiplied by the matrices still producing $\mathcal{O}(N|{\cal P}||{\cal Y}|)$ operations, and thus making the total computational cost equal to $\mathcal{O}((N|{\cal Y}|+H)|{\cal P}|)$. At this point the forward messages of the $Y_{1},\dots,Y_{N}$ output variables will be available for analysis. The following table shows the difference in computational terms determined by the introduced optimizations

[TABLE]

To provide a practical example, suppose to perform an inferential process on an LVM network with $10$ binary output variables and a single hidden variable, with an embedding space $|{\cal S}|$ of $10$. It will thus be noted that while the direct algorithm will require 1500 multiplications, the optimized one will require only 670.

Moreover, the additional memory required by the optimized algorithm is still linear with the size of the inputs, since it depends only on the temporary vectors inside the Diverter. In fact, in the previous example, the significant computational advantage obtained will correspond to an increase in memory equal to only $10$ vectors of size $10$.

6.2 Cost of the batch learning phase

For what concerns the batch learning session, assume to have $L$ different examples given in input through the backward messages of the output variables. As already mentioned, when the process begins backward messages entering the network will be saved inside the SISO blocks, requesting an amount of memory of $\mathcal{O}(L|{\cal Y}|)$ individually. After storing the vector, the process will continue by multiplying it and the probability matrix; sending the result to the Diverter. In the propagation phase the computational costs will not change with respect to what has been seen in the inference phases, but it must be remembered that the messages that will return back to the SISO blocks will be memorized again; with memory costs individually equal to $\mathcal{O}(L|{\cal P}|)$. The total cost of additional memory required is therefore $\mathcal{O}(LN(|{\cal Y}|+|{\cal P}|))$, while the computational cost is still $\mathcal{O}(L(N|{\cal Y}|+H)|{\cal P}|)$.

Once this first phase has been completed, the ML algorithm will be executed at most K times (with K fixed a priori), trying to make conditional probability matrices converge, and the whole process is then repeated $T$ times. Thus, the total computational cost of the batch learning session is equal to

[TABLE]

Again, a table is presented which facilitates the comparison between the computational costs of the direct case and those of the optimized case.

[TABLE]

Furthermore, it is easy to state that the repetition of the process for $T$ epochs, as well as the various calls of the ML algorithm, do not imply the need for any additional memory units with respect to the non-optimized case.

7 Incremental Algorithm

As noted above, the ML algorithm has many undoubted advantages, being very stable and typically converging in a few steps. Unfortunately it is a batch-type algorithm that is able to perform only on the whole training set. In order to obtain a lighter implementation we have changed the previous structure by requiring that at each epoch of the learning phase only one ML cycle (that is, $K=1$) is included. In other words, the algorithm has been made incremental obtaining the advantage of reducing the amount of memory. In fact, this approach makes it unnecessary to store backward messages within SISO blocks, since they must now be used only once. This eliminates the need for the previous storage space equal to $L(|{\cal P}|+|{\cal Y}|)$ for each SISO block.

Despite the great advantage in terms of both memory and computational costs, this type of approach has surprisingly proved to be as robust as the previous one, being less likely to provide overfitting of the data. Referring to a One-To-Many LVM structure (Figure 4.b), in which the only hidden variable has an embedding space $|{\cal S}|$ of $20$, the following tables show the accuracy of the classification for three different datasets.

[TABLE]

• Datasets composition: Breast Cancer = 10 useful variables, 699 instances, 16 missing values; Mammographic Mass = 6 variables, 961 instances, 162 missing values; Contraceptive Method = 10 variables, 1473 instances, 0 missing values.

The tables present the classification success rates, both on the training and on the test set, of three databases from the UCI repository: Wisconsin Breast Cancer Dataset (Dheeru and Taniskidou, 2017), Mammographic Mass Data Set (Elter et al., 2007) and Contraceptive Method Choice Data Set (Dheeru and Taniskidou, 2017). The values presented are obtained by making sure that the latent variable is not learned, and therefore represent the learning ability of a single layer; being a good index to represent layer-by-layer learning of more complex networks. Note that, in this particular situation, the incremental algorithm does not give results that deviate much from the values obtained with the batch learning, providing in some cases even better results. It should also be noted that the results do not improve according to the adaptability of the paradigm to the specific case, since they are better even when the realized one-layer LVM network is probably not suitable for capturing the underlying implication scheme (as seen in the case of Contraceptive Method database).

Finally, the results prove even more interesting if we consider that in both cases (batch and incremental) they are obtained using the same number of epochs (in particular equal to $20$). In fact, an incremental algorithm should typically employ many more steps to achieve the performance of a batch algorithm, whereas in the cases under examination the change seems to bring only advantages.

8 Discussion and Conclusions

In this work, an in-depth analysis of the individual elements necessary to create a Bayesian network through the FGrn paradigm was conducted, showing how it is possible to reduce memory and computational costs during implementation. The analysis led to the creation of a C++ library able to provide excellent results from a computational point of view, transforming polynomial costs into linear (respectively to the number of variables involved). The incremental use of the ML algorithm has finally demonstrated how it is further possible to reduce both the computational and memory costs of the learning phase, even improving, in some of the cases, the ability of the network to learn from the evidence. All these algorithmic choices are at the basis for extending the FGrn paradigm’s to higher-scale problems.