Title: LEARNING KOLMOGOROV MODELS FOR BINARY RANDOM VARIABLES Anonymous

Abstract: We propose a framework for learning a Kolmogorov model, for a collection of binary random variables. More specifically, we derive conditions that link (in the sense of implications in mathematical logic) outcomes of specific random variables and extract valuable relations from the data. We also propose an efficient algorithm for computing the model and show its first-order optimality, despite the combinatorial nature of the learning problem. We exemplify our general framework to recommendation systems and gene expression data. We believe that the work is a significant step toward interpretable machine learning.

1. INTRODUCTION:

Machine learning and artificial intelligence method have permeated a large number of areas (Marr, Sept 2016) . These methods are based on machine learning models, which consist of learning an input-output mapping for a given dataset. Despite the plethora of such models (e.g., matrix factorization (Koren et al., 2009) , SVD-based models (Koren, 2008) , deep neural networks (LeCun et al., 2015) , and models inspired from physics (Stark, 2016b )), they lack interpretability: they are not capable of offering insight about the data, nor the underlying process. The lack of interpretablity may have serious consequences in mission-critical systems, ethics, and validation of computeraided diagnosis (Doshi-Velez & Kim, 2017) . While there is no consensus around the definition of interpretability, causality (Lipton, 2016) is a vital component: it refers to causal relations within the data and insight about the underlying data-generating process. We thus adopt a generalized version of causality using implication in mathematical logic as our 'definition' of interpretable models.

To this end, we propose learning a so-called Kolmogorov Model (KM) associated with a set of binary Random Variables (RVs). In addition to prediction, the interpretability of the model (as defined above) enables learning causal relations. 1 We derive a sufficient condition under which the realization of one RV's outcome (causally) implies the outcome of the other. In the context of recommendation systems, causal relations identify groups of items, for which a user liking one item implies that he/she likes all other items in the group. In gene expression analysis, the same relations identify groups of DNA locations for which the expression of a gene in one of them implies its expression in all other locations. The foundation of our approach is to model binary RVs as elementary events on a Kolmogorov space, by an inner product of two vectors (formally stated in Section 2.1), which is based on established results from classical probability theory. To our best knowledge, this specific formulation is novel in the context of learning representation.

The inner product of our formulation is reminiscent of factorization methods such as, matrix factorization (Koren, 2008) , non-negative matrix factorization (Lee & Seung, 2001) , SVD (Cai et al., 2010) , and physics-inspired techniques, non-negative models (Stark, 2016a) . It should be noted that the inner product components of these methods is usually based on strong intuition about the data, provided and validated by a human expert. In contrast, the proposed KM is a fully automated approach, deeply rooted in probability theory, which finds an interpretable model (as defined above) of the data without any human intervention. Thus, both the model and the causal relations within the data have a strong mathematical basis. Moreover, in most of the existing approaches, we may need to have one pipeline for learning the representation and another one for mining these relations. However, our approach requires only a single pipeline for both tasks. This is very appealing on a conceptual level, and it can drastically simplify coding and validation. Our method also generalizes K-means and some of its variants. Detailed discussions of the relation between our proposed method and these prior works is in Appendix A.3.

In an abstract sense, we formulate a KM learning problem as a coupled combinatorial program, decompose it into two subproblems using the Block-Coordinate Descent (BCD) method, and obtain provably optimal solutions for both. For the first one, we exploit the structure of linear programs on the unit simplex and use low-complexity (yet optimal) Frank-Wolfe algorithm (Frank & Wolfe, 1956) . To bypass the inherent complexity of the second subproblem (combinatorial and NP-hard), we propose a semidefinite relaxation, and show its quasi-optimality in recovering the optimal solution of the combinatorial subproblem. Finally, we show the convergence of our algorithm to a stationary point of the original problem. We propose a simple algorithm for mining the causal relations. All the proofs and additional discussions are available in Authors (Oct 2017).

2. SYSTEM MODEL:

Notation: Lowercase letter a, uppercase bold letter A, and calligraphic letter A denote vectors, matrices, and sets, respectively. [A] i,j and A T denote element (i, j) and the transpose of A. supp(a) denotes the support set of a. The inequality x ≤ y holds element-wise. I denotes the identity matrix, 1 and 0 the all-one and all-zero vectors (of appropriate dimension). e n is the nth elementary basis,

P = {p ∈ R D + | 1 T p = 1} the unit simplex, and {n} := {1, • • • , n}. 2.1 PROBLEM FORMULATION Consider a double-indexed set of binary Random Variables (RVs), X u,i ∈ A = {1, 2}, taken from a dataset D = {(u, i) | (u, i) ∈ U × I}. Each RV is defined on a sample space Ω, consisting of elementary events Ω = {ω d | 1 ≤ d ≤ D}.

We denote by P[X u,i = z], z ∈ A, the probability that RV X u,i takes the value z ∈ A; see example in Section 2.2. Since X u,i is binary, it is fully characterized by considering one outcome. Thus, we write the Kolmogorov Model (KM) for X u,i as

EQUATION

Thus, each RV X u,i is associated with (characterized by) a Probability Mass Function (PMF) θ u , u ∈ U, and an indicator vector ψ i , i ∈ I. The model follows from established results in classical probability theory (Kolmogorov, 1957) , (Gray, 2009) 

= {(u, i) | (u, i) ∈ U × I}. 2

Given samples coming from the model in (1), we wish to infer the parameters of underlying probability distribution: find parameters of the KM, i.e., {ψ i , θ u } that best describe {p u,i | (u, i) ∈ K}. The resulting problem is a fully parametric statistical inference task. For tractability, we address it using the minimum mean-squared error as a point estimator. The corresponding optimization problem is

EQUATION

)

We present our solution to this coupled non-convex conspiratorial optimization in the next section.

The obtained solution to (Q) can be used for prediction on a test set, as well as extracting causality structures within training set (see Section 4).

Proposed Approach: While the proposal to model binary RVs as elementary events on a Kolmogorov space is based on established results, the specific learning formulation (Problem 1) is novel. Because this model is rooted in probability theory, (1) defines the outcome of a RV in the strict Kolmogorov sense, and the resulting causal relations (Section 4) also hold from a strict analytical perspective. Note that causal relations (a.k.a. association rules) may still be extracted using existing methods, e.g., (non-negative) Matrix Factorization (MF) and their variants, SVD, binary MF, and K-means. However, these relations are not based on causality and the formal relations that (mathematically) follow from the KM in (1), but rather on intuitions/heuristics, which may yield different relations. Naturally, we wish the explore further in this work the causal relations that arise from the proposed model. Additionally, unlike existing methods, the prediction and causal relations mining are done in 'one-shot', thereby simplifying the implementation/validation; see Appendix A.3.

2.2. ILLUSTRATIVE EXAMPLE: RECOMMENDATION SYSTEMS:

In this context, U and I denote the set of users and items respectively, and X u,i models the preference of user u for item i, (u, i) ∈ K. Thus, P[X u,i = 1] (or P[X u,i = 2]) is the probability that user u likes (or dislikes) item i. Moreover, θ u determines the profile/taste of user u, ψ i is related to item i (depending on genre, price, etc.), and the elementary events denote movie genres (e.g., ω 1 = "Action", ω 2 = "SciFi", etc.). The training set, consisting of an empirical probability that user u likes item i, p u,i . We obtain this probability using p u,i [R] u,i /R max , where [R] u,i ∈ N denotes the rating that user u has provided for item i, and R max the maximum rating (Stark, 2015) . For instance, if user u rates item i with a score of [R] u,i = 7 (where the maximum rating is 10), then p u,i = 0.7; Other approaches may be used to obtain p u,i depending on the specific application. As a concrete illustrative example, consider a 10-star "recommendation system", having 2 users and 2 items. We then find the D-dimensional (D = 3) KM factorization to obtain , {ψ i } 2 i=1 and {θ u } 2 u=1 . D is the size of the Kolmogorov space Ω, the number of elementary events, and the dimension of the factorization (selected via cross-validation to minimize the test error). Solving (Q) results in finding {ψ i } 2 i=1 and {θ u } 2 u=1 . An example result is given below:

0.3 1 0.1 1 {pu,i} = θ T 1 0.2 0.3 0.5 θ T 2 0.1 0.1 0.8     0 1 }Action 1 1 }SciFi 0 ψ 1 1 ψ 2 }Drama    

To showcase the model's intuition, note that p 1,1 , the probability that user 1 likes movie 1, is represented as ψ T 1 θ 1 . It is thus expressed as convex/stochastic mixture of movie genres, since elementary events are movie genres in this scenario. We underline that this high degree of interpretability is not artificially enforced but rather context-dependent. More generally, a KM represents a set of observed outcomes for RVs as (context-dependent) mixtures of elementary events. The approach consists of learning a (hidden) latent model by jointly optimizing the PMF and binary indicator vectors. Thus, interpreting elements of the indicator vectors as movie genres requires a context-dependent map, from the elements of ψ i to movie genres. Another way to find such an interpretation is when this context-dependent map is known apriori. In this setting, each item is already tagged with its movie genres, and ψ i need not be optimized; Alas, having this context-dependent map comes at the expense of a loss in training/test (see Appendix A.4). However, recall that the primary interest of this work is not this facet of interpretability, but rather by that of the causal relations.

3. PROPOSED ALGORITHM:

To approach a solution for problem (2), we use the block-coordinate descent (BCD) method to split (Q) into two sub-problems. Here, we derive our solution approach to each. We first refine the current PMF estimation θ u , and then that of the indicator ψ i . Given {ψ (n) i } at iteration n, we pose the PMF refinement, θ-step, as

EQUATION

We then pose the indicator vector refinement, ψ-step, as

(Q 2 ) : ψ (n+1) i ∈ argmin ψi∈B D g(ψ i ) ψ T i S (n+1) i := u∈U K θ (n+1) u θ (n+1) T u ψ i -2ψ T i v (n+1) i := u∈U K θ (n+1) u pu,i . (4) function [θ u ] = FW ( Qu, ru, ) for k = 1, 2, ..., IF W do d (k) u = ej , j = argmin 1≤j≤D [∇f (θ (k) u )]j θ (k+1) u = (1 -α (k) u )θ (k) u + α (k) u d (k) u Stop if θ (k+1) u -θ (k) u ≤ end for end function Table 1: θ-step solution using FW function [ ψi] = SDR ( Si, ti, M rnd ) // Repeat to approximate each ψ i , ∀i ∈ IK Solve (5) to find X (SDR) i Factorize as X (SDR) i = L T i Li for m = 1, 2, ..., M rnd do Gen. zero-mean i.i.d Gaussian vector u (m) i Compute û(m) i = sign[ L T i u (m) i ] end for Find m = argmin 1≤m≤D+1 û(m) T i Si û(m) i Compute ẑi = [u (m ) i ]1:D [u (m ) i ]D+1 Approximate ψ i , as ψi = ( ẑi + 1)/2 end function Table 2: ψ-step solution using SDR Moreover, U K and I K are defined as K = {(u, i) | u ∈ U K ⊆ U , i ∈ I K ⊆ I}.

Recall that having globally optimal solutions for both (Q 1 ) and (Q 2 ) is necessary to show the convergence of BCD (Tseng, 2001 ) -a challenging task due to the NP-hardness of (Q 2 ).

3.1. θ-STEP: REFINE PMF ESTIMATE:

We use the Frank-Wolfe (FW) algorithm (Frank & Wolfe, 1956) to solve (Q 1 ) as a succession of Linear Programs (LPs) over the unit simplex. While LP solvers generally have similar complexity as quadratic program solvers, solving an LP reduces to searching for the minimum index -which is computationally efficient, when the LP is over the unit simplex. We summarized this FW variant for solving (Q 1 ); see also Jaggi (2013) [Algorithm 1]. Here, we drop the BCD iteration index, n, and just keep the FW iteration number, k, for notation simplicity. We first determine the descent direction:

d (k) u ∈ argmin s∈P ∇f (θ (k) u ) T s.

The constraint s ∈ P greatly simplifies the above LP and yields: d

(k) u = e j , j ∈ argmin 1≤j≤D [∇f (θ (k) u )] j ; see Proposition 2 (Appendix A.6).

The solution follows from LPs over the unit probability simplex. Thus, finding the descent direction reduces to searching over the D-dimensional gradient vector (done in O(D)). Then, the current value is updated using a simple step size rule, α (k) u = k/(k + 1). Table 1 shows the θ-step solution, and Proposition 3 in Appendix A.6 characterizes its convergence.

3.2. ψ-STEP: REFINE INDICATOR ESTIMATE:

To address the NP-hard nature of (Q 2 ), we propose a solution based on Semi-Definite Relaxation and Randomization (SDR), and establish its quasi-optimality. We use the results of Ma et al. (2002) [Sec IV-C]) and a series of reformulations to rewrite (Q 2 ) in the following equivalent form (see Authors (Oct 2017) for all the derivations):

X i ∈ argmin Xi tr( Si X i ) , s. t. X i 0, [X i ] k,k = 1, ∀k , rank(X i ) = 1 where X i = x i x T i , Si = (1/4)S i -ti /2 -tT i /2 0 , x i = z i w i , z i = 2ψ i -1, w i ∈ {-1, +1} is an

auxiliary variable, and ti (v i -(1/2)S i 1). The problem is then relaxed into a convex program,

EQUATION

X (SDR) i may be solved using generic solvers. Then, a randomization procedure (Ma et al., 2002 ) extracts an approximate (binary) solution for (Q 2 ); see Table 2 . This evidently raises the issue of the sub-optimality gap for SDR. Based on the results of Tan & Rasmussen (2001) and Luo et al. (2010) , we show in Proposition 4 (see Appendix A.6) that SDR (Table 2 ) is optimal (asymptotically in D) in recovering the binary solution of (Q 2 ).

We highlight that the performance bound in Proposition 4 compares the quality of both the approximate binary solution offered by SDR (with respect to the optimal binary solution of (Q 2 )), as well as their respective cost functions. The asymptotic optimality is empirically validated in Appendix A.7.

3.3. ALGORITHM DESCRIPTION:

The BCD-based algorithm alternates between refining the indicator and PMF vectors; see Algorithm 1. Lemma 2 (Appendix A.6) shows its convergence to a stationary point of (Q).

Algorithm 1 Iterative computation of KMs // Randomly Initialize {θ

(1) u ∈ P} for n = 1, 2, ... do Compute S (n) i and t (n) i using (4) Call ψ(n) i = SDR(S (n) i , t (n) i , M rnd ), ∀i ∈ I K Compute Q (n) u and r (n) u using (3) // Initialize FW with {θ (n-1) u }, from previous iteration Call θ (n) u = FW(Q (n) u , r (n) u , ), for all u ∈ U K end for 4 INTERPRETABILITY VIA CAUSAL RELATION 4.1 CAUSAL RELATIONS

Once a solution is found using Algorithm 1, here, we propose a method to find causal relations among the RVs. More specifically, we compare the support set of each pair of RVs from the training set, X u,i and X u,j , and check if there is any 'overlap' between their support set. Intuitively, this condition means that some of the elementary events (see Section 2.1) of one RV might be 'contained' in the elementary events of another. Consequently, the RVs are mutually related by causality, and the outcome of one determines that of the other. This insight is formalized here.

Proposition 1 (Inclusion of Support Set) Consider two random variables X u,i and X u,j (belonging to the training set) whose KM are given by the model in (1). If supp(ψ j ) ⊆ supp(ψ i ), then the following two causal relations hold: X u,i = 1 implies X u,j = 1 , and X u,j = 2 implies X u,i = 2 .

(6)

Proof: See Authors (Oct 2017) for the proof. Stated plainly, when the support set condition holds, the first outcome of X u,i implies the same outcome for X u,j , and the second outcome for X u,j implies the second one for X u,i , thereby implying a mutual causal relation among them (since X u,i influences X u,j and vice-versa). Note that our above definition of causality and causal relations is different than conventional ones in Pearl (2009) [Chap 2.8]. Moreover, our definition is distinct from Granger causality, due to the mutual coupling among the RVs in question.

We present next a special case of Proposition 1. When ψ i = 1, then supp(ψ i ) = {D}, and supp(ψ j ) ⊆ supp(ψ i ) holds, for any choice of ψ j , ∀ j ∈ I K , where I K defined in Section 3.

Corollary 1 (Maximally Supported RVs) Let {ψ i , θ u } (u,i)∈K denote the KM associated with {P[X u,i = 1]} (u,i)∈K . We define M as set of RVs for which the support set of the indicator vector is one, i.e., M = {i | ψ i = 1}. Then, the condition supp(ψ j ) ⊆ supp(ψ i ) (Proposition 1) holds trivially ∀ j ∈ I K . Thus, the causal relations in (6) hold, for every i ∈ M.

For maximally supported RVs, the realization of one outcome, X u,i = 1, determines that of all RVs of the set {X u,j = 1 | ∀j ∈ I K }.

Example 1 (Causal Relations in Recommendation Systems) In addition to prediction, recommendation systems are designed to accurately mine for association rules: if a user likes item i will he/she like another item j? Thus, the causal relations of Proposition 1 and Corollary 1 will be quite powerful, as we shall see. In the example of Section 2.2, note that supp(ψ 1 ) ⊆ supp(ψ 2 ). Then, Proposition 1 yields: if user 1 (or user 2) likes movie 2 implies he/she also likes movie 1. Moreover, X 1,2 and X 2,2

X u,3 X u,1 X u,2 X u,3 X u,1 X u,2 X u,3 X u,1 X u,2 a) b) c)

a) Check support set condition between Xu,1, and Xu,2, Xu,3 (Xu,1 is a maximally supported RV here) b) Check support set condition between Xu,2, and Xu,1, Xu,3 c) Check support set condition between Xu,3, and Xu,1, Xu,2 Thus, our approach provides association rules that follow from causal relations between different RVs. Consequently, these relations are stricter (as they are rooted in Kolmogorov probability) than other methods for mining association rules, which are not based on causality.

4.2. CAUSAL RELATIONS MINING (CRM):

We provide an efficient algorithmic approach to automatically mine the above relations. In a nutshell, the algorithm does a pairwise check of the support set condition (Proposition 1), for pairs of RVs X u,i and X u,j , ∀(i, j) ∈ I K × I K . The method is illustrated in Figure 1 , for the illustrative example of Section 2.2. The above causal relations can be modeled conveniently using the adjacency matrix A ∈ B |I K |×|I K | , defined as

[A] i,j = 1, if supp(ψ j ) ⊆ supp(ψ i ), and 0 otherwise, ∀i = j . (7) Stated differently, [A] i,j = 1 if the support of X u,j is contained in that of X u,i .

Note that when A is sparse (resp. dense) is an indication for a dataset in which few (resp. many) casual relations exist. To quantify the sparsity level, we define an influence score,

β i = |I K | -1 j∈I K j =i

[A] i,j , which measures the normalized number of pairs X u,i and X u,j , satisfying the support set condition. These steps are summarized in the Algorithm 2 (illustrated in Figure 1 ). Moreover, its complexity is dominated by the pairwise search over I K (see Step 1), which comes at

O(|I K | 2 -|I K |) ≈ O(|I K | 2 ) operations.

Algorithm 2 Causal Relations Mining (CRM)

1. Check the support set condition, via a pairwise search to check for pairs ψ i and ψ j satisfying supp(ψ j ) ⊆ supp(ψ i ), ∀(i, j) ∈ I K × I K , i = j (done over the training set K) 2. Build the adjacency matrix A, in (7) , and compute the influence score β i , ∀i ∈ I K 3. Find all pairs (i, j) such that a i,j = 1. For each of these pairs it holds that (Proposition 1), X u,i and X u,j are causally related, i.e., X u,i = 1 implies X u,j = 1, and X u,j = 1 implies X u,i = 1 (8) 4. Identify (if possible) maximally supported RVs (Corollary 1), M = {i | ψ i = 1}. For each of them, the relations in ( 8) hold for all i ∈ I K

5. SPECIAL CASES AND APPLICATIONS:

We highlight relevant special cases and applications of our approach.

Special Case: Unsupervised Learning Setting. Note that Algorithm 1 is equally applicable to an unsupervised learning task (no prediction needed). Moreover, if the training set has no missing data, i.e., K = D, then an alternate solution to (Q) may be obtained using the binary matrix factorization (MF) (Slawski et al., 2013) method. However, it is not applicable when factorizing a training set K (where K ⊂ D). Unlike Algorithm 1, binary MF does not offer prediction. Application: Analysis of Gene Expression. The ability to analyze gene expression data is critical to DNA research (Zhang et al., 2009) . This task can be formulated in our proposed framework as follows: P[X u,i = 1] denotes the probability that the genes being studied are expressed in Figure 2 : Performance of proposed method sample u ∈ U at location i ∈ I of the DNA. In this setting, the set Kolmogorov elementary events, {ω 1 , • • • , ω D }, represents the various genes that might be involved. Consequently, our approach models the probability that a set of genes is expressed as a convex mixture of D various genes. We recall that this intuition follows from our model and is not enforced explicitly. More importantly, the causal relations can identify groups of DNA locations for which the expression of a gene (or its absence) in a given location implies its presence (or its absence) for all other locations in the group. We will numerically show the usefulness of these relations to DNA analysis, in Section 7.2.

6. PRACTICAL CONSIDERATIONS:

Overfitting: Regularization parameters (to mitigate overfitting) can be included without any changes to the solution method. An 2 -regularization can be included in

(Q 1 ): f (θ u ) = θ T u (Q u + λ u I D )θ u -2θ T u r u + γ u ,

where the regularizer λ u ≥ 0 is absorbed into a "new" matrix (Q u + λ u I D ). Note that an 1 -regularization for θ u would not work, since θ u ∈ P. Similarly, an

1 -regularization for (Q 2 ) is: g(ψ i ) = ψ T i S i ψ i -2(v i -(µ i /2)1) T ψ i + γ i ,

where the regularizer µ i ≥ 0 is absorbed into the linear term, since µ i ψ i 1 = µ i 1 T ψ i , for ψ i binary. Optimality Gap: We recall that the proposed SDR method was shown to be quasi-optimal in providing approximate binary solutions to (Q 2 ). Thus, the relaxation does not affect the interpretability, in the sense that Proposition 1 and Corollary 1 still hold. While the derivations pertaining to causal relations (Section 4) assume globally optimal solutions to (Q) -an NP-hard problem, Algorithm 1 guarantees locally optimal ones. Thus, a bound on the gap between these solutions is needed. We highlight this issue as an interesting topic for further investigation.

6.1. LIMITATIONS:

Computational Complexity: The computational complexity of Algorithm 1 is dominated by the solution in (5), which is ≈ O(D 4.5 ) operations per iteration of Algorithm 1 (recalling the negligible cost of the FW method). Notice that the additional complexity compared to matrix factorization (its extensions), for which the complexity is O(D 3 ) (e.g., D ≤ 16 in all numerical results). Moreover, we are already investigating complexity reduction techniques leveraging the structure of the SDP, and distributed solutions to (P ) to enable parallelization. Finally, the computational complexity of Algorithm 2 consists mainly of step 1, which has ≈ O(|I K | 2 ) operations.

Learning a non-stationary distribution: The proposed method assumes that distributions of the RVs in the training set are stationary: Indeed, scenarios with time-varying distributions are a limitation (and interesting future directions). However, in learning it is quite common to assume that the datagenerating distribution is stationary.

7. NUMERICAL RESULTS:



7.1. APPLICATION TO RECOMMENDATION SYSTEMS:

Experimental Setup: We evaluate the performance of Algorithm 1 under the following. We opted to not have a stopping criterion based on the prediction error of the algorithm (on a validation set), since we are primarily interested in the model that we learn in the training phase. Nonetheless, as we have reported in Appendix A.7, there is not loss in the predictive performance of the model. Hereafter, "Prop." refers to the proposed method. Setting 1: An artificial training dataset, where p u,i ∈ K = {U = 20} × {I = 40}, where {p u,i } are i.i.d. and uniformly chosen on the unit interval. We benchmark against a variant on Algorithm 1, where the ψ-step for Q 2 is replaced by an exhaustive search. As the data is artificial (unsupervised learning setting), we include the Binary MF in Slawski et al. (2013) [Algorithm 2]. Setting 2: The training set K, is chosen as the MovieLens 100K, with U = 943 users and I = 1682 items, split into 80% for training and 20% for testing. Let { ψi }, { θu } the output of Algorithm 1, after 5 iterations. We benchmark against matrix factorization (MF) (Koren et al., 2009) , non-negative MF (NMF) (Lee & Seung, 2001) , SVD++ (Koren, 2008) (ensuring the dimension of the factorization, k, is close to D), non-negative models (NNM) Stark (2015) , and the K-means (K-M) algorithm. The implementation and results use the MyMediaLite package (Gantner et al., 2011) . Training Performance for Supervised Learning (Setting 2): Recall that Binary MF is not applicable here, due to the supervised learning setting. The same conclusions hold when testing Algorithm 1 on the ML100K (Figure 2b ). We underline that while NNMs yield better training performance over Algorithm 1, the latter will have better test performance (since NNMs are indeed defined by relaxing KMs). In Table 4 (Appendix A.7), we empirically verify that the test performance of the proposed method outperforms the benchmarks in Setting 2.

Training

Interpretability via Causal Relations (Setting 2): We numerically evaluate the relations of Algorithm 2. The influence score for each item in the training set, β i , is shown in Figure 2c where we displayed items with 'high' influence score, β i ≥ 0.5. Note that each of these high influence items is causally related to at least half of the items in the training set. This confirms the effectiveness of Algorithm 2 for finding causal relations, and that a sparse adjacency matrix is uncommon. We next identify the set of items corresponding to maximally supported RVs, M = {119, 814, 1188, 1190, 1290, 1393, 1462, 1486, 1494, 1530, 1590, 1638}. For each of these items, a user liking one given item, implies he/she likes all other items in the training set. Interestingly, these results remain the same when D = 24, thereby suggesting that procedure for mining causal relations is quite stable.

7.2. APPLICATION TO GENE EXPRESSION:

Experimental Setup: Following the problem statement in Section 5, we show the usefulness Algorithm 2 for gene expression. We first define another setting. Setting 3: We use the REGED0 dataset3 . Element (u, i) in the input matrix, [L] u,i represents the level of gene expression for sample u at location i of the DNA, and is reported as integer between 0 and L max . Similarly to recommendation systems (Section 2.2), the training set is obtained as p u,i := [L] u,i /L max . Thus, p u,i denotes the (empirical) probability that the gene is expressed in sample u at location i. 4 After running Algorithm 1, we use Algorithm 2 to mine causal relations. As this is unsupervised learning setting, we include the binary MF (Slawski et al., 2013) 2a where both algorithms yield similar performance). This drastic degradation in the performance of binary MF compared to the small artificial of Figure 2a may be attributed to increasing the data/problem size (though we were unable to empirically verify this claim). Thus, we opted to use the solution of Algorithm 1 as Figure 3 : Training performance of proposed method for REGED0 dataset a basis for finding causal relations. We identified 10 DNA locations corresponding to the highest β i as: {813, 250, 774, 706, 380, 49, 477, 162, 740, 702}. For instance, the highest influence score of .1612 at DNA location 813 allowed to identify a set, S, of 161 different locations which are causally related to DNA location 813. More specifically, the expression (or absence) of a given gene in location 813 implies its expression (or absence) in all other DNA locations in S. While similar relations are possible using gene analysis methods, the above causal relations follow from our rigorous mathematical framework.

We also ran experiments on the ML 1M dataset (10x larger that ML 100K) and observed that main conclusions were the same. For lack of space, we instead opted to include different datasets such as the REGED0 dataset (5x larger than ML 100K) to exemplify another application of the approach. However, we note that a large-scale implementation of the method may be challenging at this stage. As mentioned earlier, we are currently investigating complexity reduction methods before running experiments involving large datasets. Rather, the current work is intended a proof of concept of the usefulness of such an approach.

8. CONCLUSION:

We have proposed a framework for learning a Kolmogorov model, associated with a collection of binary random variables. Interpretability of the model (as defined by logical implication and causality) was harnessed by deriving causal relations, i.e., by finding sufficient conditions that bind outcomes of certain random variables. We also proposed an algorithm for computing a Kolmogorov model, a combinatorial non-convex problem, and showed its convergence to a stationary point of the problem using results from block-coordinate descent. The combinatorial nature of the problem was addressed using a semi-definite relaxation. We also proposed an efficient algorithm to mine for the causal relations inherent to our model. Our results suggest that increased interpretability and improved prediction, do not cause a significant increase in complexity. We highlight several key issues for future work, e.g., complexity reduction for (5) (leveraging that its dual is piece-wise linear), and a sufficient condition for identifiability (by adapting that of Fu et al. (2017)) 

:

Hereafter, for the sake of simplicity, causal relations mean logical relations within the data.

The method for acquiring (estimates of) the empirical probabilities is application-dependent (see Section 5).

http://www.causality.inf.ethz.ch/data/REGED.html

Finding better/different mappings from [L]u,i to pu,i is an interesting direction, yet outside the scope of this our work.

Figures and Tables:

Figure 1: Algorithm for CRM (toy example)

Performance for Unsupervised Learning (Setting 1): Fig 2a shows the resulting normalized training RMSE= ( (u,i)∈K |p u,i -θT u ψi | 2 /|K|) 1/2. We observe that the monotone convergence in Lemma 2 is validated numerically, and that the training error decays with increasing model size, D. Note that the training performance of Algorithm 1 is indistinguishable from its exhaustive search variant. Moreover, Algorithm 1 converges to solution whose performance is similar to Binary MF, with a few iterations (except for D = 8 where Algorithm 1 outperforms Binary MF).

method. Results: While Fig 3a shows the training performance for several values of D, Figure 3b plots the corresponding influence score that is obtained with our method. Fig 3a reveals a huge gap (≈ 3× less) between the training error of Algorithm 1 and Binary MF (unlike Figure

Influence score for each DNA location (Setting 3)

Let p u,i denote the empirical value of P[X u,i = 1]. We assume that {p u,i } are known for elements of a training set K ⊆ D, where K

Table content: <table><tr><td>(formalized in Appendix A.2). Notice that the</td></tr><tr><td>model in (1) can approximate with arbitrarily small accuracy the measure corresponding to P[•] for a</td></tr><tr><td>large enough D.</td></tr><tr><td>Problem 1 (Problem Statement)</td></tr></table>

Peer review (iclr2019):
 The paper deals with an interesting problem. The presentation is clear the the approach intuitive. However, the reviewer has some concerns about the pertinence of the approach and the relationship with related work. It would be very helpful if the authors could contrast and compare the proposed approach (both qualitatively and quantitatively in their numerical experiments) with methods for sparse non-negative matrix factorization. These would also lend themselves to causal interpretation. The need for modeling probabilities p_u,i in the key motivating applications is questionable. Indeed in both recommendation systems and gene expression datasets the observations are not readily in the form of probabilities. For instance in the experiments, the authors normalize the level of gene expressions to make is look as a probability, which by the way is very different from the i.id uniform setup considered in setting 1. For the movielens dataset, it is unclear how the data was preprocessed to obtain the observed p_ui.