Learning Linear Temporal Properties

We present two novel algorithms for learning formulas in Linear Temporal Logic (LTL) from examples. The first learning algorithm reduces the learning task to a series of satisfiability problems in propositional Boolean logic and produces a smallest LTL formula (in terms of the number of subformulas) that is consistent with the given data. Our second learning algorithm, on the other hand, combines the SAT-based learning algorithm with classical algorithms for learning decision trees. The result is a learning algorithm that scales to real-world scenarios with hundreds of examples, but can no longer guarantee to produce minimal consistent LTL formulas. We compare both learning algorithms and demonstrate their performance on a wide range of synthetic benchmarks. Additionally, we illustrate their usefulness on the task of understanding executions of a leader election protocol.

I. INTRODUCTION
Making sense of the observed behavior of complex systems is an important problem in practice. It arises, for instance, in debugging (especially in the context of distributed systems), reverse engineering (e.g., of malware and viruses), specification mining for formal verification, and modernization of legacy systems, to name but a few examples. However, understanding a system based on examples of its execution is clearly a challenging task that can quickly become overwhelming without proper tool support.
In this paper, we address this problem and develop learningbased techniques to help engineers understand the dynamic (i.e., temporal) behavior of complex systems. More precisely, we solve the problem of learning formulas in Linear Temporal Logic (LTL) [1], which are meant to distinguish between desirable and undesirable executions of a system (e.g., to explain the root-cause of a bug). The particular choice of LTL in this work is motivated by two observations: first, logical formulas often provide concise descriptions of the observed behavior and are relatively easy for humans to comprehend; second, LTL-together with Computational Tree Logic (CTL) [2]-is widely considered to be the de facto standard for specifying temporal properties and, hence, many engineers are familiar with its use.
The precise problem we are aiming at is the following: given a sample S consisting of two finite sets of positive and negative examples, learn an LTL formula ϕ that is consistent with S in the sense that all positive examples satisfy ϕ, whereas all negative examples violate ϕ. 1 To be as general and succinct as possible, we here consider examples to be infinite, ultimately periodic words (e.g., traces of a non-terminating system) and assume the standard syntax of LTL. However, our techniques can easily be adapted to the case of finite words and extend smoothly to arbitrary future-time temporal operators, such as "release", "weak until", and so on. We fix all necessary definitions and notations in Section II.
The main contribution of this work are two novel learning algorithms for LTL formulas from data, one based on SAT solving, the other on learning decision trees.
SAT-based learning algorithm: The idea of our first algorithm, presented in Section III, is to reduce the problem of learning an LTL formula to a series of satisfiability problems in propositional Boolean logic and to use highly-optimized SAT solvers to search for solutions. Inspired by ideas from bounded model checking [10], our learning algorithm produces a series of propositional formulas Φ S n for increasing values of n ∈ N \ {0} that depend on the sample S and have the following two properties: (1) Φ S n is satisfiable if and only if there exists an LTL formula of size n (i.e., with n subformulas) that classifies the examples correctly, and (2) a model of Φ S n contains sufficient information to construct such an LTL formula. By increasing the value of n until Φ S n becomes satisfiable, we obtain an effective algorithm that learns an LTL formula that is guaranteed to classify the examples correctly (given that the sample is non-contradictory).
By design, our SAT-based learning algorithm has three distinguished features, which we believe are essential in practice. First, our algorithm learns LTL formulas of minimal size (i.e., with the minimal number of subformulas). As we seek to learn formulas to be read by humans, the size of the learned formula is a crucial metric since larger formulas are generally harder to understand than smaller ones. Second, once an LTL formula has been learned, our algorithm can be queried for further, distinct formulas that are consistent with the sample. We believe that this feature is important in practice as it allows generating multiple explanations for the observed data. Third, our algorithm does not rely on an a priori given set of templates, which is in stark contrast to existing work on learning temporal properties (e.g., Bombara et al. [11]). To the best of our knowledge, our SAT-based algorithm is in fact the first learning algorithm that is not restricted to a fixed class of templates. However, restrictions to the shape of LTL formulas (e.g., to the popular GR(1)-fragment of LTL [12]) can easily be encoded if desired.
Learning algorithm based on decision trees: Our second learning algorithm, which we present in Section IV, trades in the guarantee of finding minimal solutions in order to attain better scalability. The key idea is to perform the learning in two phases. In the first phase, we run the SAT-based learning algorithm described above on various subsets of the examples. This results in a (small) number of LTL formulas, named "LTL primitives", that classify at least these subsets correctly.
In the second phase we use a standard learning algorithm for decision trees [13] to learn a Boolean combination of these LTL primitives that classifies the whole set of examples correctly, though it might not be minimal. Note, however, that we need to carefully choose the subsets of examples such that the resulting LTL primitives (a) separate all pairs of positive and negative examples and (b) are general enough to permit "small" decision trees. We have experimented with numerous strategies to select subsets, but in this paper we present only the two that performed best. A well known advantage of decision trees is that they are simple to comprehend due to their rule-based structure.
In Section V, we evaluate the performance of both learning algorithms on a wide range of synthetic benchmarks that reflect typical patterns of LTL formulas used in practice. Additionally, we illustrate their usefulness for understanding causes of inconsistencies in the leader election used by Zookeeper's atomic broadcast protocol [14].

Related Work
Learning of temporal properties from examples has recently attracted increasing interest, especially in the area of Signal Temporal Logic (STL) [16] and parametric STL [17]. Examples include the work by Asarin et al. [17], Kong et al. [18], [19], Vaidyanathan et al. [20], and Bartocci, Bortolussi, and Sanguinetti [21]. In contrast to our SAT-based learning algorithm, however, all of these techniques either rely on usergiven templates or can only learn formulas from very restricted syntactic fragments. Various techniques for mining LTL specifications [22], [23] and CTL specifications [24] exist as well, but these also rely on templates or restrict the class of formulas severely. To the best of our knowledge, our SAT-based algorithm is in fact the first that is capable of learning unrestricted LTL formulas without relying on usergiven templates. Nonetheless, expert knowledge in form of constraints on the syntax can easily be encoded if desired.
Our SAT-based learning algorithm is inspired by bounded model checking [10] and earlier work of the first author on learning (minimal) automata over finite words [7], [9]. However, since regular languages are strictly more expressive than LTL (the former being equivalent to monadic secondorder logic [25], while the latter being equivalent to fist-order logic [26]), automata learning techniques-including active learning algorithms [27], [28] that operate in Angluin's active learning framework [29]-are not immediately applicable. However, lifting the methods developed in this work to an active learning setup, without a detour via automata, is part of our plans for future work.
Using decision trees to learn Signal Temporal Logic (STL) formulas has been explored by Bombara et al. [11], whose main contribution is an adaptation of the classical impurity measure to account for STL formulas. However, this work still requires user-defined STL primitives to be provided, which serve as the features for the decision tree learning algorithm. By contrast, our technique uses the SAT-based learning algorithm to infer LTL primitives fully automatically.
Learning of logical formulas has also been studied in the context of probably approximately correct learning (PAC) [3]. Grohe and Ritzert [30], for instance, considered learning of first-order definable concepts over structures of small degree. Subsequently, Grohe, Löding, and Ritzert [31] studied the learning of hypotheses definable using monadic second order logic on strings. Due to the fundamental differences between PAC learning and the learning model considered here (one being approximate and the other being exact), their techniques cannot easily be applied.

II. PRELIMINARIES
In this section, we set up definitions and notations used throughout the paper.
Finite and Infinite Words: An alphabet Σ is a nonempty, finite set. The elements of this set are called symbols.
A finite word over an alphabet Σ is a sequence u = a 0 . . . a n of symbols a i ∈ Σ, i ∈ {0, . . . , n}. The empty sequence is called empty word and written as ε. The length of a finite word u is denoted by |u|, where |ε| = 0. Moreover, Σ * denotes the set of all finite words over the alphabet Σ, while Σ + = Σ * \{ε} is the set of all non-empty words.
An infinite word over Σ is an infinite sequence α = a 0 a 1 . . . of symbols a i ∈ Σ, i ∈ N. We denote the i-th symbol of an infinite word α by α(i) and the infinite suffix starting at position j by α[j, ∞). Given u ∈ Σ + , the infinite word u ω = uu . . . ∈ Σ ω is the infinite repetition of u. An infinite word α is called ultimately periodic if it is of the form α = uv ω for a u ∈ Σ * and v ∈ Σ + . Finally, Σ ω denotes the set of all infinite words over the alphabet Σ.
Propositional Boolean Logic: Let Var be a set of propositional variables, which take Boolean values from B = {0, 1} (0 representing false and 1 representing true). Formulas in propositional (Boolean) logic-which we denote by capital Greek letters-are inductively constructed as follows: • each x ∈ Var is a propositional formula; and • if Ψ and Φ are propositional formulas, so are ¬Ψ and Ψ ∨ Φ. Moreover, we add syntactic sugar and allow the formulas true, false, Ψ ∧ Φ, Ψ ⇒ Φ, and Ψ ⇔ Φ, which are defined as usual.
A propositional valuation is a mapping v : Var → B, which maps propositional variables to Boolean values. The semantics of propositional logic is given by a satisfaction relation |= that is inductively defined as follows: v |= x if and only if v(x) = 1, v |= ¬Ψ if and only if v |= Ψ, and v |= Ψ ∨ Φ if and only if v |= Ψ or v |= Φ. In the case that v |= Φ, we say that v satisfies Φ and call it a model of Φ. A propositional formula Φ is satisfiable if there exists a model v of Φ. The size of a formula is the number of its subformulas (as defined in the usual way).
The satisfiability problem of propositional logic is the problem to decide whether a given formula is satisfiable. Although this problem is well-known to be NP-complete [32], modern SAT solvers implement optimized decision procedures that can check satisfiability of formulas with millions of variables [33]. Moreover, SAT solvers also return a model if the input-formula is satisfiable.
Linear Temporal Logic: Linear Temporal Logic (LTL) [1] is an extension of propositional Boolean logic with modalities that allow expressing temporal properties. Starting with a finite, nonempty set P of atomic propositions, formulas in LTLwhich we denote by small Greek letters-are inductively defined as follows: • each atomic proposition p ∈ P is an LTL formula; • if ψ and ϕ are LTL formulas, so are ¬ψ, ψ ∨ ϕ, X ψ ("next"), and ψ U ϕ ("until"). Again, we add syntactic sugar and allow the formulas true := p ∨ ¬p for some p ∈ P, false := ¬true, as well as ψ ∧ ϕ and ψ → ϕ, which are defined as usual. Moreover, we allow the additional temporal formulas F ψ := true U ψ ("finally") and G ψ := ¬ F ¬ψ ("globally"). The size of an LTL formula ϕ, which we denote by |ϕ|, is the number of its subformulas. Finally, let C = {∧, ∨, ¬, →, F, G, U, X} be the set of LTL operators.
LTL formulas are interpreted over infinite words α ∈ (2 P ) ω , though there exist various semantics for LTL over finite words and our techniques smoothly extend to these situations. For the sake of a simpler presentation, we define the semantics of LTL in a slightly non-standard way by means of a valuation function V . This functions maps pairs of LTL formulas and infinite words to Boolean values and is inductively defined as follows: Our SAT-Based learning algorithm relies on a canonical syntactic representation of LTL formulas, which we call syntax DAGs. A syntax DAG is essentially a syntax tree (i.e., the unique tree labeled with atomic propositions as well as Boolean and temporal operators that is derived from the inductive definition of an LTL formula) in which common subformulas are shared. This sharing turns the syntax tree into a directed, acyclic graph (DAG), whose number of nodes coincides with the number of subformulas of the represented LTL formula. As an example, Figure 1b (on Page 4) depicts the (unique) syntax DAG of the formula (p U G q) ∨ (F G q), in which the subformula G q is shared; the corresponding syntax tree is depicted in Figure 1a. Note that syntactically distinct formulas have different (i.e., non-isomorphic) syntax DAGs.
Samples and Consistency: Throughout this paper, we assume that the data we learn from is given as two (potentially empty) finite, disjoint sets P, N ⊂ (2 P ) ω of ultimately periodic words. The words in P are interpreted as positive examples, while the words in N are interpreted as negative examples. We call the pair S = (P, N ) a sample. Since we want to work with the ultimately periodic words in a sample algorithmically, we assume that they are stored as pairs (u, v) of finite words u ∈ (2 P ) * and v ∈ (2 P ) + , which can be accessed individually. To measure the complexity of a sample, we define its size to be |S| = uv ω ∈P ∪N |u| + |v|.
Given an LTL formula ϕ and a sample S = (P, N ), both over a set P of atomic propositions, we call ϕ consistent with S if V (ϕ, uv ω ) = 1 for each uv ω ∈ P (i.e., all positive examples satisfy ϕ) and V (ϕ, uv ω ) = 0 for each uv ω ∈ N (i.e., all negative examples do not satisfy ϕ); in this case, we also say that ϕ separates P and N . We call ϕ minimally consistent with S if ϕ is consistent with S and no consistent LTL formula of smaller size exists.

III. A SAT-BASED LEARNING ALGORITHM
The fundamental task we solve in this section is: "given a sample S, compute an LTL formula of minimal size that is consistent with S". We call this task passive learning of LTL formulas-as opposed to active learning [29] where the learning algorithm is permitted to actively query for additional data. Note that this problem can have more than one solution as there can be multiple, non-equivalent LTL formulas that are minimally consistent with a given sample.
Before we explain our learning algorithm in detail, let us briefly comment on the minimality requirement in the definition above. On the one hand, we observe that the problem becomes simple if no restriction on the size is imposed: for α ∈ P and β ∈ N , construct a formula ϕ α,β with V (ϕ α,β , α) = 1 and V (ϕ α,β , β) = 0 that describes the first symbol where α and β differ using a sequence of X-operators and an appropriate propositional formula; then, α∈P β∈N ϕ α,β is consistent with S since we assume P and N to be disjoint. However, simply characterizing all differences between positive and negative examples is clearly overfitting the sample and, hence, arguably of little help in practice. On the other hand, we believe that small formulas are easier for humans to comprehend than large ones, which justifies spending effort on learning a smallest formula. However, we do not impose any preference amongst minimal consistent formulas (which is an interesting topic for future work).
Let us now turn to describing our learning algorithm. Its underlying idea is to reduce the construction of a minimally consistent LTL formula to a satisfiability problem in propositional logic and use a highly-optimized SAT solver to search for solutions. More precisely, given a sample S and a natural number n ∈ N \ {0}, we construct a propositional formula Φ S n of size polynomial in n and |S| that has the following two properties: 1) Φ S n is satisfiable if and only if there exists an LTL formula of size n that is consistent with S; and 2) if v is a model of Φ S n , then v contains sufficient information to construct an LTL formula ψ v of size n that is consistent with S. By increasing the value of n by one and extracting an LTL formula ψ v from a model v of Φ S n as soon as it becomes satisfiable (indeed, any model is sufficient), we obtain an effective algorithm that learns an LTL formula of minimal size that is consistent with S. This idea is shown in pseudo code as Algorithm 1. In fact, the existence of a trivial solution for the passive LTL learning task (as sketched at the beginning of this section) shows that Algorithm 1 is guaranteed to terminate, and the size of this solution provides an upper bound on the value of n. Construct and solve Φ S n ; 5 until Φ S n is satisfiable, say with model v; 6 Construct and return ψ v ; The key idea of the formula Φ S n is to encode the syntax DAG of an (unknown) LTL formula ϕ ⋆ with n subformulas and then constrain the variables of Φ S n such that ϕ ⋆ is consistent with the sample S. To simplify our encoding, we assign to each node of this syntax DAG a unique identifier i ∈ {1, . . . , n} such that (a) the identifier of the root is n and (b) if the identifier of an inner node is i, then the identifiers of its children are less than i. Note that such a numbering scheme is not unique for a given syntax DAG, but it entails that the root always has identifier n and the node with identifier 1 is always labeled with an atomic proposition. We refer the reader to Figures 1b and 1c for an example.
We encode a syntax DAG using three types of propositional variables: . . , n} and λ ∈ P ∪ C; • l i,j where i ∈ {2, . . . , n} and j ∈ {1, . . . , i − 1}; and • r i,j where i ∈ {2, . . . , n} and j ∈ {1, . . . , i − 1}. Intuitively, the variables x i,λ encode a labeling of the syntax DAG in the sense that if a variable x i,λ is set to true, then node i is labeled with λ (recall that each node is labeled with either an atomic proposition from P or an operator from C). The variables l i,j and r i,j , on the other hand, encode the structure of the syntax DAG (i.e., the left and/or right child of inner nodes): if variable l i,j (r i,j ) is set to true, then j is the identifier of the left (right) child of node i. By convention, we ignore the variables r i,j if node i of the syntax DAG is labeled with an unary operator; similarly, we ignore both l i,j and r i,j if node i is labeled with an atomic proposition. Note that in the case of l i,j and r i,j , the identifier i ranges from 2 to n because node 1 is always labeled with an atomic proposition and, hence, cannot have children. Moreover, j ranges from 1 to i − 1 to reflect the fact that identifier of children have to be smaller than the identifier of the current node.
To enforce that the variables x i,λ , l i,j , and r i,j in fact encode a syntax DAG, we impose the constraints listed in Table I. Formula (1) ensures that each node is labeled with exactly one label. Similarly, Formulas (2) and (3) enforce that each node (except for node 1) has exactly one left and exactly one right child (although we ignore certain children if the node represents an unary operator or an atomic predicate). Finally, Formula (4) makes sure that node 1 is labeled with an atomic proposition.
Let Φ DAG n now be the conjunction of Formulas (1) to (4). Then, one can construct a syntax DAG from a model v of Φ DAG n in a straightforward manner: simply label node i with the unique label λ such that v(x i,λ ) = 1, designate node n as the root, and arrange the nodes of the DAG as uniquely described by v(l i,j ) and v(r i,j ). Moreover, we can easily derive an LTL formula from this syntax DAG, which we denote by ψ v . Note, however, that ψ v is not yet related to the sample S and, thus, might or might not be consistent with it.
To enforce that ψ v is indeed consistent with S, we now constrain the variables x i,λ , l i,j , and r i,j further. More precisely, we add for each ultimately periodic word uv ω in S a propositional formula Φ u,v n that tracks the valuation of the LTL formula encoded by Φ DAG n (and all its subformulas) on uv ω . The observation that enables us to do this is the following.
Intuitively, Observation 1 states that the suffixes of a word uv ω eventually repeat periodically. As a consequence, the valuation of an LTL formula on a word uv ω can be determined based only on the finite prefix uv (recall that the semantics of temporal operators only depend on the suffixes of a word). To illustrate this claim, consider the LTL formula X ϕ and assume that we want to determine the valuation V (X ϕ, uv ω [|uv| − 1, ∞)) (i.e., X ϕ is evaluated at the end of the prefix uv). Then, Observation 1 permits us to compute this valuation based on V (ϕ, uv ω [|u|, ∞)), as opposed to the original semantics of the X-operator, which recurs to V (ϕ, uv ω [|uv|, ∞)) (i.e., the valuation at the next position). Note that similar, though more involved ideas can be applied to all other temporal operators.
Each formula Φ u,v n is built over an auxiliary set of propositional variables y u,v i,t where i ∈ {1, . . . , n} is a node in the syntax DAG and t ∈ {0, . . . , |uv| − 1} is a position in the finite word uv. The meaning of these variables is that the value of y u,v i,t corresponds to the valuation V (ϕ i , uv ω [t, ∞)) of the LTL subformula ϕ i that is rooted at node i. Note that the set of variables for two distinct words from the sample must be disjoint.
To obtain this desired meaning of the variables y u,v i,t , we impose the constraints listed in Table II, which are inspired by bounded model checking [10]. Formula (5) implements the LTL semantics of atomic propositions and ensures that if node i is labeled with p ∈ P, then y u,v i,t is set to 1 if and only if p ∈ uv(t). Next, Formulas (6) and (7) implement the semantics of negation and disjunction, respectively: if node i is labeled with ¬ and node j is its left child, then y u,v i,t is the negation of y u,v j,t ; on the other hand, if node i is labeled with ∨, node j is its left child, and node j ′ is its right child, then y u,v i,t is the disjunction of y u,v j,t and y u,v j ′ ,t . Moreover, Formula (8) implements the semantics of the X-operator, following the idea of "returning to the beginning of the periodic part v" as sketched above. Finally, Formula (9) implements the semantics of the U-operator. More precisely, the first conjunction in the consequent covers the positions t ∈ {0, . . . , |u| − 1} in the initial part u, while the second conjunct covers the positions t ∈ {|u|, . . . , |uv| − 1} in the periodic part v. Thereby, the second conjunct relies on an auxiliary set t u,v t ′ defined by which contains all positions in v "between t and t ′ ". To avoid cluttering this section too much, we have omitted the description of the missing operators ∧, →, F, G and the 1<i≤n 1≤j<i constants true and false, which are implemented analogously. Moreover, our SAT encoding is extensible, and additional LTL operators such as weak until or weak and strong release can easily be added. For each uv ω ∈ P ∪ N , let Φ u,v n now be the conjunction of Formulas (5) to (9). Then, we define Note that the subformula Φ u,v n ∧y u,v n,0 makes sure that uv ω ∈ P satisfies the prospective LTL formula (more concretely, uv ω starting from position 0 satisfies the LTL formula at the root of the syntax DAG), while Φ u,v n ∧¬y u,v n,0 ensures that uv ω ∈ N does not satisfy it.
To prove the correctness of our learning algorithm, we first establish that the formula Φ S n has in fact the desired properties. Lemma 1: Let S = (P, N ) be a sample, n ∈ N\{0}, and Φ S n the propositional formula defined above. Then, the following holds: 1) If an LTL formula of size n that is consistent with S exists, then the propositional formula Φ S n is satisfiable. 2) If v |= Φ S n , then ψ v is an LTL formula of size n that is consistent with S. Termination and correctness of Algorithm 1 then follow from Lemma 1.
Theorem 1: Given a sample S, Algorithm 1 terminates eventually and outputs an LTL formula of minimal size that is consistent with S. Proof: Since there exists a consistent LTL formula for every non-contradictory sample, Part 1 of Lemma 1 guarantees that Algorithm 1 terminates. Moreover, Part 2 ensures that the output is indeed an LTL formula that is consistent with S. Since n is increased by one in every iteration of the loop until Φ S n becomes satisfiable, the output of Algorithm 1 is a consistent LTL formula of minimal size.
It is important to emphasize that the size of Φ S n and, hence, the performance of Algorithm 1 depends on the size of a sample S = (P, N ), as summarized next.
Remark 1: The formula Φ S n ranges over O(n 2 + n|S|) variables and is of size O(n 2 + n 3 uv ω ∈P ∪N |uv| 3 ). Finally. we conclude this section with a remark on incorporating expert knowledge into the learning process.
Remark 2: By adding constraints to the variables x i,λ , l i,j , and r i,j , one can easily incorporate expert knowledge (e.g., syntactic templates) into the learning process.

IV. A DECISION TREE BASED LEARNING ALGORITHM
The SAT-based algorithm described in Section III is an elegant, out-of-the-box way to discover minimal LTL formulas describing a sample. Even though it scales well beyond toy examples, its running time seems too prohibitive for real-world examples (as discussed in Section V). That is why we now present a learning algorithm based on a combination of SAT solving and decision tree learning.
Our second algorithm proceeds in two phases, outlined in Algorithm 2. In the first phase, we run Algorithm 1 on small subsets of P and N . This is repeated until we obtain a set Π of LTL formulas (we call them LTL primitives) that separate all pairs of words from P and N . In the second phase, formulas from Π are used as features for a standard decision tree learning algorithm [13]. The resulting decision tree is a Boolean combination of LTL formulas ϕ i ∈ Π that is consistent with the sample.

Algorithm 2: Learning algorithm based on decision trees
Input: a sample S 1 Run Algorithm 1 on small subsets of P and N to construct a set Π = {ϕ 1 , . . . , ϕ n } of LTL formulas such that for each pair u 1 v ω 1 ∈ P and u 2 v ω 2 ∈ N there exists a ϕ i ∈ Π with V (ϕ i , u 1 v ω 1 ) = 1 and V (ϕ i , u 2 v ω 2 ) = 0; 2 Learn a decision tree t with LTL primitives from Π as features and return the resulting Boolean combination ψ t of LTL primitives (which is consistent with S); Note that this relaxes the problem addressed in Section III: we can no longer guarantee finding a formula of minimal size. However, decision trees are among the structures that are the easiest to interpret by end-users. That makes them suitable for our use-case, and the minimality of formulas is replaced by structural simplicity of decision trees.
Learning Decision Trees: We assume familiarity with decision tree learning and refer the reader to a standard textbook for further details [5]. As illustrated in Figure 3, the decision trees we seek to learn are tree-shaped structures whose inner nodes are labeled with LTL formulas from Π and whose leaves are labeled with either true or false. The LTL formula represented by such a tree t is given by ψ t := ρ∈P ϕ∈ρ ϕ where P is the set of all paths from the root to a leaf labeled with true and ϕ ∈ ρ denotes that ϕ occurs on ρ (negated if the path follows a dashed edge).
To learn a decision tree over LTL primitives, we perform a preprocessing step and modify the sample as follows. For each word uv ω ∈ P ∪N , we use the LTL primitives as features and create a Boolean vector of size |Π| with the i-th entry set to V (ϕ i , uv ω ); this vector is then labeled with true if uv ω ∈ P or with false if uv ω ∈ N . In the second step, we apply a standard learning algorithm for decision trees to this modified sample (we used Gini impurity [34] as split heuristic in our experiments). Since we are interested in a tree that classifies our sample correctly, we disable heuristics such as pruning.
Obtaining LTL Primitives: Meaningful features are essential for a successful classification using decision trees. In our algorithm, features are generated from the set of LTL primitives Π. We used two different strategies, called Strategy α and Strategy β, for obtaining Π.
Strategy α iteratively chooses subsets P ′ ⊂ P and N ′ ⊂ N of size k according to probability distributions prob P and prob N on P and N , respectively. After a formula ϕ separating P ′ and N ′ is found using Algorithm 1 and added to Π, prob P and prob N are updated to increase the likelihood of any word that is not yet classified correctly by any of the ϕ ∈ Π to be selected. This process is repeated until all pairs of positive and negative examples are separated by some LTL primitive or restarted after a user-given number of iterations. Although this strategy is, in general, not guaranteed to terminate due to its probabilistic nature, it always did in our experiments.
Strategy β computes LTL primitives in a more aggressive way. Starting with the set S = P × N , it uniformly at random selects k pairs from S and uses Algorithm 1 to compute an LTL primitive ϕ that separates those pairs. Then, it removes all pairs separated by ϕ from S and repeats the process until S becomes empty (i.e., all pairs of examples are separated).
We refer to the extended version of this paper [35] for a detailed explanation of both strategies.
Correctness: The correctness of Algorithm 2 is formalized below.
Theorem 2: Given a sample S, Algorithm 2 learns a (not necessarily minimal) formula ψ t that is consistent with S. Theorem 2 follows from the fact that Step 1 of Algorithm 2 constructs a set of LTL primitives that allows separating any pair of positive and negative examples. Once such a set is constructed, any decision tree learner produces a decision tree t that is guaranteed to classify the examples correctly. The resulting LTL formula ψ t , hence, is consistent with S.

V. EVALUATION
In this section, we answer questions that arise naturally: how performant is Algorithm 1 and what is the performance gain of Algorithm 2. Furthermore, what is the complexity of the learned decision trees in terms of the number of decision nodes, and, finally, how do different parameters influence the performance of Algorithm 2. After answering these questions with experiments performed on synthetic data, we demonstrate the usefulness of our algorithms for understanding executions of a leader-election algorithm. We implemented both learning algorithms in a Python tool 2 using Microsoft Z3 [36]. All experiments were conducted on Debian machines with Intel Xeon E7-8857 CPUs at 3 GHz, using up to 5 GB of RAM.
Performance on Synthetic Data: To simulate real-world use-cases, we generated samples based on common LTL patterns [37], which are shown in Table III. Starting from a pattern formula ψ, we generated sets of random words and separated them into P and N depending on whether they are a model of ψ or not. Thereby, we fixed |u| + |v| = 10 for all words in the sample and added noise in form of one additional atomic proposition that is not constrained by the pattern formula. The size of the generated samples ranges between 50 and 5000. In total, we generated 192 samples. Figure 2 compares the running times of Algorithm 1 and Algorithm 2 (using Strategy α and k = 3) on samples of varying sizes. (So as not to clutter the presentation too much, we selected four LTL patterns that showed a typical behavior of our learning algorithms. The complete results are available in the technical report [35].) Overall, Algorithm 1 produces minimal formulas consistent with a sample. It does so even for samples of considerable size, but if the sample size grows beyond 2000 (varies over samples), the SAT-based learner (Algorithm 1) frequently times out. When Algorithm 2 (using decision tree learning) is applied to these samples-as shown on the right-hand-side of Figure 2-none of the computations timed out and the running times significantly improved.
What kind of trees does Algorithm 2 produce? An example output of the algorithm is shown in Figure 3. Moreover, as Table IV illustrates, Algorithm 2 learns small trees, often with less than five inner nodes. Upon closer inspection, we noticed that it often happens that one of the LTL primitives was the specified formula itself. This suggests that small subsets already characterize our samples completely.
To be able to compare decision trees to the formulas learned by Algorithm 2, we measure the size of a tree t in terms of the size of the formula ψ t this tree encodes. In our experiments, the formulas learned by Algorithm 2 were on average 1.41 times larger than those learned by Algorithm 1. However, there are outlier trees that are four times bigger than the one learned by Algorithm 1. Nonetheless, about 70 % are of the same size. Even for the outliers, as emphasized previously, the 2 Our tool is publicly available at https://github.com/gergia/samples2LTL.
false If a node evaluates to true, the outgoing full edge should be followed (and the dashed edge otherwise). The tree on the left represents the formula readability does not degrade completely because the rule-based structure of decision trees is known to be easily understandable by humans. Note that the runtime and size of decision trees depends on the parameters of Algorithm 2, which we discuss next.
Tuning the Decision Tree-Based Algorithm: As described in Section IV, Algorithm 2 can be tuned by various parameters (sampling strategy for obtaining LTL primitives, size of sample subsets, probability increase rate, and number of repetitions inside a single sampling). In this subsection, we explore how those parameters affect the performance of the algorithm.  Table IV shows the performance of Algorithm 2 for different parameters, averaged over all 192 benchmarks. As the table indicates, the less aggressive method of separating sets, Strategy α, performs better. It seems that if the subset sizes are increased, or Strategy β is used, the sampled subsets already describe the specified formula completely. Finally, we chose Strategy α and k = 3 to be our default parameters. Varying the probability decrease rate and the number of repetitions inside a single sampling did not influence the performance much.
Explaining Executions of a Leader Election Protocol: A number of methods exist for finding errors or reproducing certain behavior in distributed systems through systematic testing [38], [39]. However, finding an execution and a corresponding schedule is only a first step towards understanding an issue. In the following, we demonstrate how to apply our technique in order to obtain a minimal LTL description of a specific inconsistency in a leader election protocol.
The leader election protocol we consider is the Fast Leader Election algorithm [14], [40] used by Apache Zookeeper. In this protocol, every node has a unique ID and initially tries to become the leader. To this end, every node sends messages to all other nodes proclaiming its leadership. Upon receiving a message by an aspirant leader with a higher ID, a node gives up its claim and acknowledges its support for the aspirant. If a node learns that an aspirant node has a support of a majority of all nodes, it commits (after waiting for a constant time for new messages) to the aspirant as the leader. Once committed, the node never again changes its decision and informs any other node of its commitment (one example is the message depicted by the dotted arrow in Figure 5). If a node has not committed and learns about another node that has committed, it commits to the same leader.  Figure 4 shows an example of a successful leader election with three nodes in an UML-style message sequence chart. The messages exchanged between nodes are proposing the leader i (P i ) and node j acknowledging the claim of a leader (A j ). The arrows indicate exchanged messages and imply a precedence of events. Note that not all messages are shown in the figures, but only the ones important for understanding the protocol.
In Figure 4 all the nodes have committed to the same leader. On the other hand, Figure 5 shows a schedule that ends up in an inconsistent state where nodes committed to different leaders. This schedule was discovered by the PCTCP algorithm [41], which systematically explores the space of possible executions of distributed algorithms. The situation in Figure 5 is caused by the asynchronous communication: for performance reasons, nodes commit as quickly as possible and then discard any messages, which otherwise would have changed their commitment (indicated as a dashed line in Figure 5). Note, however, that this is not a bug in Zookeeper's broadcast algorithm, as a leader without a quorum will not be allowed to perform any action in the later phase.
To better understand how this inconsistent state arises, our goal is to generate an LTL formula that describes the difference between the schedules in Figures 4 and 5. To this end, we constructed a sample by generating 20 linearizations of the schedule from Figure 4 and 20 linearizations of the schedule from Figure 5. Since we seek an explanation for the inconsistent behavior, the former (with consistent outcomes) correspond to negative examples (set N ), and the latter (with inconsistent outcomes) correspond to positive examples (set P ). The set of atomic propositions used to construct the examples contains twelve elements: recv(i, j) for i, j ∈ {1, 2, 3} (meaning that node j received a message from node i) and comm(i) for i ∈ {1, 2, 3} (meaning that node i committed to a leader). 3 Finally, we ran Algorithm 1 on this sample. The result was the formula ¬recv(2, 1) U comm (1). Intuitively, node 1 did not receive a message from node 2 before it committed to a leader. That is exactly the difference between the schedules in Figures 4 and 5. Also, it hints at a specific reason for the inconsistency in Figure 5, thus potentially helping the engineers improve the system. Note, however, that this experiment still required a significant amount of manual effort. In order to apply the technique in practice, more automation is needed.
Summary: Algorithm 2 significantly improves upon the performance of Algorithm 1, though with a small increase in the size of the formula. The original motivation of getting readable explanations for the behavior of a system is preserved due to the fact that decision-trees are easy to comprehend. Algorithm 2 works the best using Strategy α and subsets of size k = 3. Finally, our techniques are able to give interesting insight into real-world systems.

VI. CONCLUSION
We have presented two novel algorithms for learning LTL formulas from examples. Our first algorithm is based on SAT solving, while the second algorithm extends the first with techniques for learning decision trees. We have shown that both algorithms are able to learn LTL formulas for a comprehensive set of benchmarks that we have derived from common LTL patterns. Moreover, we have demonstrated how our methods can help understand distributed algorithms.
Interesting directions of future work include the integration of LTL past-time operators, lifting our techniques to an active learning setup [29], as well as the development of similar learning algorithms for CTL. Furthermore, we plan to investigate the use of maximum-margin classifiers, such as support vector machines. To this end, one needs to develop a notion of distance between temporal formulas and words, which is clearly of independent, theoretical interest as well.


An approach for automating the design of convolutional neural networks

Image recognition is an independent field of the computer science nowadays. Image classification is one of its main domains, in which investigated objects can be represented by an image or a video stream. The objective of the image classification is correct assigning of objects to corresponding classes, and there exist many effective approaches for solving this problem. One of the most popular approaches is artificial neural networks, which are a method from the field of machine learning. Despite the fact that neural networks cover a wide range of machine learning problems, they are also able to solve the problem of the image classification. However, there is one more specific approach for neural networks-based images classification that applies the deep learning conception. The best-known deep learning algorithm is called the convolutional neural network (CNN). The CNN uses a principle of using the same parts of a neural network to manipulate with different local parts of an input image. As well as the standard neural network architecture, the convolutional neural network should be fine-tuned for solving a certain problem. Because of the CNN’s depth and complexity, the tuning process usually is complex and needs huge computational efforts. In this study, we have proposed an approach for creating ensembles of previously trained convolutional neural networks. The approach allows to increase the performance of the image classification. The results of experiments for image classification problems are presented and discussed. The experiments show that the proposed approach is able to outperform the standard perceptron and single convolutional neural network.

Introduction
Image recognition is one of the most complicated problem for machine learning techniques. For many real-world applications a solution to the problem cannot be described as a set of instructions for a machine learning algorithm. Instead of that it is needed to design a model which will be able to estimate the instances (input data) by itself. For many machine learning techniques, both the structure of the model and its parameters are needed to be fine-tuned. In a case of simple models, the tuning process can be done manually by a human-expert, but for complex image recognition problems simple models do not perform well, thus an automated designing of more complicated structures is preferable.
Starting with the LeNet-5 [1], convolutional neural networks (CNNs) usually had the standard structure included lots of layers ( Figure 1). This structure shows perfect results working with MNIST, CIFAR and ImageNet benchmark problems [2,3]. For many real-world problems, there is no stand-alone algorithm which can achieve a high performance for the given problem. In this case one of the most powerful technique for obtaining the best possible results for the problem is creating an ensemble model. Ensemble models consist of pooling together the predictions of different algorithms, to produce better predictions [4]. There are some popular methods of designing ensemble models, such as the simple voting and the weighted voting. In this study, we will use the simple voting approach in order to take into consideration all predictions from all CNN models. The output of the ensemble will be a label, which polled the majority of the votes.
In this study we will automate the process of CNNs designing by defining some special filters which are used as kernels in the CNN structure. Moreover, we will use an ensemble of previously trained CNN with the identic subsampling layers. The results of experiments are compared with fully connected NN with the same structure as in the fully-connected layers in CNNs.
The rest of the paper is organized as follows. Section 2 describes related work. Section 3 describes the proposed approach. In Section 4 the results of numerical experiments are discussed. In the Conclusion the results and further research are discussed.

Related work
The term "neural network" has its origins in attempts to find mathematical representations of information processing in biological systems (McCulloch and Pitts, 1943;Widrow and Hoff, 1960;Rosenblatt, 1962;Rumelhart et al., 1986). Indeed, it has been used very broadly to cover a wide range of different models [5].
Before deep learning appeared, David Hubel and Torsten Wiesel performed one of the most important research about striate cortex [6]. They explored parts of the visual cortex which are responsible for detecting edges. They found that some neurons activated when there were vertical lines, other neurons activated when there were horizontal lines. There were also neurons which were responsible for detecting other particular angles. A few later scientists explored that visual cortex is organized as a layered structure. Exactly these ideas affected on designing convolutional neural nets.
Deep Learning and its subfield CNNs appeared not so long ago on a scale of machine learning. Now we can claim that the first CNN was called Neocognitron and it was designed by Fukushima K. in 1979Fukushima K. in -1980. Today there is designed a wide range of CNNs architectures. And despite the fact that they show excellent results, the majority of them uses the same ideas. On the other hand, some other approaches become more popular. And one of these approaches is ensemble modeling. Ensemble models usually consist only of CNNs. There are some researches that proves its higher efficiency on various image recognition problems, such as MNIST [9] or FMOW [10].

Proposed Approach
Each CNN has three types of layers: convolutional layers, subsampling layers and fully-connected layers. First two layers differ this algorithm from the others. Convolutional layer represents the list of kernels, where each kernel is a special filter applied to the input image. These kernels aimed to find some graphical patterns in original image and in the images, which were processed by the previous layers. The output of this layer is represented by the feature maps. Subsampling layer created to reduce the dimensionality of the problem and the feature maps. This layer pools closely-spaced pixels. Let's discuss each layer in details.

Convolutional layer.
As it was mentioned, one of the most important part of convolutional layer is kernel. Kernel represents the process of filter applying to the input image. In practice, various types of kernels are used. What is more, some training techniques are implemented to tune the kernels. Usually the adapted backpropagation algorithm is used. Each filter follows its own purpose. However, the majority of them is used to increase the quality of pixel image. In this paper the following filters were considered:  Gradient filters (directions: north, south, west, and east). Belong to difference linear filters. These filters are also called edge detection filters.  Laplacians (kernel = 3x3 pixels, 5x5 pixels). These filters commonly used to increase the sharpness and applied to the images which were previously blurred in case to reduce the noise.  Line detection filters (incrementation of an angle degree = 45°).  Sobel operator (lines: horizontal, vertical). Discrete differentiation operator used to create an image with emphasizing edges.  Sharpening filters (2 levels of sharpness). These filters are used to highlight fine detail in an image or to enhance detail that has been blurred.  Smoothing filter. These filters are used to clean noised images. Kernel is just a special kind linear operator applied to the input data. Let's define the following variables: -feature map of -th layer, kernel size , -weight matrix, where = . The result of convolution with kernel of size can be represented as follows: , where -the result of convolution at the -th layer, and -the output of the whole the previous layer. In other words, to get the component of the next layer, linear operator is applied to the previous layer, so scalar multiplication of pixel values and kernel vector is performed. The example of filter applying is shown in Figure 2.

Subsampling layer.
In recognition problems it is possible to pool feature pixels in some cases. First of all, the information about some features existing is more valuable knowledge than the concrete location of feature border pixels. What is more, the dimensionality of some problems is too high and it is needed to reduce the dimensionality. The example of pooling operator is shown in Figure 3. There are some functions can be used, such as finding maximum, minimum, median, etc. And one more parameter which needed to be defined is step. In many papers maximum function is used. For example, in object detection [11] and human pose detection [12] problems this operator shows better results. Moreover, some researches explored that max-pooling is more efficient function [13].
After applying listed layers and training fully-connected layers, we need to design ensemble models. For choosing the general prediction simple voting is used. There is not a concrete strategy for uniting different CNNs, because each filter determines unique features, and it is hard to predict which combination of extracted features works better. What it more, as we shall see throughout the paper, the CNNs combinations are different in all problems.

Experimental Settings and Results
To estimate the proposed approach performance, we have used 3 image recognition problems.

CIFAR -10.
It is an object recognition problem. Number of examples = 60000. Dimensionality of the problem: 1024 (each image has a resolution = 32x32 pixels). Number of classes = 10. Each class has 6000 images, so it is a balanced dataset. The classes are following: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, truck. The classes are completely mutually exclusive and there is no overlap between automobiles and trucks, for example. Color model: RGB.

The ORL Database of Faces.
It is a face recognition problem. This dataset contains a set of face images taken between April 1992 and April 1994 at the lab. There are ten different images of each of 40 distinct subjects. For some subjects, the images were taken at different times, varying the lighting, facial expressions (open / closed eyes, smiling / not smiling) and facial details (glasses / no glasses). All the images were taken against a dark homogeneous background with the subjects in an upright, frontal position (with tolerance for some side movement). Dimensionality of the problem: 10304 (each image has a resolution = 92x112 pixels). Number of classes = 40. Color model: grayscale.

Experimental settings.
Firstly, all the filters listed in Section 3 were used to create and train different CNNs (each CNN had a one specific filter). Results of some filters applying are shown in Figures 5 -7. After that, the ensemble models were created, based on previously trained CNNs. Highest possible number of CNNs in one ensemble model equals 3.

Experimental Results.
For each experiment 10 independent runs were computed. Performance shows mean value of these runs. Performance of fully-connected NN for all considered problems is presented in Table 1 Performance of the best single CNNs using one filter is presented in Tables 2-4.   Performance of the best ensemble models is presented in Tables 5 -7.    The experimental results show that creating ensemble models from previously trained CNNs can improve the performance, comparing with a single NN. Ensemble models also outperforms single CNNs, designed only with one filter and give more reliability. What is more, there is not a statistical significance between these algorithms.

Conclusions
In this study, an automation of creating an ensemble of various CNNs is considered. All CNNs were designed using different types of image filters for increasing the quality of bitmap images and for detecting graphical patterns in images. The feedforward NN stage in all CNNs is the same for providing the equal conditions of the experiments. The experimental results have demonstrated that the performance of the proposed approach is comparable with the state-of-the-art techniques. Furthermore, the dimensionality of input data has been decreased, thus it is required smaller computational capabilities for implementing the CNN algorithm.
In further work, we will implement the proposed approach using the CUDA parallel computation framework for improving the computational performance and for carrying out more experiments. We will explore a wider range of existing image filters and will try to design new ones. More complicated structures of CNNs and ensemble models will be investigated with various image recognition problems.


extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Business

text:
Graph Theory Based Approach to Users Grouping and Downlink Scheduling in FDD Massive MIMO

Massive MIMO is considered as one of the key enablers of the next generation 5G networks.With a high number of antennas at the BS, both spectral and energy efficiencies can be improved. Unfortunately, the downlink channel estimation overhead scales linearly with the number of antenna. This does not create complications in Time Division Duplex (TDD) systems since the channel estimate of the uplink direction can be directly utilized for link adaptation in the downlink direction. However, this channel reciprocity is unfeasible for the Frequency Division Duplex (FDD) systems where different physical transmission channels are existent for the uplink and downlink. In the aim of reducing the amount of Channel State Information (CSI) feedback for FDD systems, the promising method of two stage beamforming transmission was introduced. The performance of this transmission scheme is however highly influenced by the users grouping and selection mechanisms. In this paper, we first introduce a new similarity measure coupled with a novel clustering technique to achieve the appropriate users partitioning. We also use graph theory to develop a low complexity groups scheduling scheme that outperforms currently existing methods in both sum-rate and throughput fairness. This performance gain is demonstrated through computer simulations.

I. INTRODUCTION
M OBILE traffic demand has never been as high as it is today due to proliferation of smart-phones, video streaming, and other data-hungry applications. The next generation mobile networks should, therefore, cope with the necessity of higher throughput. One of the promising technologies to enable this higher throughput is Massive MIMO [1]. Compared to the currently deployed multi-user MIMO systems, massive MIMO incorporates a much larger number of antennas. This technology was shown to provide better performance in terms of energy efficiency and overall capacity [2] which made massive MIMO a hot research topic and a key component of future standards.
However, the high number of antennas complicates the channel estimation and feedback. The downlink channel estimation overhead burden comes from the fact that it scales linearly with the number of antennas [3]. This is alleviated in TDD systems by exploiting the channel reciprocity since the channel estimate of the uplink direction can be directly utilized for the downlink direction [3]. This is not possible for FDD systems that still represent the vast majority of currently deployed cellular networks.
To deal with this difficulty, the authors in [4] proposed Joint-Spatial-and-Division-Multiplexing (JSDM), an approach to multiuser MIMO downlink that is considered one of the most promising candidates for FDD massive MIMO. It works by partitioning users with the same second order downlink channel statistics into groups and splitting the downlink beamforming into two stages: an outer precoder, that depends on the channel statistics, and an inner precoder that depends on the instantaneous effective channel realizations. The role of the precoders being to suppress inter-group and intra-group interference respectively. The dimensions of the effective channel is significantly less than the number of antennas, thanks to the outer precoder projection. Even with this reduction in CSIT feedback, the authors in [4] showed that JSDM achieves the same sum capacity of the corresponding MU-MIMO downlink channel if the eigenspaces of groups are mutually orthogonal, a condition they called "tall unitary".
In realistic scenarios, users might have similar but not necessarily identical second order downlink channel statistics. This dictates the incorporation of a clustering algorithm to partition users into groups with sufficiently similar covariance eigenspaces. On top of that, with a high number of users uniformly distributed across the cell, the eigenspaces of the groups are far from meeting the tall unitary condition and a reduction of the number of simultaneously served groups is required. These issues inspired the work in [5] where Kmeans clustering algorithm was adopted and a greedy sum-rate maximization scheduling algorithm was proposed.
Due to the greediness nature of the scheduling scheme proposed in [5] and to simplify users grouping, recent work [6] adopted a hierarchical clustering algorithm which mixes both target number of clusters and chordal distance threshold to reach an appropriate users clustering. Average Signal-to-Leakage-plus-Noise-Ratio (SLNR) based scheduling approach was also proposed and was shown to outperform in terms of sum-rate all the previous methods in the literature [6].
In our paper, we deal with the issues that are still present in the previous approaches. First, the target number of clusters is not known beforehand and choosing an arbitrary number of clusters can have severe impact on the performance of JSDM. On top of that, appropriate chordal distance thresholds are hard to predict in the clustering process. This inspired us to propose a novel similarity measure along with a new clustering scheme where the number of clusters is not required to be known. Taking into account the lack of orthogonality between the eigenspaces of groups, adopting a scheduling scheme is of great importance to enhance the overall performance of JSDM. However, groups scheduling that aims to improve the average SLNR as in [6] does not necessarily translate into a higher sum-rate. Also, fairness is an issue that was not previously addressed. For instance, when adopting the approach in [6], some groups suffer from starvation. We therefore develop, by using graph theory tools, a low complexity scheduling scheme that outperforms all currently proposed methods in both sumrate and throughput fairness.
The paper is organized as follows: Section II describes the system model. Section III presents the newly proposed metric and clustering method. Section IV includes the outer precoder design and the development of our scheduling scheme. Section V provides numerical results that demonstrate the performance of our method while Section VI concludes the paper.

II. SYSTEM MODEL
We consider a single cell downlink multiuser MIMO system with N t antennas at the base station and K single-antenna users. Let y ∈ C K×1 be the received signal by the users: where x ∈ C Nt×1 is the transmitted signal vector, z ∈ C K×1 denotes the additive white Gaussian noise vector and H ∈ C Nt×K is the channel matrix. The transmitted signal vector is actually a precoded version of the data vector x = V d where V ∈ C Nt×S is the precoder and d ∈ C S×1 is the data vector. The dimension S is equal to the number of total independent streams and is upper bounded by min{N t , K} [4]. For the sake of simplicity, we adopt the approach of [4] with equal power allocation i.e. E(dd H ) = P S I S where P is the total downlink power budget. We suppose that z ∼ CN (0, I K ). Assuming no line-of-sight propagation, we have h k ∼ CN (0, R k ) where R k is a positive semi-definite covariance matrix. Let the eigenvalue decomposition (EVD) of R k be the following: where Λ k is an r k × r k diagonal eigenvalues matrix with r k being the rank of R k and U k ∈ C Nt×r k being the set of eigenvectors corresponding to the non-zero eigenvalues. We adopt the two-stage precoders approach proposed in [4] where V = BP . Based on the similarity of their channel covariance, users are partitioned into G groups with each containing K g users such as K = G g=1 K g . The outerprecoder B, of dimensions N t × b, is designed in a way to minimize the inter-group interference based on the channel statistics, which is supposed to be known at the base station as adopted in [4] 1 . The inner precoder P , of dimensions b × S, depends on the instantaneous channel realizations and is intended to suppress intra-group interference. By taking into account the partitioning of users into G groups, we have 1 The channel statistics vary at a much slower rate than the channel coherence time and therefore can be assumed to be locally stationary and easily tracked by methods cited in [4] the following: H g = [h g1 , . . . , h ], P = diag{P 1 , . . . , P G } and we define the effective channel H = B H H. It is straightforward that the effective channel is of dimension b×K with b = G g=1 b g and b g N t . This drastically reduces the amount of CSI feedback in the case when each user g k feedback his effective channel h g k ∈ C bg×1 rather than h g k ∈ C Nt×1 . We will refer to this approach as per-group processing (PGP). The received signal by group g can be therefore written as: where d g ∈ C Sg×1 with S g being the number of independent streams intended for group g. By adopting the PGP approach and assuming perfect effective CSI at the BS, a zeroforcing (ZF) inner precoder can be calculated in the following manner: with ζ g being a normalization factor to ensure that the power budget constraint is satisfied:

III. CORRELATION CLUSTERING
In this section, we deal with the fact that users might have similar but not necessarily identical covariance matrices and appropriate grouping of users is essential for JSDM. The research papers that investigated this clustering problem presented two approaches: K-means clustering [5] and a hierarchical clustering [6]. Both of these approaches used the chordal distance as a metric. The downside of such a metric is the fact that prediction of any threshold involved in the clustering algorithm is a difficult task. Motivated by this, we propose in the following subsection a new similarity measure suitable for our problem.

A. Similarity Measure
Herdin et al. [7] introduced a novel metric named Correlation Matrix Distance (CMD). It was used to track the changes of spatial structures of the channel in non-stationary MIMO. The use of this metric has been extended to many different research work. For example, the authors in [8] used it in the context of Grassmannian subspace packing. The same metric was also adopted by the authors in [9] to study the effect of subspace alignment in multi-user MIMO. The fact that CMD is normalized makes it more sensitive to differences in the correlation structure and present an opportunity in terms of threshold design. In the previous literature [5][6], the similarity between user 1 and user 2 was solely taken based on their covariance's eigenstructures (U 1 U H 1 , U 2 U H 2 ) without taking into account the energy of the modes. In our case, we will be applying our similarity measure on the whole covariance matrices (R 1 , R 2 ). The motivation behind this is that differences in the eigenstructures of weak modes should contribute less than the one's of strong modes. Based on CMD, we can define the new similarity measure as follows: This similarity measure is lowerbounded by 0 and upperbounded by 1. A value of 0 corresponds to the case where R 1 and R 2 are orthogonal while a value of 1 takes place when R 1 and R 2 are collinear. This proposed measure can be regarded as an extension of the cosine similarity of vectors (which is a widely used metric in clustering algorithms see, e.g., [10]) to matrices and therefore can now be considered as what we will call Degree of OverLap (DOL) between the two spaces. To our knowledge, this is the first time it has been used in the context of users clustering for FDD massive MIMO.

B. Clustering Algorithm
Unlike the previously proposed approaches, we seek to use a clustering algorithm without passing the target number of clusters as a parameter. To do so, we take advantage of the ease of threshold design presented by our proposed similarity metric. An interesting way to do so is by choosing DOL th high enough such as if d s (R k , R k ) ≥ DOL th then users k and k can be thought to be laying in the same correlation space. Unlike other similarity metrics, this threshold is easily determined. One can simply say if the degree of overlap between the two spaces is above 0.95 then consider them as highly similar and are preferred to be assigned to the same cluster. Based on this, we can construct what we will call a complete advice graph G c = (V c , E c ) where each vertex represents a user and an edge would have a +1 label to signal that these two users are preferred to be in the same cluster while −1 label refers to the opposite case.

Fig. 1: Clustering Algorithm
Our goal now is therefore to produce a partition of the graph's vertices in a way that agrees as much as possible with the edge labels. To do so, we propose a cost function J as the total disagreements of our resulting partitioned graph. The total disagreements cost is defined as the overall negative weights inside a cluster added to the positive weights between clusters. Our partitioning problem can be hence formulated as follows: The constraints take into account the symmetry of x uv and the triangular inequality 2 satisfied by these variables. What makes this clustering formulation interesting is that there is no need to specify the target number of clusters. Instead, the resulting optimal number of clusters could be any value from 1 to K depending on what fits our graph the most. Our problem in (7) turns out to be the same as the one studied in [11] [12]. In general, solving (7) and finding the optimal clustering is NP-hard, as proven by the authors of [11] using a reduction from Exact Cover by 3-Sets (X3C) which is one of Karp's 21 NP-complete problems. One way to deal with this hardness is to turn the problem into a Linear Program LP by relaxing the binary condition and replacing it by The LP is then solved in polynomial time by any desired standard LP solvers followed by appropriate rounding of the fractional values. The question that remains is how to round the fractional values? The literature is rich with rounding techniques that achieves a good approximation ratio compared to the optimal solution. The most recent work in [12] proposed a new randomized technique based on Pivoting. Pivoting works by treating the fractional solution of (7) as a probability to put the two vertices in different clusters. The algorithm that was proposed in [12] is to apply the following function on the solution of (7) before proceeding to the pivoting phase: respectively. This rounding technique is guaranteed to achieve an expected (2.06-)-approximation for a = 0.19, b = 0.5095 and a constant such as 0 < < 0.01. A derandomized version of the algorithm was also considered but we omit it for the sake of space and we refer the readers to [12]. Overall, the clustering algorithm can be summarized in Algorithm 1.
Pick a pivot w t ∈ V t randomly and let S t = w t 10: ∀ u ∈ V t , add u to S t with probability 1 − p wu 11: let V t+1 = V t \ S t , let t = t + 1 12: end while 13: end procedure 14: Output the clusters S 0 , . . . , S F inal IV. DOWNLINK SCHEDULING After grouping users with similar second order channel statistics, we can now deal with the orthogonality aspect of JSDM. In realistic scenarios, groups do not lay in mutual orthogonal channel covariance spaces and inter-group interference can therefore limit the overall performance. One can seek to reduce this interference by applying appropriate outer precoding techniques but it is insufficient as will be proven shortly. Therefore, adopting a scheduling scheme is of paramount importance for the overall performance of JSDM.

A. Problem Formulation
We define the centroid that would be taken as a representative of each group's equivalent covariance: where r g is the rank of R g . This is where the clustering algorithm's effect is highlighted. Due to our clustering algorithm and with a high DOL th , we know that the covariances of users in each group are really similar which means the centroid is a good representative of each cluster. However this is not necesarrily true in the case of other clustering algorithms with pre-determined target number of clusters as in [5][6] since it is hard to predict beforehand the right number of clusters for which each group's representative is a good one. As proposed in [4] and adopted in [6], one can seek to eliminate inter-group interference by simply building up the inter-group interference matrix for group g and projecting the intended signal space on the interference matrix's orthogonal space, a method that was given the name Approximate Block Diagonalization. The interference matrix seen by group g is Ξ g = [U * 1 , . . . , U * g−1 , U * g+1 , . . . , U * G ] which is based on the eigenspace of other active groups with U * g ∈ C Nt×r * g where r * g is a design parameter which denotes the number of the channel's strongest modes 3 taken into account. Let [E (1) g , E (0) g ] denotes the set of left eigenvectors of Ξ g where E (0) g of dimensions N t ×(N t − g =g r * g ) form a unitary basis for Span ⊥ (U * g : g = g). The projected channel covariance matrix is given by: The next step would be to match the b g strongest eigenmodes of our projected channel.
Overall, the outer precoder is B g = E (0) g G (1) g . One would say that we can simply put r * g = r g i.e. includes all modes and we will lay in an interference free scenario. But by construction, we have that the effective channel dimension b g ≤ rank( R g ) = min(r g , N t − g =g r * g ), which means including more modes would shrink our dimensionality and leads to a dimensional bottleneck. Keeping in mind that S g ≤ b g , the dimensionality bottleneck is a serious matter since we have a certain number of independent streams that we are obligated to send out for each group. Since not necessarily all modes are included (r * g ≤ r g ), inter group interference would still be inevitable and we therefore include downlink scheduling to enhance the performance of the system. To formulate our scheduling problem, we have to adopt a certain network utility function.
The first building block of any network utility is the expression of the rate R g k achieved by each user g k in the network. The expression of the rate can be concluded from the SINR expression by applying R g k (SIN R) = log 2 (1+SIN R g k ) [4].
Large scale analysis: We focus in our paper on the case of multi-user massive MIMO where the number of antenna and users N t , K −→ +∞. In this case, Random Matrix Theory (RMT) tools come in handy. The authors of JSDM [4] made use of the work in [13] to propose a deterministic equivalent for the SIN R expression in JSDM. Motivated by the fact that this deterministic equivalence was shown to be accurate for realistic values of (N t , K) [4][13], we take it as a basis of our analysis. Details concerning these equations can be found in [4] which for our case reduce to: where x g is a binary variable that denotes if group g is scheduled. ζ 2 g = m g b g , Υ g,g and m g are the results of fixed point equations with R g = B H g R g B g : What makes those expressions interesting is the fact that the effect of small-scale fading is averaged out as predicted in [2]. This is quite convenient since the second order statistics change at a much smaller rate than the channel coherence time as previously pointed out in Section II.
With the rate expression approximation dealt with, the goal now is to schedule these clusters in a way to get the highest utility, while preserving fairness between clusters and ensuring a certain quality of link for users in each group. Motivated by facilitating the design of the quality of link tolerance, we consider as a criterion the Signal-to-Interference Ratio SIR g k experienced by scheduled users which is a widely used criterion in power control for wireless cellular networks [14].
The SIR g k is forced to be above a certain threshold beyond which the link to user g k is supposed to be good. Since the SIR g k of user g k depends solely on the group index g, we will drop the sub-index g k and work with SIR g experienced by the group g. Putting it all together, and taking the weighted sum-rate as utility, we can formulate our scheduling problem: The previous studies in this area have not considered fairness between users (See, e.g., [6]). Introducing this weight w g k allows us to incorporate fairness in our scheduling scheme. A special case would be the stable policy max-weight scheduling [15] where w g k is chosen to be the length of the queue Q g k .

B. Scheduling Scheme
In order to solve our problem in (17), we propose a 2steps scheme based on graph theory. The first step deals with both the SIR constraint and a combinatorial difficulty faced in our problem. The second step aims to find the appropriate combination of groups to be able to solve (17). From the SIR g expression in (16), we can construct a weighted directed graph G L = (V, E) where V is the set of groups and in which the weight of the edge e(g , g) corresponds to what we will call the normalized interference from group g to group g: e(g , g) = ζ 2 g Υ g,g ζ 2 g (18) One thing to point out is that due to the fact that the outer precoder B g depends on the eigenspace of all active groups, the weight of the edges in this graph are not constant and depend on the activity of all groups as seen from the fixed point equations (12)- (15). This complicates our problem since in normal graphs, the weight of the edges is fixed and do not change as you manipulate the graph. To alleviate this issue and deal with the SIR constraint in (17), we first proceed to what we will call the "Elimination" phase. This phase also allows us to convert our problem into a vertex coloring problem, as will be detailed in the "Grouping" phase. An example of 4 groups scenario will be presented in successive figures to fully demonstrate the mechanisms of the scheduling scheme. Fig. 2: Weighted Directed Graph G L 1) Elimination: We can picture each vertex in the graph as a sink of interference that undergoes successive iterations. In the first iteration, all groups are considered to be active. The outerprecoder B g of each group is calculated as detailed in Section IV-A. The fixed points equations (12)-(15) are then solved and e 1 (g , g) ∀(g , g) ∈ V 2 are calculated based on (18) where the sub-index "1" refers to the iteration number. For each vertex g ∈ V , the SIR g condition of (17) is tested.
If it is violated, the edge e 1 (g , g) with the highest weight is eliminated. In other words, the group that interfere most with g is chosen to be eliminated. At the next iteration, we have a new graph due to the edges removal from the previous iteration and therefore the outer precoders are to be recalculated. This time however, the outer precoder B g of each vertex g ∈ V is calculated based on the eigenspace of neighboring 4 vertices only. We repeat the same procedures of the first step: the fixed points equations are solved again (12)- (15) and the weight of the edges of neighboring vertices only are recalculated using (18). The stopping criteria would be if an iteration resulted in no new deleted edges. In other words, if simultaneous scheduling of neighboring vertices in the resulting graph will not vioalate the SIR condition of each of them. An example of the above procedure is given in Fig. 3, the first iteration resulted in four deleted edges. The edges of neighboring vertices are then updated for the second iteration. The second iteration did not result in any deleted edges and the algorithm finishes. Once the algorithm finishes, we turn our graph into an undirected graph G u = (V, E u ) by simultaneous agreements from both sides i.e. if e(g, g ) and e(g , g) are both not eliminated in G L then an edge e u (g, g ) = 1 exist in G u and e u (g, g ) = 0 otherwise. In our new undirected graph G u , an edge exist between two vertices if scheduling them together will not violate their respective SIR conditions. We can now tackle another aspect of our problem: Which groups of those that are allowed to transmit simultaneously should we schedule in order to maximize our utility?

Fig. 3: Elimination Process
2) Grouping: After proceeding with the elimination step, our SIR constraint can be replaced by making sure that two simultaneously scheduled groups should have an edge between them in G u . Therefore, our problem in (17) is turned into: To solve this new problem, we recall that for a well chosen α g ∀g, an edge exist between two vertices in G u = (V, E u ) only if they barely interfere and hence scheduling them together would normally increase their sum utility. Our goal is therefore to find combinations of groups that are adjacent one to the other in G u while covering the whole vertex set V . For this purpose, we define a clique in an undirected graph as a subset of vertices such as every two distinct vertices in the clique are adjacent. We seek to find the smallest number of cliques that cover V , where we emphasize "smallest" to ensure that each clique have the largest number of groups possible inside. Essentially, we are trying to solve the minimal clique vertex cover problem. The minimal clique vertex cover problem is known to be equivalent to vertex coloring on the complement graphḠ u , a well known NP-Complete problem.
Knowing that vertex coloring seeks to partition the set of vertices into the smallest number of independent sets, one can see the connection between the two problems since a subset of vertices is a clique in G u if and only if it is an independent set inḠ u . We will use a simple yet effective maximal independent set based vertex coloring algorithm that achieves a O( n log(n) )approximation ratio [16] presented in Algorithm 2 and apply it onḠ u . After applying Algorithm 2, each group g will be assigned a color Col(g). Groups that are assigned the same color represent a subset of groups that are allowed to transmit simultaneously. We can therefore replace the edges constraint in (19) by ensuring that the color assignments are respected: The problem in (20) is indeed simple to solve. One can simply form "Schedules", each made of groups that belong to the same color. At the start of each coherence time T c , the schedule (i.e. the color) that leads to the largest utility is selected. Let R = V \ S the set of remaining vertices 5: while R = ∅ do 6: Pick a vertex w ∈ R with lowest degree randomly and let C i = C i ∪ {w} Output the color C i and i ← i + 1 11: end while

V. NUMERICAL RESULTS
We consider a base station with a 120 • sector centered around the x-axis consisting of a ULA with N t = 128 antennas and serving K = 80 users arbitrarily distributed in the sector. For ease of correlation entries calculation, we adopt the onering model [4]. We consider a user terminal (UT) at an azimuth angle θ and angular spread ∆. The correlation entry is then calculated for 1 m, p N t using the following formula: where k(α) = −2π λ (cos(α), sin(α)) T denotes the wave vector for a planar wave with angle of arrival α, λ is the wavelength and u m , u p ∈ R 2 are the position vectors of the BS antennas in the 2D-coordinate system. For our scenario, we consider that all users have the same angular spread of ∆ = 5 • . We set the clustering high threshold to DOL th = 0.95. We suppose that S g = K g and set r * g = K g as adopted in [6]. The goal of the simulations is to compare our method to the recently proposed SLNR based scheduling scheme [6]. We take the SLNR based scheduling scheme as a benchmark since it was shown to outperform all proposed methods in the litterature in terms of sum-rate [6]. We also expose the necessity of a scheduling policy by simulating JSDM without any scheduling just as in [4]. Due to the fact that both the SLNR based scheduling and our scheduling scheme require a threshold tolerance, we iterate over a wide range of thresholds and choose the one that led to the highest sum-rate as a representative of each methods for a fair comparison. To be able to compare our method to theirs, we create our schedules and switch between them in Symmetrical Round-Robin manner and take their average sum-rate as a representative while baring in mind that higher sum-rate can be achieved for our scheme by simply choosing the schedule with the highest sum-rate. Fig. 5 shows how the JSDM method without scheduling performs poorly due to limitations in terms of inter-group interference which is a proof of the necessity of adopting a scheduling scheme. In addition, we can see how our method was able to outperform the SLNR based method over the whole SN R range. This is because we work on the inter-ference itself to improve the sum-rate. On the other hand, increasing the average SLNR of a system as in [6] does not necessarily translate into a higher sum-rate. Fig. 5: Comparison of sum spectral efficiency vs. SNR For the comparaison in terms of throughput fairness between the methods, the well known Jain's fairness index [17] was taken as a metric. It is defined as follows: The values of this index range from 1/K (the case where only one user acquire the channel while the others end up starving) to 1 (the case where resources are shared equally between users). Fig. 6 shows that the SLNR based scheduling scored the worst fairness index due to the fact that after successive elimination of groups with low SLNR, the users inside these groups end up starving. The throughput fairness of JSDM with no scheduling is high but not perfect since users suffer different interference conditions and therefore asymmetric throughput. Our method scored almost perfect throughput fairness due to two reasons: the first being that by construction, SIR of each group was chosen to be lower bounded by the same well chosen tolerance and the second being that symmetrical Round-Robin was adopted between schedules. Overall, our proposed scheme was able to outperform the SLNR based method in sum-rate while providing a huge gain in terms of throughput fairness. VI. CONCLUSION In this paper, we tackled the problem of users clustering and scheduling in the promising technique of two-stage beamforming for the downlink of FDD massive MIMO. We introduced a new similarity metric coupled with a clustering method that are characterized by ease of design and good performance. We also developed using graph theory tools an efficient groups scheduling scheme that outperforms the current methods proposed in the literature in both sum-rate and throughput fairness as was shown in the simulations.