Year: 2019

Abstract:
Most deep learning models rely on expressive high-dimensional representations to achieve good performance on tasks such as classification. However, the high dimensionality of these representations makes them difficult to interpret and prone to over-fitting. We propose a simple, intuitive and scalable dimension reduction framework that takes into account the soft probabilistic interpretation of standard deep models for classification. When applying our framework to visualization, our representations more accurately reflect inter-class distances than standard visualization techniques such as t-SNE. We experimentally show that our framework improves generalization performance to unseen categories in zero-shot learning. We also provide a finite sample error upper bound guarantee for the method.

heading: 1 INTRODUCTION
text: Dimensionality reduction is an important problem in machine learning tasks to increase classification performance of learned models, improve computational efficiency, or perform visualization. In the context of visualization, high-dimensional representations are typically converted to two or threedimensional representations so that the underlying relations between data points can be observed and interpreted from a scatterplot. Currently, a major source of high-dimensional representations that machine learning practitioners have trouble understanding are those generated by deep neural networks. Techniques such as PCA or t-SNE (Van der Maaten, 2014) are typically used to visualize them, e.g., in Law et al. (2017); Snell et al. (2017). Moreover, Schulz et al. (2015) proposed a visualization technique that represents examples based on their predicted category only. However, none of these techniques exploit the fact that deep models have soft probabilistic interpretations. For instance, the output of deep classifiers typically employs softmax regression, which optimizes classification scores across categories by minimizing cross entropy. This results in soft probabilistic representations that reflect the confidence of the model in assigning examples to the different categories. Many other deep learning tasks such as semantic segmentation (Long et al., 2015) or boundary/skeleton detection (Xie & Tu, 2015) also optimize for probability distributions. In this paper, we experimentally demonstrate that the soft probability representations learned by a neural network reveal key structure about the learned model. To this end, we propose a dimensionality reduction framework that transforms probability representations into a low-dimensional space for easy visualization.
Furthermore, our approach improves generalization. In the context of zero-shot learning where novel categories are added at test time, deep learning approaches often learn high-dimensional representations that over-fit to training categories. By learning low-dimensional representations that match the classification scores of a high-dimensional pre-trained model, our approach takes into account inter-class similarities and generalizes better to unseen categories than standard approaches.
Proposed approach: We propose to exploit as input representations the probability scores generated by a high-dimensional pre-trained model, called the teacher model or target, in order to train a
lower-dimensional representation, called the student. In detail, our approach learns low-dimensional student representations of examples such that, when applying a specific soft clustering algorithm on the student representations, the predicted clustering scores are similar to the target probability scores.
Contributions: This paper makes the following contributions: (1) We propose the first dimensionality reduction approach optimized to consider some soft target probabilistic representations as input. (2) By exploiting the probability representations generated by a pre-trained model, our approach reflects the learned semantic structure better than standard visualization approaches. (3) We experimentally show that our approach improves generalization performance in zero-shot learning. (4) We theoretically analyze the statistical properties of the approach and provide a finite sample error upper bound guarantee for it.

heading: 2 DIMENSIONALITY REDUCTION FOR PROBABILISTIC REPRESENTATIONS
text: (DRPR)
Our method, called Dimensionality Reduction of Probabilistic Representations (DRPR, pronounced “Derper”), is given probability distribution representations generated from high-dimensional data as target. Its goal is to learn a low-dimensional representation such that the soft clustering scores predicted by a soft clustering algorithm are similar to the target. If the targets are probability distributions generated by a pre-trained classifier, we want the low-dimensional space to reflect the relationships between categories interpreted by the classifier. The position of each example in the low-dimensional space should then reflect the ambiguity of the classifier for the example (see Fig. 1). We summarize in Section 2.1 the soft clustering algorithm that is used by DRPR in the low-dimensional space, the algorithm is detailed in Banerjee et al. (2005). The general learning algorithm of DRPR is introduced in Section 2.2.

heading: 2.1 FAMILY OF EFFICIENT SOFT CLUSTERING ALGORITHMS
text: Probability density: We consider that we are given a set of n vectors f1, · · · , fn ∈ V concatenated into a single matrixF = [f1, · · · , fn]> ∈ Vn. In the following, we consider V = Rd, and Vn = Rn×d. The goal is to partition n examples into k soft clusters. Each cluster Cc with c ∈ {1, · · · , k} has a center µ̂c ∈ V and its corresponding probability density is pc(fi) = exp(−d(fi, µ̂c)) b(fi), where d is a regular Bregman divergence (Banerjee et al., 2005; Bregman, 1967) and b : V → R+ is a uniquely determined function that depends on d and ensures that the integral of the density over V is 1 (e.g., b(fi) = √ 1/(2π)d if d is the squared Euclidean distance). For simplicity, we consider that d is the squared Euclidean distance. The density pc(fi) decreases as the divergence between the example fi and the center µ̂c increases.
Bregman Soft Clustering problem (BSCP): The BSCP (Banerjee et al., 2005) is defined as that of learning the maximum likelihood parameters Γ = {µ̂c, π̂c}kc=1 of a mixture model p(fi|Γ) =∑k c=1 π̂c exp(−d(fi, µ̂c)) b(fi) where π̂c is the prior probability that fi is generated by Cc. To partition the n examples into k clusters, Banerjee et al. (2005) apply the EM algorithm to maximize the likelihood parameter estimation problem for mixture models formulated as: maxΓ ∏n i=1 p(fi|Γ)
Assignment matrix: Partitioning the n observations in F into k soft clusters is equivalent to determining some soft assignment matrix Ŷ ∈ (0, 1)n×k in the set Yn×k of matrices whose rows are positive and sum to 1. Formally, Yn×k is written Yn×k := {Ŷ ∈ (0, 1)n×k : Ŷ 1k = 1n} where 1k ∈ {1}k is the k-dimensional vector containing only 1. For a given value of Γ, the element Ŷic = p(Cc|fi) ∈ (0, 1) is the posterior probability, or responsibility of Cc for fi. The higher the value of Ŷic, the more likely fi belongs to cluster Cc.
Local maximum condition: Once the BSCP has converged to a local maximum of maxΓ ∏n i=1 p(fi|Γ), the following equations are all satisfied:
(E-step) ∀i, c, p(Cc|fi) = Ŷic = π̂c exp(−d(fi, µ̂c))∑k
m=1 π̂m exp(−d(fi, µ̂m)) (1)
(M-step) ∀c, µ̂c = ∑n i=1 Ŷicfi∑n i=1 Ŷic
and ∀c, π̂c = 1
n n∑ i=1 Ŷic (2)
Eq. (1) corresponds to the E-step of the EM algorithm that computes p(Cc|fi) = Ŷic when the parameters F and Γ are given. Eq. (2) corresponds to the M-step, which has a simple form since the likelihood is a regular exponential family function (Banerjee et al., 2005). The M-step may be computationally expensive for other types of exponential family distributions (Banerjee et al., 2005, Section 5.2). It is worth noting that these optimality conditions do not depend on the function b used to define pc(fi), so b can be ignored.

heading: 2.2 LEARNING THE MAPPING TO THE LOW-DIMENSIONAL SPACE
text: Section 2.1 explains how to perform soft clustering on some fixed representation F .
We now describe how to learn F so that the soft clustering scores predicted by the BSCP match those of the target. We assume that we are given the probability representation yi ∈ [0, 1]k as target for each training example fi. These representations are concatenated into a single matrix Y = [y1, · · · ,yn]> ∈ Yn×k which is the target of our method for F . DRPR learns the representation of F ∈ Vn such that the soft assignment matrix obtained from applying the BSCP on F is close to Y . We first give the formulation of our prediction function in Eq. (3). Our dimensionality reduction problem is given in Eq. (4).
Prediction function: Let us assume that we are given the dataset matrixF = [f1, · · · , fn]> ∈ Vn = Rn×d, the centersM = [µ1, · · · ,µk]> ∈ Vk = Rk×d and the priors π = [π1, · · · , πk]> ∈ (0, 1)k. As in Eq. (1), we define our prediction function ψ(F,M,π) = Ψ ∈ Yn×k as the soft assignment matrix predicted by the BSCP given F , M and π. The elements of the matrix Ψ are then computed as follows:
∀i, c, Ψic = πc exp(−d(fi,µc))∑k
m=1 πm exp(−d(fi,µm)) (3)
Optimization problem: DRPR learns the representation of F so that the predicted assignment matrix ψ(F,M,π) = Ψ is similar to Y . Given the optimal condition properties of the BSCP stated in Section 2.1, the optimal values of M and π also depend on Ψ and are therefore variables of our dimensionality reduction problem that we formulate:
min F,M,π ∆n (ψ(F,M,π), Y ) (4)
The function ∆n(Ψ, Y ) is an empirical discrepancy loss between the predicted assignment matrix Ψ and the target assignment matrix Y . Since the rows of Ψ and Y represent probability distributions, we formulate ∆n as the average KL-divergence between the rows of Y and the rows of Ψ. Let ψ>i and y>i be the i-th rows of Ψ and Y , respectively, we formulate ∆n(Ψ, Y ) = 1 n ∑n i=1DKL(yi‖ψi). Note that the choice of the discrepancy loss ∆n is independent of the chosen Bregman divergence d. Moreover, the number of classes k has an impact on the number of clusters in the low-dimensional space but is not related to the dimensionality d of the model. DRPR considers that each class c ∈ {1, · · · , k} is represented by one cluster prototype µc ∈ Rd. Terminology: In our experiments, the target (or teacher) is the assignment matrix Y ∈ Yn×k that contains the probability scores generated by a pre-trained neural network. It corresponds to the output of a classifier trained with softmax regression in the visualization experiments, and to the matrices Y1 and Y2 described in Section 4.2. The goal is then to learn F , M and π so that (1) F , M and π reach the BSCP optimality conditions given in Section 2.1 and (2) Ψ = ψ(F,M,π) is similar to Y .
Visualization: DRPR can be used for visualization since many models (e.g., usual neural networks) have probabilistic interpretations w.r.t. the k training categories. In our visualization task, the matrices M and π are not provided, whereas the target matrix Y is given. By using the optimality conditions
Algorithm 1 Dimensionality Reduction of Probabilistic Representations (DRPR) input : Set of training examples (e.g., images) in X and their target probability scores (e.g., classification scores
w.r.t. k training categories), nonlinear mapping gθ parameterized by parameters θ, number of iterations t 1: for iteration 1 to t do 2: Randomly sample n training examples x1, · · · ,xn ∈ X and create target assignment matrix Y ∈ Yn×k containing the target probability scores y1, · · · ,yn (i.e., Y = [y1, · · · ,yn]> ∈ Yn×k) 3: Create matrix F ← [f1, · · · , fn]> ∈ Vn such that ∀i, fi = gθ(xi) 4: Create matrix of centers M ← diag(Y >1n)−1Y >F and prior vector π ← 1nY
>1n 5: Update the parameters θ by performing a gradient descent iteration of ∆n (ψ(F,M,π), Y ) (i.e., Eq. (4)) 6: end for
output : nonlinear mapping gθ
in Eq. (2), we can write the desired values of M and π as a function of F and Y : at each iteration, for some current value of F , the optimal values M = diag(Y >1n)−1Y >F and π = 1nY
>1n are computed and F is updated via gradient descent. DRPR is illustrated in Algorithm 1 in the case where F is the output of a model gθ parameterized by θ, e.g., gθ can be a neural network. However, we represent F as non-parametric embeddings in our visualization experiments to have an equitable comparison to non-parametric baselines. The learning algorithm then modifies the matrix F at each iteration. If the priors πc are all equal, then the priors are updated in step 4 as follows: π ← 1k1k. Zero-shot learning: DRPR can be used to improve zero-shot learning generalization since highdimensional models may overfit to training categories and the goal of zero-shot learning is to generalize to novel categories. In the considered zero-shot learning task, the variable F concatenates image representations (outputs of a neural network) in the same way as step 3 of Algorithm 1, and the variable M concatenates category representations extracted from text (outputs of another neural network). Both F and M are of different nature and are therefore computed as concatenating the outputs of two distinct neural networks taking different sources as input. To optimize Eq. (4), both neural networks are trained jointly by fixing the other neural network during backpropagation.
Choice of divergence: In our experiments, we consider the squared Euclidean distance d(fi, µ̂c) = ‖fi − µ̂c‖22. However, DRPR can be used with any regular Bregman divergence (Banerjee et al., 2005). The algorithm is then identical, with the exception of the chosen divergence d to compute the prediction in Eq. (3).
Convergence and scalability: Although our problem has 3 variables (F , M and π), we use the optimal properties of the BSCP to write them as a function of each other (see step 4 of Algo 1). Eq. (4) is then an optimization problem w.r.t. only one variable F . Since the problem is differentiable wrt F and unconstrained, it is easy to optimize by gradient descent (e.g., back-propagation when training neural networks). Moreover, our loss is nonnegative, hence lower-bounded. It is worth noting that we do not apply multiple iterations of the EM algorithm at each gradient descent iteration, as we first use the optimal properties of BSCP to obtain closed-form formulations of M and π, and then compute ψ(F,M,π) = Ψ to minimize our problem in Eq. (4). Unlike t-SNE and many iterative DR problems, the complexity of DRPR is linear in n (instead of quadratic) and linear in k, which makes it efficient and scalable. Our visualization experiments take less than 5 minutes to do 105 iterations while t-SNE takes 1 hour to do 1000 iterations. PCA, which has an efficient closed-form solution, is still much faster. Our approach is thus simple to optimize, hence scalable, and it generalizes to a large family of Bregman divergences.
Gradient interpretation: We now interpret the gradient of our optimization problem w.r.t. examples. To simplify its formulation, we consider that all the priors π1, · · · , πk are equal and the matrix M does not depend on F , which is the case in our zero-shot learning task.
When d is the squared Euclidean distance, all the priors are equal and M does not depend on F , the gradient of Eq. (4) w.r.t. fi is:
2
n k∑ c=1 YicΨic[ k∑ m=1 exp(‖fi − µc‖22 − ‖fi − µm‖22)][−µc + k∑ m=1 Ψimµm] (5)
One can observe that its magnitude depends on both the target scores Yic ∈ (0, 1) and the predicted responsibilities Ψic. The gradient tries to make the vector fi closer to each centroid µc while separating it from all the centroids µm depending on their predicted scores Ψim ∈ (0, 1). Statistical guarantee: We analyze the statistical property of the algorithm in Appendix A. Theorem 1 provides a finite sample upper bound guarantee on the quality of the minimizer of the empirical discrepancy of Eq. (4). We show that it is upper bounded by the minimum of the true expected discrepancy, and an estimation error term of O(n−1/2) (under certain conditions and with high probability). We defer the detail to the appendix.

heading: 3 RELATED WORK
text: This paper introduces a dimensionality reduction method that represents the relations in a dataset that has probabilistic interpretation. It can be seen as a metric learning approach.
Metric Learning: The most related approaches try to solve the “supervised” hard clustering problem (Law et al., 2016) . During training, they are given the target hard assignment matrix Y ∈ {0, 1}n×k where Yic = 1 if the training example fi has to belong to Cc, and 0 otherwise. The goal is to learn a representation such that applying a hard clustering algorithm (e.g., kmeans) on the training dataset will return the desired assignment matrix Y . These approaches can be decomposed in 2 groups: (1) the methods that optimize some regression problem (Lajugie et al., 2014; Law et al., 2016; 2017) and exploit orthogonal projection properties that hold only in the hard clustering context, and (2) the methods that use exponential families (Mensink et al., 2012; Snell et al., 2017) to describe some probability score. In the latter case, the learning problem is written as a multiclass logistic regression problem where the probability of a category Cc given an example fi is p(Cc|fi) = exp(−δW (fi,µc))∑k
e=1 exp(−δW (fi,µe)) and δW is the learned dissimilarity function. DRPR generalizes
those approaches and can also be used for hard clustering learning.1 For instance, Snell et al. (2017) is a special case of DRPR where π = 1k1k, Y ∈ {0, 1}
n×k is a hard assignment matrix and ∆n is the same ∆n as ours (i.e., ∆n(Ψ, Y ) = 1n ∑n i=1DKL(yi‖ψi) by using the convention 0 log 0 = 0). Moreover, the optimal value of M is not implicitly written as a function of F and Y in Snell et al. (2017). The approach in Mensink et al. (2012) is similar to Snell et al. (2017) but only considers linear models. In summary, both Snell et al. (2017) and Mensink et al. (2012) do not exploit the BSCP formulation to its full potential as they consider some restricted hard clustering case. DRPR generalizes these approaches to the soft clustering context at no additional algorithmic complexity.
Dimensionality reduction: Learning models by exploiting soft probability scores predicted by another pre-trained model as supervision was proposed in Ba & Caruana (2014) for classification. It was experimentally observed in Ba & Caruana (2014) that using the output of a large pre-trained neural network as supervision, instead of ground truth labels, improves classification performance of small neural networks. However, in Ba & Caruana (2014), each dimension of the student representation describes the confidence score of the model for one training category, which is problematic in contexts such as zero-shot learning where categories can be added or removed at test time. Our approach can learn a representation with dimensionality different from the number of training categories. It can therefore be used in zero-shot learning. Dimensionality reduction with neural networks has been proposed in unsupervised contexts (e.g., to maximize variance in the latent space (Hinton & Salakhutdinov, 2006)) and in supervised contexts (e.g., using ground truth labels as supervision (Salakhutdinov & Hinton, 2007)). Instead, our approach exploits probability scores generated by a teacher pre-trained model. These scores may be very different from ground truth labels if the pretrained model does not generalize well on the dataset. Our approach can then help understand what was learned by the teacher by visualizing groups of examples the teacher has trouble distinguishing.
1We experimentally show in the supplementary material that DRPR can be learned for hard supervised clustering when the centers are implicit (see Section C.3).

heading: 4 EXPERIMENTS
text: We evaluate the relevance of our method in two types of experiments. The first learns low-dimensional representations for visualization to better interpret pre-trained deep models. The second experiment exploits the probability scores generated by a pre-trained classifier in the zero-shot learning context; these probability scores are used as supervision to improve performance on novel categories.

heading: 4.1 VISUALIZATION
text: Interpreting deep models is a difficult task, and one of the most common tools to solve that task is visualization. Representations are most often visualized with t-SNE which does not account for the probabilistic interpretation of the learned models. We propose to exploit probability classification scores as input of our dimensionality reduction framework. In the visualization experiments of this section, DRPR learns non-parametric low-dimensional embeddings (i.e. our representations are not outputs of neural networks but vectors) as done by non-parametric baselines. Nonetheless, DRPR can also be learned with neural networks (e.g., as done in Section 4.2).
Toy dataset: As an illustrative toy experiment, we compare the behavior of t-SNE and DRPR when applied to reduce a simple artificial 3-dimensional dataset as 2D representations. The 3D dataset illustrated in Fig. 2 (a) contains k = 6 clusters, each generated by a Gaussian model and containing 1,000 artificial points. To generate target soft clustering probability scores in the 3D space, we compute the relative distances of the examples to the different cluster centers and normalize them to obtain (probability) responsibilities as done in Eq. (3). In detail, let us note the original 3-dimensional dataset X = [x1, · · · , xn]> ∈ Rn×3 (where n is the number of examples) plotted in Fig. 2 (a). The target soft assignment matrix Y ∈ Yn×k where k = 6 is constructed by computing Yic =
exp(−‖xi−cc‖2)∑k e=1 exp(−‖xi−ce‖2)
where cc ∈ R3 is the center of the c-th Gaussian model (defined by its color) in the original space, and priors are equal. We plot in Fig. 2 the 2D representation of our model when using Y as input/target, and the t-SNE representations obtained when using the original dataset X as input in Fig. 2 (b) and Y as input in Fig. 2 (c). Two observations can be made from Fig. 2: (1) the global structure of the original dataset is better preserved with DRPR than with t-SNE; (2) DRPR satisfies the relative distances between the different clusters better than t-SNE since DRPR tries to preserve the relative responsibilities of the different clusters. We quantitatively evaluate these two observations in the following. It is also worth noting that it is known that distances between clusters obtained with t-SNE may not be meaningful (Wattenberg et al., 2016) as t-SNE preserves local neighborhood instead of global similarities.
In the following (i.e. in Fig. 4 and tables of results), we only consider the case where t-SNE takes as input the logits (i.e., classification scores before the softmax operation) instead of probability scores since the latter case returns bad artifacts such as the one in Fig. 2 (c). Other examples of bad artifacts obtained with t-SNE exploiting probability scores are provided in the supplementary material. We also provide in the supplementary material the visualizations obtained by t-SNE when replacing the `2-norm in the input space by the KL-divergence and the Jensen-Shannon divergence to compare probabilistic representations that DRPR uses as input. These lead to worse visualizations than Fig 2 (b) in terms of preserving original clusters and their inter-cluster distances.
Quantitative evaluation metrics: We quantify the relevance of our approach with the two following evaluation metrics: (1) an adaptation of the Neighborhood Preservation Ratio (NPR) (Van der Maaten & Hinton, 2012): for each image i, it counts the ratio of κ nearest neighbors of i (i.e. that have the closest probability scores w.r.t. the KL divergence) that are also in the set of κ nearest neighbors in the learned low-dimensional space (w.r.t. the Euclidean
distance). It is averaged over all images i. This metric evaluates how much images that have similar probability scores are close to each other with the student representation. (2) Clustering Distance Preservation Ratio (CDPR): we randomly sample 105 triplets of images (i, i+, i−) such that the 3 images all belong to different categories and i has closer probability score to i+ than to i− w.r.t. the KL divergence. The metric counts the percentage of times that the learned representation of i is closer to i+ than to i− w.r.t. the Euclidean distance in the low-dimensional representation. This evaluates how well inter-cluster distances are preserved.
We evaluate our approach on the test sets of the MNIST (LeCun et al., 1998), STL (Coates et al., 2011), CIFAR 10 and CIFAR 100 (Krizhevsky & Hinton, 2009) datasets with pre-trained models that are publicly available and optimized for cross entropy.2 The dimensionality of the high-dimensional representations is equal to the number of categories in the respective datasets (i.e. 10 except for CIFAR 100 that contains 100 categories). Our goal is to visualize the teacher representations with 2-dimensional representations by using their probability scores as target, not the ground truth labels. Quantitative comparisons with standard visualization techniques such as t-SNE, ISOMAP (Tenenbaum et al., 2000), Locally Linear Embedding (LLE) (Roweis & Saul, 2000) and PCA using the 2 leading dimensions are provided in Tables 1 and 2. We also report the scores obtained with the logit representations which are not dimensionality reduction representations but provide an estimate of the behavior of the original dataset. DRPR outperforms the dimensionality reduction baselines w.r.t. both evaluation metrics and is competitive with the logit representation. Examples that have similar probability-based representations are closer with our approach than with other dimensionality reduction baselines. DRPR also better preserves inter-cluster distances. It is worth noting that DRPR exploits as much supervision as the “unsupervised” visualization baselines. Indeed, all the compared methods use as input the same source of supervision which is included in the (classifier output) representations given as input.
Qualitative results: Visualizations of pre-trained models obtained with DRPR and t-SNE are illustrated in Fig. 4 for MNIST and STL. The visualizations for CIFAR 10 and 100 are in the supplementary material. DRPR representations contain spiked groups at the corners to better reflect examples that have high confidence scores for one category. Indeed, an example in a spike at the corner of a figure has a soft assignment score w.r.t. its closest center close to 1. This means that the pre-trained model has very high confidence to assign the example to the corresponding category (see
2We use the pre-trained models available at https://github.com/aaron-xichen/ pytorch-playground.
illustration in Fig. 3). Examples that are between multiple clusters (usually in the middle of figures) are those that are harder to classify by the model. Detailed visualizations and analysis are provided in the supplementary material.
One can observe that representations obtained with DRPR reflect the semantic structure between categories. On MNIST, categories that contain a curve at the bottom of the digit (i.e., 0, 3, 5, 6, 8 and 9) are in the bottom of Fig. 4 (left); some pairs of digits that are often hard to differentiate by classifiers (i.e., 4 and 9, 1 and 7, 3 and 8) are also adjacent. On STL and CIFAR 10, animal categories are illustrated on the right whereas machines are on the left. Semantically close categories such as airplane and bird, or car and truck are also adjacent in the figures. One main difference between the DRPR and t-SNE representations for STL is the distance between the clusters ship and airplane. These categories are ac-
tually hard for the model to differentiate since they contain blue backgrounds and relatively similar objects. In particular, the STL airplane category contains many images of seaplanes lying on the sea and can then be mistaken for ships. This ambiguity between both categories is not observed on the t-SNE representation.
Due to lack of space, a detailed analysis for the CIFAR 100 and STL datasets is available in the supplementary material. A summary of the results is that categories that belong to the same superclass (e.g., categories hamster, mouse, rabbit, shrew, squirrel are part of the superclass small mammals) are grouped together with DRPR. The DRPR visualization also reflects some semantical structure: plants and insects are on the top left; animals are on the bottom left and categories on the right are outdoor categories. Medium mammals are also represented between small mammals and large carnivores.
In conclusion, the quantitative results show that the representations of DRPR are meaningful since they better preserve the cluster structure and allow observation of ambiguities between categories.

heading: 4.2 ZERO-SHOT LEARNING
text: We consider the same zero-shot learning scenario as Reed et al. (2016) and Snell et al. (2017). In particular, we test our approach on the same datasets and splits as them. The main goal is to learn two mappings, one for image representations and one for category representations, in a common space V . The latter mapping takes some category descriptions as input (e.g., from text description or visual attributes). Image representations are then learned so that they are closer to the representative of their category than to the representative of any other category. At test time, categories that were unseen during training are considered, and their representative is obtained by using the second mapping. An image is assigned to the category with closest representative.
Training datasets: We use the medium-scaled Caltech-UCSD Birds (CUB) dataset (Welinder et al., 2010) and Oxford Flowers-102 (Flowers) dataset (Nilsback & Zisserman, 2008). CUB contains 11,788 bird images from 200 different species categories split into disjoint sets: 100 categories for training, 50 for validation and 50 for test. Flowers contains 8,189 flower images from 102 different species categories: 62 categories are used for training, 20 for validation and 20 for test. To represent images, Reed et al. (2016) train a GoogLeNet (Szegedy et al., 2015) model whose output dimensionality is 1,024. For each category, Reed et al. (2016) extract some text annotations from
Table 3: Test accuracy on CUB
Method Accuracy
TMV-HLP Oquab et al. (2014) 47.9 % SJE Akata et al. (2015) 50.1% DS-SJE Reed et al. (2016) (Bag-of-words) 44.1 % DS-SJE Reed et al. (2016) (Char CNN-RNN) 54.0 % Ziming & Saligrama Zhang & Saligrama (2016) 55.3 % DS-SJE Reed et al. (2016) (Word CNN-RNN) 56.8 % Prototypical Networks Snell et al. (2017) 58.3 %
Ours – using DS-SJE (Char CNN-RNN) as supervision 57.7 % Ours – using Prototypical Networks as supervision 60.3 %
Table 4: Test accuracy on Flowers dataset
Method Accuracy
DS-SJE Reed et al. (2016) (Char CNN) 47.3 % DS-SJE Reed et al. (2016) (Bag-of-words) 57.7 % DS-SJE Reed et al. (2016) (Word CNN) 60.7 % DS-SJE Reed et al. (2016) (Char CNN-RNN) [reported] 63.7 % DS-SJE Reed et al. (2016) (Char CNN-RNN) [publicly available] 59.6 % DS-SJE Reed et al. (2016) (Word CNN-RNN) 65.6 % Prototypical Networks Snell et al. (2017) 63.9 %
Ours – using DS-SJE (Char CNN-RNN) [publicly available] 62.4 % Ours – using Prototypical Networks as supervision 68.2 %
which they learn a representative vector (e.g., based on Char CNN-RNN (Zhang et al., 2015)). The image representations of examples and the text representations of categories are learned jointly so that each image is more similar to the representative vector of its own category than to any other.
Supervision: We now describe how we generate the supervision/target of our model from the models pre-trained by (Reed et al., 2016; Snell et al., 2017) and provided by their respective authors.
Once its training is over, Prototypical Network (Snell et al., 2017) represents each image i by some vector f̃i and each category c by some vector µ̃c. By concatenating the different vectors into matrices F̃ = [f̃1, · · · , f̃n]> and M̃ = [µ̃1, · · · , µ̃k]> and formulating π̃ = 1k1k, DRPR considers the soft assignment matrix Y1 = ψ(F̃ , M̃ , π̃) ∈ Yn×k as target (i.e. supervision). In the case of Reed et al. (2016), we consider the same preprocessing as Snell et al. (2017). Each image i is represented by some vector f̌i, each category c is represented by some vector µ̌c (that is `2-normalized in Snell et al. 2017). The target soft assignment matrix of DRPR is then Y2 = ψ(F̌ , M̌ , 1k1k) in this case.
Trained models: In the model that Snell et al. (2017) provide and that obtains 58.3% accuracy on CUB, ProtoNet trains two models that take as arguments the representations learned by Reed et al. (2016). They train one model g̃θ̃1 for images such that ∀i, f̃i = g̃θ̃1(f̌i), and one model g̃θ̃2 for text representative vectors such that ∀c, µ̃c = g̃θ̃2(µ̌c). Following Snell et al. (2017), we train two (neural network) models: gθ1 for images, and gθ2 for categories. Both of them take as input the image and category representations used to create the target soft assignment matrix (i.e., we take the representations learned by Reed et al. (2016) when its probability scores Y2 are used as supervision, and the representations learned by Snell et al. (2017) otherwise). In this context, we alternately optimize gθ1 by fixing M (which depends on gθ2) and optimize gθ2 by fixing F (which depends on gθ1 ).
Implementation details: We consider that the learned models gθ1 and gθ2 have the same architecture and are multilayer perceptrons (MLP) with tanh activation functions. The number of hidden layers λ ∈ {0, 1, 2} and output dimensionality d are hyperparameters cross-validated from the accuracy on the validation set. More details on their architecture can be found in the supplementary material.
Results: We report the performance of our approach on the test categories of the CUB and Flowers datasets in Tables 3 and 4, respectively. The performance is measured as the average classification accuracy across all unseen classes. We use DS-SJE (Char CNN-RNN) and Prototypical Networks as supervision for our model because they are the only approaches whose pre-trained models are publicly available. Our approach obtains state-of-the-art results on both CUB and Flowers datasets by significantly improving the classification performance of the different classifiers. For instance, it improves the scores of 63.9% obtained by ProtoNet of Flowers up to 68.2%. In general, it improves zero-shot learning performance of the different classifiers by 2% to 4.3%.
Impact of dimensionality: We report in the supplementary material the performance of our model on both the validation and test sets using different numbers of hidden layers, and ranging the output dimensionality d from 16 to the dimensionality e of the input representations. Except for linear models (i.e. λ = 0), reducing the dimensionality improves generalization. This shows that the zeroshot learning performance of a given model can be significantly improved by taking its prediction scores as supervision of our model. To study the impact of the dimensionality reduction generated DRPR, we also ran the codes of Reed et al. (2016); Snell et al. (2017) by learning representations
with dimensionality smaller than those provided (using the same ranges as those in the tables of the supp. material). This decreased their generalization performance. Therefore, directly learning a low-dimensional representation is not a sufficient condition to generalize well. Our framework that learns representations so that examples with similar ambiguities (i.e. similar teacher predictions) are close to each other acts as a semantic regularizer. This is suggested by the fact that test accuracy is improved with DRPR even when e = d (as long as the MLPs contain hidden layers).
It is worth mentioning that one test category of the CUB dataset (Indigo bunting) belongs to the ImageNet dataset (Deng et al., 2009) that was used to pretrain GoogLeNet. By using the train/val/test category splits proposed by Xian et al., we did not observe a change of performance of the different models on CUB.

heading: 5 CONCLUSION
text: We have proposed a dimensionality reduction approach such that the soft clustering scores obtained in the low-dimensional space are similar to those given as input. We experimentally show that our approach improves generalization performance in zero-shot learning on challenging datasets. It can also be used to complement t-SNE, as a visualization tool to better understand learned models. In particular, we show that we can give a soft clustering interpretation to models that have probabilistic interpretations. Real-world applications that can be used with DRPR include distillation. For instance, when the teacher model is too large to store on a device with small memory (e.g., mobile phone), the student model which has a smaller memory footprint is used instead. Low-dimensional representations can also speed up retrieval tasks. Future work includes applying our approach to the task of distillation in the standard classification task where training categories are also test categories.

heading: 6 ACKNOWLEDGMENTS
text: We thank Fartash Faghri and the anonymous reviewers for their helpful comments on early versions of this manuscript. This work was supported by Samsung and the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (DoI/IBC) contract number D16PC00003. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. AMF acknowledges funding from the Canada CIFAR AI Chairs Program.
Disclaimer: The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoI/IBC, or the U.S. Government.

heading: A THEORETICAL ANALYSIS
text: This section provides a statistical guarantee for the algorithm presented in Section 2.2. We show that the minimizer of the empirical discrepancy (4), which uses a finite number of data points x1, · · · ,xn ∈ X , is close to the minimizer of the true expected discrepancy measure, to be defined shortly. Let us define the setup here.
We are given data points Dn = {X1, . . . , Xn}.3 We suppose that each Xi ∈ X is independent and identically distributed (i.i.d.) with the distribution ν ∈ M(X ), where M(X ) is the space of all probability distributions defined over X . The teacher is a fixed function φ = [φ(·; 1), . . . , φ(·; k)] that maps points in X to a k-dimensional simplex, and provides the target probability distributions. That is, the target yi for Xi is computed as yi = φ(Xi).
Consider a function space G whose domain is X and its range is a subset of Rd. This function space might be represented by a DNN, but we do not make such an assumption in our statistical analysis. Given a function g ∈ G (called gθ in the main article) and the number of clusters k, we define ψg(x) = [ψg(x; 1), . . . , ψg(x; k)] as
ψg(x; c) = πc exp
( −‖g(x)− µc(g)‖22 ) ∑k b=1 πb exp ( −‖g(x)− µb(g)‖22
) , with the cluster centres
µc(g) = E [φc(X)g(X)]
E [φc(X)] = ∫ φc(x)g(x)dν(x)∫ φc(x)dν(x) ,
and the priors
πc = ∫ φc(x)dν(x),
for c ∈ {1, . . . , k} (cf. Section 2.1). Note that ψg(x) defines a k-dimensional probability distribution, and µc(g) is a mapping of the function g to a point in Rd.
Similarly, given Dn, we define the empirical cluster centres
µ̂c(g) = 1 n
∑n i=1 φc(Xi)g(Xi)
1 n ∑n i=1 φc(Xi) ,
and their corresponding ψ̂g(x) = [ψ̂g(x; 1), . . . , ψ̂g(x; k)]
ψ̂g(x; c) = πc exp
( −‖g(x)− µ̂c(g)‖22 ) ∑k b=1 πb exp ( −‖g(x)− µ̂b(g)‖22
) . Note that here for simplicity of analysis, we assume that the priors πc are exact, and not estimated from data.
The student’s goal is to find a g such that ψg is close to φ. The closeness of the student to the teacher is defined based on their KL divergence. Specifically, Algorithm 1 minimizes the distorted empirical discrepancy (4), which can be written as4
∆n(ψ̂, φ) = 1
n n∑ i=1 KL(φ(Xi)||ψ̂(Xi)),
where Xis are from dataset Dn. Notice that the distorted empirical discrepancy ∆n(ψ̂, φ) is defined based on ψ̂, which uses the empirical centres µ̂c, instead of the expected centres µc. We also define an empirical discrepancy w.r.t. the true µc as ∆n(ψ, φ). This quantity is not accessible to the algorithm.
3Here we use Xi instead of xi in order to emphasize their randomness. 4Algorithm 1 does not explicitly specify what function generates y; we need to specify it here, so we use
φ(x) for that purpose.
We evaluate the quality of g, and its corresponding ψg , based on how well, in average, it performs on new points x ∈ X . We consider the expected KL divergence between ψ and φ w.r.t. distribution ν as the measure of performance. Therefore, the discrepancy is
∆(ψ, φ) = EX∼ν [KL(φ(X)||ψ(X))] = ∫ KL(φ(x)||ψ(x))dν(x).
The output of Algorithm 1 is the minimizer5 of ∆n(ψ̂, φ), which we denote by ĝ
ĝ ← argmin g∈G ∆n(ψ̂g, φ). (6)
We also define the minimizer of the discrepancy by g∗:
g∗ ← argmin g∈G ∆(ψg, φ).
We would like to compare the performance of ĝ when evaluated according to discrepancy, that is ∆(ψĝ, φ), and compare it with ∆(ψg∗ , φ).
Before stating our results, we enlist our assumptions. We shall remark on them as we introduce.
Assumption A1 (Samples) The dataset Dn = {Xi}ni=1 consists of i.i.d. samples drawn from ν(X ).
The i.i.d. assumption simplifies the analysis. With extra effort, one can provide similar results for some classes of dependent processes too. For example, if the dependent process comes from a time series and it gradually “forgets” its past, one may still obtain similar statistical guarantees. Forgetting can be formalized through the notion of “mixing” of the underlying stochastic process (Doukhan, 1994). One can then provide statistical guarantees for learning algorithms under various mixing conditions (Yu, 1994; Meir, 2000; Steinwart & Christmann, 2009; Mohri & Rostamizadeh, 2009; 2010; Farahmand & Szepesvári, 2012).
Assumption A2 (Bounded Function Space) Functions in G are L-bounded for some L > 0, i.e., g(x) ∈ [−L,L]d for any g ∈ G and x ∈ X .
This is a mild and realistic assumption on the function space G, and is mainly here to simplify some steps of the analysis.
Assumption A3 (Teacher) Part I) The output of the teacher φ is a probability distribution, i.e., for any x ∈ X and c = {1, . . . , k}, we have φ(x; c) ≥ 0 and ∑k c=1 φ(x; c) = 1. Part II) We assume that πc = E [φc(X)] is bounded away from zero for all c. We set πmin = minc πc.
This first part of the assumption explicitly expresses the fact that the algorithm expects to receive a probability distribution from the teacher. If it does not, for example if φ(x; c) is negative for some x and c but we treat it as a probability in the calculation of the KL divergence, the algorithm would not be well-defined.
The second part of this assumption requires that prior probability for each cluster is bounded away from zero, and has the probability at least πmin. This is a technical assumption used by the proof technique; it might be possible that one can relax this assumption.
We need to make some assumptions about the function space G and its complexity, i.e., capacity. We use covering number (and its logarithm, i.e., metric entropy) as the characterizer of the complexity. The covering number at resolution ε is the minimum number of balls with radius ε required to cover the spaceM according to a particular metric. We use N (ε,G, ‖·‖) to denote the covering number of G w.r.t. the norm ‖·‖, which we shall explicitly specify. As ε decreases, the covering number increases (or more accurately, the covering number is non-decreasing). For example, the covering number for a p-dimensional linear function approximator with constraint on the magnitude of its
5Since our focus is on the statistical analysis, we ignore the issue of the quality of the optimization procedure, and the fact that one may not find a global minimizer of ∆n(ψ̂g, φ) within G for a function space described by a DNN.
functions behaves like O( 1εp ). A similar result holds when the subgraphs of a function space has a VC-dimension p. Function spaces whose covering number grows faster are more complex, and estimating a function within them is more difficult. This leads to larger estimation error. On the other hand, those function spaces are often (but not always) show better function approximation properties too. For more detail and many examples, refer to Györfi et al. (2002); van de Geer (2000); Steinwart & Christmann (2008); Giné & Nickl (2015).
Let us define a norm for the function space G: ‖g‖∞,2 = sup
x∈X ‖g(x)‖2
This is a mixed-norm where we compute the `2-norm for each g(x) ∈ Rd, and then take the supremum norm over the resulting `2-norm. We use this norm to characterize the covering number of G.
Assumption A4 (Metric Entropy) There exists constants B > 0 and 0 < α < 1 such that for any ε, the following covering number (i.e., metric entropy) condition is satisfied:
logN ( ε,G, ‖·‖∞,2 ) ≤ ( B
ε
)2α .
The logarithm of the covering number of G is O( 1ε2α ). It grows much faster than the metric entropy of linear models, which is O(p log( 1ε )). This behaviour is suitable to capture the complexity of large function spaces such as the Sobolev space Wk(Rd) and many reproducing kernel Hilbert spaces (RKHS).6 Note that we use a mixed-norm to define the covering number. The use of supremum norm in the definition might be considered conservative. Using a more relaxed norm, for example based on the empirical Lp(PX1:n)-norm for some 1 ≤ p <∞, is an open technical question for future work. Finally let us define the pointwise loss function
l(x; g) = KL(φ(x)||ψg(x)) = k∑ c=1 φ(x; c) log φ(x; c) ψg(x; c) . (7)
Notice that ∆n(ψg, φ) = 1n ∑n i=1 l(Xi; g) and ∆(ψg, φ) = Eν [l(X; g)]. We define the following function space: L = {x 7→ l(x, g) : g ∈ G } . (8)
We also define the entropy integral (Dudley integral) J (ε) = ∫ ε
0
√ log 2N ( u,G, ‖·‖∞,2 ) du. (9)
We are now ready to state the main theoretical result of this section. Theorem 1. Suppose that Assumptions A1, A2, and A3 hold. Consider ĝ obtained by solving (6).
There exists a finite c1 > 0 such that for any δ > 0, with probability at least 1− δ, we have
∆(ψĝ, φ) ≤min g∈G
∆(ψg, φ) + 64 √ 2L √ dJ
( 16dL2+2 log2 k
16L √ d ) √ n + c1 ( dL2 + log(k) )√ log(1/δ) n +
c2L √ d
πmin
[ dJ (2L)√
n + L
√ d ln(k/δ)
n
] .
Furthermore, if the metric entropy satisfies Assumption A4, there exist constants c2, c4 > 0 and a function c3(α) > 0, which depends only on α, such that for any δ > 0, with probability at least 1− δ, we have
∆(ψĝ, φ) ≤min g∈G ∆(ψg, φ) + c2
[ 1 + c3(α) ( BL √ d
dL2 + log k
)α] ( dL2 + log(k) )√ log(1/δ) n +
c4 πmin
[ L2−αBαd3/2
1− α + L2d
]√ log(k/δ)
n .
6 For Wk(Rd) = Wk,2(Rd), the Sobolev space defined w.r.t. the L2-norm of the weak derivatives, we can set α = d
2k , see e.g., Lemma 20.6 of Györfi et al. (2002).
This theorem provides a finite sample error upper bound on the true discrepancy of ĝ, and relates it to the expressiveness and complexity of the function space G, the number of samples n, and some other properties. The term ming∈G ∆(ψg, φ) = ∆(ψg∗ , φ) is the function approximation error, and reflects the expressiveness of G. This is the minimum achievable error given the function space G. The other terms in the upper bound correspond to the estimation error caused by having a finite number of data points.
Let us focus on the second part of the theorem, which is under the particular choice of covering number according to Assumption A4. In that case, the estimation error shows n−1/2 dependence on the number of samples, and hence decreases as we have more training data. We observe that the upper bound increases as the range L of the function space G, the dimension d of the low-dimensional space, and the number of clusters k increases.
The effect of using the distorted empirical discrepancy ∆n(ψ̂g, φ) instead of ∆n(ψg, φ) shows itself in the last term, i.e., the term with the constant multiplier of c2L √ d
πmin in the first part of the theorem and
c4 πmin in the second part. Note the appearance of πmin in the denominator of the last term. This causes a large error if one of the clusters is improbable. Also observe that if we have k clusters, πmin ≤ 1k (equality under the uniform distribution over classes), so we have at least a linear dependence on k in the upper bound. This dependence might be due to our proof technique; it remains to be seen whether this can be improved.
Proof. To simplify the notation, we denote ∆n(g) = ∆n(ψg, φ), ∆̂n(g) = ∆n(ψ̂g, φ), and ∆(g) = ∆(ψg, φ). We want to relate ∆(ĝ) to ∆(g∗), the supremum of the empirical process ∆(g)−∆n(g) and the supremum of the distortion of the empirical loss ∆n(g) − ∆̂n(g). We have the following relations:
∆(ĝ) = ∆n(ĝ) + (∆(ĝ)−∆n(ĝ)) = ∆̂n(ĝ) + ( ∆n(ĝ)− ∆̂n(ĝ) ) + (∆(ĝ)−∆n(ĝ))
≤ ∆̂n(g∗) + ( ∆n(ĝ)− ∆̂n(ĝ) ) + (∆(ĝ)−∆n(ĝ))
= ∆n(g ∗) + ( ∆̂n(g ∗)−∆n(g∗) ) + ( ∆n(ĝ)− ∆̂n(ĝ) ) + (∆(ĝ)−∆n(ĝ))
= ∆(g∗) + (∆n(g ∗)−∆(g∗)) + ( ∆̂n(g ∗)−∆n(g∗) ) + ( ∆n(ĝ)− ∆̂n(ĝ) ) + (∆(ĝ)−∆n(ĝ))
≤ ∆(g∗) + 2 sup g∈G |∆(g)−∆n(g)|+ 2 sup g∈G ∣∣∣∆n(g)− ∆̂n(g)∣∣∣ . (10) The first inequality is because of the the optimizer property of ĝ, i.e., ∆̂n(ĝ) ≤ ∆̂n(g) for any g ∈ G, including g∗.
We need to upper bound the supremum of the empirical process, that is supg∈G |∆(g) − ∆n(g)|, and the supremum of the distortion caused by using ψ̂ in minimizing ∆n instead of ψ, that is supg∈G |∆n(g)− ∆̂n(g)|, cf (6).
Upper Bounding supg∈G |∆(g)−∆n(g)|. We use Lemma 7 in Appendix B.3 in order to upper bound the supremum of the empirical process, supg∈G |∆(g)−∆n(g)|, which is equivalent to supl∈L ∣∣ 1 n ∑n i=1 l(Xi)− E [l(X)] ∣∣. That lemma, which is originally Theorem 2.1 of Bartlett et al. (2005), relates the supremum of the empirical process to the Rademacher complexity of L, defined in the same appendix.
To apply the lemma, we first provide upper bound on l(x) and Var [l(X)]. By assumption, g(x) is within [−L,L]d. It is not difficult to see that µc(g) is also within [−L,L]d; for example, see the proof leading to (19). So
‖g(x)− µc(g)‖22 ≤ 4dL 2.
We evoke Proposition 4 with the choice of f(x; c) = ‖g(x)− µc‖22 and L′ = 4dL2 to obtain that
l(x; g) ≤ 8dL2 + log2 k. (11)
Since l(x; g) is bounded, we have Var [l(X; g)] ≤ E [ l(X; g)2 ] ≤ ( 8dL2 + log2 k )2 .
By the choice of β = 1, B = (8dL2 + log2 k), and r = (8dL 2 + log2 k) 2 in Lemma 7, we get that for any δ1 > 0,
sup l∈L ∣∣∣∣∣ 1n n∑ i=1 l(Xi)− E [l(X)] ∣∣∣∣∣ ≤4E [Rn(L)] + (8dL2 + log2 k) √ 2 log(2/δ1) n +
2(8dL2 + log2 k) 5
6
log(2/δ1)
n , (12)
with probability at least 1− δ1. It remains to provide an upper bound on E [Rn(L)]. This can be done by using Dudley’s integral to relate the Rademacher complexity of L to the covering number of L. Afterwards, we use Lemma 3 in Appendix A.1 to relate the covering number of L to the covering number of G. We have
E [Rn(L)] ≤ 4 √
2√ n E [∫ diam(L) 0 √ log 2N (u,L, L2(PX1:n)) du ]
≤ 4 √
2√ n E [∫ 16dL2+2 log2 k 0 √ log 2N (u,L, ‖·‖∞) du ]
≤ 4 √
2√ n E [∫ 16dL2+2 log2 k 0 √ log 2N ( u 16L √ d ,G, ‖·‖∞,2 ) du ]
= 4 √ 2(8L √ d)J
( 16dL2+2 log2 k
16L √ d ) √ n . (13)
In the second inequality, we benefit from two observations: first, we use l(x) ≤ 8dL2 + log2 k for any l ∈ L (11) to upper bound diam(L); second, the covering number w.r.t. L2(PX1:n) can be upper bounded by the covering number w.r.t. the supremum norm.7

heading: Upper Bounding supg∈G
text: ∣∣∣∆n(g)− ∆̂n(g)∣∣∣. We use Proposition 5 in Appendix A.2, which states
that for any δ2 > 0,
sup g∈G ∣∣∣∆n(g)− ∆̂n(g)∣∣∣ ≤ c1L√d πmin [ dJ (2L)√ n + L (√ d ln(k/δ2) n + ln(dk/δ2) n )]
≤ c2L √ d
πmin
[ dJ (2L)√
n + L
√ d ln(k/δ2)
n
] , (14)
with probability at least 1− δ2. Plugging (12) and (14) in (10) and using the entropy integral upper bound (13) lead to the desired result of the first part.
To prove the second part of the theorem, we use logN ( ε,G, ‖·‖∞,2 ) ≤ ( B ε )2α to calculate J (ε), which results in
J (ε) ≤ B αε1−α
1− α .
By plugging in ε = 16dL2 + 2 log2 k, we get that
E [Rn(L)] ≤ 32 √ 2(BL √ d)α
( 8dL2 + log2 k )1−α (1− α) √ n .
7This version of Dudley’s integral is from Theorem 2.3.7 of Giné & Nickl (2015). We use it similar to the argument in their proof of Theorem 3.5.1 and the comments thereafter. One could also use Theorem A.7 by Bartlett et al. (2005) (originally from Dudley) and note that the upper bound of the integral does not need to go up to infinity when the function space is bounded in the norm.
We upper bound J (2L) in (14) to obtain:
J (2L) ≤ 2 1−αBαL1−α
1− α .
After some simplifications these lead to the desired result of the second part.

heading: A.1 SOME TECHNICAL TOOLS
text: We develop some technical tools required in the proof of Theorem 1. Proposition 2 provides Lipschitz constant for some functions that are used in the proof of Lemma 3 to relate the covering number of the function space L, defined in (8), to that of G. This was a key step of the proof of the theorem. Proposition 4 provides an upper bound on the magnitude of l(x; f), shortly defined in (16).
We introduce a few more notations to reduce the clutter. Let dg(x; c) = ‖g(x)− µc(g)‖22, so we can write
ψg(x; c) = πc exp (−dg(x; c))∑k b=1 πb exp (−dg(x; b)) . c ∈ {1, . . . , k}
For a function f : X × {1, . . . , k} → R, we define
pf (x; c) = πc exp (−f(x; c))∑k b=1 πb exp (−f(x; b)) . c ∈ {1, . . . , k} (15)
We overload the pointwise loss function l(x; g), defined in (7), and define a similar definition for l(x; f) as follows
l(x; f) = k∑ c=1 φ(x; c) log φ(x; c) pf (x; c) . (16)
It is clear that with the choice of f = dg , the probability distribution pf is the same as ψg .
The following proposition specifies the Lipschitz properties of dg and l(x; f). Proposition 2. Suppose that Assumption A3 hold.
Part I) Consider g1, g2 ∈ G and let Assumption A2 hold. For any x ∈ X and c ∈ {1, . . . , k}, we have
|dg1(x; c)− dg2(x; c)| ≤ 4L √ d [ ‖g1(x)− g2(x)‖2 + sup
x′∈X ‖g1(x′)− g2(x′)‖2
] .
Part II) Consider two functions f1, f2 : X × {1, . . . , k} → R. For any x ∈ X we have |l(x; f1)− l(x; f2)| ≤ 2 ‖f1(x; ·)− f2(x; ·)‖∞ .
Proof. Part I) First notice that for any two vectors u and v, by the Cauchy-Schwarz inequality we have ∣∣∣‖u‖22 − ‖v‖22∣∣∣ = ∣∣∣∣∣∑ i (u2i − v2i ) ∣∣∣∣∣ = ∣∣∣∣∣∑ i (ui − vi)(ui + vi) ∣∣∣∣∣ ≤ ∣∣∣∣∣∣ √∑
i
(ui − vi)2 √∑
i
(ui + vi)2 ∣∣∣∣∣∣ = ‖u− v‖2 ‖u+ v‖2 . Consider g1, g2 ∈ G, and their corresponding dg1 and dg2 . We have
|dg1(x; c)− dg2(x; c)| = ∣∣∣‖g1(x)− µc(g1)‖22 ‖g2(x)− µc(g2)‖22∣∣∣ ≤‖(g1(x)− µc(g1))− (g2(x)− µc(g2))‖2× ‖(g1(x)− µc(g1)) + (g2(x)− µc(g2))‖2 ≤ [‖g1(x)− g2(x)‖2 + ‖µc(g1)− µc(g2)‖2]× ‖(g1(x) + g2(x))− (µc(g1) + µc(g2))‖2 . (17)
Note that ‖µc(g1)− µc(g2)‖2 = ∥∥∫ φc(x)(g1(x)− g2(x))dν(x)∥∥2∫ φc(x)dν(x) ≤ ∫ φc(x) ‖g1(x)− g2(x)‖2 dν(x)∫ φc(x)dν(x)
≤ supx ‖g1(x)− g2(x)‖2 ∫ φc(x)dν(x)∫
φc(x)dν(x) = sup x ‖g1(x)− g2(x)‖2 , (18)
where we used Jensen’s inequality and the fact that φc(x) ≥ 0. As µc(0) = 0, we also obtain that ‖µc(g)‖2 = ‖µc(g)− µc(0)‖2 ≤ sup
x ‖g(x)‖2 . (19)
As g(x) ∈ [−L,L]d, it holds that ‖g(x)‖2 ≤ L √ d. Therefore by (17), (18), and (19),
|dg1(x; c)− dg2(x; c)| ≤ 4L √ d [‖g1(x)− g2(x)‖2 + ‖µc(g1)− µc(g2)‖2]
≤ 4L √ d [ ‖g1(x)− g2(x)‖2 + sup
x′ ‖g1(x′)− g2(x′)‖2
] .
Part II) For functions f1, f2, using the definition of l(x; f) (16), we get that
l(x; f1)− l(x; f2) = ∑ c φ(x; c) log pf2(x; c) pf1(x; c) . (20)
By substituting the definition (15) and some simplifications, we get
log Pf2(x; c)
Pf1(x; c) = log  πc exp(−f2(x;c))∑kb=1 πb exp(−f2(x;b)) πc exp(−f1(x;c))∑k b=1 πb exp(−f1(x;b))  = [f1(x; c)− f2(x; c)] + [ρ(−f2(x; ·))− ρ(−f1(x; ·))] , (21)
with
ρ(u) = log (∑ b πb exp(ub) ) ,
where u ∈ Rk. We study the Lipschitz property of ρ(u) as a function of u in order to upper bound the second term on the right-hand side (RHS).
We take the derivative of ρ(u) w.r.t. the c-th component of u to obtain that
∂ρ(u)
∂uc = πc exp(uc)∑ b πb exp(ub) , qc(u).
Notice that qc(u) is a probability distribution. We denote (q1(u), . . . , qk(c)) by q(u).
By Taylor’s theorem, for u, u′ ∈ Rk, we have ρ(u′) = ρ(u) + 〈u′ − u , ∇uρ(ũ) 〉 = ρ(u) + 〈u′ − u , q(ũ) 〉 , for some ũ = (1 − λ)u + λu′ with 0 ≤ λ ≤ 1. By Hölder inequality, for any Hölder conjugate 1 r + 1 s = 1 (r, s ∈ [1,∞]), we get
|ρ(u′)− ρ(u)| ≤ ‖u′ − u‖r maxu≤ũ≤u′ ‖q(ũ)‖s ,
where max over u ≤ ũ ≤ u′ should be understood as the maximum over the line segment between u and u′.
In particular,
|ρ(u′)− ρ(u)| ≤ ‖u′ − u‖∞ maxu≤ũ≤u′ ‖q(ũ)‖1 = ‖u ′ − u‖∞ .
Here we used the fact that qc(u) is a probability distribution and its sum is equal to 1, for any choice of ũ.
We substitute (21) in (20) and use the upper bound |ρ(−f2(x; ·)) − ρ(−f1(x; ·))| ≤ ‖f1(x; ·)− f2(x; ·)‖∞, which is just shown, to get that
|l(x; f1)− l(x; f2)| ≤ 2 ‖f1(x; ·)− f2(x; ·)‖∞ , as desired.
The following lemma relates the covering number of the function space L to the covering number of the function space G. Lemma 3. Consider the function space G and its induced function space L (8). Let Assumptions A2 and A3 hold. For any ε > 0, we have
N (ε,L, ‖·‖∞) ≤ N ( ε
16L √ d ,G, ‖·‖∞,2
) .
Proof. By choosing f(x; c) = dg(x; c) = ‖g(x)− µc(g)‖22 and using both parts of Proposition 2, we get that
|l(x; g1)− l(x; g2)| ≤ 2 ‖dg1(x, ·)− dg2(x, ·)‖∞
≤ 8L √ d [ ‖g1(x)− g2(x)‖2 + sup
x′∈X ‖g1(x′)− g2(x′)‖2
] .
If we have an ε-cover of G w.r.t. ‖g‖∞,2 = supx∈X ‖g(x)‖2, it induces a 16L √ dε-cover on L w.r.t. the supremum norm.
The following proposition upper bounds the magnitude of l(x; f). Proposition 4. Suppose that Assumption A3 hold and |f(x; c)| ≤ L′ for any x ∈ X and c ∈ {1, . . . , k}. Consider pf defined in (15). It holds that
k∑ c=1 φ(x; c) log φ(x; c) pf (x; c) ≤ 2L′ + log2 k.
Proof. For simplicity, we ignore the dependence on x. We use the definition of pf (15) to get k∑ c=1 φ(c) log φ(c) pf (c) = ∑ c φ(c) log φ(c)− ∑ c φ(c) log πc + ∑ c φ(c)f(c) +
∑ c φ(c) log (∑ b πb exp(−f(b)) ) .
Let us consider each term on the RHS.
• As log φ(c) ≤ 0, we have ∑ c φ(c) log φ(c) ≤ 0.
• The priors (π1, . . . , πk) indeed defines a probability distribution, as each πc is nonnegative and ∑ c πc = ∑ c ∫ φc(x)dν(x) = ∫ ∑ c φc(x)dν(x) = ∫ 1 × dν(x) = 1. So
− ∑ c φ(c) log πc is the entropy of a probability distribution over an alphabet with size k, which is at most log2 k.
• The summation ∑ c φ(c)f(c) is upper bounded by L
′ because of the boundedness of f(c) ≤ L′ and the fact that φ is a probability distribution and sums to one.
• Consider the term ∑ c φ(c) log ( ∑ b πb exp(−f(b))). By the boundedness of f(b), we have∑
b πb exp(−f(b)) ≤ ∑ b πb exp(+L
′) ≤ exp(L′). Therefore, the term is upper bounded by ∑ c φ(c) log(exp(L ′)) ≤ L′.
Collecting all these terms leads to the upper bound of 2L′ + log2 k.

heading: A.2 CONTROLLING THE DISTORTED LOSS
text: This section provides an upper bound on the distortion of the empirical discrepancy, i.e., supg∈G |∆n(g)− ∆̂n(g)|. Proposition 5. Suppose that Assumptions A1, A2, and A3 hold. For any δ > 0, there exists a constant c1 > 0 such that with probability at least 1− δ, we have
sup g∈G ∣∣∣∆n(g)− ∆̂n(g)∣∣∣ ≤ c1L√d πmin [ dJ (2L)√ n + L (√ d ln(k/δ) n + ln(dk/δ) n )] .
Proof. Recall that dg(x; c) = ‖g(x)− µc(g)‖22, ∆n(g) = ∆n(ψg, φ), and ∆̂n(g) = ∆n(ψ̂g, φ).
Let us define d̂g(x; c) = ‖g(x)− µ̂c(g)‖22. Also from the definition of l(x, f) (16), we can write ∆n(ψg, φ) = 1 n ∑n i=1 l(Xi; dg). Similarly we have ∆̂n(g) = ∆n(ψ̂g, φ) = 1 n ∑n i=1 l(Xi; d̂g). Therefore, for any g ∈ G and by the application of Proposition 2 (Part II), we have∣∣∣∆n(g)− ∆̂n(g)∣∣∣ = ∣∣∣∣∣ 1n n∑ i=1 l(Xi, dg)− l(Xi, d̂g)
∣∣∣∣∣ ≤ 2 n n∑ i=1 ∥∥∥dg(Xi; ·)− d̂g(Xi; ·)∥∥∥ ∞ , (22)
where the supremum norm is taken over the centres c = {1, . . . , k}. For any x ∈ X and any c = {1, . . . , k}, we have∣∣∣dg(x; c)− d̂g(x; c)∣∣∣ = ∣∣∣‖g(x)− µc(g)‖22 − ‖g(x)− µ̂c(g)‖22∣∣∣
= ‖µc(g)− µ̂c(g)‖2 ‖2g(x)− (µc(g) + µ̂c(g))‖2 . (23)
It was shown in (19) that ‖µc(x)‖2 is upper bounded by supx∈X ‖g(x)‖2. One may similarly show the same for ‖µ̂c(x)‖2:
‖µ̂c(g)‖2 = ∥∥∥∥ 1n∑ni=1 φc(Xi)g(Xi)1
n ∑n i=1 φc(Xi) ∥∥∥∥ 2 ≤ 1 n ∑n i=1 φc(Xi) ‖g(Xi)‖2 1 n ∑n i=1 φc(Xi)
≤ supx∈X ‖g(x)‖2 1 n
∑n i=1 φc(Xi)
1 n ∑n i=1 φc(Xi)
= sup x∈X ‖g(x)‖2 .
As g(x) ∈ [−L,+L]d, we have that ‖g(x)‖2, ‖µc(x)‖2, and ‖µ̂c(x)‖2 are all L √ d-bounded, so∣∣∣dg(x; c)− d̂g(x; c)∣∣∣ ≤ 4L√d ‖µc(g)− µ̂c(g)‖2 .
This together with (22) and (23) show that∣∣∣∆n(g)− ∆̂n(g)∣∣∣ ≤ 8L√dmax c ‖µc(g)− µ̂c(g)‖2 . (24)
Proposition 6, which we prove soon, upper bounds supg∈G ‖µc(g)− µ̂c(g)‖2. By a union bound argument over c ∈ {1, . . . , k}, we get that for any fixed δ > 0, there exists a constant c1 > 0 such that
sup g∈G ∣∣∣∆n(g)− ∆̂n(g)∣∣∣ ≤ c1(8L√d) πmin [ dJ (2L)√ n + L (√ d ln(k/δ) n + ln(dk/δ) n )] ,
with probability at least 1− δ.
This proposition upper bounds the supremum of the `2 distance between cluster centres.
Proposition 6. Suppose that Assumptions A1, A2, and A3 hold. Consider a fixed c ∈ {1, . . . , k}. For any δ > 0, there exists a constant c1 > 0 such that with probability at least 1− δ, we have
sup g∈G ‖µc(g)− µ̂c(g)‖2 ≤ c1 πmin
[ dJ (2L)√
n + L
(√ d ln(1/δ)
n +
ln(d/δ)
n
)] .
Proof. To shorten the formulae, we use the notation En [f(Xi)] = 1n ∑n i=1 f(Xi) to denote the empirical expectation. We can decompose µc(g)− µ̂c(g) as follows µc(g)− µ̂c(g) = E [φc(X)g(X)]
E [φc(X)] − En [φc(Xi)g(Xi)] En [φc(Xi)]
= E [φc(X)g(X)]− En [φc(Xi)g(Xi)] + En [φc(Xi)g(Xi)] E [φc(X)] − En [φc(Xi)g(Xi)] En [φc(Xi)]
= En [φc(Xi)g(Xi)] [
1 E [φc(X)] − 1 En [φc(Xi)] ] ︸ ︷︷ ︸
(I)
+ E [φc(X)g(X)]− En [φc(Xi)g(Xi)]
E [φc(X)]︸ ︷︷ ︸ (I) .
In the rest of the proof, we provide an upper bound for Term (I) and Term (II).
Term (I): Fix δ1 > 0. We denote 1d as the d-dimensional vector with all components equal to 1. As g(x) ∈ [−L,+L]d and φc(x) ≥ 0, we have
|En [φc(Xi)g(Xi)]| ≤ 1
n d∑ i=1 |φc(Xi)||g(Xi)| ≤ L1d 1 n d∑ i=1 |φc(Xi)| = L1dEn [φc(Xi)] ,
where the comparison is dimension-wise. Therefore,∣∣∣∣En [φc(Xi)g(Xi)] (En [φc(Xi)]− E [φc(X)])En [φc(Xi)]E [φc(X)] ∣∣∣∣ ≤ L1dEn [φc(Xi)] |En [φc(Xi)]− E [φc(X)]|En [φc(Xi)]E [φc(X)]
= L1d
E [φc(X)] |En [φc(Xi)]− E [φc(X)]| . (25)
We provide a probabilistic upper bound on |En [φc(Xi)]− E [φc(X)]|. Since φc is a fixed 1-bounded function, we use Hoeffding’s inequality to upper bound it. After some manipulations, we obtain that
|En [φc(Xi)]− E [φc(X)]| ≤ √ ln(2/δ1)
2n ,
with probability at least 1− δ1. We use this along (25) and the assumption that E [φc(X)] ≥ πmin to get
‖(I)‖2 ≤ L √ d
πmin
√ ln(2/δ1)
2n , (26)
with probability at least 1− δ1. Term (II): We want to upper bound the supremum of the norm of second term (II). Fix δ2 > 0. To simplify the notation, let us first define function f(x) = φc(x)g(x) corresponding to a function g. Functions f are mapping from X to Rd. Furthermore, we define the function space F = {x 7→ φc(x)g(x) : g ∈ G }. For each dimension s ∈ {1, . . . , d}, we also define Fs = {x 7→ φc(x)gs(x) : g ∈ G }, where gs is the s-th dimension of g. With this notation, we can write
sup g∈G ∥∥∥∥E [φc(X)g(X)]− En [φc(Xi)g(Xi)]E [φc(X)] ∥∥∥∥
2
= 1
E [φc(X)] sup f∈F ‖En [f(Xi)]− E [f(X)]‖2 .
Let us focus on the supremum over F term. We have
sup f∈F ‖En [f(Xi)]− E [f(X)]‖22 = sup f∈F d∑ s=1 |En [fs(Xi)]− E [fs(X)]|2
≤ d∑ s=1 sup fs∈Fs |En [fs(Xi)]− E [fs(X)]|2 .
We take the square root of both sides and use the fact that √ a21 + . . .+ a 2 d ≤ a1 + . . . + ad (for non-negative ai) to get
sup f∈F ‖En [f(Xi)]− E [f(X)]‖2 ≤ d∑ s=1 sup fs∈Fs |En [fs(Xi)]− E [fs(X)]| .
Notice that as φc(x) is bounded by 1 and g is L-bounded dimension-wise, each fs(x) is also Lbounded. We use Lemma 7 (Appendix B.3) on each term to obtain a high probability upper bound. We choose the parameters of that lemma as β = 1, B = L, r = L2, and δ = δ2/d. This leads to
sup f∈F ‖En [f(Xi)]− E [f(X)]‖2 ≤ d∑ s=1
[ 4E [Rn(Fs)] + L √ 2 ln(2d/δ2)
n +
8L
3
ln(2d/δ2)
n
] , (27)
with probability at least 1− δ2. In order to control E [Rn(Fs)], we relate it to the covering number of Fs. We see that the covering number of Fs can be upper bounded by the covering number of F , which in turn can be upper bounded by the covering number of G. First, we use Dudley’s integral to get
E [Rn(Fs)] ≤ 4 √
2√ n E [∫ diam(Fs) 0 √ log 2N (u,Fs, L2(PX1:n)) du ] . (28)
The covering number argument is as follows. Choose any functions fs, f ′s ∈ Fs. For any sequence x1:n, the squared of the empirical `2-norm is
1
n n∑ i=1 |fs(xi)− f ′s(xi)| 2 ≤ 1 n n∑ i=1 d∑ r=1 |fr(xi)− f ′r(xi)| 2 = 1 n n∑ i=1 ‖f(xi)− f ′(xi)‖ 2 2
= 1
n n∑ i=1 ‖φc(xi)g(xi)− φc(xi)g′(xi)‖ 2 2 ≤ 1 n n∑ i=1 ‖g(xi)− g′(xi)‖ 2 2 ,
where in the last step we used the fact that φc(x) ≤ 1. An ε-cover of G w.r.t. ‖·‖∞,2 induces an ε-cover of G w.r.t. L2(Px1:n), which in turn induces an ε-cover on Fs w.r.t. L2(Px1:n), as the inequality above shows. Notice that this holds for any choice of x1:n, including X1:n appearing in the Dudley’s integral (28). Therefore, by (28) and setting diam(Fs) = 2L, we have
E [Rn(Fs)] ≤ 4 √
2√ n E [∫ 2L 0 √ log 2N ( u,G, ‖·‖∞,2 ) du ] = 4 √ 2J (2L)√ n , (29)
where we used the definition of the entropy integral of G (9). After plugging this upper bound in (27) and use the upper bound (26), which was obtained for Term (I), we see that
sup g∈G ‖µc(g)− µ̂c(g)‖2 ≤‖(I)‖2 +
1
E [φc(X)] sup f∈F ‖En [f(Xi)]− E [f(X)]‖2
≤L √ d
πmin
√ ln(2/δ1)
2n +
d
πmin
[ 16 √
2J (2L)√ n + L
√ 2 ln(2d/δ2)
n +
8L
3
ln(2d/δ2)
n
] ,
with probability at least 1− (δ1 + δ2). We set δ1 = δ2 = δ/2 and simplify the upper bound to obtain the desired result.

heading: B AUXILIARY RESULTS
text: For the convenience of the reader, we collect some auxiliary definitions and results that are used in the our proofs.

heading: B.1 FUNCTION SPACE AND NORMS
text: Here we briefly define some of the notations that we use throughout the paper. Consider the domain X and a function space F : X → R. We do not deal with measure theoretic considerations, so we just assume that all functions involved are measurable w.r.t. an appropriate σ-algebra for that space. We useM(X ) to refer to the set of all probability distributions defined over that space. We use symbols such as ν ∈M(X ) to refer to probability distributions defined over that space. We use ‖f‖p,ν to denote the Lp(ν)-norm (1 ≤ p <∞) of a measurable function f : X → R, i.e.,
‖f‖pp,ν , ∫ X |f(x)|pdν(x).
The supremum norm is defined as ‖f‖∞ = sup
x∈X |f(x)|
Let x1, . . . , xn be a sequence of points in X . We use x1:n to refer to this sequence. The empirical measure Pn is the probability measure that puts a mass of 1n at each xi, i.e.,
Pn = Px1:n = 1
n n∑ i=1 δxi ,
where δx is the Dirac’s delta function. For Dn = z1:x, the empirical L2(Pn)-norm of function f : X → R is
‖f‖2Dn = ‖f‖ 2 Pn
= 1
n n∑ i=1 |f(xi)|2 .
We can also define other Lp(Pn)-norms similarly.8 When there is no chance of confusion about Dn, we may denote the empirical norm simply by ‖f‖n.

heading: B.2 THE METRIC ENTROPY AND THE COVERING NUMBER
text: We quote the definition of the covering number from Györfi et al. (2002). Definition 1 (Definition 9.3 of Györfi et al. 2002). Let ε > 0, F be a set of real-valued functions defined on X , and νX be a probability measure on X . Every finite collection of Nε = {f1, . . . , fNε} defined on X with the property that for every f ∈ F , there is a function f ′ ∈ Nε such that ‖f − f ′‖p,νX < ε is called an ε-cover of F w.r.t. ‖·‖p,νX .
Let N (ε,F , ‖·‖p,νX ) be the size of the smallest ε-cover of F w.r.t. ‖·‖p,νX . If no finite ε-cover exists, take N (ε,F , ‖·‖p,νX ) =∞. Then N (ε,F , ‖·‖p,νX ) is called an ε-covering number of F and logN (ε,F , ‖·‖p,νX ) is called the metric entropy of F w.r.t. the same norm.
Given a x1:n = (x1, . . . , xn) ⊂ X and its corresponding empirical measure Pn = Px1:n , we can define the empirical covering number of F w.r.t. the empirical norm ‖·‖p,x1:n and is denoted by Np(ε,F , x1:n) = N (ε,F , Lp(x1:n)).

heading: B.3 RADEMACHER COMPLEXITY
text: We define Rademacher complexity and quote a result from Bartlett & Mendelson (2002). For more information about Rademacher complexity, we refer the reader to Bartlett & Mendelson (2002); Bartlett et al. (2005).
Let σ1, . . . , σn be independent random variables with P {σi = 1} = P {σi = −1} = 1/2. For a function space F : X → R, define RnF = supf∈F 1n ∑n i=1 σif(Xi) with Xi ∼ ν. The Rademacher complexity (or average) of F is E [RnG], in which the expectation is w.r.t. both σ and Xi. Rademacher complexity appears in the analysis of the supremum of an empirical process right after the application of the symmetrization technique. As such, its behaviour is closely related to the
8Or maybe more clearly, Pn(A) = 1n ∑n i=1 I{xi ∈ A} for any measurable subset A ⊂ X .
behaviour of the empirical process. One may interpret the Rademacher complexity as a complexity measure that quantifies the extent that a function from F can fit a noise sequence of length n (Bartlett & Mendelson, 2002).
The following is a simplified (and slightly reworded) version of Theorem 2.1 of Bartlett et al. (2005).
Lemma 7. Let F : X → R be a measurable function space with B-bounded functions. Let X1, . . . , Xn ∈ X be independent random variables. Assume that for some r > 0, Var [f(Xi)] ≤ r for every f ∈ F . Then for every δ > 0, with probability at least 1− δ,
sup f∈F ∣∣∣∣∣E [f(X)]− 1n n∑ i=1 f(Xi) ∣∣∣∣∣ ≤ inf β>0 { 2(1 + β)E [Rn(F)] + √ 2r ln(2/δ) n + 2B ( 1 3 + 1 β ) log(2/δ) n } .

heading: C EXPERIMENTAL RESULTS
text: We give implementation details about the experiments of our submitted paper in Section C.1. We report the detailed performance of our model as a function of the output dimensionality and the number of hidden layers in the zero-shot learning context in Section C.2. We show in Section C.3 that our method can be generalized to hard clustering by using implicit centers. We give additional visualization results in Section C.4.
C.1 IMPLEMENTATION DETAILS
We coded our method in PyTorch and ran all our experiments on a single Nvidia GeForce GTX 1060 which has 6GB of RAM.
PyTorch automatically calculates the gradient w.r.t. the mini-batch representation F . Nonetheless, it is worth mentioning that both the first and second arguments of our prediction function ψ(F,M,π) depend on F in the case where the centers are implicit (i.e., when we write M = diag(Y >1n)−1Y >F ). In this case, the gradient of our loss function w.r.t. F depends on both the first and second arguments of ψ.

heading: C.1.1 ZERO-SHOT LEARNING EXPERIMENTS
text: We now give details specific to the zero-shot experiments.
In the zero-shot learning experiment where F and M are computed from different sources (i.e., images and text) and are the output of two different networks, the optimization is performed by alternately optimizing one variable while fixing the other.
Mini-batch size: The training datasets of CUB and Flowers contain 5894 and 5878 images, respectively. In order to fit into memory, we set our mini-batch sizes as 421 (= 5894/14) and 735 (≈ 5878/8) for CUB and Flowers, respectively. Reed et al. (2016) and Snell et al. (2017) use 10 different views per image during training (middle, upper left, upper right, lower left and lower right crops for the original and horizontally-flipped image), and 1 view for test which is the middle crop of the original image. On the other hand, we use only 1 view per image (i.e., the middle crop of the original image) during training and test.
Optimizer: We use the Adam optimizer with a learning rate of 10−5 to train both models ϕθ1 and gθ2 .
Target soft assignment matrix: When we use the representations provided by Snell et al. (2017), our target soft assignment matrix is Y1 = ψ(F̃ , M̃ , π̃) = Ψ̃ ∈ Yn×k where the matrices F̃ and M̃ are provided and π̃ = 1k1k. The elements of Ψ̃ are written Ψ̃ic = π̃c exp(−d(f̃i,µ̃c))∑k e=1 π̃e exp(−d(f̃i,µ̃e))
. In this case, we formulate:
d(f̃i, µ̃c) = 1
10 ‖f̃i − µ̃c‖2
Using a temperature of 10 made the optimization more stable as it avoided gradients with high values.
We use a temperature of 2 when using the representations provided by Reed et al. (2016).
Initial temperature of our model: To make our optimization framework stable, we start with a temperature of 50. We then formulate our Bregman divergence as:
d(fi,µc) = 1
50 ‖fi − µc‖2
where fi and µc are the representations learned by our model. We decrease our temperature by 10% (i.e., tempt+1 = 0.9tempt) every 3000 epochs until the algorithm stops training.
We stop training at 10k epochs on CUB and 1k epochs on Flowers.
C.1.2 VISUALIZATION EXPERIMENTS
We now give details specific to the visualization experiments.
Dataset size: To be comparable to t-SNE, we directly learn two-dimensional embeddings instead of neural networks. Our mini-batch size is the size of the test set (i.e., the number of examples is n = 104 for most datasets except STL that contains n = 8000 test examples).
Optimizer: We use the RMSprop optimizer with a learning rate of 10−3, α = 0.99, = 10−6, the weight decay and momentum are both 0, and the data is not centered. We also we formulate the empirical discrepancy loss ∆n(Ψ, Y ) = ∑n i=1DKL(yi‖ψi) instead of ∆n(Ψ, Y ) = 1 n ∑n i=1DKL(yi‖ψi).
Target soft assignment matrix: Let us note zi = [zi,1, · · · , zi,k]> ∈ Rk the vector containing the logits of the learned representation of the i-th test example. We formulate yi = [yi,1, · · · , yi,k]> ∈ Rk our target assignment vector for the i-th test example as follows:
yi,c = exp( zi,c τ )∑k
e=1 exp( zi,e τ )
(30)
where τ = 5 for all the dataset except CIFAR-100 for which τ = 4.
We report the quantitative scores for τ = 1.
Initial temperature of our model: We learned our representation by using a fixed temperature of 1 (i.e., using the standard squared Euclidean distance).
We stop the algorithm after 8000 iterations.
Tuning t-SNE: we tested different ranges of scaling (1/1, 1/10, 1/100) and perplexity (i.e., 1, 10, 30 (default) and 100) and reported the representations that obtained the best quantitative results.
C.2 IMPACT OF DIMENSIONALITY IN ZERO-SHOT LEARNING
Let e ∈ N be the dimensionality of the representations taken as input and d the output dimensionality of the models gθ1 and gθ2 , the architecture of the models is e-d and e-e-d in the 1 and 2 hidden layer cases, respectively. The hyperparameter d ∈ {16, 32, 64, · · · , e} is also a hyperparameter cross-validated on the validation set.
We give the detailed accuracy performance of our model in Tables 5 to 9.
• Table 5 reports the test performance of our model on CUB when using the features provided by Reed et al. (2016) as supervision for different numbers of hidden layers and values of output dimensionality of our model.
• Table 6 (resp. Table 7) reports the validation (resp. test) performance of our model on CUB when using the features provided by Snell et al. (2017) as supervision.
• Table 8 (resp. Table 9) reports the validation (resp. test) performance of our model on Flowers when using the features provided by Snell et al. (2017) as supervision.
Dimensionality reduction improves performance, though the optimal dimensionality is dataset specific. In general, increasing the number of hidden layers also helps.

heading: C.3 GENERALIZATION TO HARD CLUSTERING
text: We validate that DRPR can be used to perform hard clustering as in Snell et al. (2017) but with implicit centers. To this end, we train a neural network with 2 convolutional layers on MNIST (LeCun et al., 1998) followed by a fully connected layer. Its output dimensionality is d = 2 or d = 3, the mini-batch size is n = 1000, the number of categories is k = 10 and the target hard assignment matrix Y ∈ {0, 1}n×k contains category membership information (i.e., Yic is 1 if the example fi belongs to category c, 0 otherwise). We train the model on the training set of MNIST and plot in Fig. 6 the representations of the test set. By assigning each test example to the category with closest centroid (obtained from the training set), the model obtains 98% (resp. 99%) accuracy when d = 2 (resp. d = 3). DRPR can then be learned for hard clustering when the centers are implicitly written as a function of the mini-batch matrix representation F and the target hard assignment matrix Y .
C.4 VISUALIZATION RESULTS
We now present visualization results.

heading: C.4.1 ARTIFACTS WITH T-SNE
text: Fig. 5 illustrates the CIFAR 100 representation learned by t-SNE when its input data is the target probability distribution that we give as supervision/input of our algorithm.
Following the recommendations mentioned in https://lvdmaaten.github.io/tsne/ when the representation form a strange ball with uniformly distributed points and obtaining very low error, we decreased the perplexity from 30 (which is the value by default) to 10 and 1 and divided our data by 10, 100 and 1000. Nonetheless, we still obtained the same type of representation as in Fig. 5.
This kind of artifact is the reason why we only report results obtained with logits.

heading: C.4.2 ADAPTING T-SNE WITH OTHER DIVERGENCES
text: We plot in Fig. 7 the visualization obtained by t-SNE when using the KL or JS divergences to compare pairs of probability distribution representations. The representations obtained in this case are still worse than using the original 3-dimensional representations as the cluster structures are not preserved, nor the inter-cluster distances. This suggests that comparing pairs of examples, as done by t-SNE, is less appropriate than our method that considers similarities between examples and the different k = 6 clusters.

heading: C.4.3 ADDITIONAL RESULTS
text: Fig. 8 illustrates the DRPR and t-SNE representations of CIFAR 10. Animal categories are illustrated on the right whereas machines are on the left.
Fig. 9 illustrates the DRPR and t-SNE representations of CIFAR 100. We learned the representations by exploiting 100 clusters but plot only 20 colors (one for each superclass of CIFAR 100), which is why multiple spikes have the same color. Groups with same colors are better defined with our approach than with t-SNE, this means that different categories from the same superclass (e.g., hamster, mouse, rabbit, shrew, squirrel which are small mammals) are grouped together with DRPR. One can observe a semantic structure in the 2D representation of DRPR: plants and insects are on the top left; animals are on the bottom left and categories on the right are outdoor categories. Medium mammals are also represented between small mammals and large carnivores.
Figures 10, 11, 12 and 13 illustrate the representations learned by our model for the STL, MNIST, CIFAR 100 and CIFAR 10 datasets, respectively. Instead of using colors that represent the categories of the embeddings as done in the submitted paper, we directly plot the images.
In general, we observe that images towards the end of spikes consist of a clearly visible object in a standard viewpoint on a simple background. Those closer to the center often have objects with a non-standard viewpoint or have a complex textured background. At a high-level, the classes appear to be organized by their backgrounds.
Taking the STL-10 visualization as an example, deer and horses are close together since they both tend to be found in the presence of green vegetation. These classes are far from boats and planes, which often have solid blue backgrounds. Looking more closely, the ordering of classes is sensible. Planes are neighbors to both boats (similar background) and birds (similar silhouette). And trucks neighbor both cars (similar background) and horses, which appear visually similar, particularly for images in which the horse is pulling a cart.
Taking the MNIST visualization as another example, one can observe that written characters in spikes are easy to recognize as they correspond to examples for which the learned model has high confidence in its scores. On the other hand, ambiguous examples are between multiple spikes (e.g., the characters 0 and 6 between spikes are more ambiguous than their neighbors in spikes).

Review summaries:
Authors propose a method of embedding training data examples into low-dimensional spaces such that mixture probabilities from a mixture model on these points are close to probability predictions from the original model in terms of KL divergence. Authors suggest two use-cases of such an approach: 1) data visualization, and 2) zero-shot learning. For the visualization use-case, authors compare against other dimensionality reduction methods with qualitative analysis on a synthetic problem, as well as evaluation metrics such as Neighborhood-Preservation Ratio and Clustering Distance Preservation Ratio. For zero-shot use-case, they take pre-trained models on two zero-shot tasks, and improve the accuracy by using probability outputs from pre-trained models as target.

This paper introduces a new supervised dimensionality reduction model. Supervision is provided in the form of class probabilities and the learning algorithm learns low-dimensional representations such that posterior cluster assignment probabilities given the representations match the observed class probabilities. The representations can be learned directly or the parameters of a neural network can be learned which maps inputs to the lower-dimensional space. The authors provide an extensive theoretical analysis of the proposed method and evaluate it on dimensionality reduction, visualization, and zero-shot learning tasks.