Title: DIMENSIONALITY REDUCTION FOR REPRESENTING THE KNOWLEDGE OF PROBABILISTIC MODELS

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.

1. INTRODUCTION:

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.

2. DIMENSIONALITY REDUCTION FOR PROBABILISTIC REPRESENTATIONS (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.

2.1. FAMILY OF EFFICIENT SOFT CLUSTERING ALGORITHMS:

Probability density: We consider that we are given a set of n vectors f 1 ,

• • • , f n ∈ V concatenated into a single matrix F = [f 1 , • • • , f n ] ∈ V n .

In the following, we consider V = R d , and

V n = R n×d .

The goal is to partition n examples into k soft clusters. Each cluster C c with c ∈ {1, • • • , k} has a center μc ∈ V and its corresponding probability density is p c (f i ) = exp(-d(f i , μc )) b(f i ), 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(f i ) = 1/(2π) d if d is the squared Euclidean distance). For simplicity, we consider that d is the squared Euclidean distance. The density p c (f i ) decreases as the divergence between the example f i 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 } k c=1 of a mixture model p(

f i |Γ) = k c=1 πc exp(-d(f i , μc )) b(f i )

where πc is the prior probability that f i is generated by C c . 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(f i |Γ) Assignment matrix: Partitioning the n observations in F into k soft clusters is equivalent to determining some soft assignment matrix Ŷ ∈ (0, 1) n×k in the set Y n×k of matrices whose rows are positive and sum to 1. Formally, Y n×k is written Y n×k := { Ŷ ∈ (0, 1) n×k : Ŷ 1 k = 1 n } where 1 k ∈ {1} k is the k-dimensional vector containing only 1. For a given value of Γ, the element Ŷic = p(C c |f i ) ∈ (0, 1) is the posterior probability, or responsibility of C c for f i . The higher the value of Ŷic , the more likely f i belongs to cluster C c .

Local maximum condition::

Once the BSCP has converged to a local maximum of max Γ n i=1 p(f i |Γ), the following equations are all satisfied:

(E-step) ∀i, c, p(C c |f i ) = Ŷic = πc exp(-d(f i , μc )) k m=1 πm exp(-d(f i , μm )) (1) (M-step) ∀c, μc = n i=1 Ŷic f i n i=1 Ŷic and ∀c, πc = 1 n n i=1

Ŷic:

(2)

Eq. ( 1) corresponds to the E-step of the EM algorithm that computes p(C c |f i ) = Ŷ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 p c (f i ), so b can be ignored. 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 y i ∈ [0, 1] k as target for each training example f i . These representations are concatenated into a single ma-

trix Y = [y 1 , • • • , y n ] ∈ Y n×k which

is the target of our method for F . DRPR learns the representation of F ∈ V n 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 matrix ∀i, c,

F = [f 1 , • • • , f n ] ∈ V n = R n×d , the centers M = [µ 1 , • • • , µ k ] ∈ V k = R

EQUATION

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:

EQUATION

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=1 D KL (y i ψ 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 ∈ R d .

Terminology::

In our experiments, the target (or teacher) is the assignment matrix Y ∈ Y n×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 Y 1 and Y 2 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 

(i.e., Y = [y1, • • • , yn] ∈ Y n×k ) 3: Create matrix F ← [f1, • • • , fn] ∈ V n such that ∀i, fi = g θ (xi) 4: Create matrix of centers M ← diag(Y 1n) -1 Y F and prior vector π ← 1

n Y 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 1 n ) -1 Y F and π = 1

n Y 1 n 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: π ← 1 k 1 k . 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(f i , μc ) = f i -μc 2 2 . 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 10 5 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. f i is:

EQUATION

One can observe that its magnitude depends on both the target scores Y ic ∈ (0, 1) and the predicted responsibilities Ψ ic . The gradient tries to make the vector f i 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.

3. RELATED WORK:

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 Y ic = 1 if the training example f i has to belong to C c , 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 C c given an example

f i is p(C c |f i ) = 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 π = 1 k 1 k , Y ∈ {0, 1} n×k is a hard assignment matrix and ∆ n is the same ∆ n as ours (i.e., ∆ n (Ψ, Y ) = 1 n n i=1 D KL (y i ψ 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. representation obtained by our method by exploiting using the same supervision as (c). The relative inter-cluster distances of the original dataset are preserved with our approach, unlike t-SNE.

4. EXPERIMENTS:

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.

4.1. VISUALIZATION:

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 = [x 1 , • • • , x n ] ∈ R n×3 (

Y ic = exp(-xi-cc 2 ) k e=1 exp(-xi-ce 2 )

where c c ∈ R 3 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 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 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 actually 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.

4.2. ZERO-SHOT LEARNING:

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) 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 fi and each category c by some vector μc . By concatenating the different vectors into matrices

F = [ f1 , • • • , fn ] and M = [ μ1 , • • • , μk ] and formulating π = 1 k 1 k , DRPR considers the soft assignment matrix Y 1 = ψ( F , M , π) ∈ Y n×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 fi , 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

Y 2 = ψ( F , M , 1 k 1 k ) 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, fi = gθ 1 ( fi ), 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 Y 2 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.

5. CONCLUSION:

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.

A THEORETICAL ANALYSIS:

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 x 1 , • • • , x n ∈ 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 D n = {X 1 , . . . , X n }. 3 We suppose that each X i ∈ 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 y i for X i is computed as y i = φ(X i ).

Consider a function space G whose domain is X and its range is a subset of R d . 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) 2 2 k b=1 π b exp -g(x) -µ b (g) 2 2

, 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 R d .

Similarly, given D n , we define the empirical cluster centres

μc (g) = 1 n n i=1 φ c (X i )g(X i ) 1 n n i=1 φ c (X i ) ,

and their corresponding ψg (x) = [ ψg (x; 1), . . . , ψg (x; k)]

ψg (x; c) = π c exp -g(x) -μc (g) 2 2 k b=1 π b exp -g(x) -μb (g) 2 2

.

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(φ(X i )|| ψ(X i )),

where X i s are from dataset D n . 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.

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

∆(ψ, φ) = E X∼ν [KL(φ(X)||ψ(X))] = KL(φ(x)||ψ(x))dν(x).

The output of Algorithm 1 is the minimizer5 of ∆ n ( ψ, φ), which we denote by

EQUATION

We also define the minimizer of the discrepancy by g * :

g * ← argmin g∈G ∆(ψ g , φ).

We would like to compare the performance of ĝ when evaluated according to discrepancy, that is ∆(ψ ĝ , φ), 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 D n = {X i } n i=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 = min c π 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 space M 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 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) ∈ R d , 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:

log N ε, 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 W k (R d ) 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 L p (P X1: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 , φ) = 1 n n i=1 l(X i ; g) and ∆(ψ g , φ) = E ν [l(X; g)].

We define the following function space:

L = { x → l(x, g) : g ∈ G } .

(8) We also define the entropy integral (Dudley integral)

EQUATION

)

We are now ready to state the main theoretical result of this section. Theorem 1. Suppose that Assumptions A1, A2, and A3 hold. Consider ĝ obtained by solving (6).

There exists a finite c 1 > 0 such that for any δ > 0, with probability at least 1 -δ, we have

∆(ψ ĝ , φ) ≤ min g∈G ∆(ψ g , φ) + 64 √ 2L √ d J 16dL 2 +2 log 2 k 16L √ d √ n + c 1 dL 2 + log(k) log(1/δ) n + c 2 L √ d π min dJ (2L) √ n + L d ln(k/δ) n .

Furthermore, if the metric entropy satisfies Assumption A4, there exist constants c 2 , c 4 > 0 and a function c 3 (α) > 0, which depends only on α, such that for any δ > 0, with probability at least 1 -δ, we have

∆(ψ ĝ , φ) ≤ min g∈G ∆(ψ g , φ) + c 2 1 + c 3 (α) BL √ d dL 2 + log k α dL 2 + log(k) log(1/δ) n + c 4 π min L 2-α B α d 3/2 1 -α + L 2 d log(k/δ) n .

This theorem provides a finite sample error upper bound on the true discrepancy of ĝ, and relates it to the expressiveness and complexity of the function space G, the number of samples n, and some other properties. The term min g∈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 ≤ 1 k (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 ∆(ĝ) 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:

EQUATION

The first inequality is because of the the optimizer property of ĝ, i.e., ∆n (ĝ) ≤ ∆n (g) for any g ∈ G, including g * .

We need to upper bound the supremum of the empirical process, that is

sup g∈G |∆(g) -∆ n (g)|,

and the supremum of the distortion caused by using ψ in minimizing ∆ n instead of ψ, that is

sup g∈G |∆ n (g) -∆n (g)|, cf (6).

Upper Bounding sup g∈G |∆(g) -∆ n (g)|. We use Lemma 7 in Appendix B.3 in order to upper bound the supremum of the empirical process,

sup g∈G |∆(g) -∆ n (g)|, which is equivalent to sup l∈L 1 n n i=1 l(X i ) -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) 2 2 ≤ 4dL 2 .

We evoke Proposition 4 with the choice of f (x; c) = g(x) -µ c 2 2 and L = 4dL 2 to obtain that l(x; g) ≤ 8dL 2 + log 2 k.

(11) Since l(x; g) is bounded, we have

Var [l(X; g)] ≤ E l(X; g) 2 ≤ 8dL 2 + log 2 k 2 .

By the choice of β = 1, B = (8dL 2 + log 2 k), and r = (8dL 2 + log 2 k) 2 in Lemma 7, we get that for any δ 1 > 0,

EQUATION

with probability at least 1 -δ 1 .

It remains to provide an upper bound on E [R n (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

EQUATION

In the second inequality, we benefit from two observations: first, we use l(x) ≤ 8dL 2 + log 2 k for any l ∈ L (11) to upper bound diam(L); second, the covering number w.r.t. L 2 (P X1:n ) can be upper bounded by the covering number w.r.t. the supremum norm.7 

Upper Bounding sup g∈G ∆ n (g) -∆n (g) . We use Proposition 5 in Appendix A.2, which states that for any δ 2 > 0,

EQUATION

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 log

N ε, G, • ∞,2 ≤ B ε 2α to calculate J (ε),

which results in

J (ε) ≤ B α ε 1-α 1 -α .

By plugging in ε = 16dL 2 + 2 log 2 k, we get that

E [R n (L)] ≤ 32 √ 2(BL √ d) α 8dL 2 + log 2 k 1-α (1 -α) √ n .

We upper bound J (2L) in ( 14) to obtain:

J (2L) ≤ 2 1-α B α L 1-α 1 -α .

After some simplifications these lead to the desired result of the second part.

A.1 SOME TECHNICAL TOOLS

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

d g (x; c) = g(x) -µ c (g) 2 2

, so we can write

ψ g (x; c) = π c exp (-d g (x; c)) k b=1 π b exp (-d g (x; b)) . c ∈ {1, . . . , k}

For a function f : X × {1, . . . , k} → R, we define

EQUATION

We overload the pointwise loss function l(x; g), defined in ( 7), and define a similar definition for l(x; f ) as follows

EQUATION

It is clear that with the choice of f = d g , the probability distribution p f is the same as ψ g .

The following proposition specifies the Lipschitz properties of d g and l(x; f ). Proposition 2. Suppose that Assumption A3 hold.

Part I) Consider g 1 , g 2 ∈ G and let Assumption A2 hold. For any x ∈ X and c ∈ {1, . . . , k}, we have

|d g1 (x; c) -d g2 (x; c)| ≤ 4L √ d g 1 (x) -g 2 (x) 2 + sup x ∈X g 1 (x ) -g 2 (x ) 2 .

Part II) Consider two functions f 1 , f 2 : X × {1, . . . , k} → R. For any x ∈ X we have

|l(x; f 1 ) -l(x; f 2 )| ≤ 2 f 1 (x; •) -f 2 (x; •) ∞ .

Proof. Part I) First notice that for any two vectors u and v, by the Cauchy-Schwarz inequality we have

u 2 2 -v 2 2 = i (u 2 i -v 2 i ) = i (u i -v i )(u i + v i ) ≤ i (u i -v i ) 2 i (u i + v i ) 2 = u -v 2 u + v 2 .

Consider g 1 , g 2 ∈ G, and their corresponding d g1 and d g2 . We have

EQUATION

)

Note that

EQUATION

where we used Jensen's inequality and the fact that φ c (x) ≥ 0. As µ c (0) = 0, we also obtain that

EQUATION

) As g(x) ∈ [-L, L] d , it holds that g(x) 2 ≤ L √ d.

Therefore by ( 17), ( 18), and ( 19),

|d g1 (x; c) -d g2 (x; c)| ≤ 4L √ d [ g 1 (x) -g 2 (x) 2 + µ c (g 1 ) -µ c (g 2 ) 2 ] ≤ 4L √ d g 1 (x) -g 2 (x) 2 + sup x g 1 (x ) -g 2 (x ) 2 .

Part II) For functions f 1 , f 2 , using the definition of l(x; f ) (16), we get that

EQUATION

)

By substituting the definition (15) and some simplifications, we get

EQUATION

with

ρ(u) = log b π b exp(u b ) ,

where u ∈ R k .

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) ∂u c = π c exp(u c ) b π b exp(u b ) q c (u).

Notice that q c (u) is a probability distribution. We denote (q 1 (u), . . . , q k (c)) by q(u).

By Taylor's theorem, for u, u ∈ R k , we have

ρ(u ) = ρ(u) + u -u , ∇ u ρ(ũ) = ρ(u) + u -u , q(ũ) ,

for some ũ = (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 max u≤ũ≤u q(ũ) s ,

where max over u ≤ ũ ≤ u should be understood as the maximum over the line segment between u and u .

In particular,

|ρ(u ) -ρ(u)| ≤ u -u ∞ max u≤ũ≤u q(ũ) 1 = u -u ∞ .

Here we used the fact that q c (u) is a probability distribution and its sum is equal to 1, for any choice of ũ.

We substitute ( 21) in ( 20) and use the upper bound |ρ(-

f 2 (x; •)) -ρ(-f 1 (x; •))| ≤ f 1 (x; •) -f 2 (x; •) ∞ , which is just shown, to get that |l(x; f 1 ) -l(x; f 2 )| ≤ 2 f 1 (x; •) -f 2 (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) = d g (x; c) = g(x) -µ c (g) 2 2

and using both parts of Proposition 2, we get that

|l(x; g 1 ) -l(x; g 2 )| ≤ 2 d g1 (x, •) -d g2 (x, •) ∞ ≤ 8L √ d g 1 (x) -g 2 (x) 2 + sup x ∈X g 1 (x ) -g 2 (x ) 2 .

If we have an ε-cover of G w.r.t. g ∞,2 = sup x∈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 ). Proof. For simplicity, we ignore the dependence on x. We use the definition of p f (15) to get

k c=1 φ(c) log φ(c) p f (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 log 2 k. Collecting all these terms leads to the upper bound of 2L + log 2 k.

A.2 CONTROLLING THE DISTORTED LOSS:

This section provides an upper bound on the distortion of the empirical discrepancy, i.e., sup g∈G |∆ n (g) -∆n (g)|.

Proposition 5. Suppose that Assumptions A1, A2, and A3 hold. For any δ > 0, there exists a constant c 1 > 0 such that with probability at least 1 -δ, we have φ) , and ∆n (g) = ∆ n ( ψg , φ). Let us define dg (x; c) = g(x) -μc (g) 2 2 . Also from the definition of l(x, f ) (16), we can write

sup g∈G ∆ n (g) -∆n (g) ≤ c 1 L √ d π min dJ (2L) √ n + L d ln(k/δ) n + ln(dk/δ) n . Proof. Recall that d g (x; c) = g(x) -µ c (g) 2 2 , ∆ n (g) = ∆ n (ψ g ,

∆ n (ψ g , φ) = 1 n n i=1 l(X i ; d g ). Similarly we have ∆n (g) = ∆ n ( ψg , φ) = 1 n n i=1 l(X i ; dg ).

Therefore, for any g ∈ G and by the application of Proposition 2 (Part II), we have

EQUATION

where the supremum norm is taken over the centres c = {1, . . . , k}. For any x ∈ X and any c = {1, . . . , k}, we have

EQUATION

)

It was shown in (19) that µ c (x) 2 is upper bounded by sup x∈X g(x) 2 . One may similarly show the same for μc (x) 2 :

μc (g) 2 = 1 n n i=1 φ c (X i )g(X i ) 1 n n i=1 φ c (X i ) 2 ≤ 1 n n i=1 φ c (X i ) g(X i ) 2 1 n n i=1 φ c (X i ) ≤ sup x∈X g(x) 2 1 n n i=1 φ c (X i ) 1 n n i=1 φ c (X i ) = 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

d g (x; c) -dg (x; c) ≤ 4L √ d µ c (g) -μc (g) 2 .

This together with ( 22) and ( 23) show that

EQUATION

)

Proposition 6, which we prove soon, upper bounds sup g∈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 c 1 > 0 such that

sup g∈G ∆ n (g) -∆n (g) ≤ c 1 (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 c 1 > 0 such that with probability at least 1 -δ, we have

sup g∈G µ c (g) -μc (g) 2 ≤ c 1 π min dJ (2L) √ n + L d ln(1/δ) n + ln(d/δ) n .

Proof. To shorten the formulae, we use the notation E n [f (X i )] = 1 n n i=1 f (X i ) 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)] - E n [φ c (X i )g(X i )] E n [φ c (X i )] = E [φ c (X)g(X)] -E n [φ c (X i )g(X i )] + E n [φ c (X i )g(X i )] E [φ c (X)] - E n [φ c (X i )g(X i )] E n [φ c (X i )] = E n [φ c (X i )g(X i )] 1 E [φ c (X)] - 1 E n [φ c (X i )] (I) + E [φ c (X)g(X)] -E n [φ c (X i )g(X i )] 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 1 d as the d-dimensional vector with all components equal to 1. As g(x) ∈ [-L, +L] d and φ c (x) ≥ 0, we have

|E n [φ c (X i )g(X i )]| ≤ 1 n d i=1 |φ c (X i )||g(X i )| ≤ L1 d 1 n d i=1 |φ c (X i )| = L1 d E n [φ c (X i )] ,

where the comparison is dimension-wise.

Therefore,

EQUATION

)

We provide a probabilistic upper bound on

|E n [φ c (X i )] -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

|E n [φ c (X i )] -E [φ c (X)]| ≤ ln(2/δ 1 ) 2n ,

with probability at least 1 -δ 1 . We use this along (25) and the assumption that E

EQUATION

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 R d . Furthermore, we define the function space

F = { x → φ c (x)g(x) : g ∈ G }. For each dimension s ∈ {1, . . . , d}, we also define F s = { x → φ c (x)g s (x) : g ∈ G }

, where g s is the s-th dimension of g. With this notation, we can write

sup g∈G E [φ c (X)g(X)] -E n [φ c (X i )g(X i )] E [φ c (X)] 2 = 1 E [φ c (X)] sup f ∈F E n [f (X i )] -E [f (X)] 2 .

Let us focus on the supremum over F term. We have

sup f ∈F E n [f (X i )] -E [f (X)] 2 2 = sup f ∈F d s=1 |E n [f s (X i )] -E [f s (X)]| 2 ≤ d s=1 sup fs∈Fs |E n [f s (X i )] -E [f s (X)]| 2 .

We take the square root of both sides and use the fact that a

2 1 + . . . + a 2 d ≤ a 1 + . . . + a d (for non-negative a i ) to get sup f ∈F E n [f (X i )] -E [f (X)] 2 ≤ d s=1 sup fs∈Fs |E n [f s (X i )] -E [f s (X)]| .

Notice that as φ c (x) is bounded by 1 and g is L-bounded dimension-wise, each f s (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 = L 2 , and δ = δ 2 /d. This leads to

EQUATION

with probability at least 1 -δ 2 .

In order to control E [R n (F s )], we relate it to the covering number of F s . We see that the covering number of F s 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

EQUATION

The covering number argument is as follows. Choose any functions f s , f s ∈ F s . For any sequence x 1:n , the squared of the empirical 2 -norm is

1 n n i=1 |f s (x i ) -f s (x i )| 2 ≤ 1 n n i=1 d r=1 |f r (x i ) -f r (x i )| 2 = 1 n n i=1 f (x i ) -f (x i ) 2 2 = 1 n n i=1 φ c (x i )g(x i ) -φ c (x i )g (x i ) 2 2 ≤ 1 n n i=1 g(x i ) -g (x i ) 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. L 2 (P x1:n ), which in turn induces an ε-cover on F s w.r.t. L 2 (P x1:n ), as the inequality above shows. Notice that this holds for any choice of x 1:n , including X 1:n appearing in the Dudley's integral (28). Therefore, by (28) and setting diam(F s ) = 2L, we have

EQUATION

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 E n [f (X i )] -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.

B AUXILIARY RESULTS:

For the convenience of the reader, we collect some auxiliary definitions and results that are used in the our proofs.

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 X 1 , . . . , X n ∈ X be independent random variables. Assume that for some r > 0, Var [f (X i )] ≤ r for every f ∈ F. Then for every δ > 0, with probability at least 1 -δ,

sup f ∈F E [f (X)] - 1 n n i=1 f (X i ) ≤ inf β>0 2(1 + β)E [R n (F)] + 2r ln(2/δ) n + 2B 1 3 + 1 β log(2/δ) n .

C EXPERIMENTAL RESULTS:

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 1 n ) -1 Y F ).

In this case, the gradient of our loss function w.r.t. F depends on both the first and second arguments of ψ.

C.1.1 ZERO-SHOT LEARNING EXPERIMENTS:

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. 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. 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). 

:

We experimentally show in the supplementary material that DRPR can be learned for hard supervised clustering when the centers are implicit (see Section C.3).

-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.

We use the pre-trained models available at https://github.com/aaron-xichen/ pytorch-playground.

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.

Here we use Xi instead of xi in order to emphasize their randomness.

Algorithm 1 does not explicitly specify what function generates y; we need to specify it here, so we use φ(x) for that purpose.

Since 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.

For W k (R d ) = W k,2 (R d ),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 ofGyörfi et al. (2002).

This version of Dudley's integral is from Theorem 2.3.7 ofGiné & 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 byBartlett 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.

Figures and Tables:

Figure 1: DRPR learns low-dimensional representations that reflect the uncertainties of a pre-trained classifier (here for categories blue and yellow represented by their centers)

k×d and the priors π = [π 1 , • • • , π k ] ∈ (0, 1) k . As in Eq. (1), we define our prediction function ψ(F, M, π) = Ψ ∈ Y n×k as the soft assignment matrix predicted by the BSCP given F , M and π. The elements of the matrix Ψ are then computed as follows:

Figure 2: (a) Original 3-dimensional dataset containing 6 clusters (one color per cluster); (b) 2D representation obtained with t-SNE by exploiting the original 3D representation as input; (c) 2D representation obtained with t-SNE by exploiting soft probability scores w.r.t. the 6 clusters; (d) 2Drepresentation obtained by our method by exploiting using the same supervision as (c). The relative inter-cluster distances of the original dataset are preserved with our approach, unlike t-SNE.

where n is the number of examples) plotted in Fig. 2 (a). The target soft assignment matrix Y ∈ Y n×k where k = 6 is constructed by computing

Figure 4: MNIST (left half) and STL (right half) representations of DRPR (left) and t-SNE (right)

Figure 3: Level sets representing responsibilities of the green cluster whose center is the green circle. Grey circles are the centers of other clusters. The responsibility of the green cluster increases for the examples that are located in the spike at the right side of the figure.

Suppose that Assumption A3 hold and |f (x; c)| ≤ L for any x ∈ X and c ∈ {1, . . . , k}. Consider p f defined in (15). It holds that k c=1 φ(x; c) log φ(x; c) p f (x; c) ≤ 2L + log 2 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 .

Reed et al. (2016) andSnell 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 bySnell et al. (2017), our target soft assignment matrix isY 1 = ψ( F , M , π) = Ψ ∈ Y n×kwhere the matrices F and M are provided and π = 1 k 1 k . The elements of Ψ are written Ψic = πc exp(-d( fi, μc )) k e=1 πe exp(-d( fi, μe )) . In this case, we formulate:d( fi , μc ) = 1 10 fiμc

Figure 7: (left) Visualization obtained with t-SNE on the toy dataset when replacing the 2 distance in the original space by the KL divergence to compare probability distribution representations. (right) Visualization obtained with t-SNE on the toy dataset when replacing the 2 distance in the original space by the Jensen Shannon divergence.

Figure 8: CIFAR 10 representation of DRPR and t-SNE

Figure 9: CIFAR 100 representation of DRPR and t-SNE

Figure 11: MNIST representation learned by DRPR.

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 ∈ Y n×k containing the target probability scores y1, • • • , yn

Table content: <table/>

NPR performance (in %) for different values of κ

Table content: <table><tr><td/><td/><td>κ</td><td>1</td><td>2</td><td>5</td><td>10</td><td>50</td><td>100</td><td/><td>κ</td><td>1</td><td>2</td><td>5</td><td>10</td><td>50</td><td>100</td></tr><tr><td/><td/><td>logits</td><td>6.3</td><td>8.7</td><td>12.4</td><td>16.0</td><td>28.2</td><td>36.5</td><td/><td>logits</td><td>1.7</td><td>2.8</td><td>4.8</td><td>6.9</td><td>17.0</td><td>24.7</td></tr><tr><td>MNIST</td><td>dim reduction</td><td>PCA LLE ISOMAP t-SNE DRPR (Ours)</td><td>0 0 0.5 5.2 7.0</td><td>0.1 0.1 1.0 7.4 9.4</td><td>1.3 0.2 2.0 9.8 13.4</td><td>2.4 1.8 3.3 12.1 16.3</td><td>8.5 4.1 10.4 20.4 24.5</td><td>14.1 8.6 16.7 27.7 30.3</td><td>CIFAR10</td><td>PCA LLE ISOMAP t-SNE DRPR</td><td>0.1 1.3 0.2 1.9 8.4</td><td>0.2 1.6 0.3 2.7 10.2</td><td>0.3 3.4 0.6 4.1 13.2</td><td>0.8 4.2 1.0 5.6 15.2</td><td>3.2 9.8 4.1 12.8 21.1</td><td>5.9 15.6 7.5 19.5 26.6</td></tr><tr><td/><td/><td>logits</td><td>5.9</td><td>8.4</td><td>11.8</td><td>14.8</td><td>25.5</td><td>32.4</td><td/><td>logits</td><td>2.5</td><td>4.1</td><td>6.3</td><td>8.5</td><td>15.9</td><td>19.3</td></tr><tr><td>STL</td><td>dim reduction</td><td>PCA LLE ISOMAP t-SNE DRPR (Ours)</td><td>0.1 0.4 0.1 5.9 13.9</td><td>0.2 0.8 0.4 8.0 18.9</td><td>0.7 1.6 1.0 10.6 25.9</td><td>1.3 2.6 1.7 12.1 31.6</td><td>4.9 6.9 6.1 19.6 45.3</td><td>8.4 11.3 10.0 26.7 52.0</td><td>CIFAR100</td><td>PCA LLE ISOMAP t-SNE DRPR</td><td>0 0 0 2.4 9.1</td><td>0.1 0 0.1 3.7 12.6</td><td>0.2 0 0.4 5.0 17.9</td><td>0.4 0.1 0.6 6.5 24.0</td><td>1.6 0.2 2.3 13.4 52.7</td><td>2.8 1.4 3.9 16.1 58.6</td></tr><tr><td colspan="9">Quantitative evaluation metrics: We quantify the</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="9">relevance of our approach with the two following</td><td/><td colspan="6">Table 2: CDPR performance (in %)</td><td/></tr><tr><td colspan="9">evaluation metrics: (1) an adaptation of the Neigh-</td><td/><td>model</td><td/><td>MNIST</td><td>STL</td><td>CIFAR10</td><td colspan="2">CIFAR100</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>logits</td><td/><td>81.5</td><td>80.6</td><td>79.9</td><td colspan="2">67.9</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>dim reduction</td><td>PCA LLE ISOMAP t-SNE DRPR (ours)</td><td/><td>69.1 53.9 67.9 52.8 76.7</td><td>66.9 55.3 70.1 65.9 70.5</td><td>67.4 56.5 69.5 66.5 70.0</td><td colspan="2">56.4 53.8 58.7 60.4 69.0</td></tr></table>

Test accuracy on CUB

Table content: <table><tr><td>Method</td><td>Accuracy</td></tr><tr><td>TMV-HLP Oquab et al. (2014)</td><td>47.9 %</td></tr><tr><td>SJE Akata et al. (2015)</td><td>50.1%</td></tr><tr><td>DS-SJE Reed et al. (2016) (Bag-of-words)</td><td>44.1 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Char CNN-RNN)</td><td>54.0 %</td></tr><tr><td>Ziming &amp; Saligrama Zhang &amp; Saligrama (2016)</td><td>55.3 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Word CNN-RNN)</td><td>56.8 %</td></tr><tr><td>Prototypical Networks Snell et al. (2017)</td><td>58.3 %</td></tr><tr><td>Ours -using DS-SJE (Char CNN-RNN) as supervision</td><td>57.7 %</td></tr><tr><td>Ours -using Prototypical Networks as supervision</td><td>60.3 %</td></tr></table>

Test accuracy on Flowers dataset

Table content: <table><tr><td>Method</td><td>Accuracy</td></tr><tr><td>DS-SJE Reed et al. (2016) (Char CNN)</td><td>47.3 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Bag-of-words)</td><td>57.7 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Word CNN)</td><td>60.7 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Char CNN-RNN) [reported]</td><td>63.7 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Char CNN-RNN) [publicly available]</td><td>59.6 %</td></tr><tr><td>DS-SJE Reed et al. (2016) (Word CNN-RNN)</td><td>65.6 %</td></tr><tr><td>Prototypical Networks Snell et al. (2017)</td><td>63.9 %</td></tr><tr><td>Ours -using DS-SJE (Char CNN-RNN) [publicly available]</td><td>62.4 %</td></tr><tr><td>Ours -using Prototypical Networks as supervision</td><td>68.2 %</td></tr></table>

Peer review (iclr2019):
 The paper deals with a problem formulation adjacent to that of the sufficient dimension reduction: given training set of pairs (x_i,y_i), how to reduce the dimension of the first element, i.e. map x_i --> f(x_i), so that f(x_i)*s still have all the information to recover y_i*s. In the paper, the output y_i is a probability distribution over k labels that softly describes inclusion of example i into k classes. They consider a nonlinear case, i.e. the mapping f is taken from a prespecified set of mappings, parameterized by Theta (e.g. neural network). Then by "recovering y_i" they mean that EM algorithm for {f(x_i)} will result in a clustering of the data into k soft clusters similar to given {y_i}. The algorithm that is presented is quite natural, though no guarantees that it will converge to something relevant were given. Theoretical analysis deals with a question --- how far the empirical discrepancy could be from the true expected one. Especially, easiness of substitution of .bar{Y}_{ic} with Y_{ic} in the algorithm is unclear (roughly speaking, the latter means that E-step is omitted in EM). If matrix Y in algorithm is fixed, why we need to compute .pi in the loop? Isn*t it going to be the same? Does this algorithm really minimizes the discrepancy?