S U M M A R Y

For any Internet service provider or network operator, it is crucial to quickly and efficiently diagnose the problems that occur on the network. The benefits of a good fault diagnosis system are mainly to minimize the costs of network and service operations and to enhance the customer's quality of experience. One major challenge for any diagnosis system concerns the discovery of new faults, that are unknown to the current version of the diagnosis system. The exploratory process for finding new faults can prove to be expensive and time consuming for internet service providers.

In this thesis, we explore an alternative approach based on learning methods, in order to build learning-based diagnosis systems. Our study explores Probabilistic Graphical Models that are capable of clustering patterns of faults in an unsupervised and a semi-supervised manner. We demonstrate the efficiency of our models on real use-cases of large scale data, extracted from Fiber-to-the Home (FTTH) services based on Gigabit-capable Passive Optical Networks.

Our contributions are as follows:

• We propose the Infinite Categorical Mixture Model (Echraibi et al., 2020b) for unsupervised clustering of faults from categorical network data, with simultaneous learning of the number of clusters. We show that our model achieves competitive results, and allows for the identification and interpretation of clusters of faults.

• We propose the Infinite Mixed-type semi-supervised Mixture Model, extending the categorical model to treat continuous and categorical data. We also introduce semi-supervision in the form of labels provided by domain experts to guide the clustering process towards relevant fault clusters.

• We propose the Deep Infinite Mixture Model for semi-supervised discovery of faults from noisy large scale data of GPON-FTTH networks. The model leverages the power of deep neural networks to learn hidden representations relevant to the fault diagnosis task, and the exploration and clustering capability of infinite mixtures to identify clusters of faults.

• We introduce a novel approach to train deep probabilistic graphical models, namely Generalized Stochastic Backpropagation (Echraibi et al., 2020a). Our approach gives a theoretical framework for computing low variance estimates of gradients often arising in probabilistic graphical models.

• We introduce a general end-to-end infinite deep probabilistic mixture model, entitled the Dirichlet Process Deep Latent Mixture Model (Echraibi et al., 2020c). We demonstrate that the model is capable of learning representative hidden features of clusters of faults, while simultaneously clustering them without a prior knowledge on the number of clusters.

• Finally, we demonstrate that shallow infinite mixture models can prove interpretable, in the sense where identified clusters can be traced to network variables (or root-causes). We also give a general method of interpreting clustering results of the deep models, using standard approaches such as decision trees.

L I S T O F F I G U R E S

Summary of the structure of the document . . . . . . Figure 2 An example of a PGM. . . . . . . . . . . . . . . . . . . Figure 3 A generative model of an image. . . . . . . . . . . . . Figure 4 Graphical representation of a mixture model. . . . . . Figure 5 Graphical representation of the model in plate notation. Figure 6 Expert Bayesian Network of the FTTH GPON [START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. Colored nodes are possible root causes whereas non-colored nodes are observations. . . . . . Figure 7 Network scope for fault pattern extraction from real operational data . . . . . . . . . . . . . . . . . . . . . . Figure 8 evidence lower bound of the model across variational inference iterations, mean and standard deviation over 10 runs. . . . . . . . . . . . . . . . . . . . . . . . . Figure 9 Convergence time as a function of the number of instances for the variational inference and Gibbs sampling approaches. . . . . . . . . . . . . . . . . . . . . . Figure 10 Comparison of the clustering accuracy across iterations for the variational inference and Gibbs sampling approaches. . . . . . . . . . . . . . . . . . . . . . . . . . Figure 11 Assignments defining each cluster from most occurring to least occurring . . . . . . . . . . . . . . . . . . . Figure 12 Graphical representation of the mixed-type infinite mixture model in plate notation. . . . . . . . . . . . . Figure 13 semi-supervised evidence lower bound across variational inference iterations. . . . . . . . . . . . . . . . . Figure 14 Accuracy of the model across variational inference iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 15 Overview of the diagnosis system. . . . . . . . . . . . Figure 16 Architecture of the feature extractor followed by the softmax unit. . . . . . . . . . . . . . . . . . . . . . . . . Figure 17 2-D representation of the feature manifold across epochs of training. . . . . . . . . . . . . . . . . . . . . . Figure 18 (Left) 2-D representation in the feature space. (Right) 2-D representation of the raw data. . . . . . . . . . . . Figure 19 Loss function on the test set across iterations of gradient descent. Accuracy on the test set across epochs.

The clustering process of the infinite Gaussian mixture model on the latent feature manifold across iterations of variational inference. . . . . . . . . . . . . . . Figure 20 Log probability of mixture across iterations. . . . . . . Figure 22 Training loss and log variance of the gradients for the different estimators for f (z) = ∑ d j=1 (z j -) 2 for d ∈ {1, 10, 100}. . . . . . . . . . . . . . . . . . . . . . . Figure 23 Training loss and log variance of the gradients for the different estimators for f (z) = ∑ d j=1 exp(-z j ) for d ∈ {1, 10, 100}. . . . . . . . . . . . . . . . . . . . . . . Figure 24 Bayesian Logistic Regression with Laplacian priors . Figure 25 (Left) Bias and variance of the gradient for different values of the truncation level at a fixed parameter value. (Right-top) Mean square error between the Laplace gradient estimator and the score function and reparameterization estimators across iterations.

(Right-bottom) Norm of the gradient estimators. . . Figure 26 The graphical representation of the generative process of the model, with the convention x = h (0) . . . . . . . Figure 27 (Left) t-SNE plot of the second stochastic hidden layer on the MNIST test set for the semi-supervised (10% labels) version of the DP-DLM. (Right) t-SNE plot of the second stochastic hidden layer on the MNIST test set for the unsupervised version of the DP-DLM. . . . Figure 28 Generated samples from the DP-DLM model for the unsupervised version on the MNIST dataset (left) and the semi-supervised version on the SVHN dataset (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 29 (Left) Evidence lower bound across learning iterations. (Right) Comparison of the accuracy between the Semi-supervised mixed type mixture and the Dirichlet process deep latent mixture. . . . . . . . . . . . . . Figure 30 2D representation of the last hidden layer for each class, and clusters discovered. . . . . . . . . . . . . . . Figure 31 Interpretation of fault cluster 40, corresponding to an account suspension of the customer on the network. . Figure 32 Posterior distributions of the visible variables given the cluster assignment. . . . . . . . . . . . . . . . . . . Figure 33 An example of a decision rule leading to cluster 4. . . Figure 34 Clear decision leading to cluster 3. . . . . . . . . . . . Figure 35 Graphical representation of the different processes in the neural variational model. . . . . . . . . . . . . . .

L I S T O F TA B L E S

Table 1 Description of each dataset. . . . . . . . . . . . . . . . Table 2 Confusion Matrix of the Infinite Categorical Mixture Model on the synthetic dataset. . . . . . . . . . . . . . Table 3 Comparison between our model and the baselines. . . Table 4 Comparison between our model and the baselines. . . Table 5 Description of the different classes in the dataset. . . Table 6 Overview of the network variables in the dataset. . . Table 7 The confusion matrix on the known classes. . . . . . . Table 8 Semi-supervised classification error (%) on the MNIST test set with 10 % labelisation. Comparison with (Nalisnick and Smyth, 2016) . . . . . . . . . . . . . . . Table 9 Test likelihood for the Bernoulli stochastic backpropagation (BSB) estimator compared to the REBAR and RELAX estimators. We report the mean and standard deviation over 5 runs. . . . . . . . . . . . . . . . . . . . Table 10 Test likelihood for the categorical stochastic backpropagation (CSB) estimator, compared to the Gumbelsoftmax estimator. We report the mean and standard deviation over 5 runs. . . . . . . . . . . . . . . . . . . . Table 11 Execution Identifying faults in networks is of crucial importance to any Internet service provider. An accurate diagnosis system capable of identifying faults, can improve the customer's experience, decrease the down time of network operations, and more importantly reduce the costs of an Internet service provider in terms of interventions to repair the network.

A considerable amount of research has been done to construct and maintain operational fault diagnosis systems in all domains, ranging from network operations to industrial plants and applications [START_REF] Mosterman | Diagnosis of continuous valued systems in transient operating regions[END_REF][START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. However, the research scope for solving the diagnosis task is often specific to the domain, and can even vary in definition.

In the network management community, the diagnosis task is divided into three sub-tasks. The first is the detection step, which can be done proactively, it aims at deciding if a customer experiences a problem and if further investigations are required. The second is the isolation step, where the goal is to identify the root cause of the problem given the current technical status or available data. Third, the mitigation step covers all the actions required to fix the problem [START_REF] Ayoubi | Machine Learning for Cognitive Network Management[END_REF][START_REF] Boutaba | A comprehensive survey on machine learning for networking: evolution, applications and research opportunities[END_REF][START_REF] Cherrared | A Survey of Fault Management in Network Virtualization Environments: Challenges and Solutions[END_REF][START_REF] Steinder | A survey of fault localization techniques in computer networks[END_REF].

The first step, consisting in identifying if a fault occurred on the network, is often associated with detecting anomalous behavior. The methods often adopted for this task are Anomaly Detection approaches (Chandola, Banerjee, and Kumar, 2009a). For the specific application that we study across the thesis, i.e. fault diagnosis of broadband access services, these methods are not considered, since fault detection is a direct outcome of a customer call to the Internet service provider's hotline. Thus, we are rather interested in the isolation step, where we know that a fault has occurred on the network, and we must identify the cause or the nature of the fault. The scope of our application is fault diagnosis of Fiber-to-the Home (FTTH) services based on Gigabit-capable Passive Optical Networks (GPON) [START_REF] Chanclou | Overview of the optical broadband access evolution: a joint article by operators in the IST network of excellence e-Photon/One[END_REF]; Gigabit-capable passive optical networks (G-PON): ONT management and control interface specification 2008; Gigabit-capable passive optical networks (GPON): General characteristics 2008).

With the various modern applications of machine learning, the objective is to build models capable of classifying faults from GPON-FTTH network data derived from operations of Orange networks. Furthermore, the challenge is to build adaptable systems. In the sense that, if a new fault appears on the network (due to changes or other phenomena) the model is capable of clustering the new fault, given that a signature of the fault exists in the data encountered.

. In order to introduce the network diagnosis task as a machine learning problem, we need to first discuss the current status of diagnosis for an Internet service provider such as Orange, and for the research community in general. Currently, in most real-world operational diagnosis systems, the preferred method is Expert Systems. Specifically for Orange, the expert system is called DELC (Diagnostic Expert de la Ligne Client) or expert diagnostic for the customer's line.

An Expert System is a rule-based system, where given a certain deterministic condition on the network variables an alarm is triggered to alert if a fault happened. This has its advantages and its drawbacks. The first advantage is that these rules can be constructed directly given enough expert knowledge, and they are interpretable and clear to other network or domain experts. However, the drawbacks are numerous. The first if the network is large the expertise and time required to build the expert system are enormous, maintaining it is even more challenging and requires tremendous resources. Second, if a new fault occurs due to changes in the network, the system has to be revisited, which can prove even more daunting. Third, since the rules are deterministic, the expert system cannot always conclude when some variables are missing.

The first approach that has been considerably researched for the diagnosis task are Bayesian Networks [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF][START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. A Bayesian Network is a probabilistic model of the network variables with dependencies. The model can be constructed by an expert, the dependencies represent causal relations between variables of the network. The objective of a bayesian network is to model the propagation of faults through variables of the network, thus by reversing the process, we can identify the root cause. The advantage of these models is that they can be learned from data (conditional probabilities estimated from data). However, as expert systems, bayesian networks are highly specific to the application: if a new fault or a change occurs on the real network, the bayesian network has to be revisited which is often an expensive task.

Classical machine learning methods have also been explored for the diagnosis task. One approach is to consider the traditional supervised learning scheme, where a classifier is trained using data labeled according to the different known types of failures [START_REF] Adda | Comparative analysis of clustering techniques in network traffic faults classification[END_REF][START_REF] Baras | Automated network fault management[END_REF][START_REF] Chen | Pinpoint: Problem determination in large, dynamic internet services[END_REF][START_REF] Chen | Failure diagnosis using decision trees[END_REF]. These classifiers are numerous: they can range from simple models such as Naive-Bayes classifiers [START_REF] Zhang | Exploring conditions for the optimality of naive Bayes[END_REF], to Support-Vector Machines (SVMs) [START_REF] Cortes | Support-vector networks[END_REF], and even Neural Networks [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF].

Although these methods are appealing, they require a retraining of the classification model with new features and new labeled data for new types of faults, which generally are not available, especially when the fault is very recent. With the Bayesian network framework, taking into account new types of failures, new data or new equipment would mean updating the Bayesian graph by adding and/or removing nodes and edges as well as changing conditional probabilities of the edges. This would require a high network expertise and may be unrealistic if the network is too complex. Though it is possible to derive the Bayesian graph from data [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF], it becomes unfeasible in practice when dealing with large datasets composed of hundreds or thousands of features.

Furthermore, in order to discover new faults, one needs to perform some data exploration and to find clusters corresponding to new types of faults. The current standard approach is an expert investigation of the data "by hand" in order to identify these clusters. Unsupervised clustering methods can be used to reach good performance. However, using such approaches is often dependent on model assumptions that could hurt the results, such as distance metrics in high dimensional spaces [START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF][START_REF] Kriegel | Clustering high-dimensional data: A survey on subspace clustering, pattern-based clustering, and correlation clustering[END_REF]. In addition, the clusters may be relevant from a technical point of view (e.g. grouping individuals because they use the same type of equipment), but this type of clustering does not provide any new information that could be used to highlight a new fault. The obvious way to bypass this issue is to resort to cumbersome feature selection and preprocessing of the data by hand, which requires deep expert knowledge of the environment that is constantly changing, and can be expensive time-wise.

In this thesis, we present a general framework based on Probabilistic Graphical Models to learn and identify clusters of faults. We argue for a specific range of probabilistic graphical models namely, Dirichlet Process Mixtures or Dirichlet Process Deep Mixtures. These models we propose are capable of clustering faults without knowing the number of clusters (or faults), they can identify hidden features relevant to the diagnosis task, and more importantly, they include expert knowledge in terms of supervision (or semi-supervision) in the form of labels in order to guide and improve fault identification from data.

. This thesis is divided into two parts: in Part i, we introduce shallow probabilistic graphical models. Specifically, Infinite Mixture Models that are capable of identifying clusters of faults from network data, in an unsupervised and a semi-supervised way. In Part ii, we evolve these models to include deep hidden representations. The goal of the hidden representations is to construct features relevant for the diagnosis task. These features are constructed to deal with noisy and complex large scale data. Combined with the infinite mixtures, we can use these representations to reach better results on such datasets.

In order to give a guide of the manuscript, the contributions made in Part i are an attempt to answer the following research questions:

Question 1: Can we develop models capable of identifying clusters adapted to network data, specifically, categorical data, without knowing a priori the number of clusters?

Our contribution to answer this question is the Infinite Categorical Mixture Model (Echraibi et al., 2020b). The model is introduced in chapter 3. We extend the definition of Dirichlet Process mixtures to the categorical variable case and we show how to use mean-field variational inference to learn the model efficiently. We show that our model is capable of identifying clusters of faults without knowing a priori the number of clusters and we demonstrate how the model assumptions lead to better performance in terms of accuracy.

Question 2: Can we extend the Dirichlet process mixtures to deal with continuous and categorical variables simultaneously, and also take into account expert knowledge in terms of supervision?

In this contribution, we extend the infinite categorical mixture model to a Mixed-type model, that can deal with continuous and categorical variables. The model also takes into account semi-supervision from experts in the form of labels. We evaluate the model on network data and we show its competitiveness with state-of-the art baselines. This work is presented in chapter 4 and is based partially on work from [START_REF] Echraibi | Bayesian Mixture Models For Semi-Supervised Clustering[END_REF].

We have introduced shallow probabilistic graphical models, in the form of infinite mixture models, that could be applied to mixed-type data. A problem remains, these models cannot create complex combinations of data in order to discover clusters. In other words, feature creation is not possible. However, in most real-world applications, data is noisy and of large dimensionality. Thus, feature creation is necessary, in Part ii, we develop models and methods in order to accomplish this goal, and the main research question addressed are the following:

Question 3: Can we construct highly relevant representations or features for the diagnosis task, in order to improve performance?

In order to answer this question, in chapter 5, we combine deep learning representations and infinite mixture models to perform fault discovery in a partially labeled large scale dataset. Our model is called Deep Infinite Mixtures (Echraibi et al., 2021b). We train a deep neural network to recognize labeled faults and build low dimensional hidden representations relevant to the diagnosis task. Then, an infinite Gaussian mixture model is applied on the unlabeled part in order to recognize clusters of faults based on the hidden representation. We evaluate the model on a large scale real-world GPON-FTTH dataset, we show that the model can identify classes of known faults with high accuracy, and that the hidden representation on the unlabeled faults lead to recognizable clusters.

Deep learning and infinite mixtures live in two different universes of methods. Their combination is not natural in the sense we do not have an end-to-end global model, capable of learning features and clustering. Thus, the following question naturally emerges:

Question 4: Can we construct an end-to-end deep probabilistic model capable of learning hidden representations and of performing clustering? It turns out this is possible through deep generative models, which are a class of probabilistic graphical models. The problem of deep generative models is the training process that involves computing estimates of gradients that are often of large variance. In our contributions Generalized Stochastic Backpropagation, and Stochastic Backpropagation through Fourier Transforms (Echraibi et al., 2020a(Echraibi et al., , 2021a)), we present a method to derive low-variance estimates of these gradients. This work is presented in chapter 6.

In chapter 7, we use generalized stochastic backpropagation to train our proposed end-to-end Dirichlet Process Deep Latent Mixture Model (Echraibi et al., 2020c). We evaluate our model on open datasets and on real-world large scale GPON-FTTH data for the fault diagnosis and discovery task.

Finally, in chapter 8, we discuss the interpretability of our models. We show that shallow models are naturally interpretable. Deep models, however, are more challenging in terms of interpretability. We discuss how to use decision trees to retrace the decision leading to a fault cluster, making it understandable for a domain expert.

The previous chapters discussed constitute the body of work done during the thesis. In chapter 2, we additionally introduce some background notions necessary for reading this manuscript. Following the body of the manuscript, we give a conclusion and directions of interesting future research in chapter 9.

In the following figure we give a guide for reading this thesis, specifically the dependencies between the chapters. The background chapter is necessary for the comprehension of practically all chapters. The variational inference appendix is used in all chapters except chapter 6. Chapter 6 is a theoretical contribution of the thesis, that stands on its own, and is used in chapter 7. In this chapter we present background subjects relevant to the overall manuscript, as well as the mathematical notations used throughout the manuscript. In the case where some additional materiel is relevant, it will be introduced in the relevant chapter. The current chapter is organized as follows: in section 2.1 we introduce important notations. In section 2.2 we present Probabilistic Graphical Models on which the whole manuscript is based. In sections 2.3 and 2.4, we present Markov Chain Monte Carlo and Variational Inference methods, required for inference and learning in probabilistic graphical models. In section 2.5, we give an introduction of deep learning methods required for the second part of the manuscript.

.

The following table gathers a non exhaustive list of mathematical notations and conventions which will be used throughout the manuscript chapters:

Notation Definition

x We denote by a lowercase letter a scalar or a scalar valued random variable.

x

We denote by a bold lowercase letter a vector or a vector valued random variable.

X

We denote by an uppercase bold letter a tensor or a tensor valued random variable.

x i denotes the i th scalar of vector x.

x -i denotes the vector x without the i th entry x i .

f θ (•)
N (x; µ, Σ) = 1 (2π) d det(Σ) e -1 2 (x-µ) T Σ -1 (x-µ) ,
where d represents the dimension of x and det(•) denotes the determinant operator.

Beta(•; γ 1 , γ 2 ) denotes the beta distribution of parameters γ 1 and γ 2 , defined by:

Beta(β; γ 1 , γ 2 ) = Γ(γ 1 + γ 2 ) Γ(γ 1 )Γ(γ 2 ) β γ 1 -1 (1 -β) γ 2 -1 ,
where Γ represents the gamma function.

Cat(•; π) denotes the categorical distribution of parameters π, defined by:

Cat(z; π) = K ∏ k=1 π 1[z=k] k s.t. K ∑ k=1 π k = 1,
where z takes values in the discrete set {1, ..., K}.

Dir(•; α, K) represents the Dirichlet distribution of parameters α and K, defined by:

Dir(x; α, K) = 1 B(α) K ∏ k=1 x α k -1 k s.t. K ∑ k=1 x k = 1
where B(•) represents the multivariate beta function.

∇ x f

represents the gradient of a scalar function f w.r.t. to the variable x, defined by:

[∇ x f ] i = ∂ f ∂x i .
Notation Definition

E x∼p

Expectation taken w.r.t. the random variable x following the distribution p(•).

H[p] or H[p(x)]

represents the entropy of probability distribution p(•) defined by:

H[p] = -p(x) log p(x)dx D KL [q||p] or D KL [q(x)||p(x)
] represents the Kullback-Leibler divergence between two probability distributions q(•) and p(•) defined by:

D KL [q||p] = q(x) log q(x) p(x) dx
.

The previous notations are not exhaustive, if necessary additional notations will be introduced in the relevant chapters.

. Probabilistic Graphical Models (PGMs) are a class of models aimed at specifying the dependency structure over random variables of the problem treated. The advantages of such representations are numerous, the most important of which are:

• The ability to visualize the implications and properties of the model (dependency between interacting random variables).

• PGMs motivate the design and discovery of new models, generalizing previous models.

• They provide an approach to manipulate complex probabilistic computations through graphs.

PGMs are a central topic throughout the thesis. For a brief introduction, consider two vector random variables h and x, where the h random variable is hidden or unobserved (i.e. no samples of data exist for such variables), and the x random variable is observed, thus samples of this random variable are available. Let us further assume that the random variable x depends on the random variable h. This provides us with a probabilistic graphical model specified by the previous assumptions. One can visualize this model by sketching the dependencies as presented in figure 2. Observable variables are denoted by gray circles, hidden variables are denoted by clear circles. Arrows or directed edges between the circles represent causality or conditional dependency, for example, in our case x depends on h.

Probabilistic graphical models provide us with a clear modeling paradigm for interacting parts in an environment, thus making them highly interpretable and subject to inspection. In addition, PGMs can be seen as generative models, where the observable quantities are generated based on hidden or unobserved factors (hidden variables). This approach is very useful for modeling causality in real-world applications. A classical example is given in figure 3, where an image is considered as an observable random variable. One can represent the process by which an image is generated based on the orientation and position of objects in it, by supposing that the object, position, and orientation are hidden factors on which the image depends.

Maximum Likelihood (ML) based learning

Let us reconsider the previous toy model of figure 2. Creating a model in the sense of PGMs is the task of specifying a joint probability distribution over the random variables p θ (x, h) often parameterized by some learnable parameters θ. Given N samples of the observed random variable x, independent and identically distributed (i.i.d.): x (1:N)

= {x (n) } N
n=1 , the objective is to maximize the likelihood w.r.t. the parameters θ defined by:

l(θ) = N ∑ n=1 log p θ (x (n) ).
(2.1)

Maximizing the function l is done mainly through the computation of the gradient ∇ θ l, either by setting it to zero and deriving a fixed point equation, or through gradient ascent algorithms. However, computing the gradient requires marginalization over the hidden variables h:

p θ (x) = p θ (x, h)dh. (2.2)
Exact marginalization is often intractable, except for simple models. We will revisit different approaches to bypass this issue in later sections.

Inference of hidden variables

In our modeling assumptions, hidden variables are generally what interests us, they often represent an interesting feature of the problem that has to be inferred from the data, for example, a class, or the position of an object as shown in figure 3. The inference process requires the computation of the posterior of the latent or hidden variables knowing the observed random variables. Namely, it requires the computation of the quantity:

p θ (h|x) = p θ (x, h) p θ (x) (Bayes rule) = p θ (x, h) p θ (x, h)dh (2.3)
thus facing the same issue of marginalizing over the hidden variables.

For further details on probabilistic graphical models, we suggest the following books by [START_REF] Bishop | Pattern recognition and machine learning[END_REF], [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF][START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF][START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF].

. As discussed in section 2.2, the main issue with the learning and inference of probabilistic graphical models, is the exact marginalization of equation (2.2). Fortunately, Markov Chain Monte Carlo (MCMC) methods provide us with a class of approaches to bypass the issue of exact marginalization through sampling based approximations [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. The main two problems solved by such approaches are the following:

• Problem 1: generating S samples {x (s) } S s=1 from a given probability distribution p(x) = p(x 1 , ..., x d ).

• Problem 2: The estimation of quantities of the form:

E x∼p [ f (x)] = f (x)p(x)dx. (2.4)
Often what interests us is problem 2, however, by solving problem 1, problem 2 becomes trivial. Indeed, for the time being let us suppose that problem 1 is solved, i.e. we have obtained S samples {x (s) } S s=1 drawn from the probability distribution p(•). An estimator of the expectation is given by: s) ).

Î = 1 S S ∑ s=1 f (x (
(2.5)

Notice that if {x (s) } S s=1 are generated from the probability distribution p(•)

the expectation of Î is equal to E x∼p [ f (x)].
Furthermore, the variance of the estimator is proportional to S -1 , thus the more samples drawn the more accurate the estimator becomes.

The theory of Markov chains provides us with a general formalism in order to evaluate sampling algorithms, in terms of consistency of the samples w.r.t. the desired sampled distribution [START_REF] Neal | Probabilistic inference using Markov chain Monte Carlo methods[END_REF]. Starting from an initial probability distribution denoted by π (0) (x) the Markov chain simulator constructs a probability distribution at time t: π (t) (x). The objective is that at convergence: t → ∞, π (t) (x) converges to the target distribution p(x). The Markov chain simulator is defined by the initial probability distribution π (0) (x), and the transition kernel T(x , x) such that:

π (t+1) (x ) = T(x , x)π (t) (x)dx.
(2.6)

In order to guarantee uniqueness of convergence to the target distribution, some required properties have to be respected. Namely, the invariance property, i.e. the target distribution has to be invariant under the kernel:

p(x ) = T(x , x)p(x)dx (2.7)
and the ergodicity property, i.e. the Markov chain converges to the target distribution irrespective of the choice of the initial distribution π (0) (x).

Throughout the manuscript we will focus on the Gibbs sampling algorithm [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], which is an instance of Markov Chain Monte Carlo algorithms, where the transition kernel is given by:

T(x , x) = d ∏ j=1 p(x j |x -j ).
(2.8)

The key advantage of Gibbs sampling is in the fact that it utilizes the conditional distributions p(x j |x -j ), which are often easy to derive in closed form. Algorithm 1 summarizes the procedure of Gibbs sampling.

A complete review of Markov Chain Monte Carlo methods is beyond the scope of the thesis, however for more details on MCMC methods we highly recommend the book by MacKay, 2003 and the seminal papers by [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] and [START_REF] Neal | Probabilistic inference using Markov chain Monte Carlo methods[END_REF].

. Variational inference is a class of approaches that provide us with an alternative way of computing or approximating the posterior of equation 2.3, without resorting to exact marginalization. Variational methods have their origin in the early work of Euler on the calculus of variations, where the (∞) idea is to minimize a functional (a function of a function) with respect to its argument, which is also a function. The calculus of variations is closely related to standard calculus, often we only require the computation of a functional derivative, i.e. elementary changes of the functional subject to small changes in the argument function [START_REF] Feynman | [END_REF].

Algorithm 1 Gibbs sampling Initialize x (0) for t = 1...∞ do sample x (t) 1 ∼ p(x 1 |x (t-1) -1 ) sample x (t) j ∼ p(x j |x (t) 1 , ..., x (t) j-1 , x (t-1) j+1 , ..., x (t-1) d ) sample x (t) d ∼ p(x d |x (t) -d ) end for return x
In the case of probabilistic graphical models, we transform the inference problem into an optimization problem of a certain functional. The functional is derived from the likelihood of equation (2.1) and often referred to as the evidence lower bound or (ELBO), a term that we will use throughout the manuscript. Let us reconsider our toy probabilistic graphical model of figure 2, defined by the joint distribution p θ (x, h) over the x and h random variables. The evidence lower bound is derived through Jensen's inequality as follows:

l(θ) = log p θ (x (1:N) )
= log p θ (x (1:N) , h)dh = log q(h)p θ (x (1:N) , h) q(h) dh ≥ q(h) log p θ (x (1:N) , h) q(h) dh := L(θ, q).

(2.9)

The main idea of variational inference is to introduce an auxiliary probability distribution q(h) over the hidden random variables, called the variational distribution. Throughout the thesis we will reserve the letter q to denote variational distributions as opposed to true model distributions denoted by p. The variational distribution serves as an approximate or surrogate posterior over the hidden random variables, and the main objective is to find the probability distribution q that maximizes the evidence lower bound L at θ fixed. In order to show the relation between q(h) and the true posterior p θ (h|x (1:N) ), we can rewrite the evidence lower bound as: L(θ, q) = q(h) log p θ (x (1:N) , h) q(h) dh = q(h) log p θ (h|x (1:N) )p θ (x (1:N) ) q(h) dh = q(h) log p θ (h|x (1:N) ) q(h) dh + log p θ (x (1:N) ) q(h)dh

= l(θ) -D KL [q(•)||p θ (•|x (1:N) )].
(2.10)

Where D KL represents the Kullback-Leibler divergence, a measure of distance between two probability distributions. Notice that the evidence lower bound is a functional of the probability distribution q(•), hence the variational aspect of the problem.

The objective of the variational distribution is to approximate the true posterior. This can be seen through the Kullback-Leibler term, by maximizing the evidence lower bound w.r.t. q, we minimize the Kullback-Leibler term, thus minimizing the distance between the approximate posterior q and the true posterior p θ (•|x (1:N) ). The optimization problem can be written as:

max q L(θ, q) s.t. q(h)dh = 1.
(2.11)

The optimization problem of equation (2.11) is often intractable over the entire space of possible probability distributions. Hence, an often used strategy is to restrict the space of search. One of the most used restrictions is what we call the mean-field hypothesis [START_REF] Jaakkola | tutorial on variational approximation methods[END_REF][START_REF] Jordan | An introduction to variational methods for graphical models[END_REF], in which we suppose that the variational distribution respects a factorization condition: q(h) = ∏ j q(h j ).

(2.12)

Another recent popular approach is to choose a specific distribution and parameterize it through some set of parameters φ: q φ . In this case the optimization is transferred to the parameters φ, thus returning to classical optimization methods over the function L(θ, φ) with respect to its parameters [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF].

The final piece of the puzzle is to ensure that the optimization of L w.r.t. the parameters θ follows the maximum likelihood criterion. This is the case here, and it can easily be proven under the hypothesis that the optimal solution of (2.11) q * approximates the true posterior. Indeed, if: (2.13) and thus with equation (2.10), we can deduce that: L(θ, q * ) ≈ l(θ).

q * (•) ≈ p θ (•|x (1:N) ) then, D KL [q * (•)||p θ (•|x (1:N) )] ≈ 0,
(2.14)

Hence maximizing L(θ, q * ) w.r.t. θ is equivalent to maximum likelihood estimation.

For further details on Variational Inference methods for probabilistic graphical models we highly recommend the book by [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF].

. Deep learning as a field has retained considerable interest from researchers from all disciplines in recent years. The ability of deep neural networks to solve with high accuracy machine learning tasks ranging from classification and regression to reinforcement learning, is one of the most noteworthy breakthroughs in recent history. Our interest throughout the thesis in neural networks will be mainly in their ability to learn complex features relevant to the task at hand. For this purpose, we will present neural networks as generalized function approximators, defined by the application of successive linear and nonlinear functions [START_REF] Rosenblatt | Principles of neurodynamics. perceptrons and the theory of brain mechanisms[END_REF]. A neural network in our case will be defined by a function often denoted by g, such as: 1) .

g θ = f (L) • σ • f (L-1) • ... • σ • f (2) • σ • f (
(2.15)

where each function f l represents the l th application of the linear or fully connected layer, defined by:

f (l) (h) = h T W l + b l (2.16)
W l , and b l represent the weight matrix and bias vector for the l th layer. • represents the composition of functions operator. Between each application of two linear layers, a non linear function is applied denoted by σ. Throughout the manuscript this function will denote the Rectified Linear Unit [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] defined by:

σ(h) = max(0, h) (2.17)
We will be restricted to the previously defined structure of neural networks as it is sufficient for our applications. The overall parameters of the neural network are assembled in the vector denoted by θ. The parameters are optimized to minimize some cost function relevant to some task, we will often denote this function by C . The optimization is done using a stochastic gradient descent algorithm with momentum [START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF] or a close variant such as the ADAM optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF].

The computation of the gradient of the loss function ∇ θ C is done by the classical backpropagation algorithm [START_REF] Werbos | Applications of advances in nonlinear sensitivity analysis[END_REF]. The algorithm is based on the chain rule of calculus, C can be seen as the composition of two functions u and v, that is C = u(y) and y = v(x), then:

∂C ∂x i = ∑ j ∂C ∂y j ∂y j ∂x i (2.18)
For further details on Deep learning methods and neural network architectures, we highly recommend the following book by [START_REF] Goodfellow | Deep learning[END_REF]. Identifying fault patterns in large telecommunication network and service infrastructures is a difficult task. Most of the considered technical solutions rely either on rule based expert systems or hand crafted expert Bayesian networks [START_REF] Kogeda | A probabilistic approach to faults prediction in cellular networks[END_REF]Kogeda and Agbinya, 2006;[START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. Although these approaches have had tremendous success, one of the unsung drawbacks is the data processing and the expert knowledge required to build the diagnosis model or rules. For expert systems, for example, the rules created by the expert require knowledge of the existing fault, and the identification of the variables describing the fault. This process requires data processing by hand by an expert of the domain, which is an expensive and time consuming task. Also, the maintenance of the model or rules in the long run can be a significant issue for operational teams.

Part I S H A L L O W P R O B A B I L I S T I C G R A P H I C

Recently, machine learning techniques have been tremendously successful in the identification and the extraction of patterns in various domains. In the context of our problem, similar approaches can be used to identify patterns of faults from diagnosis data. This task is thus an unsupervised machine learning task, where labels are not available and obtaining them is as complex as the data processing required to construct expert rules. However, the data gathered from various devices and services in the network is often structured in the form of a table, where each variable takes some range of values.

Clustering such data, gathered from telecommunication networks and services presents many challenges. The first and the main challenge is the unknown number of clusters of faults in the data. The second challenge is the types and multivariate nature of the data. The data is multi-dimensional and can contain categorical and continuous variables. Therefore, classical clustering algorithms where the number of clusters is to be set a priori require some form of model selection. Furthermore, classical approaches such as KMeans suppose a specific probability distribution for each cluster. These modeling assumptions, including the assumption on the number of clusters, can hurt the performance of the clustering when the data do not comply with such assumptions, which is often the case when dealing with real-world applications.

In this chapter, we present an infinite multivariate categorical mixture model to identify patterns of faults in an unsupervised setting, without any prior expert knowledge, and without the requirement to know a priori the number of different fault patterns. The model is based on the Dirichlet Process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], which allows for learning the number of clusters from the data. However, the Dirichlet Process supposes an infinite number of clusters which translates to an intractable inference problem on the model. Our contributions are the following:

• We propose an infinite multivariate categorical mixture model for clustering categorical data (section 3.3).

• We show how to perform approximate inference on the model, in order to extract the clusters from the data using Variational Inference (section 3.4).

• We demonstrate how the model is able to identify root causes of faults in a synthetic dataset generated from a real-world expert Bayesian Network (section 4).

• We also demonstrate the clustering performance of the model on real operational data acquired from the Fixed Access Network and the Local Area Network (section 5). Mixture models are one of the simplest forms of probabilistic graphical models. The objective of mixture models is to model the distribution of an observable random variable x as a mixture of base distributions.

The main application goal is to identify patterns captured by the base distributions.

The assumption is that the observable random variable x depends on a latent variable denoted by z taking values in a discrete set {1, ..., K}, with probabilities π k = p(z = k). Depending on the value taken by z, the random variable x follows some base probability distribution p (z)

θ , where θ and π are parameters of the model to be learned to fit the data through maximum likelihood (2.1). Figure 4 shows a graphical representation of a mixture model.

. . The Dirichlet Process

In a classical mixture model, we must specify a number of components (or clusters) for the random variable z. Identifying the "correct" number of components is a task full of challenges, that has been tackled through model selection approaches [START_REF] Bishop | Pattern recognition and machine learning[END_REF]. The Dirichlet process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] allows us to identify directly the number of components much more naturally. The idea behind the Dirichlet process is to suppose an infinite number of components, and specify a prior on the creation of new components. The creation of new components has a prior probability defined by draws from a beta distribution, formally:

β k ∼ Beta(•; 1, η) π k = β k k-1 ∏ l=1 (1 -β l ) (3.1)
The probability of the k th component being active is governed by draws from the beta distribution of hyperparameter η. One way to think of such construction is through the stick-breaking construction [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF]. The stick breaking construction is as follows, π k represents the length of the k th piece broken from a stick of length 1. If the concentration parameter is small, π k will be large (close to 1) at the first times the stick is broken, therefore the stick will be broken a small number of times (small number of clusters). If η is large, π k will be small (close to 0) and the stick can be broken a large number of times (large number of clusters). The probability that the (n + 1) th data point belongs to a new cluster k * or to K existing clusters is [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF]:

P[z n+1 = z|z 1:n , η] = 1 η + n η1[z = k * ] + K ∑ k=1 n k 1[z = k] (3.2)
where n k is the number of data points in cluster k. If η → ∞ a new cluster would be created for each data point, and if η → 0 all data points are concentrated in the first cluster. The intuition behind the Dirichlet process is the following: as the number of data points increases, we allow the number of clusters to grow according to the concentration parameter and the clusters already assigned. New samples are assigned to existing clusters if they match, otherwise a new cluster is created for them. Therefore the Dirichlet Process allows the mixture model to cluster the data and identify automatically the number of clusters necessary to explain the data.

. . Notations

In our case, we consider that the random variable x represents the variables or features of a network such as an equipment status, an alarm, or a physical metric. We only consider categorical random variables. Continuous variables such as optical powers or temperatures are discretized, using standard methods such as equal frequency or equal width discretization. Thus, each random variable takes values in a setE i = {v 1i , ..., v |x i |i }, where |x i | is the number of modalities of variable x i . Let x (1:N) represent N samples of the vector x = [x 1 , ..., x d ] T of dimension d. We denote by z (1:N) N random variables, where z n represents the fault cluster of sample x n .

. . .

Structure and Definition

Using the stick-breaking construction introduced in the previous section for the weights π k of the k th cluster, the generative process for our proposed Infinite Categorical Mixture Model is:

z ∼ Cat(•|π) = ∞ ∏ k=1 π 1[z=k] k (3.3) B ki ∼ Dir(•; α i , |x i |) ∝ ∏ v∈E i B α iv -1 k,i,v (3.4) s.t ∑ v∈E i B k,i,v = 1 (3.5) x i |z = k, B ∼ Cat(•|B ki ) = ∏ v∈E i B 1[x i =v] k,i,v (3.6)
where π k is defined through the Dirichlet process defined in the previous section.

Figure 5 shows a graphical representation of the generative process of the model. The observed variable x depends on the cluster assignment z, and the parameters B. We adopt a full bayesian approach where we consider the parameters of the model as random variables themselves. Thus, the parameters B have a prior that we choose as a Dirichlet prior with hyperparameters α. Based on the assignment of cluster z a sample x i for dimension i is drawn based on the conditional probability for variable x i taking a certain value v ∈ E i , defined by the B tensor:

B k,i,v = p[x i = v|z = k] (3.7)
The Dirichlet prior on the parameters B is governed by the hyperparameter α called the concentration parameter. It allows us to inject prior knowledge about the modalities of variable x i . α iv thus represents a weight for modality v in the Dirichlet distribution associated with variable x i . In the case where no information is available we use an uninformative prior α iv = 1 |x i | .

. .

Learning and Inference

In order to fit the model to the data and identify the clusters representing the different types of network failures, we need to compute or approximate the posterior distribution:

p(z (1:N) , B, β|x (1:N) ) = p(z (1:N) , B, β, x (1:N) ) p(x (1:N) ) (3.8)
Markov Chain Monte Carlo (MCMC) approaches are commonly used to fit such models [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. The main idea behind these approaches is to run a Markov chain long enough (until convergence) where the stationary distribution at convergence is the posterior of interest as defined in equation (3.8).

The advantage of MCMC methods is that they are guaranteed to converge to the true posterior. However, one major drawback of such approaches arises with large datasets, in this case, MCMC methods are very expensive computationally [START_REF] Blei | Variational inference: A review for statisticians[END_REF]. Variational inference by comparison can utilize all the optimization tools to bypass this issue, such as Stochastic optimization and distributed optimization [START_REF] Robbins | A stochastic approximation method[END_REF], making the inference run faster. Furthermore, in terms of predictive accuracy it has been shown that variational inference and MCMC methods often have the same performance [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF][START_REF] Braun | Variational inference for largescale models of discrete choice[END_REF].

In the next section we detail the variational inference approach for our proposed model.

. . . The mean field approximation

In order to introduce the variational inference approach, we denote by ζ = ζ 1:m = {z (1:N) , B, β} the vector grouping all the hidden variables of the model. Solving the inference problem of the model amounts to determining p(ζ|x (1:N) ). As detailed in section 2.4, closed form solutions for this quantity are intractable. Variational inference and the mean-field approximation allow us to approximate this intractable distribution by a set of distributions for which the inference is tractable, namely the mean field family. A distribution q is said to be in the mean field family for variables ζ, if it verifies:

q(ζ) = m ∏ l=1 q(ζ l ) (3.9)
The main idea of variational inference is to approximate the intractable distribution (3.8), by finding the closest mean field family member in terms of Kullback-Leibler divergence i.e.:

q * = min q D KL q(•)||p(•|x (1:N) ) (3.10)
By exploiting the factorization in the mean field family, we can show that the solution q * verifies the following fixed point equations:

log q * (ζ l ) = const + E ζ -l ∼q * log p(ζ, x (1:N) ) ∀l (3.11)
For an explicit derivation of this criterion we refer the reader to appendix A.

In the case of our model, the mean-field family is defined as :

q(z (1:N) , B, β) = N ∏ n=1 q(z n ) d ∏ i=1 T ∏ k=1 q(B k,i ) T ∏ k=1 q(β k ) (3.12)
We also suppose that q(β T = 1) = 1 hence q(z n > T) = 0, i.e. the number of clusters is truncated to an upper bound on the true number of clusters [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF]. A noteworthy aspect of this approach is that the true posterior given by (3.8) has an infinite number of factors, however the approximate distribution q is constrained based on the previous conditions. Therefore, the true model is unchanged however the minimization problem is relaxed in order to be solved efficiently.

. . Optimal Variational Distributions and Update Equations

By applying equation (3.11) to our Infinite Categorical Mixture Model, we obtain:

log q * (z (n) ) = const + E {z -n ,β,B}∼q * log p(z (1:N) , B, β, x (1:N) )
Algorithm 2 Variational Inference for the ICMM

Input: x (1:N) , T, η Initialize φ, , γ 1 , γ 2 L (0) = -∞ for t = 1...∞ do Compute: γ (t)
φ nk ∀n return z (1:N) , φ log q * (B k,i ) = const + E {z (1:N) ,β,B -{k,i} }∼q * log p(z (1:N) , B, β, x (1:N) ) log q * (β k ) = const + E {z (1:N) ,β -k ,B}∼q * log p(z (1:N) , B, β, x (1:N) ) (3.13)
And by substituting the expression of p(z (1:N) , B, β, x (1:N) ) resulting from the graphical representation of the model (Figure 5), we then deduce the following approximating distributions:

q * (z (n) ) = Cat(z n ; φ n ) q * (B k,i ) = Dir(B k,i ; k,i , |x i |) q * (β k ) = Beta(β k ; γ 1,k , γ 2,k ) (3.14)
And the mean field fixed point equations for the parameters are the following, where ψ is the digamma function (proofs given in appendix B.1):

log φ nk = d ∑ i=1 ∑ v∈E i 1[x (n) i = v][ψ( k,i,v ) -ψ( ∑ v ∈E i k,i,v )] + ψ(γ 1,k ) -ψ(γ 1,k + γ 2,k ) (3.15) + k-1 ∑ l=1 [ψ(γ 2,l ) -ψ(γ 1,l + γ 2,l )] s.t T ∑ k=1 φ nk = 1 γ 1,k = 1 + N ∑ n=1 φ nk (3.16) γ 2,k = η + N ∑ n=1 T ∑ l=k+1 φ nl (3.17) k,i,v = α i,v + N ∑ n=1 φ nk 1[x (n) i = v] (3.18)
Thus, the identification of the optimal approximating distribution is done by iterating the previous equations. The stopping criterion is defined using the Kullback-Leibler divergence:

L(q) = -D KL q(•)||p(•, x (1:N) ) (3.19)
At convergence the value of L should reach a plateau corresponding to a lower bound on the likelihood. Algorithm 2 summarizes the procedure for learning the approximating posterior and for inferring the clusters.

.

We test our proposed model on the task of fault clustering in two datasets.

The first one is a synthetic dataset generated using an expert bayesian network, and the second one a real dataset collected from equipment deployed in a real fiber-to-the-home network. We evaluate our model in comparison with other clustering baselines. Furthermore, we also compare the Variational Inference approach with the Gibbs sampling approach.

. . Datasets

Synthetic dataset

In order to test our model in a controlled environment, we reuse the work done by [START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. The authors designed an expert bayesian network that reflects the behaviour of a real GPON-FTTH network. The bayesian network is depicted on figure 6 where the colored top nodes represent the root causes and the other nodes the observations coming from the network. Note here that the bayesian network is only used to generate a synthetic dataset similar to operational data, it is not used in the clustering method.

In order to generate the synthetic data for the experiment, we simulate a fault by activating the root cause for the fault, and we sample the visible variables. Here we sampled 6 uncorrelated faults (for more details on the expert bayesian network and remaining variables, we refer the reader to [START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]): • AltONT: highlights a problem with the power supply of the optical network termination (ONT). This hidden variable controls the current of the ONT I ONT , the Dying Gasp alarm DG, and the electric voltage of the ONT V ONT .

• AltOLT: describes a problem with the power supply of the optical line termination (OLT). It controls similar variables for the OLT, V OLT and I OLT .

• FaultyONT: denotes whether the ONT is faulty, and symbolizes the global state for an ONT. It controls the alarms TIA (Transmission Interference Alarm), DOW (Drift Of Windows), and SF (Signal Fail) or SD (Signal Degraded).

• FiberDB: represents the state of the drop optical fiber and controls Rx ONT the power at which the ONT receives the signal.

• IOS: (Image Operating System) refers to an incompatibility between the OS of an OLT and the OS of an ONT. It controls the variable SWV describing the software version alarm.

• TcOLT: represents the temperature of the OLT.

We sampled 150 data points for each fault, using the likelihood weighted sampling method, based on the conditional distribution tables of each observation given the root cause [START_REF] Fung | Weighing and integrating evidence for stochastic simulation in Bayesian networks[END_REF]. Therefore we generate for each fault a pattern of the visible variables for a specific customer equipment (ONT). The resulting dataset contains 900 samples, where each sample represents a realization of the 29 visible nodes of the Bayesian network.

Real-world dataset

The second dataset used for evaluation is collected from a real world operating GPON-FTTH and Local area network depicted in figure 7. For the fixed access network we collected data describing the status of the OLT, and the status of the ONT (olt_status, ont_status, ont_download_status), in addition to continuous metrics such as the powers of transmission and reception in each equipment. For the local area network, we collected data describing the status of the router and the different services provided such as, IPTV, VOD (Video On Demand) and VoIP. We also collected data describing the account status of the client.

The data collected corresponds to instances where some faults are known (identified) and others where the fault is not known to experts or expert systems. The unidentified fault could be the same as a known fault or entirely novel to the system. In this chapter, we focus on the unidentified faults representing 70.6 % of the dataset. We adopt a full unsupervised method, i.e. we consider that no ground truth is available and we attempt to identify patterns blindly.

Table 1 summarizes each dataset and its characteristics.

. . Baselines and Metrics

We compare our clustering approach to other state of the art clustering methods, where the number of clusters is not required and is determined automatically. The baselines chosen are the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF], the affinity propagation clustering method [START_REF] Frey | Clustering by passing messages between data points[END_REF] method [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF], and the infinite Gaussian mixture model (IGMM) [START_REF] Rasmussen | The infinite Gaussian mixture model[END_REF].

The evaluation metric is the clustering accuracy, it is similar to the classification accuracy, however, in the clustering task the clusters are not identifiable with the labels, i.e. the clusters can change from one run of the algorithm to another. Therefore we need to test all possible combinations and choose the best one. The clustering accuracy is defined as in [START_REF] Xie | Unsupervised deep embedding for clustering analysis[END_REF]:

ACC = max m∈M ∑ N n=1 1 [l n = m(z n )]
N where z n is the cluster assignment, l n the true label and M the set of all possible one-to-one mappings.

. . . Clustering results on the synthetic dataset

Here, we report the clustering results of our model on the synthetic dataset. First, we compare our model with the considered baselines. In table 3, we report the mean and standard deviation of the clustering accuracy for each model over 10 runs. A search over the hyperparameters of each model has been performed and the best performing model is selected for each baseline. As shown our model outperforms all models considerably. This

Methods

Affinity propagation Mean shift DBSCAN IGMM Ours

Accuracy 0.35 ± 0.0 0.69 ±0.0 0.29 ± 0.0 0.57 ± 0.09 0.9 ± 0.07

shows the advantage of explicitly modeling a clustering problem to match the specificity of the data. Using classical approaches blindly often leads to poor results. Second, we report the evidence lower bound of equation (3.19) across variational inference iterations, for 10 runs of the model with their standard deviation. The evidence lower bound gives us a criterion for evaluating the convergence of the model. As depicted in figure 8, the evidence lower bound increases until reaching a plateau, which indicates that the model has converged.

Finally, in table 2, we report the confusion matrix for the best model (highest elbo). As shown, the confusion matrix is diagonal with the majority of instances well clustered to their ground-truth classes of faults. Furthermore, and more importantly, the advantage of the Dirichlet process is clearly demonstrated. We start with 50 clusters, however only 6 are filled corresponding to the ground-truth classes, the others are kept empty. Therefore, the model clusters the data and automatically identifies the correct number of clusters.

. . Variational inference vs Gibbs sampling

We report the results of the comparison between the Gibbs sampling fitting of the model and the variational inference method. We compare the two approaches in terms of complexity and accuracy on the synthetic dataset.

As seen in figures 9 and 10, both methods arrive at high accuracy, they have the same performance in terms of predictive accuracy around 90 %. However, as the number of instances in the dataset increases, we see that the Gibbs sampling convergence time increases considerably. In contrast with the variational inference approach the convergence time stays stable. This is a foreseeable consequence, that stems from the fact that iterations of Gibbs sampling can only be done sequentially and are of complexity O(NKD), where N is the number of instances, K the number of classes, and D the dimension of the observable random variable x. Variational inference on the other hand, has complexity O(NK) (updates of φ, others are negligible). Furthermore, the updates of variational inference can be written in matrix form that can be run in parallel form.

. . Pattern recognition on real-world data

In this second experiment, we apply our clustering model on real operational network data. However in this setup we don't have access to a ground truth so the aim of this experiment is to demonstrate the benefit of our clustering model in a data exploration process. We obtained 6 main clusters corresponding to known problems. We report for each cluster the number of customers assigned to the cluster and the most common pattern, i.e. characteristic values that the variables take and the counts for each variable (Figure 11). Furthermore, we give an interpretation based on network expertise to understand what kind of problem the customers are facing in each cluster.

• Cluster 37: (600 clients) This cluster gathers customers with a problem on the remote PVR for the IPTV: Remote_PVR = Missing, and all other values correspond to an operational state of the client's line: (client_account_status=1, OLT_status=OK, router_status=Enabled, ONT_status=NODEFECT...).

• Cluster 30: (> 2500 clients) In this cluster the ONT is not detected: (ONT_status = Missing), where all other services are operational (router_status=Enabled, voip_status=Ok, tv_profile_status=operational ... ).

• Cluster 43: (≈ 2500 clients) This cluster corresponds to the optimal behavior of the network and services, all equipment are detected and all services are operational.

• cluster 44: (≈ 550 clients) This cluster describes a problem with the router where the router is disabled (router_status=DISABLED, router_ipv6_status=DISABLED), however all other equipments are operational and client account status is activated.

• Cluster 40: (≈ 27 clients) This cluster represents the case where an account is suspended and the VOIP service is down (client_account_status = 0, tv_profile_status=SUSPENDED, VOIP_status=KO), however all pieces of equipment are operational.

• Cluster 12: (≈ 130 clients) This cluster gathers clients with deactivated WiFi (router_WiFi_status = Down, router_status=Enabled), and all other equipment is operational (ont_status = NODEFECT, OLT_status=OK...).

This experiment shows that our clustering method can be helpful in a data exploration process where the number of clusters is not known a priori. With a simple network expertise we can see that the clusters are relevant and describe a particular behaviour on the network.

. &

Machine learning techniques for fault detection in telecommunication networks are promising. Bayesian networks have been dominant in this domain [START_REF] Kogeda | A probabilistic approach to faults prediction in cellular networks[END_REF]Kogeda and Agbinya, 2006;[START_REF] Tembo | Model-Based Probabilistic Reasoning for Self-Diagnosis of Telecommunication Networks: Application to a GPON-FTTH Access Network[END_REF]. The main advantage of Bayesian network modeling is their easy interpretable nature. Root causes of faults are modeled by variables of the Bayesian network and inference can be performed in order to determine the root cause from the observable variables. This approach however, requires expert knowledge of the specific domain of faults and a time consuming task to build the Bayesian network and all the dependencies between the variables. Recently, researchers have been interested in learning the structure of the Bayesian network from data [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF], and such approaches have been successfully used for the self-diagnosis problem [START_REF] Bengio | Self-diagnosis technique for virtual private networks combining Bayesian networks and case-based reasoning[END_REF]. However, the main problem, in such approaches, is the loss of the interpretability. The structure learned from the data can be one of a class of equivalence of optimal structures, hence, this often results in structures that are optimal but hard to interpret.

The second group of machine learning approaches, are anomaly detection based approaches [START_REF] Chandola | Anomaly Detection: A Survey[END_REF]Hood and Ji, 1997). Although these methods allow for accurate detection of anomalies Figure 11: Assignments defining each cluster from most occurring to least occurring from data, their main drawback is that all anomalies identified are grouped in a single cluster. Classifying each anomaly requires either a relabeling process by hand or a clustering process. Therefore, the natural reformulation of the self-diagnosis problem is as a clustering problem. Clustering types of faults from the data has been previously proposed [START_REF] Adda | Comparative analysis of clustering techniques in network traffic faults classification[END_REF][START_REF] Hashmi | Enabling proactive self-healing by data mining network failure logs[END_REF]. However, most of these approaches rely on classical algorithms where the number of clusters is known a priori or estimated by model selection. On the contrary, our approach allows for automatic determination of the number of clusters and requires no intervention from an expert. A post processing of the clusters by an expert is still needed in order to identify the fault present in each cluster and the related root cause. [START_REF] Chen | Failure diagnosis using decision trees[END_REF] proposed a method based on decision trees to accomplish this task.

. Clustering patterns of faults in network data can be achieved using the categorical mixture model approach as proposed in this chapter. However, two main limitations arise and are quite clear at this stage. The first, is the fact that the model can be applied only to categorical random variables. Thus, continuous variables have to be preprocessed in advance. The second limitation stems from the fact that labels are often available as discussed in subsection 3.5.1. However, the proposed infinite mixture model as formulated can only cluster the data in an unsupervised manner, without taking advantage of the existing labels available. In the next chapter, we extend the definition of the mixture models to treat categorical and continuous random variables alike. Furthermore, we show how we can use semi-supervision using the labels available to guide the clustering based on the ground-truth available.

. We have introduced the infinite categorical mixture model, for clustering faults in categorical network data. We showed how the model can be used to analyze data gathered from telecommunication networks and services and how to discover fault patterns in an unsupervised manner. The model is capable of identifying the correct number of clusters to analyze the data using the Dirichlet process. We showed how Variational Inference can be used to perform inference on the model. This approach allows inference to scale well with the dimensionality of the data, and the convergence of the model can be determined explicitly.

S E M I -S U P E R V I S I O N A N D M I X E D -T Y P E I N F I N I T E M I X T U R E M O D E L S .

In the previous chapter, we have explored the problem of clustering faults from network data in an unsupervised manner. As discussed, the previous model, i.e. the infinite categorical mixture model, was capable of identifying patterns of faults and simultaneously inferring the number of faults from the data. However, as its name suggests the model only deals with categorical random variables, meaning that continuous random variables need to be preprocessed in order for the model to be applied. It turns out our model can be generalized to continuous and categorical random variables by expanding the definition of the emission distributions of the visible random variables.

Furthermore, in most real-world applications of diagnostics some expert knowledge is available. The expert knowledge is often represented by partial labels, i.e. an expert looks at a data point and determines which class of fault it belongs too. The infinite categorical mixture model of the previous chapter does not use this extra-information of partial labels. Fortunately, semi-supervised learning methods could be applied to take advantage of both large amounts of unlabeled data and a small fraction of labeled data [START_REF] Chapelle | Semisupervised learning[END_REF].

With slight modifications to our previous infinite categorical mixture model, we could bypass these issues. Indeed, probabilistic graphical models can incorporate labels as observable random variables. This approach has been successfully applied to classical mixture models with the Expectation-Maximization algorithm [START_REF] Nigam | Text classification from labeled and unlabeled documents using EM[END_REF][START_REF] Seeger | A taxonomy for semi-supervised learning methods[END_REF]. Furthermore, continuous data can be incorporated in the probabilistic graphical model through continuous distributions such as the Gaussian distribution. The challenge that arises is the non-conjugacy between the categorical distribu-tions and the continuous ones, making classical inference methods such as Gibbs sampling non viable.

In this chapter, we generalize the previous infinite categorical mixture model to a semi-supervised mixed-type infinite mixture model capable of dealing with categorical and continuous distributions alike, and also capable of incorporating expert knowledge in the form of partial labels. Our contributions are the following:

• We start by giving a formal definition of the semi-supervised mixedtype infinite mixture model as a probabilistic graphical model.

• We show how to solve the inference problem through mean-field variational inference similarly to the previous chapter, and how to use the model to make predictions on future data samples.

• We apply our model on a real-world semi-supervised fault cluster identification problem from GPON-FTTH data, and we compare our model to classical baselines.

.

We start by defining additional notations in this section in order to rigorously introduce the semi-supervision and the continuous random variables. We keep the same definition for the visible random variables as in the previous chapter, i.e. x = [x 1 , ...., x d ]. However, within the d random variables considered, there exists two subsets, those of the continuous and categorical random variables. We denote x C = {x j , if j is a continuous dimension} the vector of continuous random variables, and we denote the vector of categorical random variables by x D = {x j , if j is a categorical dimension}.

Second, in order to introduce the semi-supervision, we consider the set of indexes of the whole dataset of N instances denoted by {1, ..., N} = D l ∪ D u where D u = {n ∈ {1, ..., N}, s.t. x (n) is unlabeled} denotes the indexes of the observed instances that are unlabeled, and, D l = {n ∈ {1, ..., N}, s.t x (n) is labeled} denotes the indexes of the observed instances that are labeled. We denote the n th label for instance x (n) by y (n) .

Finally, we introduce the Wishart distribution, that will serve as a prior in the continuous case over precision matrices. The Wishart distribution is a distribution defined on the space of p × p symmetric, non-negative, definite matrices Λ. Its density function is defined as:

W (Λ; L, ν) = det(Λ) ν-p-1 2 exp{-1 2 tr(L -1 Λ)} 2 νp 2 det(L) ν 2 Γ p ( ν 2 ) , (4.1)
where L is called the scale matrix, ν is a scalar called the degrees of freedom, tr() represents the trace function and Γ p represents the multivariate gamma function. The Wishart distribution serves as a conjugate prior for the precision matrix of a multivariate normal distribution.

.

. .

Structure and definition

The proposed model in this chapter is an extension of the infinite categorical mixture model to the continuous case, with the advantage of introducing semi-supervision. Certain aspects of the previous model remain unchanged, namely, the Dirichlet process and the cluster hidden variables, defined by the following generative process, which is the same as the previous chapter:

β k ∼ Beta(•; 1, η) ∀k ∈ N * (4.2) π k = β k k-1 ∏ l=1 (1 -β l ) (4.3) z ∼ Cat(•|π). (4.4)
Following the generation of the cluster hidden random variable z, we need to specify the emission distributions of the hidden random variables x = [x D , x C ]. We can follow the same procedure as the last chapter to model the categorical random variables, that is:

∀i ∈ D B ki ∼ Dir(•; α i , |x i |) (4.5) x i |z = k, B ∼ Cat(•|B k,i ), (4.6)
where the parameters B and the sets E i are defined in the same way as in subsection 3.3.1.

The last part of the puzzle is to specify emission distributions over the continuous observable random variables given the cluster assignment. In this model, we opt for the simplest of modeling and we will explore its limitations in later chapters. We will suppose that the continuous random vector is Gaussian with conjugate priors on the mean and precision matrix (inverse of the covariance matrix):

x C |z = k, M, Λ ∼ N (•; M k , Λ -1 k ), (4.7)
where the M = [M k,j ] and Λ = [Λ k,i,j ] are respectively the mean and precision matrices for each cluster. In this approach, we also adopt a full bayesian setting where the mean and precision matrices are considered random variables with Gaussian and Wishart conjugate priors respectively:

∀k M k |Λ k ∼ N (•; µ 0 , (κ 0 Λ k ) -1 ) (4.8) Λ k ∼ W (•; L 0 , ν 0 ) (4.9)
The previous bayesian priors are standard in classical Gaussian mixture models [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF]. {µ 0 , κ 0 , L 0 , ν 0 } are hyperparameters to be set empirically.

In Figure 12 we report a graphical representation of the mixed-type infinite mixture model and the dependencies between the random variables.

. . Semi-supervision as partially observable cluster variables

One way to incorporate the partial labels in a probabilistic graphical model is to consider the values as observations of the corresponding hidden variables, i.e. the cluster hidden variable. This approach was introduced by [START_REF] Seeger | A taxonomy for semi-supervised learning methods[END_REF]. The approach consists in supposing that if x (n) is labeled, in its generative process its cluster variable z (n) has to correspond to its label y (n) . This approach generalizes quite nicely to our probabilistic model.

In our case, the condition would be:

p(x (n) , z (n) = y (n) ) = 0 ∀n ∈ D l . (4.10)
The subtlety of this condition remains in the inference process. The previous criterion of equation (3.19) no longer holds as a variational lower bound in this case. In the next section, we show how to generalize it to the semisupervised case, and how to derive the corresponding variational update equations.

. -. . The semi-supervised variational lower bound

The semi-supervision condition of equation (4.10) requires us to revisit the likelihood of our model. Indeed, with this condition the likelihood can be decomposed into two terms:

log p(x (1:N) ) = ∑ n∈D l log p(x (n) ) + ∑ n∈D u log p(x (n) ) = ∑ n∈D l log p(x (n) , z (n) = y (n) ) + ∑ n∈D u log p(x (n) )
The last equality follows from the fact that ∀n ∈

D l : p(x (n) , z (n) = y (n) ) = 0.
Given the decomposition of the log likelihood, the corresponding variational lower bound also involves two terms for the labeled and unlabeled terms.

Let us denote by ζ = {β, z (1:N) , M, Λ, B} the set of all hidden random variables, and ζ -z (1:N) the set without the z random variables. Similarly to chapter 2, using the Jensen inequality as in equation (2.9), and introducing the variational distribution q(•) over the hidden variables, we have:

log p(x (1:N) ) ≥ ∑ n∈D l q(ζ -z (n) ) log p(x (n) , ζ -z (n) , z (n) = y (n) ) q(ζ -z (n) ) dζ -z (n) + ∑ n∈D u q(ζ) log p(x (n) , ζ) q(ζ) dζ = L ss (q) (4.11)
The new semi-supervised evidence lower bound L ss can be written in the form of a Kullback-Leibler divergence, by simply noticing that the semisupervision condition can be transferred to the variational distribution, that is: ∀n ∈ D l : q(z (n) = y (n) ) = 1. We force the labeled instances to also have a variational posterior that respects the labelisation. The semi-supervised evidence lower bound in this case can be re-written as:

L ss (q) = -D KL q(•)||p(•, x (1:N) ) s.t ∀n ∈ D l : q(z (n) = y (n) ) = 1 (4.12)
Since the semi-supervised evidence lower bound can be written as a Kullback-Leibler divergence, we can apply mean-field variational inference in order to find the optimal mean-field distribution maximizing L ss , and thus fitting the model to the observed data. In the next subsection, we apply the mean-field variational inference approach to derive the update equations of the optimal q.

. .

Mean-field variational inference

As discussed in the previous subsection, optimizing the model to fit the data amounts to finding the optimal variational distribution maximizing the semi-supervised evidence lower bound of equation (4.12). i.e.:

q * = arg max q L ss (q) s.t ∀n ∈ D l : q(z (n) = y (n) ) = 1 (4.13)
In order to apply the mean-field approach of Appendix A, we need to specify the variational family. For the discrete dimensions, we keep the same variational factorization as the previous chapter. As for the continuous dimensions, we use a standard variational family for Gaussian mixtures, where the mean and precision random variables are not independent (Murphy, 2012). The variational family of the distribution q has the following form over the set of hidden random variables {β, z (1:N) , M, Λ, B}:

q({β, z (1:N) , M, Λ, B}) = N ∏ n=1 q(z n ) T ∏ k=1 q(β k ) × ∏ i∈D q(B k,i )q(M k |Λ k )q(Λ k ) (4.14)
Given the previous variational family, solving equation (4.13) leads to the following parametric variational distributions:

q * (z n ) = Cat(z n ; φ n ) q * (B k,i ) = Dir(B k,i ; k,i , |x i |) ∀i ∈ D q * (β k ) = Beta(β k ; γ 1,k , γ 2,k ) (4.15) q * (M k |Λ k ) = N (M k ; µ k , (κ k Λ k ) -1 ) q * (Λ k ) = W (Λ k ; L k , ν k )
Note that the hidden variables corresponding to the discrete dimensions remain unchanged. The variational distribution of the continuous part are simply Gaussian and Wishart distributions of variational parameters

{L k , ν k , µ k , κ k }.
The Variational update equations of the discrete part and beta random variables remain familiar as well:

γ 1,k = 1 + N ∑ n=1 φ n,k γ 2,k = η + N ∑ n=1 T ∑ l=k+1 φ n,l (4.16) k,i,v = α i,v + N ∑ n=1 φ n,k 1[x (n) i = v] ∀i ∈ D (4.17)
Algorithm 3 Variational Inference for the mixed type infinite mixture model

Input: x (1:N) , T, η, µ 0 , κ 0 , L 0 , ν 0 , α Initialize φ L (0) ss = -∞ for t = 1...∞ do Compute: γ (t)
µ k = κ 0 µ 0 + ∑ N n=1 φ n,k x (n) C κ k (4.18) κ k = κ 0 + N ∑ n=1 φ n,k ν k = ν 0 + N ∑ n=1 φ n,k + 1 (4.19) L -1 k = L -1 0 + κ 0 (µ k -µ 0 )(µ k -µ 0 ) T + N ∑ n=1 φ n,k (x (n) C -µ k )(x (n) C -µ k ) T (4.20)
Finally, the variational parameters representing the cluster probabilities φ have a different form, they depend naturally on the continuous and discrete dimensions simultaneously:

log φ n,k = - 1 2 |C| κ k + ν k (x (n) C -µ k ) T L k (x (n) C -µ k ) + ∑ i∈D ∑ v∈E i 1[x (n) i = v][ψ( k,i,v ) -ψ( ∑ v ∈E i k,i,v )] + 1 2 |C| ∑ i=1 ψ ν k + 1 -i 2 + log det(L k ) (4.21) + ψ(γ 1,k ) -ψ(γ 1,k + γ 2,k ) + k-1 ∑ l=1 [ψ(γ 2,l ) -ψ(γ 1,l + γ 2,l )] + const s.t T ∑ k=1 φ n,k = 1 and, φ n,y n = 1 ∀n ∈ D l
A proof of the previous equations is given in appendices B. 1 and B.2. Given all the variational fixed point equations for the parameters, we need only to iterate them until convergence of L ss in order to fit the model. Algorithm 3 summarizes the process.

. .

The predictive distribution

Given a trained model through the variational inference method, one needs to make predictions of classes using such model. Formally, given a new instance x (N+1) , one needs to determine the probability that the hidden cluster variable z (N+1) is equal to a certain value k, given all the observations x (1:N+1) . Meaning, we need to compute the following probability:

p(z (N+1) = k|x (1:N+1) ) ∝ p(z (N+1) = k, x (N+1) |x (1:N) ) (4.22)
Computing the previous quantity requires a marginalization over all the hidden variables except z (N+1) . Which leads to the follwing computation:

p(z (N+1) = k|x (1:N+1) ) ∝ p(z (N+1) = k, x (N+1) , B, M, Λ, β|x (1:N) )dMdΛdβdB ∝ p(x (N+1) |z (N+1) = k, B, M, Λ) (4.23) × p(z (N+1) = k|β)p(B, M, Λ, β|x (1:N) )dMdΛdβdB
In order to compute this integral, we need to compute the posterior of the parameters p(B, M, Λ, β|x (1:N) ). Unfortunately, as discussed in previous chapters a closed form of this posterior is intractable. However, variational inference gives us an approximation of this posterior: q * . Therefore, we can use this approximation to compute the integral [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF]:

p(z (N+1) = k|x (1:N+1) ) ∝ ∼ p(x (N+1) |z (N+1) = k, B, M, Λ) × p(z (N+1) = k|β)q * (B, M, Λ, β)dMdΛdβdB (4.24) ∝ ∼ E β∼q * [π k (β)] ∏ j∈D E B k,j ∼q * [Cat(x (N+1) j |B k,j )] × E M k ,Λ k ∼q * [N (x (N+1) C ; M k , Λ k )]
Furthermore, all the expectations of the previous equation can be computed in closed form (see Appendix B), and lead to the following final form of the predictive distribution:

Algorithm 4 Cluster prediction for future samples

Input: x (N+1) , φ, L, ν, µ, κ, γ, Compute: Θ k ∀k (4.26) Approximate: p(z (N+1) = k|x (1:N+1) ) ∀k (4.25) ẑ(N+1) = arg max k p(z (N+1) = k|x (1:N+1) ) return ẑ(N+1) p(z (N+1) = k|x (1:N+1) ) ∝ ∼ γ 1,k γ 1,k + γ 2,k ∏ l<k γ 2,l γ 1,l + γ 2,l ∏ j∈D k,j,x (N+1) j ∑ v∈E j k,j,v × κ k π(κ k + 1) p det (L -1 k + Θ -1 k ) -1 ν k +1 2 Γ p ( ν k +1 2 ) det(L k ) ν k 2 Γ p ( ν k 2 )
, (4.25)

where p represents the cardinal of the set C, and Θ k is p × p symmetric non-negative definite matrix defined by:

Θ -1 k = 1 κ k + 1 (κ k µ k + x (N+1) C )(κ k µ k + x (N+1) C ) T + κ k µ k µ T k + x (N+1) C x (N+1) C T (4.26)
A detailed derivation of the previous equation is given in Appendix B. Following from the previous computation, predicting the cluster of the N + 1 sample is simply done by taking the value of k representing the maximum of the probability. Algorithm 4 summarizes the process.

.

In this section, we demonstrate the ability of the semi-supervised mixed-type mixture model to learn to classify faults in a semi-supervised fashion.

. . Experimental setup and baselines

Our experimental setup in this chapter is the following, we consider the real-world dataset presented in chapter 3. This dataset is comprised of categorical variables (Alarms, status of equipment,...), and continuous variables (transmission and reception powers, bit error rates, ...). The categorical variables are treated in the same way as the previous chapter. As for the continuous variables, we scale their values between 0 and 1 using min-max scaling. We use the labeled part of the dataset. The evaluation metrics reported are computed over a held-out 10 % of the overall dataset. We use five cross-validation runs, and we report the mean and standard deviation of each metric.

In the previous chapter we showed how we can cluster patterns of faults in an unsupervised manner. In this chapter, we evaluate the semi-supervised classification capability of our model compared to some state-of-the-art baselines. The baselines considered are some of the best performing semisupervised learning methods, namely, Label Propagation [START_REF] Zhur | Learning from labeled and unlabeled data with label propagation[END_REF] and Transductive Support Vector Machines (TSVM) [START_REF] Joachims | Transductive inference for text classification using support vector machines[END_REF].

. .

Evaluating convergence of the model Same as the previous chapter the first task is to ensure that the update equations of the variational inference are correct and lead to model convergence, we report the semi-supervised evidence lower bound of equation (4.12) across updates of the variational parameters.

As shown in Figure 13, the evidence lower bound increases monotonically and converges to a plateau, which indicates that the equations derived are correct and that the model is well fitted. Second, we evaluate the classification accuracy of our model on the held-out datasets. We report the evolution of the classification accuracy across iterations. Same as the evidence lower bound the classification accuracy should increase and converge to a stable value. As shown in Figure 14, it is the case, the accuracy converges steadily over all runs to a value of about 94 % corresponding to 6 % error rate on the test set.

. . Evaluating the classification accuracy

Finally, we evaluate our model in comparison with the considered baselines, Label Propagation and Transductive Support Vector Machines (TSVM) with Linear and Gaussian kernels. In Table 4 we report the results of the comparison. The model is competitive with label propagation, with it scoring slightly better (one percent better) than our model on the dataset. Linear and Gaussian TSVMs perform much poorly compared to our model and label propagation. This could be a result of the assumptions behind the Transductive SVM. There is no guarantee that the classes in the data considered are separable linearly, or in the space generated by the kernel.

. &

In the semi-supervised learning literature, two approaches to assign classes to samples with prior knowledge have been extensively studied. The first approach supposes that all the different classes are present in the labeled part of the dataset, thus the number of classes is known. The model then attempts to propagate the labels to the unlabeled set, by respecting some notion of similarity. The previously introduced methods Transductive Support Vector Machines [START_REF] Joachims | Transductive inference for text classification using support vector machines[END_REF] and label propagation [START_REF] Bengio | Self-diagnosis technique for virtual private networks combining Bayesian networks and case-based reasoning[END_REF] and [START_REF] Zhur | Learning from labeled and unlabeled data with label propagation[END_REF] belong to this first approach. The Transductive Support Vector Machine model leverages the unlabeled data in the learning of the decision boundary between classes unlike classical support vector machines where only the labeled data is used.

Label propagation, on the other hand, is a graph based approach where the labeled and unlabeled data points represent nodes of the graph, and edges represent similarities between the nodes. The known labels are then propagated through the graph to label the unlabeled nodes. However, in some real world applications of data mining and information retrieval, the number of possible clusters is unknown and not all clusters are partially labeled, making the utility of such methods limited.

The second approach is based on the classical KMeans clustering algorithm, where constraints are added in order to guide the clustering, and improve the performance [START_REF] Wagstaff | Constrained k-means clustering with background knowledge[END_REF]. The constraints are constructed from the partial labels, must-link constraints assure that two samples must have the same label, and cannot-link constraints assure that two samples shouldn't be labeled the same. [START_REF] Basu | Active semi-supervision for pairwise constrained clustering[END_REF] proposed a Hidden Markov Random Field (HMRF) formulation of this model and later introduced distance learning to identify relevant features to the clustering [START_REF] Basu | A probabilistic framework for semi-supervised clustering[END_REF]. In this approach it is possible to find a hidden cluster not present in the partial labels. However, because it is based on the KMeans algorithm, the number of clusters must be known or some sort of model selection needs to be applied to estimate it. For a detail review of these approaches, we refer the reader to [START_REF] Chapelle | Semisupervised learning[END_REF].

By comparison, our approach does not suppose a fixed number of clusters but rather learns it directly from the data. More importantly, the modeling assumption are in accordance with the data types and structure, categorical variables remain as such and continuous variables are modeled as such. Thus, no distance learning is required to evaluate relevant features. In addition, the model in its semi-supervised form can adjust the parameters to the known labels of faults, however, if a new fault pattern emerges the Dirichlet Process mixture part automatically creates a new cluster for it.

. Although the infinite mixed-type semi-supervised mixture model can successively model the known labels of faults, and through its Dirichlet process mixture cluster new patterns, it remains that the model only functions correctly if the necessary information is directly correlated to the cluster hidden variable. Under this assumption, there is no need to create new features from the data in order to classify or cluster the faults correctly. The previous setting leads to good results, in highly specific applications with low dimensional datasets, that can be heavily preprocessed. In most real-world applications, we can collect hundreds if not thousands of variables and even with thorough preprocessing, noise remains, making shallow models perform poorly. In order to reach good performance in this case, one needs to use models capable of creating latent features relevant to each class, which naturally leads us to deep learning based methods.

. In this chapter, we extended the Infinite Mixture Model to the semi-supervised case with mixed type variables (continuous and categorical). We showed how to learn the model using variational inference, and how to make predictions using the model. In the next part of the manuscript, we discuss how to enlarge the scope of infinite mixture models to include feature learning using deep learning methods. These methods as we will show are necessary to deal with high dimensional noisy datasets, that require good feature creation to classify and identify patterns of faults. In chapters 3 and 4, we have introduced infinite mixture models for fault clustering and classification from network data. We generalized the approach to deal with categorical and continuous hidden variables in the presence of expert knowledge in the form of labels. A major limitation however of the previous approach is its inability to create new, and possibly more informative features for the classification task. In this chapter, we discuss how to combine the infinite mixture model approach with deep learning methods in order to solve this problem. We call this approach Deep Infinite Mixture Models.

Part II D E E P P R O B A B I L I S T I C G R A P H I C

Deep learning methods have demonstrated tremendous capability of learning representative features for classification tasks, especially for unstructured data such as images [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. Thanks to the ever increasing size of datasets that are labeled, deep learning can use this massive scale of data to reach good performance on a number of tasks ranging from classical supervised learning to reinforcement learning. Extending the paradigm to network data is very appealing, network data is often labeled through expert systems, is massive in terms of dimensionality and number of instances, thus deep learning methods can be applied to extract features in this domain as well.

The challenge of the diagnosis task in telecommunication networks is the discovery or clustering task. It is not sufficient to classify correctly known faults, but we need to also discover new patterns of faults. The literature of deep learning is still lacking in this task, and the subject is still widely unexplored and an active area of research. The closest deep learning paradigm is Zero-Shot Learning (ZSL), where some classes are known and others are hidden [START_REF] Xian | Zero-shot learningthe good, the bad and the ugly[END_REF]. However, in most applications of zero-shot learning extra semantic information is available (for example, a description of an image in words accompanying the image). This makes these approaches inapplicable for the diagnosis task, due to the fact that such information is often not available [START_REF] Larochelle | Zero-data learning of new tasks[END_REF][START_REF] Socher | Zero-shot learning through cross-modal transfer[END_REF].

In this chapter, we propose an approach based on deep neural networks and infinite Gaussian Mixture Models (GMM). Our approach combines the power of neural networks, and the exploration potential of infinite mixture models. The neural network identifies useful features for the fault discovery process (based on the labeled faults), then the infinite mixture model clusters the data based on the neural network representation. Thus, the fault discovery process is done on a more representative space of the faults.

We demonstrate our approach on data derived from GPON-FTTH networks, we show that the neural network is capable of identifying relevant features for the diagnosis task. Furthermore, we show that the clustering approach automatically identifies clusters of relevant faults. The remainder of the chapter is organized as follows:

• In section 5.2, we frame the problem of discovering new types of faults as a deep learning problem and we give an overview of the overall model.

• In section 5.3, we detail the feature extraction task as a deep learning classification task.

• In section 5.4 we introduce the infinite Gaussian mixture model for clustering new faults based on the extracted features.

• Finally, in section 5.5, we show how our approach can explore an unlabeled part of a more complex GPON-FTTH dataset in order to identify new types of faults that were previously considered as unidentified faults by an expert system.

. &

. .

Problem statement

Discovering and diagnosing faults from network data is extremely challenging due to three key properties of large telecommunication networks:

• Scale: when dealing with large telecommunication networks, the scale of the data in terms of dimensionality and the number of instances adds to the complexity of the diagnosis system, and restricts further the approaches that one can adopt to solve the diagnosis problem.

• Novelty: a second issue is the adaptability of the diagnosis system. Given a dynamic corpus of data gathered from the network, the system needs to identify previously known faults if they occur, in addition to clustering and identifying new patterns of faults.

• Noise: when dealing with network data, not all features and data dimensions will be relevant for a specific fault, hence adding noise to the relevant information. Therefore, an efficient diagnosis system needs to filter out the noise in some way in order to identify relevant features for the fault in question.

A learning-based efficient diagnosis system thus needs to identify and learn previously known patterns of faults, discover new patterns unknown to the experts or the expert system, and finally identify patterns relevant to actual faults on the network.

. .

Data specification

The first challenge encountered when dealing with network data is the multi-typed nature of variables. For example, when treating data gathered from the optical network, one finds continuous variables such as power levels, and variables taking discrete, categorical, or binary values such as alarms of some specific device. In order to create the data corpus, we start by normalizing the continuous variables by subtracting the minimum and then scaling the values (dividing by the maximum minus the minimum value). For the categorical variables, we use the one-hot encoding i.e. if a variable X i takes values in a set {c 1 , ..., c K }, then we create a representation in {0, 1} K , for which the k th column is equal to 1 if the variable X i is equal to c k . The resulting training data X ∈ R N×D contains N lines, each representing a customer's data vector, and D columns, each representing a variable of the network. Given the large amount of data gathered D is often very high (D >> 1, with often hundreds of variables). In addition to the data vector, for some customers, expert systems have identified the corresponding fault, thus providing us with a label. We denote these labels y ∈ R N , each line represents the label given to a customer. A label y n takes values in {-1, 1, ..., K}, where {1, ..., K} represent the known classes, and the value -1 represents an unknown class, as a convention. The key problem is to identify the labels of the class -1. This label could belong to the set of existing known faults (i.e. to the set of known classes {1, ..., K}) or a completely new set of faults occurring in the network. For what follows, we will denote x n ∈ R D the data vector corresponding to the n th customer, and suppose that the samples are independent and identically distributed. . .

Overview of the model

As shown in Figure 15, the diagnosis model is composed of three main blocks. The first one is the feature extractor, which will be a multi-layer perceptron or a neural network, responsible for compressing the data into a low dimensional feature vector containing the relevant information. The second one is the softmax unit responsible for classifying the known faults, its output is fed to the loss function in order to train the feature extractor. The final block is the clustering model responsible for identifying new patterns of faults based on the feature vector extracted.

The feature extractor represents the backbone of the system. Learning a good feature extractor is crucial for the performance of the diagnosis system. The feature extractor allows us to solve the scale and noise challenges simultaneously. The scale problem is solved by compressing the high dimensional raw data vector into a compact low dimensional manifold. The learned low dimensional manifold holds the information to correctly classify the known classes, thus creating noiseless features containing only the important information for each class.

The novelty challenge is tackled by the clustering model, which is in our case an infinite Gaussian mixture model. Based on learned features, the clustering model tries to identify new patterns or clusters of customer features, each cluster corresponding to a new unknown class. Then a final process of validation is done by an expert to identify relevant or correct new fault clusters. These correct clusters are then treated as labeled classes, the -1 labels being changed into positive values representing the new known classes. These new known classes are then re-fed to the feature extractor. In the following sections, we detail the model architectures of each block. . &

. .

Feature extractor architecture

The feature extractor in our case is a deep feed forward neural network, that takes as input a data vector x n ∈ R D and outputs a low dimensional feature vector denoted f n ∈ R p , where p << D represents the dimension of the feature manifold. As a remainder, we can think of a neural network as a nonlinear function, built by successive applications of linear layers and nonlinear activation functions. Formally, we denote the overall parameters of the neural network θ, and:

f n = g θ (x n ) = σ • f (L-1) • ... • σ • f (2) • σ • f (1) (x n ) (5.1)
where, each function f (l) represents the l th application of a linear or fully connected layer, defined by:

f (l) (h) = h T W l + b l (5.2)
W l , and b l represent the weight matrix and bias vector for the l th layer. • represents the composition of functions operator. Between each application of two linear layers, a non linear function is applied, denoted by σ, representing the ReLU activation function defined by

σ(h) = max(0, h) (5.3)
where the max is taken term by term.

For our application, we consider a neural network with 4 hidden layers (L -1 = 4), where the dimension of each layer is about half of that of the previous layer. The previous choices are based on common practices in the deep learning community. The choice of the number of layers is done by cross-validation, where we choose the lowest number of hidden layers for which the classification error rate does not change. An overview of the architecture of the feature extractor is given in Figure 16. We detail the loss function and the learning process in the next subsection.

. . Learning features by classification

In order to learn the parameters of the feature extractor, we define a classification task based on the K known classes. Given the feature vector of a customer denoted by f n , the softmax unit outputs the probability of each customer being in class k as:

P[y n = k|f n ] = e f T n W L,k +b L,k ∑ K j=1 e f T n W L,j +b L,j .
(5.4)

where W L , b L are the weights and biases of the last layer of the whole neural network (feature extractor followed by the softmax unit). The optimization of the parameters of the neural network is done by minimising the crossentropy loss function defined as:

C (θ) = - N ∑ n=1 K ∑ k=1 1[y n = k] log P θ [y n = k|f n ],
(5.5)

where θ represents the set of weights of the neural network. Learning θ is done by gradient descent and the standard backpropagation algorithm [START_REF] Goodfellow | Deep learning[END_REF]. Following the previous learning process, we learn a feature extractor and a classifier of the known classes {1, ..., K}. In the next section, we show how to identify outlier clusters representing new patterns of faults.

.

In the clustering process, the goal is to identify new patterns or clusters of unseen faults in unlabeled data (i.e. data with -1 label). The key trick in our method consists in applying the clustering process to the feature manifold.

The assumption is that the feature manifold is a more representative space in terms of filtered noise and relevant features than the original data space. Furthermore, applying the clustering process (step Given that the neural network learned highly expressive features for the known classes, it is highly likely that all information not relevant to the diagnosis task has been filtered. Furthermore, the low dimensional space of the last hidden layer is separated linearly with respect to each known class. Thus, outliers and new clusters should exist in a region outside the class regions.

3
The previous assumptions and deductions give guidance for the choice of the clustering model. The linearly separable nature of the feature space suggests that a simple Gaussian mixture model should be sufficient in order to identify clusters in it. Coupled with a Dirichlet process the Infinite Gaussian Mixture model can also learn the number of clusters directly without specifying it in advance, in the same manner as in the previous chapters.

. . Dirichlet Process Gaussian Mixture Model

We adopt the same formulation of the infinite Gaussian mixture as in chapter 4. The only exception is that the observable random variables are no longer the data samples, but the feature vector representation given by the neural network. The generative process is thus the following:

∀k β k ∼ Beta(1, η) π k = β k k-1 ∏ l=1 (1 -β l ) Λ k ∼ W (•; L 0 , ν 0 ) µ k |Λ k ∼ N (•|m 0 , (κ 0 Λ k ) -1 ) (5.6) z|π ∼ Cat(•|π) i.e. P[z = k] = π k f n |z n = k, µ, Λ ∼ N (•|µ k , Λ -1 k )
Where, {z 1:N , µ, Λ, β} represent the hidden random variables to infer, and f 1:N represent the feature vector of the observable random variables.

. . Variational Inference for the DPGMM

In the same manner as in the previous chapters, we can apply mean-field variational inference to perform inference of the hidden random variables given the observations, using the following minimization criterion:

q * = arg min q D KL [q||p(•, f 1:N )]
(5.7)

The optimal variational distributions have a familiar form, same as the last chapter given by:

q * (β k ) = Beta(β k ; γ 1,k , γ 2,k ) q * (µ k |Λ k ) = N (µ k ; m k , (κ k Λ k ) -1 ) q * (Λ k ) = W (Λ k ; L k , ν k ) q * (z n ) = Cat(z n ; φ n ) (5.8)
Where the variational parameters are updated iteratively using the following equations:

γ 1,k = 1 + N ∑ n=1 φ nk γ 2,k = η + N ∑ n=1 T ∑ l=k+1 φ nl
(5.9)

κ k = κ 0 + N ∑ n=1 φ nk ν k = ν 0 + N ∑ n=1 φ nk + 1 (5.10) m k = κ 0 m 0 + ∑ N n=1 φ nk f n κ k (5.11) L -1 k = L -1 0 + κ 0 (m k -m 0 )(m k -m 0 ) T + N ∑ n=1 φ nk (f n -m k )(f n -m k ) T (5.12) log φ nk = - 1 2 d κ k + ν k (f n -m k ) T L k (f n -m k ) + 1 2 d ∑ i=1 ψ ν k + 1 -i 2 + log det(L k ) + ψ(γ 1,k ) -ψ(γ 1,k + γ 2,k ) + k-1 ∑ l=1 [ψ(γ 2,l ) -ψ(γ 1,l + γ 2,l )] (5.13)
The details of the computation of q * are given in appendix A. A complete overview of the learning process is given in Algorithm 5. Table 5: Description of the different classes in the dataset.

We evaluate our approach on real network data. The dataset is built from the technical data for 86241 fiber customers.

In one part of the dataset, thanks to a legacy diagnosis system (in our case an expert system based on deterministic rules), the customers have been classified into 8 known types of faults (including a normal behavior) as described in For each customer, we collected 1824 features coming from different parts of the network that describe the state of the customer's line. These variables mainly characterize the properties of the FTTH GPON optical system, the Home Gateway (HG), some services used by the customer (e.g. TV, VoIP) and the Internet session (see Table 6). For instance, we have a set of variables describing the GPON properties such as the optical powers received (Rx) or transmitted (Tx) by the Optical Line Termination (OLT) at the central office, and by the Optical Network Termination (ONT) at the customer side, in addition to the properties of the ONT and OLT (temperature, voltage, version etc.), the alarms encountered by the OLT, and so on.

The dataset is thus composed of a mix of categorical variables and numerical variables. The numerical features are normalized in [0, 1] using min-max scaling. The categorical variables are encoded using the one-hot encoding, as explained in section 5.2, in order to obtain a numerical dataset. Thus the final dataset is composed of N = 86241 instances and D = 8297 variables.

In order to evaluate the classification task on the known classes, the labeled dataset is split into a training set and a test set, where each represent 80% and 20% respectively. The metrics for evaluations are reported over 5 cross validation runs, we report the mean and the standard deviation for each metric. The clustering evaluation is done on the entire unknown fault dataset in a single run.

. . Evaluation of the Feature Extractor

The first evaluation is that of the neural network representing the feature extractor. Extracting relevant features for the known classes is crucial for a good classification score, and for efficiently filtering noisy features for the subsequent clustering task. The key desirable feature of the hidden data manifold represented by the last hidden layer of the feature extractor is to be separable for each known class. In order to represent this data manifold, we compute the hidden features on the test set after each epoch (full pass of stochastic gradient descent on the training set), we then transform the features into a 2D representation using kernel principal component analysis. This procedure shows plots of the hidden feature points in a 2D representation based on a Gaussian proximity measure.

As shown in Figure 17, the hidden representation is disorganised at the first epochs of training the neural network, the classes are not well separated based on the similarity measure. However, as we converge to the optimum parameters, the hidden feature manifold becomes more separable with respect to each class (epoch 8). This suggests that the features extracted are representative of each class modeled. In addition, the structure is in the form of "lumps" which is coherent with Gaussian model assumptions. In order to demonstrate the effectiveness of the feature extractor further, we compare the raw two dimensional representation on the unknown faults dataset, and the two dimensional representation of the features extracted using the same dataset. As shown in Figure 18, the 2D representation on the extracted features is much more organized than the 2D representation on the raw data. Furthermore, the 2D representation of the features shows clear patterns of clusters unlike the raw data representations. This suggests that the feature extractor is capable of identifying features that permit us to separate and identify different clusters in the unknown fault dataset.

. . Evaluation of the Softmax Unit

The softmax unit classifies the known faults. As a standard evaluation process we report the loss function (5.5) on the test set across training iterations. We report the mean and standard deviation over 5 different runs. As shown in Figure 19 (left) the loss function converges to the set of parameters representing of the minimum of (5.5).

In addition, we report the classification accuracy on the test set. As shown in Figure 19 (right) the classification accuracy increases across training iterations and reaches a high value, showing that the classifier is well trained and robust to errors. This is demonstrated also by the confusion matrix between the known classes, reported in table 7. The standard deviation across the different validation runs is very small, suggesting that the classifier generalizes very well on unseen test set samples. In order to evaluate the clustering model, we begin by evaluating the fitting of the model to the unknown fault data. The latent infinite Gaussian mixture model is a probabilistic generative model, a standard approach of evaluating such models is to compute the log likelihood or the log probability of the data under the model across iterations. As shown in Figure 20, the log probability increases across iterations of updates of variational inference, and converges smoothly to a plateau, which suggests the convergence of the model.

The objective of the infinite latent Gaussian mixture model is to cluster patterns in the latent representation of unknown faults. In order to evaluate this clustering process, we report the clusters identified by the model in the 2D representation, across iterations of variational inference updates. As shown in Figure 21, the model starts with random clusters that do not fit the patterns in the unknown faults dataset. Across training iterations the ellipsoids representing the clusters fit with more accuracy the patterns seen in the 2D representation.

Furthermore, using the Dirichlet Process which learns the number of clusters, we see that across iterations of training of the clustering model, the number of clusters changes, in order to evaluate which number of clusters better fits the data. Thus, the model automatically identifies the number of clusters needed to converge to the optimal fit. . We have discussed how to combine the power of deep learning representations and the exploration capacity of Infinite Mixture Models, in order to discover faults in large noisy datasets. There are obvious limitations of the proposed combination, first, in the case where new faults are discovered a retraining process has to be re-implemented on the neural network side, which could prove expensive. Second, the unlabeled part only serves in the clustering part not the training of the deep model, thus wasting this unlabeled information that could prove useful for better performance and generalization. And last, the model is naturally sequential thus even the clustering process has to be repeated each time new labels are introduced. A natural answer to all these concerns, would be to consider deep generative models, where the feature learning and clustering are simultaneously done, making the model trainable end-to-end.

. In this chapter, we have introduced the Deep Infinite Mixture Model combining deep learning representations and Infinite mixture models in order to explore and identify faults in noisy datasets. The proposed model is sequential and composed of multiple bricks, making training and inference expensive. A natural extension can be proposed using deep generative models, where intermediate hidden random variables would represent the high level features necessary for the clustering and classification task. However, a major challenge in training deep generative models using variational inference, lies in the difficulty of backpropagating gradients through multiple layers of hidden random variables, with nonlinear activations. Specifically, the estimation of gradients of functions of the form L(θ) := E z∼q θ [log p(z)] is very difficult in deep generative models. Classical estimators often lead to high variance and make the learning of these models unstable. In the next chapter, we present an original approach developed in the thesis, to compute low variance estimates of these gradients.

G E N E R A L I Z E D S T O C H A S T I C B A C K P R O PA G AT I O N .

In the previous chapter, we discussed how to combine neural networks and mixture models to learn and cluster faults based on highly representative features of neural networks. A natural extension is to consider an end-toend deep generative model that models deep hidden features and cluster assignments simultaneously. The general formalism for these models are what is called deep belief networks [START_REF] Mnih | Neural variational inference and learning in belief networks[END_REF][START_REF] Neal | Connectionist learning of belief networks[END_REF], which are mainly deep neural networks with stochastic neurons instead of deterministic neurons. These models have become crucial in multiple domains, such as generative modeling [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Mnih | Neural variational inference and learning in belief networks[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF], deep reinforcement learning [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF], and attention mechanisms [START_REF] Mnih | Recurrent models of visual attention[END_REF].

Training these models amounts to performing gradient descent on functions of the form L(θ) := E z∼q θ [ f (z)] with respect to the parameters θ. The difficulty encountered in training these models arises in the computation of gradients of L(θ), which requires to backpropagate the gradient through the random variable z. One of the first and most commonly used methods is the score function or reinforce method [START_REF] Itu-T. Glynn | Optimization of stochastic systems via simulation[END_REF][START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF], that requires the computation and estimation of the derivative of the log probability function. For applications with high dimensionality however, it has been noted that reinforce gradients have high variance, making the training process unstable [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF].

Recently, significant progress has been made in tackling the variance problem. The first class of approaches that deals with continuous random variables are reparameterization tricks. In that case a standardization function is introduced, that separates the stochasticity from the dependency on the parameters θ. This allows to transport the derivative inside the expectation and to sample from a fixed distribution, resulting in low variance gradient [START_REF] Figurnov | Implicit reparameterization gradients[END_REF][START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Naesseth | Reparameterization gradients through acceptance-rejection sampling algorithms[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF][START_REF] Ruiz | The generalized reparameterization gradient[END_REF][START_REF] Titsias | Doubly stochastic variational Bayes for non-conjugate inference[END_REF]. The second class of approaches concerns discrete random variables, for which a direct reparameterization is not known. The first solution uses the score function gradient with control variate methods to reduce its variance [START_REF] Gu | MuProp: Unbiased Backpropagation for Stochastic Neural Networks[END_REF][START_REF] Mnih | Neural variational inference and learning in belief networks[END_REF]. The second consists in introducing a continuous relaxation admitting a reparameterization trick of the discrete random variable, thus being able to backpropagate low-variance reparameterized gradients by sampling from the concrete distribution [START_REF] Grathwohl | Backpropagation through the Void: Optimizing control variates for black-box gradient estimation[END_REF][START_REF] Jang | Categorical reparameterization with gumbel-softmax[END_REF][START_REF] Maddison | The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables[END_REF][START_REF] Tucker | Rebar: Low-variance, unbiased gradient estimates for discrete latent variable models[END_REF].

Although recent developments have advanced the state-of-the-art in terms of variance reduction and performance, stochastic backpropagation (i.e. computing gradients through random variables) still lacks theoretical foundation. In particular, the following questions remain open: How to develop stochastic backpropagation rules, where the derivative is transferred explicitly to the function f for a broader range of distributions? And can the discrete and deterministic cases be interpreted in the sense of stochastic backpropagation? In this chapter, we provide a new method to address these questions, and our main contributions are the following:

• We present a theoretical framework based on the link between the multivariate Fourier transform and the characteristic function, that provides a standard method for deriving stochastic backpropagation rules, for any distribution of practical interest discrete or continuous.

• We show that deterministic backpropagation can be interpreted as a special case of stochastic backpropagation, where the probability distribution q θ is a Dirac delta distribution, and that the discrete case can also be interpreted as backpropagating a discrete derivative.

• We generalize previously known estimators, and provide new stochastic backpropagation rules for the special cases of the Laplace, gamma, beta, and Dirichlet distributions.

• We demonstrate experimentally that the resulting new estimators are competitive with state-of-the art methods on simple tasks.

. &

Let (E, λ) be a d-dimensional measure space equipped with the standard inner product, and f be a square summable positive real valued function on E, that is,

f : E → R + , with E | f (z)| 2 λ(dz) < ∞.
Let q θ be an arbitrary parameterized probability density on the space E. We denote by ϕ θ its characteristic function, defined as: ϕ θ (ω) := E z∼q θ [e iω T z ]. We denote by f the Fourier transform of the function f defined as:

f (ω) := F { f }(ω) = E f (z)e -iω T z λ(dz). (6.1)
The inverse Fourier transform is given in this case by:

f (z) := F -1 { f }(z) = R d f (ω)e iω T z µ(dω), (6.2)
where µ(dω) represents the measure in the Fourier domain. In this chapter we treat the cases where E = R d for which µ(dω) = dω (2π) d , and the case where E is a discrete set, for which the measure µ is defined as:

µ(dω) = 1[ω ∈ [-π, π] d ] dω (2π) d .
Throughout the chapter, we reserve the letter i to denote the imaginary unit: i 2 = -1. To denote higher order derivatives of the function f , we use the multi-index notation [START_REF] Raymond | Elementary introduction to the theory of pseudodifferential operators[END_REF]. For a multi-index n = (n 1 , ..., n d ) ∈ N d , we define: To clarify the multi-index notation, let us consider the example where d = 3, and n = (1, 0, 2), in this case:

∂ n z := ∂ |n| ∂z n 1 1 ...∂z
∂ n z = ∂ 3 ∂z 1 ∂z 2 3 and, ω n = ω 1 ω 2 3 .
The objective is to derive stochastic backpropagation rules, similar to that of [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF], for functions of the form: L(θ) := E z∼q θ [ f (z)], for any arbitrary distribution q θ , discrete or continuous.

. Stochastic backpropagation rules similar to that of [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF] can in fact be derived for any continuous distribution, under certain conditions on the characteristic function. In the following theorem we present the main result of this chapter concerning the derivation of Generalized stochastic backpropagation rules.

Theorem 1. (Continuous Stochastic Backpropagation) Let f ∈ C ∞ (R d , R + )
, under the condition that ∇ θ log ϕ θ is a holomorphic function of iω, then there exists a unique set of real numbers {a n (θ)} n∈N d such that:

∇ θ L = ∑ |n|≥0 a n (θ)E z∼q θ [∂ n z f (z)] . (6.3)
Where {a n (θ)} n∈N d are the Taylor expansion coefficients of ∇ θ log ϕ θ (ω):

∇ θ log ϕ θ (ω) = ∑ |n|≥0 a n (θ)(iω) n . (6.4)
Proof. Let us rewrite L in terms of f :

L(θ) = q θ (z) f (z)λ(dz) = q θ (z)F -1 [ f ](z)λ(dz) = R d f (ω) E q θ (z)e iω T z λ(dz)µ(dω) Fubini's theorem = R d f (ω)ϕ θ (ω)µ(dω).
(6.5)

By introducing the derivative under the integral sign, and using the reinforce trick [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF] applied to ϕ θ , where ∇ θ ϕ θ (ω) = ϕ θ (ω)∇ θ log ϕ θ (ω), (6.5) becomes:

∇ θ L = R d f (ω)ϕ θ (ω)∇ θ log ϕ θ (ω)µ(dω). (6.6)
Under analyticity conditions of the gradient of the log characteristic function, we can expand the gradient term ∇ θ log ϕ θ (ω), in terms of Taylor series around zero as:

∇ θ log ϕ θ (ω) = ∑ |n|≥0 a n (θ)(iω) n .
(6.7)

Putting everything together, and replacing the characteristic function by its expression, the gradient of L becomes:

∇ θ L = R d f (ω) E q θ (z)e iω T z ∑ |n|≥0 a n (θ)(iω) n µ(dω)λ(dz). (6.8)
By rearranging the sums using Fubini's theorem a second time, we obtain the following expression for the gradient:

∇ θ L = E z∼q θ   F -1    ω → ∑ |n|≥0 a n (θ)(iω) n f (ω)    (z)   = ∑ |n|≥0 a n (θ)E z∼q θ F -1 ω → (iω) n f (ω) (z) = ∑ |n|≥0 a n (θ)E z∼q θ [∂ n z f (z)] .
(6.9)

Identically, we can follow the same procedure for discrete random variables. We suppose that q θ factorizes over disjoint cliques of the dependency graph, where each dimension z j takes values in a discrete space V al(z j ). In theorem 2 we derive the result concerning the discrete case.

Theorem 2. (Discrete Stochastic Backpropagation) Let E be a discrete space: E = ∏ d j=1 V al(z j ), and C the set of disjoint cliques of the dependency graph over z, that is,

q θ (z) = ∏ c∈C q θ (z c ) then, ∇ θ L = ∑ c∈C ∑ z c =z * c ∇ θ q θ (z c )E z -c ∼q θ [D f (z -c , z c )] .
(6.10)

Where:

• z * c : represents the normalizing assignment q θ (z

* c ) = 1 -∑ z c =z * c q θ (z c ). • D f (z -c , z c ) := f (z -c , z c ) -f (z -c , z * c ). (6.11)
Proof. The characteristic function for the factored distribution is given by:

ϕ θ (ω) = ∏ c∈C ϕ (c) θ (ω c ), where, ϕ (c) θ (ω c ) = ∑ z c =z * c q θ (z c )e iω T c z c + 1 -∑ z c =z * c q θ (z c ) e iω T c z * c . (6.12)
Thus the gradient of the log characteristic function becomes: .13) By plugging this expression to (6.6), we obtain:

∇ θ log ϕ (c) θ (ω c ) = ∑ z c =z * c ∇ θ q θ (z c ) e iω T c z c -e iω T c z * c ϕ (c) θ (ω c ) . ( 6
∇ θ L = ∑ c∈C ∑ z c =z * c ∇ θ q θ (z c ) ∏ c =c ϕ (c ) θ (ω c ) e iω T c z c -e iω T c z * c f (ω)µ(dω) = ∑ c∈C ∑ z c =z * c ∇ θ q θ (z c )E z -c ∼q θ [D f (z -c , z c )] . (6.14)
The estimator of (6.10) has been derived in the literature through Rao Blackwellization of the score function gradient, and it has been known under different names [START_REF] Asadi | Mean actor critic[END_REF][START_REF] Cong | GO gradient for expectation-based objectives[END_REF][START_REF] Titsias | Local expectation gradients for black box variational inference[END_REF]. Theorem 2 shows that the discrete case can also be seen as backpropagating a derivative of the function f , in this case a discrete derivative given by (6.11).

. -

Following from the previous section, we derive the stochastic backpropagation estimators for certain commonly used distributions.

The multivariate Gaussian distribution: In this case q θ (z) = N (z; µ θ , Σ θ ).

The log characteristic function is given by: log ϕ θ (ω) = iµ T θ ω + 1 2 Tr Σ θ i 2 ωω T . Thus by applying theorem 1, we recover the stochastic backpropagation rule of [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]:

∇ θ L = E z∼q θ ∂µ θ ∂θ T ∇ z f (z) + 1 2 Tr ∂Σ θ ∂θ ∇ 2 z f (z) , (6.15)
where, ∇ z and ∇ 2 z , represent the gradient and hessian operators. The multivariate Dirac distribution: q θ (z) = δ a θ (z), the log characteristic function of the Dirac distribution is given by: log ϕ θ (ω) = iω T a θ . Thus the stochastic backpropagation rule of the Dirac is given by:

∇ θ L = ∂a θ ∂θ T E z∼δ a θ [∇ z f (z)] = ∂a θ ∂θ T ∇ z f (a θ ), (6.16)
resulting in the classical backpropagation rule. In other words, the deterministic backpropagation rule is a special case of stochastic backpropagation where the distribution is a Dirac delta distribution. This result provides a link between probabilistic graphical models and classical neural networks.

The multivariate Bernoulli:

q θ (z) = ∏ d j=1 B(z j ; π (j) θ )
, where π (j)

θ = P[z j = 1]
. By applying theorem 2, we obtain the local expectation gradient of [START_REF] Titsias | Local expectation gradients for black box variational inference[END_REF]:

∇ θ L = d ∑ j=1 ∂π (j) θ ∂θ E z -j ∼q θ f (z -j , 1) -f (z -j , 0) . (6.17)
The multivariate categorical:

q θ (z) = ∏ d j=1 cat(z j ; π (j) θ )
, where the dimensions are independent and take values in the set {1, ..., K}. Similarly to the Bernoulli case, we obtain the following stochastic backpropagation rule:

∇ θ L = d ∑ j=1 K-1 ∑ k=1 ∂π θ (j) k ∂θ E z -j ∼q θ D f (z -j , k) . (6.18)
The Laplace distribution: q θ (z) = L(z; µ θ , b θ ), in this case the log characteristic function is the following: log ϕ θ (ω) = iµ θ ωlog(1 + b 2 θ ω 2 ), using the Taylor series expansion for the function x → 1 1-x , we get the following stochastic backpropagation rule for the Laplace distribution:

∇ θ L = ∂µ θ ∂θ E z d f dz (z) + 1 b 2 θ ∂b 2 θ ∂θ ∞ ∑ n=1 b 2n θ E z d 2n f dz 2n (z) . (6.19)
The gamma distribution: q θ (z) = Γ(z; k θ , µ θ ), the log characteristic function of the Gamma distribution is given by: log ϕ θ (ω) = -k θ log(1iµ θ ω). By expanding it using Taylor series of the logarithm function, we obtain the following stochastic backpropagation rule:

∇ θ L = ∞ ∑ n=1 1 n ∂k θ ∂θ + k θ µ θ ∂µ θ ∂θ µ n θ E z∼q θ d n f dz n (z) . (6.20)
The estimator of (6.20) gives a stochastic backpropagation rule for the gamma distribution and, hence also applies by extension to the special cases of the exponential, Erlang, and chi-squared distributions.

The beta distribution: q θ (z) = Beta(z; α θ , β θ ), in this case the characteristic function is the confluent hypergeometric function:

ϕ θ (ω) = 1 F 1 (α θ ; α θ + β θ ; iω).
A series expansion of the gradient of the log of this function is not trivial to derive. However, we can use the parameterization linking the gamma and beta distributions to derive a stochastic backpropagation rule.

Indeed, if ζ 1 ∼ Γ(α θ , 1) and ζ 2 ∼ Γ(β θ , 1) , then z = g(ζ 1 , ζ 2 ) = ζ 1 ζ 1 +ζ 2 ∼ Beta(α θ , β θ ).
By substituting in the gamma stochastic backpropagation rule, we obtain:

∇ θ L = ∞ ∑ n=1 1 n ∂α θ ∂θ E ζ 1 ,ζ 2 ∂ n f ∂ζ n 1 ζ 1 ζ 1 + ζ 2 + ∂β θ ∂θ E ζ 1 ,ζ 2 ∂ n f ∂ζ n 2 ζ 1 ζ 1 + ζ 2 .
(6.21) The Dirichlet distribution: q θ (z) = Dir(z; K, α θ ), following the same proce- dure, as for the beta distribution and using the following parameterization:

z k = ζ k ∑ K j=1 ζ j with, ζ k ∼ Γ(α (k)
θ , 1), we obtain: .22) .

∇ θ L = ∞ ∑ n=1 1 n K ∑ k=1 ∂α (k) θ ∂θ E ζ j ∀j ∂ n f ∂ζ n k ζ 1 ∑ K j=1 ζ j , ..., ζ K ∑ K j=1 ζ j . ( 6
The Generalized stochastic backpropagation gradient as presented in previous sections presents two major computational bottlenecks for non-trivial distributions. The first is the computation of infinite series, and the second is evaluating higher order derivatives of the function f . Depending on the application, the function f could be chosen in order to bypass the computational bottlenecks. A trivial example, is if the higher order derivatives of the function f vanish at a certain order: ∂ n z f = 0. Another example, is the exponential function f (z) = exp( T z). From the fact that it obeys the following partial differential equation ∂ f ∂z j (z) = j f (z), one can deduce that the stochastic backpropagation rule reduces in this case to:

∇ θ L = ∇ θ log ϕ θ (-i ) E z∼q θ [ f (z)] (6.23)
In most real world applications however, the infinite sum will not often reduce to a tractable expression such as that of the exponential. An example of this case is the evidence lower bound of a generative model with Bernoulli observations. In this case, the natural solution is to truncate the sum up to a finite order. The assumption (although it might be wrong), is that the components associated to higher frequencies of the gradient of the log characteristic function, do not contribute as much. In this case the gradient of the log characteristic function of (6.7), is truncated and becomes:

∇ θ log ϕ θ (ω) = ∑ n≤N a n (θ)(iω) n + o((iω) N ).
(6.24)

.

In our experimental evaluations, we test the stochastic backpropagation estimators of equations 6.19 and 6.20 for the gamma and Laplace distributions. In the case of the gamma estimator, we use toy examples where we can derive exact stochastic backpropagation rules without truncating the infinite sum. As for the Laplace stochastic backpropagation rule, we test the estimator in the case of Bayesian logistic regression with Laplacian priors and variational posteriors on the weights. We compare our estimators with the pathwise [START_REF] Jankowiak | Pathwise Derivatives for Multivariate Distributions[END_REF][START_REF] Jankowiak | Pathwise Derivatives Beyond the Reparameterization Trick[END_REF], and score function estimators, in addition to the weak reparameterization estimator in the gamma case [START_REF] Ruiz | The generalized reparameterization gradient[END_REF]. We do not use control variates in our setup, the goal is to verify the exactness of the proposed infinite series estimators and how they compare to current state-of-the-art methods in simple settings. In all our experiments, we use the Adam optimizer to update the weights [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with a standard learning rate of 10 -3 . In all the curves, we report the mean and standard deviation for all the metrics considered over 5 iterations.

. . Toy problems

In the toy problem setting, we test the gamma stochastic backpropagation rule following the same procedure as [START_REF] Mohamed | Monte carlo gradient estimation in machine learning[END_REF]. We consider the following cases: Toy problem 1: L(θ) = E z∼q θ ||z -|| 2 , where q θ (z) = ∏ d j=1 Γ(z j ; k j , µ j ), θ = {k, µ}, and = .49. In this case, we only need to compute the first and second order derivatives of the function f .

Toy problem 2: L(θ) = E z∼q θ ∑ d j=1 exp(z j ) , in this case, the infinite sum transfers to , which results in the following estimator:

∇ θ L = ∇ θ log ϕ θ ( i) E z∼q θ [ f (z)].
In figures 22 and 23 we report the training loss and log variance of the gradient across iterations of gradient descent for different values of the dimension d ∈ {1, 10, 100}. The stochastic backpropagation estimator converges to the minimal value in all cases faster than the other estimators and the variance of the gradient is competitive with the pathwise gradient in the learning regime (number of iterations < 10000).

. . Bayesian logistic regression with Laplacian Priors

We evaluate the Laplace stochastic backpropagation estimator using a Bayesian logistic regression model [START_REF] Jaakkola | A variational approach to Bayesian logistic regression models and their extensions[END_REF], similarly to [START_REF] Mohamed | Monte carlo gradient estimation in machine learning[END_REF]. In our case, we substitute the normal prior and posterior on the weights with Laplace priors and posteriors. We adopt the same notations of [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF], where the data, target and weight variables are respectively: andw. The probabilistic model in our case is the following:

x n ∈ R d , y n ∈ {-1, 1},
p(w) = d ∏ j=1 L(w j , 0, 1) p(y|x, w) = σ(yx T w), (6.25)
where σ represents the sigmoid function. We consider Laplacian variational posteriors of the form q θ (w) = ∏ d j=1 L(w j , µ j , b j ), with θ = {µ, b}. The evidence lower bound of a single sample is given by: 6.26) where the Kullback-Leibler divergence between the two Laplace distributions is the following:

L(x n , y n ; θ) = E w∼q θ log σ(y n x T n w) -D KL [q θ ||p], (
D KL [q θ ||p] = d ∑ j=1 |µ j | + b j e - |µ j | b j -log b j -1 . (6.27)
We test the model on the UCI women's breast cancer dataset [START_REF] Dua | UCI Machine Learning Repository[END_REF], with a batch size of 64 and 50 samples from the posterior to evaluate the expectation. In the case of the stochastic backpropagation estimator we truncate the infinite series for the scale parameter b of equation 6.19 to N = 4 and N = 8. In figure 24, we report the training evidence lower bound, the log variance of the gradient, and the accuracy computed on the entire dataset for the different estimators. The stochastic backpropagation estimator converges faster than the considered estimators and the variance is significantly lower. We also notice that the truncation level of the infinite series for the scale parameter has little effect on the outcome. In figure 25, we report the bias and variance of the estimator at different values of the truncation level, for a fixed parameter value during the training phase (epoch=100). The bias and variance do not vary much, with the truncation level in this case. This result confirms the intuition of neglecting higher frequencies presented in section 6.5. In addition, we compare the mean squared error between the Laplace stochastic backpropagation estimator and the score function and reparameterization estimators. As shown in figure 25 the mean squared error is small, thus the values of the gradients across iterations are close. However, the reparameterization gradient is closer to our estimator than the score function gradient, probably due to the fact that the reparameterization gradient is more stable and has lower variance.

. &

Computing gradients through stochastic computation graphs has received considerable attention from the community, due to its application in many fields. The first general approach that provides a closed form solution for any probability distribution is the score function method [START_REF] Itu-T. Glynn | Optimization of stochastic systems via simulation[END_REF][START_REF] Schulman | Gradient estimation using stochastic computation graphs[END_REF][START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF][START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF]. The main inconvenience of this approach, is that it results in high variance gradients when the dimension of the random variable becomes high. In order to bypass this issue, the second approach consisted in designing control variates to reduce the variance of the score function estimator [START_REF] Mnih | Neural variational inference and learning in belief networks[END_REF][START_REF] Paisley | Variational Bayesian inference with stochastic search[END_REF][START_REF] Ranganath | Black box variational inference[END_REF][START_REF] Tokui | Evaluating the variance of likelihoodratio gradient estimators[END_REF][START_REF] Weaver | The optimal reward baseline for gradientbased reinforcement learning[END_REF]. In addition to the score function gradient, it was proposed to use an importance weighted estimator instead of the classical score function with a multi-sample objective [START_REF] Burda | Importance weighted autoencoders[END_REF][START_REF] Mnih | Variational inference for monte carlo objectives[END_REF].

The second class of approaches is that concerning reparameterization tricks [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. Through the decoupling of the computation of the gradient from the expectation, reparameterization tricks have shown that they provide low-variance gradients using often a single sample. The issue for these methods is the necessity to find a reparameterization for each probability distribution. Certain distributions such as the Gaussian are easy to reparametrize but others like the gamma are not. In addition, discrete random variables do not admit an easy reparameterization as well. Recently, these issues has been partially solved through implicit reparameterization, the generalized reparameterization gradient, and the pathwise gradient [START_REF] Figurnov | Implicit reparameterization gradients[END_REF][START_REF] Jankowiak | Pathwise Derivatives Beyond the Reparameterization Trick[END_REF][START_REF] Ruiz | The generalized reparameterization gradient[END_REF]. For the discrete case, continuous relaxations that are reparameterizable have been proposed and combined with control variate methods [START_REF] Grathwohl | Backpropagation through the Void: Optimizing control variates for black-box gradient estimation[END_REF][START_REF] Gu | MuProp: Unbiased Backpropagation for Stochastic Neural Networks[END_REF][START_REF] Jang | Categorical reparameterization with gumbel-softmax[END_REF][START_REF] Maddison | The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables[END_REF][START_REF] Tucker | Rebar: Low-variance, unbiased gradient estimates for discrete latent variable models[END_REF].

Our approach, in contrast provides a new broad family of stochastic backpropagation rules derived using the Fourier transform. One interesting aspect of our approach is the fact that the weighting a n is separated from the expectation of the higher order derivatives of the function f . Thus the sampled variable does not intervene in the weighting in contrast to other methods such as reparameterization and pathwise gradients. In addition, by applying the derivative to the function f with respect to some variable z j , terms that do not depend on this variable are eliminated. Thus the other random variables in them do not need to be sampled, which results in lower variance.

It is worth noting that deriving stochastic backpropagation rules using the Fourier transform has been proposed for the Gaussian case [START_REF] Fellows | Fourier policy gradients[END_REF]. In our work, we extend it to non Gaussian distributions by way of the characteristic function, and exploiting the invariance of the functional inner product under Fourier transformation (Parseval's theorem).

. In this chapter, we presented a new method to compute gradients through random variables for any probability distribution, by explicitly transferring the derivative to the random variable using the Fourier transform. Our approach gives a framework to be applied for any distribution, where the gradient of the log characteristic function is analytic, resulting in a new broad family of stochastic backpropagation rules, that are unique for each distribution.

In the next chapter, we present the Infinite Latent Gaussian Model, an endto-end model to discover clusters of faults by learning representative deep hidden features. We show how Generalized Stochastic Backpropagation is used to perform variational inference to learn the model parameters.

D I R I C H L E T P R O C E S S D E E P L AT E N T M I X T U R E S F O R FA U LT I D E N T I F I C AT I O N .

In the previous chapter, we have showed how to backpropagate gradients through random variables efficiently. In what follows, we introduce an overall deep generative model capable of learning the hidden features and cluster assignment simultaneously. We show how to use Generalized Stochastic Backpropagation to learn these models, through variational inference. Our model is a generalization of the classical infinite mixture of chapters 3 and 4, with random variables representing hidden features and a Dirichlet process Mixture on the last hidden layer of random variables.

As discussed in previous chapters, nonparametric Bayesian priors, such as the Dirichlet Process (DP), have been widely adopted in the probabilistic graphical community. Their ability to generate an infinite amount of probability distributions using a discrete latent variable makes them ideally suited for automatic model selection. The most famous applications of the DP have been however limited to classical probabilistic graphical models such as Dirichlet Process Mixture Models and Hierarchical Dirichlet Process Hidden Markov Models [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF][START_REF] Fox | An HDP-HMM for systems with state persistence[END_REF][START_REF] Zhang | Stochastic variational inference for the HDP-HMM[END_REF].

Recently, deep generative models such as Deep Latent Gaussian Models (DLGMs) and Variational AutoEncoders (VAEs) [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF] have shown huge success in modeling and generating complex data structures such as images. Various proposals to generalize these models to the mixture and nonparametric mixture cases have been made [START_REF] Dilokthanakul | Deep unsupervised clustering with gaussian mixture variational autoencoders[END_REF][START_REF] Jiang | Variational deep embedding: An unsupervised and generative approach to clustering[END_REF][START_REF] Nalisnick | Approximate inference for deep latent gaussian mixtures[END_REF]Nalisnick and Smyth, 2016). Introducing such priors on top of the deep generative model can improve its generative capabilities, preserve class structure in the latent representation space, and offer a nonparametric way of performing model selection with respect to the size of the generative model.

The main challenge posed by such models lies in the inference process. Deep generative models with continuous latent variables owe their success mainly to the reparameterization trick [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. This approach provides an efficient and scalable method for obtaining low variance estimates of the gradient of the variational lower bound with respect to variational posterior parameters. Applying this approach directly to the variational posterior of the DP is not straight-forward, due to the fact that a reparameterization trick for the beta distributions is hard to obtain [START_REF] Ruiz | The generalized reparameterization gradient[END_REF]. One approach to bypass this issue have been proposed by (Nalisnick and Smyth, 2016), where the authors used the Kumaraswamy distribution [START_REF] Kumaraswamy | A generalized probability density function for double-bounded random processes[END_REF] as a higher entropy alternative for the beta distribution in the variational posterior. However, by deriving the nature of the variational posterior directly from the variational lower bound, we can show that the appropriate distribution is in fact the beta distribution.

In this chapter we provide an alternative treatment of the variational posterior of Dirichlet Process Deep Latent Mixtures (DP-DLM), where we combine classical variational inference to derive the variational posteriors of the beta distributions and cluster hidden variables, and neural variational inference for the hidden variables of the latent Gaussian model. This leads to gradient ascent updates over the parameters present in nonlinear transformations where Generalized Stochastic Backpropagation can be applied knowing the cluster assignment. As for the remaining parameters, closed-form solutions can be obtained by maximization of the evidence lower bound. The remainder of the chapter is organized as follows:

• In section 7.2, we introduce our proposition of Dirichlet process deep latent mixtures, with nonlinear stochastic hidden layers.

• In section 7.3, we show how to perform variational inference on the model, using Generalized Stochastic Backpropagation of chapter 6.

• In section 7.5, we demonstrate our model on the MNIST (LeCun and [START_REF] Lecun | MNIST handwritten digit database[END_REF] open source dataset, and more importantly, in the context of fault identification of our noisy dataset of chapter 5.

.

Generalizing deep latent models to the Dirichlet process mixture case can be obtained by adding a Dirichlet process prior on the hidden cluster assignments. We denote these cluster assignments by z. Following the assignment of a cluster hidden variable, a deep latent model is defined for

β z n h (l+1) n1 h (l+1) np l+1 h (l) n1 h (l) np l l ∈ {0, ..., L -1} n ∈ {1, ..., N} Figure 26:
The graphical representation of the generative process of the model, with the convention x = h (0) .

the assigned cluster similar to [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. We adopt the stick-breaking construction of the Dirichlet Process [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF]. The generative process of the model (Figure 26) is given by:

β k ∼ Beta(•; 1, η) π k = β k k-1 ∏ l=1 (1 -β l ) z n |π ∼ Cat(•|π) h (L) n |z n = k ∼ p(•|θ, z n = k) h (l) n |h (l+1) n , z n = k ∼ p(•|g (l) k (h (l+1) n )) ∀1 ≤ l ≤ L -1 x n |h (1) n , z n = k ∼ p(•|g (0) k (h (1) n ))
where h (l) n ∈ R p l is the l th layer hidden representation constructed by sampling a base distribution p(•|g

(l) k (h (l+1) n
)) with parameters computed using a nonlinear transformation g

(l) k (h (l+1) n ) = σ(W (l) k h (l+1) n + b (l) k ). W (l) k and b (l)
k are the weight and bias parameters for the l th hidden layer and k th cluster assignment. The base distribution in this formulation is general, but it can be specified to any classical distribution, for example a Gaussian, with parameters {m

(l) k , C (l) k } = g (l) k (h (l+1) n ), that is p(h (l) n |g (l) k (h (l+1) n )) = N (h (l) n ; m (l) k , C (l) k ).
The last hidden layer is generated using a prior p(•|θ, z n = k) with parameters θ. In the Gaussian example, we have p(•|θ,

z n = k) = N (h (L) n ; m (L) k , C (L) k ), with θ k = {m (L) k , C (L) k }.
We denote by η the concentration parameter of the Dirichlet process, which is a hyperparameter to be tuned manually. The emission distribution of the observable x n is usually chosen to be a normal distribution for continuous variables or the Bernoulli distribution for binary variables. We denote the parameters of the generative model by:

Θ = {θ 1:∞ , W (0:L-1) 1:∞ , b (1:L-1) 1:∞ }
The model thus has an infinite number of parameters due to the Dirichlet process prior. Furthermore, the posterior distribution of the hidden variables cannot be computed in closed-form. In order to perform inference on the model we need to use approximate methods such as Markov Chain Monte Carlo (MCMC) or Variational Inference. As discussed in previous chapters, MCMC methods are not suitable for high dimensional models, convergence of the Markov chain to the true posterior can prove to be slow and hard to diagnose [START_REF] Blei | Variational inference: A review for statisticians[END_REF].

In the next section, we show how to use variational inference and Generalized Stochastic Backpropagation to learn the model. We show that by choosing a suitable structure for the variational posterior, closed-form solutions can be obtained for the updates of the truncated variational posteriors of the beta distributions, and the variational posteriors of the cluster hidden variables. As for the remaining parameters, gradient ascent coupled with Generalized Stochastic Backpropagation are used to learn the variational posteriors of the hidden layers and model parameters.

. For a brief review of variational methods, we denote by x 1:N the N samples present in the dataset supposed to be independent and identically distributed. The log-likelihood of the model is intractable due to the required marginalization of all the hidden variables. In order to bypass this marginalization, we introduce an approximate distribution q Φ and use Jensen's inequality to obtain a lower bound [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF]:

l(Θ) = ln p Θ (x 1:N ) = ln ∑ z 1:N p Θ (x 1:N , z 1:N , h (1:L) 1:N , β)dh (1:L) 1:N dβ ≥ E z 1:N ,h (1:L) 1:N ,β∼q Φ ln p Θ (x 1:N , z 1:N , h (1:L) 1:N , β) q Φ (z 1:N , h (1:L) 1:N , β|x 1:N ) = L(Θ, Φ). (7.1)
We have shown that if the distribution q Φ is a good approximation of the true posterior, maximizing the evidence lower bound (ELBO) with respect to the model parameters Θ is equivalent to maximizing the log-likelihood (cf, chapter 2). For deep generative models, most state-of-the-art methods use inference networks to construct the posterior distribution (Nalisnick and Smyth, 2016;[START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. For deep mixture models with discrete latent variables, this approach leads to a mixture density variational posterior where the reparameterization trick requires additional investigation [START_REF] Graves | Stochastic backpropagation through mixture density distributions[END_REF]. Our approach combines standard variational Bayes and neural variational inference. We approximate the true posterior using the following structured variational posterior:

q Φ (z 1:N , h (1:L) 1:N , β|x 1:N ) = N ∏ n=1 L ∏ l=1 q ψ (l) z n (h (l) n |x n , z n ) × q φ n (z n |x n ) T ∏ t=1 q γ t (β t |x 1:N ), (7.2)
where T is a truncation level for the variational posterior of the beta distributions obtained by supposing that q(β T = 1) = 1 [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF].

We assume a factorized posterior over the hidden layers h

(1:L) n
, where the intra-layer dependencies are conserved.

. .

Deriving the variational posteriors q φ n and q γ t Deriving the nature of the posterior distributions of the hidden layers h

(1:L) n using the variational approach is intractable due to the nonlinearities present in the model. Thus, we take a similar approach to [START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF], and we assume that the variational posterior is specified by an inference network, where the parameters of the distribution are the outputs of auxiliary nonlinear transformations of parameters ψ (l)

t . For example, in the Gaussian case:

q ψ (l) t (h (l) n |x n , z n = t) = N h (l) n ; µ ψ (l) t (x n ), Σ ψ (l) t (x n ) , (7.3)
where

µ ψ (l) t (x n ), Σ ψ (l) t (x n ) are outputs of a one layer neural network of pa- rameters ψ (l)
t (Further details provided in Appendix A.2). In contrast to the hidden layers, we can use the proposed variational posterior of equation ( 7.2) to derive closed-form solutions for q φ n and q γ t . Let us consider the Kullback-Leibler definition of the ELBO L:

L(Θ, Φ) = -D KL [q Φ (•|x 1:N )||p Θ (•, x 1:N )] .
By plugging the variational posterior and isolating β t terms and z n terms, we can analytically derive the optimal distributions q γ t and q φ n maximizing L (see Appendix B):

q * γ t (β t |x 1:N ) = Beta(β t ; γ 1,t , γ 2,t ) q * φ n (z n |x n ) = Cat(z n ; φ n ),
where the fixed point equations for the variational parameters φ n and γ t are:

γ 1,t = 1 + N ∑ n=1 φ n,t γ 2,t = η + N ∑ n=1 T ∑ r=t+1 φ n,r (7.4) log φ n,t = const + E β∼q [log π t ] + E h (1:L) n ∼q ψ (1:L) t log p(x n , h (1:L) n |z n = t) + ∑ l H q ψ (l) t (•|z n = t, x n ) s.t. T ∑ t=1 φ n,t = 1
The fixed point equation of φ n,t , requires the evaluation of the expectation over the hidden layers, this can be performed by sampling from the variational posterior of each hidden layer and then forwarding the sample using the generative model:

E h (1:L) n ∼q ψ (1:L) t log p(x n , h (1:L) n |z n = t) ≈ 1 S S ∑ s=1 log p x n , h (1:
L)(s) n,t |z n = t (7.5) where: h (l)(s) n,t ∼ q ψ
(l) t (•|x n , z n = t).
A key insight here is the following: if a cluster t is incapable of reconstructing a sample x n from the variational posterior, this will reinforce the belief that x n should not be assigned to that cluster.

. . Generalized Stochastic Backpropagation

We next show how to perform stochastic backpropagation in order to maximize L with respect to the parameters Θ = {W

(0:L-1) 1:T , b (0:L-1) 1:T , θ 1:T } and Ψ = ψ (1:L) 1:T .
Similarly to the previous section, we isolate the terms in the evidence lower bound dependent on Ψ and Θ. We have:

L(Ψ, Θ) = const + ∑ n,t φ n,t ∑ l H q ψ (l) t (h (l) n |z n = t, x n ) + E h (1:L) n ∼q ψ (1:L) t log p Θ (x n , h (1:L) n |z n = t) . (7.6)
By taking the expectation over the hidden cluster variables z n , we obtain conditional expectations over the hidden layers h

(1:L) n knowing the cluster Algorithm 6 Variational Inference for the DP-DLM Input: x 1:N , T, η, α Initialize φ, Θ, Ψ while not converged do update: γ t ∀t (7.4) for each epoch do Θ ← Θ + α∇ Θ L Ψ ← Ψ + α∇ Ψ L end for update: φ n,t ∀n, ∀t (7.5) end while assignment. In order to apply Generalized Stochastic Backpropagation of chapter 6, let f (h

(1:L) n ) = log p(x n , h (1:L) n |z n = t).
Using Theorem 1, we have:

∇ Ψ L = ∑ |m|≥0 a m (Ψ)E h (1:L) n ∼q Ψ ∂ m h (1:L) n f (h (1:L) n ) (7.7)
In the Gaussian case, for example, that is:

q ψ ( l) t (h (l)
n |x n , z n = t) = N h (l) n ; µ ψ (l) t (x n ), Σ ψ (l) t (x n ) ,
we have:

∇ ψ (l) t L = E h (1:L) n ∼q Ψ        ∂µ ψ (l) t (x n ) ∂ψ (l) t   T ∇ h (l) n f (h (1:L) n ) + 1 2 Tr     ∂Σ ψ (l) t (x n ) ∂ψ (l) t   ∇ 2 h (l) n f (h (1:L) n )      . (7.8)
The gradient and hessian of the function f can be computed using standard automatic differentiation packages such as PyTorch [START_REF] Baydin | Automatic differentiation in machine learning: a survey[END_REF][START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF]. The expectation can be estimated using a standard Monte Carlo estimator, the same as equation (7.5). The entropy term of equation (7.6) can be computed analytically, thus gradients can be obtained using standard backpropagation (2.18). As for the Θ parameters, the expectation is independent of these parameters thus:

∇ Θ L = ∑ n,t φ n,t E h (1:L) n ∼q ψ (1:L) t ∇ Θ log p Θ (x n , h (1:L) n |z n = t) , (7.9)
which could be also estimated using standard backpropagation and a Monte Carlo estimator. Given the estimated values of these gradients, we can use standard gradient ascent to fit the parameters of the model in parallel with the cluster assignments. Algorithm 6 summarizes the process.

.

-( )

. .

SSL using the DP-DLM

In this section, similarly to chapter 4 we consider a partially labeled dataset x 1:N = D l ∪ D u , where D l = {x n , y n } n is the labeled part, y n represents the label of the sample x n , and D u represents the unlabeled part. The log likelihood can be divided for the labeled and unlabeled parts as:

l(Θ) = log p Θ (x 1:N ) = ∑ x n ∈D l log p Θ (x n ) + ∑ x n ∈D u log p Θ (x n ) = ∑ x n ∈D l log p Θ (x n , z n = y n ) + ∑ x n ∈D u log p Θ (x n ).
The last equation follows from the fact that p Θ (x n |z n = y n ) = 0. By dividing the labeled and unlabeled parts of the dataset, we can follow the same approach presented in section 7.3 in order to derive a variational inference algorithm. In this case, the fixed point updates and the gradient ascent steps remain unchanged if we set φ n,y n = 1 for a labeled x n sample.

. . The predictive distribution

In order to make predictions using the model, we need to evaluate the predictive distribution. Given a new sample x N+1 , the objective is to evaluate the following quantity p(z N+1 = k|x 1:N+1 ). In the same way as in chapter 4, we can use the variational posterior to approximate the true posterior, which in turn leads to simpler expectation terms:

p(z N+1 = k|x 1:N+1 ) ∝ p(z N+1 = k, x N+1 |x 1:N ) ∝ ∼ E β∼q [π k (β)] ×E h (1:L) N+1 ∼q ψ (1:L) k p x N+1 |g Θ (h (1:L) N+1 ), z n = k (7.10)
where g Θ (•) represents the forward pass over the generative model. The expectation with respect to the beta terms can be computed in closed-form as a product of expectations over the beta posteriors. The second expectation can be evaluated using a Monte-Carlo estimator same as equation (7.5).

kNN (k=5) DGM SB-DGM Ours 6.13 ± .13 4.86 ± .14 3.95 ± .15 2.90 ± .17 .

. -

In the previous section, we proposed a novel deep generative model. In order to demonstrate our model, we start by evaluating it on public opensource data. Given the general form of our model and its deep structure, we chose to evaluate it on open-source image data, following the standards of the generative modeling community.

. . Evaluation of the semi-supervised classification

We evaluate the semi-supervised classification capabilities of the model. We train a two hidden layer DP-DLM model (L = 2) on the MNIST dataset (Le-Cun and Cortes, 2010) with train-valid-test splits equal to {45000, 5000, 10000} similarly to (Nalisnick and Smyth, 2016), with 10 % labelisation randomly drawn. We train our algorithm until converenge using 5 cross-validation runs, and we evaluate our model on the test set. We report the mean and standard deviation of the classification error in percentages in Table 8. Our method produces a competitive score with existing state-of-the art methods: Deep Generative Models (DGM) (Kingma et al., 2014) and Stick-Breaking Deep Generative Models (SB-DGM) (Nalisnick and Smyth, 2016). Unlike the previous approaches, the loss was not up-weighted for the labeled samples. Figure 27 (left) shows the t-SNE projections (Maaten and Hinton, 2008a) obtained with 10 % of the labels provided. We notice that by introducing a small fraction of labels the class structure was highly preserved in the latent space.

. . Data generation and visualization

To further test our model, we generate samples for each cluster from the models trained on both the MNIST and SVHN [START_REF] Netzer | Reading digits in natural images with unsupervised feature learning[END_REF] datasets.

The MNIST model is trained in an unsupervised manner, and the SVHN model is trained with semi-supervision where we provide 1000 randomly generated labels. The samples obtained are represented in figure 28. For the unsupervised model, we notice that the clusters are representative of the shape of each digit. We plot the t-SNE (Maaten and Hinton, 2008b) projections of the MNIST test set of the unsupervised model in Figure 27 (right). We notice that the digits belonging to the same true class tend to group with each other. However, two groups of the same class can be very separated in the embedding space. The interpretation we can draw from this effect is that the DP-DLM tends to separate the latent space in order to distinguish between the variations of hidden representations of the same class. The clusters obtained are not always representative of the true classes which is a common effect with infinite mixture models. In a full unsupervised setting, data can be explained by multiple correct clustering results. This effect can simply be countered by adding a small supervision (figure 27 left). .

-

In this section of experiments, we fit our Dirichlet Process Deep Latent Mixture model in order to model clusters of faults in the noisy GPON-FTTH dataset of chapter 5. Our goal is to demonstrate the representational power in the hidden feature space, as well as, the classification of known faults and exploration of new ones.

. . Experimental setup

In this setting, we consider all known faults of the noisy dataset as presented in Table 5, in addition to the unknown faults. In training, we up-sample low represented classes to get a balanced dataset for each class. We use our Dirichlet process deep latent mixture to model the known classes and to cluster simultaneously the unknown fault data.

We consider four hidden layers in our probabilistic graphical model (L = 4), with Gaussian priors and variational posteriors. The visible variables x follow a Gaussian distribution if the dimension is continuous and a Bernoulli distribution with one-hot encoding if the dimension is categorical.

We implement Algorithm 6, with the following hyperparameters. The learning rate α is chosen to be equal to .001, decreasing exponentially across iterations. The parameters are learned using stochastic gradient descent with batch size equal to 256, with 100 epochs (passes over the dataset). The truncation level is set to T = 50 as an upper bound on the number of clusters. And the concentration parameter is set to η = 1/T (uniformative prior). . .

Model convergence

Following the same procedure as for previous chapters, in order to assess model convergence, we report the evidence lower bound of equation (7.1).

In Figure 29 (left), we can see that the evidence lower bound increases monotonically across iterations reaching a plateau. This indicates that the model is well-fitted and that the optimal variational parameters Ψ and Θ have been identified.

. . Semi-supervised classification results

In order to demonstrate the benefits of using deep representations, we compare the Dirichlet process deep latent mixture model and the infinite semi-supervised mixed-type mixture model of chapter 4. We evaluate the classification accuracy on the known labels of the noisy dataset for both models. We use a held out dataset unseen during training, containing a .1 fraction of the whole dataset. We run the experiment for five cross-validation runs, and we report the mean and standard deviation for the accuracy.

In Figure 29 (right), we report the classification accuracy across training epochs (passes over the dataset), for the eight known labels. As it can be noticed the DP-DLM outperforms considerably the mixed-type mixture. Furthermore, the shallow semi-supervised mixture model has really poor performance on this dataset. This result is to be expected seeing how the dataset contains many dimensions of noise (not relevant for faults). The DP-DLM has the capacity to compress and filter these dimensions into a more informative representation in the hidden layers. This capability is not available with classical shallow mixtures.

. . Cluster identification & visualisation

A second useful characteristic of our generative model is its capacity to give us a two dimensional representation of the hidden features. This provides us with a way of visualizing discovered clusters and learned classes. This is done by fixing the dimension of the last hidden layer to two (p L = 2).

In order to demonstrate the usefulness of this 2D representation, we report the evolution of clusters across learning iterations in Figure 30. For visualisation purposes we use three known classes, and two discovered clusters. As it can be seen the clusters are well separated in the feature space, which suggests that the model has learned useful features to separate these clusters.

. In this chapter, we introduced a novel model, the Dirichlet Process Deep Latent Mixture Model, for semi-supervised cluster discovery. We have shown how to use variational inference coupled with Generalized Stochastic Backpropagation to learn the model. We demonstrated the model on public open-source data, and GPON-FTTH data for the fault discovery task. We have shown that our model is competitive of state-of-the art methods in terms of semi-supervised classification. We have shown as well, that the model can be used to learn clusters of faults from network data.

The remaining issue after the cluster discovery process is the cluster interpretation problem. Interpretability of machine learning algorithms is still an open research field, with many questions. In the next and last chapter of the thesis, we will discuss the interpretability of the models introduced. We will show which models are interpretable by construction, and for the non-interpretable models, we will give a standard approach to interpret the clusters of faults obtained, so that network experts can fully exploit them.

FA U LT C L U S T E R I N T E R P R E TAT I O N F O R D O M A I N E X P E R T S

In the previous chapters, we have presented multiple models for fault diagnosis in GPON-FTTH data. We have tested their accuracy on known faults, and presented their clustering results on unknown and unlabeled instances corresponding to new faults. The clustering results present clear patterns in the data either using the shallow models or the deep models.

However, a problem remains, in order to exploit and make use of these clusters, we need to present them in a comprehensible way to network experts. In other words, either we need to specify directly a pattern of reoccurring variables for each fault, or retrace the fault using deterministic decisions of the sort "if-else" on the network variables.

Interpretability in machine learning is still considered as an open question, with many sub-areas of interesting research avenues. The multitude of models that perform predictions in a black-box way has moved this question to the center of the field, especially if we want to ensure reliability of our algorithms. The majority of the research focuses on making decisions or predictions of a model interpretable [START_REF] Fong | Interpretable explanations of black boxes by meaningful perturbation[END_REF][START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF][START_REF] Ribeiro | Why should i trust you?" Explaining the predictions of any classifier[END_REF]. Others, focus on specific models and mainly deep neural networks given their multiple applications [START_REF] Agarwal | Neural additive models: Interpretable machine learning with neural nets[END_REF][START_REF] Fong | Understanding deep networks via extremal perturbations and smooth masks[END_REF][START_REF] Selvaraju | Grad-cam: Visual explanations from deep networks via gradient-based localization[END_REF][START_REF] Yosinski | Understanding neural networks through deep visualization[END_REF].

Although these methods are widely popular in the machine learning community, their application for the network diagnosis task could prove limited. This is due to the fact, that in the machine learning community we would like to investigate how the model reached a decision, but in the network community we would like to determine why the model reached a decision. Investigating gradients and visualising hidden layers is useful for the machine learners but not for the network experts.

Therefore, in this chapter, we will go through each of the models presented and their results, and we will detail which are naturally interpretable, and which are not. For the deep models that are not easily interpretable, we will give a method based on decision trees to investigate the clustering results of these models.

.

The infinite mixture models introduced in chapters 3 and 4 directly model clusters from the data available. Thus in order to interpret the cluster, it suffices to inspect the emission distributions of the observable variables for each cluster on the fitted model. In what follows, we detail how to perform this task for each model.

. . The Infinite Categorical Mixture Model

Recall that the emission distribution of the Infinite Categorical Mixture model is defined as:

x i |z = k, B ∼ Cat(•|B k,i ) = ∏ v∈E i B 1[x i =v] k,i,v (8.1)
Meaning, for a certain value v ∈ E i , B k,i,v represents the probability of dimension i taking value v, supposing the cluster assignment is k:

B k,i,v = p[x i = v|z = k] (8.2)
Thus, in order to interpret cluster z = k, one only needs to use the posterior of B (in this case, the variational posterior q * ) to sample from the generative model by fixing z = k. Then, by counting the values occurring for each dimension, we get which occurrences contribute the most.

As an example on the dataset considered for experiments in chapter 3, we report the analysis for a clearly identified cluster in Figure 31. The cluster represents customers for which the account is suspended, this is clear by noticing that the values taken by the variables client account status equal to 0, and tv profile status being suspended. For more examples of explained clusters, we refer the reader back to chapter 3, subsection 3.6.3.

. . Continuous Dimensions and the Mixed-type Model

The Mixed-type Mixture model adds continuous dimensions to the categorical ones. The categorical ones thus can be treated in the same manner as the previous subsection, for the continuous dimensions however, we recall that:

x C |z = k, M, Λ ∼ N (•; M k , Λ -1 k ) (8.3)
where x C is the observable vector of continuous variables, following a Gaussian distribution of mean and covariance matrix M k , and Λ -1 k respectively, under the assumption that the cluster assignment z = k. In this case, if in the model assumptions we suppose that the dimensions are uncorrelated, that is, Λ -1 k = σ 2 k I, we can plot the mixture posterior distributions directly for each cluster and dimension. Thus, we have an interpretation of the mean values taken by each fault cluster.

As demonstrated in Figure 32, we plot the posterior distributions for three variables: Optical line termination transmission powers (OLTTx), Optical network termination transmission powers (ONTTx), and The down bit error count (BEC). All values are normalized between 0 and 1, -1 represents the case where an equipment is not responding thus providing no value. We plot the distributions for two recognizable clusters z = 1, and z = 3. As it can be seen, the cluster z = 1 represents the case where there is a problem with the OLT, where values are centered around -1. As for cluster z = 3, it represents normal behavior with (BEC equal to 0, and ONT equal to 1), however the OLT transmission power is slightly low ( approximately 0.6). . Deep models presented in chapters 5 and 7 are much more challenging to interpret. The introduction of non-linear hidden layers, either stochastic or deterministic, creates new features that help performance but heart interpretability. This is due to the fact that there is no way to reverse the process back to the original data. However, these models do provide us with cluster representing new patterns as demonstrated in previous chapters. Furthermore, the data is structured and can be read by experts of the domain. The challenge is to be able to link the clusters identified to patterns in the data.

In order to accomplish this task, we propose to train a decision tree classifier [START_REF] Breiman | Classification and regression trees[END_REF] on the unknown fault data with the cluster labels as the true labels of the classifier. Thus, the resulting decision tree classifier gives us a path of logical statements on the features leading to each cluster discovered. This provides a readable interpretable rule that can be exploited by experts of the network domain to identify the fault.

As an example, we revisit the clusters identified using the deep models in previous chapters, and we apply our method to them. In Figure 33, we give an example of such a decision path for the forth cluster discovered in the unknown fault dataset of chapter5. As it can be seen the decision rule mainly involves transmission and reception powers of the optical network termination of the GPON-FTTH network. This suggests a problem with Optical Network Termination (ONT). Another example where the decision is quite clear is given in Figure 34, where the issue is the link between the livebox (lb) or home router of the client and the Optical Network Termination (ONT).

C O N C L U S I O N

.

In this thesis, we have presented a class of probabilistic graphical models, based on Dirichlet process mixtures and deep generative models, in order to cluster faults from network data. Our models are capable of identifying clusters of faults without knowing a priori the number of clusters using the Dirichlet process. In addition, we presented how to extend the Dirichlet process mixtures to include latent hidden representations. These representations can learn features relevant to the fault identification task.

Furthermore, we introduced a novel method called Generalized Stochastic Backpropagation to compute low-variance gradients of the evidence lower bound with respect to the variational distribution. Our approach was used to train the proposed end-to-end Dirichlet Process Deep Latent Mixture Model for fault classification and discovery. Our approach also generalizes previous estimators and proposes new ones to other distributions.

We validated the presented models in the context of the fault diagnosis and discovery task. However, the models can be used in larger sense to other tasks in various domains. The theoretical results presented can be applied to other fields of machine learning ranging from reinforcement learning [START_REF] Sutton | Reinforcement learning: An introduction[END_REF] to attention mechanisms [START_REF] Mnih | Recurrent models of visual attention[END_REF][START_REF] Vaswani | Attention is all you need[END_REF]. The experimental results show that our models are able to classify with high accuracy known faults as well as clustering new patterns of faults.

The interpretation of new faults is straight-forward for shallow models. In this case, we can compute statistics relevant to each fault class, and investigate the variables occurring for each fault using these models. We can thus deduce the root cause of a fault cluster. As for deep models, where the interpretation is not straight-forward, we can use the clusters found as labels, to train highly interpretable models, such as decision trees, in order to retrace the path of decision leading to the cluster. This leads to an interpretation even for the deep models.

.

In this thesis we have discussed Dirichlet process mixtures and Dirichlet process deep mixtures for the fault classification and clustering problem. The extensions of these models and their application to different domains provide us with many interesting avenues for future research perspectives.

The first item treated is the clustering of faults using Dirichlet process mixtures. As discussed in chapters 3 and 4, we can see two main benefits with these models. The first is the learning of the number of clusters simultaneously with the cluster assignments. The second is their ability to treat continuous and categorical variables at the same time. However, the Dirichlet process is not the only non-parametric prior that can perform clustering and identify clusters. Other priors such as the mixture of finite mixtures prior (Miller and Harrison, 2015) can identify the number of clusters as well. Furthermore, the assumption behind the Dirichlet Process is that the size of clusters decreases as the number of clusters grows. This assumption follows intuition when dealing with clusters of faults, yet further experimental evaluations need to be performed to validate it. In the case where new faults keep occurring with homogeneous sizes of clusters, other priors might be better suited for the infinite models. Furthermore, we have assumed in the shallow infinite mixtures that the visible variables are independent. This assumption is clearly violated when dealing with network variables. These variables are often correlated and can be described by a Bayesian network with dependencies. A natural extension of the models would be to include these dependencies in the definition of the mixture models, and to re-derive the mean-field equations to include them. This could improve performance of the models in terms of accuracy and generalization. This conjecture could be verified in future work.

In part ii of the manuscript, we have discussed the extension of the infinite mixture models to include neural networks as feature extractors. Although their efficacy has been proven in chapters 5 and 7, we opted for simple fully connected architectures with a fixed number of hidden layers. The effect of these hyperparameters can be significant in terms of the effectiveness of the feature extractor. Furthermore, given that network data can easily be reconstructed into data originating from graphs, we can explore new architectures that are specifically invented to handle such data, such as, Graph Neural Networks [START_REF] Gori | A new model for learning in graph domains[END_REF][START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF][START_REF] Scarselli | The graph neural network model[END_REF]. Such models can prove effective in learning features specific to graph data, that are not accessible to classic fully connected layers, where the relationship between the elements of the layers is not taken into account.

Moreover we have introduced a new method of computing low-variance gradients through random variables in chapter 6. These computations arise in multiple domains ranging from deep learning, probabilistic models to reinforcement learning. However, we have not extensively tested our approach in different domains as is usually done in the machine learning community. This aspect still needs some investigation and could be another avenue of future research. This work may include building a generic library for generalized stochastic backpropagation that could be used in multiple applications. Also, for nontrivial distributions the approach requires computations of infinite sums of higher derivatives, which is often of great complexity. Thus further experimental and theoretical work has to be performed to deal with these aspects.

Finally, in this work we have focused mainly on the isolation step of the diagnosis problem, i.e. identifying a fault cluster and its cause. However, the detection step consisting of identifying anomalies can also benefit from the work done in this thesis. Indeed, the models introduced can be adapted to be applied to anomaly detection problems, for example, we can construct a probabilistic model that assigns normal data a high probability, and thus anomolous data is naturally assigned a low probability. In addition, the mitigation step can be addressed using approaches such as Markov Decision Processes and Reinforcement learning [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], where the theoretical methods introduced in this thesis could be of use. And last, all that is discussed can be applied to other areas of network management applications, such as Intrusion Detection Systems (IDS), repair agents, and more. All these are interesting areas of future research, that could widen the application of probabilistic graphical models to solve network management problems.

Part III A P P E N D I C E S

A V A R I AT I O N A L I N F E R E N C E .
In this appendix, we provide a review of mean field variational inference specifically the problem of solving equation (3.10):

q * = arg min q D KL q||p(•|x (1:N) ) . (A.1)
Let ζ = ζ 1:M represent the hidden random variables to infer, and x (1:N) represent the data vector of the observable random variable x. We have:

D KL q||p(•|x (1:N) ) = ∑ ζ q(ζ) log q(ζ) p(ζ|x (1:N) ) (A.2) = ∑ ζ q(ζ) log q(ζ) p(ζ, x (1:N) ) + log p(x (1:N) ) (A.3) = D KL q||p(•, x (1:N) ) + log p(x (1:N) ) independent of q (A.4)
Thus the minimization criterion of equation A.1 is equivalent to the following:

q * = arg min q D KL q||p(•, x (1:N) ) . (A.5)
Notice that the last minimizing criterion uses the joint distribution that is known, in contrast to the posterior distribution that is intractable. Let us denote by L the function to be minimized, we have .6) where :

L(q) = D KL q(ζ)||p(ζ, x (1:N) ) = ∑ ζ q(ζ) log q(ζ) p(ζ, x (1:N) ) = ∑ ζ m ∑ ζ -m q(ζ m )q(ζ -m ) eq (3.9) Factorisation log q(ζ m )q(ζ -m ) p(ζ, x (1:N) ) = ∑ ζ m q(ζ m )[-∑ ζ -m q(ζ -m ) log p(ζ, x (1:N) ) + log q(ζ m )] -∑ ζ m q(ζ m ) =1 H[q(ζ -m )] = ∑ ζ m q(ζ m ) -E ζ -m ∼q log p(ζ, x (1:N) ) + log q(ζ m ) -H[q(ζ -m )] = D KL [q(ζ m )|| f (ζ m )] -H[q(ζ -m ] (A
f (ζ m ) ∝ exp E ζ -m ∼q log p(ζ, x (1:N) ) (A.7)
Given the decomposition of L(q) to a term depending on q(ζ m ) and a term depending on q(ζ -m ), we can minimize with respect to each q(ζ m ), ∀m:

q * (ζ m ) = arg min q(ζ m ) L(q(ζ m )q * (ζ -m )) = arg min q(ζ m ) [D KL [q(ζ m )|| f * (ζ m )] -H[q * (ζ -m ]] = arg min q(ζ m ) D KL [q(ζ m )|| f * (ζ m )] where f * (ζ m ) ∝ exp E ζ -m ∼q * log p(ζ, x (1:N) ) = f * (ζ m )
Therefore, the solution to equation (A.1) is in the form of fixed point equations often referred to as the mean field update equations of the form:

log q * (ζ m ) = const + E ζ -m ∼q * log p(ζ, x (1:N) ) ∀m (A.8)
.

In the previous section we discussed mean field variational inference. In this section, we provide a review of neural variational inference for deep generative models, such as the model presented in chapter 7. As a remainder the evidence lower bound in this case is that of equation 7.1. The neural variational posterior is that of the hidden layers h

(1:L) n

, for the l + 1 layer in the Gaussian case we have:

q ψ (l+1) t (h (l+1) n |x n , z n = t) = N h (l+1) n ; µ ψ (l+1) t (x n ), Σ ψ (l+1) t (x n ) . (A.
{µ ψ (l+1) t (x n ), Σ ψ (l+1) t (x n )} = σ [1, x T n ]ψ (l+1) t (A.10)
further details of calculations

The main trick of the mean-field equation calculation is to recognize that all expectations and terms not dependent on z (n) are constant values, hence can be integrated in the constant term. Thus, we get:

log q * (z (n) ) = const + E B∼q * log p(x (n) |z (n) , B) + E β∼q * log p(z (n) |β) = const + E B∼q * ∑ i log Cat(x (n) i |B z (n) ,i,: ) + E β∼q * log π z (n) (β)
Given that z (n) is a categorical random variable, we deduce that:

q * (z (n) ) = Cat(z n ; φ n ) with: log φ n,k = const + E B∼q * ∑ i log Cat(x (n) i |B k,i,: ) + E β∼q * [log π k (β)] s.t T ∑ k=1 φ nk = 1 . . Computing q * (B k,i ):
Using the same principle, we have:

log q * (B k,i ) = const + E {z (1:N) ,β,B -k,i }∼q * log p(z (1:N) , B, β, x (1:N) ) = const + log p(B k,i ) + ∑ n E {z (n) ,B -k,i }∼q * log p(x (n) |z (n) , B) = const + log Dir(B k,i ; |x i |; α i ) + ∑ n φ n,k log Cat(x (n) i |B k,i ) = const + ∑ v∈E i α iv + ∑ n φ n,k 1[x (n) i = v] -1 log B k,i = log Dir(B k,i ; |x i |; i ) Where: i,v = α iv + ∑ n φ n,k 1[x (n) i = v] ∀v ∈ E i . . computing q * (β k ):
Following the same procedure for the variables β k , we get:

log q * (β k ) = const + E {z (1:N) ,β -k ,B}∼q * log p(z (1:N) , B, β, x (1:N) ) = const + log p(β k ) + ∑ n E {z (n) ,β -k }∼q * log p(z (n) |β) = const + log Beta(β k ; 1, η) + ∑ n T ∑ l=1 φ n,l E β -k ∼q * [log π l (β)] = const + (η -1) log(1 -β k ) + ∑ n T ∑ l=1 φ n,l E β -k ∼q * log β l + ∑ t<l log(1 -β t ) = const + (η -1) log(1 -β k ) + ∑ n φ n,k log β k + ∑ n T ∑ l=1 ∑ t<l φ n,l E β -k ∼q * [log(1 -β t )] = const + (η -1) log(1 -β k ) + ∑ n φ n,k log β k + ∑ n T ∑ t=k+1 φ n,t log(1 -β k ) = log Beta(β k ; γ 1,k , γ 2,k ) With: γ 1,k = 1 + N ∑ n=1 φ nk γ 2,k = η + N ∑ n=1 T ∑ l=k+1 φ nl
Finally, given all the variational distributions we need to compute the expectations for the φ equation, using standard results of expectations of Dirichlet and beta distributions we have:

E B∼q * ∑ i log Cat(x (n) i |B z (n) ,i,: ) = d ∑ i=1 ∑ v∈E i 1[x (n) i = v][ψ( k,i,v ) -ψ( ∑ v ∈E i k,i,v )]
and,

E β∼q * [log π k (β)] = ψ(γ 1,k ) -ψ(γ 1,k + γ 2,k ) + k-1 ∑ l=1 [ψ(γ 2,l ) -ψ(γ 1,l + γ 2,l )]
Where ψ represents the digamma function.

. -

The Mixed-Type model is an extension of the Infinite Categorical Mixture Model to incorporate continuous random variables and semi-supervision. The variational inference updates for the random variables B, β remain the same, new variational updates need to be computed for M, Λ, and the variational equation for the random variables z needs to be updated to incorporate the continuous part of the data.

. . Computing q * (z n ):

Using the same procedure as before, but for the new model, we have:

log q * (z (n) ) = const + N ∑ m=1 E z -n ,B,β,M,Λ∼q * log p(z (m) , B, β, M, Λ, x (m) ) = const + E B∼q * ∑ i∈D log Cat(x (n) i |B z (n) ,i,: ) + E β∼q * log π z (n) (β) + E M,Λ∼q * [log N (x (n) C ; M z (n) , Λ z (n) )]
. . Computing the variational posterior of the mean and precision matrix:

In order to compute q * (M k , Λ k ), we apply the same principle except in this case to the couple {M k , Λ k }, we have:

log q * (M k , Λ k ) = const + ∑ m E z (1:N) ,B,β,M -k ,Λ -k ∼q * log p(z (m) , B, β, M, Λ, x (m) ) = const + ∑ n,k φ n,k log N (x (n) C ; M k , Λ k ) + log N (M k , µ 0 , κ 0 Λ k ) + log W (Λ k ; L 0 ; ν 0 )
Using results of conjugate Gaussian and Wishart distributions, we can rewrite the last equation as:

log q * (M k , Λ k ) = log N (M k ; µ k , κ k Λ k ) + log W (Λ k ; L k ; ν k )
By correspondence between the terms, we can deduce, the following equations for the parameters:

µ k = κ 0 µ 0 + ∑ N n=1 φ n,k x (n) C κ k (B.1) κ k = κ 0 + N ∑ n=1 φ n,k ν k = ν 0 + N ∑ n=1 φ n,k + 1 (B.2) L -1 k = L -1 0 + κ 0 (µ k -µ 0 )(µ k -µ 0 ) T + N ∑ n=1 φ n,k (x (n) C -µ k )(x (n) C -µ k ) T (B.3)
Finally, in order to give a closed-form for the parameters φ, we must use the last distribution to compute the {M k , Λ k } expectation:

E M,Λ∼q * [log N (x (n) C ; M k , Λ k )] = - 1 2 |C| κ k + ν k (x (n) C -µ k ) T L k (x (n) C -µ k ) + 1 2 |C| ∑ i=1 ψ ν k + 1 -i 2 + log det(L k )
.

-

As a remainder, the predictive distribution of the mixed-type model is approximately proportional to :

p(z (N+1) = k|x (1:N+1) ) ∝ ∼ E β∼q * [π k (β)] ∏ j∈D E B k,j ∼q * [Cat(x (N+1) j |B k,j )] × E M k ,Λ k ∼q * [N (x (N+1) C ; M k , Λ k )]
In order to derive this quantity, we need to compute the different expectations, this is possible using results of conjugacy.

. . Computing the beta expectation:

E β∼q * [π k (β)] = E β t ∼Beta(•;γ 1,t ,γ 2,t )∀t [β k ∏ l<k (1 -β l )] = γ 1,k γ 1,k + γ 2,k ∏ l<k γ 2,l γ 1,l + γ 2,l . .
Computing the B expectation:

E B k,j ∼q * [Cat(x (N+1) i |B k,i )] = Dir(B k,i ; k,i , |x i |)Cat(x (N+1) i |B k,i )dB k,i = Γ(∑ v k,i,v ) ∏ v Γ( k,i,v ) ∏ v B k,i,v -1+1[x (N+1) i =v] k,i,v dB k,i,v = Γ(∑ v k,i,v ) ∏ v Γ( k,i,v ) ∏ v Γ( k,i,v + 1[x (N+1) i = v]) Γ(∑ v ( k,i,v + 1[x (N+1) i = v])) = k,i,x (N+1) j ∑ v∈E j k,i,v
. . Computing the normal expectation:

E M k ,Λ k ∼q * [N (x (N+1) C ; M k , Λ k )] = N (x (N+1) C ; M k , Λ k ) N (M k ; µ k , (κ k Λ k ) -1 )W (Λ k ; L k , ν k )dM k dΛ k = I(Λ k )W (Λ k ; L k , ν k )dΛ k
Where,

I(Λ k ) = N (M k ; µ k , (κ k Λ k ) -1 )N (x (N+1) C ; M k , Λ k )dM k = κ k 2π(κ k + 1) p det(Λ k ) exp - 1 2 tr Λ k Θ -1 k with, Θ -1 k = 1 κ k + 1 (κ k µ k + x (N+1) C )(κ k µ k + x (N+1) C ) T + κ k µ k µ T k + x (N+1) C x (N+1) C
T Thus we have:

E M k ,Λ k ∼q * [N (x (N+1) C ; M k , Λ k )] = κ k 2π(κ k + 1) p det(Λ k ) ν k +1-p-1 2 exp -1 2 tr[(L -1 k + Θ -1 k )Λ k ] dΛ k 2 ν k p 2 det(L k ) ν k 2 Γ p ( ν k 2 ) = κ k 2π(κ k + 1) p 2 (ν k +1)p 2 det (L -1 k + Θ -1 k ) -1 ν k +1 2 Γ p ( ν k +1 2 ) 2 ν k p 2 det(L k ) ν k 2 Γ p ( ν k 2 ) = κ k π(κ k + 1) p det (L -1 k + Θ -1 k ) -1 ν k +1 2 Γ p ( ν k +1 2 ) det(L k ) ν k 2 Γ p ( ν k 2 )
.

-

The evidence lower bound of equation 7.1 can be written as :

L(Θ, Φ) = -D KL [q Φ (•|x 1:N )||p Θ (•, x 1:N )] ,
in the following we develop this equation in order to derive the nature of the variational posteriors and their fixed point updates.

= const + ∑ z n q(z n ) q(h (1:L) n |x n , z n ) log p(x n |h (1:L) n )p(h (1:L) n |z n ) q(h (1:L) n |x n , z n ) dh (1:L) n + E β∼q [log p(z n |β)] -log q(z n ) = const -D KL [q(z n )||Cat(z n ; φ n )]
Hence, the optimal distribution q * φ n (z n |x n ) maximizing L satisfies:

q * φ n (z n |x n ) = Cat(z n ; φ n ) where, log φ n,k = const + q(h (1:L) n |x n , z n = k) log p(x n , h (1:L) n |z n = k) q(h (1:L) n |x n , z n = k) dh (1:L) n + E β∼q [log π k ] = const + E h (1:L) n ∼q ψ (1:L) k log p(x n , h (1:L) n |z n = k) + E β∼q [log π k ] + ∑ l H q ψ (l) k (•|z n = k, x n ) C A D D I T I O N A L T H E O R E T I C A L C O N T R I B U T I O N S . -
In the reinforcement learning setting, we would like to optimize the long term expected reward. We follow the notation of [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF] and we denote this expected reward ρ(π), where the policy π depends on parameters φ. The policy gradient theorem states:

∇ φ ρ = ∑ s d π (s)∇ φ E a∼π(•|s;φ) [Q π (s, a)] (C.1)
where Q π represents the action-value function, and d π represents the stationary distribution, or the discounted weighting of probabilities under policy π, depending on the problem formulation. In both cases, we can apply the results of generalized stochastic backpropagation to the expectation inside the sum. For discrete RL, where the action space is A = a (1) , ..., a (K) , we obtain the following expression for the gradient:

∇ φ ρ = ∑ s d π (s) K-1 ∑ k=1 ∂π(a (k) |s; φ) ∂φ Q π (s, a (k) ) -Q π (s, a (K) ) (C.2)
In the special case of episodic tasks, we would like to optimize the expected sum of rewards ρ(π) = E τ∼q φ [R(τ)]. The function R maps a trajectory τ = (s 1:T , a 1:T ) to the sum of rewards obtained from the environment: R(τ) = ∑ T t=1 r t (s t , a t ). The dynamics are encoded via the probability distribution q φ , where for a set of parameters φ, we have:

q φ (s 1:T , a 1:T ) = q(s 1 ) T ∏ t=1 π(a t |s t ; φ)q(s t+1 |s t , a t ) (C.3) additional theoretical contributions
In this case, by exploiting the dependency structure of the distribution q φ , on which the expectation is taken, equation (C.2) transforms into:

∇ φ ρ = T ∑ t=1 K-1 ∑ k=1 E s 1:t ,a 1:t-1 ∂π(a (k) |s t ; φ) ∂φ E s t+1:T ,a t+1:T T ∑ t =t r t (s t , a t )|a t = a (k) -E s t+1:T ,a t+1:T T ∑ t =t r t (s t , a t )|a t = a (K) (C.4)
In contrast to the Actor-Critic estimator [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF], the stochastic backpropagation estimator only needs the actor neural network of parameters φ to output the values of the policy. However, it does require looking ahead to get the expected cumulative reward of each action at time step t, which is an expensive task. In the supplementary materiel we provide a detailed derivation of equation (C.4).

. -

Stochastic backpropagation through mixture density posteriors have been studied recently [START_REF] Graves | Stochastic backpropagation through mixture density distributions[END_REF][START_REF] Roeder | Sticking the landing: Simple, lower-variance gradient estimators for variational inference[END_REF]. Following from the last section, we provide an alternative way of deriving stochastic backpropagation rules for mixture density distributions, as long as the base distributions of the mixture admits a stochastic backpropagtion rule. Let us consider a mixture density distribution of the form q φ (z) = ∑ K k=1 π φ k q (k) φ (z), where the component base distributions are q (k) φ for all k. In this case, we can rewrite the function L as :

L(φ) = E h∼cat(π φ ) E z∼q (h) φ [ f (z)] (C.5)
Following from equation (C.5), we can apply the categorical stochastic backpropagation rule on the first expectation, thus obtaining the following expression for the gradient:

∇ φ L = K-1 ∑ k=1 ∂π φ k ∂φ [g(k) -g(K)] + E h∼cat(π φ ) ∇ φ E z∼q (h) φ [ f (z)] (C.6) where g(k) := E z∼q (k) φ [ f (z)].
For the second term, we can re-apply stochastic backpropagation to derive a stochastic gradient for each base component and simply take the expectation over the class random variable.

. : -

x h (1) . . .

h (L) y

Figure 36: A hidden variable probabilistic model, where the observed variables are the data x and target y, with L hidden stochastic layers h (1:L) .

In this section, we explore the connection between neural networks and probabilistic graphical models following from the stochastic backpropagation rule of the Dirac delta distribution. To this end, let us consider the probabilistic graphical model of figure 36. The observed random variables in this model are denoted x and y representing the data and target variables. We place the analysis in a supervised learning context, but the argument is valid for unsupervised models as well. As usual the goal is to maximize the log likelihood for the data samples (x, y), which is intractable, given that we need to integrate over the hidden variables. However using variational inference, we can maximize an evidence lower bound of the form:

L(θ; x, y) = E h (1:L) ∼q θ (•|x) log p(y, h (1:L) , x) + H[q θ (•|x)] (C.7)
As suggested in the Dirac stochastic backpropagation rule, let us assume that the variational posteriors and priors are Dirac delta distribution of the form:

∀l : q θ (h (l+1) |h (l) ) = p(h (l+1) |h (l) ) = δ a (l+1) (W (l+1) T h (l) +b (l+1) ) (h (l+1) ) (C.8)
where, the a (l) , W (l) , and b (l) represent respectively the activation functions, the weights and biases for layer l, with the convention x := h (0) . Under these assumptions, the Kullback-Leibler divergence term is equal to zero, and the evidence lower bound reduces to the the log-likelihood of a classic neural network: L(θ; x, y) = log p(y|g θ (x)), with, g θ (x) = a (L) (W (L) (....a (1) (W (1) x + b (1) )...).

Thus, when using neural networks we are indirectly using a probabilistic graphical model and making the strong assumption that the hidden layers follow a parameterized Dirac distribution knowing the previous layer. We evaluate the Bernoulli and Categorical Stochastic Backpropagation estimators (BSB and CSB) of equations 6.17 and 6.18 on standard generative modeling benchmark tasks, using the MNIST and Omniglot datasets [START_REF] Lake | Oneshot learning by inverting a compositional causal process[END_REF][START_REF] Lecun | MNIST handwritten digit database[END_REF]. We use the REBAR, RELAX, and Gumbel-softmax (or Concrete) estimators as baselines for our comparison [START_REF] Grathwohl | Backpropagation through the Void: Optimizing control variates for black-box gradient estimation[END_REF][START_REF] Jang | Categorical reparameterization with gumbel-softmax[END_REF][START_REF] Maddison | The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables[END_REF][START_REF] Tucker | Rebar: Low-variance, unbiased gradient estimates for discrete latent variable models[END_REF]. The Bernoulli stochastic backpropagation is compared to the REBAR and RELAX estimators for three models: the sigmoid belief network of one and two stochastic hidden layers [START_REF] Neal | Connectionist learning of belief networks[END_REF] and the variational autoencoder. In this case, we adopt the same architectures as [START_REF] Grathwohl | Backpropagation through the Void: Optimizing control variates for black-box gradient estimation[END_REF]. The categorical stochastic backpropagation estimator is compared to the Gumbel-softmax estimator [START_REF] Jang | Categorical reparameterization with gumbel-softmax[END_REF][START_REF] Maddison | The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables[END_REF] using two models: a variational autoencoder and a single layer belief network with categorical priors. In this case, we set the dimension of the hidden layer to d = 20 and the number of modalities for each dimension to K = 10.

All models are trained using the ADAM optimizer (Kingma and Ba, 2014) using a standard learning rate α = 10 -4 and batch size of 100. We train the models for 500 epochs on the MNIST dataset and 100 epochs on the Omniglot dataset, longer learning epochs leads to overfitting and lower performance on the test sets for all estimators and models. We perform 5 iterations of training in all experiments and we report the mean and standard deviation of each performance metric considered.

For all models and estimators, we report the mean marginal test likelihood in tables 9 and 10 for both datasets. The test likelihood is estimated via importance sampling using 200 samples from the variational posterior. In all cases the stochastic backpropagation estimator, a control variate free method outperforms the baselines. In the case of the one layer sigmoid belief network the BSB estimator exhibits an increase of performance of about 4 nats in the case of the MNIST dataset and 10 nats in the case of the Omniglot dataset. We also vary the number of samples used to estimate the expectation in the stochastic backpropagation rule S ∈ {1, 5, 10}. We notice that using a single sample estimate does not hurt performance and leads to a faster training process.

We estimate the mean variance of the gradients w.r.t the parameters of the models using exponential moving averages of the first and second moments computed by the ADAM optimizer. The BSB estimator significantly outperforms the REBAR, RELAX estimators in terms of variance reduction with a difference of about 2 nats in the case of sigmoid belief networks, and 1 nat in the case of the categorical variational autoencoder on the mnist dataset, leading to a more stable training process as shown in figures 37 and 38. critic neural network considered also has two hidden layers of the same dimensionality with ReLU activation functions.

The actor and policy networks are trained using the ADAM optimizer, with a learning rate α = 10 -3 . We compute the cumulative reward at the completion of each episode of the environment, then we update the weights using a single step of gradient descent.

As shown in figure 39, the discrete stochastic backpropagation estimator outperforms the A2C estimator in the sense that it converges faster to the optimal policy. We also estimate the mean variance of the gradient of the policy network. As demonstrated at convergence the stochastic backpropagation estimator has a lower variance of about 1 nat compared to the A2C estimator, leading to more stable training process. Compte tenu des enjeux économiques, de nombreux travaux de recherche ont été menés pour construire et maintenir des systèmes de diagnostic de pannes opérationnels dans tous les domaines, allant des opérations de réseau aux installations et applications industrielles. Cependant, le champ de recherche pour résoudre la tâche de diagnostic est souvent spécifique au domaine, et peut même varier dans sa définition. Afin de donner une définition formelle de la tâche de diagnostic telle qu'elle est couramment admise dans la communauté de gestion de réseau, nous devons la diviser en trois étapes. La première étape est l'étape de détection, qui peut être effectuée de manière proactive, elle vise à décider si un client rencontre un problème et si des investigations supplémentaires sont nécessaires. La seconde est l'étape d'isolement, où l'objectif est d'identifier la cause racine du problème compte tenu de l'état technique actuel ou des données disponibles. Troisièmement, l'étape de réparation couvre toutes les actions nécessaires pour résoudre le problème.

La première étape, consistant à identifier si une panne s'est produite sur le réseau, est souvent associée à la détection d'un comportement anormal. Les méthodes souvent associées à cette tâche sont des approches de détection d'anomalies. Pour l'application spécifique que nous étudions tout au long de la thèse, à savoir le diagnostic de pannes des services d'accès à large bande, ces méthodes ne sont pas considérées : en effet, une panne est déclarée détectée à la suite d'un appel client.

Ainsi, nous nous intéressons plutôt à l'étape d'isolement, où l'on sait qu'une panne s'est produite sur le réseau, et il faut identifier la cause ou la nature de la panne. Nous nous intéressons plus spécifiquement au diagnostic des pannes des services FTTH (Fiber to the home). Considérant les différentes applications modernes de l'apprentissage automatique, l'objectif est de construire des modèles capables de classifier les pannes à partir des données du réseau GPON-FTTH issues du fonctionnement des réseaux Orange. De plus, l'enjeu est de construire des systèmes adaptables, dans le sens où, si une nouvelle panne apparaît sur le réseau (due à changements logiciels ou techniques) le modèle est capable de regrouper ensemble les instances de la nouvelle panne, étant donné qu'une signature de la panne existe dans les données rencontrées.

' -

Afin d'introduire la tâche de diagnostic de réseau comme un problème d'apprentissage automatique, nous devons d'abord discuter de l'état actuel du diagnostic pour Orange en particulier, et pour la communauté de recherche en général. Actuellement, dans la plupart des systèmes de diagnostic opérationnels du monde réel, la méthode préférée est les systèmes experts. Spécifiquement pour Orange, le système expert est appelé DELC (Diagnostic Expert de la Ligne Client).

Un système expert est un système basé sur des règles, où, étant donné une certaine condition sur les variables du réseau, une alarme est déclenchée pour alerter si une panne s'est produite et indiquer le cas échéant la cause racine de la panne. Cela a ses avantages et ses inconvénients. Le premier avantage est que ces règles peuvent être construites directement avec suffisamment de connaissances spécialisées, et elles sont interprétables et claires pour les autres experts du réseau ou du domaine. Cependant, les inconvénients sont nombreux. Le premier si le réseau est suffisamment grand, l'expertise et le temps requis pour construire le système expert sont énormes, le maintenir est encore plus difficile. Deuxièmement, si une nouvelle panne se produit en raison de changements dans le réseau, le système doit être revisité, ce qui peut s'avérer encore plus complexe. Troisièmement, puisque les règles sont déterministes, le système expert ne peut pas toujours conclure lorsque certaines variables sont manquantes.

La première approche qui a fait l'objet de nombreuses recherches pour la tâche de diagnostic fait appel aux réseaux bayésiens. Un réseau bayésien est un modèle probabiliste des variables de réseau avec dépendances. Le modèle peut être construit par un expert, les dépendances représentent des relations causales entre les variables du réseau. L'objectif d'un réseau bayésien est de modéliser la propagation des pannes à travers les variables du réseau bayésien, ainsi en inversant le processus, on peut identifier la cause racine. L'avantage de ces modèles est qu'ils peuvent être appris à partir de données (probabilités conditionnelles estimées à partir de données). Cependant, l'inconvénient reste le même que les systèmes experts : ils sont très spécifiques à l'application, si une nouvelle panne ou un changement survient sur le réseau réel, le réseau bayésien doit être revisité ce qui est souvent une tâche coûteuse.

Des méthodes classiques d'apprentissage automatique ont également été explorées pour la tâche de diagnostic. Une approche consiste à considérer le schéma d'apprentissage supervisé traditionnel, où un classificateur est appris à l'aide de données étiquetées selon les différents types de pannes connus. Ces classificateurs sont nombreux et peuvent aller de modèles simples tels que les classificateurs Naive-Bayes, aux Support-Vector Machines (SVM), et même aux réseaux de neurones.

Bien que ces méthodes soient attrayantes, en utilisant ces approches, nous aurions besoin de réapprendre notre modèle de classification avec de nouvelles données étiquetées pour de nouveaux types de pannes, qui ne sont généralement pas disponibles, surtout lorsque la panne est très récente. En utilisant les réseaux bayésiens, la prise en compte de nouveaux types de pannes, de nouvelles données ou de nouveaux équipements reviendrait à mettre à jour le graphe bayésien en ajoutant et/ou en supprimant des noeuds et des arêtes ainsi qu'en modifiant les probabilités conditionnelles des arêtes. Cela nécessiterait une grande expertise du réseau et peut être irréaliste si le réseau est trop complexe. Bien qu'il soit possible de dériver le graphe bayésien à partir de données, cela devient impossible en pratique lorsqu'il s'agit de grands ensembles de données composés de centaines ou de milliers d'entités. D'autre part, afin de découvrir de nouvelles pannes, il faut effectuer une exploration des données et trouver des clusters correspondant à de nouveaux types de pannes. L'approche standard actuellement, consiste pour un expert à étudier les données "à la main" afin d'identifier ces clusters. Les méthodes de clustering non supervisées peuvent être utilisées pour obtenir de bonnes performances. Cependant, l'utilisation de telles méthodes dépend souvent d'hypothèses de modèle qui pourraient nuire aux performances, telles que les métriques de distance dans les espaces de grande dimension. De plus, dans les méthodes non supervisées, les clusters peuvent être pertinents d'un point de vue technique (par exemple, regrouper des individus parce qu'ils utilisent le même type d'équipement), mais ce type de clustering n'apporte aucune information nouvelle qui pourrait être utilisée pour mettre en évidence une nouvelle panne. Le moyen évident de contourner ce problème est de recourir à une sélection fastidieuse de variables caractéristiques et à un prétraitement manuel des données, ce qui nécessite une connaissance approfondie de l'environnement en constante évolution et peut être coûteux en temps.

Dans cette thèse, nous présentons un cadre général basé sur des modèles probabilistes graphiques pour apprendre et identifier des clusters de pannes. Nous plaidons en faveur d'un type spécifique de modèles probabilistes graphiques, à savoir les mélanges à base de processus de Dirichlet ou les mélanges profonds à base de processus de Dirichlet. Ces modèles que nous proposons sont capables de clusteriser les pannes sans connaître le nombre de clusters (ou de pannes), ils peuvent identifier des signatures cachées pertinentes pour la tâche de diagnostic, et plus important encore, ils incluent des connaissances expertes en termes de supervision (ou semi-supervision) sous forme d'étiquettes pour guider et améliorer l'identification des pannes à partir des données.

Cette thèse est divisée en deux parties: dans la partie I, nous introduisons des modèles probabilistes graphiques superficiels. Plus précisément, des modèles de mélange infini capables d'identifier des clusters de pannes à partir de données de réseau, de manière non supervisée et semi-supervisée. Dans la partie II, nous faisons évoluer ces modèles pour inclure des représentations cachées profondes. Le but des représentations cachées est de construire des signatures pertinentes pour la tâche de diagnostic. Ces fonctionnalités sont conçues pour traiter des données à grande échelle bruitées et complexes. Combiné avec les mélanges infinis, nous pouvons utiliser ces représentations pour obtenir de meilleurs résultats sur de tels ensembles de données. Afin de guider la lecture du manuscrit, les contributions apportées dans la partie I tentent de répondre aux questions de recherche suivantes:

Question 1 : Peut-on développer des modèles capables d'identifier des clusters adaptés aux données du réseau, notamment des données catégorielles, sans connaître a priori le nombre de clusters ? Notre contribution pour répondre à cette question est le modèle de mélange infini catégoriel. Le modèle est présenté au chapitre 3. Nous étendons la définition des mélanges à base de processus de Dirichlet au cas des variables catégorielles multidimensionnelles et nous montrons comment utiliser l'inférence variationnelle de type champ-moyen pour apprendre le modèle efficacement. Nous montrons que notre modèle est capable d'identifier des clusters de pannes sans connaître a priori le nombre de clusters et nous montrons comment les hypothèses du modèle conduisent à de meilleures performances en termes de précision.

Question 2 : Peut-on étendre les mélanges à base de processus de Dirichlet pour traiter simultanément des variables continues et catégorielles, et prendre également en compte les connaissances expertes en termes de supervision ?

Dans cette contribution, nous étendons le modèle de mélange catégoriel infini à un modèle de type mixte, qui peut traiter des variables continues et catégorielles. Le modèle prend également en compte la semi-supervision d'experts sous forme de labels. Nous évaluons le modèle sur les données du réseau et nous montrons ses très bonnes performances comparées à des modèles références de l'état de l'art. Ce travail est présenté dans le chapitre 4.

Nous avons introduit des modèles probabilistes graphiques superficiels, sous la forme de modèles de mélanges infinis, qui pourraient être appliqués à des données de type mixte. Un problème demeure, ces modèles ne peuvent pas créer de combinaisons complexes de données qui seraient nécessaires pour découvrir des clusters. En d'autres termes, aucune création de signatures n'est possible. Cependant, dans la plupart des applications du monde réel, les données sont bruitées et de grande dimensionnalité. Ainsi, la création de signatures ou "features" est nécessaire. Dans la partie II, nous développons des modèles et des méthodes afin d'atteindre cet objectif, et la principale question de recherche abordée est la suivante: Question 3: Pouvons-nous construire des représentations ou des signatures très pertinentes pour la tâche de diagnostic, afin d'améliorer les performances ? Afin de répondre à cette question, dans le chapitre 5, nous combinons des représentations de l'apprentissage profond et des modèles de mélange infini pour effectuer la découverte de pannes dans un ensemble de données à grande échelle partiellement étiqueté. Nous parlons pour cette approche de modèles de mélanges infinis profonds. Nous créons un réseau de neurones profond pour reconnaître les pannes étiquetées et construire des représentations cachées de faible dimension pertinentes pour la tâche de diagnostic. Ensuite, un modèle de mélange gaussien infini est appliqué sur la partie non étiquetée afin de reconnaître les clusters de pannes sur la base de la représentation cachée. Nous évaluons le modèle sur un ensemble de données GPON-FTTH du monde réel à grande échelle, nous montrons que le modèle peut identifier des classes de pannes connues avec une grande précision, et que la représentation cachée sur les pannes non étiquetées conduit à des clusters reconnaissables.

Le deep learning et les mélanges infinis vivent dans deux univers de méthodes différents. Leur combinaison n'est pas naturelle, dans le sens où nous n'avons pas de modèle global de bout en bout, capable d'apprendre des signatures et de faire du clustering. Ainsi se pose naturellement la question suivante:

Question 4: Peut-on construire un modèle probabiliste profond de bout en bout capable d'apprendre des représentations cachées et d'effectuer du clustering ? Il s'avère que cela est possible grâce aux modèles génératifs profonds, qui sont une classe de modèles probabilistes graphiques. Le problème des modèles génératifs profonds est le processus d'apprentissage qui implique le calcul des estimateurs de gradients qui sont souvent de grande variance. Dans nos contributions Rétropropagation stochastique généralisée et Rétropropagation stochastique par transformées de Fourier, nous présentons une méthode pour dériver des estimateurs à faible variance de ces gradients. Ce travail est présenté au chapitre 6.

Dans le chapitre 7, nous utilisons la rétropropagation stochastique généralisée pour former notre modèle de mélange latent profond à base de processus de Dirichlet de bout en bout. Nous évaluons notre modèle sur des ensembles de données ouverts et sur des données GPON-FTTH à grande échelle du monde réel pour la tâche de diagnostic et de découverte des pannes.

Enfin, au chapitre 8, nous discutons de l'interprétabilité de nos modèles. Nous montrons que les modèles superficiels sont naturellement interprétables. Les modèles profonds, cependant, sont plus difficiles en termes d'interprétabilité. Nous discutons de la façon d'utiliser les arbres de décision pour retracer la décision menant à un cluster de pannes, la rendant compréhensible pour un expert du domaine.

A C K N O W L E D G E M E N T S

First and foremost, I would like to thank Sandrine Vaton and my team at Orange Labs for giving me the opportunity to pursue a PhD under their supervision. I am tremendously grateful for their support, and for the continuous encouragement to pursue challenging ideas. Special thanks also to my Orange Labs advisers Joachim Flocon-Cholet and Stéphane Gosselin for their support, encouragement, and valuable discussions and insights through the years. My time at Orange Labs has truly been a wonderful and fulfilling experience thanks to them all. I would also like to thank my family who has supported me through the years across my academic journey, and especially during the unusual circumstances under which the thesis took place. Without the many welcomed distractions, calls, jokes, the stress of the situation at the time would have surely gotten the better of me. I would like to thank my sister Meriem for her unwavering optimism and support, without her I would be a different person today. Finally, I would like to thank my parents, Omar and Khadija for their support, encouragements, and the sacrifices that they had to endure for me to be where I am today. A debt that I would never be able to repay! Sincerely, Amine Echraibi

- 3 . 1

 31 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Mixture Models as Probabilistic Graphical Models . . . 3.2.2 The Dirichlet Process . . . . . . . . . . . . . . . . . . . . 3.2.3 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Infinite Categorical Mixture Models . . . . . . . . . . . . . . . 3.3.1 Structure and Definition . . . . . . . . . . . . . . . . . . 3.3.2 Learning and Inference . . . . . . . . . . . . . . . . . . . 3.4 Variational Inference for the ICMM . . . . . . . . . . . . . . . . 3.4.1 The mean field approximation . . . . . . . . . . . . . . 3.4.2 Optimal Variational Distributions and Update Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Baselines and Metrics . . . . . . . . . . . . . . . . . . . 3.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Clustering results on the synthetic dataset . . . . . . . 3.6.2 Variational inference vs Gibbs sampling . . . . . . . . . 3.6.3 Pattern recognition on real-world data . . . . . . . . . . 3.7 Discussion & Prior Work . . . . . . . . . . . . . . . . . . . . . . 3.8 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -ix x 8.1.2 Continuous Dimensions and the Mixed-type Model . 91 8.2 Interpretability of Deep Models . . . . . . . . . . . . . . . . . . 92 93 9.1 Thesis work summary . . . . . . . . . . . . . . . . . . . . . . . 93 9.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 97 99 a.1 Mean Field Variational Inference . . . . . . . . . . . . . . . . . 99 a.2 Neural Variational Inference . . . . . . . . . . . . . . . . . . . . 100 103 b.1 Variational Inference for the Infinite Categorical Mixture Model103 b.1.1 Computing q * (z n ): . . . . . . . . . . . . . . . . . . . . . 103 b.1.2 Computing q * (B k,i ): . . . . . . . . . . . . . . . . . . . . 104 b.1.3 computing q * (β k ): . . . . . . . . . . . . . . . . . . . . . 104 b.2 Variational Inference for the Mixed-Type Model . . . . . . . . 105 b.2.1 Computing q * (z n ): . . . . . . . . . . . . . . . . . . . . . 106 xii contents b.2.2 Computing the variational posterior of the mean and precision matrix: . . . . . . . . . . . . . . . . . . . . . . b.3 Predictive distribution of the mixed-type model . . . . . . . . b.3.1 Computing the beta expectation: . . . . . . . . . . . . . b.3.2 Computing the B expectation: . . . . . . . . . . . . . . . b.3.3 Computing the normal expectation: . . . . . . . . . . . b.4 Variational inference for the Dirichlet Process Deep Latent Mixture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . b.4.1 Computing q γ t (β t |x 1:N ): . . . . . . . . . . . . . . . . . . b.4.2 Computing q φ n (z n |x n ): . . . . . . . . . . . . . . . . . . . c.1 Stochastic Backpropagation for Reinforcement Learning . . . c.2 Stochastic Backpropagation for mixture density distributions c.3 The Dirac distribution: the link between neural networks and PGMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . d.1 Experiments using discrete stochastic backpropagation . . . . d.2 Reinforcement Learning Application . . . . . . . . . . . . . . .

Figure A hidden variable

  FigureA hidden variable probabilistic model, where the observed variables are the data x and target y, with L hidden stochastic layers h(1:L) . . . . . . . . . . . . . . . 113 Figure The training evidence lower bound on the MNIST training set (top) and the log variance of the gradient (bottom) over 5 runs. Comparison with the REBAR and RELAX estimators. . . . . . . . . . . . . . . . . . . 116 Figure Training evidence lower bound and the log variance of the gradient for the categorical VAE on the MNIST dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Figure Mean and standard deviation of the cumulative reward (left) and the log variance of the gradient for the actor network (right) over 5 iterations of the training process for the cart-pole task. . . . . . . . . . . . . . . 118

1 I

 1 time of one epoch of training on the mnist dataset per estimator, per model. . . . . . . . . . . . . xvi

Figure 1 :B

 1 Figure 1: Summary of the structure of the document

Figure 2 :Figure 3 :

 23 Figure 2: An example of a PGM.

Figure 4 :

 4 Figure 4: Graphical representation of a mixture model.

Figure 5 :

 5 Figure 5: Graphical representation of the model in plate notation.

Figure 6 :

 6 Figure 6: Expert Bayesian Network of the FTTH GPON (Tembo et al., 2017). Colored nodes are possible root causes whereas non-colored nodes are observations.

Figure 7 :

 7 Figure 7: Network scope for fault pattern extraction from real operational data

  , the mean shift clustering

Figure 8 :

 8 Figure 8: evidence lower bound of the model across variational inference iterations, mean and standard deviation over 10 runs.

Figure 9 :Figure 10 :

 910 Figure 9: Convergence time as a function of the number of instances for the variational inference and Gibbs sampling approaches.

Figure 12 :

 12 Figure 12: Graphical representation of the mixed-type infinite mixture model in plate notation.

  Figure 13: semi-supervised evidence lower bound across variational inference iterations.

Figure 14 :

 14 Figure 14: Accuracy of the model across variational inference iterations.

Figure 15 :

 15 Figure 15: Overview of the diagnosis system.

Figure 16 :

 16 Figure16: Architecture of the feature extractor followed by the softmax unit.

  of Algorithm 5) in the Algorithm 5 Deep Infinite Gaussian Mixture Model Input: X, y Initialize parameters # Step 1: Training of the feature extractor on labeled data compute θ * by minimizing (5.5) #Step 2: Feature computation of unlabeled data compute f n using the last hidden layer # Step 3: Clustering process on unlabeled data compute q * by minimizing (5.7) return ẑn = arg max k q * (z n = k) ∀n low dimensional feature manifold reduces considerably the complexity and execution time of the approach.

Figure 17

 17 Figure 17: 2-D representation of the feature manifold across epochs of training.

Figure 18 :

 18 Figure 18: (Left) 2-D representation in the feature space. (Right) 2-D representation of the raw data.

Figure 19 :

 19 Figure 19: Loss function on the test set across iterations of gradient descent. Accuracy on the test set across epochs.

Figure 21 :Figure 20 :

 2120 Figure 21: The clustering process of the infinite Gaussian mixture model on the latent feature manifold across iterations of variational inference.

Figure 22 :Figure 23 :

 2223 Figure 22: Training loss and log variance of the gradients for the different estimators for f (z) = ∑ d j=1 (z j -) 2 for d ∈ {1, 10, 100}.

Figure 25 :

 25 Figure 24: Bayesian Logistic Regression with Laplacian priors

Figure 27 :

 27 Figure 27: (Left) t-SNE plot of the second stochastic hidden layer on the MNIST test set for the semi-supervised (10% labels) version of the DP-DLM. (Right) t-SNE plot of the second stochastic hidden layer on the MNIST test set for the unsupervised version of the DP-DLM.

Figure 28 :

 28 Figure 28: Generated samples from the DP-DLM model for the unsupervised version on the MNIST dataset (left) and the semi-supervised version on the SVHN dataset (right).

Figure 29 :

 29 Figure 29: (Left) Evidence lower bound across learning iterations. (Right) Comparison of the accuracy between the Semi-supervised mixed type mixture and the Dirichlet process deep latent mixture.

Figure 30 :

 30 Figure 30: 2D representation of the last hidden layer for each class, and clusters discovered.

Figure 31 :

 31 Figure 31: Interpretation of fault cluster 40, corresponding to an account suspension of the customer on the network.

Figure 32 :

 32 Figure 32: Posterior distributions of the visible variables given the cluster assignment.

Figure 33 :

 33 Figure 33: An example of a decision rule leading to cluster 4.

Figure 34 :

 34 Figure 34: Clear decision leading to cluster 3.

Figure 37 :

 37 Figure 37: The training evidence lower bound on the MNIST training set (top) and the log variance of the gradient (bottom) over 5 runs. Comparison with the REBAR and RELAX estimators.

Figure 39 :

 39 Figure 39: Mean and standard deviation of the cumulative reward (left) and the log variance of the gradient for the actor network (right) over 5 iterations of the training process for the cart-pole task.

  des pannes dans les réseaux est d'une importance cruciale pour tout fournisseur de services Internet. Un système de diagnostic précis, capable d'identifier les pannes, peut améliorer l'expérience client, diminuer l'indisponibilité du réseau, et surtout réduire les coûts d'interventions sur le terrain pour le fournisseur d'accès Internet.

Table 1 :

 1 Description of each dataset.

Table 2 :

 2 Confusion Matrix of the Infinite Categorical Mixture Model on the synthetic dataset.

Table 3 :

 3 Comparison between our model and the baselines.

  For the variational parameters corresponding to the continuous dimensions, i.e. {L k , ν k , µ k , κ k }, the update equations are the following:

		1,k , γ	(t) 2,k , κ	(t)	(t) k	∀k (4.16)(4.19)
	Compute: Compute: µ k , L k ∀k (4.18)(4.20) (t) k,i,v ∀v, ∀k, ∀i (4.17)
	Compute: φ	(t) n,k ∀n, ∀k (4.21)
	if L	(t) ss changes by less than 10 -6 then
	break
	end if
	end for
	return φ, L, ν, µ, κ, γ,

k , ν

Table 4 :

 4 Comparison between our model and the baselines.

		Label
	Methods	Propagation TSVMs (Linear) TSVMs (RBF)	Ours
	Accuracy 0.95 ± 0.01	0.59 ± 0.02	0.62 ± 0.04 0.94 ± 0.01

Table 5 .

 5 This

	labelled sub-dataset corre-
	sponds to 64279 customers. Another
	part of the dataset (21962 customers)
	cannot be classified by the legacy diag-
	nosis system, and the corresponding customers thus fall into the "unknown
	faults" category.
	Network scope # variables Description
	Wide Area Network 652	GPON (OLT, ONT...)
	Home LAN 446	Gateway, home devices
	Internet Session 256	DHCP
	Services 204	VoIP, TV
	Offer 41	Customer offer profile
	Miscellaneous 225	-
	TOTAL 1824
	8297	(after One-Hot Encoding)

Table 6 :

 6 Overview of the network variables in the dataset.

Table 7 :

 7 The confusion matrix on the known classes.

								Predicted classes
						1	2	3	4		5	6	7	8
					1 1288 96	10	39		63	13	14	4
				True classes	2 47 3619 3 0 0 4 0 1 5 1 1	0 6024 5 8	0 3 18590 1	0 12 0 18258	45 5 304 621 14 2 7 6	7 4 83 0
					6	1	5	9	226		0	6580	74	11
					7	2	0	0	56		1	10 6710	3
					8	1	0	3	22		0	2	4	1364
		1.7
	Loss function	1.4 1.5 1.6							Accuracy	0.8 0.9
		1.3	0	10000 20000 30000 40000 50000 60000 iterations		0.7	0	20	40	epochs	60	80	100

Table 8 :

 8 Semi-supervised classification error (%) on the MNIST test set with 10 % labelisation. Comparison with(Nalisnick and Smyth, 2016)

  -114.14 ± 0.44 -114.55 ± 0.48 -110.87 ± 0.2 -110.70 ± 0.11 -110.59 ± 0.08 two layer SBN -101.33 ± 0.04 -101.09 ± 0.07 -99.74 ± 0.3 -100.44 ± 0.28 -100.66 ± 0.21 Bern. VAE -127.76 ± 0.84 -128.06 ± 2.66 -107.4 ± 1.47 -108.46 ± 0.37 -109.19 ± 1.31 Omniglot one layer SBN -123.66 ± 0.05 -123.82 ± 0.17 -113.53 ± 0.21 -114.34 ± 0.16 -114.37 ± 0.19 two layer SBN -117.81 ± 0.17 -117.89 ± 0.04 -102.05 ± 0.19 -102.16 ± 0.09 -102.29 ± 0.14 Bern. VAE -136.83 ± 0.31 -136.53 ± 0.32 -126.94 ± 0.81 -128.69 ± 0.34 -129.48 ± 0.38

	Dataset	Model	REBAR	RELAX	BSB (S=1)	BSB (S=5)	BSB (S=10)
		one layer SBN
	MNIST

Table 9 :

 9 Test likelihood for the Bernoulli stochastic backpropagation (BSB) estimator compared to the REBAR and RELAX estimators. We report the mean and standard deviation over 5 runs.

In order to optimize for the variational distribution, it suffices to optimize the parameters ψ (l+1) t

.

Given the outputs of the decoder defining the variational distribution, we sample h (l+1) n from q:

Then we forward it through the probabilistic model to construct the mean and covariance for the next layer as outputs of another neural network called the decoder: {m

Then the loss function minimization calibrates the parameters in order to accurately reconstruct the data x n . A summary of this process is given in Figure 35.

. -

In order to derive the mean-field equations of the Infinite Categorical Mixture model, we must compute the conditional expectations of equation ( 3.13), we start by specifying the joint distribution of the model as described in chapter 3:

In order to determine the optimal variational distributions and the update equations, we need to compute the conditional expectations of the mean-field equations for each hidden random variable.

. . Computing q * (z n ):

. . Computing q γ t (β t |x 1:N ):

The distribution q * (β t ) maximizing L, minimizes the kullback-leibler term. Hence:

By isolating the terms dependent on z n in L, we obtain: Finally, we evaluate the computational overhead of the categorical stochastic backpropagation estimator compared the Gumbel-softmax estimator. We compare the two estimators in terms of execution time of one epoch of training. The comparison is done using GPU implementations on a NVIDIA GeForce RTX 2080 Ti GPU, where the stochastic backpropagation rule of equation (6.18) is vectorized, thus leveraging the parallel batch treatment of the GPU. As shown in table 11, the Gumbel-softmax method is faster than stochastic backpropagation (S = 1) by a difference of about 3 seconds per training epoch. This is due to the forward passes performed to compute each of the terms in equation (6.18). The variance reduction and the increase in performance outweigh however the computational cost.

. In the reinforcement learning setting, we compare the stochastic backpropagation estimator of equation (C.4) to the famous actor-critic baseline [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF]. We compare the two methods on the classical cart-pole benchmark of openAI gym [START_REF] Brockman | OpenAI Gym[END_REF]. The Actor network in both cases is a two layer neural network with dimensions equal to [64,32] and hyperbolic tangent activation functions. The A2C method requires the critic network which serves as the baseline in the A2C gradient estimator. The Abstract : For any Internet service provider or network operator, it is crucial to quickly and efficiently diagnose the problems that occur on the network. The benefits of a good fault diagnosis system are mainly to minimize the costs of network and service operations and to enhance the customer's quality of experience. One major challenge for any diagnosis system concerns the discovery of new faults, that are unknown to the current version of the diagnosis system. The exploratory process for finding new faults can prove to be expensive and time consuming for internet service providers.

In this thesis, we explore an alternative approach based on learning methods, in order to build learning-based diagnosis systems. Our study explores Probabilistic Graphical Models that are capable of clustering patterns of faults in an unsupervised and a semi-supervised manner. We demonstrate the efficiency of our models on real use-cases of large scale data, extracted from Fiber-to-the Home (FTTH) services based on Gigabit-capable Passive Optical Networks.