Andreas Ar- Gyriou

Nicolas Baskiotis

Rodroigo Fernandez

Stéphane Gaïffas

Vianney Perchet

First of all I have to thank my uncle Mehrdad Shahshahani who contaminated me with mathematics. Then I would like to thank my advisors Theodoros Evgeniou and Nicolas Vayatis for their patience, insights, and for the profusion of their rich ideas which

Abstract

Predicting connections among objects, based either on a noisy observation or on a sequence of observations, is a problem of interest for numerous applications such as recommender systems for e-commerce and social networks, and also in system biology, for inferring interaction patterns among proteins. This work presents formulations of the graph prediction problem, in both dynamic and static scenarios, as regularization problems. In the static scenario we encode the mixture of two different kinds of structural assumptions in a convex penalty involving the ℓ 1 and the trace norm. In the dynamic setting we assume that certain graph features, such as the node degree, follow a vector autoregressive model and we propose to use this information to improve the accuracy of prediction. The solutions of the optimization problems are studied both from an algorithmic and statistical point of view. Empirical evidences on synthetic and real data are presented showing the benefit of using the suggested methods.

Overview and contributions

Industrial context and modeling choices

The current work is the result of close collaboration between the web-advertisement and marketing company 1000mercis and the the CMLA academic research laboratory at Ecole Normale Supérieure de Cachan. The collaboration was done in the context of a CIFRE (Conventions Industrielles de Formation par la REcherche) contract between the company and the laboratory. The research project was motivated by data mining challenges 1000mercis had faced. We provide algorithmic tools to solve some of these problems that are frequently met in the context of internet applications and also in other fields.

Predictive algorithms for graph data are needed for the conception of recommender systems, the valuation of network agents, prediction of trends and discovery of hidden patterns in complex real-world networks arising from various application domains such as internet data and gene-expression networks in biology. The purpose of this thesis is to introduce a methodology and related practical implementation techniques for exploiting the temporal and static patterns of complex networks in a predictive perspective. The need for robust and scalable recommender systems and valuation methodology for relational customer relationship management (CRM) data sets such as those arising from purchase histories in customer-to-customer (C-to-C) contexts (see [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF]), or social networks has motivated and defined the scope of the research project. The contribution of this work, in addition to the descriptive presentation and the analysis of the data provided by the company, is to formulate the problem of prediction in graph sequences as a regularization problem and explore the properties of the solutions both empirically and theoretically. In fact even though standard link prediction heuristics may perform well over some databases, the need for more accurate systems, robust measurements and functionality guarantees reveal the importance of developing more powerful tools. We chose to formulate the problem of prediction in graphs with the vocabulary of regularization methods. These methods have received a lot of interest from signal processing, imaging, computer vision, machine learning and high-dimensional statistics communities during the recent years. We can summarize the main motivation of the modeling choices made throughout this work as follows:

• Complexity measure. The crucial ingredient for formulating the prediction problem in graphs as a regularization problem is to define adequate complexity measures for graphs.

A very basic and insightful observation is that link prediction heuristics can be expressed as spectral functions of the adjacency matrix that flatten the spectrum of it (see Figure 1.5, page 35). Hence spectral measures such as the rank or a convex surrogate, as the nuclear norm, seem to be promising candidates. These spectral measures are closely related to the community structures of the graphs. In fact block-diagonal and overlapping blockdiagonal matrices have low rank and they belong to sets of adjacency matrices representing graphs containing sharp clusters or overlapping clusters. The low-rank prior and its convex surrogates [START_REF] Jameson | Summing and Nuclear Norms in Banach Space Theory[END_REF] (the trace norm [START_REF] Bach | Consistency of trace norm minimization[END_REF], the max norm [LRS + 10] for instance) have been studied by Srebro [Sre04] in the collaborative filtering context. We adapt them to the framework of highly clustered graph data where low-rank and sparsity effects are mixed.

• Simultaneous prediction of local and global effects. Another important direction of investigation is the multi-scale structure of the data: as in many complex, self-organized systems, global and local effects highly interact. In case of dynamic graphs, the idea is that global effects such as life cycles, trends or communities and local connections among nodes of the networks are inseparable phenomena. The estimation method should therefore take into account the nature of mutual interactions among these objects.

• Convex optimization tools. It is a practical choice for guaranteeing the existence and possibly the uniqueness of the solution of the formulated problem to use a convex optimization framework. Such a framework allows in many cases to design algorithms converging to the true or approximate solution in polynomial time, even if the initial problem of interest is non-polynomial. In fact, we know from compressed sensing theory [START_REF] Candès | Decoding by linear programming[END_REF] that, in some regime, convex relaxation gives the same solution as the corresponding NP hard problem, leading to no loss of information. The very rich literature on optimization provides algorithms for a large body of convex functions that include our functionals. We will provide methodological tools for implementing efficient algorithms, for tuning the parameters, and will provide empirical evidence of the algorithms performance. In matrix factorization problems, the convex formulation leads to computationally intensive iterative algorithms because of the computation of a singular value decomposition at each iteration step. We will argue in a specific matrix factorization task how to replace the convex formulation with a partially convex but more scalable one, and will provide empirical evidence of success of the obtained method. Theoretical guarantees of convergence for these algorithms are missing and constitute a future direction of research.

• Generative models. We also introduce generative models of temporal and static graphs that mimic significant properties of the real-world data, and illustrate cases where all the assumptions required for the models to fit the data hold true. After arguing their interesting properties and motivating the modeling choices, we provide theoretical analysis of the predictive models over the synthetic generative graphs. In particular we prove oracle bounds for the regression problem subject to a mixed ℓ 1 -trace norm penalty we designed and deduce related bounds for the prediction problem in an autoregressive framework.

In terms of methodology, introducing such models allows one to study the impact of parameters on the model separately and to understand their interactions.

Contributions

The main contribution of this work is the formulations of the graph prediction problem in the static and dynamic scenarios as convex optimization problems. We provide theoretical and empirical evidence in support of this approach.

1. Estimation of sparse low-rank matrices and graph denoising. Previous works in highdimensional statistics had tackled the problem of estimating a vector having a few nonzero elements. In the same spirit, the estimation of matrices that have low rank has been investigated in an undersampled setting. It is known that the estimation is inaccurate in the first case when the few active variables are grouped into sets of highly correlated variables.

On the other hand, the standard matrix completion procedure fails when the coefficients of the low-rank matrix are mostly equal to zero. It turns out that in many applications the objects of interest are sparse and low-rank simultaneousely. Such an object cannot be estimated using each of the existing methodologies alone, as it has the two undesirable properties. Namely, only a few nonzero coefficients appear in grouped structures. If the groups are known in advance recent works [START_REF] Bach | Structured sparsity through convex optimization[END_REF][START_REF] Baldassarre | A general framework for structured sparsity via proximal optimization[END_REF] have been devoted to the estimation of the coefficients. In many applications the groups themselves are unknown making the task even harder. Block-diagonal matrices are examples of such objects that naturally arise from clustering or community detection applications. We develop a methodology to estimate, under several scenarios, such matrices. Denote by X a candidate adjacency matrix with coefficients X i,j . We introduce a regularizer (see [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF]) composed by the sum of two well-studied terms as X → τ X * + γ X 1 . We will argue that the mixed norm combines the benefits of each of the two norms and is suitable for estimating sparse low-rank matrices. This penalty can be written as τ i σ i (X) + γ i,j |X i,j | where σ i (X) represents the ith singular value of X or the square root of the ith largest eigenvalue of X ⊤ X. The set of sparse low-rank matrices obtained by minimizing objectives penalized by such a term turns out to contain matrices that are -up to permutation of rows and columns-block-diagonal and overlapping block-diagonal forms. These matrices can be interpreted as adjacency matrices of networks containing highly connected groups of nodes and therefore are of particular interest for prediction and denoising applications in graph data, and in covariance estimation. Assume we are given a loss function ℓ(•, Y ) : R n×n → R that measures for instance the fit to the observed data Y . We use the latter norm in order to relax the hard optimization problem

min X∈S ℓ(X, Y )
subject to X 0 ≤ q and rank(X) ≤ r to its convex surrogate

X = arg min X∈S {ℓ(X, Y ) + γ X 1 + τ X * } .
We will provide algorithms to minimize the objective over the convex set of admissible solutions S for different loss functions ℓ and also provide oracle bounds for the regression loss. In particular we study the problem that we call graph denoising where we consider an ℓ 0 loss ℓ(X, Y ) = X -Y 0 . It consists of finding the fewest edges to flip (turn off if on, and on if off) making the new graph adjacency matrix low-rank and sparse.

We will provide algorithms to solve the convex surrogate of this problem that consists of minimizing the convex nonsmooth objective X → X -A 1 + τ X * + γ X 1 .

2. Prediction of graphs with autoregressive features. We also study the graph prediction problem in a dynamic setting [START_REF] Richard | Link discovery using graph feature tracking[END_REF][START_REF] Richard | Link prediction in graphs with autoregressive features[END_REF]. In this case we need to assume that a set of graph descriptors such as degrees or any linear function of the graph evolve smoothly over time. The objective function we design for this problem has two variables: W and A. The variable W is the VAR parameters set that handles the prediction of the features by minimizing the least squares fit to past data term W → X T -X T -1 W 2 F where T is the number of observed snapshots and X T , X T -1 represent collections of observed features from time t = 1, • • • T and t = 0, • • • T -1 respectively. The second variable A represents the adjacency matrix of the graph to be predicted. We argue that the prediction of A can be done through the prediction of the features of A that we compute through a mapping ω : R n×n → R d . We refer to this mapping as the feature map. For link prediction heuristics, the stability of the distribution of some features explains the success of the methods. We try to generalize this mechanism by designing a methodology where the slow evolution of a set of well-chosen features makes their prediction easier than the prediction of each individual edge directly. Then estimating the graph given a linear observation of it is very similar to matrix completion, or regression with spectral penalty. We use the mixed ℓ 1 -trace norm to this end. The objective can be written in the following

A 0 , A 1 , • • • A T A T +1 Adjacency matrices observed ∈ R n×n ↓ ω ↓ ω ↓ ω ↑ ω(A 0 ) , ω(A 1 ) , • • • ω(A T ) → W ω(A T +1 ) Feature vector representations ∈ R d
Table 1: General scheme of the method suggested for prediction in dynamic graph sequences through a feature map ω.

form:

L(A, W ) . = 1 dT X T -X T -1 W 2 F + κ W 1 + 1 d ω(A) -W ⊤ ω(A T ) 2 2 + τ A * + γ A 1
This objective function contains a smoothly differentiable convex least squares term interpreted as a loss (A, W )

→ 1 dT X T -X T -1 W 2 F + 1 d ω(A) -W ⊤ ω(A T ) 2 2
, and a penalty term acting on the joint variable (A, W ). The penalty term acts differently on the A subset and W subset of variables. The penalty on W is a basic element-wise ℓ 1 term that encourages sparsity of the VAR matrix W . This penalty can be replaced by an ℓ 2 term if the sparsity assumption is not relevant for the domain. On A we incorporate a mixed ℓ 1 plus trace norm penalty encouraging simultaneously sparsity and low rank of A. It is motivated by the block-diagonal and overlapping block-diagonal matrices. Table 1 summarizes the methodology in a scheme where the symbols ↓ ω represent the feature extraction procedure through the map ω : R n×n → R d . The prediction in the feature space is represented by → W , and is handled in practice by the least squares regression on W , and finally the symbol ↑ mapping the predicted feature vector ω(A T +1 ) to A T +1 represents the inverse problem that is to be solved through the regression penalized by the mixed penalty. It is important to keep in mind that the mapping ω that is relevant for real world applications may be nonlinear but in the current work we only study linear approximations to such maps. For instance an interesting set of features of graphs is obtained by counting the number of length k (k = 3, 4, • • • ) cycles passing through each node, these numbers are obtained algebraically by taking the diagonal elements of A k where A represents the adjacency matrix of the graph. These mappings are obviously nonlinear. We leave their studies for future work.

3. Algorithmic and empirical evaluation of the models. The algorithms introduced for solving the convex but nonsmooth optimization problems defined in both static and dynamic contexts arise from the family of forward-backward algorithms [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] and their generalizations. Specifically, minimizing an objective function containing a smoothly differentiable convex loss plus the mixed penalty τ X * + γ X 1 can be done through the generalized forward-backward algorithm [START_REF] Raguet | Generalized forward-backward splitting[END_REF], and the same algorithmic scheme can be used in the dynamical context. A Douglas-Rachford scheme using a generalization of the soft-thresholding operator to piecewise affine functions handles the case where the loss function is the graph denoising loss ℓ : X → X -A 1 . The need for scalability and the high computational cost of the shrinkage operator through a soft-thresholded singular value decomposition suggests the use of the variational characterization of the trace norm X * = 1 2 min U V ⊤ =X U 2 F + V 2 F in the objective function rather that its closedform expression. Hence we also wrote alternate minimization algorithms for optimizing surrogate objectives. These algorithms are known to perform well in practice [START_REF] Koren | Factorization meets the neighborhood: a multifaceted collaborative filtering model[END_REF], to be scalable and implementable in a parallel architecture [START_REF] Recht | Parallel stochastic gradient algorithms for large-scale matrix completion[END_REF] but the lack of convexity makes their analysis more challenging. We used the MATLAB implementations of the algorithms for evaluating the suggested models on both real and synthetic data and report empirical results illustrating the comparative superiority of the suggested methods relative to the standard ones and also describing regimes of noise and regularity where accurate predictions can be expected as phase transition diagrams.

4. Theoretical guarantees. The formulation of the problems as convex optimization procedures allows one to benefit from the large body of work on the analysis of these classes of problems. In particular the oracle bounds that we were eager to obtain have been derived by using the tools introduced for studying the Lasso (ℓ 1 penalized least squares regression) and the trace-norm penalized regression. In our static setting we consider a low dimensional linear measurement ω(A 0 ) of the target matrix A 0 to be observed after corruption by noise: y = ω(A 0 ) + ǫ. Let A denote the estimator obtained by minimizing the convex surrogate objective. In this setting these inequalities bound the prediction error of the estimator ω( A)ω(A 0 ) 2 by the the infimum over the set of admissible solutions of the squared error plus a complexity term involving both the rank and sparsity index of the candidate. The inequality shows in particular that optimizing the convex surrogate leads us to the best solution, if we were able to afford the high cost of minimizing the NP-hard real problem containing a rank and a sparsity index constraint. Technically, instead of directly decomposing the noise into the estimator support and the orthogonal, we first split the noise into two parts and perform the projections on each part separately.

The splitting constant α ∈ (0, 1) appears in the propositions' assumptions as a trade-off constant between the two tuning parameters τ and γ. Using this mechanism, all the inequalities obtained can be interpreted as smooth interpolation between the corresponding inequalities for the Lasso and for the trace-norm penalized matrix completion, and shows that our estimation procedure performs at least as well as the best of the two. We also introduce the restricted isometry constant, the cone of restriction and restricted eigenvalue constants that extend the existing definitions to our framework.

In a nutshell, the oracle bounds we provide for the time-dependent setup hold the following message: the prediction error and the estimation error can be simultaneously bounded by the sum of three terms that involve homogeneously (a) the sparsity, (b) the rank of the adjacency matrix A T +1 , and (c) the sparsity of the VAR model matrix W 0 . The tight bounds we obtain are similar to the bounds of the Lasso and are upper bounded by:

C 1 log d T d 2 W 0 0 + C 2 log n d A T +1 0 + C 3 log n d rank A T +1 .
The positive constants C 1 , C 2 , C 3 are proportional to the noise level σ. The interplay between the rank and sparsity constraints on A T +1 are reflected in the observation that the values of C 2 and C 3 can be changed as long as their sum remains constant. The bounds obtained for regression with the mixed penalty are similar to the latter results for C 1 = 0.

Organization

This dissertation is organized as follows. The Chapter 1 is a general introduction. In Chapter 2 we introduce tools to explore dynamic graph data sets. We provide definitions of some topological features of graphs that allow one to measure quantities of interest and illustrate their distributions over several real data sets, discussing implications on various properties of the data we can read from the plots and data. We point out the connections between specific distribution of graph features and heuristics approaches used for predicting links in networks. We also review standard random graph generative models before introducing some random graph models that reproduce the properties of real-world graphs we are interested in. In particular, these models allow us to test the performance of our algorithm in a framework where the working hypotheses are rigorously fulfilled. This is a standard way of validating an estimation methodology before testing the validity of the working assumptions. In Chapter 3 we introduce a new regularization tool for estimating sparse low-rank matrices. This regularizer is a norm formed by taking the positive linear combination of two norms, each being used separately in previous applications. The ℓ 1 norm has been used for estimating sparse vectors independently from the dependence structure among variables, and the trace norm, or the sum of singular values, is used for estimating low-rank matrices.

Finally in Chapter 4 we formulate the problem of prediction in a sequence of graphs as a regression problem with a mixed penalty term. We provide optimization algorithms for joint prediction and VAR estimation and prove theoretical bounds guaranteeing the quality of the estimator under natural assumptions on the noise process. We point out that several simplifications had to be made in the models we study in order to formulate it using the tools of convex optimization. This formulation has the benefit of relating the work to a large body of theoretical and numerical analysis. The current work should be seen as a first attempt to develop prediction models for dynamic graphs using convex estimation tools.

Chapter 1

Introduction

Motivations

Prediction for anticipation, discovery and serendipity

The methods and algorithms studied in the current thesis are devoted to explore some of the strategic research and development axes identified by the company 1000mercis that is our industrial partner. 1000mercis, a pioneer in interactive marketing and advertisement founded in Paris in 2000, provides acquisition and retention tools through different channels: emailing, mobile and viral games. The main challenges motivating our collaboration with 1000mercis can be summarized as follows.

1. Reducing promotional pressure. Emailing has no cost, therefore the immediate danger many emailing service providers face is to overwhelm consumers. The mission of a responsible company is to struggle against such abusive practices by trading the amount of advertisements sent with their quality. The obvious benefit of such a sustainable strategy is to avoid long-term fatigue of users and reducing churn rate, i.e. the ratio of users who unsubscribe from a promotional campaign.

2. Serendipity. Web 2.0 places the customer at the heart of its development. By providing socio-demographic as well as browsing and purchase data of users to the websites, many potential services involving personalized advertisement have been made possible. The recommender systems, that are one of the main application of the current work, aim precisely at using past users data in order to provide the users information about products they may like and do not know about.

3. Optimize immediate efficiency. The last but not the least point is the immediate benefit generated by using powerful tools of modern statistics. Some preliminary tests in 1000mercis had shown that the efficiency of advertisement campains can be improved by a factor of 5 by tuning the parameters of a basic algorithm.

The purpose of storing exponentially growing amounts of data is to use the past data in order to make better decisions. We argue that the key to the three before-mentionned challenges is to predict several types of phenomena by forecasting users activity levels and trends in dynamic relational databases.

• Forecasting users activity levels and trends one can adjust promotional pressure at individual and global levels, anticipate important variations in the supply-chain, and push fragile tendencies

• Predicting users preferences filters information for adapting the content of advertisement to individuals tastes and also to share similar users preferences making one discover new items.

• Fitting predictive models to the market data allows to build valuation tools at micro and macro levels, and hence attribute the promotional assets optimally.

In statistics and machine learning the graphs have classically been used for modeling dependence among variables. Numerous examples can be given from biology applications where the huge number of variables imposes to find sparse dependency patterns among them. The new type of problems involving these objects uses them in a very different fashion. In fact the growth of internet has produced considerable amounts of relational data (social networks, purchase histories, communication networks, hyperlinked websites) where the observations are modeled directly as graphs. Graphs are in this case not anymore artifacts used for better visualizing or understanding dependence structures among variables, but rather constitute the observation and the output variable of interest.

Growing graphs of the world wide web

The data produced nowadays by humans every year overpass by an order of magnitude the amount of data produced up to 2005. The cost of storing 1Tera-Byte of data in 1980 was $14M and it is not more than $30 in 2010. What The Economist calls the 'The Data Deluge' [The10] opens the door to a wide spectrum of blue oceans of scientific research and economic development. These new range of data, generated by various types of applications including governments records, scientific data, internet websites, raise many original challenges for statisticians, computer scientists and mathematicians. An important portion of the data generated by internet can be viewed as network data evolving over time. In fact they contain pairs of connected objects such as websites connected by hyperlinks, users linked to websites through theirs surfing and clicks, users associated to products through purchases or users tied through social networks. The link structure among these entities, and its evolution dynamics contains highly relevant information. The value of this information can be highlighted by noticing the success of those who were able to use it in a proper way. For instance • Google's pagerank [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF] that is the basic ingredient of the search engine aims at quantifying the relevance of each website given the hyperlink structure of its surrounding.

• Amazon was the first to use purchase history as a mirror for predicting future purchases. It has used its valuable database for earning up to 35%of its turnover through the recommender engine [START_REF] Linden | com recommendations : Item-to-item collaborative filtering[END_REF].

• Facebook Open Graph is a protocol1 that aims at superimposing various internet graphs (mainly user-products and user-users) in order to use the social link structure for more efficient advertisement targeting.

The emergence of internet as a new communication medium has fundamentally changed the nature of advertisement by placing the user at the heart of the communication chain. For the advertiser, besides making communication almost cost-free and instantaneous, two major revolutionary assets of internet as a communication medium are the measurability of the users feedbacks to ads and the availability of inter-users communication traces. In fact, in comparison with the standard offline Business-to-Customer, the clicks or purchases initiated by an advertisement can be interpreted as the success of the ad or the promoted product and are a golden key for the marketers in their product design. On the other hand the growing popularity of user generated websites and social networks has made the management of viral marketing2 possible. In fact the former non-measurable word-of-mouth advertisement has turned out to be observable and measurable at individual and aggregated levels, and given its tremendous efficiency, marketers run after tools allowing them to master related technologies. Note that despite being rich in content and widely available, the data collected from internet are not always smoothly useable raw materials for advertisers. In fact most declarative data as questionnaires, turn out to contain many missing and noisy values, and the spontaneous surfing and behavior data is high-dimensional, noisy and incomplete. A common point among many types of data collected from internet, such as purchase history, surfing data, communication, collaboration and social media data, are that they can be modeled as connections among users, or connections of users and items (webpages or products), or connections among objects. The technical application of the thesis is to suggest methodologies to handle these types of data for predictive goals.

The complexity and the novelty of graph data make their statistical analysis, engineering and retrieval non-trivial. As graph data can not be easily embedded onto a usual vector space, the standard tools of statistics and machine learning do not apply to them when one needs to achieve clustering, classification or prediction. The goal of the current thesis is to introduce a set of mathematical tools enabling the formulation of the problem of prediction in dynamic graphs such as those produced massively by web applications as a regularization problem, and suggest methods to solve them efficiently.

Graphs of multivariate data

When observations are modeled as vector data, graphs are useful objects for describing highdimensional distributions. In fact when the number of variables grows, not only the study of each individual variable is of interest, but in many cases the interaction structure among different variables turns out to be both informative and interesting by its own. Forecasting the behavior of systems with multiple responses has been a challenging problem in the context of many applications such as collaborative filtering, financial markets, or bioinformatics, where responses may be, respectively, movie ratings, stock prices, or activity of genes within a cell. Statistical modeling techniques have been widely applied for learning multivariate time series either in the multiple linear regression setting [START_REF] Breiman | Predicting multivariate responses in multiple linear regression[END_REF] or with autoregressive models [START_REF] Tsay | Analysis of Financial Time Series[END_REF]. More recently, kernel-based regularized methods have been developed for multitask learning [START_REF] Th | Learning multiple tasks with kernel methods[END_REF][START_REF] Argyriou | A spectral regularization framework for multi-task structure learning[END_REF][START_REF] Lounici | Oracle inequalities and optimal inference under group sparsity[END_REF]. These approaches share in common the use of the interactions between input variables to enhance prediction of every single output.

In some applications the interactions among variables constitute the main object of interest. For instance, a problem of major interest in functional genomics, systems biology and drug discovery is the estimation of regulatory networks in gene expression, and interaction networks among proteins that could explain genetic variability in cancer and other diseases (see [START_REF] Barabasi | Network biology: Understanding the cell's functional organization[END_REF][START_REF] Hue | Large-scale prediction of protein-protein interactions from structures[END_REF][START_REF] Jacob | Protein-ligand interaction prediction: an improved chemogenomics approach[END_REF] for explanations and the figure 1.1.3 for illustration3 ). Nodes represent genes or proteins and edges interactions. The studies in these fields aim at discovering the mechanisms leading to diseases such as the development of cancer cells. Understanding the functions of the central nervous system and the brain [VLEM01] also involves a large number of objects interacting together. Based on the observation of gene expressions, protein interactions, or cell activities statistical tools are designed to infer the correlation mechanisms, and at a higher level, measure causality among different expressions. For this goal methodologies involving Bayesian Networks [START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF] and causality networks have been developed, both in static and dynamic [START_REF] Song | Time-varying dynamic bayesian networks[END_REF] cases. Bayesian Networks model the random variables as nodes of a graph that are connected through an edge if they are dependent. If no cycle appears in the graph of variables, then the joint density function of the d variables

x 1 , • • • , x d can be written as p(x 1 , • • • , x d ) = (a,b)∈E p(a|b)
where a and b are pairs of variables and E represents the edge set that contains all the pairwise connections. The undirected Bayesian Networks are also known as Markov Random Fields.

They have the power of modeling cyclic dependency behavior and have been used in computer vision and image processing [START_REF] Komodakis | Mrf energy minimization and beyond via dual decomposition[END_REF]. Even though causality remains a very delicate notion to measure, concise definitions have attempted to recover this notion. Granger causality for instance is defined for a pair of temporal random variable (X, Y ) as follows. If the prediction of Y given past values of both X and Y outperforms the prediction of Y given only its own past values in an autoregressive way, then X is said to be Granger-causal for Y . Networks of Granger Causality describe the paths and mechanisms leading to specific expressions in living metabolisms [START_REF] Shojaie | Discovering graphical granger causality using a truncating lasso penalty[END_REF].

For studying interactions among different variables a common practice is to introduce matrices that have entries quantifying pairwise interactions. Two matrices are of particular interest in this context: the covariance matrix itself and its inverse that is also called the precision matrix. The correlation matrix can be seen as a dimensionless covariance matrix. The precision matrix entries are interpreted as partial covariance between two variables given all the others are set to fixed values. In case the data is i.i.d. multivariate gaussian, it turns out that the sparsity pattern of the precision matrix determines the pairwise dependence structure among random variables. Even though estimating the empirical covariance matrix of a sample is straightforward, when the number of variables is larger than the sample size or in presence of noise, the empirical sample covariance is a poor estimation of the true covariance matrix, and additional prior knowledge is needed to estimate the true covariance or precision matrix. Efficient procedures developed in the context of sparse model estimation mostly rely on the use of ℓ 1 -norm regularization [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF]. Natural extensions include cases where subsets of related variables are known to be active simultaneously [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. These methods are readily adapted to matrix valued data and have been applied to covariance estimation [START_REF] Karoui | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF][START_REF] Bien | Sparse estimation of a covariance matrix[END_REF] and graphical model structure learning through the estimation of the inverse covariance by using the sample covariance as the input of an estimation procedure [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF].

Context

Users behavior modeling with redundant and missing data

Modeling users behavior is of significant interest for marketers. The use of electronic devices such as internet as advertisement platforms has drastically changed the marketers job by suggesting new tools to measure users sensitivities to different stimuli. In fact, as opposite to traditional advertisement campaigns, the feedbacks of the electronic media users to advertisement campaigns can now be stored at individual level and give access to measurements of The key ingredient for designing individual level ads is to have access to users preference data, which are often incomplete or noisy. A common remedy to the missing data problem [START_REF] Marlin | Missing Data Problems in Machine Learning[END_REF] is to use the redundancies in users behaviors in order to infer the missing features. This controvertial but still rich method has a crucial benefit: it exploits the redundancy of the users behavior for considerably reducing the search space, or equivalently, multiplying data points. In mathematical terms, the manifold on which the users descriptor parameters lie has a much lower dimension than the number of apparent characteristics and the number of users. This implies that by using such prior knowledge on the domain one can infer the missing entities by detecting similarities among users and hence boost the quality of the database artificially. Techniques of the same spirit have extensively been used in multivariate regression [START_REF] Breiman | Predicting multivariate responses in multiple linear regression[END_REF], multitask learning [START_REF] Argyriou | Learning to Integrate Data from Different Sources and Tasks[END_REF] and matrix completion [START_REF] Srebro | Maximum-margin matrix factorization[END_REF].

Estimation for valuation

Customer valuation consists of scoring the clients and prospects by using their past actions for estimating their expected spendings, and allocate accordingly the promotional assets. In fact, given the unbalanced power-law distribution of consumer spendings, or the Pareto 80-20 rule that states that 80 % of the revenue is generated by 20% of the customers, it is crucial to detect or predict the small portion of the population that generates the majority of the company's income for privileging them. Standard tools evaluate the customer value by simply averaging customers past activities or tracking the corresponding curves. In many modern data sets, the network structure makes the prediction of many actions much harder due to interactions among the users that strongly impact their behavior. For highlighting the relevance of this application of our work, we show how to valuate buyers and sellers in a C-to-C (costumer-to-costumer) market seen as a bipartite network, and refer to the paper [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF] for further details on the methodology. Here is a brief description of the data and the major elements of the estimation tools used in the paper.

Seller and buyers valuation in a C-to-C book market based on sales time series. We used approximately 6 years historic data from an online C-to-C bookstore. The market is formed by 1.3M buyers, 0.3M sellers, 1.5M products tracked over 322 weeks, representing an overall 8.6M transactions, having generated 80 M euros turnover (i. We did not take the December data into account due to the particular purchase behavior of users that should be addressed independently from the rest of the year. We also did not include the first week of january as it may have less than 7 days.

Following standard econometric methodology we suggested in [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF], we estimate using an ordinary least squares regression the matrix autoregressive coefficients of the multivariate time series formed by the number of new buyers, returning buyers, new sellers, returning sellers and the turnover. The obtained matrix is then used for computing the cumulative turnover, after a period of time, resulting from shocks on each actor (new/returning, buyer/seller). The interpretation is that the unit actor at time t = 0 is seen as a perturbation to the system, and one measures the cumulative value generated by a small perturbation on each of the actors volumes. This gives an estimate of the value of each of the users of the C-to-C market, and the cumulative impulse response of market users are interpreted as valuations of the network effect. We refer to [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF] and references therein for details on the methodology, and to 1.2.2 for the cumulative impulse response curves obtained in our experiments, and see table 1.1 for the table of results obtained in [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF]. We highlight that the purpose of the current thesis is to develop an estimation methodology taking into account both time series and network data in order to capture the effects of each on the other thanks to a joint estimation.

Revenue

Network effect: Network effect: For classical data sets available for econometric measurements only the time series of some variable are available and therefore vector autoregressive (VAR) models are the most appropriate tools. In the case of e-commerce data sets, individual level data are available, and therefore the estimation method can take them into account in order to get a more accurate estimation. Using the vocabulary of dynamic graph sequences we formulate the problem of consumer and contributor valuation for the C-to-C data presented below. Consider a C-to-C market formed by Sellers (S), Products (P) and Buyers (B). We represent the market by a dynamic tripartite graph. The three sets of nodes are S,P and B. A transaction is modeled by a pair of links between a buyer and a product, and a seller to the latter product.

Let n B , n P , n S design the number of buyers, products and seller, and n = n B + n P + n S be the overal number of actors in the market.

We can represent the transactions at each date by 3 incidence matrices B P of size n B × n P , B S of size n B × n S and P S of size n P × n S . The (i, j) entry of a matrix X Y is non-zero if a We can denote by ω(A t ) ∈ R d the vector containing the set of all such scalar variables of interest. Given the observation of ω(A t ), we aim at predicting ω(A t+1 ). For this end we use the linear autoregressive process

ω(A t+1 ) = f * (ω(A t )) + ǫ t+1 (1.1)
where ǫ t+1 represents noise, that must be assumed to be i.i.d. and to have zero mean. The function f , taken linear at the first glance, can be estimated not only based on the set of vectors (ω(A t )) t , but also using the network information, namely the adjacency matrix sequence.

To summarize, in the setting of our problem we consider

• Adjacency matrices A t ∈ R n×n representing a graph sequence for t ∈ {1, 2, ..., T }

• A feature map ω : R n×n → R d . We refer to ω(A t ) as the feature vector of A t .

and we want to find

• A matrix A ∈ R n×n representing a link discovery score, or the estimate of A T +1 , • A prediction function f : R d → R d such that f (ω(A t )) ≃ ω(A t+1 )
The methodology suggested allows to not only results in more accurate estimation of the system parameter but also allows to valuate customers at individual level.

One of the contribution of the current thesis is to suggest a methodology for simultaneously estimating f and predicting A T +1 . We just discussed why the accurate knowledge of the multivariate function f is of interest for valuation. In the next section we see why the prediction of A T +1 is of interest.

Diffusion curves and network effect

The diffusion or penetration of a product refer to the different cycles a product goes across during its life, starting from the innovation, and ending by obsolescence. Studies in marketing [Bas63, Rog62, MMF90, SJT09, VFK04] have attempted to characterize the cycles of life of products. The parametric models introduced, despite being simplistic, provide interpretable schemes. Bass modeled the diffusion of innovations in networks of individual [START_REF] Bass | A dynamic model of market share and sales behavio[END_REF] by stating that "the probability of adopting by those who have not yet adopted is a linear function of those who had previously adopted". In mathematical terms, these models that aim at also including the models of epidemic propagation in networks can be stated as the following discrete time equation, involving the fraction of infected population at time t:

F (t) = F (t -1) + p(1 -F (t -1)) + q(1 -F (t -1))F (t -1)
where p is the coefficient of innovation and q the coefficient of imitation. This model results in the following differential equation

∂F ∂t = (p + qF (t))(1 -F (t)) , with F (0) = 0
which admits a closed form solution:

F (t) =
1e -(p+q)t 1 + q p e -(p+q)t It can be easily seen that when p > q, F (t) is a concave function of time, and if p < q, then the curve has an inflection point and is said to have an S-shape. In the S-curves sales volumes are said to follow over time, different phases of development, innovation, early adoption, early majority late majority and lagged sales are delimited by quantiles of a standard gaussian distribution. More recent works on the topic (e.g. [START_REF] Sood | Functional regression: A new model for predicting market penetration of new products[END_REF][START_REF] Vakratsas | The Shape of Advertising Response Functions Revisited: A Model of Dynamic Probabilistic Thresholds[END_REF]) have attempted to relax this very rigid parametric model of diffusion.

In marketing and economics, the notion of network effect is used to design the surmodular [START_REF] Topkis | Supermodularity and complementarity[END_REF] increase of value of a network of people with the population size. In fact the potential number of connections among network nodes increasing quadratically, it has been observed that the responses of a network of individuals is richer than the sum of individual response. A more data-oriented study by Leskovec and coauthors [START_REF] Leskovec | Graphs over time: Densification laws, shrinking diameters and possible explanations[END_REF] rather suggests that in many real-world graphs the evolution of the number of edges is related to the number of nodes by a densification power law:

E(t) = N (t) a
with a coefficient a ∈ [0; 1]. The network effect has been a leading concept in the study of word of mouth in marketing and recommendation [START_REF] Villanueva | The impact of marketing-induced versus word-of-mouth customer acquisition on customer equity growth[END_REF][START_REF] Leskovec | Dynamics of viral marketing[END_REF] and also C-to-C two-sided markets [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF].

Design of recommender systems

One of the motivating applications of predictive algorithms on graph data is the design of recommender systems. Predicting univariate binary outputs in user behavior, such as clicks on ads, or opening advertisement emails can be handled by standard statistical tools, and for which measurement methodologies have been studied extensively in statistics and information retrieval literature. Note that the statistical tools needed for building recommender systems can not be found in standard statistics and data-mining methods. The reason is the particular characteristics of the data. In fact a naive approach would consist of treating the recommendation problem as a multi-class classification or multi-class ranking problem for each individual taken separately, the classes being the different items. In this case the number of classes being extremely high, and the number of observations very small, the problem is highly underdetermined and the algorithm will perform very poorly. A common remedy to deal with the very small number of observation by individual is to form groups of similar individuals and to complete unobserved features of each by the observed features of the others. The rigorous definition of functional spaces that handle the task of measuring similarity among partially observed data points and completing the missing data turns out to be very challenging. The objective of the current thesis is to attempt to answer some questions risen in this direction. Developed first by Greg Liden at Amazon in early 2000's [START_REF] Linden | com recommendations : Item-to-item collaborative filtering[END_REF], recommender systems provide lists of items that users of an e-commerce website are likely to purchase, given the past collective purchase history. The philosophy behind such systems is that given the overwhelming amount of information available online, filtering the content based on the past users declaration, surfing and purchase history is a useful service. The goal being to emphasize for each user the very few items he is likely to purchase as a recommendation list rather than facing him to the whole items catalogue equally weighted. Even though personalized filtering of web information at user level sounds seducing, the use of purely predictive methodologies for such a goal has been subject to criticism [START_REF] Bodapati | Recommendation systems with purchase data[END_REF][START_REF] Fleder | Blockbuster culture's next rise or fall: The impact of recommender systems on sales diversity[END_REF][START_REF] Pariser | The filter bubble : what the Internet is hiding from you[END_REF]. A main argument that predictive methods face is that by fitting to the tastes of a user the system does not make him discover new items. In other words, if one knows that a user is likely to buy an item, there is no need to recommend him that item. As one of the major benefits of recommender systems is serendipity (discovery of unexpected good things), fitting to the users tastes is indeed undesirable. We argue that once we acquired the technical tools allowing to fit the tastes of users, making them discover new items can be handled efficiently. Therefore we emphasize that our goal in the current thesis is to explore methods that allow to fit to the evolution of dynamic graphs. Following the applications in mind, it is up to the practitioner to use the suggested methodology in the convenient way.

Related works

The problem of predicting the evolution of a graph sequence is a hard problem. Similar problems have been widely studied in different fields of Physics, Mathematics, Computer Science, Signal Processing and Statistics. This is why before stating our approach to the problem we overview related ideas in various scientific disciplines. The most classical related problem is possibly the n-body problem which aims at forecasting the trajectory of n planets. The planets interaction is subject to the gravitation whose closedform expression is known. Nevertheless solving the system of these n dynamical equations turns out to be extremely challenging for n ≥ 3. In our problem the inherent dimensionality of the problem is larger that 3, the closed form expression of the interaction is unknown and only a partially observed and noisy observation is available.

In the following we review some of the related problems and ideas. We quickly recall some vector and matrix operators needed for introducing the objects of interest. For a matrix X = (X i,j ) i,j , we set the following matrix norms:

X 1 . = i,j |X i,j | and X * . = rank(X) i=1 σ i ,
where σ i are the singular values of X and rank(X) is the rank of X. We consider the following setup. In the sequel, the projection of a matrix Z onto S is denoted by P S (Z). The matrix (M ) + is the componentwise positive part of the matrix M, and sgn(M ) is the sign matrix associated to M with the convention sgn(0) = 0. The component-wise product of matrices is denoted by •. The class S + n of matrices is the convex cone of positive semidefinite matrices in R n×n . The sparsity index of M is M 0 = |{M i,j = 0}| and the Frobenius norm of a matrix M is defined by M2 F = i,j M 2 i,j . We also use M op = sup x : x 2 =1 M x 2 and M ∞ = max |M i,j |.

Underdetermined linear systems

Using first order local approximation of smooth functions, we can approximate a large class of problems by solving linear systems. In their general form, linear systems can be formulated as Xw = y where w ∈ R d is the unknown, X ∈ R N ×d is the covariate matrix, N is the number of observations, and y ∈ R N is the output vector. In such problems, the pair of observed input-output (X, y) are used to infer the vector w. In general X ∈ R N ×d is not invertible.

In case there are more observations than unknowns, d < N , the problem is said to be overdetermined. We may therefore prefer solving the problem X ⊤ Xw = X ⊤ y instead. When the number of observations is smaller than the number of variables N < d, the problem may have infinitely many solutions. Consequently, the optimization problem

min w∈R d Xw -y 2 (1.2)
suffers from numerical instability of solutions. In 1943, for stabilizing the solutions, Tikhonov suggested to solve the following regularized problem :

min w∈R p Xw -y 2 2 + Γw 2
data-dependent. In general, when no prior knowledge on data distribution is available, datadriven calibration procedures ( e.g. cross-validation) are often used. For further reading on the calibration methods we refer to [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]. In 1996, Tibshirani [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] has introduced an ℓ 1 -type penalized problem called the Lasso

min w w 1 subject to Xw -y 2 ≤ ǫ (1.4)
which despite not having a closed form solution, has the advantage of leading to vector solutions with few nonzero components. Hence the Lasso selects variables at the same time as it stabilizes the solutions. In fact the ℓ 1 -norm, defined as w 1 = d i=1 |w i | is the seen as the smallest convex surrogate of ℓ 0 . The ℓ 0 , sometimes called sparsity index of the vector, is defined by x 0 = d i=1 1{x i = 0}. It simply counts the number of non-zero elements of a vector. Minimizing the ℓ 1 -norm has the virtue of leading to sparse solutions as minimizing ℓ 0 would. Compared to the exponential cost of minimizing the ℓ 0 directly through a combinatorial approach (see [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]), minimizing the ℓ 1 norm has the advantage of defining the regularized problem as a convex optimization problem that can be solved in polynomial time.

Analysis of the Lasso

The analysis of the Lasso relies on the Restricted Isometry Constants first introduced by Candès and Tao [START_REF] Candès | Decoding by linear programming[END_REF]. These constants are useful for quantifying the quality of sparse recovery solutions obtained by ℓ 1 norm minimization [START_REF] Candès | An introduction to compressive sampling[END_REF]. The design of a low-dimensional feature map and the use of low-rank and sparse criteria enforced by the trace norm and ℓ 1 -norm relates our approach to the field of compressed sensing (see [START_REF] Candès | Decoding by linear programming[END_REF]) that specifically focuses on techniques that under-sample high-dimensional signals and yet recover them accurately. In simple words, compressed sensing says if a signal is sparse in some basis, then the measurement basis must be as incoherent as possible to the reconstruction basis. For instance if a signal is decomposed with few coefficients in the Fourier basis, then sinusoids are not the most suitable measurement basis, but one should rather choose a measurement basis that is the most de-correlated with sinusoids. This notion of uncorrelation of basis is formalized by the definitions of incoherence and restricted isometry.

Definition 1 (Restricted isometry constant) Given a matrix X ∈ R N ×d with N < d, for any integer s ≥ 1 we define the constant δ s as the smalest constant such that for all s-sparse (having s nonzero coefficients) vectors w ∈ R n ,

(1 -δ s ) w 2 2 ≤ Xw 2 2 ≤ (1 + δ s ) w 2 2 .
(1.5)

We start by stating two theorems from [START_REF] Candès | An introduction to compressive sampling[END_REF] Theorem 1 (Noiseless recovery) Suppose δ 2s < √ 2 -1. Let w s be a copy of w 0 where all but the slargest entries (in absolute value) of w 0 are set to 0. Let ŵ be the minimizer of w 1 subject to Xw = y where y = Xw 0 . Then there are constants C 0 , C 1 > 0 such that the two following hold true 1. Bound on the ℓ 1 error ŵw 0 1 ≤ C 0 w sw 0 1 2. Bound on the ℓ 2 error

ŵ -w 0 2 ≤ C 1 √ s w s -w 0 1 .
Theorem 2 (Robust sparse recovery in presence of noise) Suppose δ 2s < √ 2-1 and y-Xw 0 2 ≤ ǫ. Let w s be a copy of w 0 where all but the s-largest entries (in absolute value) of w 0 are set to 0. Let ŵ be the minimizer of w 1 subject to Xwy ≤ ǫ. Then there are constants C 0 , C 1 > 0 such that

ŵ -w 0 2 ≤ C 0 √ s w s -w 0 1 + C 1 ǫ .
A similar but less restrictive property called restricted eigenvalue has been introduced in [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF]. Before stating the definition, we introduce some notations. For J a subset of {1, • • • , d} of cardinality |J|, and a vector w ∈ R d , we denote by w J the vector of R d that is a copy of w over the set J and the value of coordinates indexed by k / ∈ J are set to zero. We denote by J c the complimentary of J in {1, • • • , d}, and introduce similarly w J c that is a copy of w over J c and vanishes on J. The notation • p is used to denote the ℓ p norm.

Definition 2 (Cone of restriction) For a positive integer s and a fixed real value c 0 > 0, we define the cone of restriction, as the set of vectors which values are concentrated on a subset of size at most s of the entries,

C(s, c 0 ) = {w ∈ R d |∃J ⊂ {1, • • • , d}, |J| ≤ s, w J c 1 ≤ c 0 w J 1 }
We now turn to the definition of the restricted eigenvalue of a design matrix X, which imposes a less restrictive condition on the covariates.

Definition 3 (Restricted eigenvalue) Let κ(s, c 0 ) . = min w∈C(s,c 0 ) Xw 2 √ N w J 0 2 . If κ(s, c 0 ) > 0, then X verifies RE(s, c 0 ).
With this definition in mind we can state the following result.

Theorem 3 Let w * an s-sparse vector (having s ≥ 1 nonzero elements), and y = Xw * + ǫ with ǫ ∼ N (0, σ 2 I N ) zero mean gaussian noise with variance σ 2 > 0, A > 2 √ 2 and γ = Aσ log d N . Define the Lasso estimator using the unconstrained formulation ŵ(γ) . = arg min

w∈R d 1 N y -Xw 2 2 + 2γ w 1 . (1.6)
Then provided that 1 N X ⊤ X has diagonal entries equal to 1 and verifies RE(s, 3), with probability at least 1d 1-A 2 /8 , by letting κ(s) = κ(s, 3), we have 1. Bound on the ℓ 1 error of the solution

ŵ(γ) -w * 1 ≤ 16A κ(s) 2 σs log d N 2.
Bound on the ℓ 2 error of the covariate times the solution

X( ŵ(γ) -w * ) 2 2 ≤ 16A κ(s) 2 σ 2 slog d 3.
Bound on the sparsity of the solution

ŵ(γ) 0 ≤ 64 1 N X ⊤ X op κ(s) 2 s .
We refer the interested reader to [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF] for proofs, further reading and also to [START_REF] Lounici | Oracle inequalities and optimal inference under group sparsity[END_REF] for extensions to group sparsity and multitask learning.

Practical procedures for solving penalized linear systems

The Tikhonov and ridge regression have closed form solutions

ŵ = (X ⊤ X + Γ ⊤ Γ) -1 X ⊤ y if X ⊤ X + Γ ⊤ Γ ∈ R d×d is
invertible, and this is the case for usual regularizer choices. Therefore, the algorithmic aspects of such estimation procedures rely on the linear system solver.

For minimizing the least squares problem penalized with an ℓ 1 term, note that the objective

J (w) . = 1 N y -Xw 2 2 + 2γ
w 1 is convex, ensuring thus the existence of a minimizer. The nondifferentiable ℓ 1 term makes the practical minimization of J intractable with standard gradient descent. Different strategies have been considered (see [START_REF] Schmidt | Optimization methods for ℓ 1regularization[END_REF] for a survey) for minimizing such objectives. A basic method [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] consists in cyclicly picking coordinates and minimizing the objective with respect to each coordinate. The benefit of this method is that the coordinate-wise updates are exact and the closed-form expression is easy and fast to compute. Despite computational efficiency at each step, when variables are highly correlated, the number of iterations needed to converge to the solution increases dramatically. Other methods start by smoothing the ℓ 1 norm and consider the minimization of regularizers as

d i=1 w 2 i + ǫ or d i=1 w 2 i α i + α i
that are respectively smooth approximation and upper bounds of w 1 for ǫ, α i > 0.

Proximal algorithms [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] have been designed for solving precisely such problems where the objective J = ℓ + γR can be written as the sum of 1. a Lipschitz differentiable and convex term ℓ(w)

2. a convex term R(w) who is not differentiable, but which proximal is fast to compute.

In the case of the ℓ 1 norm, the proximal operator is given by

prox θ . 1 (x) . = arg min z θ z 1 + 1 2 z -x 2 2 = sign(x i ) max(|x i | -θ, 0) i .
It has been proven [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] that the sequence

w k+1 = prox θγ . 1 (w k -θ∇f (w k ))
converges after O( 1 ǫ ) steps to a ball of radius ǫ of the minimizer of J . The step size θ is usually taken of the order of magnitude of the inverse of the Lipschitz constant of ∇ℓ. An accelerated algorithm (FISTA) that reaches the optimal convergence rate [Nes05] of O( 1 √ ǫ ) can be written using an auxiliary sequence [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Tseng | On accelerated proximal gradient methods for convex-concave optimization[END_REF]. The intuition behind the design of the latter algorithms relies on 1. the linear expansion of f around the point w k 2. the quadratic term L 2 ww k 2 2 that controls the closeness of the next step point w k+1 to the latter w k :

J (w) ≃ ℓ(w k ) + ∇ℓ(w k ) ⊤ (w -w k ) + γR(w) + L 2 w -w k 2 2 = L 1 2 (w -w k ) + 1 L ∇ℓ(w k ) 2 2 - 1 2L 2 ∇f (w k ) 2 2 + 1 L ℓ(w k ) + γ L R(w) = L 1 2 w -(w k - 1 L ∇ℓ(w k )) 2 2 + γ L R(w) + constant It follows that the point w k+1 = prox γ L R (w k -1 L ∇ℓ(w k )
) is a fair approximation to the minimizer of J around w k . The detailed analysis and extensions can be found in [START_REF] Tseng | On accelerated proximal gradient methods for convex-concave optimization[END_REF].

Matrix completion

The problem of matrix completion became popular thanks to challenges suggested by DVD retailers MovieLens and Netflix. It has been extensively studied on both theoretical and algorithmic aspects [SRJ05, CT09, Gro11, RR11, KLT11]. The matrix completion setting assumes that some ratings (1 up to 5 stars) are given by some users on some movies. The goal is to fill in the missing entries. The application of this problem is to recommend to each user a movie he does not know (did not yet rate) and that he is likely to give full score score (high predicted score). Formally, the problem can be formulated by seeing the rating process as linear measurements over the hidden users × products rating matrix. Let A, B = Tr(A ⊤ B) denote the inner product for matrices. If E i,j denotes the matrix having all zeros but a one at its (i, j)-th position, the function A → A, E i,j is linear and measures the (i, j)-th entry of A. Assume Ω ⊂ {1, • • • , #users} × {1, • • • , #movies} denotes the set of indices of the available ratings, i.e. if (i, j) ∈ Ω, then the rating of user i on movie j is available. For simplicity, take

n = max(#users, #movies), d = |Ω| , and denote by Ω = {(i 1 , j 1 ), (i 2 , j 2 ), • • • (i d , j d )} the elements of Ω. Define the linear mapping ω : R n×n → R d X → A, E i 1 ,j 1 , • • • , A, E i d ,j d .
Matrix completion aims at recovering a matrix A 0 given the observation y = ω(A 0 ). In particular, the theoretical studies explain for which matrices A 0 and which sampling sets Ω one can recover the matrix A. Statistical approaches study the noisy case when the observation is corrupted by an additive noise with mean zero ξ:

y = ω(A 0 ) + ξ .

Exact and robust recovery of low-rank matrices from few observations

Theoreticians [Bac08, CT09, Gro11, KLT11] have proved using tools comparable to those used in ℓ 1 -recovery literature, that if the rank of A is low and the sampling uniform at random, then minimizing convex functionals as

1 d ω(A) -y 2 2 + τ A * (1.7)
leads to consistent and optimal solutions for the matrix completion problem. To see the analogy with the sparse recovery, one can think of the set of singular values of a matrix as a vector whose sparsity is equivalent to the low-rank-ness of the matrix. This analogy between vectors and matrices is formalized using convex calculus (subgradient computations) in [START_REF] Lewis | The convex analysis of unitarily invariant matrix norms[END_REF].

Here we state the most recent results for the exact recovery and noisy matrix completion.

Exact recovery We consider the following problem

min A A * subject to ω(A) = ω(A 0 ) (1.8)
and wish to establish under which conditions its solution is unique and equal to A 0 .

Let us recall an easy linear lower bound on the number of required observations first. Note that even if the set of rank-r matrices of size n × n is constituted of matrices having n 2 coefficients, the low-rank constraint reduces the degree of freedom of the set. The number of free parameters of the set of n × n rank r matrices is bounded by

(n -1) + (n -2) + • • • + (n -r) left singular vectors + (n -1) + (n -2) + • • • + (n -r) right s.vectors + r s. values = r(2n -r) ≤ n 2
with equality only for full rank r = n. Therefore the minimum number of observation required to complete the matrix is at least r(2nr).

In order to establish sharp lower bounds on the minimum number of observations, note that a matrix with a single constant block has low rank but can hardly be recovered from entry-wise measurements. To formalize this intuition, Candès and Recht have introduced the incoherence of a matrix with the canonical basis, that can be extended [START_REF] Gross | Recovering low-rank matrices from few coefficients in any basis[END_REF] to any basis, which measures the discrepancy of the two bases and is related to the hardness of completion subject to linear observations. This notion has been defined first by Candès and Recht [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] and was refined by Gross [START_REF] Gross | Recovering low-rank matrices from few coefficients in any basis[END_REF]. Before introducing the definition, let A = U ΣV ⊤ be the SVD of A where U, V ∈ R n×r are orthogonal and Σ = diag(σ 1 , • • • , σ r ). Let U ⊥ and V ⊥ matrices of size n × (nr) ortho-normally completing the bases of U and V , and define

P U ⊥ = U ⊥ U ⊥ ⊤ P V ⊥ = V ⊥ V ⊥ ⊤
and also the orthogonal projections

P A (B) = B -P U ⊥ BP V ⊥ and P ⊥ A (B) = P U ⊥ BP V ⊥
Definition 4 (Coherence of a subspace) Let U be a subspace of R n of dimension r and P U be the orthogonal projection onto U . Then the coherence of U (vis-à-vis the standard basis (e i )) is defined to be

µ(U ) . = n r max 1≤i≤n P U e i 2
We now state the theorem

Theorem 4 (Exact matrix completion [CCSA08]) Let A = r m=1 σ m u m v ⊤
m be the SVD of A. Let U and V denote respectively spaces respectively spanned by left and right singular vectors. Assume the two following assumptions hold :

• A0 The subspace coherences obey max(µ( U ), µ( V )) ≤ µ 0 . • A1 The maximum entry of r m=1 u m v ⊤ m is bounded in absolute value by µ 0 r n .
Then the exact recovery of A given the observation of d entries selected uniformly at random is possible if

d > Crn log 2 n max(µ 2 1 , µ 1/2 0 µ 1 , µ 0 n 1/4 ) .
(1.9)

The latter results have been improved by David Gross, who generalizes the low-rank matrix completion framework to low-rank matrix recovery from few linear observations in any basis. We first recall the definition Definition 5 (Coherence of a matrix) Let (ǫ i ) n 2 i=1 be an orthonormal basis of R n×n , let A ∈ R n×n be a symmetric matrix with eigenvalue decomposition r m=1 σ m u m u ⊤ m , and let P T denote the projection onto the span of its eigenvectors. The matrix A has incoherence µ with respect to the basis

(ǫ i ) n 2 i=1 if either max i ǫ i op 2 ≤ µ n or the two estimates max i P U ǫ i ≤ 2µ r n max i ǫ i , r m=1 sign σ i u m u ⊤ m ≤ µ r n 2 hold.
A result due to Gross [START_REF] Gross | Recovering low-rank matrices from few coefficients in any basis[END_REF] improves the dependency in µ of the bound :

Theorem 5 (Exact matrix recovery in any basis

[Gro11]) Let ω : A → ( ǫ i 1 , A , • • • , ǫ i d , A ) where the set Ω = (i 1 , • • • i d ) is chosen uniformly at random in {1, • • • , n 2 },
then the solution to the exact recovery problem is unique and equal to A 0 with probability at least

1 -n -β if d ≥ O(nrν(1 + β) log 2 n) .
Trace-norm regularized regression. Consider the linear mapping

ω : A → ( A, X 1 , • • • , A, X d )
for a set of random matrices X i . Assume a corrupted observation y = ω(A 0 ) + ξ is available, ξ denoting a zero-mean noise. Let the solution of the unconstrained optimization problem

min A 1 d d i=1 E A, X i 2 - 2 d d i=1 y i A, X i + τ A * be denoted by Ŝ(τ )
, where the expectation is taken over the distribution of X i s. The following result is proven in [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] Theorem 6 (Trace-norm regression) Assume that there exist a constant µ > 0 such that

E ω(A) 2 2 ≥ µ -2 A 2 F . If τ > 2 1 d d i=1 (y i X i -Ey i X i ) ∞ , then Â(τ ) -A 0 2 F ≤ min 2τ A 0 * , 1 + √ 2 2 2 τ 2 µ 2 rank(A 0 ) (1.10)
The expectation being taken over the law of X i s.

In order to refine the last results, we need the two following definitions that are similar to the restricted cone and restricted eigenvalue property defined for the ℓ 1 -penalized regression.

Definition 6 (Cone of restriction) For a matrix A ∈ R n×n of rank r, we define the cone of restriction as the set of matrices that projection onto the singular spaces of A dominate the orthogonal projection, or the set of matrices nearly aligned with the row and column spaces of A. For c 0 > 0 let

C(A, c 0 ) . = {B ∈ R n×n | P ⊥ A (B) * ≤ c 0 P A (B) * }
We can now define the constant κ c 0 (A) for any matrix A by

Definition 7 (Restricted eigenvalue constant κ c 0 ) Let κ(A, c 0 ) . = inf B∈C(A,c 0 )\{0} 1 d E ω(B) 2 2 P A (B) F .
Koltchinskii et al. [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] have proved the following theorem.

Theorem 7 (Trace-norm regression)

If τ > 3 1 d d i=1 (y i X i -Ey i X i ) ∞ , then 1 d E ω( Â(τ ) -A 0 ) 2 F ≤ inf A 1 d E ω(A -A 0 ) 2 F + τ 2 κ(A, c 0 ) rank(A) (1.11)

Proximal gradient descent for regression with trace norm penalty

Proximal gradient descent algorithms [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] can be used for matrix completion. The proximal operator associated with the trace norm, also called the shrinkage operator [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] is given by:

D τ (B) . = arg min A 1 2 A -B 2 F + τ A * .
The practical computation of D τ requires the singular value decomposition of

B = U diag(σ 1 , • • • , σ n )V ⊤ ,
and is equal to

D τ (B) = U diag max(σ i -τ, 0) i V ⊤ .
To optimize

ℓ(A, B) + τ A *
where ℓ is a Lipschitz differentiable function, the sequence suggested by the proximal iterative gradient descent follows the update rule

S k+1 = D τ L S k - 1 L ∇ℓ(S k )
where L is the Lipschitz constant of ℓ. The corresponding accelerated algorithms exist. Although the memory requirement of FISTA may be a limitation in case of large values of n.

Link prediction heuristics

The problem of predicting missing links of a partially observed graph has a particular interest in many applications such as predicting hyperlinks of webpages [START_REF] Taskar | Link prediction in relational data[END_REF], finding protein-protein interactions [KKY + 09], studying social networks [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF], as well as collaborative filtering and recommendations [START_REF] Hu | Collaborative filtering for implicit feedback datasets[END_REF]. The study of the linkage structure of social networks dates back at least to Katz [START_REF] Katz | A new status index derived from sociometric analysis[END_REF]. More recently, Liben-Nowell and Kleinberg [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF] have popularized the problem of link prediction by comparing several classical methods on a scientific collaboration network. More recently, several works have built Bayesian and regularization-based types of algorithms and methods for dealing with the link prediction problem and variants [RBEV10, KX11, SM06, SSG07, SCM10a, MGJ09].

The link prediction problem has been studied in two different perspectives. The static version considers that some of the edges of a graph are hidden and one wants to reconstruct the partially observed graph. The temporal or dynamic version assumes that the edges appear over time on some fixed set of nodes, and the order of arrival of the edges is relevant. Therefore predicting the future edges requires to extract some temporal patterns that guide the apparition of the edges over time. Such mechanisms often comprise diffusion-like processes on graphs that are believed to handle specific temporal patterns observed in networks [START_REF] Leskovec | Graphs over time: Densification laws, shrinking diameters and possible explanations[END_REF]. The static link prediction heuristics have been studied and compared by Liben-Nowell and Kleinberg [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF]. We review some of the heuristics compared by them as well as more recent developments in this direction. Roughly the different algorithms can be grouped using the type of graph features they use:

1. Degree-based link prediction only relies on the distribution of node degrees, and predicts creation of links among high degree nodes.

2. Path-based link prediction uses the distance among nodes to predict connections among close nodes. Despite nice algorithmic properties, these algorithms are not flexible enough for reaching optimal performances.

3. Factorization based link prediction builds over the assumption that a small number of latent features describe the nodes. The nodes that are close in the space of latent factors are likely to be connected.

Degree-based link prediction

The Preferential Attachment Model relies on the observation that degrees in social networks follow a power-law distributions. Compared to the Poisson distribution of degrees in Erdos-Renyi random graphs, unbalances power-law distributions give higher importance to high degree nodes. In the context of recommendation, such an approach consists in recommending the most popular items to all users, or to all users of an identified cluster. This algorithm is inspired by the following graph generation process that produces graphs with power law degree distribution having a pareto coefficient β ∈ [2, 3].

Definition 8 (Preferential attachment model [AB02])

Starts from an initial node, and then follows iteratively two types of steps : 1-Vertex step-with probability p, add a new node, and connect it to a vertex with probability proportional to degrees 2-Edge step-with probability 1p, connect two nodes chosen with probability proportional to their degrees

The preferential attachment model is a recursive random graph generation process. This is why given some observed graph, it suggests edges to add to a the graph. Note that even though this simple model produces graphs with the desired degree distribution, neither the "small world" [AB02] effect, nor the high clustering coefficient property of real world graphs are reproduced.

In terms of recommender systems, applying such a model to an e-commerce market history -seen as a bipartite graph-is too simplistic. The items degrees in the bipartite graph of a market are in fact nothing but the sales volumes of the items. As a consequence, the preferential attachment model used in a recommendation contexts would promote best-selling items (high degree nodes) to all the users, which is very naive. For instance the movie market is composed of dramas, horrors, science fictions etc. and recommending a drama best-seller to science-fiction fanatics is a poor recommendation.

The next algorithms presented involve features of higher order than the degree. Nearest-Neighbors, Random-Walks and Katz approaches rely on the distribution of length of paths in the graph and emphasize the probability of establishment of links among nodes that are close to each other but are not connected yet.

Path-based link prediction

"People who bought this book also bought that one" is the sentence used by Amazon to introduce the recommended items to users. Such models link vertices that are reachable through short paths on the graph. Their success is explained by the high frequency of short paths and high clustering coefficients in many social and market networks. Algebraically, computing such paths requires to compute powers of adjacency matrices. Matrix multiplication using a naive algorithm has a complexity rising with n 3 , the fastest known algorithm for dense matrices runs in O(n 2.38 ). For sparse matrices such as most adjacency matrices, the complexity is independent of the size of matrices and scales linearly with their respective number of nonzero elements, and depends also on their sparsity patterns. Here we review some basic facts about walks on graphs and the corresponding matrices.

Proposition 1 (Powers of the Adjacency Matrix)

If A is the adjacency matrix of a graph G, then for all integer k, and all indices i and j, the (i, j) entry of A k , (A k ) i,j equals the number of walks of length k joining the vertex i to j .

The score we refer to as Katz [START_REF] Katz | A new status index derived from sociometric analysis[END_REF], first normalizes the adjacency matrix by the diagonal matrix formed by inverse degrees of nodes

W = diag(1/d 1 , • • • , 1/d n )A ,
and then sums powers of this stochastic matrix with decreasing exponential weights, ensured by 0 < β < 1:

∞ k=1 β k W k = βW (I n -βW ) -1 .

Factor based link prediction

Modeling graphs as strong connections among points in the factor space having close latent factors has been considered in building geometric random graph models, but also in link prediction and link analysis [START_REF] Hoff | Multiplicative latent factor models for description and prediction of social networks[END_REF][START_REF] Sarkar | Theoretical justifiation of popular link prediction heuristics[END_REF]. Most of these approaches build on Bayesian models for the factors, especially when dealing with temporal trends in the link formation and factor evolution [SM06, SSG07, YCY + 07, FAD + 11, VAHS11].

Link prediction

Matrix completion Observed entries no infomation (all equal 1) relevant Sparsity pattern relevant no infomation (uniformly random) Table 1.2: Comparison of the information sources in link prediction and matrix completion problems.

Given a loss function ℓ(A, B) measuring the fit of the observed adjacency matrix B to the score matrix A ∈ R n×n , the link prediction problem subject to low-dimensional factor assumption can be formulated as the minimization of ℓ(U V ⊤ , B) for U, V ∈ R n×k where k is the number of factors. Regularized algorithms of this type have been developed for collaborative filtering tasks [START_REF] Srebro | Maximum-margin matrix factorization[END_REF][START_REF] Hu | Collaborative filtering for implicit feedback datasets[END_REF]. They give iterative solutions by alternately minimizing the following objective in the two variables U and V :

ℓ(U V ⊤ , B) + λ( U 2 F + V 2 F )
despite efficiency and scalability [START_REF] Recht | Parallel stochastic gradient algorithms for large-scale matrix completion[END_REF], these approaches aim at minimize non-convex objectives and therefore the theoretical study of these algorithms is quite poor. The use of the singular value decomposition in the factorization of the matrix leads to well-posed formulations. Given a fixed target rank k one can compute the rank k matrix the closest to A in Frobenius norm by first computing A's SVD A = U diag(σ i )V ⊤ and then setting to 0 all but the first k singular values of A. That is

A k = U diag(σ 1 , • • • , σ k , 0, • • • , 0)V ⊤ .
A similar method is to solve the optimization problem

min A rank(A) subject to A -B F ≤ ǫ
that has a solution in terms of the SVD of A obtained by hard-thresholding the singular values of A. Using the tightest convex relaxation of the rank, the nuclear norm (sum of the singular values) leads to the following optimization problem:

min A A * subject to A -B F ≤ ǫ .
The corresponding lagrangian is

L(A) = 1 2 A-B 2 F +τ A * .
It is a strictly convex function of A and its minimizer is obtained by soft-thresholding the singular values of B,

D τ (B) = U diag(max(σ i -τ, 0))V ⊤
and is called the shrinkage operator [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF]. The advantage of the formulation using the nuclear norm is that as long as the loss function ℓ(A, B) = 1 2 A -B 2 F is replaced by other convex loss functions and the problem remains convex. The convex formulation allows to develop efficient algorithms and provide theoretical guarantees about the estimator.

Spectral functions

An interesting property pointed out by Kunegis et al. [START_REF] Kunegis | Learning spectral graph transformations for link prediction[END_REF] is that many of the standard link prediction algorithms can be written as spectral functions of the adjacency matrix or the normalized adjacency matrix. This characterization allows to search for the best scoring rules in a very low dimensional subspace of matrix functions if one affords to pay the cost of computing an SVD and storing possibly dense n × n matrices. Even though for large values of n such a process becomes computationally challenging, this characterization motivates the use of the trace norm as a regularization term in a link-prediction objective.

The score functions obtained by walks on graphs and by low-rank approximation of adjacency matrices belong to a specific subset of the functions mapping R n×n to itself called unitary invariant functions [START_REF] Lewis | The convex analysis of unitarily invariant matrix norms[END_REF]. A characterization of these functions goes as follows, if A = U diag(σ i )V ⊤ is the SVD of A, a unitary invariant function of A can be written as U diag(f (σ i ))V ⊤ for a real-valued function f : R + → R + . The powers of A and consequently their linear combinations, the projection onto the r highest singular vector space and shrinkage are all unitary invariant functions of the adjacency matrix. Note that the unitary invariant functions form a linear space of dimension n that is much lower than the dimension of the set of functions R n×n → R n×n , that equals (n 2 ) n 2 . Chapter 2

Description and simulation of marketing graph data 2.1 Exploratory analysis of graph data

Static graph descriptors

We borrow the mathematical vocabulary of graphs for representing relations. We denote by V = {1, • • • , n} the set of individual objects called nodes or vertices, and suppose that the overall number of nodes n is fixed. We denote by E the edges, links or connections among nodes. E can formally be seen as a subset of R V ×V . The basic case of binary graphs corresponds to a subset of {0, 1} V ×V , the cases where edges have non-negative attributes is called weighted graphs. We recall the definition of a graph and a bipartite graph for the sake of completeness.

Definition 9 (Graph and adjacency matrix) Given a set of vertices V = {1, • • • , n}, we define a graph by G = (V, E) where E ⊂ R V ×V is called the edge set. The adjacency matrix A ∈ R n×n of a graph G = (V, E) is defined as follows. Each element of E of the form e : (i, j) → a is reflected in A through its (i, j)-th element : A i,j = a.

For modeling the preferences of users among a set of items, we use bipartite graphs.

Definition 10 (Bipartite graph and incidence matrix) Given two disjoint sets of vertices

V 1 = {1, • • • , n 1 } and V 2 = {1, • • • , n 2 } , we call bipartite graph G = (V 1 , V 2 , E) where E ⊂ R V 1 ×V 2 is the edge set. We associated an incidence matrix M ∈ R n 1 ×n 2 to a bipartite graph G = (V 1 , V 2 , E). Each element of E of the form e : (i, j) → a is reflected in M through its (i, j)-th element : M i,j = a.
It is straightforward to see that a bipartite graph can be seen as a unipartite one over the vertex set V 1 ∪ V 2 . The adjacency matrix of the unipartite graph is given by

A = 0 M M ⊤ 0 ∈ R (n 1 +n 2 )×(n 1 +n 2 ) .
All along the current document we will use unipartite graphs and adjacency matrices unless the contrary is explicitly mentioned. We highlight the fact that this choice is only to simplify the notations and does not hurt the generality of the methodology. We will be interested by low-dimensional representations of the graphs. We call feature map the mapping of the graph to the vector space R d whose coordinates are descriptors of the graph.

Definition 11 (Graph Feature Map)

We refer to ω : R n×n → R d as the graph feature map.

It is important to keep in mind that by feature we refer indifferently to two different types of graph descriptors.

1. Global descriptors. These graph features describe the whole graph as an entity. For instance in the parametric graph models, the set of parameters generating the graph can be seen as global graph descriptors. The power-law exponents, the average clustering coefficient or the Gini index of the degrees are other examples of global descriptors.

2. Local descriptors. Various quantities, such as the degree, the pagerank, the number of length-k cycles passing by a node, are used to describe the nodes of the graph. We refer to these as local descriptors.

We only quickly review some of the features usually used to describe graph data.

Degree of nodes, power-law distribution, Lorenz curves. The degree of a node is the number of connection it establishes with other nodes. In the case of weighted graphs, the degree is defined as the sum of weights of the edges of a node. If the graph is directed the in-degree and out-degree denote respectively the sum of weights of edges from a node to its neighbors, and from its neighbors to it. If A is the adjacency matrix of a graph, its degree vector is given by A1 n , where 1 n is the all 1 vector of length n. note that the mapping A → A1 n is linear.

The distribution of degrees in many graphs describing human activities such as social networks, e-commerce marketplace graphs, information graphs such as wikipedia etc. are known to have long tails. It has been discussed and empirically verified ,see [START_REF] Clauset | Power-law distributions in empirical data[END_REF] or Figure (2.1.3), that many of them follow Pareto laws, sometimes called power-law or scale-free. The benefits generated by different markets also concentrate around lines in log-log plots (see Figure 2.4) when plotted as a function of sales volumes. Several types of behaviors can be distinguished. For instance books that are sold very few times generate more benefits than best-sellers, as opposite to video-games or electronic items. Bottom left: concave regions (books for instance): the more popular, the less revenue-generating the items are. Convex regions (electronic devices and video games) have the opposite behavior.

Definition 12 (Power-law)

The random variable X is said to follow a power-law distribution of parameter (x min , α) if its density function is of the form

p(x) = Cx -α for x > x min where C = (α -1)x α-1
min is a normalization constant, and p(x) = 0 if x ≤ x min .

Lorenz curve and Gini index. The characterization of the unbalance of long-tailed data is also done using the Lorenz curves, especially in economics. Lorenz curves allows to obtain the famous 80%-20% rule that says "80% of wealth is on the hand of 20% of the population" or "80% of the sales are generated by 20% of the items". The Lorenz curve represents the cumulative weight of the population as a function of the fraction of the population ordered in the increasing order of weight. If the weight measured is the wealth, then the Lorenz curve is the set of points which x coordinate represents to which fraction of the population an individual belongs in terms of wealth, and the y coordinate represents the percentage of cumulative wealth of the population that is at most as rich as x. If L : [0; 1] → [0; 1] is the Lorenz function (cumulative weight of the population as a function of the fraction of the population ordered in the increasing order of weight), and F ∈ [0; 1] represents the fraction of the population considered, then the Gini index is defined by

G = 1 -2 1 0 L(F )dF . It is straightforward to see that if X follows a (x min , α) power-law, then G(X) = 1 2α-1 .
Another well-knwon fact about Pareto law is its scale-invariance: for a constant c > 0, the two probabilities p(x) and p(cx) are proportional. In fact

p(cx) = c -α p(x)
Because of this property sometimes authors refer to social graphs having power-law degree distribution as free-scale graphs, for emphasizing that microscopic phenomena are repeated at larger scales in such graphs.

Adjacency matrix powers.

If A ∈ {0, 1} n×n is the adjacency matrix of an unweighted graph, one sees thanks to the formula

(A 2 ) i,j = n k=1 A i,k A k,j that A 2 coefficient (i, j
) is the number of paths from i to j, which coincides in the undirected case with the common neighbors of i and j. One can see similarly that A k coefficients are the number of paths of length k joining pairs of nodes, and in particular A k diagonal elements count the number of self-cycles of length k. In particular the diagonal elements of A 3 , counting the number of triangles surrounding a node have been made popular for measuring the homophily in a network. We call the clustering coefficient of a node the ratio number of triangles surrounding a node number of 'V'surrounding a node

= (A 3 ) i,i d i (d i -1)/2 .
A smoothed generalization of random walks leads us to PageRank, the famous vertex ranking algorithm suggested by Brin and Page [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF] that is said to be used in the Google searching engines for evaluating the relevance of each website based on the hyper-link structure of the internet web. See Figure 2.2 for an illustration of the distribution of some basic topological features of graphs. Clustering coefficients and order 1,2,3 degree (ratios) on each side of the bipartite graphs. For each node of a bipartite graph the clustering coefficient is the ratio of number of length 4 loops over the number of length 4 paths starting from that node. The higher values belong to Music and Books markets where homophily seems to play a dominant effect compared to Video Games and Electronic Device markets. Similarly the relatively higher importance of hubs in electronic and video games market can be read in the degree features distribution.

Definition 13 (Personalized PageRank) Let G be a graph, and A the corresponding adjacency matrix, and W = D -1 A the associated stochastic matrix. The stationary distribution of the random walk with teleport probability α ∈]0; 1[ to s, also called personalized PageRank vector is the solution of

pr α (s) = αs + (1 -α)pr α (s)W Proposition 2 The matrix R α = α ∞ k=0 (1 -α) k W k satisfies pr α (s) = sR α
the convergence being ensured by α ∈]0; 1[. Linear approximation to the feature map. Some graph features are non-linear functions of the adjacency matrix. For instance the number of triangles incident to a node i is given by (A 3 ) i,i , and therefore the clustering coefficient is a non-linear feature of the graph. In the current document, unless the contrary is mentioned explicitly, the feature map ω is linear. This simplifying hypothesis allows to build convex and differentiable objective functions in the Chapter 4. These feature maps ω should be considered as the tangent map to the actual feature map of interest. Further work should treat more general non-linear features.

Dynamics of graph descriptors

When the temporal evolution of the graph data matters for the application of interest, either directly or indirectly, the dynamic behavior of the graph features gets involved in the methodology. Rich information patterns are hidden in the temporal trends of graph descriptors. In subsection 1.2.2 we provided an example where the evolution of a linear feature vector of an adjacency matrix is used for valuation of customers of a C-to-C website. A more basic example (see subsection 1.2.3) is the characterization of product diffusion in a marketplace through the shape of the time series of its sales volume. The sales volume of a product being equal to its degree in the bipartite user-product graph, the diffusion level can be characterized using the time series of that graph feature. In the case of recommender systems discussed in subsection 1.2.4 it is the future value of the graph itself that is of interest. We will see that most heuristic link prediction algorithms rely on the assumption that a set of graph descriptors evolve smoothly. The prediction of edges becomes tractable when the relevant set of features is identified and their evolution well understood. In chapter 4 we explore regularized regression types of methods in order to encode this idea for a times series of graphs.

Marketplace seen as a dynamic bipartite graph

The graph describing a marketplace is a bipartite graph formed by the set of users and the set of products as disjoint node sets. The edges of the graph model the connections observed among the users and the products. The most frequent type of connection we are interested in is the transactions: the purchase of a product by a user is modeled as an edge between the corresponding nodes. In the sequel we provide a concrete idea of the basic characteristics of typical e-commerce data sets. We review properties of the databases that have been studied and illustrated in this thesis. The properties we report are chosen in order to clarify the objects of interest in the following sections, and are not exhaustive for a descriptive study of these databases. For more descriptive studies on information or internet domain graphs one could refer to [START_REF] Leskovec | Dynamics of large networks[END_REF] for instance.

Descriptive study of e-commerce markets

In order to give a concrete idea of the domain specificities, we provide description of four real-world markets with different characteristics. We refer to them as Music, Books, Electronic Devices,Video Games; the market products being categorized by database managers. The data are collected from January 2009 until February 2011 unless the contrary is mentioned, see basic statistics in table 2.1. In the remainder of this work we model the marketplace data as dynamic bipartite graphs. In fact we see each purchase or declaration of affinity through any other channel as a link between a user and an item. The data we deal with can be modeled as a list of pairs of linked objects.

For basic statistics on the size and typical characteristics of the data see table 2.1, where we used the following notations. n 1 : number of users, n 2 : number of products, ∆E : daily number of sales, αU and αP estimated power coefficient for users and products. More imbalanced distributions appear in Electronic and Video Games markets. 2 and3). The spikes in sales series drastically increase the amount of input data but, in the same time, they introduce a bias in the collection of input data. We define the penetration[SJT09] of a product over a market as the normalized cumulative volume of sales. The shape of penetration of products in different markets has an impact on the accuracy of the prediction of future popularity. For instance on the Video Games market, the shortness of product cycles make the prediction of popularity of the corresponding products intractable if based only on their purchase history and without injecting exogenous information. Benefits (euros)

Domain

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 1 5 10
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • BOOKS MUSIC VIDEO GAMES ELECTRONICS COMPUTER HIFI WHITE BABY CLOTHING VEHICLE SPORT WINE 1e+00 1e+02 1e+04 1e+06 1e+04 5e+04 2e+05 5e+05 2e+06 5e+06
Cumulative sales volume (increasing order)

Cumulative benefits (euros)

• • • • • • • •• •• ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• ••••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • BOOKS MUSIC VIDEO GAMES
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •••••••••••• •••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••••• ••••••• ••••• • • • • • • • • • • • BOOKS MUSIC VIDEO GAMES ELECTRONICS COMPUTER HIFI WHITE BABY CLOTHING VEHICLE SPORT WINE

Social networks

We have performed experiments with the Facebook100 data set analyzed by [START_REF] Traud | Social structure of facebook networks[END_REF]. The data set comprises all friendship relations between students affiliated to a specific university, for a selection of one hundred universities.

Protein interactions

We use data from [HJB + 09], in which protein interactions in Escherichia coli bacteria are scored by strength in [0, 2]. The data is, by nature, sparse. In addition to this, it is often suggested that interactions between two proteins are governed by a small set of factors, such as surface accessible amino acid side chains [START_REF] Bock | Predicting protein-protein interactions from primary structure[END_REF], which motivates the estimation of a low-rank representation.

Synthetic data sets

In this section we introduce new generative models for random graphs. Further reading about standard generative models of static random graphs can be found in [START_REF] Kolaczyk | Statistical Analysis of Network Data[END_REF]. The motivation for developing new generative random graph models is to simulate situations where the working assumptions are rigorously verified in order to have an understanding and evaluation of the algorithms performance before testing them on the real data. In the sequel we review standard random graph models and develop new models that incorporate desired properties. In the static case we build sparse low-rank matrices and add a sparse random noise in order to

Lorenz Curves (Books)

Fraction of Products u Fraction of Sales by Lowest 100u% Products obtain a sparse noisy version of the signal matrix, and in the dynamic case we build a sequence of sparse low-rank adjacency matrices having linear feature vectors evolving in an autoregressive way. In Table (2.7) page 59 we bring together in a single table both the standard baseline models and the generative models introduced in the current thesis and also the corresponding estimation methodology.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 2002

Gini Coefficients Evolution

years Gini Coefficient • • • • • • • • • • • BOOKS MUSIC VIDEO GAMES ELECTRONICS COMPUTERS

Generative models of static random graphs

Erdos-Renyi graphs. The most basic model of random graphs has been introduced by Erdös and Renyi [START_REF] Erdos | On the evolution of random graphs[END_REF]. In this model any pair of nodes have equal probability to be linked by an edge. The only parameter that characterizes such a graph is the probability of presence of a link or equivalently the total number of edges. We denote by G (n, p) 2 ) have equivalent properties. The following facts are known about Erdos-Renyi graphs, we refer the interested reader to [START_REF] Bollobas | Random graphs[END_REF] for more details.

a

Property 1 (Degree distribution)

The degree distribution of G (n, p) follows a binomial law

P(d = k) = p k (1 -p) n-k-1 n -1 k
provided that pn → λ, the asymptotic degree distribution of a random Erdos Renyi graph follows a Poisson law with parameter λ : Preferential attachment. This generative model has been designed to reproduce the popularity distribution observed in real-world networks. Its temporal formation process mimics the rich-get-richer phenomenon. This iterative generative model starts with an initial graph G 0 = {V 0 , E 0 } and adds at each iteration a new vertex that connects to the older vertices with probability proportional to their degrees : high degree nodes are preferred to low-degree nodes. For further details, variants and links to some problems in statistical physics refer to [START_REF] Albert | Statistical mechanics of complex networks[END_REF].

P(d = k) ≈ λ k e -k k! Property
Geometric random graphs. This family of random graphs [START_REF] Penrose | Random geometric graphs[END_REF] is built as follows:

1. draw points randomly from a distribution over a subset of R r 2. define the vertices of the graph as the sampled points 3. connect pairs of vertices by an edge if the distance of the sample points is lower than a fixed threshold.

These graphs are often drawn from uniform distribution and over cubes or spheres.

Exponential random graph models (ERGM) [WP96]

. These models see the graph z as the realization of a random variable Z. The probability mass distribution of these models is assumed to have the form

P θ (A = a) = e θ ⊤ ω(a)-Φ(θ) (2.1)
where θ ∈ R d is a vector of parameters, g a function of the graph a : ω(a) ∈ R d and Ψ a normalization factor. In practice, if a set of parameters of interest θ (such as the number of triangles, cycles or other patterns of interest) is identified as being the sufficient statistics for estimating the graph, this model turns out to fit the right parameters by log-likelihood maximization easily.

Stochastic block model. This model introduced by Nowicki and Snijders [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] suggests that k blocks or clusters form the structure of the graph. This assumption is based on the observation of presence of communities in social networks. Some extensions of this work have relaxed the assumption of the existence of rigid blocks to overlapping[LBA09, ABFX08], or time-evolving blocks [START_REF] Ho | Evolving cluster mixed-membership blockmodel for time-varying networks[END_REF].

Latent factor models. This family of models introduced by Hoff et al. [START_REF] Hoff | Latent space approaches to social network analysis[END_REF] assumes that given some node latent factors z i , z j and possibly pair-wise covariate descriptors x i,j the edges can be assumed to be independently distributed :

P(A|X, Z, θ) = i =j P(A i,j |X i,j , Z i , Z j , θ) (2.2)
for a parameter θ and the n × n random variable Y denoting the graph adjacency matrix. The authors mention that a convenient choice of modeling, that relates this model to the exponential random graph model, may simplify the estimation of parameters to a log-likelihood maximization. They take the example of log-likelihood depending linearly on the euclidean distance : log odd(a i,j = 1|x i,j , z i , z j , α, β) ∝ αβx i,j + z iz j 2

Intuitively, the closer the latent factors z i and z j are, the higher the probability of observing a link between nodes i and j is. For further reading and a dynamic extension refer to [START_REF] Krivitsky | A Separable Model for Dynamic Networks[END_REF].

Linear time series analysis

We call time series, or sequence of real or vector valued random variables indexed by time.

We need to recall some basic models devoted to the study of time series before introducing generative models of dynamic random graph sequence. Predicting the future value of a time series is a challenge of interest for a wide range of applications as meteorology, economics, finance and supply chain management. Various models have been developed focusing each on time series with particular properties. We recover a family of models and methods used in the current thesis and refer to [Tsa05, BD09] for further reading. We start by a basic definition.

For simplicity, we index all time series with non-negative integers t ∈ N.

Definition 14 (Stationarity)

There are two forms of stationarity:

• If the joint distribution of any vector

(ζ t-1 , . . . , ζ t-k ),
where k < t, remains constant over time, the the time series is strongly stationary.

• The weak form only assumes that the mean µ t = E(ζ t ) and the covariance between lagged values

Σ t-k,t = Cov(ζ t-k , ζ t ),
where k < t, are time-invariant.

The weak form of stationarity is a more realistic assumption than the strong form. Linear econometric models have thus been developed for weakly stationary time series. We call linear time series a time series that can be written as

ζ t = µ + ∞ i=0 ψ i ǫ t-i
with ψ 0 = 1 and ǫ t-i stand for zero-mean random variables. Among linear models we focus on some particularly simple and interesting families.

Auto-Regressive models.

Definition 15 (AR(p) models) A time series (ζ t ) t is said to be autoregressive of order p if it follows a recursive relation of the form

ζ t = φ 0 + p i=1 φ i ζ t-i + ǫ t ,
for real coefficients (φ i ), where (ǫ t ) t are i.i.d. zeros-mean noise often called shocks, innovations, or innovative residuals.

To see the link with the definition of linear time series, consider the case of AR(1). In this case

ζ t = φ 0 + φ 1 ζ t-1 + ǫ t .
Therefore the average value of the series is µ = φ 0 +φ 1 µ or µ = φ 0 1-φ 1 . The model can be written as

ζ t = φ 0 1 -φ 1 + ∞ i=0 φ i 1 ǫ t
which is the equation of a linear time series.

The AR(p) model implies the following linear recurrent equation between the lag-i autocovariance ρ(i) = Σ t-i,t :

ρ k = p i=1 φ i ρ k-i .
The associated polynomial equation is called the characteristic equation of the model:

Q(x) = 1 - p i=1 φ i x i = 0. Let (ζ -1
i,⋆ ) 1≤i≤p denote the inverses of the roots of the polynomial Q. These complex numbers are known as characteristic roots. It can be shown that a sufficient and necessary condition ensuring the weak stationarity of the AutoRegressive process (ζ t ) is that its the characteristic roots are smaller than 1 in modulus, i.e. |ζ i,⋆ | > 1.

Moving Average. It is interesting to note that Moving-Average models can be seen as a class of AutoRegressive models of infinite order:

ζ t = φ 0 + +∞ i=1 φ i ζ t-i + ǫ t ,
which, under specific parameter constraints, reduce to the expression:

ζ t = θ 0 + q i=1 θ i ǫ t-i ,
where q is said to be the order of the Moving-Average model MA(q). Auto-Regressive Moving-Average. AutoRegressive and Moving-Average models can be refined by introducing a combination of both, called AutoRegressive Moving-Average (ARMA) models. An ARMA process of order (p, q) is given by

ζ t = φ 0 + p i=1 φ i ζ t-i + a t + q i=1 θ i a t-i .
Usually, this mixed form requires fewer parameters (than pure AR or MA models) to describe dynamic weak stationary data.

Besides, the condition for weak stationarity, computed with the AutoRegressive part of the ARMA model, is the same as it would be for this pure AR model. Therefore, ζ t is weakly stationary if the characteristic roots computed with the autoregressive part are smaller than 1 in modulus. As a consequence, the lag-i autocorrelation ρ(i) of an ARMA process converges exponentially to 0 as the lag i increases. Yet, specific time series can display a decay of the autocorrelation at a polynomial rate. These long-memory effects can be reproduced by fractionally differenced processes.

Moreover, ARMA models can be extended by allowing a fixed number of characteristic roots to equal 1. Simple examples are the random walk, defined by

ζ t = ζ t-1 + a t ,
and the random walk with drift

ζ t = µ + ζ t-1 + a t ,
where µ is a time trend of the series (ζ t ) and is called the drift of the model. More generally, AutoRegressive Integrated Moving-Average (ARIMA) models of order (p, d, q) are obtained by differencing d times the series ζ t and by modelling the residuals with an ARMA model of order (p, q).

Estimating autoregressive models by ordinary least squares. When dealing with real valued time series, a simple strategy is to fit linear models to the time series. The autoregressive models attempt at establishing a linear relation among successive values of a time series. We call lagged time series, those which values at time t depends on values on times

t -1, • • • , T -l + 1.
In the framework of autoregressive models, one would model the linear dependencies of such series as

ζ t = l k=1 ζ t-k a k + ǫ t
where ǫ t is zero mean noise. Note that the vector

z t . =      ζ t ζ t-1 . . . ζ t-l+1      ∈ R l
follows the autoregressive relation

z t =        a 1 a 2 a 3 • • • a l 1 0 0 • • • 0 0 1 0 • • • 0 0 0 . . . 0 0 0 0 • • • 1 0        z t-1 .
The natural generalization of autoregressive models to lagged series or multivariate series leads to vector autoregressive models. VAR models fit a vector time series Z t ∈ R d to its past value Z t-1 through a linear relation of the form

Z t = W ⊤ Z t-1 + E t
where E t sands for the noise vector, that we assume having zero mean and independent from its past values. Notice that by introducing the slack variables

X =    Z ⊤ T -1 . . . Z ⊤ 0    ∈ R T ×d and Y =    Z ⊤ T . . . Z ⊤ 1    ∈ R T ×d
the regression problem aims at estimating W ∈ R d×d in the relation Y = XW + E E denoting the noise. One easily gets the ordinary least squares estimate

W o.l.s. = (X ⊤ X) -1 X ⊤ Y .

Generative models of dynamic graphs

In this section we introduce some models of dynamical graph that despite being extremely simplified still incorporate some desired characteristics. They have two main properties:

• A low-dimensional vector, called latent factors vector, evolves smoothly over time and dictates the evolution of the linkage structure of the graph

• Each noiseless version of the graph is simple at any time, meaning that a community structure exists in every snapshot of the graph. This structure is reflected onto the adjacency matrix as the low-rank property of the matrix. This prior knowledge will be used later for solving the prediction problem in an inverse problem fashion.

We built the generative models inspired by the standard random graph models, especially geometric random graph models that are related to the literature on low-rank matrices, low-rank matrix completion, and factor models. We will argue why they benefit from several properties of former generative models, making them realistic candidates for time evolving graph models. The construction process of these models is closely related to the hypotheses made for designing our main objective.

Dynamical system of latent factors. This model is close in spirit to the latent factor model in that a set of r independent latent factors guide the structure of the network. Let r be an integer and lets refer to R r as the space of latent factors. Suppose a set of n points

Ψ (1) (t), • • • , Ψ (n) (t) ∈ R r
evolve over time following the same differential equation

∂Ψ (i) (t) ∂t = f ⋆ (Ψ (i) (t)) (2.3)
where f ⋆ : R r → R r is continuous. The observation of latent factors is by assumption impossible. For instance, in the case of social network data they may correspond to the psychological state of the users that evolves smoothly over time following some unknown but still stable laws.

We assume that we are also given an attraction or similarity function

k : R r × R r → R
that is not necessarily symmetric, We observe a link between nodes i and j at time t if

k(Ψ (i) (t), Ψ (j) (t)) ≥ 1 .
S-curve latent factors. Take the case of directed graphs and say

Ψ(t) = [U t , V t ] where ∀k ∈ {1, • • • , r}, U (i,k) t = 1 √ 2πσ i,k e -(t-µ i,k ) 2 2σ 2 i,k + ǫ i,k
suppose V t s have similar behaviors, and take

k(Ψ (i) (t), Ψ (j) (t)) = 1{U (i) (t) > θ}1{V (j) (t) > θ} ⊤
for some real valued threshold θ and parameters µ i,k , σ i,k > 0, additive noise ǫ ∼ N (0, η 2 i,k ). Hence the data observed can be modeled as a discrete sequence of graphs G 1 , • • • , G T where for all t ∈ {1, • • • , T }, the set of nodes is V = {1, • • • , n} and the set of edges of G t is defined by the latter equation. The main motivation of such a model comes from the connection it establishes with the marketing literature on the cycle of life of products (see for instance [START_REF] Sood | Functional regression: A new model for predicting market penetration of new products[END_REF]) having S-like curves, and the factor model and geometric random graphs. Discrete time dynamical system of latent factors. In a similar fashion we introduce a discretetime generative model for the graph sequence and the node-related features without fixing the closed-form expression of the feature values, but rather by discretising a dynamic model and fixing the expression of f ⋆ . In fact notice that if for some smooth function f ⋆ , ∂ ∂t z

(i) t = f ⋆ (z (i) t ) (2.4)
then in the discrete time domain it means that ∀i, z

(i) t -z (i) t-1 = ηf ⋆ (z (i) t-1 ) (2.5)
where η > 0 is the time step unit.

         X t = A t Ω A t = U t V ⊤ t + z t ∀i U (i) t = U (i) t-1 + ηf ⋆ (U (i) t-1 ) + u t,i V (i) t = V (i) t-1 + ηf ⋆ (V (i) t-1 ) + v t,i
where u t,i , v t,i are vectors representing noise in the latent factor space, for instance drawn from N (0, δ 2 I r ) in R r , and the entries of z t are independently drawn from a centered Gaussian with variance σ 2 and represent noise at edge level. The important assumption on f ⋆ : R r → R r is continuity.

In the synthetic datasets used in our numerical experiments we took

f ⋆ (x) = e -x-v 1 2 σ 2 1 (x -v 1 ) + e -x-v 2 σ 2 (x -v 2 ) ,
where vectors v 1 , v 2 ∈ R r + are chosen randomly and ǫ, σ 1 , σ 2 are positive constants. We used n = 100, r = 4, T = 60 and the entries of U 0 , V 0 ∈ R n×r are drawn according to a uniform distribution in (0, 1).

Low-rank matrix sequence having linear autoregressive features (LR-LinFeat).

Let V 0 ∈ R n×r be an orthogonal matrix i.e. V ⊤ 0 V 0 = I r . Take W 0 ∈ R r×r be a square matrix, and fix also U 0 ∈ R n×r . Define the sequence of matrices (A t ) t≥0 by for t = 1, 2, • • •

U t = U t-1 W 0 + N t and A t = U t V ⊤ 0
for zero-mean i.i.d noise N t , for example gaussians, (N t ) i,j ∼ N (0, σ 2 ) If we define the linear feature map ω(A) = AV 0 , note that 1. The sequence ω(A t ) t = U t t follows the linear autoregressive relation

ω(A t ) = ω(A t-1 )W 0 + N t 2.
For any time index t, the matrix A t has rank at most r the requirement not fulfilled by this matrix sequence is sparsity of A t .

Low-rank and sparse sequence with thresholded linear features. Let V 0 ∈ R n×r be a sparse matrix and V † 0 its pseudo-inverse

V † 0 V 0 = V ⊤ 0 V ⊤ † 0 = I r . Fix the two matrices W 0 ∈ R r×r and U 0 ∈ R n×r . For t = 1, 2, • • • let U t = 1{U t-1 W 0 + N t > θ} and A t = U t V ⊤ 0 + M t for white gaussian i.i.d noise (N t ) i,j ∼ N (0, σ 2 ).
Now define the linear feature map ω(A) = AV ⊤ † 0 and note that 1. The sequence ω(A t ) t = U t t does not follows linear autoregressive relation. In fact

ω(A t ) -ω(A t-1 )W 0 = 1{U t-1 W 0 + N t > θ} -U t-1 W 0
2. ∀t, A t has rank at most r 3. ∀t, A t is sparse as V ⊤ 0 and U t are both sparse The advantage of this model is to come up with sparse and low-rank matrices, but the linear features do not follow a linear autoregressive relation. One may possibly argue that in some particular cases their asymptotic behavior approaches the linear autoregressive one. Notice also that we avoided thresholding a low rank dense matrix for obtaining a sparse and low-rank matrix. In fact a basic example suffices to get convinced that thresholding low-rank matrices does not result in low-rank matrices. Take u > 1 and v = (u, • • • , u n ) ⊤ ∈ R n , and A = vv ⊤ = (u i+j ) i,j . The matrix 1{A > θ} for θ = u n is anti-triangular and has full rank n despite the minimal rank 1 of A.

The following model attempts to fulfill all the requirements.

Sparse low-rank sequence having permutation-autoregressive features. Let ν : R n×r × N be a discrete noise function, such that ν(U, σ) is a copy of U where σ pairs of elements have been chosen uniformly at random and permuted. Let

• V 0 ∈ {0, 1} n×r be a sparse matrix and V † 0 its pseudo-inverse V † 0 V 0 = V ⊤ 0 V ⊤ † 0 = I r .
V † 0 has no reason to be sparse

• W 0 ∈ S(r) a permutation matrix on r elements: W 0 ∈ {0, 1} r×r is invertible with W -1 0 ∈ {0, 1} r×r • U 0 ∈ {0, 1} n×r
Positive semi-definite low rank matrices with continuously evolving features. The motivation for introducing this model is to add PSD condition to the low-rank graphs, which appear when studying the graphical model underlying high-dimensional distributions [START_REF] Zhou | Time varying undirected graphs[END_REF][START_REF] Kolar | Estimating time-varying networks[END_REF]. Let z (i) t ∈ R r denote latent factors of each node i at time t. Suppose we observe the sequence of covariances

Σ t = z (i) t , z (j) t 1≤i,j≤n 1≤t≤T
.

Equivalently if we denote by

Z t =     z (1) ⊤ t . . . z (n) ⊤ t     ∈ R n×r
the matrix of latent factors, we have

Σ t = Z t Z ⊤ t
Suppose there exists a function f ⋆ : R r → R r such that in the continuous time domain

∂ ∂t z (i) t = f ⋆ (z (i) t )
Smoothness assumption. The function f ⋆ assumed to be smoothly differentiable. In the discrete time domain it means that ∀i, z

t-1 = ηf ⋆ (z

(i) t-1 )
where η > 0 is the time step unit.

Positive semi-definite low rank matrices with linearly evolving features. Take

Σ t = Z t Z ⊤ t + M t M ⊤ t
where for some

W 0 ∈ R d×d Z t = Z t-1 W 0 + N t
   A t = U t V 0 + M t U t = U t-1 W 0 + N t ω(A t ) = A t V ⊤ 0 V 0 orthogonal, M t , N t zero-mean i.i.d. noise 1. A t low-rank 2. ω(A t ) linear autoregressive Minimizing L for γ = 0, pen = 1 2 . 2 F SPLR-LinFeat    A t = U t V 0 + M t U t = U t-1 W 0 + N t ω(A t ) = A t V ⊤ † 0 V 0 , W 0 , sparse M t , N t zero-mean sparse i.i.d. noise 1. A t sparse 2. A t low-rank 3. ω(A t ) linear autoregressive Minimizing L for pen = . 1 DynSysLF              A t = U t V ⊤ t + z t ∀i          U (i) t = U (i) t-1 + ηf ⋆ (U (i) t-1 ) + u t,i V (i) t = V (i) t-1 + ηf ⋆ (V (i) t-1 ) + v t,i f ⋆ smooth 1. A t low

Chapter 3

Estimation of sparse low-rank matrices

Context

Matrix estimation is at the center of many modern applications and theoretical advances in the field of high dimensional statistics. The key element which differentiates this problem from standard high dimensional vector estimation lies in the structural assumptions which are formulated in this context. Indeed, the notion of sparsity assumption has been transposed into the concept of low-rank matrices and opened the way to numerous achievements (see for instance [START_REF] Srebro | Learning with matrix factorizations[END_REF][START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF]). In this chapter, we argue that the low-rank-ness is not only an equivalent of sparsity for matrices but that being low-rank and sparse can actually be seen as two orthogonal concepts. The underlying structure we have in mind is that of a block diagonal matrix. This situation occurs for instance in covariance matrix estimation in the case of groups of highly correlated variables or when denoising/clustering social graphs. Efficient procedures developed in the context of sparse model estimation mostly rely on the use of ℓ 1 -norm regularization [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Natural extensions include cases where subsets of related variables are known to be active simultaneously [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. These methods are readily adapted to matrix valued data and have been applied to covariance estimation [EK09, BT10] and graphical model structure learning [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. In the low-rank matrix completion problem, the standard relaxation approach leads to the use of the trace norm as the main regularizer within the optimization procedures [START_REF] Srebro | Maximum-margin matrix factorization[END_REF][START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] and their resolution can either be obtained in closed form (loss measured in terms of Frobenius norm) or through iterative proximal solutions [CP11, BT09] (for general classes of losses). However, solutions of low-rank estimation problems are in general not sparse at all, while denoising and variable selection on matrix-valued data are blind to the global structure of the matrix and process each variable independently.

In this chapter, we study the benefits of using the sum of ℓ 1 and trace-norms as regularizer. This sum of penalties on the same object allows to benefit from the virtues of both of them, in the same way as the elastic-net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] combines the sparsity-inducing property of the ℓ 1 norm with the smoothness of the quadratic regularizer. Trace norm and ℓ 1 penalties have already been combined in a different context. In Robust PCA [START_REF] Candès | Robust principal component analysis[END_REF] and related literature, the signal X is assumed to have an additive decomposition X = S+L where S is sparse and L lowrank. Note that X is not in general sparse nor low-rank and that this decomposition is subject to identifiability issues, as analyzed, e.g., in [START_REF] Chandrasekaran | Rank-sparsity incoherence for matrix decomposition[END_REF]. The decomposition is recovered by using ℓ 1 -norm regularization over X and trace norm regularization over Y . This technique has been successfully applied to background substraction in image sequences, to graph clustering [START_REF] Jalali | Clustering partially observed graphs via convex optimization[END_REF] and covariance estimation [START_REF] Luo | High dimensional low rank and sparse covariance matrix estimation via convex minimization[END_REF].

Here, we consider the different situation where the matrix S is sparse and low-rank at the same time. We demonstrate the applicability of our mixed penalty on different problems. We develop proximal methods to solve these convex optimization problems and we provide numerical evidence as well as theoretical arguments which illustrate the trade-off which can be achieved with the suggested method.

A leading insight for inferring graph structures coming from social activities is that the informative graph is highly structured [START_REF] Girvan | Community structure in social and biological networks[END_REF] while the noise can be modeled as a parsimonious random pattern, such as an Erdős-Rényi random graph. A variety of approaches have been developed for the link prediction problem from a snapshot or a sequence of snapshots of a graph using matrix factorization [START_REF] Koren | Factorization meets the neighborhood: a multifaceted collaborative filtering model[END_REF], probabilistic methods [START_REF] Taskar | Link prediction in relational data[END_REF] or heuristics [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF][START_REF] Sarkar | Theoretical justification of popular link prediction heuristics[END_REF]. We treat the temporal case in chapter 4 of this thesis. Most methods assume, however, that the observed links are reliable and try to predict links that are either missing or, in a time-dependent setting, that are likely to appear in the future. As a consequence, these approaches can not remove irrelevant links and can be sensitive to their presence. Another feature of these relational databases is that both the informative component and the noise are sparse. In terms of matrix recovery, the setup corresponds to the case where the observed adjacency matrix is a sum of a simultaneously sparse and low-rank matrix (highly structured part) and a sparse matrix (noise) that is potentially full-rank with an almost flat spectrum [START_REF] Van Mieghem | Graph spectra for complex networks[END_REF].

Recently, Candès et al. suggested in [START_REF] Candès | Robust principal component analysis[END_REF] the robust principal component analysis (RPCA) approach for decomposing a matrix as the sum of a sparse matrix and of a low-rank matrix. A closely related method has been applied to graph clustering [START_REF] Jalali | Clustering partially observed graphs via convex optimization[END_REF]. In spite of computationally more attractive variants [LGW + 09, ZT11], RPCA and its derivatives have high computational and memory requirements which may be unsuitable for the analysis of massive datasets. More important, this approach has severe limitations when the sparse component turns out to be low-rank or when the low-rank component is itself sparse. This intuition has been formalized as the rank-sparsity uncertainty property in [START_REF] Chandrasekaran | Rank-sparsity incoherence for matrix decomposition[END_REF]. The setup considered in this chapter corresponds exactly to the case where the RPCA approach fails to perform accurate recovery. We propose a method, explicitly designed to recover matrices that are simultaneously sparse and low-rank when corrupted by sparse noise. Recovering simultaneously sparse and low-rank matrices was addressed recently in the case of the squared loss in [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF] for well concentrated and dense noise. In the latter case, the natural loss for the estimation error is the ℓ 1 -norm [START_REF] Bassett | Asymptotic theory of least absolute error regression[END_REF][START_REF] Gao | Asymptotic analysis of high-dimensional lad regression with lasso[END_REF]. We develop two algorithms to solve relaxations of the graph denoising problem and demonstrate their efficiency. The first algorithm that we propose is based on a convex formulation. It enforces both sparsity and low-rank on the estimator under a sparsity inducing loss function. We also present a scalable method based on matrix factorization, where we perform rank-one updates with a flavor similar to sparse principal component analysis (SPCA, see, e.g., [START_REF] D'aspremont | A direct formulation for sparse pca using semidefinite programming[END_REF]).

The remainder of the chapter is organized as follows. After motivating the problem in section 3.1, the regularized approach is described in Section 3.2. The oracle bounds for regression with the introduced penalty are provided in Section 3.3. The algorithms for the optimization of both the differentiable and non-differentiable loss are introduced in Section 3.4. Section 3.5 contains numerical results which illustrate the potential of these algorithms.

Setup and motivations

Notations

We first set some notations. For a matrix X = (X i,j ) i,j , the Frobenius norm is X 2

F . = i,j X 2 i,j , the matrix ℓ 1 -norm is X 1 . = i,j |X i,j | and the nuclear norm is defined by X * . = rank(X) i=1 σ i ,
where σ i are the singular values of X and rank(X) is the rank of X. The sparsity index of X or ℓ 0 pseudo-norm is nnz(X) = X 0 = |{X i,j = 0}|. We shall also use X op . = sup w : w 2 =1 Xw 2 and X ∞ . = max |X i,j |. The inner product of two matrices X and Y of the same size is defined by X, Y = Tr(X ⊤ Y ) . The component wise product of matrices is denoted by •. The class S + n of matrices is the convex cone of positive semidefinite matrices in R n×n . We denote the projection of a matrix X onto a convex set S by P S (X). The matrix (X) + is the componentwise positive part of the matrix X, and sign(X) is the sign matrix associated to X with the convention sign(0) = 0.

Motivations

We consider a matrix regression model with observations of the form

y i = Ω i , X 0 + ε i , i = 1, . . . , d (3.1)
where ε = (ε i ) is i.i.d. noise, X 0 ∈ R n×n is an unknown matrix to be estimated and

Ω 1 , . . . , Ω d ∈ R n×n
are fixed (deterministic) design matrices. As a shorthand, define the linear map ω :

R n×n → R d as ω(X) = ( Ω 1 , X , • • • , Ω d , X ) .
Shortly, we have y = ω(X 0 ) + ε.

For S ⊂ R n×n some convex admissible set, we consider the estimation procedure

X = arg min X∈S {ℓ(ω(X), y) + τ X * + γ X 1 } (3.2)
and in the following, we study properties of the estimator X. The underlying assumption in this work is that the unknown matrix to be recovered is simultaneously sparse and low-rank. This is the case of matrices that are block-diagonal in a certain basis, but they are not the only sparse low-rank matrices.

In the following we provide examples of applications where target matrices to be estimated are sparse low-rank. We write the loss functions desirable for each case and provide basic properties of each. In the next section we study the case of the quadratic loss by providing oracle bounds.

Multivariate regression

The problem (3.1) is a multivariate regression problem (see [START_REF] Breiman | Predicting multivariate responses in multiple linear regression[END_REF]). Assuming that the multiple outputs share a common structure, i.e. they can be expressed using the same small set of linear combinations of initial features, then the regression matrix X 0 to be estimated is lowrank. Moreover, in a setup where the number of descriptors is potentially high we expect most of them to be irrelevant for the prediction task, hence the same low-rank matrix X 0 is also sparse. A special case of regression involving sparse low-rank matrices arises when predicting linear features of a graph sequence in an autoregressive fashion. Such a setting has been considered in problems such as customer valuation and network effect measurement [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF] and link prediction in a temporal setting (e.g. [START_REF] Richard | Link discovery using graph feature tracking[END_REF], [START_REF] Vu | Continuous-time regression models for longitudinal networks[END_REF]), where the observations consist of a set of linear features of the graph (such as the degrees) and the prior knowledge on the domain is that the graph to be predicted contains highly connected communities. The community structure of the graph is translated by a joint rank and sparsity constraint, that must be imposed on the linear regression.

In section 3.3 we provide oracle bounds that guarantee the quality of the estimator defined by (3.2) in the case of least-squares regression.

Covariance matrix estimation

Many interesting examples arise from the simplest denoising case, where the design matrices correspond to the canonical basis of R n×n , or equivalently ω(X) = X. In this setting the matrix A represents a noisy estimate of the true covariance matrix. It can be obtained for instance by estimating the sample covariance with few (less than n) observations which may be noisy themselves. The search space in this case is S = X + n , the cone of positive semidefinite matrices. For this application the noise matrix can be a full matrix, theregore the squared norm ℓ(X, A) = X -A 2 F is often used in matrix covariance estimation applications. Corollaries of oracle bounds presented in section 3.3 apply to the current case ensuring the quality of the estimated covariance matrix.

Least absolute deviation and graph denoising

We now turn to the case of sparse noise, in the following sense: for each i = 1, . . . , d, ε i is zero with probability p close to one. Hence, the expected number of nonzero entries of ε is d p. In a regression setup such a situation is often treated by using an ℓ 1 loss ℓ(X, y) = ω(X)y 1 , and is called least absolute deviation (see [START_REF] Bassett | Asymptotic theory of least absolute error regression[END_REF]), which has been studied in a penalized framework by [START_REF] Gao | Asymptotic analysis of high-dimensional lad regression with lasso[END_REF]. In the present work, for shortness and as our principal goal is to study the virtues of the mixed penalty, we study only the simple denoising case ω(A) = A.

A problem which arises in many relational databases is the presence of irrelevant relations between objects, while many relevant links are often missing. Representing these data as graphs, both missing and undesirable links can be modeled as noisy edges. These noisy links can be due to measurement bias or to the specific interpretation on what type of interaction should be encoded through these links. In the case of social networks, a link is supposed to indicate friendship between two members. However, relations may be declared after randomly meeting at a party, while close friends or family may see no reason to establish digital friendship through the social network. A similar example comes from purchase data on e-commerce websites where users may purchase items for their relatives and produce irrelevant links in the bipartite graph of users and products. Recommender systems based on these relational data face the challenge of removing those noisy links, since they are meaningless in terms of affinity between members of the network. Recovering structure in large noisy graphs presents a major interest both in collaborative filtering and exploratory web mining.

Our objective is to explore the following question: "how can one obtain a sparse and lowrank matrix by adding or removing as few edges as possible from a given sparse adjacency matrix? " In the presence of a sparse unstructured noise, the problem we are interested in is to estimate a matrix simultaneously sparse and low-rank by adding or removing only a few edges from an observed matrix A = S 0 + ε. This can be formulated as follows: min X X -A 0 subject to X 0 ≤ q and rank(X) ≤ r . This is, however, a hard combinatiorial problem. Instead, we consider the relaxation obtained by replacing . 0 by its tightest convex upper bound, the matrix ℓ 1 -norm. Similarly, we relax the rank constraint using the nuclear norm. This motivates the following penalized estimator, for some parameters τ, γ > 0 (set by cross-validation in practice) :

X ∈ arg min X∈S { X -A 1 + τ X * + γ X 1 } . (3.3)
The objective function is convex but neither smooth, nor strictly convex, due to the presence of the ℓ 1 -loss. This particular loss has been used in regression under the name of least absolute deviations as an alternative to ordinary least squares [START_REF] Bassett | Asymptotic theory of least absolute error regression[END_REF], or in a high-dimensional setting using ℓ 1 -norm penalization [START_REF] Gao | Asymptotic analysis of high-dimensional lad regression with lasso[END_REF]. The introduction of the ℓ 1 -loss is generally motivated by the presence of heavy-tailed noise and outliers, as it is more robust. The combination of ℓ 1norm and nuclear norm penalty on the same matrix to obtain a simultaneously sparse and low-rank estimator has been recently proposed in [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF] using a squared loss. As it turns out, the non differentiable ℓ 1 loss is significantly more challenging to tackle, but also much more suited to the problem of inferring links in graphs in the presence of sparse noise.

This approach encompasses various other methods. In particular, RPCA corresponds to the problem min

X,Y ∈R n×n { Y 1 + τ X * } s.t. X + Y = A
which can be rewritten using

X = A -Y min X∈R n×n { X -A 1 + τ X * } .
Hence, RPCA corresponds to γ = 0 in our formulation. Another standard approach for estimating the low-rank structure of the adjacency matrix consists of minimizing the convex objective with squared loss min

X∈R n×n { X -A 2 F + τ X * } .

Recovering a partially observed graph

Assume the matrix A is the adjacency matrix of a partially observed graph: entries are 0 for both not-existing and undiscovered links. The search space is unrestricted as before and the matrix X contains the scores for link prediction; the ideal loss function is the empirical average of the zero-one loss for each coefficient

ℓ E (X, A) = 1 |E| (i,j)∈E 1{(A i,j -1/2) • X i,j ≤ 0} ,
where E is the set of edges in A. However, as in classification, practical algorithms should use a convex surrogate (e.g., the hinge loss). Given a subset E of observed edges from a graph adjacency matrix A ∈ {0, 1} n×n , we set out to predict unobserved links by finding a sparse rank r predictor X ∈ R n×n with small zero-one loss

ℓ(X, A) = 1 n 2 (i,j)∈{1,...,n} 2 1{(A i,j -1/2) • X i,j ≤ 0}
by minimizing the empirical zero-one loss ℓ E (X, A). The objective of a generalization bound is to relate ℓ(X, A) with ℓ E (X, A). In the case of the sole rank constraint, Srebro [START_REF] Srebro | Learning with matrix factorizations[END_REF] remarked that all low-rank matrices with the same sign pattern are equivalent in terms of loss and applied a standard argument for generalization in classes of finite cardinality. In the work of Srebro, an argument is used to upper bound the number of distinct sign configurations for predictors of rank r s lr (n, r) = |{sign(X) | X ∈ R n×n , rank(X) ≤ r}| leading to the following generalization performance: for δ > 0, A ∈ {0, 1} n×n and with probability 1δ over choosing a subset E of entries in {1, . . . , n} 2 uniformly among all subsets of |E| entries, we have for any matrix X of rank at most r and ∆(n, r)

= 8en r 2nr ℓ(X, A) < ℓ E (X, A) + log ∆(n, r) -log δ 2|E| . (3.4)
We consider the class of sparse rank r predictors

M(n, r, s) = {U V T | U, V ∈ R n×r , ||U || 0 + ||V || 0 ≤ s}
and let s splr (n, r, s) be the number of sign configurations for the set M(n, r, s). By upper bounding the number of sign configurations for a fixed sparsity pattern in (U, V ) using an argument similar to [START_REF] Srebro | Learning with matrix factorizations[END_REF], a union bound gives

ℓ(X, A) < ℓ E (X, A) + log Γ(n, r, s) -log δ 2|E| . (3.5)
In general, bound (3.5) is tighter than (3.4) for sufficiently large values of n as shown in the next proposition. The two bounds coincide when s = 2nr, that is, when (U, V ) is dense and there is no sparsity constraint.

Proposition 1 For r n = nβ with β ∈]0, 1] and

s n = nα with α ≤ 2β, ∆(n, r n ) Γ(n, r n , s n ) = Ω 8en(βn -α) (βn) 2 2n 2 β ,
where Ω denotes bounded bellow asymptotically.

The result follows from the application of Stirling's formula. The expression on the right hand side of the equation, giving the asymptotic ratio between the two bounds diverges when n goes to infinity. In particular it shows that the bound for sparse low-rank matrices 3.5 is asymptotically tighter than the bound 3.4 for low-rank (potentially full) matrices.

Oracle bounds for the estimation of sparse low-rank matrices with quadratic loss

In this section, we present theoretical guarantees for the estimation of sparse low-rank matrices in the presence of noises and using quadratic loss.

In this section we consider the case of linear regression. That is, we consider

y = ω(X 0 ) + ǫ ∈ R d
with ǫ ∈ R d having i.i.d zero mean entries and ω : R n×n → R d linear. Define the objective

L(X) = 1 d ω(X) -y 2 2 + τ X * + γ X 1 , (3.6)
consider the following estimation procedure

X = arg min X∈S L(X) . (3.7)
We aim at studying properties of the estimator X. Define the matrix M = d i=1 ǫ i Ω i . We begin by stating a simple oracle inequality.

Proposition 3 Fix a real number α ∈ [0, 1]. Under assumption τ ≥ 2α d M op , and γ ≥ 2(1-α) d M ∞ if X = arg min S L, we have 1 d ω( X -X 0 ) 2 2 ≤ inf X∈S 1 d ω(X -X 0 ) 2 2 + 2τ X * + 2γ X 1 (3.8)
and in particular, when X 0 ∈ S, taking X = X 0 gives

1 d ω( X -X 0 ) 2 2 ≤ 2τ X 0 * + 2γ X 0 1 . (3.9)
In order to state recovery results on X directly, we introduce an assumption that is similar to the restricted isometry property introduced in [START_REF] Candès | Decoding by linear programming[END_REF].

Assumption 1 (Restricted Isometry Propery)

There exists µ > 0 such that for any X 1 , X 2 ∈ S

1 d ω(X 1 -X 2 ) 2 2 ≥ µ -2 X 1 -X 2 2 F .
(3.10)

We can now state the following result, where we set c 0 = √ 2+1 2 . Proposition 4 Under assumption 1, for any real number α ∈ [0, 1], for constants τ ≥ 2α d M op , and

γ ≥ 2(1-α) d M ∞ , we have 1 d ω( X -X 0 ) 2 2 ≤ inf X∈S 1 d ω(X -X 0 ) 2 2 + µ 2 d c 0 τ rank(X) + γ X 0 2 (3.11)
and in particular if X 0 ∈ S, taking X = X 0 results in

1 d ω( X -X 0 ) 2 2 ≤ min µ 2 d c 0 τ rank(X 0 ) + γ X 0 0 2 , 2τ X 0 * + 2γ X 0 1 (3.12)
Notice in particular that the bound 3.11 containing the term c 0 τ rank(X 0 ) + γ X 0 0 shows that minimizing the convex surrogate objective L leads to the optimal solution obtained by minimizing the nonconvex (and exponentially hard to optimize) objective

1 d ω(X -X 0 ) 2 2 + µ 2 d c 0 τ rank(X) + γ X 0 2 .
An immediate corollary of Propositions 3 and 4 is the following.

Corollary 1 Under the assumptions of Proposition 4,

X -X 0 2 F ≤ µ 2 min µ 2 d c 0 τ rank(X 0 ) + γ X 0 0 2 , 2τ X 0 * + 2γ X 0 1 . (3.13)
In the particular case of denoising ω(X) = X, the constant µ equals n = √ d, and therefore we have the following result for the minimizer of J : S → X -A 2 F + τ X * + γ X 1 (pay attention to the change of normalization compared to L).

Corollary 2 (Denoising) Let X be the minimizer of J , and α ∈ [0, 1] a real number, for constants τ ≥ 2α M op , and γ ≥ 2(1α) M ∞ , we have

X -X 0 2 F ≤ min c 0 τ rank(X 0 ) + γ X 0 0 2 , 2τ X 0 * + 2γ X 0 1 (3.14)
Assumption 1 can be restrictive. In order to state bounds that hold under less stringent conditions, we introduce the following assumption similar to restricted eigenvalue conditions related to those introduced in [BRT09] and [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] that holds in more realistic scenarios than Assumption 1. We first define a cone that contains matrices such that projections onto the singular spaces and onto the sparsity pattern of a matrix S dominate the orthogonal projections, with a trade-off ratio β between the two constraints.

Definition 16 (Cone of restriction) For a matrix X ∈ S ⊂ R n×n and any β, κ > 0, letting P X be the projection onto the singular space of X and Θ X = sign(X) the sign pattern of X,

C(X, κ, β) = B ∈ S | P ⊥ X (B) * + β Θ ⊥ X • B 1 ≤ κ ( P X (B) * + β Θ X • B 1 )
is the cone of dominant coordinates.

We now relax Assumption 1 by only requiring that it holds over the cone.

Definition 17 (Restricted eigenvalue constant) For a linear map ω :

R n×n → R d , µ κ,β (X) = inf µ ′ > 0 : max( P X (B) F , Θ X • B F ) ≤ µ ′ d ω(B) 2 ∀B ∈ C(X, κ, β) .
Assumption 2 (Restricted eigenvalue) The map ω : R n×n → R d verifies RE assumption κ, β at point X if µ κ,β (X) is bounded away from zero.

We have the following proposition, that holds true under weaker assumptions on ω.

Proposition 5 Let α ∈ [0, 1]. With τ > 3α d M op and γ > 3(1-α) d M ∞ , we have under Assumption 2 1 d ω( X -X 0 ) 2 2 ≤ inf X∈S 1 d ω(X -X 0 ) 2 2 + µ 5, γ τ (X) 2 d c 1 τ rank(X) + c 2 γ X 0 2 (3.15)
and, in particular,

ω( X -X 0 ) 2 2 ≤ µ 5, γ τ (X 0 ) 2 c 1 τ rank(X 0 ) + c 2 γ X 0 0 2 with c 1 = 3+2 √ 2 6
, c 2 = 5 6 .

The techniques used in the proof (see the Appendix) are very similar to those introduced in [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF]. Note that the upper bound interpolates between the results known for trace-norm penalization and the Lasso. In fact, for α = 0, τ can be set to zero, and we get a sharp bound for the Lasso, while the trace-norm regression bounds of [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] are obtained for α = 1.

In the next Theorem 8, we obtain convergence rates for the procedure 3.7 by combining Proposition 5 with controls on the noise process, by using techniques introduced in [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF], see Proposition 8 in appendix. For this, we need to make the assumption that noise entries ǫ i s are independent, zero-mean and subgaussian i.e. there exists σ > 0 such that for any positive real number λ > 0,

Ee λǫ i ≤ e σ 2 λ 2 /2 .
We introduce two quantities of interest:

v 2 Ω,op = 1 d d j=1 Ω ⊤ j Ω j op ∨ 1 d d j=1 Ω j Ω ⊤ j op , v 2 Ω,∞ = 1 d d j=1 Ω j • Ω j ∞
which are the (observable) variance terms that naturally appear in the controls of the noise processes.

Theorem 8 Consider the procedure X given by (3.7) with smoothing parameters given by

τ = 3ασv Ω,op 2(t + log(2n)) d , γ = 3(1 -α)σv Ω,∞ 2(t + 2 log n) d
for some α ∈ (0, 1) and fix a confidence level t > 0. Then, we have

ω( X -X 0 ) 2 2 ≤ inf X∈S ω(X -X 0 ) 2 2 + 25µ(X) 2 rank(X)α 2 σ 2 v 2 Ω,op 2(t + log(2n)) d + 25µ(X) 2 X 0 (1 -α) 2 σ 2 v 2 Ω,∞ 2(t + 2 log n) d
with a probability larger than 1 -3e -t , where µ(X) = µ 5, γ τ (X).

Optimization methods

Preliminary results on proximal operators

We now present how to solve the optimization problem with mixed penalties presented in subsection 2. We consider a loss function ℓ(X, A) convex in X. When ℓ is also differentiable in X, we assume that its gradient is Lipschitz with constant L and can be efficiently computed. This is, in particular, the case for the squared Frobenius norm previously mentioned and for other classical choices such as the hinge loss.

We encode the presence of a constraint set S using the indicator function 1 S (X) that is zero when X ∈ S and +∞ otherwise, leading to

X = arg min X∈R n×n {ℓ(X, A) + γ||X|| * + τ ||X|| 1 + 1 S (X)} .
This formulation involves a sum of a convex differentiable loss and of convex nondifferentiable regularizers which renders the problem non-trivial. A string of algorithms have been developed for the case where the optimal solution is easy to compute when each regularizer is considered in isolation. Formally, this corresponds to cases where the proximal operator is easy to compute for each regularizer taken separately. We first recall the definition of a proximal operator (Xee [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] for a broad overview of proximal methods).

Definition 18 Let f : R d → R a convex function, then the proximal operator of f is defined for x ∈ R d by

prox f (x) = arg min z∈R d f (z) + 1 2 z -x 2 2 .
The proximal operator of the indicator function is simply the projection onto S, which justifies the alternate denomination of generalized projection operator for prox R . The proximal operator for the trace norm is given by the shrinkage operation as follows [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

If Z = U diag(σ 1 , • • • , σ n )V T is the singular value decomposition of Z, SHR τ (Z) := prox τ ||.|| * (Z) = U diag((σ i -τ ) + ) i V T .
Similarly, the proximal operator for the ℓ 1 -norm is the soft thresholding operator

ST γ (Z) := prox γ||.|| 1 = sgn(Z) • (|Z| -γ) + .
We now generalize this last result to more complex combinations of ℓ 1 -norms.

Definition 19 Let α 1 < . . . < α k a sorted sequence of distinct real numbers and (β 1 , . . . , β k ) nonnegative numbers. We note ξ α, β (x) the proximal operator of the function h defined for x ∈ R by

h(x, α, β) = k i=1 β i |x -α i | .
We provide in Lemma 1 below a closed-form expression for ξ α, β . In particular, ξ α, β is continuous and nondecreasing piecewise affine which extends classical results on the proximal operator of the ℓ 1 -norm.

Lemma 1 (Proximal operator of h) Let α 1 < . . . < α k a sorted sequence of distinct real numbers, (β 1 , . . . , β k ) nonnegative numbers, and h defined for x ∈ R as in Definition 19 Noting

B i = i j=1 β j - k j=i+1 β j , for i ∈ {0, . . . , k} , the proximal operator of h at x ∈ R is ξ α,β (x) =          x -B k if x > B k + α k x -B 0 if x < B 0 + α 1 α i if there exists i s.t. x ∈ [B i-1 + α i ; B i + α i ] x -B i if there exists i s.t. x ∈ [B i + α i ; B i + α i+1 ] . (3.16) Note that since α 1 + B 0 < α 1 + B 1 < α 2 + B 1 < α 2 + B 2 < .
. ., there actually exists a single case that applies.

Splitting methods

The family of Forward-Backward splitting methods are iterative algorithms applicable when there is only one nondifferentiable regularizer. These methods alternately take a gradient step over the Lipschitz-differentiable ℓ and and a proximal step over the convex and possibly nonsmooth R, leading to updates of the form

X k+1 = prox θR (X k -θ grad X ℓ(X, A)) .
In particular, this corresponds to projected gradient descent when R is the indicator function of a convex set.

Douglas-Rachford splitting

We now turn to the question of minimizing

L(X) = X -A 1 + τ X * + γ X 1 .
This can be carried out using proximal methods, which are popular methods for nonsmooth convex analysis [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Douglas-Rachford splitting allows to minimize a sum of two nonsmooth convex functions as long as the proximal operator of each of them can be computed [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Extensions of the method have been proposed to extend the algorithm to more than two functions but usually at the cost of introducing extra variables [START_REF] Spingarn | Partial inverse of a monotone operator[END_REF], which results in a higher computational and memory complexity. Instead of using the proximal of each of the three components of the objective, we propose to apply Douglas-Rachford splitting using the decomposition

L(X) = f (X) + g(X) with f (X) = X -A 1 + γ X 1 , g(X) = τ X * .
Since f (X) is separable with respect to the coefficients (X i,j ), the proximal operator of f can be computed coordinate-by-coordinate as prox f (X) = (prox f (X i,j )) i,j ∈ R n×n . For each coordinate, this reduces to computing the proximal operator of x → γ|x| + |x -a| which can easily be computed using Lemma 1. In particular, this corresponds to the proximal operator of ξ α, β with α = (0, a), β = (γ, 1) when a ≥ 0, and to α = (a, 0), β = (1, γ) when a < 0 .

For convenience in the exposition of the algorithm, we define for X ∈ R n×n the reflexive proximal of a function h as rprox h (X) = 2 prox h (X) -X and P S denotes orthogonal projection onto the convex set S. The pseudocode oh the algorithm based on the Douglas-Rachford scheme is shown in Algorithm 1.

Algorithm 1 Douglas-Rachford Scheme for Graph Denoising Input A, τ, γ > 0, θ ∈ (0, 2) Initialization:

Z (0) = A for k = 1, 2, . . . , K do Z (k) = 1 -θ 2 Z (k-1) + θ 2 rprox τ . * • rprox f (Z (k-1) ) end for X = prox τ . * (Z (K) )

Generalized Forward-Backward splitting

A generalization of the two setups (Lipschitz-differentiable + easy proximal) and (easy proximal + easy proximal) has been recently proposed in [START_REF] Raguet | Generalized forward-backward splitting[END_REF] under the name of Generalized Forward-Backward, which we specialize to our problem in Algorithm 2. The proximal operators are applied in parallel, and the resulting (Z 1 , Z 2 , Z 3 ) is projected onto the constraint that Z 1 = Z 2 = Z 3 which is given by the mean. The auxiliary variable Z 3 can be simply dropped when S = R n×n . The algorithm converges under very mild conditions when the step size θ is smaller than 2 L .

Algorithm 2 Generalized Forward-Backward Initialize X, Z 1 , Z 2 , Z 3 = A, q = 3 repeat Compute G = ∇ X ℓ(X, A). Compute Z 1 = prox qθτ ||.|| * (2X -Z 1 -θG) Compute Z 2 = prox qθγ||.|| 1 (2X -Z 2 -θG) Compute Z 3 = P S (2X -Z 3 -θG) Set X = 1 q q k=1 Z k until convergence return X

Incremental Proximal Descent

Although Algorithm 2 performs well in practice, the O(n 2 ) memory footprint with a large leading constant due to the parallel updates can be a drawback in some cases. As a consequence, we mention a matching serial algorithm (Algorithm 3) introduced in [Ber11] that has a flavor similar to multi-pass stochastic gradient descent, and has a convergence rate in O( 1 ǫ 2 ) to an approximate solution in a radius of ǫ > 0. This convergence rate is slow compared to the standard proximal gradient descent that converges in O( 1 ǫ ) , or O( 1 √ ǫ ) for the accelerated version. We present here a version where updates are performed according to a cyclic order, although random selection of the order of the updates is also possible.

Algorithm 3 Incremental Proximal Descent Initialize X = A repeat Set X = X -θ∇ X ℓ(X, A) Set X = prox θτ ||.|| * (X) Set X = prox θγ||.|| 1 (X) Set X = P S (X) until convergence return S

PSD constraint

For any positive semidefinite matrix, we have ||Z|| * = Tr(Z). The simple form of the trace norm allows to take into account the positive semidefinite constraint at no additional cost, as the shrinkage operation and the projection onto the convex cone of positive semidefinite matrices can be combined into a single operation.

Lemma 2 For τ ≥ 0 and X ∈ R n×n , if P S + n denotes the projection onto the PSD cone,

prox τ ||.|| * +1 S + n (X) = arg min Z 0 1 2 ||Z -X|| 2 F + τ ||Z|| * = P S + n (X -τ I n ) .
The solution of this problem is given by soft-thresholding the singular values of A and is known as the shrinkage operator (e.g., [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]). This is similar to the best rank-r approximation of a matrix in Eucliean norm that is obtained through hard-thresholding of the singular values.

Fast methods using factorization-based updates

The method proposed in Algorithm 1 is slow due to full SVD computations at every iteration in order to compute the proximal operator of the nuclear norm. Given the size of real data sets, there is a real need for exploring more scalable algorithms. It is well known [START_REF] Srebro | Maximum-margin matrix factorization[END_REF] that the nuclear norm can be characterized as

X * = 1 2 min X=U V ⊤ , U,V ∈R n×n U 2 F + V 2 F .
Motivated by this characterization, number of approaches have been developed [START_REF] Srebro | Maximum-margin matrix factorization[END_REF] for matrix estimation subject to rank constraint without explicitly penalizing the objective with the nuclear norm. Most of these approaches optimize for factorizations of the form X = U V ⊤ and address the case of the squared loss

U V T -A 2 F with penalization proportional to U 2 F + V 2 F .
The motivation for such a factorization rather the convex methods is 1. Good empirical results. See for instance [START_REF] Koren | Factorization meets the neighborhood: a multifaceted collaborative filtering model[END_REF].

2. Attracting scalability. In fact, as clearly argued in [START_REF] Recht | Parallel stochastic gradient algorithms for large-scale matrix completion[END_REF], optimizing factors leads to scalable methods both in terms of memory requirement and computational complexity, and may be parallelized easily.

3.

Interpretability. The factors U and V can be interpreted as the latent features describing individuals nodes. Nonnegative matrix factorization (see for instance [START_REF] Kim | Sparse nonnegative matrix factorization for clustering[END_REF]) encourages further interpretability by imposing nonnegativity to the entries of the factors.

Nevertheless it should be pointed out that this method leads to nonconvex objectives. Alternating minimization schemes are usually developed, where one alternates minimization over U with minimization over V . This method has no theoretical guarantee to converge to a global optimum, but empirical success has been reported in many experiments.

Sparse factorizations

Sparse matrix factorizations have been proposed in sparse principal component analysis (SPCA) to extend the classical principal component analysis and find sparse directions that maximize the variance of the projection by penalizing the sparsity of the principal components [START_REF] Zou | Sparse principal component analysis[END_REF].

Similarly, [START_REF] Kim | Sparse nonnegative matrix factorization for clustering[END_REF] have proposed in the context of nonnegative matrix factorization to encourage sparsity of the factors U and V by penalizing U 1 and V 1 instead of directly enforcing the sparsity of X by penalizing by X 1 . A simple baseline for graph denoising is thus the following squared-loss objective for sparse factorization:

L 2 (U, V ) = U V ⊤ -A 2 F + τ 2 U 2 F + V 2 F + γ ( U 1 + V 1 ) . (3.17)
As we already outlined, the squared loss is not adequate for sparse noise. We propose to use a factorization method along these lines to find an approximate solution to the optimization problem defined at Equation (3.3). The surrogate objective we consider is hence

L U V (U, V ) = U V ⊤ -A 1 + τ 2 U 2 F + V 2 F + γ ( U 1 + V 1 ) .
The penalty on U and V can also be interpreted as the elastic-net penalty [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] over the factors U and V . Alternating minimization of L U V is more complex than with traditional alternating minimization approaches in the presence of the squared loss, and we will now present a procedure to optimize L U V in U, V ∈ R n×r for some fixed r.

Alternating minimization with rank one updates

Recently, [DHM12, SSGO11] have suggested to recover low-rank matrices by successively finding rank-one updates minimizing a given loss, while [START_REF] D'aspremont | A direct formulation for sparse pca using semidefinite programming[END_REF] introduced a method for SPCA using successive sparse rank-one updates. In a similar fashion, we propose to use successive sparse rank one updates minimizing the following objective in (u, v) ∈ R n :

ℓ(u, v) = uv ⊤ -A 1 + τ 2 ( u 2 2 + v 2 2 ) + γ( u 1 + v 1 ) .
We then deflate the original A by uv T and repeat the procedure r times where r is a fixed target maximum rank. Note that for τ > 0 the functional ℓ is strictly convex in each variable u or v, the other one being kept fixed. Although ℓ is also jointly convex in (u, v) for τ > 2, the joint minimization is hard to carry out, and we suggest to perform a regularized altering minimization in u and v instead. More precisely, we consider sequences (u (k) , v (k) ) k defined by

u (k+1) = arg min u∈R n uv (k) ⊤ -A 1 + γ u 1 + 1 2θ u -(1 -θτ )u (k) 2 2 v (k+1) = arg min v∈R n u (k+1) v ⊤ -A 1 + γ v 1 + 1 2θ v -(1 -θτ )v (k) 2 2 .
for some θ > 0. Each update corresponds to one step of proximal gradient descent [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] with respect to each variable, while the other one is kept fixed. Notice that for θ = 1/τ , there is no extra regularization and this is simply a standard alternate minimization of ℓ. In practice, we choose θ ∈ (0, 1/τ ) as this exhibited empirical superiority to the standard alternate minimization in numerical experiments, as shown in Table 3.3. In addition to this, we initialize the alternate minimization using the left/right singular vectors corresponding to the highest singular value. The full procedure is detailed in Algorithm 4.

In the following, we sketch the computations for the minimization with respect to u and point out that equivalent formulas can be obtained by swapping u and v and replacing A by

Algorithm 4 Factorization-based Graph Denoising with Rank-1 updates

Input A, τ, γ > 0, θ < 1/τ, r ≤ n, T : number of iterations Initialization:

Z (0) = 0 for k = 1, 2, . . . , r do Initialize u (k) , v (k) for i = 1, • • • , T do u (k) = arg min u∈R n ℓ(u, v (k) ) v (k) = arg min v∈R n ℓ(u (k) , v) end for Z (k) = Z (k-1) + u (k) v (k) T A = A -u (k) v (k) T end for X = Z (r)
A ⊤ . For computing u (k+1) , first notice that

uv (k) ⊤ -A 1 + γ u 1 + 1 2θ u -(1 -θτ )u (k) 2 2 = n i=1   n j=1 |u i v (k) j -A i,j | + γ|u i | + 1 2θ (u i -(1 -θτ )u (k) i ) 2   = n i=1    j | v (k) j =0 |v (k) j ||u i - A i,j v (k) j | + γ|u i | + 1 2θ (u i -(1 -θτ )u (k) i ) 2    + C = 1 θ n i=1    j | v (k) j =0 θ|v (k) j ||u i - A i,j v (k) j | + θγ|u i | + 1 2 (u i -(1 -θτ )u (k) i ) 2    + C . The term C = j | v (k) j =0 |A i,j |
does not depend on u. We have a sum (over i) of terms of the form

k β k |u i -α k | + 1 2θ (u i -δ) 2
, and we can thus use Lemma 1 to find the minimizer of this expression with respect to u i . The vector u (k) can hence be updated by applying Lemma 1 for each i at the point

x = (1 -θτ )u (k) i
and using

{α j } = {A i,j /v (k) j , j s.t. v (k) j = 0} ∪ {0}, {β j } = {θ|v (k) j |, j s.t. v (k) j = 0} ∪ {θγ}
where we used a set notation as in practice, for each i, the coefficients of α i and β i must be sorted by increasing value of A i,j /v (k) j , that is, we must find an injection

π : {1, • • • , q} → {1, • • • , n} by sorting the q non-zeros elements of v (k) such that A i,π(1) /v (k) π(1) < • • • < A i,π(q) /v (k)
π(q) . Each iteration thus has a cost of O( u 0 v 0 log v 0 + v 0 u 0 log u 0 ), which leads to an efficient procedure for sparse u and v in O r T u 0 v 0 log v 0 + v 0 u 0 log u 0 .

The terms . 0 should be read as the average sparsity over the iterations. The runtime of the algorithm in practice highly depends on the sparsity of the initial vectors u (k) and v (k) . It can considerably be reduced by initializing for each k, u (k) and v (k) with sparse vectors. In our experiments we initialized them with thresholded first singular vectors or with normalized hard-thresholded degrees of the nodes. Further experiments are needed to understand the dependency of the solution on the initialization.

Remark 1 (Penalizing the product u v rather than the sum u + v ) It is known [START_REF] Jameson | Summing and Nuclear Norms in Banach Space Theory[END_REF] that the function of X defined by the value

inf ∞ i=1 u i L v i R u i , v i ∈ R n s.t. X = ∞ i=1 u i v ⊤ i
for any pair of norms . L and . R is a norm. If both the norms are the ℓ 1 norm, then we obtain the standard element-wise ℓ 1 norm of the matrix:

inf ∞ i=1 u i 1 v i 1 u i , v i ∈ R n s.t. X = ∞ i=1 u i v ⊤ i = i,j |X i,j | .
In order to use this observation, rather than minimizing ℓ, one may rather want to minimize the nonconvex functional

ℓ(u, v) = uv ⊤ -A 1 + τ 2 ( u 2 2 + v 2 2 ) + γ u 1 v 1 .
Note that the same Algorithm 4 applies by simply replacing the inside loop iterations to

u (k) = arg min u∈R n ℓ(u, v (k) ) and v (k) = arg min v∈R n ℓ(u (k) , v) .
The minimizer of ℓ can be obtained from the minimizer of ℓ (see Lemma 1) by replacing γ by γ v 1 . As we are using only r < ∞ factors the optimized functional does not have a convex equivalent in the same way as the trace norm is the equivalent of using sum of squares of Frobenius norms. This does not cure the nonconvexity of the function

µ(x) = inf u 1 + v 1 x = uv ⊤ .
To see why it is nonconvex, notice that for n = 1 one gets the square root |x|. In the numerical experiments we have not yet seen a significant difference between the results of the two methods. An issue with these methods, due to nonconvexity, is the existence of non-optimal stationary points such as (0, 0).

Remark 2 (Rank r updates) Instead of the objective ℓ(u, v) of two vector variables, one can use a similar function L U V (U, V ) of matrices U, V ∈ R n×r . Minimizing following each single coordinate variable u i,j requires the same tools as minimizing ℓ(u, v), so one can perform coordinate descent on L U V as well. The comparative study of several competing methods for sparse matrix factorization is left to future work.

Numerical experiments

We present numerical experiments to highlight the benefits of our method. For efficiency reasons, we use the incremental proximal descent algorithm (Algorithm 3).

Covariance estimation and graph denoising with Frobenius norm loss

Synthetic data

Covariance matrix estimation. We draw N vectors x i ∼ N (0, Σ) for a block diagonal covariance matrix Σ ∈ R n×n . We use r blocks of random sizes and of the form vv ⊤ where the entries of v are drawn i.i.d. from the uniform distribution on [-1, 1]. Finally, we add gaussian noise N (0, σ 2 ) on each entry. In our experiments r = 5, N = 20, n = 100, σ = 0.6. We apply our method (SPLR), as well as trace norm regularization (LR) and ℓ 1 norm regularization (SP) to the empirical covariance matrix, and report average results over ten runs. Figure 3.1 shows the RMSE normalized by the norm of Σ for different values of τ and γ. Note that the effect of the mixed penalty is visible as the minimum RMSE is reached inside the (τ, γ) region. We perform, on the same data, separate cross-validations on (τ, γ) for SPLR, on τ for LR and on γ for SP. We show in Figure 3.2 the supports recovered by each algorithm, the output matrix of LR being thresholded in absolute value. The support recovery demonstrates how our approach discovers the underlying patterns despite the noise and the small number of observations. low-rank representation. Representing the data as a weighted graph, we filter to retain only the 10% of all 4394 proteins that exhibit the most interactions as measured by weighted degree. We corrupt 10% of entries of the adjacency matrix selected uniformly at random by uniform noise in [0, η]. Parameters are selected by cross-validation and algorithms are evaluated using mean RMSE between estimated and original adjacency matrices over 25 runs. RMSE scores are shown in Table 3.1 and show the empirical superiority of our approach (SPLR). Social Networks. We have performed experiments with the Facebook100 data set analyzed by [START_REF] Traud | Social structure of facebook networks[END_REF]. The data set comprises all friendship relations between students affiliated to a specific university, for a selection of one hundred universities. We select a single university with 41554 users and filter as in the previous case to keep only the 10% users with highest degrees. In this case, entries are corrupted by impulse noise: a fixed fraction σ of randomly chosen edges are flipped, thus introducing noisy friendship relations and masking some existing relations. The task is to discover the noisy relations and recover masked relations. We compare our method to standard baselines in link prediction [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF]. Nearest Neighbors (NN) relies on the number of common friends between each pair of users, which is given by A 2 when A is the noisy graph adjacency matrix. Katz's coefficient connects a pair of nodes according to a score based on the number of paths connecting them, emphasizing short paths. Results are reported in Table 3.2 using the area under the ROC curve (AUC). SPLR outperforms LR but also NN and Katz which do not directly seek a low-rank representation.

Graph denoising with ℓ 1 norm loss

Regularization paths

The model parameters τ, γ can be tuned using cross-validation in a supervised setting. In an exploratory perspective, one may rather look at the sparsity and rank of the solution as a function of τ and γ. In particular, exploring the solutions of the optimization problem defined by Equation (3.3) for various parameters clearly shows the benefit of our method compared to RPCA. In the left panel of Figure 3.3, we show the regularization path of RPCA for various values of the parameter τ . We insist on the fact that using this unique parameter does not allow to achieve different levels of sparsity and rank independently. The right panel of Figure 3.3 shows on an example that by choosing nonzero values for both τ and γ, we achieve better support recovery: top left shows the support of A 0 in blue and of the noise in red, while bottom left corresponds to RPCA and bottom right to Convex Graph Denoising. We point out that for τ = 0, nonzero values of γ do not lead to any interesting solution (top right) as separating sparse noise from sparse signal without any additional prior is meaningless.

Phase transition

We analyzed the influence of the characteristics of the noise and of A 0 on the performance of Factorization-based Graph Denoising. For different values of the sparsity of the noise, of the sparsity of A 0 and of the rank of A 0 , we report the mean relative ℓ 1 error after 50 iterations of cross-validation on the parameters of the algorithm. and the sparsity of the noise vary. Unsurprisingly, the recovery is all the more accurate than the noise is sparse and A 0 low-rank. Level sets take the form of diagonal lines which indicate the following tradeoff: to maintain an equal recovery error when the density of the noise increases, the rank of A 0 must decrease, or equivalently, A 0 must become more structured.

Empirical evaluation

Synthetic data. We generate synthetic data A = A 0 +N , where A 0 = U 0 V ⊤ 0 . The two matrices V 0 ∈ R n×r and N ∈ R n×n have sparsity patterns chosen uniformly at random and their entries are drawn i.i.d. from the uniform distribution U (-1, 1). The matrix U 0 has a one on each row at a randomly chosen position. We point out that using this construction, A 0 is both sparse and low-rank, N is sparse and has high rank, and both have entries of with the same order of magnitude.

In our experiments, the number of nodes is n = 100, r = 6, N 0 = 200, A 0 0 = 350. We performed 10-fold cross validation for choosing the values of the parameters τ and γ, which not surprisingly turn out to be the very close for all the competitors. We then evaluated the algorithms through 150 test iterations on data with the same characteritics. We report in Table 3.3 the performance of the algorithms in terms of relative ℓ 1 error A -A 0 1 / A 0 1 (mean absolute deviation), relative ℓ 2 error A -A 0 F / A 0 F and relative ℓ 0 error A -A 0 0 / A 0 0 .

We present results with Convex Graph Denoising (Algorithm 1), Factorization-based Graph Denoising (Algorithm 4) with and without regularization, Robust Principal Component Analysis (RPCA, [START_REF] Candès | Robust principal component analysis[END_REF]) implemented using Douglas-Rachford splitting, a dual implementation of RPCA ( [LGW + 09]), sparse matrix factorization corresponding to Equation (3.17), and shrinkage.

Convex Graph Denoising achieves the best results both in ℓ 1 -and ℓ 2 -norm. Factorizationbased Graph Denoising achieves slightly inferior results but still performs second to best at a smaller computation cost (see Figure 3.5). In addition to this, Factorization-based Graph Denoising achieves a better recovery of the support. Due to the lack of explicit regularization of the sparsity of X, RPCA and shrinkage perform poorly. The sparse factorization method of Equation ((3.17)) also achieves poor results due to the use of the squared loss. For small values of n, the overhead due to the alternate minimization makes the factorization approach more expensive than the convex formulation. For larger values of n, however, the factorization approach gradually becomes more and more efficient.

Real data. We performed empirical evaluation of our Douglas-Rachford algorithm for graph denoising over a subset of the e-commerce network data for comparing its predictive performance with baseline link prediction approaches. The data set contains the purchase history of the 1000 most popular music items by the 1000 most active users. Their purchases are tracked over a period of 10 months that constitute the training set and the 11th month purchases where subject to prediction. We compared the output of the algorithm with the link prediction scores computed by nearest neighbors, Katz algorithm and a low-rank approximation of the training adjacency matrix. Figure 3.6 contains the precision and recall scores of the different algorithms and shows empirical superiority of our algorithm on this recommender system task.

Discussion

• Optimization. Other optimization techniques can be considered for future work. A trace norm constraint alone can be taken into account without projection or relaxation into a penalized form by casting the problem as a SDP as proposed in [START_REF] Jaggi | Convex optimization without projection steps[END_REF]. The special form of this SDP can be leveraged to use the efficient resolution technique from [START_REF] Hazan | Sparse approximate solutions to semidefinite programs[END_REF]. This method applies to a differentiable objective whose curvature determines the performances. Extending these methods with projection onto the ℓ 1 ball or a sparsity-inducing penalty could lead to interesting developments.

• Generalization of the approach -The methods presented in this chapter can be seamlessly extended to nonsquare matrices, which can arise, for instance, from adjacency matrices of bipartite graphs. Our work also applies to a wide range of other losses. A useful example that links our work to the matrix completion framework is when linear measurements of the target matrix or graph are available, or can be predicted as in [START_REF] Richard | Link discovery using graph feature tracking[END_REF][START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF]. In this case, the loss can be defined in the feature space, and we treat this in Chapter 4 of this manuscript. Due to the low-rank assumption, our method does not directly apply to the estimation of precision matrices often used for gaussian graphical model structure learning [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], and the applications of conditional independence structures generated by low-rank and possibly sparse models is to be discussed. The suggested approach can be seen as a method to estimate sparse matrices whose nonzero coordinates are structured in specific low-rank patters. When the groups of simultaneously active or when hierarchy of active variables are known in advance, a family of approaches called structured sparsity [START_REF] Bach | Structured sparsity through convex optimization[END_REF][START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] generalize the Lasso by penalizing the sum of ℓ 2 norm of groups. The first approach dedicated to estimate sparse vectors under strong variable correlation assumption was the Elastic Net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] that uses the hybrid ℓ 1 + ǫℓ 2 penalty. The trace-Lasso [START_REF] Grave | Trace lasso: a trace norm regularization for correlated designs[END_REF] builds a more powerful regularizer for sparse regression under highly correlated designs by inserting the covariate matrix X into the regularizer, and penalizes (X ⊤ X) 1/2 diag(w) * for the variable w. In some applications we have a general prior knowledge of how the variables will be organized, but do not know in advance any specific labeling / ordering / grouping / hierarchy of the variables. In the applications discussed in this chapter, as community detection or clustering, we saw how the block-dioagonal in some permuted basis assumption was translated onto a low-rank assumption. We enforced this structural assumption on the estimation procedure by overlapping a trace-norm over the ℓ 1 of the same matrix variable. It would be very interesting to explore methods for estimating vectors / matrices under other types of prior knowledge. For instance estimating adjacency matrices of graphs containing tree-like hierarchical structures has potential interesting applications.

Chapter 4

Graph prediction in a temporal setting

Context

Forecasting systems behavior with multiple responses has been a challenging issue in many contexts of applications such as collaborative filtering, financial markets, or bioinformatics, where responses can be, respectively, movie ratings, stock prices, or activity of genes within a cell. Statistical modeling techniques have been widely investigated in the context of multivariate time series either in the multiple linear regression setup [START_REF] Breiman | Predicting multivariate responses in multiple linear regression[END_REF] or with autoregressive models [START_REF] Tsay | Analysis of Financial Time Series[END_REF]. More recently, kernel-based regularized methods have been developed for multitask learning [START_REF] Th | Learning multiple tasks with kernel methods[END_REF][START_REF] Argyriou | A spectral regularization framework for multi-task structure learning[END_REF]. These approaches share the use of the correlation structure among input variables to enrich the prediction on every single output. Often, the correlation structure is assumed to be given or it is estimated separately. A discrete encoding of correlations between variables can be modeled as a graph so that learning the dependence structure amounts to performing graph inference through the discovery of uncovered edges on the graph. The latter problem is interesting per se and it is known as the problem of link prediction where it is assumed that only a part of the graph is actually observed [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF][START_REF] Kolar | On time varying undirected graphs[END_REF]. This situation occurs in various applications such as recommender systems, social networks, or proteomics, and the appropriate tools can be found among matrix completion techniques [SRJ05, CT09, ABEV09]. In the realistic setup of a time-evolving graph, matrix completion was also used and adapted to take into account the dynamics of the features of the graph [START_REF] Richard | Link discovery using graph feature tracking[END_REF].

The estimation of a VAR model for node degrees (that are linear graph features) has been considered by [START_REF] Zhang | Content contributor management and network effects in a ugc environment[END_REF], and successfully applied to customer valuation, and to measure network effect in user generated content market places. In this work, we study the prediction problem where the observation is a sequence of graphs adjacency matrices (A t ) 0≤t≤T and the goal is to predict A T +1 . This type of problem arises in applications such as recommender systems where, given information on purchases made by some users, one would like to predict future purchases. In this context, users and products can be modeled as the nodes of a bipartite graph, while purchases or clicks are modeled as edges. In functional genomics and systems biology, estimating regulatory networks in gene expression can be performed by modeling the data as graphs and fitting predictive models is a natural way for estimating evolving networks in these contexts. A large variety of methods for link prediction only consider predicting from a single static snapshot of the graph -this includes heuristics [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF][START_REF] Sarkar | Theoretical justification of popular link prediction heuristics[END_REF], matrix factorization [START_REF] Koren | Factorization meets the neighborhood: a multifaceted collaborative filtering model[END_REF], diffusion [START_REF] Myers | On the convexity of latent social network inference[END_REF], or probabilistic methods [START_REF] Taskar | Link prediction in relational data[END_REF]. More recently, some works have investigated using sequences of observations of the graph to improve the prediction, such as using regression on features extracted from the graphs [RBEV10], using matrix factorization [START_REF] Koren | Collaborative filtering with temporal dynamics[END_REF], continuous-time regression [START_REF] Vu | Continuous-time regression models for longitudinal networks[END_REF] or non-parametric models [START_REF] Sarkar | Nonparametric link prediction in dynamic networks[END_REF]. In fact link prediction can also be seen as a special case of matrix completion where the goal is to estimate the missing entries of the adjacency matrix of the graph where the entries can be only "0s" and "1s". The ratings of matrix completion provide a richer information than the binary case (existence or absence of a link). Consider for instance the issue of link prediction in recommender systems. In that case, we consider a bipartite graph for which the vertices represent products and users, and the edges connect users with the products they have purchased in the past. The setup we consider in the present work corresponds to the binary case where we only observe purchase data, say the presence of a link in the graph, without any score or feedback on the product for a given user. The fundamental different in this case is that no negative feedback is observed. By definition, the non-observation of a link can either mean negative label or unknown label. Hence, we will deal here with the situation where the components of snapshots of the adjacency matrix only consist of "1s" and missing values.

As opposite to matrix completion where the sparsity pattern of observation is supposed to by chosen at random and contain no information, all the information in the case of link prediction remains in the sparsity pattern and the value of the entries (all equal to 1) are not relevant. Our main assumption is that the network effect is a cause and a symptom at the same time, and therefore, the edges and the graph features should be estimated simultaneously. We propose a regularized approach to predict the uncovered links and the evolution of the graph features simultaneously. We provide oracle bounds under the assumption that the noise sequence has subgaussian tails and we prove that our procedure achieves a trade-off in the calibration of smoothing parameters which adjust with the sparsity and the rank of the unknown adjacency matrix. The rest of this chapter is organized as follows. In Section 2, we describe the general setup of our work with the main assumptions and we formulate a regularized optimization problem which aims at jointly estimating the autoregression parameters and predicting the graph. In Section 3, we provide technical results with oracle inequalities and other theoretical guarantees on the joint estimation-prediction. Section 4 is devoted to the description of the numerical simulations which illustrate our approach. We also provide an efficient algorithm for solving the optimization problem and show empirical results. The proof of the theoretical results are provided as supplementary material in a separate document.

Estimation of low-rank graphs with autoregressive features

Our approach is based on the asumption that features can explain most of the information contained in the graph, and that these features are evolving with time. We make the following assumptions about the sequence (A t ) t≥0 of adjacency matrices of the graphs sequence.

Low-Rank. We assume that the matrices A t have low-rank. This reflects the presence of highly connected groups of nodes such as communities in social networks, or product categories and groups of loyal/fanatic users in a market place data, and is sometimes motivated by the small number of factors that explain nodes interactions.

Autoregressive linear features. We assume to be given a linear map ω : R n×n → R d defined by

ω(A) = Ω 1 , A , • • • , Ω d , A , (4.1)
where (Ω i ) 1≤i≤d is a set of n×n matrices. These matrices can be either deterministic or random in our theoretical analysis, but we take them deterministic for the sake of simplicity. The vector time series (ω(A t )) t≥0 has autoregressive dynamics, given by a VAR (Vector Auto-Regressive) model:

ω(A t+1 ) = W ⊤ 0 ω(A t ) + N t+1 ,
where W 0 ∈ R d×d is a unknown sparse matrix and (N t ) t≥0 is a sequence of noise vectors in R d . An example of linear features is the degree (i.e. number of edges connected to each node, or the sum of their weights if the edges are weighted), which is a measure of popularity in social and commerce networks. Introducing X T -1 = (ω(A 0 ), . . . , ω(A T -1 )) ⊤ and X T = (ω(A 1 ), . . . , ω(A T )) ⊤ , which are both T × d matrices, we can write this model in a matrix form:

X T = X T -1 W 0 + N T , (4.2)
where N T = (N 1 , . . . , N T ) ⊤ . This assumes that the noise is driven by time-series dynamics (a martingale increment), where each coordinates are independent (meaning that features are independently corrupted by noise), with a sub-gaussian tail and variance uniformly bounded by a constant σ 2 . In particular, no independence assumption between the N t is required here.

Notations. The notations • F , • p , • ∞ , • * and • op stand, respectively, for the Frobenius norm, entry-wise ℓ p norm, entry-wise ℓ ∞ norm, trace-norm (or nuclear norm, given by the sum of the singular values) and operator norm (the largest singular value). We denote by A, B = tr(A ⊤ B) the Euclidean matrix product. A vector in R d is always understood as a d × 1 matrix. We denote by A 0 the number of non-zero elements of A. The product A • B between two matrices with matching dimensions stands for the Hadamard or entrywise product between A and B. The matrix |A| contains the absolute values of entries of A. The matrix (M ) + is the componentwise positive part of the matrix M, and sign(M ) is the sign matrix associated to M with the convention sign(0) = 0

If A is a n × n matrix with rank r, we write its SVD as A = U ΣV ⊤ = r j=1 σ j u j v ⊤ j where Σ = diag(σ 1 , . . . , σ r ) is a r × r diagonal matrix containing the non-zero singular values of A in decreasing order, and U = [u 1 , . . . , u r ], V = [v 1 , . . . , v r ] are n × r matrices with columns given by the left and right singular vectors of A. The projection matrix onto the space spanned by the columns (resp. rows) of A is given by P U = U U ⊤ (resp. P V = V V ⊤ ). The operator P A : R n×n → R n×n given by P A (B) = P U B + BP V -P U BP V is the projector onto the linear space spanned by the matrices u k x ⊤ and yv ⊤ k for 1 ≤ j, k ≤ r and x, y ∈ R n . The projector onto the orthogonal space is given by P ⊥ A (B) = (I -P U )B(I -P V ). We also use the notation a ∨ b = max(a, b).

Joint prediction-estimation through penalized optimization

In order to reflect the autoregressive dynamics of the features, we use a least-squares goodnessof-fit criterion that encourages the similarity between two feature vectors at successive time steps. In order to induce sparsity in the estimator of W 0 , we penalize this criterion using the ℓ 1 norm. This leads to the following penalized objective function:

J 1 (W ) = 1 dT X T -X T -1 W 2 F + κ W 1 ,
where κ > 0 is a smoothing parameter. Now, for the prediction of A T +1 , we propose to minimize a least-squares criterion penalized by the combination of an ℓ 1 norm and a trace-norm. This mixture of norms induces sparsity and a low-rank of the adjacency matrix. Such a combination of ℓ 1 and trace-norm was already studied in [START_REF] Gaiffas | Sharp oracle inequalities for high-dimensional matrix prediction[END_REF] for the matrix regression model, and in [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF] for the prediction of an adjacency matrix.

The objective function defined below exploits the fact that if W is close to W 0 , then the features of the next graph ω(A T +1 ) should be close to W ⊤ ω(A T ). Therefore, we consider

J 2 (A, W ) = 1 d ω(A) -W ⊤ ω(A T ) 2 F + τ A * + γ A 1 ,
where τ, γ > 0 are smoothing parameters. The overall objective function is the sum of the two partial objectives J 1 and J 2 , which is jointly convex with respect to A and W :

L(A, W ) . = 1 dT X T -X T -1 W 2 F + κ W 1 + 1 d ω(A) -W ⊤ ω(A T ) 2 2 + τ A * + γ A 1 , (4.3)
If we choose convex cones A ⊂ R n×n and W ⊂ R d×d , our joint estimation-prediction procedure is defined by ( Â, Ŵ ) ∈ arg min

(A,W )∈A×W L(A, W ). (4.4)
It is natural to take W = R d×d and A = (R + ) n×n since there is no a priori on the values of the feature matrix W 0 , while the entries of the matrix A T +1 must be positive.

In the next section we propose oracle inequalities which prove that this procedure can estimate W 0 and predict A T +1 at the same time.

Main result

The central contribution of our work is to bound the prediction error with high probability under the following natural hypothesis on the noise process.

Assumption 3 We assume that (N t ) t≥0 satisfies E[N t |F t-1 ] = 0 for any t ≥ 1 and that there is σ > 0 such that for any λ ∈ R and j = 1, . . . , d and t ≥ 0:

E[e λ(Nt) j |F t-1 ] ≤ e σ 2 λ 2 /2 .
Moreover, we assume that for each t ≥ 0, the coordinates (N t ) 1 , . . . , (N t ) d are independent.

The main result (Theorem 12) can be summarized as follows. The prediction error and the estimation error can be simultaneously bounded by the sum of three terms that involve homogeneously (a) the sparsity, (b) the rank of the adjacency matrix A T +1 , and (c) the sparsity of the VAR model matrix W 0 . The tight bounds we obtain are similar to the bounds of the Lasso and are upper bounded by:

C 1 log d T d 2 W 0 0 + C 2 log n d A T +1 0 + C 3 log n d rank A T +1 .
The positive constants C 1 , C 2 , C 3 are proportional to the noise level σ. The interplay between the rank and sparsity constraints on A T +1 are reflected in the observation that the values of C 2 and C 3 can be changed as long as their sum remains constant.

Oracle inequalities

In this section we give oracle inequalities for the mixed prediction-estimation error which is given, for any A ∈ R n×n and W ∈ R d×d , by

E(A, W ) 2 . = 1 d (W -W 0 ) ⊤ ω(A T ) -ω(A -A T +1 ) 2 2 + 1 dT X T -1 (W -W 0 ) 2 F . (4.5)
It is important to have in mind that an upper-bound on E implies upper-bounds on each of its two components. It entails in particular an upper-bound on the feature estimation error

X T -1 ( W -W 0 ) F that makes ( W -W 0 ) ⊤ ω(A T
) 2 smaller and consequently controls the prediction error over the graph edges through ω( A -A T +1 ) 2 .

The upper bounds on E given below exhibit the dependence of the accuracy of estimation and prediction on the number of features d, the number of edges n and the number T of observed graphs in the sequence.

Let us recall N T = (N 1 , . . . , N T ) ⊤ and introduce the noise processes

M = - d j=1 (N T +1 ) j Ω j and Ξ = T +1 t=1 ω(A t-1 )N ⊤ t ,
which are, respectively, n × n and d × d random matrices. The source of randomness comes from the noise sequence (N t ) t≥0 , see Assumption 3. If these noise processes are controlled correctly, we can prove the following oracle inequalities for procedure (4.4). The next result is an oracle inequality of slow type (see for instance [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF]), that holds in full generality.

Theorem 9 Let ( Â, Ŵ ) be given by (4.4) and suppose that

τ ≥ 2α d M op , γ ≥ 2(1 -α) d M ∞ and κ ≥ 2 dT Ξ ∞ (4.6)
for some α ∈ (0, 1). Then, we have

E( A, W ) 2 ≤ inf (A,W )∈A×W E(A, W ) 2 + 2τ A * + 2γ A 1 + 2κ W 1 .
For the proof of oracle inequalities of fast type, the restricted eigenvalue (RE) condition introduced in [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF] and [START_REF] Koltchinskii | The Dantzig selector and sparsity oracle inequalities[END_REF][START_REF] Koltchinskii | Sparsity in penalized empirical risk minimization[END_REF] is of importance. Restricted eigenvalue conditions are implied by, and in general weaker than, the so-called incoherence or RIP (Restricted isometry property, [START_REF] Candès | Decoding by linear programming[END_REF]) assumptions, which excludes, for instance, strong correlations between covariates in a linear regression model. This condition is acknowledged to be one of the weakest to derive fast rates for the Lasso (see [START_REF] Van De Geer | On the conditions used to prove oracle results for the Lasso[END_REF] for a comparison of conditions).

Matrix version of these assumptions are introduced in [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF]. Below is a version of the RE assumption that fits in our context. First, we need to introduce the two restriction cones.

The first cone is related to the W 1 term used in procedure (4.4). If W ∈ R d×d , we denote by Θ W = sign(W ) ∈ {0, ±1} d×d the signed sparsity pattern of W and by Θ ⊥ W ∈ {0, 1} d×d the orthogonal sparsity pattern. For a fixed matrix W ∈ R d×d and c > 0, we introduce the cone

C 1 (W, c) . = W ′ ∈ W : Θ ⊥ W • W ′ 1 ≤ c Θ W • W ′ 1 .
This cone contains the matrices W ′ that have their largest entries in the sparsity pattern of W .

The second cone is related to mixture of the terms A * and A 1 in procedure (4.4). Before defining it, we need further notations and definitions.

For a fixed A ∈ R n×n and c, β > 0, we introduce the cone

C 2 (A, c, β) . = A ′ ∈ A : P ⊥ A (A ′ ) * + β Θ ⊥ A • A ′ 1 ≤ c P A (A ′ ) * + β Θ A • A ′ 1 .
This cone consist of the matrices A ′ with large entries close to that of A and that are "almost aligned" with the row and column spaces of A. The parameter β quantifies the interplay between these too notions.

Definition 20 (Restricted Eigenvalue (RE))

For W ∈ W and c > 0, we introduce

µ 1 (W, c) = inf µ > 0 : Θ W • W ′ F ≤ µ √ dT X T +1 W ′ F , ∀W ′ ∈ C 1 (W, c) .
For A ∈ A and c, β > 0, we introduce

µ 2 (A, W, c, β) = inf µ > 0 : P A (A ′ ) F ∨ Θ A • A ′ F ≤ µ √ d W ′⊤ ω(A T ) -ω(A ′ ) 2 , ∀W ′ ∈ C 1 (W, c), ∀A ′ ∈ C 2 (A, c, β) .
The RE assumption consists of assuming that the constants µ 1 and µ 2 are non-zero. Now we can state the following Theorem that gives a fast oracle inequality for our procedure using RE.

Theorem 10 Let ( Â, Ŵ ) be given by (4.4) and suppose that

τ ≥ 3α d M op , γ ≥ 3(1 -α) d M ∞ and κ ≥ 3 dT Ξ ∞ (4.7)
for some α ∈ (0, 1). Then, we have

E( A, W ) 2 ≤ inf (A,W )∈A×W E(A, W ) 2 + 25 18 µ 2 (A, W ) 2 rank(A)τ 2 + A 0 γ 2 ) + 25 36 µ 1 (W ) 2 W 0 κ 2 , where µ 1 (W ) = µ 1 (W, 5) and µ 2 (A, W ) = µ 2 (A, W, 5, γ/τ ) (see Definition 20).
The proofs of Theorems 9 and 10 use tools introduced in [KLT11] and [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF]. Note that the residual term from this oracle inequality mixes the notions of sparsity of A and W via the terms rank(A), A 0 and W 0 . It says that our mixed penalization procedure provides an optimal trade-off between fitting the data and complexity, measured by both sparsity and low-rank. This is the first result of this nature to be found in literature.

In the next Theorem 11, we obtain convergence rates for the procedure (4.4) by combining Theorem 10 with controls on the noise processes. We introduce

v 2 Ω,op = 1 d d j=1 Ω ⊤ j Ω j op ∨ 1 d d j=1 Ω j Ω ⊤ j op , v 2 Ω,∞ = 1 d d j=1 Ω j • Ω j ∞ , σ 2 ω = max j=1,...,d 1 T + 1 T +1 t=1 ω j (A t-1 ) 2 ,
which are the (observable) variance terms that naturally appear in the controls of the noise processes. We introduce also

ℓ T = 2 max j=1,...,d log log T +1 t=1 ω j (A t-1 ) 2 T + 1 ∨ T + 1 T +1 t=1 ω j (A t-1 ) 2 ∨ e ,
which is a small (observable) technical term that comes out of our analysis of the noise process Ξ. This term is a small price to pay for the fact that no independence assumption is required on the noise sequence (N t ) t≥0 , but only a martingale increment structure with sub-gaussian tails.

Theorem 11 Consider the procedure ( Â, Ŵ ) given by (4.4) with smoothing parameters given by

τ = 3ασv Ω,op 2(x + log(2n)) d , γ = 3(1 -α)σv Ω,∞ 2(x + 2 log n) d , κ = 6σσ ω 1 d 2e(x + 2 log d + ℓ T ) T + 1
for some α ∈ (0, 1) and fix a confidence level x > 0. Then, we have

E( A, W ) 2 ≤ inf (A,W )∈A×W E(A, W ) 2 + 25µ 2 (A) 2 rank(A)α 2 σ 2 v 2 Ω,op 2(x + log(2n)) d + 25µ 2 (A) 2 A 0 (1 -α) 2 σ 2 v 2 Ω,∞ 2(x + 2 log n) d + 25µ 1 (W ) 2 W 0 σ 2 σ 2 ω 2e(x + 2 log d + ℓ T ) d 2 (T + 1)
with a probability larger than 1 -17e -x , where µ 1 and µ 2 are the same as in Theorem 10.

The proof of Theorem 11 follows directly from Theorem 10 and basic noise control results, see Proposition 8 in appendix. In the next Theorem, we propose more explicit upper bounds for both the indivivual estimation of W 0 and the prediction of A T +1 .

Theorem 12 Under the same assumptions as in Theorem 11, for any x > 0 the following inequalities hold with a probability larger than 1 -17e -x :

1 dT X T ( Ŵ -W 0 ) 2 F ≤ inf A∈A 1 d ω(A) -ω(A T +1 ) 2 F + 25 18 µ 2 (A, W ) 2 rank(A)τ 2 + A 0 γ 2 ) + 25 36 µ 1 (W 0 ) 2 W 0 0 κ 2 (4.8) Ŵ -W 0 1 ≤ 5µ 1 (W 0 ) 2 W 0 0 κ + 6 W 0 0 µ 1 (W 0 ) inf A∈A 1 d ω(A) -ω(A T +1 ) 2 F + 25 18 µ 2 (A, W ) 2 rank(A)τ 2 + A 0 γ 2 ) (4.9) Â -A T +1 * ≤ 5µ 1 (W 0 ) 2 W 0 0 κ + (6 rank A T +1 + 5β A T +1 0 )µ 2 (A T +1 ) × inf A∈A 1 d ω(A) -ω(A T +1 ) 2 F + 25 18 µ 2 (A, W ) 2 rank(A)τ 2 + A 0 γ 2 ) . (4.10)

Algorithms and numerical experiments 4.4.1 Generalized forward-backward algorithm for minimizing L

We use the algorithm designed in [START_REF] Raguet | Generalized forward-backward splitting[END_REF] for minimizing our objective function. Note that this algorithm is inherently superior to the method introduced in [RBEV10] as it directly minimizes L jointly in (S, W ) rather than alternately minimize in W and S. Moreover we use the novel joint penalty from [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF] that is more suited for estimating graphs. The proximal operator for the trace norm is given by the shrinkage operation, if

Z = U diag(σ 1 , • • • , σ n )V T is the singular value decomposition of Z, prox τ ||.|| * (Z) = U diag((σ i -τ ) + ) i V T . Algorithm 6 Minimize J Initialize U, V, W repeat Compute (G U , G V , G W ) = ∇ U,V,W Φ(U V ⊤ , W ). Set U = prox θγ||.|| 1 (U -θG U ) Set V = prox θγ||.|| 1 (V -θG V ) Set W = prox θκ||.|| 1 (W -θG W ) until convergence return (U, V, W ) minimizing J

Laplacian regularizer for node features attracted following graph edges

Consider the case where the extracted features and descriptors are node features and node descriptors. In this case ω is matrix valued, ω(A) ∈ R n×q , each row of it coding the features of the corresponding node. We denote by f (i) the i-th component of f . The n prediction functions f (i) are assumed to belong to the same Hilbert space H, equipped with the norm

• H . The norm of f ∈ H n is defined by f H n := n i=1 f (i) 2
H . We develop here a methodology to relate the graph and the predictors through an additional term that acts as a common regularizer. Before introducing this term we need to make the hypotheses

• A t ∈ S n ≥0 , t ∈ {0, 1, 2, .
.., T }, where S n ≥0 denotes the set of n × n symmetric matrices with nonnegative entries.

• We assume also that, the weights of the edges are nondecreasing functions of time, which in the unweighted case means edges do not disappear but new ones may appear over time.

For simplicity of presentation, we consider the case of H being an RKHS corresponding to a linear kernel. In addition, each f (i) is a linear function represented by a q × q matrix W (i) , so that f (i) (ω(A t ) (i) ) = ω(A t ) (i) W (i) . Then we may assume f (i) H = W (i) F and the graph regularity term can be expressed in terms of the graph Laplacian of A. Recall that the Laplacian is defined as the operator Λ such that Λ(A) = D -S where D = diag(d 1 , . . . , d n ) and d i = n j=1 A i,j are the degrees of the graph. A standard computation gives the following expression for J 4 :

Tr A ⊤ ∆(f ) = 1≤i,j≤n A i,j W (i) -W (j) 2 F where W = (W (1) , . . . , W (n) ) ∈ R n×q×q . We use the standard extension of the Frobenius norm to 3-tensors, W 2 F = n i=1 W (i) 2
F and the notations:

Q(W, Λ, V ) := 1≤i,j≤n Λ ij Tr(W (i) ⊤ V (j) ) , ∆(W ) := W (i) -W (j) 2 F n i,j=1
.

We also assume in this work that the node features are linear functions of adjacency matrices : ω(A t ) = A t Ω for some Ω ∈ R n×q . Note that degree, inter/intra cluster degrees and projection onto specific linear subspaces are such features, but not clustering coefficients or statistics on path lengths. We explain later how we chose Ω in our experiments.

The Laplacian-like term that ensures joint regularity of graph and predictors can be written as follows

Tr A ⊤ ∆(f ) = 1 2 i,j A i,j f (i) -f (j) 2

H

where we have used

∆(f ) = f (i) -f (j) 2 2 i,j ∈ R n×n
This term relates the contributions of f and A enforcing regularity of the predictors over the estimated graph through A. The objective to minimize in this case contains the additional Laplacian-like term enforcing predictors regularity over A:

L Lap (A, f ) . = 1 d y -f (X) 2 F + κ pen(f ) + 1 2 ω(A) -f (ω(A T )) 2 2 + τ A * + γ A 1 + λ 2 i,j A i,j f (i) -f (j) 2

H

We only point out that the objective function L Lap (A, f ) can be optimized inside the region it is convex in using projected gradient algorithms, future investigation is required on both empirical and theoretical properties of the resulting estimator.

In order to determine the convexity region of L Lap we use intuition from what happens with a similar optimization problem in dimension 2. Indeed, consider the minimization of (w, a) → aw 2 + αa 2 + βw 2 which is the simplified functional in the degenerate case where n = 1. The eigenvectors of the Hessian are positive iff w 2 ≤ αβ, so this condition defines the convexity region. By adding the quadratic terms in A and W to Q(W, Λ(A), W ), let us define

Q(W, Λ, V ) := 1≤i,j≤n Λ ij Tr(W (i) ⊤ V (j) ) and Ψ(A, W ) := κ 2 W 2 F + ν 2 A -A T 2 F + λQ(W, Λ(A), W
). If we suppose that the entries of Z = A -A T are nonnegative, then the minimizer of Ψ is the trivial solution A = A T , W = 0. In fact,

Ψ(A, W ) = κ 2 W 2 F + ν 2 Z 2 F + λ Q(W, Λ(A T ), W ) + Q(W, Λ(Z), W ) ≥ 0 (4.11)
and for Z = 0, W = 0, Ψ(A, W ) = 0. We prove that L Lap is convex over a set E around the minimizer of Ψ and we will ensure henceforth that the descent algorithm takes place inside this convex domain.

Proposition 2

The function Ψ is convex in the interior of the set:

E = A ∈ S n ≥0 , W ∈ R n×d×q W F ≤ √ νκ 2λ √ n + 1 .

Bias of the measurement methodology

The objective of this thesis was to develop estimation methods for prediction purposes in temporal sequences of graphs. A natural way of testing the validity of the working hypotheses and the suggested algorithms is backtesting. It consists of using past data for both training and testing. The stability of the estimators can be measured also by taking several training and test periods in the past. One of the major criticism to such a methodology is that given the unbalanced degree distributions in real-world data, the outputs of the algorithms fitting to real world data are often over-represent the popular items. Figure 4.1 shows for illustrating this statement, the linear dependency of magnitude of standard scores on the sales volumes illustrates the bias due to the heavy-tailed distributions that hurts these scoring rules. Therefore they can suffer from poor recommendation diversity, by either recommending the same best-selling items to all the users or by fitting to the users tastes and not letting them discover new items.

Empirical evaluation of heuristics on real data sets

It is a clear industrial requirement that the algorithms developed for recommender systems should scale to millions of users and millions of products database size. However, the storage and manipulation of such data would require tera-bytes of cash and RAM memory. Given the fast growing volumes of data that need to be processed, one could clearly not handle quadratic or cubic algorithms as those studied in the current thesis and wait for computers to have sufficient power to deal with these data. It turned obvious after a short while of studying recommender systems that local algorithms should be preferred to global algorithms. In fact in all the studied approaches, the trace-norm term, or the low-rank criterion was the key to handle the collaborative aspect of the learning process by generalizing the characteristics / tastes of a node / user to the others. Somehow it is a fact that a minor change in one entry of the input matrix will affect the spectrum of it and consequently affect the trace-norm. This is why a first conclusive remark would be to point out the limitation of the approach due to the presence of this term as the only generalizer. As a case study for our industrial partner, and also by the desire of better understanding the nature of the data-mining problem, whatever statistical relevance this had, we also developed some very basic algorithms using an item-to-item approach. The idea is to compute the similarity among pairs of products and store only the top-K most similar corresponding to each. The Figure 4.2 shows how to use a cross-validation procedure for tuning the length of lists to merge for books market: K = 7 maximizes the average Hit Rate in this case. The suggested list of items to a user would be a list obtained by mixing the corresponding lists of his last purchases. Our numerical studies on several databases have proven empirically that such an approach is outperforming very easily its full-matrix scoring rivals by de-biasing the scores to the popularity of the recommended items. The second asset of these approaches was that they were able to provide recommendations lists in sub-linear time, as the heavy computation is done in advance and the products are stored in tables. The Figure 4.3 shows that except for Electronic products market where the notion of similarity/vicinity of products seems meaningless we see the benefits of both steps of our approach on the predictions. The reason why we did not push the study of these algorithms further is that we wanted to build statistical models for better understanding and capturing evolution patterns of the data. The immediate applicability and use of the methods was not the major concern. Given the directions shown by the mixed regularizers design we think there is clearly a room for improvement by designing and implementing local versions of these regularized algorithms. The deep question remains in the interpretation one makes from the term local. In fact the metric of the spaces in which our work was developed is a euclidean metric that can be easily generalized to Hilbert spaces as Reproducing Kernel Hilbert Spaces by using a PSD kernel. Somehow recent research has suggested [START_REF] Chen | On the hyperbolicity of smallworld networks and tree-like graphs[END_REF] that the hyperbolic geometry may be more adapted to the study of graphs of human activity. Even though at the first glance such an assertion may sound more as a non-pragmatic intellectual-game idea, a closer look at the question reveals deep insights which may give a clue to demystifying the nature of these random graphs.

Empirical evaluation of the regularization method on synthetic data sets

We tested the presented methods on synthetic data generated as in section (4.4.4). In our experiments the noise matrices M t and N t where built by soft-thresholding i.i.d. noise N (0, σ 2 ). We took as input T = 10 successive graph snapshots on n = 50 nodes graphs of rank r = 5.

We used d = 10 linear features, and finally the noise level was set to σ = .5. We compare our methods to standard baselines in link prediction. We use the area under the ROC curve as the measure of performance and report empirical results averaged over 50 runs with the corresponding confidence intervals in Figure 4.5. The competitor methods are the nearest neighbors (NN) and static sparse and low-rank estimation, that is the link prediction algorithm suggested in [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF]. The algorithm NN scores pairs of nodes with the number of common friends between them, which is given by A 2 when A is the cumulative graph adjacency matrix A T = T t=0 A t and the static sparse and low-rank estimation is obtained by minimizing the objective X -A T 2 F + τ X * + γ X 1 , and can be seen as the closest static version of our method. The two methods autoregressive low-rank and static low-rank are regularized using only the trace-norm, (corresponding to forcing γ = 0) and are slightly inferior to their sparse and low-rank rivals. Since the matrix V 0 defining the linear map ω is unknown we consider the feature map ω(A) = SV where A T = U ΣV ⊤ is the SVD of A T . The parameters τ and γ are chosen by 10-fold cross validation for each of the methods separately.

Discussion

1. Comparison with the baselines. This experiment sharply shows the benefit of using a temporal approach when one can handle the feature extraction task. The left-hand plot shows that if few snapshots are available (T ≤ 4 in these experiments), then static approaches are to be preferred, whereas feature autoregressive approaches outperform as soon as sufficient number T graph snapshots are available (see phase transition). The decreasing performance of static algorithms can be explained by the fact that they use as input a mixture of graphs observed at different time steps. Knowing that at each time step the nodes have specific latent factors, despite the slow evolution of the factors, adding the resulting graphs leads to confuse the factors.

2. Phase transition. The right-hand figure is a phase transition diagram showing in which part of rank and time domain the estimation is accurate and illustrates the interplay between these two domain parameters.

3. Choice of the feature map ω. In the current work we used the projection onto the vector space of the top-r singular vectors of the cumulative adjacency matrix as the linear map ω, and this choice has shown empirical superiority to other choices. The question of choosing the best measurement to summarize graph information as in compress sensing seems to have both theoretical and application potential. Moreover, a deeper understanding of the connections of our problem with compressed sensing, for the construction and theoretical validation of the features mapping, is an important point that needs several developments. One possible approach is based on multi-kernel learning, that should be considered in a future work. An extension to nonlinear graph features such as the distribution of triangles or other nonlinear patterns of interest is also to be considered.

4. Generalization of the method. In this chapter we consider only an autoregressive process of order 1. For better prediction accuracy, one could consider mode general models, such as vector ARMA models, and use model-selection techniques for the choice of the orders of the model. A general modelling based on state-space model could be developed as well. We presented a procedure for predicting graphs having linear autoregressive features.

Our approach can easily be generalized to non-linear prediction through kernel-based methods.

5.

A matrix covariate multi-armed bandits problem As highlighted several times in the current thesis, seeing the recommendation problem as a prediction problem quickly reaches its limits. In fact recommending very likely items, even though practitioners argue that it avoids customers outflow to competitor retailers, is not necessarily the best business decision for a long-term contact with the customer [START_REF] Bodapati | Recommendation systems with purchase data[END_REF]. The trade-off risen by this problem is an exploitation-exploration trade-off that can easily be modeled using the multi-armed bandit formalism. We were eager to formulate this problem, but not able to solve it so far, a closely related problem has been solved meanwhile by Hazan and coauthors [START_REF] Hazan | Near-optimal algorithms for online matrix prediction[END_REF]. The formulation of the problem goes as follows.

The Bandit Completion problem is an interesting problem that arised when trying to optimally allocate advertisement ressources in typical CRM databases containing N ≃ 10 7 users and used by around K ≥ 10 3 advertisers. The concrete application requires to max-imize the click-rate (reward) for a limited number of advertisements T . In the following we refer to the matrix containing the expected click rate of a user i on an advertisment sent by advertiser j as the expected reward matrix.

Take N bandits each bandit having K arms to play. We want to study the following scenari :

(a) Clusters of Bandits : suppose the bandits can be partitionned in r clusters, r ≪ K, N . The reward laws being identical for the bandits belonging to the same cluster. We aim at performing simultaneousely the clustering of the bandit population and also finding the best arm for each cluster.

(b) Low rank rewards matrix : if A ∈ R N ×K is the rewards matrix, i.e. A i,j represents the expected reward of arm j for bandit i, we assume A has rank r ≪ K, N .

At round t the playing bandit and the arm (n t , k t ) ∈ {1, • • • , N } × {1, • • • , K} are chosen and reward R(n t , k t ) is obtained. We want to find the optimal strategy for choosing (n t , k t ) in order to maximize the cumulated reward after T rounds T t=1 R(n t , k t ), or equivalently if R * = max (n,k)∈{1,••• ,N }×{1,••• ,K} ER(n, k) is the maximal expected reward, we want to minimize the regret

R T = T R * - T t=1
R(n t , k t ) 6. Convexity. This work has been devoted to develop methods for estimating interactions among objects. The data domain having motivated this work is e-commerce data, while the methods introduced can be used for measuring interactions in other domains such as in biology and drug discovery. We met several nonconvex objective functions. Examples of nonconvex terms are

• the laplacian inner product that pushes the node features to be close when the covariate is high and vice versa:

(X, w) → i,j X i,j w iw j 2 2

• by replacing the distance w iw j 2 2 by an inner product we obtain (X, w) → i,j X i,j w i , w j

• in the context of matrix factorization, the standard loss function takes the form

(U, V ) → U V ⊤ -A
• the squared error term in matrix feature estimation may be nonconvex if the feature map is not linear X → ω(X)ω(A 0 ) where ω is not linear and we saw that very basic graph characteristics such as the number of triangles surrounding a node, or the clustering coefficient are nonlinear and nonconvex. As the purpose of using regularization methods is to come up with fast algorithms leading to accurate estimators/predictors, dealing with nonconvex objects is hopeless. Nevertheless these terms are natural candidates for describing real-world effects. The classical issue for dealing with these terms is either to ignore their nonconvexity by restricting the search space to the region of convexity (see Proposition 2 for instance) or to replace the problem with a convex one [START_REF] Argyriou | Multi-task feature learning[END_REF] which has the same solution, or to reformulate the problem with use SVD-related tools that ensure the uniqueness of the solution for the high computational cost of spectral computation. A future direction of work is to explore the information relying in such nonconvex terms, possibly for other goals than estimating accurately a vector of parameters.

7. Detecting hierarchical structures. The use of a mixed ℓ 1 -trace norm penalty has been motivated by the presence of highly connected groups of nodes in real world graphs. In fact such a penalty allows to estimate matrices containing such structures. Another very important characteristic of real-worl networks is the presence of hierarchical structures. A promising direction of work is to develop a methodology to estimate such structures. A series of previous work (see [START_REF] Chen | On the hyperbolicity of smallworld networks and tree-like graphs[END_REF] and references therein), relying on the properties of hyperbolic spaces have suggested to embed the nodes onto a hyperbolic space rather than onto euclidean space. The Figure 4.6 illustrates the Poincaré hyperbolic disk where the multiscale, hierarchical structure of the space is visible. The challenge will be to develop efficient algorithms in such a framework where the expression of the basic objects such as inner products and distances is more complicated than in euclidean case.

2. Computing the SVD of a matrix M ∈ R n×m is equivalent to finding the eigenvalue decomposition of the symmetric matrix A = 0 M M T 0 ∈ R (n+m)×(n+m) .

For obtaining the eigen-decomposition of symmetric matrices most numerical solvers first tridiagonalize the matrix using a Householder transformation for a cost of O(mn 2 ) flops (where without loss of generality m ≥ n). The next step consists of computing the eigenvectors using a power method. Each iteration of the power method costs O(n).

Remark 11 (Computational Cost of SVD on Dense Matrices) If m > n, the computational cost of SVD is O(mn 2 ). Lanczos If the matrix A is very large and sparse (as adjacency matrices of real-world graphs are often), then a Lanczos algorithm can be applied to more efficiently use the storage memory. Before further details we need the following Definition 29 (Krylov spaces) Given a vector v and a square matrix A, the Krylov subspace of order r ∈ N is the subspace

K r (A, v) = span{v, Av, A 2 v, • • • , A r-1 v}
A key idea is that singular values of the restriction of a matrix to a Kyrlov space are accurate approximations to its singular vectors. Orthogonal bases of Kyrlov space allow to restrict the heavy part of the computations to tridiagonal matrices, and thereby reduce the computational cost of the decomposition.

There is a family of algorithms named Arnoli algorithms used for building an orthogonal basis of the Krylov subspace K r . Among Arnoldi algorithms, Lanczos algorithms are of special interest when A is symmetric. These algorithm calculates two matrices T r and V r which are respectively symmetric-tridiagonal and orthogonal and T r = V T m AV m . The non-zero coefficients of T r (α j = t j,j and β j = t j-1,j = t j,j-1 ) are obtained during the following Algorithm 7 Lanczos v 0 ← 0 v 1 ← unit norm random vector β 0 ← 0 for j = 1, 2, . . . , r do w j ← Av jβ j v j-1 α j ← w j , v j w j ← w jα j v j β j+1 ← w j v j+1 ← w j /β j+1 end for

The families of matrices T r and V r are then used to compute an approximate singular value decomposition of A.

Randomization Recently randomized methods have been introduced for making the SVD computation possible in very large scale cases where the singular values decay is fast and only a few number k of singular elements are required. We quickly review the main sketch of these algorithms, further details can be found in [START_REF] Halko | Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions[END_REF] we have

E(A, W ) 2 = 1 d Φ(A -A T +1 , W -W 0 ) 2 F .
Let us introduce also the empirical risk

R n (A, W ) = 1 d 1 √ T X T , 0 -Φ(A, W ) 2 F .
The proofs of Theorem 9 and 10 are based on tools developped in [START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy matrix completion[END_REF] and [START_REF] Bickel | Simultaneous analysis of lasso and dantzig selector[END_REF]. However, the context considered here is very different from the setting considered in these papers, so our proofs require a different scheme.

A.3.1 Proof of Theorem 9

First, note that

R n ( Â, Ŵ ) -R n (A, W ) = 1 d Φ( Â, Ŵ ) 2 F -Φ(A, W ) 2 F -2 ( 1 √ T X T , 0), Φ( Â -A, Ŵ -W ) . Since 1 d Φ( Â, Ŵ ) 2 F -Φ(A, W ) 2 F = E( Â, Ŵ ) 2 -E(A, W ) 2 + 2 d Φ( Â -A, Ŵ -W ), Φ(A T +1 , W 0 ) , we have R n ( Â, Ŵ ) -R n (A, W ) = E( Â, Ŵ ) 2 -E(A, W ) 2 + 2 d Φ( Â -A, Ŵ -W ), Φ(A T +1 , W 0 ) -( 1 √ T X T , 0) = E( Â, Ŵ ) 2 -E(A, W ) 2 + 2 d Φ( Â -A, Ŵ -W ), (- 1 √ T N T , N T +1 ) .
The next Lemma will come in handy several times in the proofs.

Lemma 3 For any A ∈ R n×n and W ∈ R d×d we have

( 1 √ T N T , -N T +1 ), Φ(A, W ) = (M, 1 T Ξ), (A, W ) = 1 T W, Ξ + A, M .
This Lemma follows from a direct computation, and the proof is thus omitted. This Lemma entails, together with (4.4), that

E( Â, Ŵ ) 2 ≤ E(A, W ) 2 + 2 dT Ŵ -W, Ξ + 2 d  -A, M + τ ( A * -A * ) + γ( A 1 -A 1 ) + κ( W 1 -W 1 ).
Now, using Hölder's inequality and the triangle inequality, and introducing α ∈ (0, 1), we obtain

E( Â, Ŵ ) 2 ≤ E(A, W ) 2 + 2α d M op -τ Â * + 2α d M op + τ A * + 2(1 -α) d M ∞ -γ Â 1 + 2(1 -α) d M ∞ + γ A 1 + 2 dT Ξ ∞ -κ Ŵ 1 + 2 dT Ξ ∞ + κ W 1 ,
which concludes the proof of Theorem 9, using (4.6).

where stands for the semidefinite order on self-adjoint matrices. Using Corollary 3.7 from [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF], this leads for any x > 0 to , which leads easily to (A.17). Inequality (A.18) comes from the following standard bound on the sum of independent sub-gaussian random variables:

P 1 d d j=1
(N T +1 ) j (Ω j ) k,l ≥ x ≤ 2 exp -x 2 2σ 2 (Ω j ) 2 k,l together with an union bound on 1 ≤ k, l ≤ n.

Inequality (A.19) is based on a classical martingale exponential argument together with a peeling argument. We denote by ω j (A t ) the coordinates of ω(A t ) ∈ R d and by N t,k those of N t , so that

T +1 t=1 ω(A t-1 )N ⊤ t j,k = T +1 t=1 ω j (A t-1 )N t,k .
We fix j, k and denote for short ε t = N t,k and x t = ω j (A t ). Since E[exp(λε t )|F t-1 ] ≤ e σ 2 λ 2 /2 for any λ ∈ R, we obtain by a recursive conditioning with respect to F T -1 , F T -2 , . . . , F 0 , that

E exp θ T +1 t=1 ε t x t-1 - σ 2 θ 2 2 T +1 t=1 x 2 t-1 ≤ 1.
Hence, using Markov's inequality, we obtain for any v > 0:

P T +1 t=1 ε t x t-1 ≥ x, T +1 t=1 x 2 t-1 ≤ v ≤ inf θ>0 exp(-θx + σ 2 θ 2 v/2) = exp - x 2 2σ 2 v ,
that we rewrite in the following way:

P T +1 t=1 ε t x t-1 ≥ σ √ 2vx, T +1 t=1 x 2 t-1 ≤ v ≤ e -x .
Let us denote for short V T = T +1 t=1 x 2 t-1 and S T = T +1 t=1 ε t x t-1 . We want to replace v by V T from the previous deviation inequality, and to remove the event {V T ≤ v}. To do so, we use a peeling argument. We take v = T + 1 and introduce v k = ve k so that the event {V T > v} is decomposed into the union of the disjoint sets {v k < V T ≤ v k+1 }. We introduce also ℓ T = 2 log log On {V T ≤ v} the proof is the same: we decompose onto the disjoint sets {v k+1 < V T ≤ v k } where this time v k = ve -k , and we arrive at t=1 ω j (A t-1 ) 2 ∨ e .

P S T ≥ σ 2eV T (x + ℓ T ), V T ≤ v ≤ 3.47e -x .

This leads to

The conclusion follows from an union bound on 1 ≤ j, k ≤ d. This concludes the proof of Proposition 8.

Figure 1 . 1 :

 11 Figure 1.1: These two representations of the same graph represent a communication network in a viral marketing situation.

Figure 1 . 2 :

 12 Figure 1.2: Example of regulatory network.

Figure 1

 1 Figure 1.4: Cartoon of the innovation curve of a product.

∝

  Figure 1.5: Some spectral functions of the adjacency matrix used for link preiction.

  Figure 2.1: Zones of validity of the power-law fit on different markets.

-

Figure 2

 2 Figure 2.3: Music,Books,Video Games and Computers contingency heat maps. Sum of each row/column equals 10%. The top-left corner (most frequent user and rarest products) is much more dense than in the cultural market.

Figure 2 . 4 :

 24 Figure2.4: The benefits generated by different markets.

Figure 2 . 5 :Figure 2

 252 Figure 2.5: Evolution of the daily number of sales on different markets. The seasonality in Electronic device transaction seem less important, peaks in book sales are in September and before Christmas for Music and Video Games items.

Figure 2 . 7 :

 27 Figure 2.7: The evolution of products diversity over 10 years in several products category: books sold on the website seem to concentrate around popular items.

  graph having |V | = n vertices, and each pair of vertices independently connected with probability p. G (n, |E|) is a graph having |V | = n vertices, and |E| edges are chosen at random among the n 2 possibilities. For n → ∞, provided that n 2 p → ∞, G (n, p) and G (n, |E| = p n

  2 (Number of cliques) The expected number of k-vertex cliques in G (n, p) is p k ν k Property 3 (Clustering coefficient) The clustering coefficient, or the probability that 3 connected vertices of G (n, p) belong to a triangle, is p ≈ 2|E| n(n-1) .

Figure 3 . 1 :Figure 3

 313 Figure 3.1: Covariance estimation. Cross-validation: normalized RMSE scores (SPLR)

  0854 ±0.012 0.1487 ±0.02 0.1023 ±0.02 0.2 0.2073 ± 0.03 0.2673 ± 0.3 0.2484 ± 0.03 0.3 0.3105 ± 0.03 0.3728 ± 0.03 0.3104 ± 0.02Table 3.1: Prediction of interactions in Escherichia coli. Mean normalized RMSE and standard deviations.

Figure 3 . 3 :

 33 Figure 3.3: Regularization path for RPCA (left), sparsity patterns of solutions for different values of the parameters (right)

Figure 3 . 4 :

 34 Figure 3.4: Influence of A 0 and N on Factorization-based Graph Denoising

  Figure 3.5 shows the computation time for Convex Graph Denoising and Factorization-

Figure 3

 3 Figure 3.5: Computation time

Figure 3 . 6 :

 36 Figure 3.6: Graph denoising for recommender system application: precision-recall curves.

Figure 4 . 1 :

 41 Figure 4.1: Dependency of magnitude of standard scores on the sales volumes.

Figure 4

 4 Figure 4.2: Cross-validation for the parameter K, the length of intermediary lists.

  .S. std NN std LR std Katz list NN list LR list Katz list-.S. std NN std LR std Katz list NN list LR list Katz list-.S. std NN std LR std Katz list NN list LR list Katzlist-.S. std NN std LR std Katz list NN list LR list Katzlist-.S. std NN std LR std Katz list NN list LR list Katz list-.S. std NN std LR std Katz list NN list LR list Katz list-.S. std NN std LR std Katz list NN list LR list Katz list-.S. std NN std LR std Katz list NN list LR list Katz list-

Figure 4 . 10 Figure 4 . 4 :

 41044 Figure 4.3: NDCG@N in the left column and Hit Rate@N at right. Row 1 : Music; Row 2 : Books; Row 3 : Electronic Devices; Row 4 : Video Games.

Figure 4

 4 Figure 4.5: Left: performance of algorithms in terms of Area Under the ROC Curve, average and confidence intervals over 50 runs. Right: Phase transition diagram.

Figure 4 . 6 :

 46 Figure 4.6: Poincaré hyperbolic disk.

P λ max d j=1 (

 j=1 N T +1 ) j X j ≥ x ≤ n exp -x 2 2v 2 , where v 2 = σ 2 [Tro10], we consider the dilation operator L : R n×n → R 2n×2n given by +1 ) j Ω j op = λ max L d j=1 (N T +1 ) j Ω j = λ max d j=1 (N T +1 ) j L(Ω j )and an easy computation gives (A.20) with the self-adjoint X j = L(Ω j ) givesP d j=1 (N T +1 ) j Ω j op ≥ x ≤ 2n exp -x 22v 2 where v 2 = σ 2

  P S T ≥ σ 2eV T (x + ℓ T ), V T > v = k≥0 P S T ≥ σ 2eV T (x + ℓ T ), v k < V T ≤ v k+1 = k≥0 P S T ≥ σ 2v k+1 (x + 2 log log(e k ∨ e)), v k < V T ≤ v k+1 ≤ e -x (1 + k≥1 k -2 ) ≤ 3.47e -x .

  A t-1 )N t,k ≥ σ 2e T +1 t=1 ω j (A t-1 ) 2 (x + ℓ T,j ) 1/2 ≤ 7e -xfor any 1 ≤ j, k ≤ d, where we introduced ℓ T,j = 2 log logT +1 t=1 ω j (A t-1 ) 2 T + 1 ∨ T + 1 T +1

Table 1

 1

		contribution	new sellers	new buyers
	New seller	11	8	11
	New buyer	12	11	12
	Returning seller	6	6	5
	Returning buyer	1	1	1

.1: Revenue contribution and network effects measured by cumulative impulse responses at 95 % of the cumulative levels.

Formulation as a prediction problem in evolving graph data.

It is very natural to use network data for valuation purposes because the revenue generated by an e-commerce network is a simple linear function of the market graph. Therefore we argue that fitting to the graph dynamics allows to extract wealth generating patterns in the graph and hence establish the proper value of each of the market actors.

6 0.8 Cumulative Impulse Response to New Product

	4
	2
	0
	-2
	-4
	-6
	-8
	-10	0	5	10	15	20	25
	Time (week)

Cumulative Impulse Response to Returning Product

  link relates the i-th user/item of set X to the j-th user/item of set Y . Following the application we may either take binary matrices or weighted adjacency matrices where the weight of an edge equals the amount of transaction or commission generated by the transaction. Note that X ⊤ Y = Y X . The n-by-n symmetric adjacency matrix

	A t =	 	0 B P B S B ⊤ B ⊤ S S 0 P ⊤ P 0 P S	  =	 	S B S P 0 B P B S 0 P B 0 P S	 
	summarizes all the transaction informations at time t.
	Note that the following quantities are linear functions of the adjacency matrices :
	1. The overall turnover / commission generated by the transactions
	2. Number of Buyers/Products/Sellers participating to a transaction at time t
	3. Number of Buyers/Products/Sellers participating to a transaction up to time t if we either
	consider cumulated adjacencies or a function of the set of matrices A 1 , • • • , A t
	4. Individual user activity / product popularity indicators : number or amount of transac-
	tions at time t
	Sales Volume
	Commission
	New Buyer
	New Product
	New Seller
	Returning Buyer
	Returning Product
	Returning Seller
							Time (week)
	Figure 1.3: Top: time series of the number of different types of actors of a C-to-C market and
	the corresponding commission generated over weeks.. Bottom: impulse response curves of
	different times series of the marketplace.

transaction

Table 2

 2 .5 show the contingency table for each market. From these data, we can notice that in all markets very frequent users seem to buy very rare products, occasional buyers seem to be attracted by the popular products.

			n 1	n 2	∆E	αU	αP
		Music	0.4M 60 k	2k	2.8 2.5
		Books	1.2M 1.7M 18k	2.8 2.7
		Electronic	0.5M 60k	2k	2.9 2.5
		Video Games 0.9M 0.2M 9k	3.0 1.6
		Table 2.1: Basic statistics on our data sets.
	Very Popular (> 80) Popular(>25) Rare (>7) Very Rare The tables 2.2, 2.3, 2.4 and 2Buyer \ Item Very Frequent (>36) 4.77 6.30 6.87 7.05
	Frequent (>11)	6.39		6.26	6.15	6.18
	Occasional (>4)	7.14		6.29	5.87	5.69
	Very Occasional	6.69		6.13	6.09	6.06
	Very Popular (> 90) Popular(> 25) Rare (> 7) Very Rare Very Frequent (> 43) 3.72 Buyer \ Item 5.97 7.48 7.80
	Frequent(> 12)	5.40		6.42	6.55	6.61
	Occasional (> 4)	7.34		6.43	5.70	5.50
	Very Occasional	8.51		6.16	5.24	5.06

.2: Books contingency table. Sum of each row/column equals 25% of the total amount of purchases (edges of the bipartite graph).

Table 2

 2

.3: Music contingency table. Sum of each row/column equals 25%.

The dynamics of the market graph

We are interested in a sequence of graphs evolving over time. It is a choice of modelling to keep the vertex set V fixed and let the edges change over time. We formally represent a graph sequence as a series indexed by time (G t ) t = (V, E t ) t where E t ⊂ R V ×V . We choose to discretize the time domain to the set {1, • • • , T }. Temporal evolution of product sales. Users have different habits and behaviors on different market categories. The sales explode seasonally (before Christmas or in September for books for example) on some markets (see figure

Table 2

 2

										Users side
						1
					Density	0.4 0.6 0.8			Video Games Music Electronic Devices Books
						0.2
						8 0	-7	-6	-5	-4		-3	-2	-1	0
								log(Clustering Coefficient)
									Products side
						1
					Density	0.4 0.6 0.8			Video Games Music Electronic Devices Books
						0.2
						-8 0	-7	-6	-5	-4		-3	-2	-1	0
								log(Clustering Coefficient)
						user side							product side
		0.9									1
								Video Games								Video Games
		0.8						Music									Music
		0.7						Electronic Books			0.8						Electronic Books
		0.6
	Density	0.4 0.5									0.4 0.6 Density
		0.3
		0.2									0.2
		0.1
		0	8	9		10	11	12	13		0	7	8	9	10	11	12	13
						log(degree)							log(degree)
						user side							product side
		0.5									0.45
	Density	0.2 0.3 0.4						Video Games Music Books Electronic		Density	0.4 0.2 0.25 0.3 0.35						Video Games Music Electronic Books
											0.15
		0.1									0.1
											0.05
		0	7	8	9	10	11	12	13		0	7	8	9	10	11	12	13
					log(d (2) /degree)							log(d (2) /degree)
						user side							product side
		0.5									0.45
								Video Games Music		0.4						Video Games Music
		0.4						Electronic Books			0.35						Electronic Books
											0.3
	Density	0.2 0.3								Density	0.2 0.25
											0.15
		0.1									0.1
											0.05
		0	7	8	9	10	11	12	13		0	7	8	9	10	11	12	13
						log(d (3) /d (2) )							log(d (3) /d (2) )
			Figure 2.2: Distribution of topological network features.
	Very Popular(>1180) Popular(>355) Rare (>75) Very Rare Very Frequent (>32) 4.28 Buyer \ Item 5.10 6.20 9.41
	Frequent (>11)			6.19					6.44			6.45	5.914
	Occasional (>4)			6.90					6.81			6.35	4.92
	Very Occasional			7.62					6.63			6.00	4.74

.4: Video Games contingency table. Sum of each row/column equals 25%.

Table 3 .

 3 2: Facebook denoising data. Mean AUC over 10 simulation runs. All standard deviations are lower than 3 • 10 -4 .

	σ	SPLR LR	NN	Katz
	5 %	0.9293 0.9291 0.7680 0.9298
	10 % 0.9221 0.9174 0.7620 0.9189
	15 % 0.9117 0.9024 0.7555 0.9068
	20 % 0.8997 0.8853 0.7482 0.8941

http://developers.facebook.com/docs/reference/api/

Viral marketing, also called viral advertisement or buzz marketing, is a marketing technique based on wordof-mouth type of exchanges among users, leading to a spread of information about a brand that is similar to virus propagation in a social network.

3 Source : http://www.bioquicknews.com/ Vast New Regulatory Network Discovered in Mammalian Cells, 10/14/2011

(1.3) for a well-chosen matrix Γ. The particular case Γ = √ λI d is known as Ridge Regression[START_REF] Hastie | The Elements of Statistical Learning[END_REF], and has the practical virtue of reducing noise and leading to smooth solutions. The choice of the parameter λ that trades the smoothness of the solution with fit to the data is

Acknowledgments

• Fix an integer σ < nr and take for t = 1, 2, • • •

Now define the linear feature map ω(A) = AV ⊤ † 0 and note that

follows the linear autoregressive relation

• ∀t, A t has rank at most r • ∀t, A t is sparse as V ⊤ 0 and U t are both sparse

Low-rank and sparse matrix sequence having linear autoregressive features (SPLR-LinFeat).

The difference with the non-sparse low-rank case is that the noise is supposed to be sparse.

The construction goes as follows. Let V 0 ∈ R n×r be a sparse matrix V † 0 its pseudo-inverse

Fix the two sparse matrices W 0 ∈ R r×r and U 0 ∈ R n×r . Now define the sequence of matrices (A t ) t≥0 by for t = 1, 2, • • •

A t = U t V ⊤ 0 + M t for i.i.d sparse noise matrices N t and M t , which means that for any pair of indices (i, j), with high probability (N t ) i,j = 0 and (M t ) i,j = 0.

If we define the linear feature map ω(A) = AV ⊤ † 0 , note that 1. The sequence ω(A t )

follows the linear autoregressive relation

2. For any time index t, the matrix A t is close to U t V 0 that has rank at most r 3. A t is sparse, and furthermore U t is sparse

The two following models are designed for generating sequences of SPD matrices, used frequently in the study of gaussian graphical models. Similarly, the proximal operator for the ℓ 1 -norm is the soft thresholding operator defined by using the entry-wise product of matrices denoted by •:

Method

The algorithm converges under very mild conditions when the step size θ is smaller than 2 L , where L is the operator norm of the joint quadratic loss:

Nonconvex factorization method

An alternative method to the estimation of low-rank and sparse matrices by penalizing a mixed penalty of the form τ A * + γ A 1 as in [START_REF] Richard | Estimation of simultaneously sparse and low-rank matrices[END_REF] is to factorize A = U V ⊤ where U, V ∈ R n×r are sparse matrices, and penalize γ( U 1 + V 1 ). The objective function to be minimized is

which is a non-convex function of the joint variable (U, V, W ), making the theoretical analysis more difficult. Given that the objective is convex in a neighborhood of the solution, by initializing the variables adequately, we can write an algorithm inspired by proximal gradient descent for minimizing it. The main motivation for directly optimizing the factors U and V of the matrix A rather than penalizing its nuclear norm come from the parallelization and stochastic optimization perspective such an approach offers. We leave further analysis for future work and settle for reporting empirical success over convex methods over simulated data in table 4.1.

Proof. We introduce the slack variable Z = A-A T and isolate the quadratic part of Ψ(Z 0 +Z, W 0 +W ) for some (Z 0 , W 0 ) :

Thanks to Cauchy-Schwarz and the basic norm property

We have Λ(Z) = D -Z. We get, again by Cauchy-Schwarz

On the other hand

νκ -4λ 2 √ n + 1 2 w 2 0 ≥ 0, and this completes the proof.

A generative model for graphs having linearly autoregressive features

Let V 0 ∈ R n×r be a sparse matrix, V † 0 its pseudo-inverse such, that

noise matrices N t and M t , which means that for any pair of indices (i, j), with high probability (N t ) i,j = 0 and (M t ) i,j = 0. We define the linear feature map ω(A) = AV ⊤ † 0 , and point out that

follows the linear autoregressive relation

2. For any time index t, the matrix A t is close to U t V 0 that has rank at most r 3. The matrices A t and U t are both sparse by construction.

Evaluation metrics

Various prediction tasks could be studied in our setup with specific criteria for optimization and evaluation. We focus here on regression (for feature prediction) and graph completion (for link prediction), but classification of vertices may also be a useful task. Here we denote by A the prediction we made for A T +1 given the past observations. 1. Regression. When the effective value of an asset in the future is of interest, a squared error metric is appropriate, leading to the use of ω(A T +1 ) -A F for evalution. The reported values are the relative errors

Our method aims simultaneousely at the prediction of ω(A T +1 ) and A T +1 . We measure the quality of prediction of A T +1 by A -A T +1 F , and the relative values reported are

3. Classification. Specific patterns may appear in the time series ω t which we may wish to predict. For instance, predicting top-selling items in a given market is definitely of interest. This problem can be formulated as follows: the matrix ω t represents the sales volumes over a market (each component corresponds to a product), and we assign a ±1 label depending on the order of magnitude of the increase of sales volumes over a specific time window.

Ranking

In cases where the ranking of all the instances are of interest for the application, the standard tool to measure the quality of a scoring rule is the use of Receiver O Characteristics (ROC) curves. The ROC curves is the plot of the True Positive Rate as a function of the False Positive Rate. The curves that are situated above others correspond to the best scoring rules. The area under the ROC curve is sometimes used to evaluate the quality of a scoring rule.

5. Top-k ranking In the recommender system applications for each user, a list of k = 10 recommendations are computed. We evaluate the performance with two different measures commonly used in Information Retrieval:

• the Hit Rate -evaluates if there is a correct hit in the list:

Hit Rate@N = # correct hits n ,

• the Normalized Discounted Cumulative Gain [KK02] -sensitive to the ranking of the item in the list:

and y i = 1 if the i-th recommendation was effectively purchased and y i = 0 otherwise.

Note that even though these two metrics seem very similar, they measure slightly different things. In fact NDCG is sensitive to the order in which the N recommendations are presented, while Hit Rate -that is more easily interpretable-considers this list as an unordered set.

Appendix A

Appendices

A.1 Mathematical tools and notations

A.1.1 Basic notions in graph theory

We present in this section usual notions and statistics used to capture topology and characteristics of a graph.

Standard graph features

Definition 21 (Degree) The degree of a node is the sum of the weights of its adjacent edges. If the graph is directed, we define the in degree and out degree for the sum of weights of edges directed to and from a vertex. If the graph is not weighted, we apply the definition by taking all the weights equal to 1.

Definition 22 (Path)

A path (or walk) is a sequence of vertices (v 1 . . . v n ), {v i } ∈ V such that two consecutive vertices are linked by an edge : ∀k ∈ {1, n -1}, (v k , v k+1 ) ∈ E. The length of a path is the number of edges n -1.

Definition 23 (Neighborhood)

The neighborhood of a vertex v, noted N (v), is the set of neighbors of v that is vertices adjacent to v not including v itself : N (x) = {y|(x, y) ∈ E}.

Definition 24 (Geodesic distance)

The geodesic distance between two nodes is the shortest path between the two nodes.

Definition 25 (Clique) A clique of a graph G = (V, E) is a subset C of the vertex set such that for every two vertices in C there exists an edge connecting the two :

Definition 26 (Clustering coefficient)

The clustering coefficient of a graph is the probability that 3 nodes belong to a triangle.

A "high" value of C i translates a disposition to establish direct links with the elements of of distance 2. By the same kind of consideration, one can quantify the inclination to attach to distance k elements.

Matrices related to a graph

Another way of normalizing the adjacency matrix is Remark 3 (Symmetric normalized adjacency) The symmetric normalized adjacency matrix is

Remark 4 (Graph laplacian) The unnormalized laplacian of a graph is the matrix

The normalized laplacian of a graph is the matrix

The following definition and the related remarks are specific to the bipartite case.

Definition 27 (Incidence matrix of a bipartite graph)

Remark 5 (Adjacency matrix of a bipartite graph)

The adjacency matrix of a bipartite graph is

Remark 6 (Projected graph adjacency matrices) The adjacency matrices of the projected graphs G 1 and G 2 are respectively

M T G M G Remark 7 (Symmetric normalized incidence) The symmetric normalized incidence matrix is

Think of a particle walking on the graph and define a probability map on the transitions (edges) as follows : given that the particle is located at the node i at time t, it has a uniform probability of moving to each of its neighbors at time t+1. It is easy to check that the stochastic matrix of such a random walk is given by the following Remark 8 (Random walk stochastic adjacency) The random walk stochastic matrix of a graph is D -1 A where D is an n × n diagonal matrix containing the degrees of the vertices.

Remark 9 (Random walk stochastic incidences)

The random walk stochastic matrix of the whole graph (related to the adjacency matrix A), is D -1 A where D is an (n + m) × (n + m) diagonal matrix containing the degrees of the vertices. The corresponding random walk incidence matrices are respectively

where D 1 and D 2 are respectively n × n and m × m diagonal matrices containing degrees of vertices on their diagonal.

Remark 10 (Random walk and symmetric normalized incidence)

and more generally

A.1.2 Singular Value Decomposition

The singular value decomposition of a matrix is a golden mine of information on its structure that has many powerful properties both from a mathematical and data-mining viewpoint. As we use it frequently in analysis and algorithms it is more than useful to remind some basic theoretical and computational aspects of it for a better understanding.

Theory

Theorem 13 (Singular Value Decomposition) For any M ∈ K n×m , with K = C or R, there exists a unique decomposition of the form

We refer to Σ i,i = σ i 's as singular values and to columns of U and V as singular vectors u i and v i . These notations allow to write the singular value decomposition using vector notations :

Proposition 6 (Rank r approximation) The closest rank r < n, m matrix to M in the sense of Frobenius norm, i.e. the solution to

is U ΣV * where Σ is a copy of Σ where the singular values are set to 0 for i > r : Σi,i = 0 for i > r.

The following shrinkage operator and the minimization problem it solves are of particular importance in our applications.

Definition 28 (Shrinkage Operator)

, the authors prove that

Computational issues

There are two main ways to simplify the SVD computation to a symmetric eigen-decomposition problem.

1. One can also obtain singular values of M by taking the square roots of eigen values of M M T or M T M . The singular vectors correspond to eigenvectors of the symmetric matrices M M T and M T M . Due to numerical error, matrix multiplication may hurt the smallest singular values precision. Hence for some applications, despite the higher storage requirement, one might prefer the following transformation.

Algorithm 8 Single-vector Lanczos recursion [LASVD]

Use a Lanczos algorithm to generate a family of real symmetric tridiagonal T r

for some k ≤ r compute the k relevant eigenvalues of T k that are approximations to the eigenvalues of A for each λ, compute a corresponding unit vector z such that T k z = λz, and compute the corresponding approximate eigenvector of A : y = V r z.

A.2 Proof of propositions: oracle bounds for regression with sparselow-rank matrices

We recall that P U denotes the projection associated with the orthogonal matrix U , P U = U U ⊤ . The projection P X is defined for a matrix X having a singular value decomposition X = U ΣV ⊤ by P X (B) = P U B + BP V -P U BP V .

It is the projector onto the linear space spanned by the matrices u k x ⊤ and yv ⊤ k for 1 ≤ j, k ≤ r and x, y ∈ R n . The projector onto the orthogonal space is given by P ⊥ A (B) = (I -P U )B(I -P V ). The notation Θ X is used for the sign pattern of the matrix X:

and Θ ⊥ X is the matrix having 1s where X i,j = 0 and 0s elsewhere.

A.2.1 Proof of proposition 3

We have for any X ∈ S, by optimality of X:

Thanks to trace-duality and ℓ 1 -duality we have M, X ≤ M ∞ X 1 and M, X ≤ M op X * for any X, so for any α ∈ [0, 1]:

now using assumptions τ ≥ 2α d M op , and γ ≥ 2(1-α) d M ∞ , and then triangle inequality

A.2.2 Proof of proposition 4.

We introduce the the convex function g : X → τ X * + γ X 1 , and let

denote an element of the subgradient of g, where W * op ≤ 1 and W 1 ∞ ≤ 1. There exist two matrices W 1 and W * such that

By optimality, an element of the subgradient of L at X belongs to the normal cone of S at X, we have ∂L( X), X -X ≤ 0 and on the other hand, by the monotonicity of the subgradient of the convex function g we have X -X, Z -Z ≥ 0. Therefore we can deduce

The inequality (A.1) can be written as

We need to bound ǫ, ω( X -X) . For this, note that by definition of M ,

We can write the following decomposition, that holds true for any real number α ∈ [0, 1],

We get by applying triangle inequality, Cauchy-Schwarz, Hölder inequality written for the trace-norm and ℓ 1 -norm

By using rank and support inequalities obtained by applying Cauchy-Schwarz inequality to P X (M ) F and Θ X • M F respectively, we get

we can rewrite the inequality (A.3) as follows:

the last inequality being due to the assumption τ ≥ 2α d M op , and γ ≥ 2(1-α) d M ∞ . So finally, and by using again this assumption,

The special case result and corollary results follows by using

and setting X = X 0 :

A.2.3 Proof of proposition 5

We recall the identity

then the bound trivially holds. Therefore lets assume ω( X -X 0 ), ω( X -X) > 0. In this case the inequality (A.3) implies

On the other hand for any α

By (A.7) and (A.8) we get

M ∞ , by using the fact that for x ≥ 3, x+2 x-2 ≥ 5, the following holds true

where β = γ τ . This proves that X -X ∈ C X,5,β , and therefore by definition of µ,

By using the inequality (A.4) we have

After noticing that by writing µ 5,β (X) = µ(X),

and similarly for the

So by bxx 2 ≤ b 2 2 we finally get

A.3 Proof of propositions: oracle bounds for prediction of feature autoregressive graphs

From now on, we use the notation (A, a) 2

Let us introduce the linear mapping Φ : R n×n × R d×d → R T ×d × R d given by

Using this mapping, the objective (4.3) can be written in the following reduced way:

Recalling that the error writes, for any A and W :

A.3.2 Proof of Theorem 10

Let A ∈ R n×n and W ∈ R d×d be fixed, and let A = U diag(σ 1 , . . . , σ r )V ⊤ be the SVD of A. Recalling that • is the entry-wise product, we have A = Θ A •|A|+Θ ⊥ A •A, where Θ A ∈ {0, ±1} n×n is the entry-wise sign matrix of A and Θ ⊥ A ∈ {0, 1} n×n is the orthogonal sparsity pattern of A. The definition (4.4) of ( Â, Ŵ ) is equivalent to the fact that one can find Ĝ ∈ ∂L( Â, Ŵ ) (an element of the subgradient of L at ( Â, Ŵ )) that belongs to the normal cone of A × W at ( Â, Ŵ ). This means that for such a Ĝ, and any A ∈ A and W ∈ W, we have Ĝ, ( Â -A, Ŵ -W ) ≤ 0.

(A.11)

Any subgradient of the function g(A) = τ A * + γ A 1 writes

). So, if Ẑ ∈ ∂g( Â), we have, by monotonicity of the sub-differential, that for any Z ∈ ∂g(A)

and, by duality, we can find Z such that

By using the same argument with the function W → W 1 and by computing the gradient of the empirical risk (A, W ) → R n (A, W ), Equation (A.11) entails that

(A.12) Using Pythagora's theorem, we have

Using Hölder's inequality, we obtain

and the same is done for | Θ W , Ŵ -W | ≤ Θ W • ( Ŵ -W ) 1 . So, when (A.14) holds, we obtain by rearranging the terms of (A.12):

Using Lemma 3, together with Hölder's inequality, we have for any α ∈ (0, 1):

(A.16) Now, using (A.15) together with (A.16), we obtain

which proves, using (4.7), that

In the same way, using (A.15) with A = Â together with (A.16), we obtain that Ŵ -W ∈ C 1 (W, 5). Now, using together (A.12), (A.13) and (A.16) , and the fact that the Cauchy-Schwarz inequality entails

and similarly for Ŵ -W , we arrive at

Since  -A ∈ C 2 (A, 5, γ/τ ) and Ŵ -W ∈ C 1 (W, 5), we obtain using Assumption 20 and ab ≤ (a 2 + b 2 )/2:

which concludes the proof of Theorem 10.

A.3.3 Proof of Theorem 12

For the proof of (4.8), we simply use the fact that 1 dT X T -1 ( Ŵ -W 0 ) 2 F ≤ E( Â, Ŵ ) 2 and use Theorem 11. Then we take W = W 0 in the infimum over A, W .

For (4.9), we use the fact that since Ŵ -W 0 ∈ C 1 (W 0 , 5), we have (see the Proof of Theorem 10),

and then use again Theorem 11. The proof of (4.10) follows exactly the same scheme.

A.3.4 Concentration inequalities for the noise processes

The control of the noise terms M and Ξ is based on recent developments on concentration inequalities for random matrices, see for instance [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF]. Moreover, the assumption on the dynamics of the features's noise vector (N t ) t≥0 is quite general, since we only assumed that this process is a martingale increment. Therefore, our control of the noise Ξ rely in particular on martingale theory.

Proposition 8 Under Assumption 3, the following inequalities hold for any x > 0. We have

(N T +1 ) j Ω j op ≤ σv Ω,op 2(x + log(2n)) d (A.17) with a probability larger than 1e -x . We have

with a probability larger than 1 -2e -x , and finally

with a probability larger than 1 -14e -x , where ℓ T = 2 max j=1,...,d log log

t=1 ω j (A t-1 ) 2 ∨ e .

Proof. For the proofs of Inequalities (A.17) and (A.18), we use the fact that (N T +1 ) 1 , . . . , (N T +1 ) d are independent (scalar) subgaussian random variables.

From Assumption 3, we have for any n × n deterministic self-adjoint matrices X j that E[exp(λ(N T +1 ) j X j )] exp(σ 2 λ 2 X 2 j /2) ,