1 Introduction
Probabilistic models play a central role in modern machine learning. They offer a powerful framework for learning from data, and they have found applications in a variety of scientific fields beyond machine learning. A longstanding goal in machine learning and statistics is to achieve a separation between modeling and inference, so that users of these tools may focus on specifying models without having to implement new inference algorithms every time the models change. Recently, work in probabilistic programming has taken up this challenge, seeking to unify probabilistic modeling with computer programming in order to dramatically increase the effectiveness of machine learning experts
(DARPA, 2013) and to equip nonexperts with effective tools for specifying models and performing inference. We anticipate that continued success toward these goals will decrease the reliance of machine learning practitioners on triedandtrue models and will shift the community toward a paradigm grounded in flexible tools for rapidly prototyping and designing new models (Bishop, 2013).There are many benefits to this increased flexibility. Scientists will be able to swap one model for another without rewriting the inference algorithms. People will be able to use domainspecific knowledge to build highlyexpressive models for their data instead of engineering complex features to pass into a generic algorithm.
In addition to these benefits, the proliferation of new models introduces some serious challenges. Foremost perhaps, is the fact that the available inference algorithms do not work well in all models. One common cause of failure, which is the focus of this work, is symmetry (nonidentifiability) in the parameterization of a model. Model symmetries, which we define formally in Section 2, can greatly affect the quality of inference and the interpretability of the output. We will elaborate on the problems caused by symmetries in Section 3. Symmetries have been identified in common models used in fields from machine learning to econometrics to political science (see Section 5), and mechanisms for breaking these symmetries have been discussed in depth. These mechanisms depend on the specific type of symmetry under consideration, but common techniques include reparameterizing the problem or imposing additional constraints on the parameter values.
If we expect machine learning to be done by experts using wellunderstood models, then it is feasible to maintain a catalog of the symmetries that occur in each common model. If, on the other hand, we anticipate a wide range of people designing completely original models, then it is important to automatically detect the presence of model symmetries so that practitioners can be informed of the problem or so that the model symmetries can automatically be broken.
To our knowledge, with the exception of some work specific on permutation symmetries (see Section 5), this is the first paper to consider the problem of detecting symmetries in probabilistic models.
2 Background
We will represent probabilistic models using factor graphs defined by the pair of variables and factors . The factor graph includes all observations, e.g. all data points in a clustering problem. Some of the variables in represent parameters shared across data points (like the cluster parameters), while others may be perpoint hidden variables (like cluster assignments), and others are the observations themselves (which are fixed but still considered part of ). We will use to denote the space of variable values, so . The (unnormalized) posterior distribution over the variables given the data can be expressed as the product of the factors
(1) 
A factor need not depend on all parameters, but we write the product as above for simplicity. Inference tasks typically require us to manipulate the posterior distribution. For instance, we may wish to compute some expectation with respect to the posterior distribution or to compute the posterior marginal distribution of some variable.
A symmetry is a measurable function with a measurable inverse satisfying
(2) 
where the product is taken over nonprior terms and the two sides of the equation are viewed as functions of . We also require that the symmetries fix the components of corresponding to observed variables.
From now on, when writing the factors , it will be taken for granted that the factors we consider are not priors. As an aside, the division of factors into likelihood factors and prior factors makes sense in directed models. For arbitrary factor graphs, we can can try to divide the factors into those that get repeated with more data and those that do not repeat. We are interested in symmetries in the factors that repeat, since those are the ones that prevent convergence to a peak.
A statistical model can often be parameterized in different ways and can be expressed with or without hidden variables. For example, a mixture model can be written with explicit variables that select which cluster a data point came from, or it can be written as a summation over clusters, i.e. with these selectors integrated out. These different options all lead to different factor graphs, and potentially to different symmetries. Thus it doesn’t make sense to say that a model (in the abstract sense) has symmetries, only that a particular factor graph encoding of that model has symmetries.
3 The Problem with Parameter Symmetries
In Bayesian approaches to machine learning, it is common to make predictions by manipulating the posterior distribution over the model parameters as expressed in Equation 1. In variational approaches to inference, the exact posterior is approximated by a simpler function , which is often factorized and unimodal. This approximation is effective because the posterior tends to be increasingly peaked as we receive more and more data. For a precise statement of this result, see van der Vaart (2000, Chapter 10.2).
Unfortunately, in many real applications, the likelihood function (which is given by the product of a subset of the factors in the factor graph) has symmetries in . For instance, the likelihood may be invariant to some permutation of the components of (as in mixture models), or the likelihood may be invariant to adding a constant to certain parameters (as in ranking models). These symmetries mean that is no longer identifiable from the data, and therefore the posterior does not become a simple peak as we get more data. Instead, the symmetries of the likelihood are preserved in the posterior regardless of how much data we receive.
As a result, variational approximations can be inaccurate, even with large quantities of data. For example, suppose the likelihood in some model depends only on . We may observe a lot of data informing us that . If the prior is for , then the marginal distribution for either parameter is . If we use a factorized approach that takes the product of these marginals as the approximate posterior and then try to make a prediction for , we arrive at the very uncertain answer of . In this example, the problem in the parameterization is easy to see, but in a more complicated model, the symmetries can be difficult to find.
Our motivation stems from problems like this one, which is specific to variational inference. However, there are other reasons to be concerned about symmetries.
In particular, parameter symmetries impact the reporting of results. A common practice is to report the mean and variance of each parameter in isolation. But in the presence of symmetries, these can be inaccurate or uninformative (in a mixture model, all of the mixture components may have the same marginals). Visualizing or analyzing the distance between latent features can be misleading if the different features are invariant to scaling. For instance, in a collaborative filtering model based on matrix factorization, any component of the user features can be scaled as long as the corresponding component of the item features is scaled in the opposite direction.
More subtle problems can arise from the interpretation of nonidentifiable parameters. As shown in Tsiatis (1975)
, the joint distribution of the
potential survival times in the competing risks model is nonidentifiable from observations of and . Overlooking this property may lead people to mistakenly believe that the ’s are independent when they are not. Symmetries can also lead to confusion when attempting to compute model evidence, as discussed in Neal (1999).Symmetries do not affect the correctness of sampling algorithms, though the presence of strong correlation between variables and of multiple modes can be difficult for many sampling algorithms to handle.
4 Uses for Symmetry Detection
As evidenced by the many papers discussing the symmetries present in specific models (see Section 5), symmetry detection has a variety of uses. Though this paper does not seek to innovate on this front, we nevertheless briefly review several possible uses of automatic symmetry detection. These uses include providing warnings to users of probabilistic programming languages, sanity checking models, and improving inference by automatically breaking model symmetries.
We do not expect all users to be familiar with the problem of model symmetries. By automatically detecting model symmetries, we may help users avoid misinterpreting the results of inference. The warnings that we display can be tailored to the class of model symmetries present as well as to the particular inference algorithm being used. For example, someone running expectation propagation (Minka, 2001)
on a model with permutation symmetries may like to know that the resulting approximate posterior may cover multiple symmetric modes and so its mean and variance may be uninformative. Meanfield approximations may tend to seek individual modes, underestimating posterior variance, and therefore producing reasonable point estimates (even in the presence of symmetries). In such situations, users may try to interpret latent parameters that are in fact nonidentifiable from the data. If the model in question possesses translation or scaling symmetries, we may like to alert the user to the existence of a continuum of comparable point estimates.
When developing a complex model, reporting symmetries may be a useful sanity check. For example, if a parameter meant to denote “age” participates in a signflip symmetry, something has probably gone wrong in the model specification. More generally, if any variable intended to represent some property of the data is in fact nonidentifiable from the data, there may be a problem.
Automatic symmetry breaking is the most obvious use of automatic symmetry detection. It is an interesting and important topic, but it is not the focus of this work. Many approaches have been discussed for breaking symmetry in specific models with known symmetries. Many of the techniques used to deal with specific symmetries (such as constraining or reparameterizing the model) generalize to broad classes of model symmetries. If used in conjunction, we could automatically detect symmetries, then automatically break them, then run inference. Inference algorithms may be better able to summarize the posterior after the symmetries have been broken. Even sampling algorithms may perform better as symmetry breaking can reduce the multimodality and correlation in the posterior. See Steyvers et al. (2009) for an example of symmetry breaking within Gibbs sampling. We can also use symmetry breaking as a postprocessing step. Someone doing inference with sampling may wish to project samples from the symmetric posterior onto a constrained subspace so as to give more meaningful marginal distributions to the parameters.
There are a number of suggestions from the literature on how to break symmetry. Several examples are discussed in Gelman and Hill (2007, Chapters 5.3, 5.8, 14.3). Other examples appear in Bafumi et al. (2005); Nobile (1998); Stephens (2000); Buesing et al. (2012); Palomo et al. (2007); Erosheva and Curtis (2011); Lopes and West (2004). But before we can build an automatic symmetry breaker, we need to understand which symmetries we can detect.
5 Related Work
A variety of probabilistic programming languages currently exist or are under development. These languages include BLOG (Milch et al., 2007), BUGS (Lunn et al., 2000), Church (Goodman et al., 2008), Infer.NET (Minka et al., 2012), and Stan (Stan Development Team, 2013). Our implementations of the algorithms in this paper operate on models specified in Infer.NET, but the algorithms are general and can be adapted to other model specification languages.
Much effort has been invested in detecting symmetries in constraintsatisfaction problems (Puget, 2005) and mathematical programs (Liberti, 2012). These methods are most closely related to the permutation symmetries that we discuss in Section 8 and are typically handled by reducing symmetry detection to a graph automorphism problem. These algorithms can easily be adapted to find permutation symmetries in a probabilistic model.
The lifted inference community is also interested in finding permutation symmetries in models. In this context, symmetries arise from the interplay between a particular model and a particular inference algorithm, and they are viewed as an asset which can be exploited to speed up inference algorithms such as belief propagation (Kersting et al., 2009). In our context, the symmetries that we consider exist independent of how inference is performed. Furthermore, we wish to find and remove symmetries in order to improve the performance of inference algorithms and to improve the interpretability of the results.
One way to avoid parameter symmetries is to directly approximate the predictive distribution, instead of the posterior, as described in (Snelson and Ghahramani, 2005). Unfortunately, the procedure suggested in their paper is quite expensive. Until a more efficient approach is found, approximating the posterior is still the most practical option.
Symmetries have been identified in specific models arising from problems in a variety of fields. The Rasch model is used in political science for ideal point estimation (Rasch, 1960), but as it is commonly implemented, it possesses scaling, signflip, and translation symmetries (Bafumi et al., 2005). In econometrics, the multinomial probit model is used to model choice behavior (Geweke et al., 1994). It too has scaling and translation symmetries (Nobile, 1998). Hierarchical models and mixture models are widely used in the social sciences (Draper, 1995) as well as in the biological sciences (Pritchard et al., 2000). Such models are rife with symmetries (Gelman and Hill, 2007; Stephens, 2000). Samaniego (2010) considers a material reliability problem in which the stress and strength parameters exhibit scaling symmetry.
Neal (1998)
pointed out that softmax classifiers have a parameter symmetry. He did inference by sampling without breaking the symmetry. Later work adopted the same model in the context of variational inference without breaking the symmetry
(Girolami and Rogers, 2006; Girolami and Zhong, 2007; Kim and Ghahramani, 2006). Buesing et al. (2012) considered parameter symmetries in linear dynamical systems and proposed a constraint to break them. Triggs et al. (2000) borrowed terminology from physics and described a parameter symmetry as a gauge freedom, and symmetry breaking as gauge fixing. However, they considered optimization problems in which gauge fixing affected only the computational complexity and not the final answer.6 Local Parameter Symmetries
Though we are interested in all symmetries that preserve the likelihood of the data, we do not imagine that all conceivable symmetries can be found efficiently. Therefore, we are motivated to consider classes of symmetries for which efficient detection algorithms can be constructed. A local symmetry is a symmetry that preserves the likelihood at each nonprior factor. In contrast with Equation 2, a local symmetry satisfies
(3) 
for all nonprior factors . This formulation has two advantages. First, we believe that it encompasses many of the symmetries that occur in practice. Second, to compute the local symmetries of a model, we can consider each factor on its own and ignore higherorder interactions.
6.1 Nonlocal Symmetries
Not all symmetries are local. Permutation symmetries, which we discuss in Section 8, are quite common and do not fall in this category. Less common and more complex symmetries can arise from mathematical identities. To illustrate this phenomenon, consider the following model. We draw from a prior over and from a prior over . We then draw
(4)  
(5) 
independently. Then we observe the data
(6)  
(7) 
This model possesses the symmetry given by
(8)  
(9)  
(10) 
This symmetry preserves the product of the nonprior factors. In particular, it preserves the equalities given in Equation 6 and Equation 7, and it preserves the product
(11) 
but it does not preserve the individual Gaussian factors.
We believe that symmetries like this one are less common than the local symmetries that we focus on, and we suspect that their presence is more difficult to detect. For these reasons, we do not give an algorithm for detecting more general types of symmetries.
7 Automatically Detecting Local Symmetries
Suppose that we have identified some specific class of transformations of the model parameters (the elements of need not be symmetries) and that we would like to compute the subset of local symmetries contained in defined by
(12) 
We can compute as
(13) 
In order to compute the class of symmetries , we require each factor to be annotated with the constraints that it imposes on transformations . It is important to note that these annotations need only be made once for each factor by the creators of the model specification language. There are currently about factors in Infer.NET (Minka et al., 2012), most of which do not require any annotation (as long as we interpret the lack of any annotation on a factor to mean that the arguments of that factor may not be transformed by any symmetry in ). In fact, the effort required to add these annotations is minimal. We will then be able to compute by aggregating the constraints from each factor and solving the resulting system of equations. The ease with which this formulation gives rise to an algorithm for computing depends on the nature of , but in several useful cases,
will be a vector space and the constraints from Equation
3 will be linear. In such situations, will also be a vector space. We now demonstrate this procedure for several specific classes of symmetries.7.1 Scaling Symmetries
A scaling symmetry is a symmetry which multiplies pointwise by a vector
(14) 
where each and the exponent (the case of is covered by a combination of scaling symmetry and signflip symmetry). We will represent this symmetry by the vector of exponents . In this case, the underlying space of multiplicative transformations is represented by the vector space of exponents.
In order to detect scaling symmetries, we will require each factor to be annotated with the constraints that it imposes on the scaling of its arguments. To illustrate this concept, consider the factor encoding the constraint . Any scaling symmetry which hopes to preserve this constraint must scale , , and by the same amount. In other words, we must have . A list of example factors and the constraints they impose is shown in Figure 1.
factor  constraints 

no constraints 
The constraints imposed on by some factor are typically linear. Indeed, any integer linear combination of scaling symmetries is again a scaling symmetry. It is possible that a particular factor will impose a nonlinear constraint on (for instance, a factor enforcing the constraint that be a power of ), but such factors are rare and are typically not implemented in probabilistic programming languages. Though this framework can accommodate such factors, we will require factors to impose linear constraints on scaling symmetries so that the resulting set of scaling symmetries is a vector space.
Now, if we construct the matrix by stacking together all of the constraints (as row vectors), the scaling symmetries of the model are the vectors in the null space of . We can compute the scaling symmetries of the model by constructing the matrix and finding its null space using an algorithm such as Gaussian elimination.
To illustrate the algorithm, we work through an example. Consider the simple model in Figure 2. Our algorithm ignores the Gaussian priors. It then collects the constraints imposed by each factor on potential scaling symmetries. The “plus” factor imposes the constraint . The “multiply” factor imposes the constraint . The observation adds the requirement that .
The scaling symmetries of the model are then the solutions to the equation
(15) 
The rows of the matrix are labeled with the constraints from which they are derived. The solution to Equation 15 is a onedimensional vector space given by , which translates into the space of scaling symmetries given by
(16)  
(17)  
(18)  
(19)  
(20) 
Remark 1
We mentioned earlier that symmetries are not allowed to change the values of observed variables. However, one of the most common observed values is , which scaling does not affect. It is often convenient to allow the symmetries we search for to scale observed values that equal (doing so allows us to to find the scaling symmetries in the factor and in the factor for instance). The same is true for signflip symmetries.
7.2 SignFlip Symmetries
A signflip symmetry is a symmetry in which is multiplied pointwise by a vector
(21) 
with . As with scaling symmetries, we will represent this symmetry by the binary vector of exponents .
Our algorithm for computing signflip symmetries is the same as our algorithm for computing scaling symmetries except that some factors impose different constraints and the computations must be done modulo . As before, each factor must be annotated with the constraints that it imposes on signflip transformations. Examples are shown in Figure 3.
factor  constraints 

and  
Once again, the signflip symmetries of the model form a vector space (over the finite field with elements), the constraints imposed by the factors are linear, and the signflip symmetries of the model correspond to the null space of the matrix formed by aggregating the constraints from each factor. As before, we can compute the null space of using an algorithm such as Gaussian elimination.
7.3 Translation Symmetries
We could have defined a translation symmetry as a symmetry which translates by some vector . Then it would follow, as in the case of scaling symmetries and signflip symmetries, that each factor must impose linear constraints on potential translation symmetries and that the translation symmetries of the model are the vectors in the null space of the matrix formed by aggregating all of the individual constraints.
Unfortunately, this definition does not capture the full range of symmetries that we feel ought to fall in this category. To illustrate a symmetry that this definition misses, consider the Bayes point machine (Girolami and Rogers, 2006). This model maintains latent vectors for each of classes, and it labels a data point with label by computing a noisy score
(22) 
for each class (where represents the noise), and setting
(23) 
Given any vector of the appropriate dimension, the translation
(24) 
does not change the likelihood of the data. So far, this translation appears to satisfy our preliminary definition of a translation symmetry. However, the vectors are not the only parameters in the model. For the purpose of running inference, the model that we construct will also maintain latent variables for the scores . When we translate each by , we must translate each by . If we wish to categorize such symmetries as translation symmetries, we must allow the vector by which we translate to depend on . Because of this dependence, translation symmetries differ from scaling and signflip symmetries.
With this example in mind, we define a translation symmetry to be a symmetry satisfying
(25) 
where such that can only depend on the parameters not being translated (i.e. on the such that ). To motivate this requirement, we note that it leads to the result that the existence of a translation symmetry as above implies the existence of a family of translation symmetries
(26) 
satisfying properties such as
(27) 
at least for integer . The existence of this family is consistent with an intuitive notion of the properties that a translation symmetry ought to satisfy.
As before, each factor must be annotated with the constraints that it imposes on potential translation symmetries. Examples are shown in Figure 4.
factor  constraints 

and  
and  
Therefore, we can compute the translation symmetries of the model by aggregating all of the constraints and solving the resulting (nonlinear) system of equations. The translation symmetries of a model can be written as a union of vector spaces, and we can compute this set by solving several separate linear systems. Since we allow the constraints to depend on , the coefficients in the resulting system of equations may be symbolic.
To illustrate the way in which translation symmetries may arise, consider the simple model shown in Figure 5. In the model, , , and are drawn from Gaussian priors, and and are observed. Our algorithm ignores the Gaussian priors and finds the translation symmetries
(28)  
(29)  
(30)  
(31)  
(32)  
(33)  
(34) 
for all , and
(35)  
(36)  
(37)  
(38)  
(39)  
(40)  
(41) 
for all .
8 Permutation Symmetries
A permutation symmetry is a symmetry that permutes the components of . Such symmetries arise in models with multiple interchangeable units such as the mixture components in a mixture of Gaussians, the topics in Latent Dirichlet Allocation (Blei et al., 2003)
, the hidden units in a neural network, and the latent features in a collaborative filtering model based on matrix factorization.
In a mixture model with symmetric priors, all mixture components have the same marginal distributions. This property can make it difficult to summarize the output of a sampler. This is known as the labelswitching problem. Papers on Latent Dirichlet Allocation sometimes avoid this by using only one sample (Griffiths and Steyvers, 2004). Whether or not permutation symmetry poses a problem for variational inference can depend on how well separated the mixture components are. If the mixture components are wellseparated, variational inference tends to get stuck in one mode and so isn’t bothered by the symmetry. The posterior distribution is inaccurate, but not in a way that significantly impacts test results.
We give an algorithm for detecting permutation symmetries by reduction to a graph automorphism computation, a problem for which there is no known polynomial time algorithm. This approach is very similar to existing work in the constraint satisfaction and mathematical programming literature (Puget, 2005; Liberti, 2012), and we do not consider this section to be novel. In Section 9.1.2, we will give an efficient algorithm for finding many of the permutation symmetries that arise in practice by making use of arrays.
8.1 Detecting Permutation Symmetries
More precisely, we will label the edges and vertices of the factor graph so that the maps from the factor graph to itself that preserve the edge and vertex labels correspond to the permutation symmetries of the model that can be deduced without any knowledge of the mathematical properties of the factors other than symmetry of each factor with respect to its arguments (for instance, the factor encoding the constraint is symmetric with respect to and , but not with respect to and ).
Each type of factor must be annotated with its identity (for instance, all binary plus factors have the same label, but a binary plus factor and a ternary plus factor have different labels). Each factor must also associate a label with each of its arguments such that two arguments have the same label if and only if the factor is symmetric with respect to those two arguments (we will refer to the information contained in a factor’s argument labels as the “symmetry structure” of that factor). Now, in the factor graph, label each factor with its identity, and for each edge connecting variable to a factor , label the edge with the label that factor associates with the argument .
The permutation symmetries that we can deduce from these annotations (factor labels and factor symmetry structures) are the permutation symmetries such that
(42) 
for all factor labels . These permutation symmetries correspond to the automorphisms of the factor graph that preserve the factor and edge labels.
Now, we can take our factor graph and compute its automorphism group (McKay, 1981; Darga et al., 2004). As is common with these computations, the group of permutation symmetries that we compute can be compactly represented by a set of generators (symmetries that can be composed to produce the entire set). Despite the fact that the problem of detecting the existence of nontrivial graph automorphisms is not known to be solvable in polynomial time (Köbler et al., 1993), we expect it to perform well in practice because graph automorphism can be tested in linear time for almost all graphs (Babai et al., 1980). However, as we describe in Section 9.1.2, we will ultimately compress the factor graphs using arrays (plate notation) so that the graphs whose automorphisms we attempt to compute are very small (possessing on the order of several dozen vertices).
Note that, with this approach, it is important to write models so that the symmetries are not hidden because this approach will not examine the interactions between factors. This approach will find the symmetry of , but it will not find the symmetry of , where the first refers to an addition factor with three arguments and the second refers to the sequential use of two addition factors. Therefore, when summing more than two variables, it will be important to use an ary version of the sum factor (we point out that it is not possible to find all of the symmetries of using only the information that we considered, namely the commutativity of addition, in this case, associativity is necessary as well, which we did not take into account). This problem does not arise with the other forms of symmetry that we have discussed in this paper.
9 Arrays, For Loops, and If Statements
Probabilistic programs often include support for conventional programming constructs such as arrays, for loops, and if statements, so it is important for our algorithms to be compatible with models specified using these objects.
9.1 Arrays
The detection methods, as we described them so far, operate on the full factor graph. A matrix factorization problem with millions of users, thousands of items, and hundreds of underlying features will have an unmanageably large factor graph. However, such large models do not have billions of different variables, each with their own unique interpretations. These models typically have billions of variables many of which share similar modeling roles (indeed, specifying the model in code requires on the order of one hundred lines as opposed to one billion). When defining a model in code, it is generally convenient to declare all of these similar variables at once using an array (this is analogous to drawing a graphical model using plate notation).
Intuitively, the symmetries of a model should not depend on the sizes of the arrays involved (though it is possible to contrive examples of symmetries that appear only for arrays of certain sizes). For instance, the translation symmetries in the Bayes point machine persist when we add or remove data and when we add or remove classes. For this reason, we ought to be able to formulate an algorithm for detecting the symmetries of our model with complexity independent of the sizes of the arrays involved. If we can adapt our algorithms to enable us to treat arrays of variables as individual variables, then we will have achieved this goal.
Note that the performance enhancement obtained from treating arrays as individual variables arises from the presumed repetitive structure that arrays represent. This assumption holds whenever arrays are used to describe the repetition in a graphical model arising from plate notation. However, a programmer could choose to create a model with many variables (each with its own distinct role) and store them all in an array. In such a situation, the presumed structure does not exist, and there is no computational benefit to be gained from the array. In order to be precise about the guarantees of our arrayboosted algorithms, we must either allow our algorithms to default to operating on the full factor graph or we must consider the possibility that we will not detect all symmetries. Both options have their advantages, and it may be the case that some combination of the two is ideal. However, due to the potential size of the full factor graph, we choose to devise an algorithm that never unrolls the factor graph and that instead bounds the set of symmetries between a subset and a superset.
We must adapt our definitions of the various classes of symmetries at the coarser resolution of arrays. We do this by treating arrays as individual variables. In the case of translation symmetries, the restriction we place on applies at the level of arrays (if a variable in an array is translated, then no variable in that array can affect the extent to which another variable is translated). In the case of permutation symmetries, we allow the permutation of entire arrays and the permutation of the indices within an array, but we do not allow permutations that would send some variables in the same array to different arrays.
9.1.1 Arrays with Scaling, SignFlip, and Translation Symmetries
Let be a representation of a scaling, signflip, or translation symmetry of the model. For concreteness, we will imagine that the vector consists of the exponents representing a scaling symmetry. If the variable is an array, then will be an array of exponents. Suppose that the arrays in the model are collectively indexed by the indices , where for each . Then we will use to indicate the element of indexed by . If is not indexed by some (or any) of the , then we simply ignore those indices. We detect symmetries in this context as follows.
The Simple Case: Suppose that for some model and some class of symmetries, the set of constraints imposed by the factors can be partitioned into subsets corresponding to the functions , where the set of constraints corresponding to is given by
(43) 
for all . We proceed by solving the system of equations given by
(44) 
Importantly, Equation 44 treats arrays of variables as individual variables and can be solved without unrolling the full factor graph. Now we can compactly represent the space of symmetries as
(45) 
Partial Solution to the General Case: Now suppose that we cannot partition the constraints into subsets of the form described in Equation 43. For instance, some constraint may single out some specific index in one of the ’s, or some constraint may relate multiple elements of one of the ’s. Let be the subvector of consisting of the arrays for which some constraint singles out an individual element or relates multiple elements. Requiring to be constant for each (in the sense that does not change as we vary ) allows us to write the constraints as in Equation 43 and places us back in the simple case. Solving for the space as before, we can find and represent all symmetries subject to the constraint that is constant for each .
This algorithm imposes additional constraints and therefore finds a subset of the symmetries. We can find a superset of the symmetries by following the same procedure but ignoring the constraints that relate multiple elements of the same array (or index specific elements). Though we have not implemented this suggestion, we speculate that we can then tighten the superset by attempting to prune the resulting solutions that are inconsistent with the ignored constraints (it may make sense to invest effort in doing this for problematic factors that occur frequently).
9.1.2 Arrays with Permutation Symmetries
Permutation symmetries interchange variables with similar roles, and variables with similar roles are often declared using arrays. Arrays in Infer.NET (Minka et al., 2012) are indexed by “range” objects (for instance, in a neural network example, one range may index the data set, one range may index the components of a data point, and one range may index the hidden units). In many common mixture model examples, permutation symmetries correspond to the permutation of the indices of a single range (the hidden units in this case).
A given range may index multiple arrays. For instance, the range corresponding to hidden units must index the array storing the weights into the hidden units, the array storing the weights out of the hidden units, and the array storing the biases of the hidden units (at least). In general, the indices of a range will be permutable if each factor is symmetric with respect to permutations of the indices of that range (this is often the case, although it would not be the case for some factors like a cumulative sum factor or a factor that accesses a specific element of the array) and if the range does not index an array containing observed values. If an array is indexed by a random variable, we check to see if any factor prevents the permutation of the values taken on by the random index (and if so, the permutation of the corresponding range is disallowed).
Now, we restrict our search for permutation symmetries to two separate searches, one for permutations of variables as before (but now treating arrays of variables as individual variables) and one for permutations of the indices of a range.
Note that when we treat an array of variables as a single variable and then proceed to search for permutation symmetries, we exclude the possibility of finding permutation symmetries that send elements of the same array to different arrays (something that makes little semantic modeling sense). In return, the complexity of our permutation detection algorithm becomes independent of the sizes of the arrays. We believe that this tradeoff is a reasonable one.
9.2 If Statements
Gate notation can be used to describe contextsensitive independence in factor graphs (Minka and Winn, 2008) and are one way of implementing if statements and switch statements in probabilistic programs. All of the symmetry detection algorithms that we have described are compatible with gates. Consider the example
bool b = Bernoulli(0.5) if (b) x = y else x = z
Local symmetries must be compatible with the constraints imposed by the factors in both the if block and the else block. In this example, in the case of scaling symmetries, the if block imposes the constraint , and the else block imposes the constraint . The ifelse block could equivalently be rewritten as a pair of factors, and by doing so, we would arrive at the same constraints.
This procedure will produce no false positives, but there are situations in which it will be overly restrictive. For instance, if the condition of some if block is always false, there is no need to impose the constraints from the factors in that block. If an if block and the corresponding else block both use a temporary variable with the same name, and if this variable is never accessed outside of the ifelse block, then the two instances of the variable ought to be treated as separate variables. Consider the following example.
bool b = Bernoulli(0.5) double x if (b) x = Gaussian(0, 1) y1 = x + 1 else x = Gaussian(0, 1) y2 = x ^ 2
The if block prevents any symmetry from scaling x, and the else block prevents any symmetry from translating x. As a consequence, no symmetries will be able to scale or translate y1 or y2, the variable of interest. But if x is a temporary variable that is never used outside of the ifelse block, then the two instances of x ought to be treated as different variables. To solve this problem, we could rewrite the code as below.
bool b = Bernoulli(0.5) double x1, x2 if (b) x1 = Gaussian(0, 1) y1 = x1 + 1 else x2 = Gaussian(0, 1) y2 = x2 ^ 2
Doing so would enable us to find symmetries that translate y1 and scale y2.
In a probabilistic programming language that allows users to create very general programs, these kinds of edge cases are bound to arise. Such edge cases will not lead our algorithms to produce any false positives, but they can lead to the omission of some symmetries. This is an argument for running symmetry detection algorithms after performing various program canonicalizations (such as pruning unreachable code and renaming variables). In fact, these program transformations are often already done within inference engines for other reasons.
10 Some Demonstrations
We have evaluated our algorithms on a number of models, including a neural network, Latent Dirichlet Allocation, and a mixture of Gaussians model. Our algorithms discovered all of the symmetries in these models. In this section, we give a detailed description of three more models and the symmetries that our algorithms detect in these models. These examples showcase the capacity of our algorithms to find symmetries in complicated models and also highlight some of the limitations of our algorithms.
10.1 Multinomial probit classifier
The multinomial probit classifier, which is described in Girolami and Rogers (2006), maintains latent vectors for each of classes. It assigns a data point to class by computing a noisy score via
(46)  
(47) 
for each class, and setting
(48) 
The factor graph for this model is shown in Figure 6.
This model has no permutation, scaling, or signflip symmetries. However, our algorithm finds the translation symmetries given by
(49)  
(50)  
(51) 
for all vectors of the same dimension as the data.
As pointed out in section 2, the way that a model is expressed as a factor graph can affect its symmetries. We illustrate this by writing the multinomial probit model differently, using an additional set of hidden variables representing the additive noise
(52)  
(53)  
(54) 
If we consider Equation 53 to be a prior factor, then we now have a scaling symmetry
(55)  
(56)  
(57)  
(58) 
and the larger set of translation symmetries
(59)  
(60)  
(61)  
(62) 
So which of these factor graphs is best? Since we are interested in symmetries that affect inference, it ultimately depends on what inference algorithm is being used. If the algorithm treats as a separate variable with its own approximate distribution, then we should include it in the factor graph and the extra symmetries above will be genuine symmetries of the algorithm. On the other hand, if the inference algorithm analytically marginalizes out (as most would), then should not appear in the factor graph and these extra symmetries will not occur.
10.2 Difficulty Versus Ability
The difficulty versus ability model is a generative model for how students perform on multiplechoice tests (Bachrach et al., 2012). Each question has a latent difficulty , a level of discrimination , and a correct answer . Each student has a latent ability . The advantage of student over question is . A noisy version of the advantage is computed as
(63) 
Student then answers question correctly if . Otherwise, the students response is random. The factor graph for this model is shown in Figure 7.
This model has no permutation or signflip symmetries. Our algorithm finds the scaling symmetries
(64)  
(65)  
(66)  
(67)  
(68) 
for any and the translation symmetries
(69)  
(70) 
for any .
10.3 Collaborative Filtering
The collaborative filtering model that we consider maintains latent features for each of users and latent features for each of items. There is also an dimensional vector of user biases and an dimensional vector of item biases . The affinity of user for item is given by
(71) 
Noisy versions of some of the affinities, given by
(72) 
are observed. The factor graph for this model is shown in Figure 8.
Our algorithm finds the permutation symmetries in which we permute any of the underlying features indexed by . Our algorithm finds the scaling symmetries given by
(73)  
(74) 
for all , such that for all . Similarly, we find the signflip symmetries
(75)  
(76) 
for all , such that for all .
The model’s translation symmetries are a bit more complex. Note that the summation factor introduces the constraint
(77) 
which introduces constraints relating the translation constants of elements of the same array. In the first case, we apply the most restrictive form of our algorithm in which the are constrained to all be the same and the are constrained to all be the same. In this situation, we find the translation symmetries
(78)  
(79) 
for all . These symmetries are a subset of the model’s translation symmetries. This misses the set of translation symmetries
(80)  
(81) 
and
(82)  
(83) 
By ignoring the constraint imposed by the summation factor and solving the remaining system of equations, we can find a superset of the translation symmetries that includes these omitted symmetries. However, there will be some spurious solutions, and there will be remaining work to do in pruning away the solutions that are inconsistent with the omitted constraint. For any particular factor, (such as the ary sum), we may be able to automate this pruning procedure.
Collaborative filtering is an interesting example because this model in fact contains a much broader form of rotational symmetry. For any invertible matrix , the transformation
(84)  
(85) 
is a symmetry of the model. The permutation, scaling, and signflip symmetries that we find are special cases of this matrixmultiplication symmetry. These matrix multiplication symmetries are local, and we could devise a class of transformations encompassing them (for instance, by allowing the components of the transformations to be linear combinations of the variables in the model). We would even be able to represent the constraints imposed by the factors as a disjunction of linear constraints. Doing this would require introducing a large number of auxiliary variables, and solving the system would potentially be expensive, but it could be done. However, the presence of such a class of symmetries will imply the presence of various permutation, scaling, and signflip symmetries as in the collaborative filtering example, all of which are much simpler to detect.
11 Discussion
Parameter symmetries are ubiquitous, and the presence of symmetries in a model can degrade the quality and interpretability of inference. Therefore, symmetries present an obstacle to the success of generalpurpose probabilistic programming. In the context in which many nonexpert machine learning practitioners seek to create novel models to describe their specific research problems, we will not be able to rely on every practitioner to be familiar with the problem of parameter symmetries or to know how to find them. As such, it will be important to have mechanisms in place for automatically detecting the presence of parameter symmetries and modifying the model so as to break the symmetries.
With this goal in mind, we introduced the concept of local parameter symmetries and described a procedure for constructing algorithms to automatically detect the presence of these symmetries. We illustrated this procedure by deriving algorithms for automatically detecting scaling, signflip, and translation symmetries. We discussed more general types of symmetries, the most common of which are permutation symmetries arising in mixture models. Our implementation of these algorithms works for models specified in Infer.NET (Minka et al., 2012), but the algorithms are general and can be easily adapted to any model specification language that arises in the future.
Graphical models can be massive, especially when we have large quantities of data or when the models are specified with plate notation. In these situations, we described how our algorithms can accommodate models specified with for loops and arrays and how the complexity of our symmetry detection algorithms is independent of the sizes of the arrays in the model. However, the fact that our algorithms do not depend on the sizes of the arrays in the model means that we compromise the exactness of our algorithms in some settings. As we seek to represent symmetries by maintaining only a single value for each array, our algorithms have trouble finding symmetries that require specifying the interactions between different elements of the same array. This is primarily a problem when it comes to representing translation symmetries in models with that use the ary sum factor. In such cases, we return a subset of the translation symmetries of the model. In total, however, our methods are powerful enough to detect most of the symmetries that we have found to arise in practice, as is demonstrated in our empirical study.
One area of future work will be to restore the exactness of our algorithms in the presence of arrays. Ideally, our algorithms would scale with the length of the model specification (regardless of whether or not the specification uses arrays). In such cases, when the arrays used in the models do not correspond to repetition in the model, our algorithms would first have to unroll the factor graph to the point where sufficient repetition exists for our algorithms to behave exactly.
The most important area of future work is combining the methods of this paper with automated techniques for breaking symmetries. In conjunction with specific techniques from the literature on breaking model symmetries, these results can enable the automatic detection and removal of parameter symmetries without any effort on the part of the modeler. Though much work has been done on breaking symmetry in specific models, there is work to be done in automating the procedure.
Acknowledgements
We would like to thank Chris Maddison and Boris Yangel for insightful discussions.
References
 Babai et al. (1980) L. Babai, P. Erdős, and S. M. Selkow. Random graph isomorphism. SIAM Journal on Computing, 9(3):628–635, 1980.
 Bachrach et al. (2012) Y. Bachrach, T. Minka, J. Guiver, and T. Graepel. How to grade a test without knowing the answers – a Bayesian graphical model for adaptive crowdsourcing and aptitude testing. In Proceedings of the 29nd international Conference on Machine Learning, pages 1183–1190, 2012.
 Bafumi et al. (2005) J. Bafumi, A. Gelman, D. K. Park, and N. Kaplan. Practical issues in implementing and understanding Bayesian ideal point estimation. Political Analysis, 13(2):171–187, 2005.
 Bishop (2013) C. M. Bishop. Modelbased machine learning. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984), 2013.
 Blei et al. (2003) D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent Dirichlet allocation. The Journal of Machine Learning Research, 3:993–1022, 2003.
 Buesing et al. (2012) L. Buesing, J. H. Macke, and M. Sahani. Learning stable, regularised latent models of neural population dynamics. Network: Computation in Neural Systems, 23(12):24–47, 2012.
 Darga et al. (2004) P. T. Darga, M. H. Liffiton, K. A. Sakallah, and I. L. Markov. Exploiting structure in symmetry detection for CNF. In Proceedings of the 41st annual Design Automation Conference, pages 530–534. ACM, 2004.
 DARPA (2013) DARPA. Broad agency announcement: Probabilistic programming for advanced machine learning, 2013.
 Draper (1995) D. Draper. Inference and hierarchical modeling in the social sciences. Journal of Educational and Behavioral Statistics, 20(2):115–147, 1995.
 Erosheva and Curtis (2011) E. A. Erosheva and S. M. Curtis. Dealing with rotational invariance in Bayesian confirmatory factor analysis. Technical report, University of Washington, 2011.
 Gelman and Hill (2007) A. Gelman and J. Hill. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, 2007.
 Geweke et al. (1994) J. Geweke, M. Keane, and D. Runkle. Alternative computational approaches to inference in the multinomial probit model. The review of economics and statistics, pages 609–632, 1994.
 Girolami and Rogers (2006) M. Girolami and S. Rogers. Variational Bayesian multinomial probit regression with Gaussian process priors. Neural Computation, 18(8):1790–1817, 2006.
 Girolami and Zhong (2007) M. Girolami and M. Zhong. Data integration for classification problems employing Gaussian process priors. In Advances in Neural Information Processing Systems: Proceedings of the 2006 Conference, volume 19, pages 465–472, 2007.

Goodman et al. (2008)
N. Goodman, V. Mansinghka, D. Roy, K. Bonawitz, and J. Tenenbaum.
Church: a language for generative models.
In
Proceedings of the 24th Conference on Uncertainty in Artificial Intelligence
, 2008.  Griffiths and Steyvers (2004) T. L. Griffiths and M. Steyvers. Finding scientific topics. Proceedings of the National academy of Sciences, 101(Suppl 1):5228–5235, 2004.
 Kersting et al. (2009) K. Kersting, B. Ahmadi, and S. Natarajan. Counting belief propagation. In Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence, pages 277–284, 2009.
 Kim and Ghahramani (2006) H.C. Kim and Z. Ghahramani. Bayesian Gaussian process classification with the EMEP algorithm. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 28(12):1948–1959, 2006.
 Köbler et al. (1993) J. Köbler, U. Schöning, and J. Torán. The Graph Isomorphism Problem: Its Structural Complexity. Birkhäuser, 1993.
 Liberti (2012) L. Liberti. Reformulations in mathematical programming: automatic symmetry detection and exploitation. Mathematical Programming, 131:273–304, 2012.
 Lopes and West (2004) H. F. Lopes and M. West. Bayesian model assessment in factor analysis. Statistica Sinica, 14:41–67, 2004.
 Lunn et al. (2000) D. J. Lunn, A. Thomas, N. Best, and D. Spiegelhalter. WinBUGSa Bayesian modelling framework: concepts, structure, and extensibility. Statistics and Computing, 10(4):325–337, 2000.
 McKay (1981) B. D. McKay. Practical graph isomorphism. Congressus Numerantium, 30:45–87, 1981.
 Milch et al. (2007) B. Milch, B. Marthi, S. Russell, D. Sontag, D. L. Ong, and A. Kolobov. BLOG: Probabilistic models with unknown objects. Statistical relational learning, page 373, 2007.
 Minka and Winn (2008) T. Minka and J. Winn. Gates. In Advances in Neural Information Processing Systems, pages 1073–1080, 2008.
 Minka et al. (2012) T. Minka, J. Winn, J. Guiver, and D. Knowles. Infer.NET 2.5, 2012. Microsoft Research Cambridge. http://research.microsoft.com/infernet.

Minka (2001)
T. P. Minka.
Expectation propagation for approximate Bayesian inference.
In Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence, pages 362–369, 2001.  Neal (1998) R. M. Neal. Regression and classification using Gaussian process priors. Bayesian Statistics, 6:475–501, 1998.
 Neal (1999) R. M. Neal. Erroneous results in“marginal likelihood from the Gibbs output”. 1999.

Nobile (1998)
A. Nobile.
A hybrid Markov chain for the Bayesian analysis of the multinomial probit model.
Statistics and Computing, 8(3):229–242, 1998.  Palomo et al. (2007) J. Palomo, D. B. Dunson, and K. Bollen. Bayesian structural equation modeling. Handbook of Latent Variable and Related Models, pages 163–188, 2007.
 Pritchard et al. (2000) J. K. Pritchard, M. Stephens, and P. Donnelly. Inference of population structure using multilocus genotype data. Genetics, 155(2):945–959, 2000.
 Puget (2005) J.F. Puget. Automatic detection of variable and value symmetries. In Principles and Practice of Constraint Programming  CP 2005, pages 475–489. Springer, 2005.
 Rasch (1960) G. Rasch. Studies in mathematical psychology: I. probabilistic models for some intelligence and attainment tests. 1960.
 Samaniego (2010) F. J. Samaniego. A comparison of the Bayesian and frequentist approaches to estimation. Springer, 2010.
 Snelson and Ghahramani (2005) E. Snelson and Z. Ghahramani. Compact approximations to Bayesian predictive distributions. In Proceedings of the 22nd international conference on Machine learning, pages 840–847. ACM, 2005.
 Stan Development Team (2013) Stan Development Team. Stan: A c++ library for probability and sampling, version 1.3, 2013. URL http://mcstan.org/.
 Stephens (2000) M. Stephens. Dealing with label switching in mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(4):795–809, 2000.
 Steyvers et al. (2009) M. Steyvers, B. Miller, P. Hemmer, and M. D. Lee. The wisdom of crowds in the recollection of order information. In Advances in Neural Information Processing Systems, pages 1785–1793, 2009.
 Triggs et al. (2000) B. Triggs, P. F. McLauchlan, R. I. Hartley, and A. W. Fitzgibbon. Bundle adjustment—a modern synthesis. In Vision Algorithms: Theory and Practice, pages 298–372. Springer, 2000.
 Tsiatis (1975) A. Tsiatis. A nonidentifiability aspect of the problem of competing risks. Proceedings of the National Academy of Sciences, 72(1):20–22, 1975.
 van der Vaart (2000) A. W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 2000.
Comments
There are no comments yet.