## Thurstonian Boltzmann Machines: Learning from Multiple Inequalities

Truyen Tran †‡ , Dinh Phung , Svetha Venkatesh † † † Pattern Recognition and Data Analytics, Deakin University, Australia ‡ Department of Computing, Curtin University, Australia { truyen.tran,dinh.phung,svetha.venkatesh } @deakin.edu.au

## Abstract

We introduce Thurstonian Boltzmann Machines (TBM), a unified architecture that can naturally incorporate a wide range of data inputs at the same time. Our motivation rests in the Thurstonian view that many discrete data types can be considered as being generated from a subset of underlying latent continuous variables, and in the observation that each realisation of a discrete type imposes certain inequalities on those variables. Thus learning and inference in TBM reduce to making sense of a set of inequalities. Our proposed TBM naturally supports the following types: Gaussian, intervals, censored, binary, categorical, muticategorical, ordinal, (in)-complete rank with and without ties. We demonstrate the versatility and capacity of the proposed model on three applications of very different natures; namely handwritten digit recognition, collaborative filtering and complex social survey analysis.

## 1 Introduction

Restricted Boltzmann machines (RBMs) have proved to be a versatile tool for a wide variety of machine learning tasks and as a building block for deep architectures [12, 24, 28]. The original proposals mainly handle binary visible and hidden units. Whilst binary hidden units are broadly applicable as feature detectors, non-binary visible data requires different designs. Recent extensions to other data types result in type-dependent models: the Gaussian for continuous inputs [12], Beta for bounded continuous inputs [16], Poisson for count data [9], multinomial for unordered categories [25], and ordinal models for ordered categories [37, 35].

The Boltzmann distribution permits several types to be jointly modelled, thus making the RBM a good tool for multimodal and complex social survey analysis. The work of [20, 29, 40] combines continuous (e.g., visual and audio) and discrete modalities (e.g., words). The work of [34] extends the idea further to incorporate ordinal and rank data. However, there are conceptual drawbacks: First, conditioned on the hidden layer, they are still separate type-specific models; second, handling ordered categories and ranks is not natural; and third, specifying direct correlation between these types remains difficult.

The main thesis of this paper is that many data types can be captured in one unified model. The key observations are that (i) type-specific properties can be modelled using

one or several underlying continuous variables , in the spirit of Thurstonian models 1 [31], and (ii) evidences be expressed in the form of one or several inequalities of these underlying variables. For example, a binary visible unit is turned on if the underlying variable is beyond a threshold; and a category is chosen if its utility is the largest among all those of competing categories. The use of underlying variables is desirable when we want to explicitly model the generative mechanism of the data. In psychology and economics, for example, it gives much better interpretation on why a particular choice is made given the perceived utilities [2]. Further, it is natural to model the correlation among type-specific inputs using a covariance structure on the underlying variables.

The inequality observation is interesting in its own right: Instead of learning from assigned values, we learn from the inequality expression of evidences , which can be much more relaxed than the value assignments. This class of evidences indeed covers a wide range of practical situations, many of which have not been studied in the context of Boltzmann machines, as we shall see throughout the paper.

To this end, we propose a novel class of models called Thurstonian Boltzmann Machine (TBM). The TBM utilises the Gaussian restricted Boltzmann machine (GRBM): The top layer consists of binary hidden units as in standard RBMs; the bottom layer contains a collection of Gaussian variable groups, one per input type. The main difference is that TBM does not require valued assignments for the bottom layer but a set of inequalities expressing the constraints imposed by the evidences . Except for a limiting case of point assignments where the inequalities are strictly equalities, the Gaussian layer is never fully observed. The TBM supports more data types in a unified manner than ever before: For any combination of the point assignments, intervals, censored values, binary, unordered categories, multi-categories, ordered categories, (in)-complete ranks with and without ties, all we need to do is to supply relevant subset of inequalities.

We evaluate the proposed model on three applications of very different natures: handwritten digit recognitions, collaborative filtering and complex survey analysis. For the first two applications, the performance is competitive against methods designed for those data types. On the last application, we believe we are among the first to propose a scalable and generic machinery for handle those complex data types.

## 2 Gaussian RBM

Let x = ( x , x 1 2 , ..., x N ) glyph[latticetop] ∈ R N be a vector of input variables. Let h = ( h , h 1 2 , ..., h K ) glyph[latticetop] ∈ { 0 1 , } K be a set of hidden factors which are designed to capture the variations in the observations. The input layer and the hidden layer form an undirected bipartite graph, i.e., only cross-layer connections are allowed. The model admits the Boltzmann distribution

$$P ( x, h ) = \frac { 1 } { Z } \exp \left \{ - E ( x, h ) \right \}$$

where Z = ∑ h GLYPH&lt;1&gt; exp {-E ( x h , ) } d x is the normalising constant and E ( x h , ) is the state energy. The energy is decomposed as

$$E ( x, h ) = \sum _ { i } \left ( \frac { x _ { i } ^ { 2 } } { 2 } - ( \alpha _ { i } + W _ { i } \bullet h ) \, x _ { i } \right ) - \gamma ^ { \top } h$$

1 Whilst Thurstonian models often refer to human's judgment of discrete choices, we use the term 'Thurstonian' more freely without the notion of human's decision.

where { α i } N i =1 , W = { W ik } , γ = { γ k } are free parameters and W i · denotes the i -th row. Given the input x , the posterior has a simple form

$$P ( h \, | \, x ) \ & = \ \prod _ { k } P ( h _ { k } \, | \, x ) \\ P ( h _ { k } = 1 \, | \, x ) \ & = \ \frac { 1 } { 1 + e ^ { - \gamma _ { k } - W ^ { \prime } _ { \bullet k } x } }$$

where W · k denotes the k -th column. Similarly, the generative process given the binary factor h is also factorisable

$$P ( x \, | \, h ) \ & = \ \prod _ { i } P ( x _ { i } \, | \, h ) \\ P ( x _ { i } \, | \, h ) \ & = \ \mathcal { N } ( \alpha _ { i } + W _ { i }, h, 1 )$$

where N ( µ, 1) is the normal distribution of mean µ and unit deviation.

## 3 Thurstonian Boltzmann Machines

We now generalise the Gaussian RBM into the Thurstonian Boltzmann Machine (TBM). Denote by e an observed evidence of x . Standard evidences are the point assignment of x to some specific real-valued vector, i.e., x = e . Generalised evidences can be expressed using inequality constraints

$$b \leq _ { \odot } A x \leq _ { \odot } c$$

for some transform matrix A ∈ R M N × and vectors b c , ∈ R M , where ≤ glyph[circledot] denotes elementwise inequalities. Thus an evidence can be completely realised by specifying the triple 〈 A, b c , 〉 . For example, for the point assignment, 〈 A = I , b = e c , = e 〉 , where I is the identity matrix. In what follows, we will detail other useful popular realisations of these quantities.

## 3.1 Boxed Constraints

This refers to the case where input variables are independently constrained, i.e., A = I , and thus we need only to specify the pair 〈 b c , 〉 .

Censored observations. This refers to situation where we only know the continuous observation beyond a certain point, i.e., b = e and c = + ∞ . For example, in survival analysis, the life expectancy of a person might be observed up to a certain age, and we have no further information afterward.

Interval observations . When the measurements are imprecise, it may be better to specify the range of possible observations with greater confidence rather than a singe point, i.e., b = e -δ and c = e + δ for some pair ( e δ , ). For instance, missile tracking may estimate the position of the target with certain precision.

Binary observations . A binary observation e i can be thought as a result of clipping x i by a threshold θ i , that is e i = 1 if x i ≥ θ i and e i = 0 otherwise. The boundaries in Eq. (3.2) become:

$$\langle b _ { i }, c _ { i } \rangle = \begin{cases} \langle - \infty, \theta _ { i } \rangle & \ e _ { i } = 0 \\ \langle \theta _ { i }, + \infty \rangle & \ e _ { i } = 1 \end{cases}$$

Thus, this model offers an alternative 2 to standard binary RBMs of [28, 8].

Ordinal observations. Denote by e = ( e , e 1 2 , ..., e N ) the set of ordinal observations, where each e i is drawn from an ordered set { 1 2 , , .., L } . The common assumption is that the ordinal level e i = l is observed given x i ∈ [ θ l -1 , θ l ] for some thresholds θ 1 ≤ θ 2 ≤ ...θ L -1 . The boundaries thus read

$$\langle b _ { i }, c _ { i } \rangle = \begin{cases} \langle - \infty, \theta _ { 1 } \rangle & l = 1 \\ \langle \theta _ { l - 1 }, \theta _ { l } \rangle & l = 2, 3,.. L - 1 \\ \langle \theta _ { L - 1 }, + \infty \rangle & \text{otherwise} \end{cases}$$

This offers an alternative 3 to the ordinal RBMs of [37].

## 3.2 Inequality Constraints

glyph[negationslash]

Categorical observations . This refers to the situation where out of an unordered set of categories, we observe only one category at a time. This can be formulated as follows. Each category is associated with a 'utility' variable. The category l is observed (i.e., e i = m ) if it has the largest utility, that is x il ≥ max m l = x im . Thus, x il is the upperthreshold for all other utilities. On the other hand, max m l = x im is the lower-threshold for x il . This suggests an EM-style procedure: ( ) fix i x il (or treat it as a threshold) and learn the model under the intervals x im ≤ x il for all m = , and ( l ii ) fix all categories other than l , learn the model under the interval x il ≥ max m l = x im . This offers an alternative 4 to the multinomial logit treatment in [26].

glyph[negationslash]

glyph[negationslash]

To illustrate the point, suppose there are only four variables z , z 1 2 , z 3 , z 4 , and z 1 is observed, then we have z 1 ≥ max { z , z 2 3 , z 4 } . This can be expressed as z 1 -z 2 ≥ 0; z 1 -z 3 ≥ 0 and z 1 -z 4 ≥ 0. These are equivalent to

$$\left \langle A = \left [ \begin{array} { c c c } 1 & - 1 & 0 & 0 \\ 1 & 0 & - 1 & 0 \\ 1 & 0 & 0 & - 1 \end{array} \right ] ; \ \ b = 0 ; \ \ c = + \infty \right \rangle$$

Imprecise categorical observations . This generalises the categorical case : The observation is a subset of a set, where any member of the subset can be a possible observation . 5 For example, when asked to choose the best sport team of interest, a person

2 To be consistent with the statistical literature, we can call it the probit RBM , which we will study in Section 6.1.

3 This can be called ordered probit RBM .

4 We can call this model multinomial probit RBM.

5 This is different from saying that all the members of the subset must be observed.

glyph[negationslash]

may pick two teams without saying which is the best. For instance, suppose the subset is { z , z 1 2 } , then min { z , z 1 2 } ≥ max { z , z 3 4 } , which can be expressed as z 1 -z 3 ≥ 0; z 2 -z 3 ≥ 0, z 1 -z 4 ≥ 0 and z 2 -z 4 ≥ 0. This translates to the following triple

$$\left \langle A = \left [ \begin{array} { c c c c } 1 & 0 & - 1 & 0 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \\ 0 & 1 & 0 & - 1 \end{array} \right ] ; \ \ b = 0 ; \ \ c = + \infty \right \rangle$$

Rank (with Ties) observations . This generalises the imprecise categorical cases : Here we have a (partially) ranked set of categories. Assume that the rank is produced in a stagewise manner as follows: The best category subset is selected out of all categories, the second best is selected out of all categories except for the best one, and so on. Thus, at each stage we have an imprecise categorical setting, but now the utilities of middle categories are constrained from both sides - the previous utilities as the upper-bound, and the next utilities as the lower-bound.

As an illustration, suppose there are four variables z , z 1 2 , z 3 , z 4 and a particular rank (with ties) imposes that min { z , z 1 2 } ≥ z 3 ≥ z 4 . This be rewritten as z 1 ≥ z 3 ; z 2 ≥ z 3 ; z 3 ≥ z 4 , which is equivalent to

$$\left \langle A = \left [ \begin{array} { c c c } 1 & 0 & - 1 & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & 0 & 1 & - 1 \end{array} \right ] ; \ \ b = 0 ; \ \ c = + \infty \right \rangle$$

## 4 Inference Under Linear Constraints

Under the TBM, MCMC-based inference without evidences is simple: we alternate between P ( h | x ) and P ( x | h ). This is efficient because of the factorisations in Eqs. (3,4). Inference with inequality-based evidence e is, however, much more involved except for the limiting case of point assignments.

Denote by Ω ( e ) = { x | b ≤ glyph[circledot] A x ≤ glyph[circledot] c } the constrained domain of x defined by the evidence e . Now we need to specify and sample from the constrained distribution P ( x h , | e ) defined on Ω ( e ). Sampling P ( h | x ) remains unchanged, and in what follows we focus on sampling from P ( x | h e , ).

## 4.1 Inference under Boxed Constraints

For boxed constraints (Section 3.1), due to the conditional independence, we still enjoy the factorisation P ( x | h e , ) = ∏ i P x ( i | h e , ). We further have

$$P ( x _ { i } \, | \, h, e ) \ = \ \frac { P ( x _ { i } \, | \, h ) } { \Phi ( c _ { i } \, | \, h ) - \Phi ( b _ { i } \, | \, h ) }$$

where Φ( · | h ) is the normal cumulative distribution function of P x ( i | h ). Now P x ( i | h e , ) is a truncated normal distribution, from which we can sample using the simple rejection method, or more advanced methods such as those in [23].

## 4.2 Inference under Inequality Constraints

For general inequality constraints (Section 3.2), the input variables are interdependent due to the linear transform A . However, we can specify the conditional distribution P x ( i | x ¬ i , h e , ) (here x ¬ i = x \ x i ) by realising that glyph[negationslash]

$$b _ { m } - \sum _ { j \neq i } A _ { m j } x _ { j } \leq A _ { m i } x _ { i } \leq c _ { m } - \sum _ { j \neq i } A _ { m j } x _ { j }$$

glyph[negationslash]

glyph[negationslash]

where A mi = 0 for m = 1 2 , , ..., M . In other words, x i is conditionally box-constrained given other variables.

glyph[negationslash]

This suggests a Gibbs procedure by looping through x , x 1 2 , ..., x N . With some abuse of notation, let ˜ b mi = ( b m -∑ j = i A mj x j ) /A mi and ˜ c mi = ( c m -∑ j = i A mj x j ) /A mi . The constraints can be summarised as glyph[negationslash]

$$x _ { i } \ & \in \ \cap _ { m = 1 } ^ { M } \left [ \min \left \{ \tilde { b } _ { m i }, \tilde { c } _ { m i } \right \}, \max \left \{ \tilde { b } _ { m i }, \tilde { c } _ { m i } \right \} \right ] \\ & = \ \left [ \max _ { m } \min \left \{ \tilde { b } _ { m i }, \tilde { c } _ { m i } \right \}, \min _ { m } \max \left \{ \tilde { b } _ { m i }, \tilde { c } _ { m i } \right \} \right ]$$

The min and max operators are needed to handle change in inequality direction with the sign of A mi , and the join operator is due to multiple constraints.

For more sophisticated Gibbs procedures, we refer to the work in [10].

## 4.3 Estimating the Binary Posteriors

We are often interested in the posteriors P ( h | e ), e.g., for further processing. Unlike the standard RBMs, the binary latent variables here are coupled through the unknown Gaussians and thus there are no exact solutions unless the evidences are all point assignments. The MCMC-based techniques described above offer an approximate estimation by averaging the samples { h ( s ) } S s =1 . For the case of boxed constraints, mean-field offers an alternative approach which may be numerically faster. In particular, the mean-field updates are recursive:

$$Q _ { i }$$

$$\begin{matrix} \\ \end{matrix}$$

$$Q _ { k } \quad & \leftarrow \quad \frac { 1 } { 1 + \exp \left \{ - \gamma _ { k } - \sum _ { i } W _ { i k } \hat { \mu } _ { i } \right \} } \\ \mu _ { i } \quad & \leftarrow \quad \alpha _ { i } + \sum _ { k } W _ { i k } Q _ { k } \\ \hat { \mu } _ { i } \quad & = \quad \mu _ { i } + \frac { \phi ( b _ { i } - \mu _ { i } ) - \phi ( c _ { i } - \mu _ { i } ) } { \Phi ( c _ { i } - \mu _ { i } ) - \Phi ( b _ { i } - \mu _ { i } ) } \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \end{array}$$

where Q k is the probability of the unit k being activated, ˆ µ i is the mean of the normal distribution truncated in the interval [ b , c i i ], φ z ( ) is the probability density function, and Φ( ) is normal cumulative distribution function. z Interested readers are referred to the Supplement 6 for more details.

6 Supplement material will be available at: http://truyen.vietlabs.com

## 4.4 Estimating Probability of Evidence Generation

Given the hidden states h we want to estimate the probability that hidden states generate a particular evidence e

$$P ( e \, | \, h ) = \int _ { \Omega ( e ) } P ( x \, | \, h ) d x$$

For boxed constraints, analytic solution is available since the Gaussian variables are decoupled, i.e., P e ( i | h ) = Φ( c i -µ i ) -Φ( b i -µ i ), where µ i = α i + ∑ k W h ik k . For general inequality constraints, however, these variables are coupled by the inequalities. The general strategy is to sample from P ( x | h ) and compute the portion of samples falling into the constrained domain Ω ( e ). For certain classes of inequalities we can approximate the Gaussian by appropriate distributions from which the integration has the closed form. In particular, those inequalities imposed by the categorical and rank evidences can be dealt with by using the extreme value distributions . The integration will give the logit form on distribution of categories and Plackett-Luce distribution of ranks. For details, we refer to the Supplement.

## 5 Stochastic Gradient Learning with Persistent Markov Chains

Learning is based on maximising the evidence likelihood

$$\mathcal { L } \ = \ \log P ( e ) = \log \sum _ { h } \int _ { \Omega ( e ) } P ( h, x ) d x$$

where P ( h x , ) is defined in Eq. (1). Let Z ( e ) = ∑ h GLYPH&lt;1&gt; Ω ( e ) exp {-E ( x h , ) } d x , then L = log Z ( e ) -log Z . The gradient w.r.t. the mapping parameter reads

$$\partial _ { W _ { i k } } \mathcal { L } \ = \ \mathbb { E } _ { P ( x _ { i }, h _ { k } | e ) } \left [ x _ { i } h _ { k } \right ] - \mathbb { E } _ { P ( x _ { i }, h _ { k } ) } \left [ x _ { i } h _ { k } \right ]$$

The derivation is left to the Supplement.

## 5.1 Estimating Data Statistics

The data-dependent statistics E P x ,h ( i k | e ) [ x h i k ] and the data-independent statistics E P x ,h ( i k ) [ x h i k ] are not tractable to compute in general, and thus approximations are needed.

Data-dependent statistics. Under the box constraints, the mean-field technique (Section 4.3) can be employed as follows

$$\mathbb { E } _ { P ( x _ { i }, h _ { k } | e ) } \left [ x _ { i } h _ { k } \right ] \approx \hat { \mu } _ { i } Q _ { k }$$

For general cases, sampling methods are applicable. In particular, we maintain one persistent Markov chain [42, 32] per data instance and estimate the statistics after a very short run. This would explore the space of the data-dependent distribution P ( x h , | e ) by alternating between P ( h | x ) and P ( x | h e , ) using techniques described in Section 4.

Data-independent statistics. Mean-field distributions are not appropriate for exploring the entire state space because they tend to fit into one mode. One practical solution is based on the idea of Hinton's Contrastive Divergence (CD), where we create another Markov chain on-the-fly starting from the latest state of the clamped chain. This chain will be discarded after each parameter update. This is particular useful when the models are instance-specific, e.g., in collaborative filtering, it is much cheaper to build one model per user, all share the same parameters. If it is not the case, then we can maintain a moderate set of parallel chains and collect the samples after a short run at every updating step [42, 32].

## 5.2 Learning the Box Boundaries

In the case of boxed constraints, sometimes it is helpful to learn the boundaries 〈 b , c i i 〉 themselves. The gradient of the log-likelihood w.r.t. the lower boundaries reads

$$\text{s. The gradient of the log-likelihood w.r.t. the lower boundaries read} \\ \partial _ { b _ { i } } \mathcal { L } \ & = \ \frac { 1 } { Z ( e ) } \sum _ { h } \partial _ { b _ { i } } \int _ { \Omega ( e ) } \exp \left \{ - E ( x, h ) \right \} d x \\ & = \ \sum _ { h } \partial _ { b _ { i } } \int _ { \Omega ( e ) } P \left ( x, h \, | \, e \right ) d x \\ & \approx \ \frac { 1 } { S } \partial _ { b _ { i } } \int _ { b _ { i } } ^ { c _ { i } } \int _ { \Omega ( e ) \left \lfloor b _ { i }, c _ { i } \right \rfloor } P \left ( x \, | \, h ^ { ( s ) }, e \right ) d x _ { - i } d x _ { i } \\ & = \ \frac { 1 } { S } \sum _ { h } P \left ( x _ { i } = b _ { i } \, | \, h ^ { ( s ) }, e \right ) \\ ( s ) \right \} S \ &, \quad \,, \quad \,, \quad \,, \quad \,, \quad \,, \quad \,.$$

where { h ( s ) } S s =1 are samples collected during the MCMC procedure running on the datadependent distribution P ( x h , | e ). Similarly we would have the gradient w.r.t. the upper boundaries:

$$\partial _ { c _ { i } } \mathcal { L } = \frac { 1 } { S } \sum _ { h } P \left ( x = c _ { i } \, | \, h ^ { ( s ) }, e \right ).$$

## 6 Applications

In this section, we describe applications of the TBM for three realistic domains, namely handwritten digit recognition, collaborative filtering and worldwide survey analysis. Before going to the details, let us first address key implementation issues (see Supplement for more details).

One observed difficulty in training the TBM is that the hidden samples can get stuck in one of the two ends and thus learning cannot progress. The reasons might be the large mapping parameters or the unbounded nature of the underlying Gaussian variables, which can saturate the hidden units. We can control the norm of the mapping parameters, either by using the standard glyph[lscript] 2 -norm regularisation, or by rescaling the norm of the parameter vector for each hidden unit. To deal with the non-boundedness of the Gaussian variables, then we can restrict their range, making the model bounded.

Another effective solution is to impose a constraint on the posteriors by adding the

regularising term to the log-likelihood, e.g.,

$$\lambda \left \{ \sum _ { k } \left [ \rho \log P ( h _ { k } ^ { 1 } \, | \, e ) + ( 1 - \rho ) \log \left ( 1 - P ( h _ { k } ^ { 1 } \, | \, e ) \right ) \right ] \right \}$$

where ρ ∈ (0 1) is , the expected probability that a hidden unit will turn on given the evidence and λ &gt; 0 is the regularisation weight. Maximising this quantity is essentially minimising the Kullback-Leibler divergence between the expected posteriors and the true posteriors. In our experiments, we found ρ ∈ (0 1 . , 0 3) and . λ ∈ (0 1 . , 1) gave satisfying results.

The main technical issue is that P h ( 1 k | e ) does not have a simple form due to the integration over all the constrained Gaussian variables. Approximation is thus needed. The use of mean-field methods will lead to the simple sigmoid form, but it is only applicable for boxed constraints since it breaks down deterministic constraints among variables (Section 4.3). However, we can estimate the 'mean' truncated Gaussian ˆ µ i by averaging the recent samples of the Gaussian variables in the data-dependent phase.

Once these safeguards are in place, learning can greatly benefit from quite large learning rate and small batches as it appears to quickly get the samples out off the local energy traps by significantly distorting the energy landscape. Depending on the problem sizes, we vary the batch sizes in the range [100 , 1000].

## 6.1 Probit RBM for Handwritten Digits

We use the name Probit RBM to denote the special case of TBM where the observations are binary (i.e., boxed constraints, see Section 3.1). The threshold θ i for each visible unit i is chosen so that under the zero mean, the probability of generating a binary evidence equals the empirical probability, i.e., 1 -Φ( θ i ) = ¯ , and thus e i θ i = Φ -1 (1 -¯ ). e i Since any mismatch in thresholds can be corrected by shifting the corresponding biases, we do not need to update the thresholds further.

We report here the result of the mean-field method for computing data-dependent statistics, which are averaged over a random batch of 500 images. For the data-independent statistics, 500 persistent chains are run in parallel with samples collected after every 5 Gibbs steps. The sparsity level ρ is set to 0 3 and the sparseness weight . λ is set to 0 5. . Once the model has been learned, mean-field is used to estimate the hidden posteriors. Typically this mean-field is quite fast as it converges in a few steps.

We take the data from MNIST and binarize the images using a mid-intensity threshold. The learned representation is shown in Figure 2. Most digits are well separated in 2D except for digits 4 and 9. The learned representation can be used for classifications, e.g., by feeding to the multiclass logistic classifier. For 500 hidden units, the Probit RBM achieves the error rate of 3 28%, comparable with those obtained by the RBM trained . with CD-1 (3 02%), and much better than the raw pixels (8 46%). The features discovered . . by the Probit RBM and RBM with CD-1 are very different (Figure 1), and this is expected because they operate on different input representations. The energy surface learned by the Probit RBM is smooth enough to allow efficient exploration of modes, as shown in Figure 3.

Figure 1: 100 MNIST feature weights (one image per hidden unit) learned by Probit RBM (left) and RBM with CD-1 (right).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-3c8cbaf6-e610-450a-8f25-8362c78d361d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-f50cddf7-4caf-44f1-bbeb-9376252c4aa5.jpg)

Figure 2: t-SNE visualisation of the learned representations on MNIST (left) and random samples of two confused digits (4,9) (right). Best viewed in colours.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-87079b2e-ca88-4347-b3e8-a93a1d00017c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-59ed65f1-6285-405c-b1af-cab899855f55.jpg)

Figure 3: Samples generated from the Probit RBM (left) and RBM with CD-1 (right) by mode seeking procedures (see Supplement).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-120a30a4-24e0-4ce1-b648-18cb8b46394c.jpg)

## 6.2 Rank Evidences for Collaborative Filtering

In collaborative filtering, one of the main goals is to produce a personalized ranked list of items. Until very recently, the majority in the area, on the other hand, focused on predicting the ratings, which are then used for ranking items. It can be arguably more efficient to learn a rank model directly instead of going through the intermediate steps.

We build one TBM for ranking with ties (i.e., inequality constraints, see Section 3.2) per user due to the variation in item choices but all the TBMs share the same parameter set. The handling of ties is necessary because during training, many items share the same rating. Unseen items are simply not accounted for in each model: We only need to compare the utilities between the items seen by each person. The result is that the models are very sparse and fast to run. For the data-dependent statistics, we maintain one Markov chain per user. Since there is no single model for all data instances, the dataindependent statistics cannot be estimated from a small set of Markov chains. Rather we also maintain a data-independent chain per data instance, which can be persistent on their own, or restarted from the data-dependent chains after every parameter updating step. The latter case, which is reported here, is in the spirit of the Hinton's Contrastive Divergence, where the data-independent chain is just a few steps away from the datadependent chain.

Once the model has been trained, the hidden posterior vector ˆ h = ( ˆ h , h 1 ˆ 2 , ..., h ˆ K ) , where ˆ h k = P h ( k = 1 | e u ), is used as the new representation of the tastes of user u . The rank of unseen movies is the mode of the distribution P ( e ∗ | ˆ ), h where e ∗ are the rankbased evidences (see Section 4.4). For fast computation of P ( e ∗ | ˆ ), we approximate the h Gaussian by a Gumbel distribution, which leads to a simple way of ranking movies using the mean 'utility' µ ui = α i + W i · ˆ h u for user u (see the Supplement for more details).

The data used in this experiment is the MovieLens, which contains 1M ratings by approximately 6K users on 4K movies. To encourage diversity in the rank lists, we remove the top 10% most popular movies. We then remove users with less than 30 ratings on the

Table 1: Item ranking results on MovieLens - the higher the better ( K = 50). Here N@T is a shorthand for NDCG@T .

|                 |   ERR |   N@1 |   N@5 |   N@10 |
|-----------------|-------|-------|-------|--------|
| Item popularity | 0.587 | 0.56  | 0.68  |  0.835 |
| ListRank.MF     | 0.653 | 0.673 | 0.751 |  0.873 |
| PMOP            | 0.648 | 0.664 | 0.747 |  0.871 |
| TBM             | 0.678 | 0.722 | 0.792 |  0.893 |

remaining movies. The most recently rated 10 movies per user are held out for testing, the next most recent 5 movies are used for tuning hyper-parameters, and the rest for training.

For comparison, we implement a simple baseline using item popularity for ranking, and thus offering a naive non-personalized solution. For personalized alternatives, we implement two recent rank-based matrix factorisation methods, namely ListRank.MF [27] and PMOP [36]. Two ranking metrics from the information retrieval literature are used: the ERR [3] and the NDCG@T [13]. These metrics place more emphasis on the top-ranked items. Table 1 reports the movie ranking results on test subset (each user is presented with a ranked list of unseen movies), demonstrating that the TBM is a clear winner in all metrics.

## 6.3 Mixed Evidences for World Attitude Analysis

Finally, we demonstrate the TBM on mixed evidences. The data is from the survey analysis domain, which mostly consists of multiple questions of different natures such as basic facts (e.g., ages and genders) and opinions (e.g., binary choices, single choices, multiple choices, ordinal judgments, preferences and ranks). The standard approach to deal with such heterogeneity is to perform the so-called 'coding', which converts types into some numerical representations (e.g., ordinal scales into stars, ranks into multiple pairwise comparisons) so that standard processing tools can handle. However, this coding process breaks the structure in the data and thus significant information will be lost. Thus our TBM offers a scalable and generic machinery to process the data in its native format and then convert the mixed types into a more homogeneous posterior vector.

We use the global attitude survey dataset collected by the PewResearch Centre 7 . The survey was conducted on 24 717 people from 24 countries during the period of March 17 , - April 21, 2008 on a variety of topics concerning people's life, opinions on issues in their countries and around the world as well as future expectations. There are 52 binary, 124 categorical (of variable category sizes), 3 continuous, 165 ordinal (of variable level sizes) question types.

Like the case of collaborative filtering, we build one TBM per respondent due to the variation in questions and answers but all the TBMs share the same parameter set. Unanswered/inappropriate questions are ignored. For each respondent, we maintain 2 persistent and non-interacting Markov chains for the data-dependent statistics and the data-independent statistics, respectively.

Figure 4 shows the 2D distribution of respondents from 24 countries obtained by feeding the posteriors to the t-SNE [38] (here no explicit information of countries is used).

7 The datasets are publicly available from http://pewresearch.org/

Figure 4: Distribution of global attitudes obtained by projecting the hidden posteriors (100 hidden units) on 2D using t-SNE. A dot represents one respondent. Best viewed in colours.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0055v1-figure-ae2263fc-d833-42c7-90df-7e6a3f79530b.jpg)

It is interesting to see the cultural/social clustering and gaps between countries as opposed to the geographical distribution (e.g., between Indonesia and Egypt, Australia and UK and the relative separation of the China, Pakistan, Turkey and the US from the rest). To predict the 24 countries, we feed the posteriors into the standard multiclass logistic regression and achieve an error rate of 0 49%, suggesting that the TBM has captured the . intrastate regularities and separated the interstate variations well.

## 7 Related Work

Latent multivariate Gaussian variables have been widely studied in statistical analysis, initially to model correlated binary data 8 [1, 4] then now used for a variety of data types such as ordered categories [15], unordered categories [43], and the mixture of types [6] . Learning with the underlying Gaussian model is notoriously difficult for large-scale setting: independent sampling costs cubic time due to the need of inverting the covariance matrix, while MCMC techniques such as Gibbs sampling can be very slow if the graph is dense and the interactions between variables are strong. This can be partly overcome by adding one more layer of latent variables as in factor analysis [39, 14] and probabilistic principle component analysis [33]. The main difference from our TBM is that those models are directed with continuous factors while ours is undirected with binary factors.

Gaussian RBMs have been used for modelling continuous data such as visual features [12], where the evidences are the value assignments, and thus a limiting case of our evidence system. Some restrictions to the continuous Boltzmann machines have been studied: In [5], Gaussian variables are assumed to be non-negative, and in [41], continuous variables are bounded. However, we do not make these restrictions on the model but rather placing restrictions during the training phase only. GRBMs that handle ordinal evidences have been studied in [35], which is an instance of the boxed-constraints in our TBM.

## 8 Discussion and Conclusion

Since the underlying variables of the TBM are Gaussian, various extensions can be made without much difficulty. For example, direct correlations among variables, regardless of their types , can be readily modelled by introducing the non-identity covariance matrix [22]. This is clearly a good choice for image modelling since nearby pixels are strongly correlated. Another situation is when the input units are associated with their own attributes. Each unit can be extended naturally by adding a linear combination of attributes to the mean structure of the Gaussian.

The additive nature of the mean-structure allows the natural extension to matrix modelling (e.g., see [37, 35]). That is, we do not distinguish the role of rows and columns, and thus each row and column can be modelled using their own hidden units (the row parameters and columns parameters are different). Conditioned on the row-based hidden units, we return to the standard TBM for column vectors. Inversely, conditioned on the column-based hidden units, we have the TBM for row vectors.

To sum up, we have proposed a generic class of models called Thurstonian Boltzmann machine (TBM) to unify many type-specific modelling problems and generalise them

8 This is often known as multivariate probit models .

to the general problem of learning from multiple groups of inequalities. Our framework utilises the Gaussian restricted Boltzmann machines, but the Gaussian variables are never observed except for one limiting case. Rather, those variables are subject to inequality constraints whenever an evidence is observed. Under this representation, the TBM supports a very wide range of evidences, many of which were not possible before in the Boltzmann machine literature, without the need to specify type-specific models. In particular, the TBM supports any combination of the point assignments, intervals, censored values, binary, unordered categories, multi-categories, ordered categories, (in)-complete ranks with and without ties.

We demonstrated the TBM on three applications of very different natures, namely handwritten digit recognition, collaborative filtering and complex survey analysis. The results are satisfying and the performance is competitive with those obtained by typespecific models.

## A Supplementary Material

## A.1 Inference

## A.1.1 Estimating the Partition Function

For convenience, let us re-parameterise the distribution as follows

$$\phi _ { i } ( x _ { i } ) \ & = \ \exp \left \{ - \frac { x _ { i } ^ { 2 } } { 2 } + \alpha _ { i } x _ { i } \right \} & & ( 9 ) \\ \psi _ { i k } ( x _ { i }, h _ { k } ) \ & = \ \exp \left \{ W _ { i k } x _ { i } h _ { k } \right \} & & \\ \phi _ { k } ( h _ { k } ) \ & = \ \exp \left \{ \gamma _ { k } h _ { k } \right \}$$

The model potential is then the product of all local potentials

$$\Psi ( x, h ) = \left [ \prod _ { i } \phi _ { i } ( x _ { i } ) \right ] \left [ \prod _ { i k } \psi _ { i k } ( x _ { i }, h _ { k } ) \right ] \left [ \prod _ { k } \phi _ { k } ( h _ { k } ) \right ]$$

The partition function can be rewritten as

$$Z \ & = \ \sum _ { h } \int _ { x } \Psi ( x, h ) d x \\ & = \ \sum _ { h } \Omega ( h )$$

where Ω( h ) = GLYPH&lt;1&gt; x Ψ( x h , ) d x . We now proceed to compute Ω( h ):

$$\text{where} \, \Omega ( h ) = \int _ { x } & \Psi ( x, h ) d x. \, \text{We now proceed to compute } \Omega ( h ) \colon \\ \Omega ( h ) & \ = \ \left [ \prod _ { k } \phi _ { k } ( h _ { k } ) \right ] \int _ { x } \left [ \prod _ { i } \phi _ { i } ( x _ { i } ) \right ] \left [ \prod _ { i k } \psi _ { i k } ( x _ { i }, h _ { k } ) \right ] d x \\ & \ = \ \left [ \prod _ { k } \phi _ { k } ( h _ { k } ) \right ] \prod _ { i } \int _ { x _ { i } } \exp \left \{ - \frac { x _ { i } ^ { 2 } } { 2 } + ( \alpha _ { i } + \sum _ { k } W _ { i k } h _ { k } ) x _ { i } \right \} d x _ { i } \\ & \ = \ \left [ \prod _ { k } \phi _ { k } ( h _ { k } ) \right ] \prod _ { i } C _ { i } \int _ { x _ { i } } \exp \left \{ - \frac { ( x _ { i } - \mu _ { i } ( h ) ) ^ { 2 } } { 2 } \right \} d x _ { i } \\ & \ = \ \left [ \prod _ { k } \phi _ { k } ( h _ { k } ) \right ] \prod _ { i } C _ { i } \sqrt { 2 \pi \sigma _ { i } ^ { 2 } } \\ \text{where}$$

where

$$\mu _ { i } ( h ) \ & = \ \alpha _ { i d } + \sum _ { k = 1 } ^ { K } W _ { i k } h _ { k } \\ C _ { i } \ & = \ \exp \left \{ \frac { 1 } { 2 } \left ( \frac { \mu _ { i } ( h ) } { \sigma _ { i } } \right ) ^ { 2 } \right \}$$

Now we can define the distribution over the hidden layer as follows

$$P ( h ) = \frac { 1 } { Z } \Omega ( h )$$

Now we apply the Annealed Importance Sampling (AIS) [19]. The idea is to introduces the notion of inverse-temperature τ into the model, i.e., P ( h | τ ) ∝ Ω( h ) τ .

Let { τ s } S s =0 be the (slowly) increasing sequence of temperature, where τ 0 = 0 and τ S = 1, that is τ 0 &lt; τ ... 1 &lt; τ S . At τ 0 = 0, we have a uniform distribution, and at τ S = 1, we obtain the desired distribution. At each step s , we draw a sample h s from the distribution P ( h | τ s -1 ) (e.g. using some Metropolis-Hastings procedure). Let P ∗ ( h | τ ) be the unnormalised distribution of P ( h | τ ), that is P ( h | τ ) = P ∗ ( h | τ ) /Z τ ( ). The final weight after the annealing process is computed as

$$\omega = \frac { P ^ { * } ( h ^ { 1 } | \tau _ { 1 } ) } { P ^ { * } ( h ^ { 1 } | \tau _ { 0 } ) } \frac { P ^ { * } ( h ^ { 2 } | \tau _ { 2 } ) } { P ^ { * } ( h ^ { 2 } | \tau _ { 1 } ) } \dots \frac { P ^ { * } ( h ^ { S } | \tau _ { S } ) } { P ^ { * } ( h ^ { S } | \tau _ { S - 1 } ) }$$

The above procedure is repeated T times. Finally, the normalisation constant at τ = 1 is computed as Z (1) ≈ Z (0) ( ∑ T t =1 ω ( ) t /T ) where Z (0) = 2 K , which is the number of configurations of the hidden variables h .

## A.1.2 Estimating Posteriors using Mean-field

Recall that for evidence e we want to estimate posteriors P ( h | e ) = ∑ x ∈ Ω( ) e P Ω( ) e ( h x , | e ). Assume that the evidences can be expressed in term of boxed constraints, which lead

to the following factorisation

$$P ( x \, | \, e, h ) = \prod _ { i } P ( x _ { i } \, | \, e, h )$$

This factorisation is critical because it ensures that there are no deterministic constraints among { x i } n i =1 , which are the conditions that variational methods such as mean-fields would work well. This is because mean-field solution will generally not satisfy deterministic constraints, and thus may assign non-zeros probability to improbably areas.

To be more concrete, the mean-field approximation would be Q ( h x , ) ≈ P ( h x , | e )

$$Q ( h, x ) \ & = \ \prod _ { k } Q _ { k } ( h _ { k } ) \prod _ { i } Q _ { i } ( x _ { i } ) \\ \text{s.t.} \ \ x \ \in \ \Omega ( e )$$

The best mean-field approximation will be the minimiser of the Kullback-Leibler divergence

$$\mathcal { D } \left ( Q \| P \right ) \ & = \ \sum _ { h } \sum _ { x \in \Omega \left ( e \right ) } Q \left ( h, x \right ) \log \frac { Q \left ( h, x \right ) } { P \left ( h, x \, | \, e \right ) } \\ & = \ - H \left [ Q \left ( h, x \right ) \right ] - \sum _ { h } \sum _ { x \in \Omega \left ( e \right ) } Q \left ( h, x \right ) \log P \left ( h, x \, | \, e \right ) \quad \ \ ( 1 1 )$$

where H Q [ ( h x , )] is the entropy function. Now first, exploit the fact that Q is factorisable, and thus its entropy is decomposable, i.e.,

$$H \left [ Q ( h, x ) \right ] = \sum _ { k } H \left [ Q _ { k } ( h _ { k } ) \right ] + \sum _ { i } H \left [ Q _ { i } ( x _ { i } ) \right ]$$

Second recall from Eq. (17) that

$$P ( h, x \, | \, e ) = \frac { 1 } { Z ( e ) } \exp \left \{ - E ( x, h ) \right \}$$

and thus

$$\sum _ { h } \sum _ { x \in \Omega ( e ) } Q ($$

Since log Z ( e ) is a constraint w.r.t. Q ( h x , ), we can safely ignore it here.

Now since E ( x h , ) is decomposable (see Eq. (2)), we have

$$\begin{array} {$$

where

$$E _ { i } ( x _ { i } ) \ & = \ \frac { x _ {$$

Combining this decomposition and Eq. (12), we have completely decomposed the Kullback-Leibler divergence in Eq. (11) into local terms:

$$\mathcal { D } ( Q | | P ) \ = \ & \sum _ { i }$$

Now we wish to minimise the divergence with respect to the local distributions { Q x i ( i ) , Q k ( h k ) } for i = 1 2 , , ..., N and k = 1 2 , , ..., K knowing the proper distribution constraints

$$\int _ { x _ { i } \in \Omega ( e _ { i } ) }$$

By the method of Lagrangian multiplier, we have

$$L ( \lambda ) = \mathcal { D } \left ( Q | | P \right ) + \sum _ { i } \lambda _ { i } \left ( \int _ { x _ { i } \in \Omega ( e _ { i } ) } Q _ { i } ( x _ { i } ) - 1 \right ) + \sum _ { k } \kappa _ { k } \left ( \sum _ { h _ { k } } Q _ { k } ( h _ { k } ) - 1 \right )$$

- · Let us compute the partial derivative w.r.t. Q x i ( i ):

$$\partial _ { Q _ { i } ( x _ { i } ) } L ( \lambda ) \ = \ \log Q _ { i } ( x _ { i } ) + 1 + E _ { i } ( x _ { i } ) + \sum _ { k } Q _ { k } ( h _ { k } ^ { 1 } ) E _ { i k } ( x _ { i }, h _ { k } ^ { 1 } ) + \lambda _ { i }$$

where h 1 k is a short hand for h k = 1 and we have made use of the fact that E ik ( x , h i k = 0) = 0. Setting this gradient to zero yields

$$\begin{array} { r c l } \cdot \cdot \cdot & \cdot & \cdot & \cdot \\ Q _ { i } ( x _ { i } ) & = & \exp \left \{ - \left ( E _ { i } ( x _ { i } ) + \sum _ { k } Q _ { k } ( h _ { k } ^ { 1 } ) E _ { i k } ( x _ { i }, h _ { k } ^ { 1 } ) \right ) - 1 - \lambda _ { i } \right \} \\ & = & \exp \left \{ - \frac { 1 } { 2 } \left ( x _ { i } - \left ( \alpha _ { i } + \sum _ { k } Q _ { k } ( h _ { k } ^ { 1 } ) W _ { i k } \right ) \right ) ^ { 2 } - 1 - \lambda _ { i } \right \} \end{array}$$

for x i ∈ Ω( e i ). Normalising this distribution would lead to the truncated form of the normal distribution those the mean is

$$\mu _ { i } = \alpha _ { i } + \sum _ { k } Q _ { k } ( h _ { k } ^ { 1 } ) W _ { i k } h _ { k }$$

- · In a similar way, the partial derivative w.r.t. Q k ( h k ) would be

$$\partial _ { Q _ { k } ( h _ { k } ) } L ( \lambda ) = \log Q _ { k } ( h _ { k } ) - h _ { k } \left ( \gamma _ { k } + \sum _ { i } W _ { i k } \sum _ { x _ { i } \in \Omega ( e _ { i } ) } Q _ { i } ( x _ { i } ) x _ { i } \right ) + 1 + \kappa _ { k }$$

Equating the gradient to zero, we have

$$Q _ { k } ( h _ { k } ) \infty \exp \left \{ h _ { k } \left ( \gamma _ { k } + \sum _ { i } W _ { i k } \hat { \mu } _ { i } \right ) \right \}$$

where ˆ µ i is the mean of the truncated normal distribution

$$\hat { \mu } _ { i } = \sum _ { x _ { i } \in \Omega ( e _ { i } ) } Q _ { i } ( x _ { i } ) x _ { i }$$

Normalising Q k ( h k ) would lead to

$$Q _ { k } ( h _ { k } ^ { 1 } ) = \left [ 1 + \exp \left \{ - \gamma _ { k } - \sum _ { i } W _ { i k } \hat { \mu } _ { i } \right \} \right ] ^ { - 1 }$$

- · Finally, combining these findings in Eqs. (13,14,15), and letting Ω( e i ) = [ b , c i i ] be the boxed constraint, and using the fact that the mean of the truncated distribution is

$$\hat { \mu } _ { i } = \mu _ { i } + \frac { \phi ( b _ { i } - \mu _ { i } ) - \phi ( c _ { i } - \mu _ { i } ) } { \Phi ( c _ { i } - \mu _ { i } ) - \Phi ( b _ { i } - \mu _ { i } ) }$$

we would arrive at the three recursive equations

$$\begin{matrix} \\ \end{matrix}$$

$$\omega \, \max \, \max \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col.: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colot \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colots \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col on \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colons \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colON \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \color \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colOn \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colun \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col>. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col.. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col { \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \text{ } \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \
 \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \$$

where q k is a short hand for Q k ( h 1 k ) and φ z ( ) is the normal probability density function, and Φ( z ) is the cumulative distribution function ♣

## A.1.3 Seeking Modes and Generating Representative Samples

Once the model has been learned, samples can be generated straightforwardly by first sampling the underlying Gaussian RBM and then collect the true samples that satisfy the inequalities of interest. For example, for binary samples, if the generated Gaussian value for a visible unit is larger than the threshold, then we have an active sample. Likewise, rank samples, we only need to rank the sampled Gaussian values.

However, this may suffer from the poor mixing if we use standard Gibbs sampling, that is the Markov chain may get stuck in some energy traps. To jump out of the trap

we propose to periodically raise the temperature to a certain level (e.g., 10) and then slowly cool down to the original temperature (which is 1). In our experiment, the cooling is scheduled as follows

$$T \leftarrow \eta T$$

where η ∈ (0 1) is estimated so that for , n steps, the temperature will drop from T max to T min . That is, T min = η n T max , leading to η = ( T min /T max ) 1 /n .

To locate a basis of attraction, we can lower the temperature further (e.g., to 0 1) to . trap the particles there. Then we collect k successive samples and take the average to be the representative sample. In our experiments, k = 50.

## A.2 Learning

## A.2.1 Gradient of the Likelihood

The log-likelihood of an evidence is

$$\mathcal { L } \ & = \ \log P ( e ) = \log \sum _ { h } \int _ { \Omega ( e ) } P ( h, x ) d x \\ & = \ \log \sum _ { h } \int _ { \Omega ( e ) } \exp \left \{ - E ( x, h ) \right \} d x - \log Z \\ & = \ \log Z ( e ) - \log Z \\ \mathcal { \Phi } \quad \ \. \quad \. \quad \. \quad \.$$

where Z ( e ) = ∑ h GLYPH&lt;1&gt; Ω ( e ) exp {-E ( x h , ) } d x . The gradient of log Z ( e ) w.r.t. the mapping parameter W ik reads

$$\partial _ { W _ { i k } } \log Z ( e ) \ & = \ \frac { - 1 } { Z ( e ) } \sum _ { h } \int _ { \Omega ( e ) } \exp \left \{ - E ( x, h ) \right \} \partial _ { W _ { i k } } E ( x, h ) d x \\ & = \ - \sum _ { h } \int _ { \Omega ( e ) } P ( x, h \, | \, e ) \partial _ { W _ { i k } } E ( x, h ) d x$$

where we have moved the constant Z -1 ( e ) into the sum and integration and make use of the fact that

$$P ( x, h \, | \, e ) = \frac { P ( x, h ) } { P ( e ) } = \frac { 1 } { Z ( e ) } \exp \left \{ - E ( x, h ) \right \}$$

where the domain of the Gaussian is constrained to x ∈ Ω ( e ).

From the definition of the energy function in Eq. (2), we know that the energy is decomposable, and thus the gradient w.r.t. W ik only involves the pair ( x , h i k ). In particular

$$\partial _ { W _ { i k } } E ( x, h ) = - x _ { i } h _ { k }$$

This simplifies Eq. (16)

$$\partial _ { W _ { i k } } \log Z ( e ) \ & = \ - \sum _ { h _ { i } } \int _ { \Omega ( e _ { i } ) } P ( x _ { i }, h _ { k } \ | \ e ) \partial _ { W _ { i k } } E ( x, h ) d x _ { i } \\ & = \ \mathbb { E } _ { P ( x, h | e ) } \left [ x _ { i } h _ { k } \right ]$$

A similar process would lead to

$$\partial _ { W _ { i k } } \log Z = \mathbb { E } _ { P ( x _ { i }, h _ { k } ) } \left [ x _ { i } h _ { k } \right ]$$

and finally:

$$\partial _ { W _ { i k } } \mathcal { L } \ = \ \mathbb { E } _ { P ( x _ { i }, h _ { k } | e ) } \left [ x _ { i } h _ { k } \right ] - \mathbb { E } _ { P ( x _ { i }, h _ { k } ) } \left [ x _ { i } h _ { k } \right ] \bullet$$

## A.2.2 Regularising the Markov Chains

One undesirable feature of the MCMC chains used in learning we have experiences so far is the tendency for the binary hidden states to get stuck, i.e., after some point they do not flip their assignments as learning progresses. We conjecture that this phenomenon may be due to the saturation effect inherent in the factor posterior:

$$P ( h _ { k } = 1 \ | \ x ) = \frac { 1 } { 1 + \exp \left ( - \gamma _ { k } - \sum _ { i } W _ { i k } x _ { i } \right ) }$$

i.e., once the collected value to a node ( γ k + ∑ i W x ik i ) is too high or too low, it is very hard to over turn.

Fortunately, there is a known technique to regularise the chain: we enforce that at a time, there should be only a fraction ρ of nodes which are active, where ρ ∈ (0 1). , One way is to maximise the following objective function

$$\mathcal { L } _ { 2 } \ = \ \mathcal { L } + \lambda \int \left [ \sum _ { k } \sum _ { h _ { k } } \rho ( h _ { k } ) \log P ( h _ { k } \, | \, x ) \right ] P ( x \, | \, e ) d x$$

where λ &gt; 0 is the weighting factor, ρ h ( k ) = ρ if h k = 1 and ρ h ( k ) = 1 -ρ otherwise. The gradient with respect to g k = ( γ k + ∑ i W x ik i ) is then

$$\partial _ { g _ { k } } \mathcal { L } _ { 2 } \ & = \ \partial _ { g _ { k } } \mathcal { L } + \lambda \int \left [ \partial _ { g _ { k } } \sum _ { h _ { k } } \rho ( h _ { k } ) \log P ( h _ { k } \, | \, x ) \right ] P ( x | e ) d x \\ & = \ \partial _ { g _ { k } } \mathcal { L } + \lambda \int \left [ \rho - P ( h _ { k } ^ { 1 } \, | \, x ) \right ] P ( x | e ) d x \\ & \approx \ \partial _ { g _ { k } } \mathcal { L } + \frac { 1 } { S } \lambda \sum _ { s } \left [ \rho - P ( h _ { k } ^ { 1 } \, | \, x ^ { ( s ) } ) \right ] \\ \intertext { h o w o d o n d o n d } \, \mathcal { G } \, \intertext { h o w o d o n d o n d } \, \mathcal { D } ( h ^ { 1 } + \mathcal { G } ) \, \intertext { h o w o d o n d o n d o n d } \, \mathcal { D } ( h ^ { 1 } + \mathcal { G } ) \, \intertext { H o w o d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d }$$

where S is the number of samples and P h ( 1 k | x ) is a shorthand for P h ( k = 1 | x ). Using the chain rule, we have:

$$\partial _ { \gamma _ { k } } \mathcal { L } _ { 2 } \ \ & \approx \ \partial _ { \gamma _ { k } } \mathcal { L } + \frac { 1 } { S } \lambda \sum _ { s } \left [ \rho - P ( h _ { k } ^ { 1 } \, | \, x ^ { ( s ) } ) \right ] \\ \partial _ { W _ { i k } } \mathcal { L } _ { 2 } \ \ & \approx \ \partial _ { w _ { i k d } } \mathcal { L } + \frac { 1 } { S } \lambda \sum _ { s } x _ { i } \left [ \rho - P ( h _ { k } ^ { 1 } \, | \, x ^ { ( s ) } ) \right ]$$

## A.2.3 Online Estimation of Posteriors

For tasks such as data completion (e.g., collaborative filtering) we need the posteriors P ( h | e ) for the prediction phase. One way is to run the Markov chain or doing meanfield from scratch. Here we suggest a simple way to obtain an approximation directly from the training phase without any further cost. The idea is to update the estimated posterior ˆ at each learning step h t in an exponential smoothing fashion:

$$\hat { h } ^ { ( t ) } \leftarrow \eta \hat { h } ^ { ( t - 1 ) } + ( 1 - \eta ) \hat { h } ^ { ( t ) }$$

for some smoothing factor η ∈ (0 1) and initial , ˆ h (0) , where ¯ h ( ) t k = P h ( 1 k | x ( ) t , e ) and x ( ) t is the sampled Gaussian at time t .

As learning progresses, early samples, which are from incorrect models, will be exponentially weighted down. Typically we choose η close to 1, e.g., η = 0 9. .

## A.2.4 Monitoring the Learning Progress

It is often of practical importance to track the learning progress, either by the reconstruction errors or by the data likelihood. The data likelihood can be estimated as

$$P ( e ) \approx \frac { 1 } { S } \sum _ { s = 1 } ^ { S } \int _ { \Omega ( e ) } P ( x \, | \, h ^ { ( s ) } ) d x$$

where h ( s ) are those samples collected as learning progressed in the data-independent phase, and the integration can be carried out using the technique described in the main text.

## A.3 Extreme Value Distributions

Extreme value distributions are a class of distributions of extremal measurements [11]. Here we are concerned about the popular Gumbel's distribution.

## A.3.1 Gumbel Distribution for Categorical Choices

Let us start from the Gumbel density function

$$P ( x ) = \frac { 1 } { \sigma } \exp \left \{ - \left ( \frac { x - \mu } { \sigma } + e ^ { - \frac { x - \mu } { \sigma } } \right ) \right \}$$

where µ is the mode (location) and σ is the scale parameter.

Using Laplace's approximation (e.g., via Taylor's expansion of ( x -µ σ + e -x -µ σ ) using the second-order polynomial around µ ), we have

$$P ( x ) \in \frac { 1 } { e \sigma } \exp \left \{ - \frac { ( x - \mu ) ^ { 2 } } { 2 \sigma ^ { 2 } } \right \}$$

Renormalising this distribution, e.g., by replacing e ≈ 2 7183 by . √ 2 π ≈ 2 5066 we obtain . the standard Gaussian distribution. Thus, we can use the Gumbel as an approximation to the Gaussian distribution .

Now we turn to the categorical model using Gumbel variables. We maintain one variable per category, which plays the role of the utility for the category. Assume that all utilities share the same scale parameter σ . The existing literature [18] asserts that the probability of choosing the m -th category is

$$P ( e = c _ { m } ) = \frac { e ^ { \mu _ { m } / \sigma } } { \sum _ { l } e ^ { \mu _ { l } / \sigma } }$$

where µ l is the location of the l -th utility.

glyph[negationslash]

When we choose a the m -th category we must ensure that x m &gt; max l = m x l . Let y l = exp ( -x l -µ l σ ) or equivalently x l = µ l -σ log y l . Thus x l &lt; x m means µ l -σ log y l &lt; µ m -σ log y m , or y l &gt; y m exp ( µ l -µ m σ ) .

The CDF of the l -th Gumbel distribution is

$$\iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iia \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \julia \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \ iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota\iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \jul \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \ jul \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota\ iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \ioda \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \sigma \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \sigma \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi$$

Thus choosing category c m would mean glyph[negationslash]

glyph[negationslash]

$$P ( e = c _ { m } ) \ = \ \int P ( x _ { m } ) \left ( \prod _ { l \neq m } \int ^ { x m } P ( x _ { l } ) d x _ { l } \right ) d x _ { m } \\ = \ \int P ( x _ { m } ) \prod _ { l \neq m } F _ { l } ( x _ { m } ) d x _ { m } \\ \dots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \colon \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdot \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdclots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \end{array}$$

We rewrite the Gumbel density function by changing variable from x m to y m :

$$P ( y _ { m } ) = \frac { 1 } { \sigma } y _ { m } \exp \left \{ - y _ { m } \right \}$$

for y m ≥ 0. Thus glyph[negationslash]

$$P ( x _ { m } ) \prod _ { l \neq m } F _ { l } ( x _ { m } ) = \frac { 1 } { \sigma } y _ { m } \exp \left ( - y _ { m } \left \{ 1 + \sum _ { l \neq m } e ^ { \frac { \mu _ { l } - \mu _ { m } } { \sigma } } \right \} \right )$$

glyph[negationslash]

Now, by changing variable under the integration from x m to y m , we have glyph[negationslash]

glyph[negationslash]

$$\text{by changing variable under the integration from $x_{m}$ to $y_{m}$, we have} \\ P ( e = c _ { m } ) & \ = \ \sigma \int _ { 0 } ^ { \infty } P ( y _ { m } ) \prod _ { l \neq m } F _ { l } ( y _ { m } ) ^ { \frac { 1 } { y _ { l } } } d y _ { m } \\ & = \ \int _ { 0 } ^ { \infty } \exp \left ( - y _ { m } \left \{ 1 + \sum _ { l \neq m } e ^ { \frac { \mu _ { l } - \mu _ { m } } { \sigma } } \right \} \right ) d y _ { m } \\ & = \ \frac { 1 } { 1 + \sum _ { l \neq m } e ^ { \frac { \mu _ { l } - \mu _ { m } } { \sigma } } } \\ & = \ \frac { e ^ { \mu _ { m } / \sigma } } { \sum _ { l } e ^ { \mu _ { l } / \sigma } } \\ & \quad \text{$\in$}$$

glyph[negationslash]

## A.3.2 Gumbel Distribution for Rank

We now extend the case of categorical evidences rank evidences. Again we maintain one Gaussian variable per category. Without loss of generality, for a particular rank π we assume that we must ensure that x 1 &gt; x 2 &gt; ... &gt; x D . This is equivalent to

$$\left [ x _ { 1 } > \max _ { l > 1 } x _ { l } \right ] \cap \left [ x _ { 2 } > \max _ { l > 2 } x _ { l } \right ] \cap \dots \cap [ x _ { D = 1 } > x _ { D } ]$$

The probability of this is essentially

$$P \left ( \{ e _ { m } = m \} _ { m = 1 } ^ { D } \right ) = P ( e _ { 1 } = 1 ) \prod _ { m > 2 } P \left ( e _ { l } = m \ | \ \{ e _ { d } = l \} _ { l = 1 } ^ { m - 1 } \right )$$

In words, this offers a stagewise process to rank categories: first we pick the best category, the pick the second best from the remaining categories and so on (see also [7]). The probability of picking the best category out of a subset is already given in Appendix A.3.1:

$$P ( e _ { 1 } = 1 ) \ & = \ \frac { e ^ { \mu _ { 1 } / \sigma } } { \sum _ { l \geq 1 } e ^ { \mu _ { l } / \sigma } } \\ P \left ( e _ { l } = m \ | \ \{ e _ { d } = l \} _ { l = 1 } ^ { m - 1 } \right ) \ & = \ \frac { e ^ { \mu _ { m } / \sigma } } { \sum _ { l \geq m } e ^ { \mu _ { l } / \sigma } } \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

This gives us the Plackett-Luce model [17, 21] as mentioned in [30].

## A.4 Global Attitude: Sample Questions

- · Q4 ( Ordinal ): [...] how would you describe the current economic situation in (survey country) -{ very good, somewhat good, somewhat bad, or very bad } ?
- · Q11a ( Binary ): How do you think people in other countries of the world feel about China? -{ like, disliked } ?
- · Q35,35a ( Category-ranking ): Which one of the following, if any, is hurting the world's environment the most/second-most { India, Germany, China, Brazil, Japan, United States, Russia, Other } ?
- · Q76 ( Continuous ): How old were you at your last birthday?
- · Q85 ( Categorical ): What is your current employment situation { A list of employment categories } ?

## A.5 Other Supporting Materials

## A.5.1 Laplace Approximation

Laplace approximation is the technique using a Gaussian distribution to approximate another distribution. For the univariate case, assume that the original density distribution has the form

$$P ( x ) \in \exp \left \{ - f ( x ) \right \}$$

First we find the mode µ of P x ( ) or equivalently the minimiser of f ( x ) given it exists. Then we apply Taylor's expansion

$$f ( x ) \approx f ( \mu ) + f ^ { \prime \prime } ( \mu ) \frac { \left ( x - \mu \right ) ^ { 2 } } { 2 }$$

The Gaussian approximation has the form

$$P ^ { * } ( x ) \in \exp \left \{ - f ^ { \prime \prime } ( \mu ) \frac { \left ( x - \mu \right ) ^ { 2 } } { 2 } \right \}$$

where 1 /f ′′ ( µ ) is the new variance.

## A.5.2 Some Properties of the Truncated Normal Distribution

For a normal distribution P x µ, σ ( | ) of mean µ and standard deviation σ truncated from both sides, i.e., α &lt; x &lt; β , the new density reads

$$\bar { P } _ { [ \alpha, \beta ] } ( x \, |, \mu, \sigma ) = \frac { Q ( x ^ { * } ) } { \sigma \left [ \Phi ( \beta ^ { * } ) - \Phi ( \alpha ^ { * } ) \right ] }$$

where

$$x ^ { * } \ = \ \frac { x - \mu } { \sigma } ; \ \alpha ^ { * } = \frac { \alpha - \mu } { \sigma } ; \ \beta ^ { * } = \frac { \beta - \mu } { \sigma }$$

and Q ( ) · and Φ( ) are the probability density function and the cumulative distribution · of the standard normal distribution, respectively. In particular, we are interested in the mean of the distribution ¯ P [ α,β ] :

$$\bar { \mu } = \mu _ { i } + \sigma \frac { Q ( \alpha ^ { * } ) - Q ( \beta ^ { * } ) } { \Phi ( \beta ^ { * } ) - \Phi ( \alpha ^ { * } ) }$$

Some special cases:

- · When α = β , this distribution reduces to the Dirac's delta.
- · When α = -∞ , we have a one-sided truncation from above since Φ( α ∗ ) = 0.
- · When β = + ∞ , we obtain a one-sided truncation form below since Φ( β ∗ ) = 0.

## References

- [1] JR Ashford and RR Sowden. Multi-variate probit analysis. Biometrics , pages 535546, 1970. 7
- [2] U. B¨ ockenholt. Thurstonian-based analyses: past, present, and future utilities. Psychometrika , 71(4):615-629, 2006. 1
- [3] O. Chapelle, D. Metlzer, Y. Zhang, and P. Grinspan. Expected reciprocal rank for graded relevance. In CIKM , pages 621-630. ACM, 2009. 6.2

- [4] S. Chib and E. Greenberg. Analysis of multivariate probit models. Biometrika , 85(2):347-361, 1998. 7
- [5] O.B. Downs, D.J.C. MacKay, and D.D. Lee. The nonnegative Boltzmann machine. Advances in Neural Information Processing Systems , 12:428-434, 2000. 7
- [6] D.B. Dunson and A.H. Herring. Bayesian latent variable models for mixed discrete outcomes. Biostatistics , 6(1):11, 2005. 7
- [7] M.A. Fligner and J.S. Verducci. Multistage ranking models. Journal of the American Statistical Association , 83(403):892-901, 1988. A.3.2
- [8] Y. Freund and D. Haussler. Unsupervised learning of distributions on binary vectors using two layer networks. Advances in Neural Information Processing Systems , pages 912-919, 1993. 3.1
- [9] P.V. Gehler, A.D. Holub, and M. Welling. The rate adapting Poisson model for information retrieval and object recognition. In Proceedings of the ICML , pages 337-344, 2006. 1
- [10] J. Geweke. Efficient simulation from the multivariate normal and student-t distributions subject to linear constraints and the evaluation of constraint probabilities. In Computing science and statistics: Proceedings of the 23rd symposium on the interface , pages 571-578, 1991. 4.2
- [11] EJ Gumbel. Statistical of extremes . Columbia University Press, New York, 1958. A.3
- [12] G.E. Hinton and R.R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science , 313(5786):504-507, 2006. 1, 7
- [13] K. J¨rvelin and J. Kek¨l¨inen. a a a Cumulated gain-based evaluation of IR techniques. ACM Transactions on Information Systems (TOIS) , 20(4):446, 2002. 6.2
- [14] E. Khan, B. Marlin, and K. Murphy. Variational bounds for mixed-data factor analysis. In Proc. of Neural Information Processing Systems , 2010. 7
- [15] A. Kottas, P. M¨ller, and F. Quintana. Nonparametric Bayesian modeling for multiu variate ordinal data. Journal of Computational and Graphical Statistics , 14(3):610625, 2005. 7
- [16] N. Le Roux, N. Heess, J. Shotton, and J. Winn. Learning a generative model of images by factoring appearance and shape. Neural Computation , 23(3):593-650, 2011. 1
- [17] R.D. Luce. Individual choice behavior . Wiley New York, 1959. A.3.2
- [18] D. McFadden. Conditional logit analysis of qualitative choice behavior. Frontiers in Econometrics , pages 105-142, 1973. A.3.1
- [19] R.M. Neal. Annealed importance sampling. Statistics and Computing , 11(2):125-139, 2001. A.1.1

| [20]   | J. Ngiam, A. Khosla, M. Kim, J. Nam, H. Lee, and A.Y. Ng. Multimodal deep learning. In ICML , 2011. 1                                                                                                    |
|--------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [21]   | R.L. Plackett. The analysis of permutations. Applied Statistics , pages 193-202, 1975. A.3.2                                                                                                             |
| [22]   | M.A. Ranzato and G.E. Hinton. Modeling pixel means and covariances using fac- torized third-order Boltzmann machines. In CVPR , pages 2551-2558. IEEE, 2010. 8                                           |
| [23]   | C.P. Robert. Simulation of truncated normal variables. Statistics and computing , 5(2):121-125, 1995. 4.1                                                                                                |
| [24]   | R. Salakhutdinov and G. Hinton. Deep Boltzmann Machines. In Proceedings of 20th AISTATS , volume 5, pages 448-455, 2009. 1                                                                               |
| [25]   | R. Salakhutdinov and G. Hinton. Replicated softmax: an undirected topic model. Advances in Neural Information Processing Systems , 22, 2009. 1                                                           |
| [26]   | R. Salakhutdinov, A. Mnih, and G. Hinton. Restricted Boltzmann machines for collaborative filtering. In Proceedings of the 24th ICML , pages 791-798, 2007. 3.2                                          |
| [27]   | Y. Shi, M. Larson, and A. Hanjalic. List-wise learning to rank with matrix fac- torization for collaborative filtering. In ACM RecSys , pages 269-272. ACM, 2010. 6.2                                    |
| [28]   | P. Smolensky. Information processing in dynamical systems: Foundations of har- mony theory. Parallel distributed processing: Explorations in the microstructure of cognition , 1:194-281, 1986. 1, 3.1   |
| [29]   | N. Srivastava and R. Salakhutdinov. Multimodal learning with deep Boltzmann machines. In NIPS , pages 2231-2239, 2012. 1                                                                                 |
| [30]   | H. Stern. Models for distributions on permutations. Journal of the American Sta- tistical Association , 85(410):558-564, 1990. A.3.2                                                                     |
| [31]   | L.L. Thurstone. A law of comparative judgment. Psychological review , 34(4):273, 1927. 1                                                                                                                 |
| [32]   | T. Tieleman. Training restricted Boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th ICML , pages 1064-1071, 2008. 5.1, 5.1                                   |
| [33]   | M.E. Tipping and C.M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B , 61(3):611-622, 1999. 7                                                    |
| [34]   | T. Tran, D.Q. Phung, and S. Venkatesh. Mixed-variate restricted Boltzmann ma- chines. In Proc. of 3rd Asian Conference on Machine Learning (ACML) , Taoyuan, Taiwan, 2011. 1                             |
| [35]   | T. Tran, D.Q. Phung, and S. Venkatesh. Cumulative restricted Boltzmann machines for ordinal matrix data analysis. In Proc. of 4th Asian Conference on Machine Learning (ACML) , Singapore, 2012. 1, 7, 8 |

| [36]   | T. Truyen, D.Q Phung, and S. Venkatesh. Probabilistic models over ordered parti- tions with applications in document ranking and collaborative filtering. In Proc. of SIAM Conference on Data Mining (SDM) , Mesa, Arizona, USA, 2011. SIAM. 6.2   |
|--------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [37]   | T.T. Truyen, D.Q. Phung, and S. Venkatesh. Ordinal Boltzmann machines for collab- orative filtering. In Twenty-Fifth Conference on Uncertainty in Artificial Intelligence (UAI) , Montreal, Canada, June 2009. 1, 3.1, 8                           |
| [38]   | L. van der Maaten and G. Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research , 9(2579-2605):85, 2008. 6.3                                                                                                                   |
| [39]   | M. Wedel and W.A. Kamakura. Factor analysis with (mixed) observed and latent variables in the exponential family. Psychometrika , 66(4):515-530, 2001. 7                                                                                           |
| [40]   | E. Xing, R. Yan, and A.G. Hauptmann. Mining associated text and images with dual-wing harmoniums. In Proceedings of the 21st UAI , 2005. 1                                                                                                         |
| [41]   | M. Yasuda and K. Tanaka. Boltzmann machines with bounded continuous random variables. Interdisciplinary Information Sciences , 13(1):25-31, 2007. 7                                                                                                |
| [42]   | L. Younes. Parametric inference for imperfectly observed Gibbsian fields. Probability Theory and Related Fields , 82(4):625-645, 1989. 5.1, 5.1                                                                                                    |
| [43]   | X. Zhang, W.J. Boscardin, and T.R. Belin. Bayesian analysis of multivariate nominal measures using multivariate multinomial probit models. Computational statistics & data analysis , 52(7):3697-3708, 2008. 7                                     |
Truyen Tran Dinh Phung Svetha Venkatesh
Thurstonian Boltzmann Machines: Learning from Multiple Inequalities
We introduce Thurstonian Boltzmann Machines (TBM), a unified architecture
that can naturally incorporate a wide range of data inputs at the same time.
Our motivation rests in the Thurstonian view that many discrete data types can
be considered as being generated from a subset of underlying latent continuous
variables, and in the observation that each realisation of a discrete type
imposes certain inequalities on those variables. Thus learning and inference in
TBM reduce to making sense of a set of inequalities. Our proposed TBM naturally
supports the following types: Gaussian, intervals, censored, binary,
categorical, muticategorical, ordinal, (in)-complete rank with and without
ties. We demonstrate the versatility and capacity of the proposed model on
three applications of very different natures; namely handwritten digit
recognition, collaborative filtering and complex social survey analysis.