Melissa Ailem
email: melissa.ailem@inria.fr

Bowen Zhang

Aurelien Bellet
email: aurelien.bellet@inria.fr

Pascal Denis
email: pascal.denis@inria.fr

Fei Sha
email: feisha@usc.edu

A Probabilistic Model for Joint Learning of Word Embeddings from Texts and Images

Several recent studies have shown the benefits of combining language and perception to infer word embeddings. These multimodal approaches either simply combine pre-trained textual and visual representations (e.g. features extracted from convolutional neural networks), or use the latter to bias the learning of textual word embeddings. In this work, we propose a novel probabilistic model to formalize how linguistic and perceptual inputs can work in concert to explain the observed word-context pairs in a text corpus. Our approach learns textual and visual representations jointly: latent visual factors couple together a skip-gram model for co-occurrence in linguistic data and a generative latent variable model for visual data. Extensive experimental studies validate the proposed model. Concretely, on the tasks of assessing pairwise word similarity and image/caption retrieval, our approach attains equally competitive or stronger results when compared to other state-of-the-art multimodal models.

Introduction

Continuous-valued vector representation of words has been one of the key components in neural architectures for natural language processing [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF][START_REF] Pennington | Glove: Global vectors for word representation[END_REF][START_REF] Levy | Neural word embedding as implicit matrix factorization[END_REF]. The main idea is based on the distributional hypothesis [START_REF] Zellig | Distributional structure[END_REF], which states that words used in similar contexts have similar semantic meanings. To this end, words are mapped to points in an Euclidean space such that the displacement between their coordinates (i.e., embeddings) reflects similarity and difference in semantics [START_REF] Pennington | Glove: Global vectors for word representation[END_REF]. As such, word embeddings * On leave from U of Southern California have been shown to be useful in determining semantic and syntactic similarity between individual words [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF][START_REF] Baroni | Don't count, predict! a systematic comparison of context-counting vs. context-predicting semantic vectors[END_REF][START_REF] Levy | Improving distributional similarity with lessons learned from word embeddings[END_REF], as well as in downstream NLP tasks, e.g., sentiment analysis, question answering, and coreference resolution, just to name a few.

Most existing approaches rely solely on text corpora to infer word representations. While successful, the embeddings produced by such models do not necessarily reflect all inherent aspects of human semantic knowledge, such as the perceptual aspect [START_REF] Feng | Visual information in semantic representation[END_REF]. This has motivated many researchers to explore different ways to infuse visual information, often represented in the form of pre-computed visual features, into word embeddings [START_REF] Kiela | Learning Image Embeddings using Convolutional Neural Networks for Improved Multi-Modal Semantics[END_REF][START_REF] Silberer | Visually grounded meaning representations[END_REF][START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF][START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF]. The main theme is to take either the text embeddings, or the visual features or both as such to derive multimodal embeddings: through concatenation [START_REF] Kiela | Learning Image Embeddings using Convolutional Neural Networks for Improved Multi-Modal Semantics[END_REF], or by treating visual features as regression targets [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF][START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF].

Despite the success of these prior efforts in yielding multimodal embeddings and applying them to downstream NLP tasks, there are still several deficiencies. In particular, the visual features (as such) are not guaranteed to be suitable for the word embedding task since they are typically optimized independently for another objective (e.g., image classification). Hence, fusing pre-computed word representations and visual features may not be a good strategy.

To address the above issues, we explore a new way to integrate linguistic and perceptual information. We develop a new model which jointly learns word embeddings from text and extracts latent visual information, from pre-computed visual features, that could supplement the linguistic embeddings in modeling the co-occurrence of words and their contexts in a corpus. Instead of using pre-trained visual features as it is or as regression targets, we posit that they contain latent perceptual information that could complement text in representing words.

More specifically, the proposed model consists of two components. The visual component is an unsupervised probabilistic model for learning latent factors that generates the visual data. The linguistic component is a revised SKIP-GRAM model in which the text embeddings work in concert with the latent visual factors to explain the occurrence of wordcontext pairs in a corpus. One advantage of our joint modeling is that it allows two-way interaction. On one hand, the linguistic information can guide the extraction of latent visual factors. On the other hand, the extracted visual factors can improve the modeling of word-context co-occurrences in text data. Another appealing property of our model is its natural ability to propagate perceptual information to the embeddings of words lacking visual features (e.g., abstract words) during learning.

We conduct extensive quantitative and qualitative experiments to examine and understand the effectiveness of our approach, on the tasks of word similarity and image/caption retrieval. We show its matching or stronger performance when compared to other state-of-the-art approaches for learning multimodal embeddings.

Our Approach

We start by introducing the problem setup and notations. We then describe our model, namely PIXIE (ProbabIlistic teXtual Image Embeddings), for joint learning of word representations from text and perceptual information.

Setup and Background

We are given a corpus of H tokens (words): w 1 , . . . , w i , . . . , w H . From the corpus, we form a collection of word-context pairs δ w,c = (w, c), such that w ∈ V w , c ∈ V c , with V w and V c denoting respectively the word and context vocabularies. As in most previous work, the contexts for word w i are the words that surround it in a L-sized window. We introduce the binary indicator variables y wc , such that y wc = 1 if δ w,c appears in our collection, and y wc = 0 otherwise.

For some words with visual grounding (we will refer to them as visual words), we have access to a visual representation x w . In practice, we use con- volutional net features (see Section 4 for details).

SKIP-GRAM WITH NEGATIVE SAMPLING

(SGNS) The SGNS's objective is to learn word representations that are good at distinguishing the observed pairs (y wc = 1) from non-observed or "negative" pairs (y wc = 0), using logistic regression. Formally, SGNS maximizes the following log-likelihood:

w,c [y wc log σ(v T c e w )+(1-y wc ) log σ(-v T c e w )], (1)
where σ(•) is the sigmoid function, v c , e w denote respectively the vectors for the context c and target word w. The second term in (1) is intractable due to the large number of possible negative pairs, and is approximated by sampling N negative examples {c i } N i=1 for every observed pair of words and their contexts. This gives rise to the following objective function for each observed pair:

log σ(v T c e w ) + N i=1 log σ(-v T c i e w ), (2)
where c i is a (negative) context that does not appear in the context of w [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. In practice, criterion (2) is optimized in an online fashion, by using Stochastic Gradient Descent (SGD) over the observed pairs δ wc in the corpus. Each observed pair δ wc typically occurs several times in the corpus, therefore performing SGD over the corpus amounts to weighting equation (2) by the number of occurrence of each pair.

Joint Visual and Text Modeling

We now describe our model, namely PIXIE (Proba-bIlistic teXtual Image Embeddings) illustrated in Fig. 1, for joint learning of word representations from textual and perceptual information.

Formally, PIXIE is a probabilistic model of image features x and word-context pairs' labels y. Similar to SKIP-GRAM, PIXIE represents each word w and context c with low dimensional embeddings noted respectively e w ∈ R K and v c ∈ R K . PIXIE further assumes latent visual factors, z w ∈ R K , for each word's visual representation x w . Next we describe the two main components of PIXIE, namely textual and perceptual, in more details.

Perceptual Component Each visual vector x w is drawn conditional on its latent representation z w , i.e., x w ∼ p θ (x|z w ), with p θ (z) = N (0, I). Since x w is real valued, we let p θ (x w |z w ) be a Gaussian parameterized by a generative neural network (or decoder). That is,

p θ (x w |z w ) = N (x w |µ θ (z w ), Σ θ (z w )). (3)
For tractability purposes, Σ θ is restricted to be diagonal. Moreover, both the co-variance Σ θ (z w ) (its diagonal) and the mean µ θ are the outputs of a decoder network with parameters θ and input z w .

Textual Component To model the occurrence/absence of word-context pair δ wc in the linguistic corpus, we adopt a Bernoulli (Ber) model:

p(y wc | e w , v c , z w ) = Ber(σ[f (e w , v c , z w )]) (4)
The function f (•) defines how multimodal embeddings are fused. While many choices can be experimented, we use the simple additive model:

f (e w , v c , z w ) = v T c (e w + z w ). (5)
For words without visual representation, we simply set the corresponding latent factors z w to the zero vector. Note that, without the visual factors z w , equation (4) reduces to the Skip-Gram with negative sampling objective (1).

Joint Model

The perceptual and the textual information interact through the shared latent z w . The joint model of the above two sources of information takes the following form:

p(x w , y wc | e w , v c ) = p θ (z w )p θ (x w |z w )p(y wc | e w , v c , z w )dz w (6)
The intuition behind our joint formulation is to let the textual information guide the extraction of latent visual factors z w . Through equations ( 4) and ( 5), the model will put high probability on factors z w reflecting patterns that can supplement the linguistic embeddings e w in explaining the wordcontext co-occurrences. Thus, the extracted latent visual factors can contribute to improve the performance on predicting the occurrence of a word and its contexts in the linguistic corpus, which would encourage the model to leverage the perceptual information. The underlying assumption here is to infer visual and textual embeddings that can work in concert to represent words.

Visual Information Propagation Equation ( 5) implies that the embeddings z, e and v will affect each other during the learning process. Interestingly, if a non-visual word w 1 shares a similar context c with a visual word w 2 , then the factor z w 2 will affect e w 1 via v c . In other words, our formulation makes it possible to implicitly propagate perceptual information from one word to another through shared contexts. We illustrate this aspect in our experiments.

Approximate Inference and Learning

Training PIXIE amounts to inferring the posterior over the visual latent factors, p θ (z|x, y), as well as finding the decoder's parameters θ, the word and context embeddings, e and v, that maximize the likelihood (6). However, as in many complex probabilistic models, the likelihood (due to the integral over z) and the posterior are intractable.

We therefore resort to approximation techniques.

More precisely, we rely on Variational Inference (VI) [START_REF] David M Blei | Variational inference: A review for statisticians[END_REF]. The idea of VI is to introduce a tractable approximate posterior distribution q φ (z|x) (the variational distribution) and optimize a lower bound on the likelihood, known as Evidence Lower BOund (ELBO). The latter can be written for each word w as follows:

L w = E q [log p θ (x w |z w ) + c log p(y wc |e w , v c , z w )] -KL(q φ (z w |x w ) p θ (z w )) (7)
where KL(• •) is the Kullback-Leibler divergence.

The variational distribution is chosen to be a multivariate Gaussian parameterized by an inference network (or encoder) which takes x as input, namely

q φ (z w |x w ) = N (z w |µ φ (x w ), Σ φ (x w )), (8)
where we drop the dependency on all y wc variables to be computationally tractable. The pair of encoder and decoder neural networks gives rise to the interpretation of PIXIE's visual component (formed by z and x) as a probabilistic autoencoder. In fact, if we drop the textual part in PIXIE, namely y, e and v, then we recover the Variational Auto-Encoder (VAE) [START_REF] Diederik | Autoencoding variational bayes[END_REF]. Lastly, we approximate the intractable (i) expectation with respect to q φ (z|x) and the (ii) sum over the negative pairs in (7), by relying on a Monte Carlo estimator of L. Concretely, for (ii) we use negative sampling as in (2). Concerning (i), for every observed x w , we sample {z (j) w } J j=1 from q φ (z w |x w ) using the reparameterization trick (Kingma and Welling, 2013), i.e., z

(j) w = µ φ (x w ) + Σ φ (x w ) , with ∼ N (0, I).
Then we approximate L with:

L w = 1 J j log p θ (x w |z (j) w ) -KL(q φ (z w |x w ) p θ (z w )) + 1 J c,j y wc log σ[v T c (e w + z (j) w )] + 1 J c,j N i=1 y wc log(σ[-v T c i (e w + z (j) w )]). (9)
The last two summands correspond to the familiar conditional likelihood term in the SGNS model, augmented with latent visual factors. We optimize the objective (9) via SGD with respect to both the encoder/decoder networks parameters (θ and φ) and the embeddings (e and v). We evaluate the gradients of L, with respect to θ and φ using backpropagation. Similarly, the gradient with respect to e and v can be easily carried out using automatic differentiation tools. Our learning procedure is summarized in Algorithm 1.

Algorithm 1 Variational PIXIE Input: x, y, sample sizes B and J Steps: Randomly initialize θ, φ, e and v repeat

• Draw a minibatch W B of words: w (1) , . . . , w (B) • For each xw with w ∈ W B , sample {z

(j)
w } J j=1 from q φ (zw|xw) using the reparameterization trick. For each observed pair δwc draw

N negative examples {c i } N i=1 . • Compute the estimator L W B ← w∈W B Lw • Compute the gradient: G ← ∇ θ,φ,e,v L W B
• Use G to update θ, φ, e and v (e.g., with ADAM) until convergence return θ, φ, e and v Inference Once the parameters of the model are learned, for any given word with or without visual representation, we can compute its multimodal embedding. As a short hand, let the binary variable m w denote whether or not the word w has a visual representation. The multimodal embedding for w can be written as

s w = e w + m w µ(x w ), (10)
where µ(x w ) = E q (z w ) is the output of the encoding neural network, cf. Eq. ( 8).

In our experiments, we have also studied an alternative way to compose multimodal embeddings by concatenating the two vectors e w and µ(x w )

t w = [e w m w µ(x w )]. (11)
Note that for non-visual words, only zeros are appended to e w . One advantage of t over s is that if one uses distances to measure similarity, t can be seen as a simple summation of distances in two different spaces (in terms of e and µ respectively).

Related Work

Combining language and perception has been recently considered in various NLP tasks such as machine translation [START_REF] Calixto | Incorporating global visual features into attention-based neural machine translation[END_REF], visual question generation [START_REF] Mostafazadeh | Generating natural questions about an image[END_REF], image captioning [START_REF] Klein | Associating neural word embeddings with deep image representations using fisher vectors[END_REF], etc. In this work, we focus on learning word embeddings from images and texts. Multimodal embeddings have been studied in several recent research work. One strategy is to obtain word embeddings from linguistic data and visual data independently and then proceed with some kind of fusion steps. [START_REF] Kiela | Learning Image Embeddings using Convolutional Neural Networks for Improved Multi-Modal Semantics[END_REF] simply concatenates pre-trained linguistic word embeddings and visual features computed by convolutional nets. [START_REF] Bruni | Multimodal distributional semantics[END_REF] performs an additional step of dimensionality reduction via singular value decomposition. [START_REF] Silberer | Visually grounded meaning representations[END_REF] extend on this work by feeding the linguistic embedding and visual features into a stacked auto-encoder for nonlinear dimensionality reduction. The abovementioned approaches perform a two-stage process to derive multimodal representations (unimodal inference followed by fusion) and have been evaluated only on words for which both perceptual and textual representations are available.

A standing question is how to propagate visual information from words with visual features to words lacking them (for instance, abstract words). While the previous methods fall short on that, the recent work by [START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF] addresses this challenge by learning a mapping from language to vision, using a set of words with known linguistic embeddings and visual features. This mapping can then be used to infer visual representations for new words from their textual embeddings.

All the aforementioned methods rely on independently pre-trained linguistic embeddings and visual features. In this work, we propose a different strategy, which consists in adapting those representations so that the information can be fused in earlier stages. In this respect, the closest work to ours is [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF], which proposes to augment the SKIP-GRAM objective function with a term mapping the textual embeddings to the visual features. Crudely, the linguistic embeddings must therefore predict both the text co-occurrences and (pre-trained) visual features. We emphasize two key differences with our approach. First, instead of performing a regression or mapping from the textual embeddings to the visual features, our model learns to infer perceptual latent factors to retain only the portion of visual information that can supplement the linguistic embeddings in representing words. Second, while [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF] combines two objectives, we use a joint probabilistic model integrating both visual and text information in a principled way. Specifically, our model seeks latent factors that are good at explaining the word-context co-occurrences. For instance, a visual feature of (an image of) OCEAN often contains information about SKY and BLUE -such visual information could be beneficial to predict cooccurrence of tokens in the context of OCEAN. This desiderata further strengthens the learned embeddings to be visually grounded. In our experiments, we show that our approach tends to group concrete visually similar concepts together.

Experiments

In this section, we evaluate our model and contrast it to other competing approaches on two different tasks: word similarity and image/caption retrieval.

Setup

Text corpus We use the Text8 WIKIPEDIA corpus 1 containing over 17 million tokens. Text8 was pre-processed to contain only letters and nonconsecutive spaces. After removing infrequent words, we obtain a vocabulary of 50, 000 unique words.

Image features

We use the ImageNet dataset [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], including the fall Ima-1 http://mattmahoney.net/dc/textdata geNet 2011 release [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF]. It contains 14, 188, 125 images organized according to 21, 842 synsets of WordNet [START_REF] Fellbaum | WordNet[END_REF]. Each synset contains 600 images on average. To extract image features, we rely on the Caffe toolkit [START_REF] Jia | Caffe: Convolutional architecture for fast feature embedding[END_REF] and use the GoogLeNet convolutional neural nets [START_REF] Szegedy | Going deeper with convolutions[END_REF] pre-trained on the 1000 synsets of ILSVRC 2012. The 1024-dimensional activation of the pooling units (before the softmax layer) are then taken as our image features.

Visual representation of words For each word in the vocabulary, we recover all the synsets that it belongs to using the WordNet interface of the NLTK module (Python) [START_REF] Bird | Natural language processing with Python: analyzing text with the natural language toolkit[END_REF]. We then remove the synsets not covered by our Im-ageNet dataset. This results in 9,713 words, out of the 50,000 words in the vocabulary. For each visual word, we randomly draw 1,000 distinct images in ImageNet. If the number of images for a word is less than 1,000, we increase the coverage using images belonging to the hypernyms of the considered word's synsets, as in [START_REF] Kiela | Learning Image Embeddings using Convolutional Neural Networks for Improved Multi-Modal Semantics[END_REF]. We then take the average of these features as the word's visual representation x.

Hyper-parameter setting For all models, we set the dimension of linguistic and visual embeddings, e and z, to 100, following many previous works. In our model, the encoder/decoder neural networks are implemented as one-hidden-layer neural nets with 500 hidden units each. The dimensions of the inputs and the outputs of the decoder neural networks are 100 and 1, 024 respectively (1, 024 and 100 for the encoder). For the encoder, the hidden units are hyperbolic tangent, and the output units are linear. For the decoder, the hidden units are hyperbolic while the outputs are sigmoid. For SGNS, we set the window size L to 10 and the number of negative samples to 64. Our model is learned by Stochastic Gradient Descent using the ADAM optimizer (Kingma and Ba, 2014) with a learning rate set to 0.001.

Datasets

#word pairs MEN [START_REF] Bruni | Multimodal distributional semantics[END_REF] 3000 EN-MC (Miller and Charles, 1991) 31 EN-RG [START_REF] Rubenstein | Contextual correlates of synonymy[END_REF] 65 SimLex [START_REF] Hill | Simlex-999: Evaluating semantic models with (genuine) similarity estimation[END_REF] 999 MTurk [START_REF] Radinsky | A word at a time: computing word relatedness using temporal semantic analysis[END_REF][START_REF] Halawi | Large-scale learning of word relatedness with constraints[END_REF] 287 WORDSIM [START_REF] Finkelstein | Placing search in context: The concept revisited[END_REF] 350 REL [START_REF] Agirre | A study on similarity and relatedness using distributional and wordnet-based approaches[END_REF] 150 SIM [START_REF] Agirre | A study on similarity and relatedness using distributional and wordnet-based approaches[END_REF] 200 SEMSIM [START_REF] Silberer | Learning Grounded Meaning Representations with Autoencoders[END_REF] 5494 VISSIM [START_REF] Silberer | Learning Grounded Meaning Representations with Autoencoders[END_REF] 5494 [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF], § : [START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF] 4.2 Task 1: Word Similarity Datasets Word similarity is a common type of evaluation task for measuring the effectiveness of word embeddings. To this end, we retain 10 benchmark datasets consisting of pairs of words associated with similarity scores given by human judges. Table 1 summarizes their basic properties. There are different types of similarities being assessed. SEMSIM, SimLex, SIM, EN-RG and EN-MC focus on semantic or taxonomic similarity -e.g. CAR is similar to AUTOMOBILE. MEN, REL and MTurk consider general relatedness -e.g. CAR is related to GARAGE. VISSIM is about visual similarity -e.g. GOOSE looks like SWAN. Note that SIM and REL are the similarity and relatedness subsets of the full WORDSIM dataset [START_REF] Finkelstein | Placing search in context: The concept revisited[END_REF] respectively. VISSIM contains the same word pairs as SEMSIM.

Competing models We benchmark our model PIXIE against several strong uni-and multi-modal models listed below:

• SGNS: Skip-Gram with Negative Sampling [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. Without the visual component, PIXIE reduces to SGNS. We can thus assess the impact of the perceptual information by comparing PIXIE to SGNS. • VAE: Vatiational Auto-Encoder [START_REF] Diederik | Autoencoding variational bayes[END_REF], which corresponds to the visualspecific component of PIXIE.

• CNN: Visual features extracted from a convolutional neural net as described in Section 4.1. • CNN⊕SGNS [START_REF] Kiela | Learning Image Embeddings using Convolutional Neural Networks for Improved Multi-Modal Semantics[END_REF]: Concatenation of CNN and SKIP-GRAM embeddings. • VAE⊕SGNS: Concatenation of VAE and SKIP-GRAM embeddings.

• V-SGNS [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF]: A multimodal approach which augments SGNS with a term that treats CNN visual features as regression targets.

Comparisons with V-SGNS will allow us to evaluate the impact of our modeling assumptions. • IV-SGNS [START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF]: Learns a mapping from SGNS embeddings to CNN visual features. Due to a large degree of discrepancies in experimental setups across previously published methods and results,2 we re-implemented all the baselines and evaluate them under the same conditions.3 For [START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF], we implemented its model "A" as model "B" is comparable according to the original authors. For [START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF], we implemented both linear and nonlinear variants.

Evaluation metrics We use the cosine to measure the similarity between word representations. To assess the coherence between human ratings and models' predictions, we use the Spearman correlation coefficient.

Main results

The results across different datasets are shown in Table 2. We perform evaluations under two settings: by considering (i) word similarity between visual words only and (ii) between all words (column 100% in Table 2). For the models CNN, VAE and their concatenation with SGNS embeddings, the latter setting is not applicable. The two last rows correspond to the multimodal embeddings inferred from our model. In particular, PIXIE + (resp. PIXIE ⊕ ) represents the multimodal embeddings built using Eq. (10) (resp. Eq. ( 11)).

Overall, we note that PIXIE ⊕ offers the best performance in almost all situations. This provides strong empirical support for the proposed model. Below, we discuss the above results in more depth to better understand them and characterize the circumstances in which our model performs better.

How relevant is our formulation? Except PIXIE and V-SGNS, most of the multimodal competing methods rely on independently pre-computed linguistic embeddings. As Table 2 shows, PIXIE and V-SGNS are often the best performing multimodal models, which provides empirical evidence that accounting for perceptual information while learning word embeddings from text is beneficial. Moreover, the superior performance of PIXIE ⊕ over V-SGNS suggests that our model does a better job at combining perception and language to learn word representations.

Joint learning is beneficial PIXIE ⊕ outperforms VAE⊕SGNS in almost all cases, which demonstrates the importance of joint learning.

Where does our approach perform better? On datasets that focus on semantic/taxonomic similarity, our approach dominates all other methods.

On datasets focusing on general relatedness, our approach obtains mixed results. While dominating other approaches on MEN, it tends to perform worst than SGNS on MTurk and REL (under the 100% setting). One possible explanation is that general relatedness tends to focus more on "extrapolating" from one word to another word (such as SWAN is related to LAKE), while our approach better models more concrete relationships (such as SWAN is related to GOOSE). The low performance of CNN and VAE confirms this hypothesis.

On the VISSIM dataset focusing on visual similarity, both CNN⊕SGNS and VAE⊕SGNS perform the best, strongly suggesting that visual and linguistic data are complementary. Our approach comes very close to these two methods. Note that our learning objective is to jointly explain visual features and word-context co-occurrences. Thus, two visually similar words, which never occur within the same context, could be mapped into slightly different directions in the latent space.

Visual Propagation Here we wish to evaluate the ability of our model to propagate perceptual information to words lacking visual features. To * ) ). We then perform evaluations, under the two settings, on the datasets of Table 1 considering only pairs composed of words in the above subset of 2K words. As baselines for this experiment, we consider SGNS and the multimodal approaches which can propagate perceptual information, namely V-SGNS and IV-SGNS, as well as their outputs when the 2K visual features are ignored (denoted by V-SGNS ( * ) , and IV-SGNS ( * ) ).

The results are given in Table 3. We observe that PIXIE ⊕ ( * ) outperforms SGNS in almost all cases. Recall that, if we ignore the visual features for all words, PIXIE reduces to SGNS. We can therefore attribute the performance improvement of PIXIE ⊕ ( * ) over SGNS to the propagation of visual information to the subset of 2K words. Compared to multimodal methods, PIXIE ⊕ ( * ) (resp. PIXIE ⊕ ) performs better than V-SGNS ( * ) (resp. V-SGNS) and IV-SGNS ( * ) (resp. IV-SGNS) in almost situations. This suggests that our formulation allows perceptual information to propagate better. In Table 4, we report the cosine similarity between 5 semantically/visually coherent word pairs (from our subset of 2K words). Although the visual vectors of these words were removed during training, the PIXIE's word embeddings of each pair correlate better as compared to their SGNS counterparts, which provides further support to the propagation of visual information under PIXIE.

Qualitative analysis

Table 5 displays several qualitative examples of word similarity. We have selected 4 words: goose, brave, birthstone and savagery. The first two have visual feature representations in our training dataset and the last two do not. Furthermore, for each case, we chose one concrete and one abstract word. 4 For each word, we identify their nearest neighbors in the embedding space.

For the visual words, there is a noticeable difference between our method and others. For instance, for word goose, SGNS expresses more "general" relatedness and returns other animals like pig or shark, while our approach is more specific and tends to give visually similar neighbors by focusing on goose looks-like birds. V-SGNS's result is somewhat in between. On the abstract word brave, we observe that PIXIE ⊕ tends to select more explicit embodiments of the adjective brave than SGNS and V-SGNS.

Moving towards the non-visual words, we do not seem to find a consistent discrepancy pattern between V-SGNS and PIXIE ⊕ , though, as for visual words, both methods seem to select more explicit exemplars compared to SGNS. For instance, for the abstract word savagery, both multimodal approaches suggest cannibals and zombies. We now study the usefulness of the learned word embeddings for the tasks of image and caption retrieval. Our hypothesis is that multimodal word embeddings will perform better for downstream tasks involving multimodal information.

Experimental setup We use the Flickr30K dataset [START_REF] Young | From image descriptions to visual denotations: New similarity metrics for semantic inference over event descriptions[END_REF] containing 31,000 images and 155,000 sentences (5 captions per image). The sentences describe the images. The task is to identify the best sentence describing an image or to identify the best image depicting a sentence. We follow the data split setting provided by [START_REF] Karpathy | Deep visualsemantic alignments for generating image descriptions[END_REF], in which 1,000 images are used for validation and 1,000 for testing. The rest is used for training.

The retrieval models compute the proximity between the image features and the sentence embeddings. For image features, we use the precomputed features provided by [START_REF] Faghri | VSE++: Improving Visual-Semantic Embeddings with Hard Negatives[END_REF], which are extracted from the FC7 layer of VGG-19 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. These 4,096dimensional features are then linearly mapped to 1,024-dimensional features. For sentences, we use an one GRU-layer over the sequences of the word embeddings, resulting in 1,024-dimensional sentence embeddings.

We use a triplet loss to train the retrieval model such that the inner product between the corresponding image feature and the sentence is greater than the inner products with incorrect sentences (or im-ages) [START_REF] Kiros | Unifying visual-semantic embeddings with multimodal neural language models[END_REF]. The linear mapping and the GRU layer are then optimized to minimize the loss. We use the ADAM optimizer with the learning rate of 0.0002 and divide it by 10 every 15 epochs and we train the model for 30 epochs. Note that we do not fine-tune either the original visual feature or the word embeddings.

Results Table 6 summarizes the results. The evaluation metrics are accuracies at top-K (K=1, 5, or 10) retrieved sentences or images. Our model consistently outperforms SGNS and other competing multimodal methods, which provides additional support for the benefits of our approach.

Conclusion

We propose PIXIE, a novel probabilistic model joining textual and perceptual information to infer multimodal word embeddings. In our model, both linguistic and visual latent factors work in concert to explain the co-occurrences of words and their contexts in a corpus. Empirical results show that our model achieves equally competitive or stronger results when compared to state-of-the-art methods for multimodal embeddings.

Currently our model relies on unsupervised learning to infer visual factors. Explicit knowledge of similar and dissimilar visual categories could potentially disentangle latent factors better for alignment with linguistic data. How to incorporate visual domain knowledge more explicitly into the model would be an interesting direction for future research. While we build on skip-gram, the idea of PIXIE could be extended to other word embedding models, e.g., Glove [START_REF] Pennington | Glove: Global vectors for word representation[END_REF], ELMO (Peters et al., 2018), etc.

Figure 1 :

 1 Figure 1: Plate representation of our model PIXIE. The model consists of a generative model for visual data and a conditional model for text data. The latent visual factors z and text embedding e jointly predict the wordcontext pair's label y.

Table 1 :

 1 Datasets for the task of word similarity.

Table 2 :

 2 Results on word similarity task. Reported are the Spearman's rank order correlation between model prediction and human judgment (higher is better and bolds highlight the best methods). See text for details.

						Semantic/taxonomic similarity	General relatedness	Visual similarity REL+SIM
		Models		SEMSIM SimLex 100% 98% 100% 39% 100% 44% 100% 72% 100% 73% 100% 54% 100% 53% 100% 26% 100% 98% 100% 39% SIM EN-RG EN-MC MEN REL MTurk VISSIM WORDSIM
		CNN			-0.49 -0.41 -0.49 -0.54 -0.46 -0.54 -0.20 -0.18 -	0.53	-	0.28
		VAE			-0.65 -0.43 -0.51 -0.56 -0.55 -0.62 -0.22 -0.40 -	0.62	-	0.37
		SGNS			0.50 0.50 0.33 0.35 0.66 0.66 0.60 0.55 0.60 0.52 0.65 0.67 0.56 0.51 0.65 0.63 0.38	0.38	0.61 0.60
	CNN⊕SGNS	†	-0.67 -0.48 -0.65 -0.60 -0.55 -0.74 -0.44 -0.51 -	0.63	-	0.56
	VAE⊕SGNS	-0.70 -0.51 -0.67 -0.61 -0.60 -0.76 -0.45 -0.55 -	0.63	-	0.56
	V-SGNS	‡		0.58 0.58 0.29 0.30 0.66 0.71 0.73 0.73 0.69 0.69 0.64 0.65 0.51 0.52 0.60 0.65 0.42	0.42	0.59 0.64
	IV-SGNS	§ (LINEAR)	0.49 0.50 0.31 0.33 0.55 0.61 0.58 0.56 0.59 0.65 0.60 0.62 0.41 0.38 0.57 0.71 0.36	0.37	0.46 0.51
	IV-SGNS	§ (NONLINEAR) 0.44 0.44 0.30 0.32 0.53 0.59 0.54 0.53 0.59 0.63 0.57 0.59 0.40 0.37 0.56 0.71 0.32	0.33	0.44 0.48
		PIXIE+		0.63 0.63 0.35 0.48 0.63 0.72 0.65 0.60 0.62 0.62 0.64 0.73 0.46 0.56 0.55 0.55 0.54	0.54	0.50 0.59
		PIXIE⊕		0.71 0.71 0.39 0.53 0.68 0.71 0.73 0.73 0.69 0.71 0.68 0.76 0.52 0.59 0.60 0.59 0.60	0.61	0.58 0.65
						† :(Kiela and Bottou, 2014),

‡ :

Table 3 :

 3 Spearman's score on subsets of visual words. The symbol ( * ) indicates the visual features for the subset of 2K words have been ignored when training. Bold highlights the best performing method. Blue color highlights the best performing model under the ( * ) setting.

				Semantic/taxonomic	Relatedness Visual REL
												+SIM
	models	SEMSIM	SimLex	SIM	EN-RG	EN-MC	MEN	REL	MTurk	VISSIM	WORDSIM
	SGNS	0.55 0.40 0.67 0.57 0.52 0.68 0.53 0.75 0.41 0.63
	V-SGNS	0.61 0.32 0.72 0.71 0.69 0.69 0.52 0.75 0.44 0.66
	IV-SGNS (LINEAR) 0.50 0.40 0.60 0.59 0.60 0.66 0.41 0.72 0.37 0.54
	PIXIE⊕	0.71 0.59 0.72 0.73 0.71 0.77 0.57 0.74 0.60 0.62
	V-SGNS	( * )	0.58 0.29 0.72 0.68 0.61 0.67 0.44 0.76 0.43 0.60
	IV-SGNS	( * )	0.50 0.40 0.61 0.58 0.58 0.65 0.40 0.72 0.37 0.54
	PIXIE⊕	( * )	0.60 0.43 0.72 0.65 0.60 0.70 0.56 0.76 0.45 0.63
	this end, we randomly select a subset of 2, 000
	words for which we have visual features, and we
	train our model under two different settings: the
	visual features of the selected 2K words (i) are
	taken into account (PIXIE ⊕ ), (ii) are ignored, i.e.
	set to zero (PIXIE ⊕

Table 4 :

 4 Word pair cosine similarity computed based on SGNS and PIXIE ⊕ ( * ) embeddings.

	Word pairs	SGNS PIXIE⊕	( * )
	(chicken, turkey)	0.35	0.55
	(helicopter, jet)	0.63	0.76
	(falcon, hawk)	0.49	0.70
	(cathedral, chapel) 0.69	0.80
	(cup, mug)	0.39	0.46

Table 5 :

 5 10 nearest neighboring words to the query words in different embedding spaces generated by different methods. Only "visual" words contain direct visual representations in our dataset. The concreteness score of each query word is reported between parenthesis (see text for details).

		Query word	SGNS	V-SGNS	PIXIE⊕
	Visual	goose (4.81) brave (1.26)	quail, pig, shark, gull, smoky, sooty, bald, owl, guppy, bird courageous, young, man, fearless, heroic, thief, horrible, adventures, carefree, cowardly	fowl, quail, duck, bat, bulldog, puppy, warbler, blossom, wolfhound, crows heroes, valiant, ominous, wondrous, fearless, wanton, sabers, beast, excalibur, carefree	geese, duck, swans, swan, teal, loon, albatross, ostrich, gull, eider bodyguard, wives, benefactors, womanizer, heroes, valiant, immortal, housewives, fear-
					less, warrior
	Non Visual	birthstone (4.25) savagery (1.73)	heaviest, yeti, koala, snowfall, intrusions, mourning, amalthea, gleaming, incidentally, dolly swamp, man, crazy, madman, mysterious, thief, jeffrey, mad, hardy	emerald, lily, lavender, earthy, olive, delicacy, acacia, belladonna, flower, poppy mad, zombies, beast, fabulous, mysterious, nightmare, ghosts, alien, mayhem, bandits	beryl, emerald, lily, sandalwood, hula, guine-vere, pearls, holly, jasmine, jewels cannibals, evil, legends, werewolves, beast, haunting, ghosts, zombies, mayhem, thrillers

Table 6 :

 6 Results for image (I) ↔ sentence (S) retrieval.

	Models	I → S K=1 K=5 K=10 K=1 K=5 K=10 S → I
	SGNS	23.1 49.0 61.6 16.6 41.0 53.8
	V-SGNS	21.9 51.7 64.2 16.2 42.0 54.8
	IV-SGNS (LINEAR) 22.7 50.5 61.7 17.1 42.6 55.4
	PIXIE+	24.2 52.5 65.4 17.5 43.8 56.2
	PIXIE⊕	25.7 55.7 67.7 18.4 44.9 56.9
	4.3 Task 2: Image and Caption Retrieval

For instance,[START_REF] Lazaridou | Combining Language and Vision with a Multimodal Skip-gram Model[END_REF] reports only 5,100 visual words, nearly half of what we have defined.[START_REF] Collell | Imagined Visual Representations as Multimodal Embeddings[END_REF] used pre-trained GloVe word vectors obtained from a different corpus.

For each method we use the hyper-parameters recommended by the authors.

We rely on the concreteness ratings made available by[START_REF] Brysbaert | Concreteness ratings for 40 thousand generally known english word lemmas[END_REF], ranging from 1 to

Acknowledgements

This work was partially supported by the virtual lab Inria@SiliconValley, ANR Grant GRASP No. ANR-16-CE33-0011-01, and a grant from CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020, NSF IIS-1065243, 1451412, 1513966/ 1632803/1833137, 1208500, CCF-1139148, a Google Research Award, an Alfred P. Sloan Research Fellowship, gifts from Facebook and Netflix, and ARO# W911NF-12-1-0241 and W911NF-15-1-0484. We thank anonymous reviewers for their suggestions and comments.


Wenjie Zheng

Aurélien Bellet

Patrick Gallinari

A Distributed Frank-Wolfe Framework for Learning Low-Rank Matrices with the Trace Norm

Keywords: Frank-Wolfe algorithm, Low-rank learning, Trace norm, Distributed optimization, Multi-task learning, Multinomial logistic regression

We consider the problem of learning a high-dimensional but low-rank matrix from a large-scale dataset distributed over several machines, where lowrankness is enforced by a convex trace norm constraint. We propose DFW-Trace, a distributed Frank-Wolfe algorithm which leverages the low-rank structure of its updates to achieve efficiency in time, memory and communication usage. The step at the heart of DFW-Trace is solved approximately using a distributed version of the power method. We provide a theoretical analysis of the convergence of DFW-Trace, showing that we can ensure sublinear convergence in expectation to an optimal solution with few power iterations per epoch. We implement DFW-Trace in the Apache Spark distributed programming framework and validate the usefulness of our approach on synthetic and real data, including the ImageNet dataset with high-dimensional features extracted from a deep neural network.

Introduction

Learning low-rank matrices is a problem of great importance in machine learning, statistics and computer vision. Since rank minimization is known to be NP-hard, a principled approach consists in solving a convex relaxation of the problem where the rank is replaced by the trace norm (also known as the nuclear norm) of the matrix. This strategy is supported by a range of theoretical results showing that when the ground truth matrix is truly low-rank, one can recover it exactly (or accurately) from limited samples and under mild conditions (see [START_REF] Bach | Consistency of trace norm minimization[END_REF][START_REF] Candès | Exact matrix completion via convex optimization[END_REF][START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF][START_REF] Recht | A simpler approach to matrix completion[END_REF][START_REF] Gross | Quantum state tomography via compressed sensing[END_REF][START_REF] Gross | Recovering low-rank matrices from few coefficients in any basis[END_REF][START_REF] Koltchinskii | Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion[END_REF]Bhojanapalli et al, 2016). Trace norm minimization has led to many successful applications, among which collaborative filtering and recommender systems [START_REF] Koren | Matrix factorization techniques for recommender systems[END_REF], multi-task learning [START_REF] Argyriou | Convex multi-task feature learning[END_REF][START_REF] Pong | Trace norm regularization: Reformulations, algorithms, and multi-task learning[END_REF], multi-class and multi-label classification [START_REF] Goldberg | Transduction with matrix completion: Three birds with one stone[END_REF][START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF][START_REF] Harchaoui | Large-scale image classification with trace-norm regularization[END_REF], robust PCA [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF], phase retrieval [START_REF] Candes | Phase retrieval via matrix completion[END_REF] and video denoising [START_REF] Ji | Robust video denoising using low rank matrix completion[END_REF].

We consider the following generic formulation of the problem: min

W ∈R d×m F (W ) = n i=1 f i (W ) s.t. W * ≤ µ, (1)
where the f i 's are differentiable with Lipschitz-continuous gradient, W * = k σ k (W ) is the trace norm of W (the sum of its singular values), and µ > 0 is a regularization parameter (typically tuned by cross-validation). In a machine learning context, an important class of problems considers f i (W ) to be a loss value which is small (resp. large) when W fits well (resp. poorly) the i-th data point (see Section 2.3 for concrete examples). 1 In this work, we focus on the large-scale scenario where the quantities involved in (1) are large: typically, the matrix dimensions d and m are both in the thousands or above, and the number of functions (data points) n is in the millions or more.

Various approaches have been proposed to solve the trace norm minimization problem (1). 2 One can rely on reformulations as semi-definite programs and use out-of-the-shelf solvers such as SDPT3 [START_REF] Toh | SDPT3-a MATLAB software package for semidefinite programming, version 1.3[END_REF] or SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF], but this does not scale beyond small-size problems. To overcome this limitation, first-order methods like Singular Value Thresholding [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], Fixed Point Continuation algorithms [START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF] and more generally projected/proximal gradient algorithms [START_REF] Parikh | Proximal algorithms[END_REF] have been proposed. These approaches have two important drawbacks preventing their use when the matrix dimensions d and m are both very large: they require to compute a costly (approximate) SVD at each iteration, and their memory complexity is O(dm). In this context, Frank-Wolfe (also known as conditional gradient) algorithms [START_REF] Frank | An algorithm for quadratic programming[END_REF]) provide a significant reduction in computational and memory complexity: they only need to compute the leading eigenvector at each iteration, and they maintain compact low-rank iterates throughout the optimization [START_REF] Hazan | Sparse approximate solutions to semidefinite programs[END_REF][START_REF] Jaggi | A simple algorithm for nuclear norm regularized problems[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF][START_REF] Harchaoui | Conditional gradient algorithms for norm-regularized smooth convex optimization[END_REF]. However, as all firstorder algorithms, Frank-Wolfe requires to compute the gradient of the objective function at each iteration, which requires a full pass over the dataset and becomes a bottleneck when n is large.

The goal of this paper is to propose a distributed version of the Frank-Wolfe algorithm in order to alleviate the cost of gradient computation when solving problem (1). We focus on the Bulk Synchronous Parallel (BSP) model with a master node connected to a set of slaves (workers), each of the workers having access to a subset of the f i 's (typically corresponding to a subset of training points). Our contributions are three-fold. First, we propose DFW-Trace, a Frank-Wolfe algorithm relying on a distributed power method to approximately compute the 1 More general cases can be addressed, such as pairwise loss functions f i,j corresponding to pairs of data points.

2 Some methods consider an equivalent formulation where the trace norm appears as a penalization term in the objective function rather than as a constraint.

Algorithm 1 Centralized Frank-Wolfe algorithm to solve (2)

Input: Initial point W 0 ∈ D, number of iterations T for t = 0, . . . , T -1 do S * ← arg min S∈D S, ∇F (W t ) solve linear subproblem γ t ← 2 t+2 (or determined by line search) step size W t+1 ← (1 -γ t )W t + γ t S *
update end for Output: W T leading eigenvector with communication cost of O(d+m) per pass over the dataset (epoch). This dramatically improves upon the O(dm) cost incurred by a naive distributed approach. Second, we prove the sublinear convergence of DFW-Trace to an optimal solution in expectation, quantifying the number of power iterations needed at each epoch. This result guarantees that DFW-Trace can find low-rank matrices with small approximation error using few power iterations per epoch. Lastly, we provide a modular implementation of our approach in the Apache Spark programming framework [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF] which can be readily deployed on commodity and commercial clusters. We evaluate the practical performance of DFW-Trace by applying it to multi-task regression and multi-class classification tasks on synthetic and real-world datasets, including the ImageNet database [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF] with high-dimensional features generated by a deep neural network. The results confirm that DFW-Trace has fast convergence and outperforms competing methods. While distributed FW algorithms have been proposed for other classes of problems [START_REF] Bellet | A Distributed Frank-Wolfe Algorithm for Communication-Efficient Sparse Learning[END_REF][START_REF] Moharrer | Distributing Frank-Wolfe via Map-Reduce[END_REF][START_REF] Wang | Parallel and distributed block-coordinate Frank-Wolfe algorithms[END_REF], to the best of our knowledge our work is the first to propose, analyze and experiment with a distributed Frank-Wolfe algorithm designed specifically for trace norm minimization.

The rest of this paper is organized as follows. Section 2 introduces some background on the (centralized) Frank-Wolfe algorithm and its specialization to trace norm minimization, and reviews some applications. After presenting some baseline approaches for the distributed setting, Section 3 describes our algorithm DFW-Trace and its convergence analysis, as well as some implementation details. Section 4 discusses some related work, and Section 5 presents the experimental results.

Background

We review the centralized Frank-Wolfe algorithm in Section 2.1 and its specialization to trace norm minimization in Section 2.2. We then present some applications to multi-task learning and multi-class classification in Section 2.3.

Frank-Wolfe Algorithm

The original Frank-Wolfe (FW) algorithm dates back from the 1950s and was originally designed for quadratic programming [START_REF] Frank | An algorithm for quadratic programming[END_REF]. The scope of the algorithm was then extended to sparse greedy approximation [START_REF] Clarkson | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF] and semi-definite programming [START_REF] Hazan | Sparse approximate solutions to semidefinite programs[END_REF]. Recently, [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF] generalized the algorithm further to tackle the following generic problem:

min W ∈D F (W ) , (2)
where F is convex and continuously differentiable, and the feasible domain D is a compact convex subset of some Hilbert space with inner product •, • .

Algorithm 1 shows the generic formulation of the FW algorithm applied to (2). At each iteration t, the algorithm finds the feasible point S * ∈ D which minimizes the linearization of F at the current iterate W t . The next iterate W t+1 is then obtained by a convex combination of W t and S * , with a relative weight given by the step size γ t . By convexity of D, this ensures that W t+1 is feasible. The algorithm converges in O(1/t), as shown by the following result from [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF].

Theorem 1 [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF]

Let C F be the curvature constant of F . 3 For each t ≥ 1, the iterate W t ∈ D generated by Algorithm 1 satisfies F (W t ) -F (W * ) ≤ 2C F t+2
, where W * ∈ D is an optimal solution to (2).

Remark 1 There exist several variants of the FW algorithm, for which faster rates can sometimes be derived under additional assumptions. We refer to [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF], and Lacoste-Julien and [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF] for details.

From the algorithmic point of view, the main step in Algorithm 1 is to solve the linear subproblem over the domain D. By the linearity of the subproblem, a solution always lies at an extremal point of D, hence FW can be seen as a greedy algorithm whose iterates are convex combinations of extremal points (adding a new one at each iteration). When these extremal points have some specific structure (e.g., sparsity, low-rankness), the iterates inherit this structure and the linear subproblem can sometimes be solved very efficiently. This is the case for the trace norm constraint, our focus in this paper.

Specialization to Trace Norm Minimization

The FW algorithm applied to the trace norm minimization problem (1) must solve the following subproblem:

S * ∈ arg min S * ≤µ S, ∇F (W t ) , (3)
where W t ∈ R d×m is the iterate at time t and S ∈ R d×m . The trace norm ball is the convex hull of the rank-1 matrices, so there must exist a rank-1 solution to (3). This solution can be shown to be equal to -µu 1 v 1 , where u 1 and v 1 are the unit left and right top singular vectors of the gradient matrix ∇F (W t ) [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF].

Finding the top singular vectors of a matrix is much more efficient than computing the full SVD. This gives FW a significant computational advantage over projected and proximal gradient descent approaches when the matrix dimensions are large. Furthermore, assuming that W 0 is initialized to the zero matrix, W t can be stored in a compact form as a convex combination of t rank-1 matrices, which requires O(t(d + m)) memory instead of O(dm) to store a full rank matrix. As implied by Theorem 1, FW is thus guaranteed to find a rank-t whose approximation error is O(1/t) for any t ≥ 1. In practice, when the ground truth matrix is indeed low-rank, FW can typically recover a very accurate solution after t min(d, m) steps. We note that in the special case where the matrix W is square (d = m) and constrained to be symmetric, the gradient ∇F (W t ) can always be written as a symmetric matrix, and the solution to the linear subproblem has a simpler representation based on the leading eigenvector of the gradient, see Jaggi (2013).

Applications

We describe here two tasks where trace norm minimization has been successfully applied, which we will use to evaluate our method in Section 5.

Multi-task least square regression. This is an instance of multi-task learning [START_REF] Caruana | Multitask Learning[END_REF], where one aims to jointly learn m related tasks. Formally, let X ∈ R n×d be the feature matrix (n training points in d-dimensional space) and Y ∈ R n×m be the response matrix (each column corresponding to a task). The objective function aims to minimize the residuals of all tasks simultaneously:

F (W ) = 1 2 XW -Y 2 F = 1 2 n i=1 m j=1 (x T i w j -y ij ) 2 , ( 4
)
where • F is the Frobenius norm. Using a trace norm constraint on W allows to couple the tasks together by making the task predictors share a common subspace, which is a standard approach to multi-task learning (see e.g., [START_REF] Argyriou | Convex multi-task feature learning[END_REF][START_REF] Pong | Trace norm regularization: Reformulations, algorithms, and multi-task learning[END_REF].

Multinomial logistic regression. Consider a classification problem with m classes. Let X ∈ R n×d be the feature matrix and y ∈ {1, . . . , m} n the label vector. Multinomial logistic regression minimizes the negative log-likelihood function:

F (W ) = i log 1+ l =y i exp(w T l x i -w T y i x i ) = i -w T y i x i +log l exp(w T l x i ) .
(5) The motivation for using the trace norm is that multi-class problems with a large number of categories usually exhibit low-rank embeddings of the classes (see [START_REF] Amit | Uncovering shared structures in multiclass classification[END_REF][START_REF] Harchaoui | Large-scale image classification with trace-norm regularization[END_REF].

Distributed Frank-Wolfe for Trace Norm Minimization

We now consider a distributed master/slave architecture with N slaves (workers). The master node is connected to all workers and acts mainly as an aggregator, while most of the computation is done on the workers. The individual functions f 1 , . . . , f n in the objective (1) are partitioned across workers, so that all workers can collectively compute all functions but each worker can only compute its own subset. Recall that in a typical machine learning scenario, each function f i corresponds to the loss function computed on the i-th data point (as in the examples of Section 2.3). We will thus often refer to these functions as data points. Formally, let I j ⊆ {1, . . . , n} be the set of indices assigned to worker j, where

I 1 ∪ • • • ∪ I N = {1, . . . , n} and I 1 ∩ • • • ∩ I N = ∅.
We denote by F j = i∈I j f i the local function (dataset) associated with each worker j, and by n j = |I j | the size of this local dataset.

We follow the Bulk Synchronous Parallel (BSP) computational model: each iteration (epoch) alternates between parallel computation at the workers and communication with the master (the latter serves as a synchronization barrier).

Baseline Strategies

Before presenting our algorithm, we first introduce two baseline distributed FW strategies (each with their own merits and drawbacks).

Naive DFW. One can immediately see a naive way of running the centralized Frank-Wolfe algorithm (Algorithm 1) in the distributed setting. Starting from a common initial point W 0 , each worker j computes at each iteration t its local gradient ∇F j (W t ) and sends it to the master. The master then aggregates the messages to produce the full gradient ∇F (W t ) = N j=1 F j (W t ), solves the linear subproblem by computing the leading right/left singular vectors of ∇F (W t ) and sends the solution back to the workers, who can form the next iterate W t+1 . Naive-DFW exactly mimics the behavior of the centralized FW algorithm, but induces a communication cost of O(N dm) per epoch as in many applications (such as those presented in Section 2.3) the local gradients are dense matrices. In the large-scale setting where the matrix dimensions d and m are both large, this cost dramatically limits the efficiency of the algorithm. Singular Vector Averaging. A possible strategy to avoid this high communication cost is to ask each worker j to send to the master the rank-1 solution to the local version of the subproblem (3), in which they use their local gradient ∇F j (W t ) as an estimate of the full gradient ∇F (W t ). This reduces the communication to a much more affordable cost of O(N (d + m)). Note that averaging the rank-1 solutions would typically lead to a rank-N update, which breaks the useful rank-1 property of FW and is undesirable when N is large. Instead, the master averages the singular vectors (weighted proportionally to n j ), resolving the sign ambiguity by setting the largest entry of each singular vector to be positive and using appropriate normalization, as mentioned for instance in [START_REF] Bro | Resolving the sign ambiguity in the singular value decomposition[END_REF]. We refer to this strategy as Singular Vector Averaging (SVA). SVA is a reasonable heuristic when the individual functions are partitioned across nodes uniformly at random: in this case the local gradients can be seen as unbiased estimates of the full gradient. However the singular vector estimate itself is biased (averaging between workers only reduces its variance), and for n fixed this bias increases with the matrix dimensions d and m but also with the number of workers N (which is not a desirable property in the distributed setting). It is also expected to perform badly on arbitrary (non-uniform) partitions of functions across workers. Clearly, one cannot hope to establish strong convergence guarantees for SVA.

Algorithm 2 Our distributed algorithm DFW-Trace to solve (1)

1: Input: Initial point W 0 ∈ D, number of iterations T 2: for t = 0, . . . , T -1 do 3:
Each worker j: ∇F j (W t ) ← i∈I j ∇f i (W t ) 4:

All workers: draw the same v 0 ∈ R m uniformly on unit sphere 5:

for k = 0, . . . , K(t) -1 do distributed power method 6:

Each worker j: send u k+1,j ← ∇F j (W t )v k to master 7:

Master: broadcast u k+1 ← ( N j=1 u k+1,j )/ N j=1 u k+1,j 8:

Each worker j: send v k+1,j ← ∇F j (W t ) u k+1 to master 9:

Master: broadcast v k+1 ← ( N j=1 v k+1,j )/ N j=1 v k+1,j 10:
end for 11:

γ t ← 2 t+2 (or determined by line search) step size 12:

Each worker j:

W t+1 ← (1 -γ t )W t -γ t µu K(t) v K(t
(d + m) 1 DFW-Trace 2N K(t)(d + m) 2K(t)

Proposed Approach

We now describe our proposed approach, referred to as DFW-Trace. We will see that DFW-Trace achieves roughly the small communication cost of SVA while enjoying a similar convergence rate as Naive-DFW (and hence centralized FW).

Algorithm. The main idea of DFW-Trace (Algorithm 2) is to solve the linear subproblem of FW approximately using a distributed version of the power method applied to the matrix ∇F (W t ) F (W t ). At each outer iteration (epoch) t, the workers first generate a common random vector drawn uniformly on the unit sphere. 4 Then, for K(t) iterations, the algorithm alternates between the workers computing matrix-vector products and the master aggregating the results. At the end of this procedure, workers hold the same approximate versions of the left and right singular vectors of ∇F (W t ) and use them to generate the next iterate 1 for a comparison with baselines). It is clear that as K(t) → ∞, DFW-Trace computes the exact solution to the linear subproblems and hence has the same convergence guarantees as centralized FW. However, we would like to set K(t) min(d, m) to provide a significant improvement over the O(N dm) cost of the naive distributed algorithm. The purpose of our analysis below is to show how to set K(t) to preserve the convergence rate of the centralized algorithm.

W t+1 . The communication cost of DFW-Trace per epoch is O(N K(t)(d + m)) (see Table
Remark 2 (Other network topologies) Since any connected graph can be logically represented as a star graph by choosing a center, our method virtually works on any network (though it may incur additional communication). Depending on the topology, special care can be taken to reduce the communication overhead. An interesting case is the rooted tree network: we can adopt a hierarchical aggregation scheme which has the same communication cost of O(N K(t)(d + m)) as the star network but scales better to many workers by allowing parallel aggregations.5 For a general graph with M edges, O(M K(t)(d + m)) communication is enough to broadcast the values to all workers so they can perform the aggregation locally.

Analysis. We will establish that for some appropriate choices of K(t), DFW-Trace achieves sublinear convergence in expectation, as defined below.

Definition 1 Let δ ≥ 0 be an accuracy parameter. We say that DFW-Trace converges sublinearly in expectation if for each t ≥ 1, its iterate W t satisfies

E[F (W t )] -F (W * ) ≤ 2C F t+2 (1 + δ), (6)
where C F is the curvature constant of F .

We have the following result.

Theorem 2 (Convergence) Let F be a convex, differentiable function with curvature C F and Lipschitz constant L w.r.t. the trace norm. For any accuracy parameter δ ≥ 0, the following properties hold for DFW-Trace (Algorithm 2):

1. If m ≥ 8 and for any t ≥ 0, K(t) ≥ 1 + µL(t+2) ln m δC F

, then DFW-Trace converges sublinearly in expectation. 2. For any t ≥ 0, let σ t 1 , σ t 2 be the largest and the second largest singular values of ∇F (W t ) and assume that σ t 1 has multiplicity 1 and there exists a constant β such that

σ t 2 σ t 1 < β < 1. If K(t) ≥ max( ln(δC F )-ln[mµL(t+2)] 2 ln β + 1, K)
where K is a large enough constant, DFW-Trace converges sublinearly in expectation.

Proof (Sketch) We briefly outline the main ingredients (see Appendix A for the detailed proof). We first show that if the linear subproblem is approximately solved in expectation (to sufficient accuracy), then the FW algorithm converges sublinearly in expectation. Relying on results on the convergence of the power method [START_REF] Kuczyński | Estimating the largest eigenvalue by the power and Lanczos algorithms with a random start[END_REF] and on the Lipschitzness of F , we then derive the above results on the number of power iterations K(t) needed to ensure sufficient accuracy under different assumptions.

Theorem 2 characterizes the number of power iterations K(t) at each epoch t which is sufficient to guarantee that DFW-Trace converges sublinearly in expectation to an optimal solution. Note that there are two regimes. The first part of the theorem establishes that if K(t) scales linearly in t, the expected output of DFW-Trace after t epochs is a rank-t matrix with O(1/t) approximation error (as in centralized FW, see Theorem 1). In the large-scale setting of interest, this implies that a good low-rank approximation can be achieved by running the algorithm for t min(d, m) iterations, and with reasonable communication cost since K(t) = O(t). Remarkably, this result holds without any assumption about the spectral structure of the gradient matrices. On the other hand, in the regime where the gradient matrices are "well-behaved" (in the sense that the ratio between their two largest singular values is bounded away from 1), the second part of the theorem shows that a much lower number of power iterations K(t) = O(log t) is sufficient to ensure the sublinear convergence in expectation. In Section 5, we will see experimentally on several datasets that this is indeed sufficient in practice to achieve convergence. We conclude this part with a few remarks mentioning some additional results, for which we omit the details due to the lack of space.

Remark 3 (Convergence in probability) We can also establish the sublinear convergence of DFW-Trace in probability (which is stronger than convergence in expectation). The results are analogous to Theorem 1 but require K(t) to be quadratic in t for the first case, and linear in t for the second case.

Remark 4 (Constant number of power iterations) If we take the number of power iterations to be constant across epochs (i.e., K(t) = K for all t), DFW-Trace converges in expectation to a neighborhood of the optimal solution whose size decreases with K. We can establish this by combining results on the approximation error of the power method with Theorem 5.1 in [START_REF] Freund | New analysis and results for the Frank-Wolfe method[END_REF].

Implementation

Our algorithm DFW-Trace (Algorithm 2) can be implemented as a sequence of map-reduce steps [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF]. This allows the computation to be massively parallelized across the set of workers, while allowing a simple implementation and fast deployment on commodity and commercial clusters via existing distributed programming frameworks [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF][START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF]. Nonetheless, some special care is needed if one wants to get an efficient implementation. In particular, it is key to leverage the fundamental property of FW algorithms that the updates are rank-1. This structural property implies that it is much more efficient to compute the gradient in a recursive manner, rather than from scratch using the current parameters. We use a notion of sufficient information to denote the local quantities (maintained by each worker) that are sufficient to compute the updates. This includes the local gradient (for the reason outlined above), and sometimes some quantities precomputed from the local dataset. Depending on the objective function and the relative size of the problem parameters n, m, d and N , the memory and/or time complexity may be improved by storing (some of) the sufficient information in low-rank form. We refer the reader to Appendix B for a concrete application of these ideas to the tasks of multi-task least square regression and multinomial logistic regression used in our experiments. Based on the above principles, we developed an open-source Python implementation of DFW-Trace using the Apache Spark framework [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF]. 6 The package also implements the baseline strategies of Section 3.1, and currently uses dense representations. The code is modular and separates generic from task-specific components. In particular, the generic DFW-Trace algorithm is implemented in PySpark (Spark's Python API) in a task-agnostic fashion. On the other hand, specific tasks (objective function, gradient, etc) are implemented separately in pure Python code. This allows users to easily extend the package by adding their own tasks of interest without requiring Spark knowledge. Specifically, the task interface should implement several methods: stats (to initialize the sufficient information), update (to update the sufficient information), and optionally linesearch (to use linesearch instead of default step size) and loss (to compute the value of the objective function). In the current version, we provide such interface for multi-task least square regression and multinomial logistic regression.

Related Work

There has been a recent surge of interest for the Frank-Wolfe algorithm and its variants in the learning community. The renewed popularity of this classic algorithm, introduced by [START_REF] Frank | An algorithm for quadratic programming[END_REF], can be largely attributed to the work of [START_REF] Clarkson | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF] and more recently [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF]. They generalized its scope and showed that its strong convergence guarantees, efficient greedy updates and sparse iterates are valuable to tackle high-dimensional machine learning problems involving sparsity-inducing (non-smooth) regularization such as the L 1 norm and the trace norm. Subsequent work has extended the convergence results, for instance proving faster rates under some additional assumptions (see Lacoste-Julien and [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF][START_REF] Garber | Faster Rates for the Frank-Wolfe Method over Strongly-Convex Sets[END_REF][START_REF] Freund | New analysis and results for the Frank-Wolfe method[END_REF].

As first-order methods, FW algorithms rely on gradients. In machine learning, computing the gradient of the objective typically requires a full pass over the dataset. To alleviate this computational cost on large datasets, some distributed versions of FW algorithms have recently been proposed for various problems. Bellet et al ( 2015) introduced a communication-efficient distributed FW algorithm for a class of problems under L 1 norm and simplex constraints, and provided an MPIbased implementation. [START_REF] Tran | Distributed frank-wolfe under pipelined stale synchronous parallelism[END_REF] extend the algorithm to the Stale Synchronous Parallel (SSP) model. [START_REF] Moharrer | Distributing Frank-Wolfe via Map-Reduce[END_REF] further generalized the class of problems which can be considered (still under L 1 /simplex constraints) and proposed an efficient and modular implementation in Apache Spark (similar to what we propose in the present work for trace norm problems). [START_REF] Wang | Parallel and distributed block-coordinate Frank-Wolfe algorithms[END_REF] proposed a parallel and distributed version of the Block-Coordinate Frank-Wolfe algorithm (Lacoste-Julien et al, 2013) for problems with block-separable constraints. All these methods are designed for specific problem classes and do not apply to trace norm minimization. For general problems (including trace norm minimization), [START_REF] Wai | Decentralized Frank-Wolfe Algorithm for Convex and Non-convex Problems[END_REF] recently introduced a decentralized FW algorithm in which workers communicate over a network graph without master node. The communication steps involve local averages of iterates and gradients between neighboring workers. In the master/slave distributed setting we consider, their algorithm essentially reduces to the naive distributed FW described in Section 3.1 and hence suffers from the high communication cost induced by transmitting gradients. In contrast to the above approaches, our work proposes a communicationefficient distributed FW algorithm for trace norm minimization.

Another direction to scale up FW algorithms to large datasets is to consider stochastic variants, where the gradient is replaced by an unbiased estimate computed on a mini-batch of data points [START_REF] Hazan | Projection-free Online Learning[END_REF][START_REF] Lan | Conditional Gradient Sliding for Convex Optimization[END_REF][START_REF] Hazan | Variance-Reduced and Projection-Free Stochastic Optimization[END_REF]. The price to pay is a slower theoretical convergence rate, and in practice some instability and convergence issues have been observed (see e.g., [START_REF] Liu | Approximate Conditional Gradient Descent on Multi-Class Classification[END_REF]. The experimental results of [START_REF] Moharrer | Distributing Frank-Wolfe via Map-Reduce[END_REF] show that current stochastic FW approaches do not match the performance of their distributed counterparts. Despite these limitations, this line of work is complementary to ours: when the number of workers N is small compared to the training set size n, each worker could compute an estimate of its local gradient to further reduce the computational cost. We leave this for future work.

We conclude this section by mentioning that other kinds of distributed algorithms have been proposed for special cases of our general problem (1). In particular, for the matrix completion problem, [START_REF] Mackey | Divide-and-conquer matrix factorization[END_REF] proposed a divide-and-conquer strategy, splitting the input matrix into submatrices, solving each subproblem in parallel with an existing matrix completion algorithm, and then combining the results.

Experiments

In this section, we validate the proposed approach through experiments on two tasks: multi-task least square regression and multinomial logistic regression (see Section 2.3). We use both synthetic and real-world datasets.

Experimental Setup

Environment. We run our Spark implementation described in Section 3.3 on a cluster with 5 identical machines, with Spark 1.6 deployed in standalone mode. One machine serves as the driver (master) and the other four as executors (workers). Each machine has 2 Intel Xeon E5645 2.40GHz CPUs, each with 6 physical cores. Each physical core has 2 threads. Therefore, we have 96 logical cores available as workers. The Spark cluster is configured to use all 96 logical cores unless otherwise stated. Each machine has 64GB RAM: our Spark deployment is configured to use 60GB, hence the executors use 240GB in total. The network card has a speed of 1Gb/s. The BLAS version does not enable multi-threading.

Datasets. For multi-task least square, we experiment on synthetic data generated as follows. The ground truth W has rank 10 and trace norm equal to 1 (we thus set µ = 1 in the experiments). This is obtained by multiplying two arbitrary orthogonal matrices and a sparse diagonal matrix. X is generated randomly, with each coefficient following a Gaussian distribution, and we set Y = XW . We generate two versions of the dataset: a low-dimensional dataset (n = 10 5 samples, d = 300 features and m = 300 tasks) and a higher dimensional one (n = 10 5 , d = 1, 000 and m = 1, 000). For multinomial logistic regression, we use a synthetic and a real dataset. The synthetic dataset has n = 10 5 samples, p = 1, 000 features and m = 1, 000 classes. The generation of W and X is the same as above, with the label vector y set to the one yielding the highest score for each point. The test set has 10 5 samples. Our real-world dataset is ImageNet from ILSVRC2012 challenge [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF][START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], which has n = 1, 281, 167 training images in m = 1, 000 classes. We use the learned features of dimension p = 2048 extracted from the deep neural network ResNet50 [START_REF] He | Deep residual learning for image recognition[END_REF] provided by Keras. 7 The validation set of the competition (50, 000 images) serves as the test set.

Compared methods. We compare the following algorithms: Naive-DFW, SVA (the baselines described in Section 3.1) and three variants of our algorithm, DFW-Trace-1, DFW-Trace-2 and DFW-Trace-log (resp. using 1, 2 and O(log t) power iterations at step t). We have also experimented with DFW-Trace with K(t) = O(t), but observed empirically that far fewer power iterations are sufficient in practice to ensure good convergence. We have also used SVA as warm start to the power iterations within DFW-Trace, which marginally improves the performance of DFW-Trace. We do not show these variants on the figures for clarity.

Results

Multi-task least square. For this task, we simply set the number of power iterations of DFW-Trace-log to K(t) = 1 + log(t) . All algorithms use line search.

Figure 1 shows the results for all methods on the low and high-dimensional versions of the dataset. The performance is shown with respect to the number of epochs and runtime, and for two metrics: the value of the objective function and the estimation error (relative Frobenius distance between the current W and the ground truth). On this dataset, the estimation error behaves similarly as the objective function. As expected, Naive-DFW performs the best with respect to the number of epochs as it computes the exact solution to the linear subproblem. On the low-dimensional dataset (left panel), it also provides the fastest decrease in objective/error. SVA also performs well on this dataset. However, when the dimension grows (right panel) the accuracy of SVA drops dramatically and Naive-DFW becomes much slower due to the increased communication cost. This confirms that these baselines do not scale well with the matrix dimensions. On the other hand, all variants of DFW-Trace perform much better than the baselines on the higher-dimensional dataset. This gap is expected to widen as the matrix dimensions increase. Remarkably, only 2 power iterations are sufficient to closely match the reduction in objective function achieved by the exact solution on this task. One can see the influence of the number of power iterations on the progress per epoch (notice for instance the clear break at iteration 10 when DFW-Trace-log switches from 1 to 2 power iterations), but this has a cost in terms of runtime.

Overall, all variants of DFW-Trace reduce the objective/error at roughly the same speed. On a smaller scale version of the dataset, we verified that the gradients are well-behaved in the sense of Theorem 2: the average ratio between the two largest singular values over 100 epochs was found to be 0.86.

Multinomial logistic regression.

Here, all algorithms use a fixed step size as there is no closed-form line search. As we observed empirically that this task requires a larger number of FW iterations to converge, we set K(t) = 1+0.5 log(t) for DFW-Trace-log so that the number of power iterations does not exceed 2 as in the previous experiment. Figure 2 shows the results on the synthetic dataset for several values ofµ (the upper bound on the trace norm). They are consistent with those obtained for multi-task least square. In particular, SVA achieves converges to a suboptimal solution, while Naive-DFW converges fast in terms of epochs but its runtime is larger than DFW-Trace. DFW-Trace-2 and DFW-Tracelog perform well across all values of µ: this confirms that very few power iterations are sufficient to ensure good convergence. For more constrained problems (µ = 10), the error does not align very well with the objective function and hence optimizing the subproblems to lower accuracy with DFW-Trace-1 works best. We now turn to the ImageNet dataset. The results for µ = 30 with 24 cores are shown on Figure 3.8 Again, the DFW-Trace variants clearly outperform Naive-DFW and SVA. While DFW-Trace-2 and DFW-Trace-log reduce the objective value faster than DFW-Trace-1, the latter reduces the error slightly faster. When run until convergence, all variants converge to state-of-the-art top-5 misclassification rate with these features (around 0.13, on par with the pre-trained deep neural net provided by Keras).

We conclude these experiments by investigating the speed-ups obtained when varying the number of cores on the ImageNet dataset. As seen on the left panel of Figure 4, the time per epoch nicely decreases with the number of cores (with diminishing returns, as expected in distributed computing). The right panel of Figure 4 illustrates this effect on the convergence speed for DFW-Trace-1.

Conclusion

In this work, we introduced a distributed Frank-Wolfe algorithm for learning highdimensional low-rank matrices from large-scale datasets. Our DFW-Trace algo-rithm is communication-efficient, enjoys provable convergence rates and can be efficiently implemented in map-reduce operations. We implemented DFW-Trace as a Python toolbox relying on the Apache Spark distributed programming framework, and showed that it performs well on synthetic and real datasets.

In future work, we plan to investigate several directions. First, we would like to study whether faster theoretical convergence can be achieved under additional assumptions. Second, we wonder whether our algorithm can be deployed in GPUs and be used in neural networks with back-propagated gradients. Finally, we hope to explore how to best combine the ideas of distributed and stochastic Frank-Wolfe algorithms.

Lemma 1 Let δ ≥ 0 be an accuracy parameter. If at each step t ≥ 0, the linear subproblem is approximately solved in expectation, i.e. we find a random variable Ŝ such that

E[ Ŝ|W t ], ∇F (W t ) ≤ min S∈D S, ∇F (W t ) + 1 2 δγ t C F , (7)
then the Frank-Wolfe algorithm converges sublinearly in expectation.

Proof At any step t, given W t we set W t+1 = W t +γ t ( Ŝ -W t ) with arbitrary step size γ t ∈ [0, 1]. From the definition of the curvature constant C F [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF]:

F (W t+1 ) ≤ F (W t ) + γ t Ŝ -W t , ∇F (W t ) + (γ t ) 2 2 C F .
We can now take conditional expectation on both sides and use (7) to get

E[F (W t+1 )|W t ] ≤ F (W t ) + γ t E[ Ŝ|W t ] -W t , ∇F (W t ) + (γ t ) 2 2 C F ≤ F (W t ) + γ t min S∈D S -W t , ∇F (W t ) + (γ t ) 2 2 C F (1 + δ) ≤ F (W t ) -γ t G(W t ) + (γ t ) 2 C,
where we denote G(W ) := max S∈D W -S, ∇F (W ) and C := C F 2 (1 + δ). The function G(W ) is known as the duality gap and satisfies

F (W ) -F (W * ) ≤ G(W ) -see Jaggi (2013) for details. Denoting H(W ) := F (W ) -F (W * ), we have E[H(W t+1 )|W t ] ≤ H(W t ) -γ t G(W t ) + (γ t ) 2 C ≤ H(W t ) -γ t H(W t ) + (γ t ) 2 C = (1 -γ t )H(W t ) + (γ t ) 2 C,
where we use the duality H(x) ≤ G(x).

We shall use induction over t to prove the sublinear convergence in expectation (6), i.e., we want to show that

E[H(W t )] ≤ 4C
t+2 , for t = 1, 2, ...

We prove this for the default step size γ t = 2 t+2 (we can easily prove the same thing for the line search variant, as the resulting iterates always achieve a lower objective than with the default step size). For t = 1, we have γ

0 = 2 0+2 = 1. For any W ∈ D, we have H(W ) ≤ C F 2 < C < 4 3 C. This proves the case of t = 1. Consider now t ≥ 2, then E[H(W t+1 )] = E[E[H(W t+1 )|W t ]] ≤ (1 -γ t )E[H(W t )] + (γ t ) 2 C ≤ 1 -2 t+2 4C t+2 + 2 t+2 2 C.
Simply rearranging the terms gives

E[H(W t+1 )] ≤ 4(t+1)C (t+2) 2 < 4(t+1)C (t+1)(t+3) = 4C t+3 .
This concludes the proof.

Based on Lemma 1, in order to prove Theorem 2 we need to quantify the number of power method iterations needed to achieve the desired accuracy (7) for the linear subproblems. We will rely on some results from Kuczyński and Woźniakowski (1992, Theorem 3.1 therein), which we recall in the lemma below.

Lemma 2 [START_REF] Kuczyński | Estimating the largest eigenvalue by the power and Lanczos algorithms with a random start[END_REF] Let A ∈ R m×m be any symmetric and positive definite matrix, and b be a random vector chosen uniformly on the unit sphere (with P the corresponding probability measure). Denote by λ 1 the largest eigenvalue of A and by ξ = ξ(A, b, K) the estimate given by K power iterations. We define its average relative error e(ξ) as where π -1/2 ≤ α(m) ≤ 0.871 and, for large m, α(m) ≈ π -1/2 ≈ 0.564. Moreover, if λ has multiplicity 1, denoting the second largest eigenvalue by λ 2 , then there exists a constant K, so that for any K > K, we have

e(ξ) ≤ m λ 2 λ 1 K-1
.

We introduce a last technical lemma. According to the mean value theorem, there exists a matrix X between W and W + ∆W such that sup ∆W =0 ∇F (X), ∆W ∆W * ≤ L.

Denote the largest singular value of W by σ 1 (W ). Since the spectral norm is the dual norm of the trace norm, we have σ 1 (∇F (X)) ≤ L. Letting ∆W → 0, we get σ 1 (∇F (W )) ≤ L.

Based on the above intermediary results, we can now prove Theorem 2. For any t ≥ 0, denote A t := ∇F (W t ). The largest eigenvalue of A t A t is the square of the largest singular value of A t , denoted as σ t 1 . We estimate (σ t 1 ) 2 as v K(t) A t A t v K(t) , where v K(t) is the normalized unit vector after K(t) power iterations. We also denote u K(t) := A t v K(t) / A t v K(t) .

According to Lemma 2, we have

E v K(t) A t A t v K(t) -(σ t 1 ) 2 (σ t 1 ) 2 ≤ ln m K(t) -1 .
Therefore:

E A t v K(t) σ t 1 -1 ≤ E A t v K(t) σ t 1 -1 A t v K(t) σ t 1 + 1 = E A t v K(t) 2 (σ t 1 ) 2 -1 ≤ ln m K(t) -1 .
Let K(t) = 1 + µL(t+2) ln m δC F

, we get

E A t v K(t) σ t 1 -1 ≤ 1 2 δγ t C F µL ≤ 1 2 δγ t C F µσ t 1 ,
where the last inequality uses Lemma 3.

Removing the absolute sign and the denominator, we get

E[µ(σ t 1 -A t v K(t) )] ≤ 1 2 δγ t C F .
Rearranging the terms, we obtain

-µE A t v K(t) ≤ -µσ t 1 + 1 2 δγ t C F . (8)
On the other hand, we have

E A t v K(t) = E v K(t) A t A t v K(t) A t v K(t) = E[u K(t) A t v K(t) ] = E u K(t) v K(t) , A t , (9)
Replacing ( 9) and ( 10) into (8), we obtain (7). The first assertion of Theorem 2 thus holds by application of Lemma 1. For the second assertion, the proof is nearly identical. Indeed, by replacing ln m K(t)-1 with mβ 2K(t)-2 , we get the desired result.

Appendix B Implementation Details for Two Tasks

For the two tasks studied in this paper, we describe the sufficient information maintained by workers and how to efficiently update it. Table 2 summarizes the per-worker time and memory complexity of DFW-Trace depending on the representation used for the sufficient information. Generally, the low-rank representation is more efficient when the number of local data points n j < min(d, m).

Multi-task least square regression. Recalling the multi-task regression formulation in (4), for any worker j we will denote by X j the n j × d matrix representing the feature representation of the data points held by j. Similarly, we use Y j to denote the n j × m response matrix associated with these data points. The gradient of ( 4) is given by ∇F (W ) = X (XW -Y ). At each step t, each worker j will store (X j Y j , X j X j , X j X j W t , W t , ∇F j (W t )) as sufficient information. The quantities X j Y j and X j X j are fixed and precomputed. Given W t , W t+1 = (1-γ t )W t +γ t S t

Fig. 1

 1 Fig.1Results for multi-task least square regression. Left: low-dimensional dataset (n = 10 5 , d = 300 and m = 300). Right: higher-dimensional dataset (n = 10 5 , d = 1, 000 and m = 1, 000).

Fig. 2

 2 Fig. 2 Results for multinomial logistic regression (synthetic data) for several values of µ. Left: µ = 10. Middle: µ = 50. Right: µ = 100. The error stands for the top-5 misclassification rate.

Fig. 3 Fig. 4

 34 Fig.3Results for multinomial logistic regression (ImageNet dataset). The error stands for the top-5 misclassification rate.

  Then for any K ≥ 2 and m ≥ 8, regardless of A, we havee(ξ) ≤ α(m) ln m K -1 ,

Lemma 3

 3 If a differentiable function F is L-Lipschitz continuous w.r.t. the trace norm, then for any matrix W , all singular values of ∇F (W ) are smaller than L. Proof For any matrix W , the definition of L-Lipschitzness implies that sup ∆W =0 |F (W + ∆W ) -F (W )| ∆W * ≤ L.

Table 1

 1 Communication cost per epoch of the various algorithms. K(t) is the number of power iterations used by DFW-Trace.

	)	update

This constant is bounded above by L diam(D) 2 , where L is the Lipschitz constant of the gradient of F (see[START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF].

This can be done without communication: for instance, the workers can agree on a common random seed before running the algorithm.

In Apache Spark, this is implemented in treeReduce and treeAggregate.

https://github.com/WenjieZ/distributed-frank-wolfe

The relative performance of the methods is the same for other values of µ. We omit these detailed results due to the lack of space.

Acknowledgements This work was partially supported by ANR Pamela (grant ANR-16-CE23-0016-01) and by a grant from CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020. The first author would like to thank Ludovic Denoyer, Hubert Naacke, Mohamed-Amine Baazizi, and the engineers of LIP6 for their help during the deployment of the cluster.

Appendix A Proof of Theorem 2

Notice that our distributed version of the power method used in DFW-Trace (Algorithm 2, lines 5-10) exactly corresponds to the serial power method applied to the full gradient ∇F (W t ). Hence DFW-Trace performs the same steps as a centralized Frank-Wolfe algorithm that would use the power method to approximately solve the subproblems. We will thus abstract away the details related to the distributed setting (e.g., how the data is split, how parallel computation is organized): our analysis consists in characterizing the approximation error incurred by the power method and showing that this error is small enough to ensure that the Frank-Wolfe algorithm converges in expectation.

We start by establishing that if the linear subproblem is approximately solved in expectation (to sufficient accuracy), then the standard Frank-Wolfe algorithm converges sublinearly in expectation (in the sense of Definition 1).

is efficiently obtained by rescaling W t and adding the rank-1 matrix γ t S t . A similar update scheme is used for X j X j W t . Assuming W 0 is initialized to the zero matrix, the local gradient is initialized as ∇F j (W 0 ) = -X j Y j and can be efficiently updated using the following formula:

The same idea can be applied to perform line search, as the optimal step size at any step t is given by the following closed-form formula:

Multinomial logistic regression. We now turn to the multi-class classification problem (5). As above, for a worker j we denote by X j its local n j × d feature matrix and by Y j ∈ R n j the associated labels. The gradient of ( 5) is given by ∇F (W ) = X (P -H), where P and H are n × m matrices whose entries (i, l)

are P il = exp(w T l x i ) k exp(w T k x i ) and H il = I[y i = l] respectively. The sufficient information stored by worker j at each step t is (X j , X j H j , X j W t , ∇F j (W t )). X j H j is fixed and precomputed. Assuming that W 0 is the zero matrix, X j W t is initialized to zero and easily updated through a low-rank update. The local gradient ∇F j (W t ) = X j P j -X j H j can then be obtained by applying the softmax operator on X j W t . Note that there is no closed-form for the line search.