Performative predictions are forecasts which influence the outcomes they aim to predict, undermining the existence of correct forecasts and standard methods of elicitation and estimation. We show that conditioning forecasts on covariates that separate them from the outcome renders the target distribution forecast-invariant, guaranteeing well-posedness of the forecasting problem. However, even under this condition, classical proper scoring rules fail to elicit correct forecasts. We prove a general impossibility result and identify two solutions: (i) in decision-theoretic settings, elicitation of correct and incentive-compatible forecasts is possible if forecasts are separating; (ii) scoring with unbiased estimates of the divergence between the forecast and the induced distribution of the target variable yields correct forecasts. Applying these insights to parameter estimation, conditional forecasts and proper scoring rules enable performatively stable estimation of performatively correct parameters, resolving the issues raised by Perdomo et al. (2020). Our results expose fundamental limits of classical forecast evaluation and offer new tools for reliable and accurate forecasting in performative settings.

Reasoning about fairness through correlation-based notions is rife with pitfalls. The 1973 University of California, Berkeley graduate school admissions case from \citet{BickelHO75} is a classic example of one such pitfall, namely Simpson’s paradox. The discrepancy in admission rates among male and female applicants, in the aggregate data over all departments, vanishes when admission rates per department are examined. We reason about the Berkeley graduate school admissions case through a causal lens. In the process, we introduce a statistical test for causal hypothesis testing based on Pearl’s instrumental-variable inequalities \citep{Pearl95}. We compare different causal notions of fairness that are based on graphical, counterfactual and interventional queries on the causal model, and develop statistical tests for these notions that use only observational data. We study the logical relations between notions, and show that while notions may not be equivalent, their corresponding statistical tests coincide for the case at hand. We believe that a thorough case-based causal analysis helps develop a more principled understanding of both causal hypothesis testing and fairness.

Maximal Ancestral Graphs (MAGs) provide an abstract representation of Directed Acyclic Graphs (DAGs) with latent (selection) variables. These graphical objects encode information about ancestral relations and d-separations of the DAGs they represent. This abstract representation has been used amongst others to prove the soundness and completeness of the FCI algorithm for causal discovery, and to derive a do-calculus for its output. One significant inherent limitation of MAGs is that they rule out the possibility of cyclic causal relationships. In this work, we address that limitation. We introduce and study a class of graphical objects that we coin "$\sigma$-Maximal Ancestral Graphs" ("$\sigma$-MAGs"). We show how these graphs provide an abstract representation of (possibly cyclic) Directed Graphs (DGs) with latent (selection) variables, analogously to how MAGs represent DAGs. We study the properties of these objects and provide a characterization of their Markov equivalence classes.

When estimating a regression model, we might have data where some labels are missing, or our data might be biased by a selection mechanism. When the response or selection mechanism is ignorable (i.e., independent of the response variable given the features) one can use off-the-shelf regression methods; in the nonignorable case one typically has to adjust for bias. We observe that privileged information (i.e. information that is only available during training) might render a nonignorable selection mechanism ignorable, and we refer to this scenario as Privilegedly Missing at Random (PMAR). We propose a novel imputation-based regression method, named repeated regression, that is suitable for PMAR. We also consider an importance weighted regression method, and a doubly robust combination of the two. The proposed methods are easy to implement with most popular out-of-the-box regression algorithms. We empirically assess the performance of the proposed methods with extensive simulated experiments and on a synthetically augmented real-world dataset. We conclude that repeated regression can appropriately correct for bias, and can have considerable advantage over weighted regression, especially when extrapolating to regions of the feature space where response is never observed.

We present a new, efficient procedure to establish Markov equivalence between directed graphs that may or may not contain cycles. It is based on the Cyclic Equivalence Theorem (CET) in the seminal works on cyclic models by Thomas Richardson in the mid ’90s, but now rephrased from an ancestral perspective. The resulting characterization leads to a procedure for establishing Markov equivalence between graphs that no longer requires explicit tests for $d$-separation, leading to a significantly reduced algorithmic complexity. The conceptually simplified characterization may help to reinvigorate theoretical research towards sound and complete cyclic discovery in the presence of latent confounders.

Mathematical models of the real world are simplified representations of complex systems. A caveat to using mathematical models is that predicted causal effects and conditional independences may not be robust under model extensions, limiting applicability of such models. In this work, we consider conditions under which qualitative model predictions are preserved when two models are combined. Under mild assumptions, we show how to use the technique of causal ordering to efficiently assess the robustness of qualitative model predictions. We also characterize a large class of model extensions that preserve qualitative model predictions. For dynamical systems at equilibrium, we demonstrate how novel insights help to select appropriate model extensions and to reason about the presence of feedback loops. We illustrate our ideas with a viral infection model with immune responses.

For a continuous random variable $Z$, testing conditional independence $X\indep Y|Z$ is known to be a particularly hard problem. It constitutes a key ingredient of many constraint-based causal discovery algorithms. These algorithms are often applied to datasets containing binary variables, which indicate the ‘context’ of the observations, e.g. a control or treatment group within an experiment. In these settings, conditional independence testing with $X$ or $Y$ binary (and the other continuous) is paramount to the performance of the causal discovery algorithm. To our knowledge no nonparametric ‘mixed’ conditional independence test currently exists, and in practice tests that assume all variables to be continuous are used instead. In this paper we aim to fill this gap, as we combine elements of Holmes et al. (2015) and Teymur and Filippi (2020) to propose a novel Bayesian nonparametric conditional two-sample test. Applied to the Local Causal Discovery algorithm, we investigate its performance on both synthetic and real-world data, and compare with state-of-the-art conditional independence tests.

One of the core assumptions in causal discovery is the faithfulness assumption—i.e. assuming that independencies found in the data are due to separations in the true causal graph. This assumption can, however, be violated in many ways, including xor connections, deterministic functions or cancelling paths. In this work, we propose a weaker assumption that we call 2-adjacency faithfulness. In contrast to adjacency faithfulness, which assumes that there is no conditional independence between each pair of variables that are connected in the causal graph, we only require no conditional independence between a node and a subset of its Markov blanket that can contain up to two nodes. Equivalently, we adapt orientation faithfulness to this setting. We further propose a sound orientation rule for causal discovery that applies under weaker assumptions. As a proof of concept, we derive a modified Grow and Shrink algorithm that recovers the Markov blanket of a target node and prove its correctness under strictly weaker assumptions than the standard faithfulness assumption.

Real-world complex systems are often modelled by sets of equations with endogenous and exogenous variables. What can we say about the causal and probabilistic aspects of variables that appear in these equations without explicitly solving the equations? We make use of Simon's causal ordering algorithm (Simon, 1953) to construct a causal ordering graph and prove that it expresses the effects of soft and perfect interventions on the equations under certain unique solvability assumptions. We further construct a Markov ordering graph and prove that it encodes conditional independences in the distribution implied by the equations with independent random exogenous variables, under a similar unique solvability assumption. We discuss how this approach reveals and addresses some of the limitations of existing causal modelling frameworks, such as causal Bayesian networks and structural causal models.

While feedback loops are known to play important roles in many complex systems, their existence is ignored in a large part of the causal discovery literature, as systems are typically assumed to be acyclic from the outset. When applying causal discovery algorithms designed for the acyclic setting on data generated by a system that involves feedback, one would not expect to obtain correct results. In this work, we show that—surprisingly—the output of the Fast Causal Inference (FCI) algorithm is correct if it is applied to observational data generated by a system that involves feedback. More specifically, we prove that for observational data generated by a simple and sigma-faithful Structural Causal Model (SCM), FCI is sound and complete, and can be used to consistently estimate (i) the presence and absence of causal relations, (ii) the presence and absence of direct causal relations, (iii) the absence of confounders, and (iv) the absence of specific cycles in the causal graph of the SCM. We extend these results to constraint-based causal discovery algorithms that exploit certain forms of background knowledge, including the causally sufficient setting (e.g., the PC algorithm) and the Joint Causal Inference setting (e.g., the FCI-JCI algorithm).

The gold standard for discovering causal relations is by means of experimentation. Over the last decades, alternative methods have been proposed that can infer causal relations between variables from certain statistical patterns in purely observational data. We introduce Joint Causal Inference (JCI), a novel approach to causal discovery from multiple data sets from different contexts that elegantly unifies both approaches. JCI is a causal modeling framework rather than a specific algorithm, and it can be implemented using any causal discovery algorithm that can take into account certain background knowledge. JCI can deal with different types of interventions (e.g., perfect, imperfect, stochastic, etc.) in a unified fashion, and does not require knowledge of intervention targets or types in case of interventional data. We explain how several well-known causal discovery algorithms can be seen as addressing special cases of the JCI framework, and we also propose novel implementations that extend existing causal discovery methods for purely observational data to the JCI setting. We evaluate different JCI implementations on synthetic data and on flow cytometry protein expression data and conclude that JCI implementations can considerably outperform state-of-the-art causal discovery algorithms.

We study the performance of Local Causal Discovery (LCD) [5], a simple and efficient constraint-based method for causal discovery, in predicting causal effects in large-scale gene expression data. We construct practical estimators specific to the high-dimensional regime. Inspired by the ICP algorithm [13], we use an optional preselection method and two different statistical tests. Empirically, the resulting LCD estimator is seen to closely approach the accuracy of ICP, the state-of-the-art method, while it is algorithmically simpler and computationally more efficient.

An important goal common to domain adaptation and causal inference is to make accurate predictions when the distributions for the source (or training) domain(s) and target (or test) domain(s) differ. In many cases, these different distributions can be modeled as different contexts of a single underlying system, in which each distribution corresponds to a different perturbation of the system, or in causal terms, an intervention. We focus on a class of such causal domain adaptation problems, where data for one or more source domains are given, and the task is to predict the distribution of a certain target variable from measurements of other variables in one or more target domains. We propose an approach for solving these problems that exploits causal inference and does not rely on prior knowledge of the causal graph, the type of interventions or the intervention targets. We demonstrate our approach by evaluating a possible implementation on simulated and real world data.

Learning individual-level causal effects from observational data, such as inferring the most effective medication for a specific patient, is a problem of growing importance for policy makers. The most important aspect of inferring causal effects from observational data is the handling of confounders, factors that affect both an intervention and its outcome. A carefully designed observational study attempts to measure all important confounders. However, even if one does not have direct access to all confounders, there may exist noisy and uncertain measurement of proxies for confounders. We build on recent advances in latent variable modeling to simultaneously estimate the unknown latent space summarizing the confounders and the causal effect. Our method is based on Variational Autoencoders (VAE) which follow the causal structure of inference with proxies. We show our method is significantly more robust than existing methods, and matches the state-of-the-art on previous benchmarks focused on individual treatment effects.

Constraint-based causal discovery from limited data is a notoriously difficult challenge due to the many borderline independence test decisions. Several approaches to improve the reliability of the predictions by exploiting redundancy in the independence information have been proposed recently. Though promising, existing approaches can still be greatly improved in terms of accuracy and scalability. We present a novel method that reduces the combinatorial explosion of the search space by using a more coarse-grained representation of causal information, drastically reducing computation time. Additionally, we propose a method to score causal predictions based on their confidence. Crucially, our implementation also allows one to easily combine observational and interventional data and to incorporate various types of available background knowledge. We prove soundness and asymptotic consistency of our method and demonstrate that it can outperform the state-of-the-art on synthetic data, achieving a speedup of several orders of magnitude. We illustrate its practical feasibility by applying it on a challenging protein data set.

The discovery of causal relationships from purely observational data is a fundamental problem in science. The most elementary form of such a causal discovery problem is to decide whether $X$ causes $Y$ or, alternatively, $Y$ causes $X$, given joint observations of two variables $X,Y$. An example is to decide whether altitude causes temperature, or vice versa, given only joint measurements of both variables. Even under the simplifying assumptions of no confounding, no feedback loops, and no selection bias, such bivariate causal discovery problems are challenging. Nevertheless, several approaches for addressing those problems have been proposed in recent years. We review two families of such methods: methods based on Additive Noise Models (ANMs) and Information Geometric Causal Inference (IGCI). We present the benchmark CauseEffectPairs that consists of data for 100 different cause-effect pairs selected from 37 data sets from various domains (e.g., meteorology, biology, medicine, engineering, economy, etc.) and motivate our decisions regarding the ground truth causal directions of all pairs. We evaluate the performance of several bivariate causal discovery methods on these real-world benchmark data and in addition on artificially simulated data. Our empirical results on real-world data indicate that certain methods are indeed able to distinguish cause from effect using only purely observational data, although more benchmark data would be needed to obtain statistically significant conclusions. One of the best performing methods overall is the method based on Additive Noise Models that has originally been proposed by Hoyer et al. (2009), which obtains an accuracy of 63 $\pm$ 10 % and an AUC of 0.74 $\pm$ 0.05 on the real-world benchmark. As the main theoretical contribution of this work we prove the consistency of that method.

By aggregating data for complex traits in a biologically meaningful way, gene and gene-set analysis constitute a valuable addition to single-marker analysis. However, although various methods for gene and gene-set analysis currently exist, they generally suffer from a number of issues. Statistical power for most methods is strongly affected by linkage disequilibrium between markers, multi-marker associations are often hard to detect, and the reliance on permutation to compute p-values tends to make the analysis computationally very expensive. To address these issues we have developed MAGMA, a novel tool for gene and gene-set analysis. The gene analysis is based on a multiple regression model, to provide better statistical performance. The gene-set analysis is built as a separate layer around the gene analysis for additional flexibility. This gene-set analysis also uses a regression structure to allow generalization to analysis of continuous properties of genes and simultaneous analysis of multiple gene sets and other gene properties. Simulations and an analysis of Crohn's Disease data are used to evaluate the performance of MAGMA and to compare it to a number of other gene and gene-set analysis tools. The results show that MAGMA has significantly more power than other tools for both the gene and the gene-set analysis, identifying more genes and gene sets associated with Crohn's Disease while maintaining a correct type 1 error rate. Moreover, the MAGMA analysis of the Crohn's Disease data was found to be considerably faster as well.

We study a particular class of cyclic causal models, where each variable is a (possibly nonlinear) function of its parents and additive noise. We prove that the causal graph of such models is generically identifiable in the bivariate, Gaussian-noise case. We also propose a method to learn such models from observational data. In the acyclic case, the method reduces to ordinary regression, but in the more challenging cyclic case, an additional term arises in the loss function, which makes it a special case of nonlinear independent component analysis. We illustrate the proposed method on synthetic data.

Inference in matrix-variate Gaussian models has major applications for multi- output prediction and joint learning of row and column covariances from matrix- variate data. Here, we discuss an approach for efficient inference in such models that explicitly account for iid observation noise. Computational tractability can be retained by exploiting the Kronecker product between row and column covariance matrices. Using this framework, we show how to generalize the Graphical Lasso in order to learn a sparse inverse covariance between features while accounting for a low-rank confounding covariance between samples. We show practical utility on applications to biology, where we model covariances with more than 100,000 di- mensions. We find greater accuracy in recovering biological network structures and are able to better reconstruct the confounders.

We present a graphical model framework for decoding in the visual ERP-based speller system. The proposed framework allows researchers to build generative models from which the decoding rules are obtained in a straightforward manner. We suggest two models for generating brain signals conditioned on the stimulus events. Both models incorporate letter frequency information but assume different dependencies between brain signals and stimulus events. For both models, we derive decoding rules and perform a discriminative training. We show on real visual speller data how decoding performance improves by incorporating letter frequency information and using a more realistic graphical model for the dependencies between the brain signals and the stimulus events. Furthermore, we discuss how the standard approach to decoding can be seen as a special case of the graphical model framework. The letter also gives more insight into the discriminative approach for decoding in the visual speller system.

We propose a novel method for inferring whether X causes Y or vice versa from joint observations of X and Y. The basic idea is to model the observed data using probabilistic latent variable models, which incorporate the effects of unobserved noise. To this end, we consider the hypothetical effect variable to be a function of the hypothetical cause variable and an independent noise term (not necessarily additive). An important novel aspect of our work is that we do not restrict the model class, but instead put general non-parametric priors on this function and on the distribution of the cause. The causal direction can then be inferred by using standard Bayesian model selection. We evaluate our approach on synthetic data and real-world data and report encouraging results.

This letter introduces a nonlinear measure of independence between random variables for remote sensing supervised feature selection. The so-called Hilbert–Schmidt independence criterion (HSIC) is a kernel method for evaluating statistical dependence and it is based on computing the Hilbert–Schmidt norm of the cross-covariance operator of mapped samples in the corresponding Hilbert spaces. The HSIC empirical estimator is easy to compute and has good theoretical and practical properties. Rather than using this estimate for maximizing the dependence between the selected features and the class labels, we propose the more sensitive criterion of minimizing the associated HSIC$p$-value. Results in multispectral, hyperspectral, and SAR data feature selection for classification show the good performance of the proposed approach.

Motivated by causal inference problems, we propose a novel method for regression that minimizes the statistical dependence between regressors and residuals. The key advantage of this approach to regression is that it does not assume a particular distribution of the noise, i.e., it is non-parametric with respect to the noise distribution. We argue that the proposed regression method is well suited to the task of causal inference in additive noise models. A practical disadvantage is that the resulting optimization problem is generally non-convex and can be difficult to solve. Nevertheless, we report good results on one of the tasks of the NIPS 2008 Causality Challenge, where the goal is to distinguish causes from effects in pairs of statistically dependent variables. In addition, we propose an algorithm for efficiently inferring causal models from observational data for more than two variables. The required number of regressions and independence tests is quadratic in the number of variables, which is a significant improvement over the simple method that tests all possible DAGs.

The discovery of causal relationships between a set of observed variables is a fundamental problem in science. For continuous-valued data linear acyclic causal models are often used because these models are well understood and there are well-known methods to fit them to data. In reality, of course, many causal relationships are more or less nonlinear, raising some doubts as to the applicability and usefulness of purely linear methods. In this contribution we show that in fact the basic linear framework can be generalized to nonlinear models with additive noise. In this extended framework, nonlinearities in the data-generating process are in fact a blessing rather than a curse, as they typically provide information on the underlying causal system and allow more aspects of the true data-generating mechanisms to be identified. In addition to theoretical results we show simulations and some simple real data experiments illustrating the identification power provided by nonlinearities.

We propose a novel bound on single-variable marginal probability distributions in factor graphs with discrete variables. The bound is obtained by propagating bounds (convex sets of probability distributions) over a subtree of the factor graph, rooted in the variable of interest. By construction, the method not only bounds the exact marginal probability distribution of a variable, but also its approximate Belief Propagation marginal (``belief''). Thus, apart from providing a practical means to calculate bounds on marginals, our contribution also lies in providing a better understanding of the error made by Belief Propagation. We show that our bound outperforms the state-of-the-art on some inference problems arising in medical diagnosis.

Novel conditions are derived that guarantee convergence of the Sum-Product Algorithm (also known as Loopy Belief Propagation or simply Belief Propagation (BP)) to a unique fixed point, irrespective of the initial messages, for parallel (synchronous) updates. The computational complexity of the conditions is polynomial in the number of variables. In contrast with previously existing conditions, our results are directly applicable to arbitrary factor graphs (with discrete variables) and are shown to be valid also in the case of factors containing zeros, under some additional conditions. The conditions are compared with existing ones, numerically and, if possible, analytically. For binary variables with pairwise interactions, sufficient conditions are derived that take into account local evidence (i.e., single-variable factors) and the type of pair interactions (attractive or repulsive). It is shown empirically that this bound outperforms existing bounds.

We introduce a computationally efficient method to estimate the valid- ity of the BP method as a function of graph topology, the connectiv- ity strength, frustration and network size. We present numerical results that demonstrate the correctness of our estimates for the uniform random model and for a real-world network ("C. Elegans"). Although the method is restricted to pair-wise interactions, no local evidence (zero "biases") and binary variables, we believe that its predictions correctly capture the limitations of BP for inference and MAP estimation on arbitrary graphi- cal models. Using this approach, we find that BP always performs better than MF. Especially for large networks with broad degree distributions (such as scale-free networks) BP turns out to significantly outperform MF. 1 Introduction Loopy Belief Propagation (BP) [1] and its generalizations (such as the Cluster Variation Method [2]) are powerful methods for inference and optimization. As is well-known, BP is exact on trees, but also yields surprisingly good results for many other graphs that arise in real-world applications [3, 4]. On the other hand, for densely connected graphs with high interaction strengths the results can be quite bad or BP can simply fail to converge. Despite the fact that BP is often used in applications nowadays, a good theoretical understanding of its convergence properties and the quality of the approximation is still lacking (except for the very special case of graphs with a single loop [5]). In this article we attempt to answer the question in what way the quality of the BP re- sults depends on the topology of the underlying graph (looking at structural properties such as short cycles and large "hubs") and on the interaction potentials (i.e. strength and frus- tration). We do this for the special but interesting case of binary networks with symmetric pairwise potentials (i.e. Boltzmann machines) without local evidence. This has the practical advantage that analytical calculations are feasible and furthermore we believe that adding local evidence will only serve to extend the domain of convergence, implying this to be the worst-case scenario. We compare the results with those of the variational mean-field (MF) method. Real-world graphs are often far from uniformly random and possess structure such as clus- tering and power-law degree distributions [6]. Since we expect these structural features to arise in many applications of BP, we focus in this article on graphs modeling this kind of features. In particular, we consider Erdos-Renyi uniform random graphs [7], Barabasi- Albert "scale-free" graphs [8], and the neural network of a widely studied worm, the Caenorhabditis elegans. This paper is organized as follows. In the next section we describe the class of graphical models under investigation and explain our method to efficiently estimate the validity of BP and MF. In section 3 we give a qualitative discussion of how the connectivity strength and frustration generally govern the model behavior and discuss the relevant regimes of the model parameters. We show for uniform random graphs that our validity estimates are in very good agreement with the real behavior of the BP algorithm. In section 4 we study the influence of graph topology. Thanks to the numerical efficiency of our estimation method we are able to study very large (N 10000) networks, for which it would not be feasible to simply run BP and look what happens. We also try our method on the neural network of the worm C. Elegans and find almost perfect agreement of our predictions with observed BP behavior. We conclude that BP is always better than MF and that the difference is particularly striking for the case of large networks with broad degree distributions such as scale-free graphs. 2 Model, paramagnetic solution and stability analysis Let G = (V, B) be an undirected labelled graph without self-connections, defined by a set of nodes V = {1, . . . , N } and a set of links B {(i, j) | 1 i < j N }. The adjacency matrix corresponding to G is denoted M and defined as follows: Mij := 1 if (ij) B or (ji) B and 0 otherwise. We denote the set of neighbors of node i V by Ni := {j V | (ij) B} and its degree by di := #(Ni). We define the average degree d := 1 d N iV i and the maximum degree := maxiV di. To each node i we associate a binary random variable xi taking values in {-1, +1}. Let W be a symmetric N N -matrix defining the strength of the links between the nodes. The probability distribution over configurations x = (x1, . . . , xN ) is given by 1 1 1 M P(x) := eWijxixj = e 2 ijWijxixj (1) Z Z (ij)B i,jV with Z a normalization constant. We will take the weight matrix W to be random, with i.i.d. entries {Wij}1i For this model, instead of using the single-node and pair-wise beliefs bi(xi) resp. bij(xi, xj), it turns out to be more convenient to use the (equivalent) quantities m := {mi}iV and := {ij}(ij)B, defined by: mi := bi(+1) - bi(-1); ij := bij(+1, +1) - bij(+1, -1) - bij(-1, +1) + bij(-1, -1). We will use these throughout this paper. We call the mi magnetizations; note that the expectation values E xi vanish because of the symmetry in the probability distribution (1). As is well-known [2, 9], fixed points of BP correspond to stationary points of the Bethe free energy, which is in this case given by N 1 + mixi FBe(m, ) := - Wijij + (1 - di) 2 (ij)B i=1 xi=1 1 + mixi + mjxj + xixjij + 4 (ij)B xi,xj =1 with (x) := x log x. Note that with this parameterization all normalization and overlap constraints (i.e. b x ij (xi, xj ) = bi(xi)) are satisfied by construction [10]. We can mini- j mize the Bethe free energy analytically by setting its derivatives to zero; one then immedi- ately sees that a possible solution of the resulting equations is the paramagnetic1 solution: mi = 0 and ij = tanh Wij (for (ij) B). For this solution to be a minimum (instead of a saddle point or maximum), the Hessian of FBe at that point should be positive-definite. This condition turns out to be equivalent to the following Bethe stability matrix 2 ij (A ik Be)ij := ij 1 + - Mij (with ij = tanh Wij) (2) 1 - 2 1 - 2 kN ik ij i being positive-definite. Whether this is the case obviously depends on the values of the weights Wij and the adjacency matrix M . Since for zero weights (W = 0), the stability matrix is just the identity matrix, the paramagnetic solution is a minimum of the Bethe free energy for small values of the weights Wij. The question of what "small" exactly means in terms of J and J0 and how this relates to the graph topology will be taken on in the next two sections. First we discuss the situation for the mean-field variational method. The mean-field free energy FMF (m) only depends on m; we can set its derivatives to zero, which again yields the paramagnetic solution m = 0. The corresponding stability matrix (equal to the Hes- sian) is given by (AMF )ij := ij - WijMij and should be positive-definite for the paramagnetic solution to be stable. One can prove [11] that ABe is positive-definite whenever AMF is positive-definite. Since the exact mag- netizations are zero, we conclude that the Bethe approximation is better than the mean-field approximation for all possible choices of the weights W . As we will see later on, this dif- ference can become quite large for large networks. 3 Weight dependence The behavior of the graphical model depends critically on the parameters J0 and J. Taking the graph topology to be uniformly random (see also subsection 4.1) we recover the model known in the statistical physics community as the Viana-Bray model [12], which has been thoroughly studied and is quite well-understood. In the limit N , there are different relevant regimes ("phases") for the parameters J and J0 to be distinguished (cf. Fig. 1): The paramagnetic phase, where the magnetizations all vanish (m = 0), valid for J and J0 both small. The ferromagnetic phase, where two configurations (characterized by all magne- tizations being either positive or negative) each get half of the probability mass. This is the phase occurring for large J0. 1Throughout this article, we will use terminology from statistical physics if there is no good corresponding terminology in the field of machine learning available. BP convergence behavior Stability m=0 minimum Bethe free energy 0.4 0.4 m=0 stable (spin-glass phase) no convergence ? 0.3 0.3 marginal instability J 0.2 J 0.2 convergence m=0 stable m=0 instable convergence to ferromagnetic (paramagnetic 0.1 0.1 (ferromagnetic to m=0 solutions phase) phase) 0 0 0 0.02 0.04 0.06 0.08 0.1 0 0.02 0.04 0.06 0.08 0.1 J J (a) 0 (b) 0 Figure 1: Empirical regime boundaries for the ER graph model with N = 100 and d = 20, averaged over three instances; expectation values are shown as thick black lines, standard- deviations are indicated by the gray areas. See the main text for additional explanation. The exact location of the boundary between the spin-glass and ferromagnetic phase in the right-hand plot (indicated by the dashed line) was not calculated. The red dash-dotted line shows the stability boundary for MF. The spin-glass phase where the probability mass is distributed over exponentially (in N ) many different configurations. This phase occurs for frustrated weights, i.e. for large J . Consider now the right-hand plot in Fig. 1. Here we have plotted the different regimes con- cerning the stability of the paramagnetic solution of the Bethe approximation.2 We find that the m = 0 solution is indeed stable for J and J0 small and becomes unstable at some point when J0 increases. This signals the paramagnetic-ferromagnetic phase transition. The lo- cation is in good agreement with the known phase boundary found for the N limit by advanced statistical physics methods as we show in more detail in [11]. For comparison we have also plotted the stability boundary for MF (the red dash-dotted line). Clearly, the mean-field approximation breaks down much earlier than the Bethe approximation and is unable to capture the phase transitions occurring for large connectivity strengths. The boundary between the spin-glass phase and the paramagnetic phase is more subtle. What happens is that the Bethe stability matrix becomes marginally stable at some point when we increase J , i.e. the minimum eigenvalue of ABe approaches zero (in the limit N ). This means that the Bethe free energy becomes very flat at that point. If we go on increasing J , the m = 0 solution becomes stable again (in other words, the minimum eigenvalue of the stability matrix ABe becomes positive again). We interpret the marginal instability as signalling the onset of the spin-glass phase. Indeed it coincides with the known phase boundary for the Viana-Bray model [11, 12]. We observe a similar marginal instability for other graph topologies. Now consider the left-hand plot, Fig. 1(a). It shows the convergence behavior of the BP al- gorithm, which was determined by running BP with a fixed number of maximum iterations and slight damping. The messages were initialized randomly. We find different regimes that are separated by the boundaries shown in the plot. For small J and J0, BP converges to m = 0. For J0 large enough, BP converges to one of the two ferromagnetic solutions 2Although in Fig. 1 we show only one particular graph topology, the general appearance of these plots does not differ much for other graph topologies, especially for large N . The scale of the plots mostly depends on the network size N and the average degree d as we will show in the next section. Mean Field Bethe 2 2 1.5 1.5 1/2 1 d 1 J c 0.5 0.5 0 0 10 100 1000 10000 10 100 1000 10000 N N Figure 2: Critical values for Bethe and MF for different graph topologies ( : ER, : BA) in the dense limit with d = 0.1N as a function of network size. Note that the y-axis is rescaled by d. (which one is determined by the random initial conditions). For large J , BP does not con- verge within 1000 iterations, indicating a complex probability distribution. The boundaries coincide within statistical precision with those in the right-hand plot which were obtained by the stability analysis. The computation time necessary for producing a plot such as Fig. 1(a), showing the conver- gence behavior of BP, quickly increases with increasing N . The computation time needed for the stability analysis (Fig. 1(b)), which amounts to calculating the minimal eigenvalue of the N N stability matrix, is much less, allowing us to investigate the behavior of BP for large networks. 4 Graph topology In this section we will concentrate on the frustrated case, more precisely on the case J0 = 0 (i.e. the y-axis in the regime diagrams) and study the location of the Bethe marginal instability and of the MF instability for various graph topologies as a function of network size N and average degree d. We will denote by J Be c the critical value of J at which the Bethe paramagnetic solution becomes marginally unstable and we will refer to this as the Bethe critical value. The critical value of J where the MF solution becomes unstable will be denoted as J MF c and referred to as the MF critical value. In studying the influence of graph topology for large networks, we have to distinguish two cases, which we call the dense and sparse limits. In the dense limit, we let N and scale the average degree as d = cN for some fixed constant c. In this limit, we find that the influence of the graph topology is almost negligible. For all graph topologies that we have considered, we find the following asymptotic behavior for the critical values: 1 1 J Be , J MF c c d 2 d The constant of proportionality is approximately 1. These results are illustrated in Fig. 2 for two different graph topologies that will be discussed in more detail below. In the sparse limit, we let N but keep d fixed. In that case the resulting critical values show significant dependence on the graph topology as we will see. 4.1 Uniform random graphs (ER) The first and most elementary random graph model we will consider was introduced and studied by Erdos and Renyi [7]. The ensemble, which we denote as ER(N, p), consists of 0.5 Bethe J 0.4 c 1/2 1/d 0.3 J c MF Jc 0.2 1/(21/2) 0.1 0 10 100 1000 10000 N Figure 3: Critical values for Bethe and MF for Erdos-Renyi uniform random graphs with average degree d = 10. the graphs with N nodes; links are added between each pair of nodes independently with probability p. The resulting graphs have a degree distribution that is approximately Poisson for large N and the expected average degree is E d = p(N - 1). As was mentioned before, the resulting graphical model is known in the statistical physics literature as the Viana-Bray model (with zero "external field"). Fig. 3 shows the results for the sparse limit, where p is chosen such that the expected aver- age degree is fixed to d = 10. The Bethe critical value J Be c appears to be independent of network size and is slightly larger than 1/ d. The MF critical value J MF c does depend on network size (it looks to be proportional to 1/ instead of 1/ d); in fact it can be proven that it converges very slowly to 0 as N [11], implying that the MF approximation breaks down for very large ER networks in the sparse limit. Although this is an interesting result, one could say that for all practical purposes the MF critical value J MF c is nearly independent of network size N for uniform random graphs. 4.2 Scale-free graphs (BA) A phenomenon often observed in real-world networks is that the degree distribution be- haves like a power-law, i.e. the number of nodes with degree is proportional to - for some > 0. These graphs are also known as "scale-free" graphs. The first random graph model exhibiting this behavior is from Barabasi and Albert [8]. We will consider a slightly different model, which we will denote by BA(N, m). It is defined as a stochastic process, yielding graphs with more and more nodes as time goes on. At t = 0 one starts with the graph consisting of m nodes and no links. At each time step, one node is added; it is connected with m different already existing nodes, attaching preferably to nodes with higher degree ("rich get richer"). More specifically, we take the probability to connect to a node of degree to be proportional to + 1. The degree dis- tribution turns out to have a power-law dependence for N with exponent = 3. In Fig. 4 we illustrate some BA graphs. The difference between the maximum degree and the average degree d is rather large: whereas the average degree d converges to 2m, the maximum degree is known to scale as N . Fig. 5 shows the results of the stability analysis for BA graphs with average degree d = 10. Note that the y-axis is rescaled by to show that the MF critical value J MF c is proportional to 1/ . The Bethe critical values are seen to have a scaling behavior that lies somewhere between 1/ d and 1/ . Compared to the situation for uniform ER graphs, BP now even more significantly outperforms MF. The relatively low sensitivity to the maximum degree that BP exhibits here can be understood intuitively since BA graphs resemble forests of sparsely interconnected stars of high degree, on which BP is exact.