Keywords: indexing invariance for VRGC

Since the introduction of Google's PageRank method for Web searches in the late 1990s, graph algorithms have been part of our daily lives. In the mid 2000s, the arrival of social networks has amplified this phenomenon, creating new use-cases for these algorithms. Relationships between entities can be of multiple types: user-user symmetric relationships for Facebook or LinkedIn, follower-followee asymmetric ones for Twitter or even user-content bipartite ones for Netflix, Deezer or Amazon. They all come with their own challenges and the applications are numerous: centrality calculus for influence measurement, node clustering for knowledge discovery, node classification for recommendation or embedding for link prediction, to name a few.

In the meantime, the context in which graph algorithms are applied has rapidly become more constrained. On the one hand, the increasing size of the datasets with millions of entities, and sometimes billions of relationships, bounds the asymptotic complexity of the algorithms for industrial applications. On the other hand, as these algorithms affect our daily lives, there is a growing demand for explanability and fairness in the domain of artificial intelligence in general. Graph mining is no exception. For example, the European Union has published a set of ethics guidelines for trustworthy AI. This calls for further analysis of the current models and even new ones.

This thesis provides specific answers via a novel analysis of not only standard, but also extensions, variants, and original graph algorithms. Scalability is taken into account every step of the way. Following what the Scikit-learn project does for standard machine learning, we deem it important to make these algorithms available to as many people as possible and participate in graph mining popularization. Therefore, we have developed an open-source software, Scikitnetwork, which implements and documents the algorithms in a simple and efficient way. With this tool, we cover several areas of graph mining such as graph embedding, clustering, and semi-supervised node classification.

Resumé

Depuis l'invention du PageRank par Google pour les requêtes Web à la fin des années 1990, les algorithmes de graphe font partie de notre quotidien. Au milieu des années 2000, l'arrivée des réseaux sociaux a amplifié ce phénomène, élargissant toujours plus les cas d'usage de ces algorithmes. Les relations entre entités peuvent être de multiples sortes : relations symétriques utilisateur-utilisateur pour Facebook ou LinkedIn, relations asymétriques follower-followee pour Twitter, ou encore, relations bipartites utilisateur-contenu pour Netflix ou Amazon. Toutes soulèvent des problèmes spécifiques et les applications sont nombreuses : calcul de centralité pour la mesure d'influence, le partitionnement de noeuds pour la fouille de données, la classification de noeuds pour les recommandations ou l'embedding pour la prédiction de liens en sont quelques exemples.

En parallèle, les conditions d'utilisation des algorithmes de graphe sont devenues plus contraignantes. D'une part, les jeux de données toujours plus gros avec des millions d'entités et parfois des milliards de relations limite la complexité asymptotique des algorithmes pour les applications industrielles. D'autre part, dans la mesure où ces algorithmes influencent nos vies, les exigences d'explicabilité et d'équité dans le domaine de l'intelligence artificielle augmentent. Les algorithmes de graphe ne font pas exception à la règle. L'Union européenne a par exemple publié un guide de conduite pour une IA fiable. Ceci implique de pousser encore plus loin l'analyse des modèles actuels, voire d'en proposer de nouveaux.

Cette thèse propose des réponses ciblées via l'analyse d'algorithmes classiques, mais aussi de leurs extensions et variantes, voire d'algorithmes originaux. La capacité à passer à l'échelle restant un critère clé. Dans le sillage de ce que le projet Scikit-learn propose pour l'apprentissage automatique sur données vectorielles, nous estimons qu'il est important de rendre ces algorithmes accessibles au plus grand nombre et de démocratiser la manipulation de graphes. Nous avons donc développé un logiciel libre, Scikit-network, qui implémente et documente ces algorithmes de façon simple et efficace. Grâce à cet outil, nous pouvons explorer plusieurs tâches classiques telles que l'embedding de graphe, le partitionnement, ou encore la classification semi-supervisée.

La thèse démarre par une introduction sur les structures de graphe ainsi qu'un certain nombre de rappels d'algèbre linéaire. Une attention particulière est accordée aux décompositions matricielles ainsi qu'aux méthodes stochastiques implémentées en pratique pour permettre de traiter les jeux de données massives. Dans un deuxième temps, cette thèse présente Scikitnetwork en détaillant les choix de design, les algorithmes implémentés ainsi que les performances obtenues vis-à-vis de logiciels concurrents.

Les chapitres qui suivent présentent les travaux originaux inclus dans Scikit-network. La première contribution théorique concerne la régularisation de l'embedding spectral obtenu par décomposition du Laplacien de graphe. On démontre sur un modèle simple qu'ajouter un graphe complet à la donnée d'origine permet à l'embedding spectral d'être moins sensible aux sous-graphes disjoints de la composante connexe principale. Dans un deuxième temps, on propose une variante de cet embedding pour les graphes dirigés et bipartis. On exhibe les propriétés mathématiques qui permettent d'interpréter facilement les représentations obtenues. L'ensemble de ces résultats est illustré par un grand nombre d'expériences sur données réelles.

On s'intéresse ensuite au célèbre algorithme de partition de noeuds Louvain. On démontre que cet algorithme peut être généralisé à toute une classe de fonctions de modularité et non pas seulement celle proposée par Newman. On vérifie la pertinence des diverses modularités sur un ensemble de graphes réels de tous types : non-dirigés, dirigés et bipartis.

Dans une partie suivante, cette thèse propose une modification simple d'un célèbre algorithme de classification de noeuds semi-supervisé introduit par Zhou. On démontre sur un modèle simple en quoi l'algorithme d'origine est biaisé envers la classe majoritaire dans les seeds et que la modification proposée assure bien une classification correcte. Cette amélioration est confirmée par la suite sur données simulées et réelles.

Enfin, la thèse se termine par quelques contributions sur les algorithmes de classification de graphes. On compare une méthode basée sur la décomposition spectrale du Laplacien à une autre qui combine parcours de graphe et réseau de neurones récurrents. Les performances sont mises en perspective de l'état de l'art dans le domaine afin de démontrer la pertinence des approches proposées. 6.5.2 General benchmark . . . . . . . . . .

Chapter 1 Introduction

Graphs

Graphs are mathematical objects describing the relations between entities called vertices or nodes. If two nodes are in relation, they are said to be connected by an edge. Mathematically, we denote by G = (V, E) a graph with V the set of vertices, which we assume finite, and E ⊂ V × V, the set of edges.

In this work, we consider an extension of the previous definition where weights are associated with edges. These objects are called weighted graphs, but we simply refer to them as graphs.

Let n be the number of nodes of a graph, n = |V|.

The adjacency matrix associated with the graph is the matrix A ∈ R n×n such that

A ij = ω if (i, j) ∈ E, 0 otherwise.
ω being the weight associated to the edge (i, j).

Besides, we differentiate three main families of graphs which we illustrate in Figure 1.1. Undirected graphs In these graphs, relationships are symmetric. This means that if (i, j) ∈ E, then (j, i) ∈ E as well. This is for example the case for most social networks such as Facebook friendship graph or LinkedIn professional relationships. In particular, this implies that the associated adjacency matrix is symmetric.

Directed graphs In these graphs, the relationships are not necessarily symmetric, i.e. (i, j) ∈ E does not imply (j, i) ∈ E. This is the case for the Twitter graph in which a user can follow another one without being followed in return. Another example is the graph of hyperlinks between Web pages. In this case, A is not necessarily symmetric. We use the terms source node and target node to differentiate both ends of directed edges.

Bipartite graphs These are a special case of undirected graphs with two disjoint sets of nodes V 1 and V 2 such that relationships only exist between V 1 and V 2 . On the other hand, two nodes of V 1 cannot be connected nor can two nodes of V 2 . This is the case for relationships such as user-content in Netflix or user-item in retail. We denote by B the biadjacency matrix, possibly rectangular, such that the adjacency matrix of the graph seen as undirected is

A = 0 B B 0 .
Note that this list is not exhaustive. There are many other types of graphs such as knowledge graphs, signed graphs, hypergraphs, annotated or infinite graphs... However, they lie beyond the scope of the present work.

To conclude this section, we introduce a few notations that will help further reading. Unless otherwise specified:

• n denotes the number of nodes in a graph.

• m denotes the number of non-zero entries in an adjacency matrix. For an undirected graph without self-loops, this is twice the number of edges.

• 1 denotes a vector of ones whose dimension depends on context.

• d is the out-degree vector defined by A1 or B1. If the graph is not weighted, the out-degree d i of a given node i is simply the number of edges for which i is the source. If the edges are weighted, d i is the total weight of edges for which i is the source.

• D is the out-degree diagonal matrix defined by D = diag(d).

• f is the in-degree vector defined by A 1 or B 1. If the graph is not weighted, the indegree f i of a given node i is simply the number of edges for which i is the target. If the edges are weighted, f i is the total weight of edges for which i is the target.

• F is the in-degree diagonal matrix defined by F = diag( f ).

• w is the total weight of the graph defined by 1 A1 or 1 B1.

Graph Mining

The notion of graph mining covers a wide variety of tasks on graphs, most of which have direct industrial applications. Here, we provide a non-exhaustive list of examples.

Node embedding This task consists in mapping each node of the graph to a vector, see Figure 1.2. Usually, the embedding space has a small dimension with respect to the number of nodes in the graph. Embeddings are designed such that the distances between nodes in the vector space reflect some topological properties of the original graph. The embedding can be used to perform standard machine learning tasks such as clustering or classification or link prediction as in the famous Netflix Prize [START_REF] Bell | The bellkor solution to the netflix prize[END_REF]. Chapter 3 is dedicated to the study of a widely used embedding: the spectral embedding. Node clustering This task consists in assigning each node to a group or cluster such that nodes in the same cluster are more connected with one another than with the rest of the graph. See Figure 1.3 for an illustration and [START_REF] Schaeffer | Graph clustering[END_REF] for a review of some algorithms. Applications include load balancing for wireless sensors networks [START_REF] Younis | Node clustering in wireless sensor networks: Recent developments and deployment challenges[END_REF], anomaly detection for cyber-security or financial fraud detection [START_REF] Akoglu | Graph based anomaly detection and description: a survey[END_REF]. Node clustering is the subject of Chapter 4. Node classification This task consists in assigning a label to each node of the graph based on some partial ground-truth information, see Figure 1.4 for an illustration. It can be used to assign a polarity to tweets [START_REF] Speriosu | Twitter polarity classification with label propagation over lexical links and the follower graph[END_REF], predict potentially dangerous drug-drug interactions [START_REF] Zhang | Label propagation prediction of drug-drug interactions based on clinical side effects[END_REF] or automatically assign a genre to a movie or a song in a database. We discuss node classification algorithms in Chapter 5. Node ranking This task consists in assigning an importance score to each node of the graph with respect to a certain criterion, see Figure 1.5 for an illustration. This is particularly useful for recommender systems such as Web search engines [START_REF] Page | Method for node ranking in a linked database[END_REF] or contact recommendation in social networks. Such methods can also be used in a pipeline for node classification as we will see in Chapter 5.

Graph classification

This task is the same as the standard supervised classification in machine learning with the main difference that the samples are graphs instead of vectors. Standard applications are molecular compounds classification in biology [START_REF] Childs | Identification and classification of ncrna molecules using graph properties[END_REF][START_REF] Kudo | An application of boosting to graph classification[END_REF] or malware detection in cyber-security [START_REF] Hu | Large-scale malware indexing using function-call graphs[END_REF][START_REF] Gascon | Structural detection of android malware using embedded call graphs[END_REF][START_REF] Kinable | Malware classification based on call graph clustering[END_REF]. We discuss some graph classification algorithms in Chapter 6.

Link prediction Among other graph mining tasks, we deem it important to mention edge prediction which consists in assigning a probability of existence to unobserved relations among nodes. However, we do not tackle this task in the present work.

Data structures

In order for graph algorithms to scale, they must rely on efficient data structures. This section introduces some standard graph formats. We present some of their respective strengths and weaknesses with a special emphasis on the Compressed Sparse Row format (CSR) which is used in most of our implementations. We refer the reader to SciPy's documentation1 for a more extensive list.

Sparse matrix formats

Let G = (V, E) be a graph with n nodes and associated adjacency matrix A ∈ R n×n such that m is the number of nonzero entries in A. Storing G boils down to storing A. The main assumption we make is that A is sparse, which means that it has very few nonzero coefficients: m n 2 . In other words, the average number of neighbors of a node is small compared to n. For real large graphs, it is common to have densities, i.e. m/n 2 ratios, below 10 -5 (see Table 1.4). Hence, these adjacency matrices are made of more than 99.99% of zeros which it is useless to store explicitly in memory. Table 1.1 presents an overview of some standard formats which we detail below. format memory space linear algebra friendly flexible COO 2mM int + mM edge CSR/CSC (n + m)M int + mM edge DOK mM key Table 1.1: Some standard formats for adjacency matrices. M int , M edge and M key denote the memory space of an integer, an edge and a key of the hash table, respectively.

Coordinates (COO)

A first natural way to store the adjacency matrix of a graph is to store the corresponding edges of the graph (i, j, ω) where i is the source node, j is the target node and ω is the associated weight (possibly Boolean, float or integer). In this format, the edges are stored in no particular order which has some implications for basic graph primitives.

Compressed Sparse Row or Column (CSR, CSC)

The CSR format is composed of three arrays, an array of size m named indices, which is the concatenation of the nodes successors, an array of size m named data to store the associated edges weights and an index pointer, indptr, of size n + 1 indicating which section of indices corresponds to which node in the graph. Example 1.3.1 illustrate this format on the house graph introduced in Figure 1.1a. See that node 0, top of the roof, has two neighbors, nodes 1 and 4. The CSC format works the same way but the indices array indicates the nodes predecessors.

Example. 1.3.1: CSR format

[1]: from sknetwork.data import house A = house() print(A.indptr) print(A.indices)

[1]: [ 0 2 5 7 9 12] [1 4 0 2 4 1 3 2 4 0 1 3] Dictionary of keys (DOK) As suggested by the name, this format maps tuples of indices (i, j) to keys of a dictionary with associated edges weights (possibly Boolean). It is effective for values access or updates but the hash table is usually more expensive in memory usage than other formats.

Performance analysis

We illustrate the performance of the different graph formats for a few standard graph primitives on a real dataset, Wikivitals (WV). The graph represents hyperlinks between some pages of Wikipedia and the adjacency matrix has 10 012 rows and 792 091 Boolean entries. The pages are labeled with their Wikipedia category (Art, Science, History...) which allows to later use supervised metrics on this dataset in Chapters 4 and 5. Note that we made this dataset available online2 and via Scikit-network as explained in Chapter 2.

The test tasks are the following:

• Memory space: Make a pickle dump of the adjacency matrix and record the memory space.

• Node insertion: Add one isolated node i.e. add one extra row and one extra column with a diagonal entry to the adjacency matrix.

• Node removal: Remove the first row and first column of the matrix.

• Edge insertion: Add a link between nodes 0 and 1 (which does not exist in the original graph).

• Edge removal: Remove the previously created link between nodes 0 and 1 and eliminate the explicit 0 value in the matrix.

• Out-neighbors access: Get all successors of node 0.

• In-neighbors access: Get all predecessors of node 0.

We repeat each task involving a runtime measurement a hundred times for each format. Average values are reported in Table 1.2. If a format requires to be converted to perform a given operation, we do not record the runtime and replace it by "-".

Setups Unless otherwise specified, we always use one of the following setups for our experiments which we refer to as laptop setup and server setup:

• Laptop: Laptop equipped with an Intel(R) Core(TM) i7-7820HQ CPU @ 2.90GHz processor and 16 Go of RAM.

• Server: Computer equipped with an AMD Ryzen Threadripper 1950X 16-Core Processor and 32 Go of RAM. For deep learning experiments, the computer is equipped with an NVIDIA TitanXp GPU.

In this case, all experiments are run on the laptop setup.

The conclusion of these experiments is that, as long as the structure of the matrix does not need to be modified, the CSR and CSC formats are both lighter in memory and more efficient. On the other hand, the DOK format is relevant when edges need to be modified often (in a dynamic graph, for example). In the rest of this thesis, the default format of adjacency matrices is the CSR one, unless otherwise specified.

Linear algebra for graph mining

This section covers some standard linear algebra primitives in the context of graph mining. Part of this work has been published in [START_REF] Lara | The sparse+ low rank trick for matrix factorization-based graph algorithms[END_REF].

Transposition

Depending on whether an algorithm requires access to the successors or the predecessors of the nodes in the graph, it might need to transpose the adjacency matrix. This operation is immediate in the COO format as it only requires to switch the row and col vectors. However, this is not the case for the CSR format.

The canonical way to compute A in CSR format from A in CSR format requires an intermediate transformation into COO format. Going from CSR to COO only requires to compute the row vector from the indptr one which is an operation linear in m. However, going from COO to CSR requires a sort of the indices which is O(m log(m)) operations. In the end, the complexity of CSR matrix transposition is O(m log(m)).

Matrix-vector product

Matrix-vector product is a key primitive for many graph algorithms. It can be used to emulate discrete time random walks, perform low dimension projections, message passing... It is therefore critical to perform it efficiently.

Respective pseudocodes for COO and CSR formats are presented in Example 1.4.1. In both cases, the complexity is linear with the number of nonzero coefficients in the matrix. The main difference is that, as indices are sorted in the CSR, the processor can optimize caching operations.

Example. 1.4.1: Sparse matrix-vector products

[1]: import numpy as np def coo_dot(n, row, col, data, x): y = np.zeros(n) for i, w in enumerate(data):

y[row[i]] += w * x[col[i]] return y def csr_dot(indptr, indices, data, x): n = len(indptr) -1 y = np.zeros(n) for i in range(n): start, end = indptr[i], indptr[i+1] y[i] = data[start:end].dot(x[indices[start:end]]) return y
There are several more complex sparse matrix formats designed to optimize matrix-vector product, especially in the context of parallel computing and GPU computations. See [START_REF] Liu | Efficient sparse matrix-vector multiplication on x86-based many-core processors[END_REF], [START_REF] Bell | Efficient sparse matrix-vector multiplication on cuda[END_REF] for more details.

Eigenvalue decomposition

Let A be a symmetric matrix with real coefficients. The spectral theorem states that there exits an orthonormal matrix X and a diagonal matrix Λ with real coefficients such that

A = XΛX .
The columns of X are the eigenvectors of A while the diagonal entries of Λ are the associated eigenvalues.

Several graph algorithms rely on the computation of some of these eigenvectors and eigenvalues. In this case, we refer to this as a partial eigenvalue decomposition.

In practice, iterative power methods are used to compute the k leading eigenvectors, i.e. the eigenvectors associated with the k largest eigenvalues. The idea is to start from a random vector x 0 and compute the sequence (x t ) defined by

x t+1 = Ax t ||Ax t || .
The sequence (x t ) provably converges with probability 1 towards the leading eigenvector of A, provided that it is unique. Based on this result, the Arnoldi method [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF] uses the stabilized Gram-Schmidt process to compute the k leading eigenvectors by orthonormalization of x t , x t-1 , . . . , x t-k+1 . Note that the Arnoldi method applies to non-Hermitian matrices. In the Hermitian case, it is called the Lanczos method [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF].

Observe once again that the main primitive for this method is the matrix-vector product.

An extension of the eigenvalue decomposition is the generalized eigenvalue decomposition. Given another matrix W, the generalized eigenvalue of A is such that

AX = WXΛ.
Proposition 1. Assume that W is a diagonal matrix with positive entries and let X denote some eigenvectors for the eigenvalue decomposition of W -1/2 AW -1/2 . Then, X = W -1/2 X are eigenvectors for the generalized eigenvalue decomposition of A.

Proof. Let X be eigenvectors of W -1/2 AW -1/2 such that

W -1/2 AW -1/2 = XΛX . Then, AW -1/2 X = W 1/2 XΛ(X X), AX = WX Λ.
Such that X are generalized eigenvectors of A.

Singular value decomposition

The singular value decomposition (SVD) of a real matrix M ∈ R n 1 ×n 2 is a factorization of the form M = UΣV , where U ∈ R n 1 ×n 1 and V ∈ R n 2 ×n 2 are orthonormal matrices whose columns are respective left and right singular vectors and Σ ∈ R n 1 ×n 2 is a diagonal matrix containing the nonnegative singular values.

Note that, U and V are respective eigenvectors of MM and M M, the associated eigenvalues being the square of the singular values. Besides, a partial SVD, extracting only the k leading singular vectors couples is an optimal rank k approximation of M with respect to the Froebenius norm:

U (k) Σ (k) V (k)T ∈ arg min rank(X)=k ||M -X|| F .
As for the eigenvalue decomposition, this factorization admits a generalization. Given two matrices W 1 and W 2 , the generalized singular value decomposition of M satisfies

MV = W 1 UΣ, M U = W 2 VΣ.
Proposition 2. Assume that W 1 and W 2 are diagonal matrices with positive entries and let U, V be left and right singular vectors of W

-1/2 1 MW -1/2 2 . Then, U = W -1/2 1 U and V = W -1/2 2 V are singular vectors for the generalized eigenvalue decomposition of M. Proof. Let U, V be singular vectors of W -1/2 1 MW -1/2 2 such that W -1/2 1 MW -1/2 2 = UΣV . Then, MW -1/2 2 V = W 1/2 1 UΣ(V V), MV = W 1 U Σ.
Similar derivations lead to the second equation such that U , V are generalized singular vectors of M.

Halko's randomized method

Even with an efficient matrix-vector product, computing partial decompositions of very large matrices can be computationally prohibitive. In this section we present Halko's randomized projection method [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] which is a good solution to this issue.

The idea is to compute the decomposition of another matrix, much smaller than the original one and approximate the desired eigenvectors or singular vectors through a linear transformation.

First, one must find an orthonormal projection matrix

Q ∈ R n×p such that M ≈ QQ M.
Algorithm 1.1 offers to do so using a combination of random projection, power iteration and QR decompositions. Note that, once again, the only two actions on the sparse matrix M are matrix-vector products and transposition. In practice, p has to be chosen between k + 7 and 2k

where k is the number of desired eigenvectors or singular vectors. As k is usually very small compared to n, the dimension p of the rows of Q is also small. The power iteration parameter q has to be a small integer as well, typically between 3 and 10.

Algorithm 1.1 Randomized Subspace Iteration

Require

: input matrix M ∈ R n×n Generate a Gaussian test matrix Ω ∈ R n×p . Form Y 0 = MΩ ∈ R n×p and compute its QR factorization Y 0 = Q 0 R 0 , with Q 0 ∈ R n×p and R 0 ∈ R p×p . for i = 1, . . . , q do Form Ỹi = M Q i-1 and compute its QR factorization Ỹi = Qi Ri . Form Y i = M Qi and compute its QR factorization Y i = Q i R i . end for Q = Q q .
Once the matrix Q is computed, one can perform an eigenvalue decomposition of either

M = Q MQ ∈ R p×p for eigenvalue decomposition of M, or M = Q MM Q ∈ R p×p
for its singular value decomposition. As p n, such decompositions are much less computationally expensive than a direct decomposition of M. Now, let x ∈ R p be an eigenvector of M associated with the eigenvalue λ. Then, x = Q x ∈ R n approximates an eigenvector of M as

Mx = MQ x ≈ QQ MQ x = QM x = Qλ x = λx.
The same thing goes for ū, v left and right singular vectors of M associated with the eigenvalue

σ 2 . Let u = Q ū and v = Q v, then MM v = MM Q v ≈ QQ MM Q v = QM v = Qσ 2 v = σ 2 v.
So that v is an approximate eigenvector of MM associated with the eigenvalue σ 2 and consequently an approximate right singular vector of M associated with the singular value σ. Similar derivations show that u is an approximate left singular vector of M.

A classic use-case for Halko's algorithm is the implementation of large-scale Principal Component Analysis [START_REF] Yu | Single-pass pca of large high-dimensional data[END_REF]. We present some other ones in Section 1.4.6 and in Chapter 3.

Linear operators

A linear operator of R n is a function f such that for u, v ∈ R n and λ ∈ R, f (u + λv) = f (u) + λ f (v).
In computer science, they are defined by two primitives, namely matvec (for matrix-vector product) and transpose which would respectively be f and its conjugate f * in mathematical terms. Recall that the i th column of a matrix is simply the image of the i th vector of the canonical basis by the associated linear operator.

These structures are of great interest in graph mining for memory reasons. Indeed, they do not require to explicitly store the images of the canonical vectors. As previously stated, the matrixvector product and the transposition are enough to actually compute eigenvalue and singular value decompositions.

As a matter of fact, several matrices of interest in graph mining have an explicit "sparse + lowrank" structure, which means they are of the form

M = S + xy ,
where S is a sparse matrix and x, y are vectors. Such matrices are dense, however, storing them as linear operators requires only to store 2m + 2n coefficients (as S might have to be stored as well depending on the use-case). The matrix-vector product being computed as follows:

Mz = Sz + (y z)x.
The operator corresponding to M can be written (S, (x, y)), such that the operator corresponding to M is simply S , (y, x) .

For each one of the following examples from the literature, Table 1.3 identifies the sparse term S and the low-rank term (x, y).

matrix

S (sparse)

x y

Random surfer αD

-1 A 1-α n 1 1 Modularity 1 w A -γ w d 1 w d Laplacian I -D -1/2 α AD -1/2 α -α n √ d + α1 √ d + α1 PCA A 1 w f 1
Table 1.3: Identification of sparse and low-rank factors. (D α = D + αI.)

• Random surfer's transition matrix [Page et al., 1999a],

R ij = αA ij d i + 1 -α n ,
where α is a float between 0 and 1 and n is the number of nodes in the graph.

• The modularity matrix [START_REF] Newman | Modularity and community structure in networks[END_REF],

M ij = 1 w A ij - γ w 2 d i d j ,
where γ ≥ 0 is a resolution parameter.

• The regularized graph Laplacian [START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF],

(L α ) ij = δ ij - A ij + α/n (d i + α)(d j + α) ,
where δ is the Kronecker symbol such that δ ij = 1 if i = j and 0 otherwise and α ≥ 0.

• The principal components analysis (PCA) matrix [START_REF] Saerens | The Principal Components Analysis of a Graph, and Its Relationships to Spectral Clustering[END_REF],

Āij = A ij -f j /w.

Experiments

In this section, we illustrate the use of linear operators combined to Halko's algorithm on real datasets.

We compute the partial eigenvalue decomposition and the partial singular value decomposition of the modularity matrix of 7 undirected and 7 bipartite graphs, respectively. In the bipartite case, the modularity matrix is defined as

Q = 1 w B - γ w 2 d f ,
as later discussed in Chapter 4.

The graphs are collected from the Konect database [START_REF] Kunegis | KONECT: The Koblenz Network Collection[END_REF] so that they have a varying number of nodes and edges. In a first experiment, we extract k = 16 components in each case and record the computation times for our modified Python implementation of Halko's method. We compare our results to the computation times of the direct Lanczos method using SciPy [START_REF] Virtanen | Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. Results are displayed in Figure 1.6.

In a second experiment, we record the respective running times for the Reactome (undirected) and MovieLens user-tag (bipartite) graphs for k ∈ {8, 16, 32, 64, 128, 256}. See Figure 1.7.

Results

In the laptop setup, memory errors arise around n 2 or n 1 n 2 = 10 9 for the dense implementation. At this point, the Lanczos method is already between one and two orders of magnitude slower than Halko's. On the other hand, the "sparse + low-rank trick" enables to handle graphs with more than 3 million nodes and 234 million edges in about five minutes.

Furthermore, note that the running time for the implicit dense matrix is quite close to the running time for its sparse component. This holds for a fixed k and varying graphs as well as for fixed graphs and varying k.

Contributions

In this chapter, we have introduced graphs, their related tasks and data structures as well as some useful linear algebra that will help further reading [START_REF] Lara | The sparse+ low rank trick for matrix factorization-based graph algorithms[END_REF]. The rest of this thesis is organized as follows:

In Chapter 2, we introduce Scikit-network, a Python software for graph mining which is used in almost all our experiments.

Chapter 3 presents two distinct contributions to spectral embedding. First, a novel analysis of a standard regularization technique [Lara and Bonald, 2020]. Second, an interpretable extension to bipartite and directed graphs based on the Generalized Singular Value Decomposition [START_REF] Bonald | The forward-backward embedding of directed graphs[END_REF].

Chapter 4 is dedicated to node clustering. We propose a general framework to interpret modularity-like quality functions. Specifically, we make a novel connection between modularity maximization and block model fitting. Then, we detail how the Louvain heuristic can be generalized to find local optimums for all these quality functions.

In Chapter 5, we propose a novel algorithm for semi-supervised node classification based on heat diffusion. We highlight its efficiency with respect to other diffusion-based methods by proving its consistency on a simple model, then, we assess its performance on real graphs.

In Chapter 6, we discuss graph classification. We propose two distinct algorithms to address this problem [de [START_REF] Lara | The sparse+ low rank trick for matrix factorization-based graph algorithms[END_REF][START_REF] Lara | A Simple Baseline Algorithm for Graph Classification[END_REF]de Lara, 2019b] and benchmark them on standard datasets.

Finally, Chapter 7 concludes the thesis.

A schematic structure of the thesis is displayed in Figure 1.8.

Publications

This section enumerates the publications related to the present thesis in chronological order.

Chapter 2

Scikit-network

Scikit-network (Figure 2.1) is a Python package inspired by Scikit-learn for the analysis of large graphs [START_REF] Bonald | Scikit-network: Graph analysis in python[END_REF]. Graphs are represented by their adjacency matrix in the sparse CSR format of SciPy. The package provides state-of-the-art algorithms for ranking, clustering, classifying, embedding and visualizing the nodes of a graph. High performance is achieved through a mix of fast matrix-vector products (using SciPy), compiled code (using Cython) and parallel processing. The package is distributed under the BSD license, with dependencies limited to NumPy and SciPy. It is compatible with Python 3.6 and newer. Source code, documentation and installation instructions are available online1 . This chapter describes the version 0.18.0 of the package. Journal of Machine Learning Research.

Motivation

Scikit-learn [Pedregosa et al., 2011b] is a machine learning package based on the popular Python language. It is well-established in today's machine learning community thanks to its versatility, performance and ease of use, making it suitable for both researchers, data scientists and data engineers. Its main assets are the variety of algorithms, the performance of their implementation and their common API.

Scikit-network is a Python package inspired by Scikit-learn for graph analysis. The sparse nature of real graphs, with up to millions of nodes, prevents their representation as dense matrices and rules out most algorithms of Scikit-learn. Scikit-network takes as input a sparse matrix in the CSR format of SciPy and provides state-of-the-art algorithms for ranking, clustering, classifying, embedding and visualizing the nodes of a graph.

The design objectives of scikit-network are the same as those having made Scikit-learn a success: versatility, performance and ease of use. The result is a Python-native package, like NetworkX [START_REF] Hagberg | Exploring network structure, dynamics, and function using networkx[END_REF], that achieves the state-of-the-art performance of iGraph [START_REF] Csardi | The igraph software package for complex network research[END_REF]] and graph-tool [START_REF] Peixoto | The graph-tool python library[END_REF] (see the benchmark in section 2.5). Scikit-network uses the same API as Scikit-learn, with algorithms available as classes with the same methods (e.g., fit). It is distributed with the BSD license, with dependencies limited to NumPy [START_REF] Walt | The numpy array: a structure for efficient numerical computation[END_REF] and SciPy [START_REF] Virtanen | Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF].

Software features

The package is organized in modules with consistent API, covering various tasks:

• Data. Module for loading graphs from distant repositories, including Konect [START_REF] Kunegis | KONECT: The Koblenz Network Collection[END_REF], parsing tsv files into graphs, and generating graphs from standard models, like the stochastic block model [Airoldi et al., 2008].

• Clustering. Module for clustering graphs, including a soft version that returns a nodecluster membership matrix. It currently implements several variations of the Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], as a pipeline to apply k-means on a graph embedding or even a label propagation [START_REF] Raghavan | Near linear time algorithm to detect community structures in large-scale networks[END_REF].

• Hierarchy. Module for the hierarchical clustering of graphs, returning dendrograms in the standard format of SciPy. The module also provides various post-processing algorithms for cutting and compressing dendrograms. Implemented algorithms include Paris [Bonald et al., 2018a] and a pipeline to apply Ward's algorithm [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] on a graph embedding.

• Embedding. Module for embedding graphs in a space of low dimension. This includes spectral embedding and standard dimension reduction techniques like SVD and GSVD, with key features like regularization. This module also implements the spring layout [START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF] for visualization purposes.

• Ranking. Module for ranking the nodes of the graph by order of importance. This includes PageRank [Page et al., 1999b] and various centrality scores such as HITS [START_REF] Kleinberg | Authoritative sources in a hyperlinked environment[END_REF], diffusion based-ranking [START_REF] Chung | The heat kernel as the pagerank of a graph[END_REF], closeness or harmonic centrality [START_REF] Marchiori | Harmony in the small-world[END_REF]].

• Classification. Module for classifying the nodes of the graph based on the labels of a few nodes (semi-supervised learning). Some of the implemented algorithms are evaluated in Chapter 5.

•

Project Assets

Code quality and availability. Code quality is assessed by standard code coverage metrics. This version's coverage is greater than 98% for the whole package. Requirements are also kept up to date thanks to the PyUp tool. Scikit-network relies on TravisCI for continuous integration and cibuildwheel and manylinux for deploying on common platforms. OSX, Windows 32 or 64-bit and most Linux distributions [START_REF] Mcgibbon | A platform tag for portable linux built distributions[END_REF] are supported for Python versions 3.6 and newer.

Open-source software. The package is hosted on GitHub2 and part of SciPy kits aimed at creating open-source scientific software. Its BSD license enables maximum interoperability with other software. Guidelines for contributing are described in the package's documentation3 and guidance is provided by the GitHub-hosted Wiki. Documentation. Scikit-network is provided with a complete documentation 3 . The API reference presents the syntax while the tutorials present applications on real graphs. Algorithms are documented with relevant formulas, specifications, examples and references, when relevant.

Code readability. The source code follows the stringent PEP8 guidelines. Explicit variable naming and type hints make the code easy to read. The number of object types is kept to a minimum.

Data collection. The package offers multiple ways to fetch data. Some small graphs are embedded in the package itself for testing or teaching purposes. Part of the API makes it possible to fetch data from selected graph databases easily. Parsers are also present to enable users to import their own data and save it in a convenient format for later reuse.

Resources

Scikit-network relies on a very limited number of external dependencies for ease of installation and maintenance. Only SciPy and NumPy are required on the user side.

SciPy. Many elements from SciPy are used for both high performance and simple code. The sparse matrix representations allow for efficient manipulations of large graphs while the linear algebra solvers are used in many algorithms. Scikit-network also relies on the LinearOperator class for efficient implementation of certain algorithms.

NumPy. NumPy arrays are used through SciPy's sparse matrices for memory-efficient computations. NumPy is used throughout the package for the manipulation of arrays. Some inputs and most of the outputs are given in the NumPy array format.

Cython. In order to speed up execution times, Cython [START_REF] Behnel | Cython: The best of both worlds[END_REF] generates C++ files automatically using a Python-like syntax. Thanks to the Python wheel system, no compilation is required from the user on most platforms. Note that Cython has a built-in module for parallel computing on which scikit-network relies for some algorithms. Otherwise, it uses Python's native multiprocessing.

Performance

To show the performance of scikit-network, we compare the implementation of some representative algorithms with those of the graph softwares of Table 2.1: the Louvain clustering algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], PageRank [Page et al., 1999b], HITS [START_REF] Kleinberg | Authoritative sources in a hyperlinked environment[END_REF] and the spectral embedding [Belkin and Niyogi, 2002c]. For PageRank, the number of iterations is set to 100 when possible (that is for all packages except iGraph). For Spectral, the dimension of the embedding space is set to 16. Table 2.2 gives the running times of these algorithms on the Orkut graph of Konect [START_REF] Kunegis | KONECT: The Koblenz Network Collection[END_REF]

The NetSet repository

As mentioned in Section 2.2, Scikit-network allows to load datasets from some repositories. Among those is the NetSet4 repository which has been implemented by the Scikit-network team.

As in Scikit-learn, datasets are stored as Bunch objects with consistent attribute names. Adjacency and biadjacency matrices are stored in the compressed npz format of SciPy which is highly efficient both in memory space and loading speed. Other attributes such as labels or names of the nodes are stored in the npy format of NumPy.

The repository currently gathers some standard networks from the literature such as Cora or the 20newsgroup, but also some original ones such as Wikipedia Vitals, already mentioned in Chapter 1. See Example 2.6.1 for an illustration.

Example. 2.6.1: Loading a graph from the NetSet repository

[1]: from sknetwork.data import load_netset dataset = load_netset('cora') dataset

[1]: {'name': 'cora', 'labels': array([2, 5, 4, ..., 1, 0, 2]), 'adjacency': <2708x2708 sparse matrix of type '<class 'numpy.

→float64'>'

with 10556 stored elements in Compressed Sparse Row format>, 'meta': {'name': 'cora'}} Chapter 3

Spectral embedding

Learning from network-structured data often requires to embed the underlying graph in some Euclidian space [START_REF] Yan | Graph embedding and extensions: A general framework for dimensionality reduction[END_REF], [START_REF] Grover | node2vec: Scalable Feature Learning for Networks[END_REF], [START_REF] Bronstein | Geometric deep learning: going beyond euclidean data[END_REF]. The quality of the subsequent learning tasks (e.g., clustering, classification) then critically depends on this embedding, that must reflect the graph structure.

Usual graph embedding techniques rely on the heuristic that nodes which are close in the graph according to some similarity metric should also be close in the embedding space [START_REF] Cai | A comprehensive survey of graph embedding: Problems, techniques, and applications[END_REF]. Standard similarity metrics include first-order proximity (i.e., edges) [Belkin and Niyogi, 2002c], second-order proximity [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Dhillon | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF], Tang et al., 2015a] or higher-order proximity (i.e., random walks) [START_REF] Grover | node2vec: Scalable Feature Learning for Networks[END_REF][START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF].

The embedding is then the solution to some optimization problem, which is solved either by matrix factorization [Belkin and Niyogi, 2002c[START_REF] Dhillon | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Wold | Principal component analysis[END_REF] or by iterative methods based on skip-gram negative sampling [START_REF] Grover | node2vec: Scalable Feature Learning for Networks[END_REF][START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF], Tang et al., 2015a]. However, even for SVD-based methods [START_REF] Dhillon | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF][START_REF] Wold | Principal component analysis[END_REF], there is most of the time no clear interpretation of the distance between nodes in the embedding space.

A classical embedding of undirected graphs is based on the spectral decomposition of the Laplacian [Belkin andNiyogi, 2002c, Ng et al., 2002]; after proper normalization, the distance between nodes in the embedding space corresponds to the mean commute time of a random walk in the graph, making this embedding meaningful and easy to interpret [START_REF] Qiu | Clustering and embedding using commute times[END_REF]Hancock, 2007, Fouss et al., 2007].

• [START_REF] Bonald | The forward-backward embedding of directed graphs[END_REF]. The Forward-Backward Embedding of Directed Graphs.

Openreview.

• De Lara, N., [START_REF] Hollocou | Modularity-based Sparse Soft Graph Clustering[END_REF]. Spectral embedding of regularized block models.

International Conference on Learning Representations.

Problem statement

In Belkin and Niyogi [2002c], the embedding is presented as the minimizer of the following loss function:

X ∈    arg min L = ∑ i,j A ij ||X i -X j || 2 , X WX = I,
where W can be either I or D = diag(A1). See that the loss function L penalizes nodes which are connected but distant in the embedding space while the orthonormality constraint rules out the trivial constant embedding.

Yet, the name spectral comes from the following equivalent formulation. Given the adjacency matrix A ∈ R n×n + of the graph, the embedding is obtained by solving either the eigenvalue problem:

LX = XΛ, with X X = I, (3.1)
or the generalized eigenvalue problem:

LX = DXΛ, with X DX = I, (3.2)
where L = D -A is the Laplacian matrix of the graph, Λ ∈ R k×k is the diagonal matrix of the k smallest (generalized) eigenvalues of L and X ∈ R n×k is the corresponding matrix of (generalized) eigenvectors.

Besides, the spectral embedding can be interpreted as equilibrium states of some physical systems [START_REF] Snell | Random walks and electric networks[END_REF][START_REF] Spielman | Spectral graph theory and its applications[END_REF], Bonald et al., 2018b] as illustrated in Figure 3.1. In this virtual mechanical system, nodes are masses connected by springs corresponding to edges, sliding without friction on an axis rotating at angular speed ω. Let W denote the diagonal matrix of node masses and x the position of the nodes along the axis. The system is in equilibrium if and only if Lx = ω 2 Wx.

If all nodes have the same mass, W ∝ I and x is a solution of (3.1) associated with the eigenvalue λ = ω 2 . If nodes have masses proportional to their degree, W ∝ D and x is a solution of (3.2) for the same eigenvalue. In this chapter, unless otherwise specified, we consider the generalized eigenvalue problem. In view of proposition 1 the solution of (3.2) is given by the spectral decomposition of the normalized Laplacian matrix

L norm = I -D -1/2 AD -1/2 .
This choice will be justified by our experiments. We refer to the eigenvectors of this normalized Laplacian as the vanilla spectral embedding.

In the following sections, we discuss how to improve this vanilla embedding. Specifically, Section 3.2 presents a new interpretation of a standard graph regularization technique, published in Lara and Bonald [2020] and Section 3.3 offers an interpretable extension to directed and bipartite graphs, partially introduced in Bonald and De Lara [2018].

Graphs embeddings are traditionally combined with k-means [MacQueen et al., 1967] for clustering, Ward's algorithm [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] for hierarchical clustering or some classifier such as a nearest neighbor's one for classification. Thus, we defer some related experiments to Chapters 4 and 5, respectively.

Regularization

Vanilla spectral embedding tends to produce poor results on real datasets if applied directly on the graph [START_REF] Amini | Pseudo-likelihood methods for community detection in large sparse networks[END_REF]. One reason is that real graphs are often either disconnected or have sets of dangling nodes loosely connected to the rest of the graph due to noise or outliers in the dataset. In this section, we analyze how a simple graph regularization can help to overcome this issue.

Several regularization techniques have been proposed to improve the quality of the embedding with respect to downstream tasks like clustering. In this section, we explain on a simple block model the impact of the complete graph regularization, whereby a constant is added to all entries of the adjacency matrix. Specifically, we show that the regularization forces the spectral embedding to focus on the largest blocks, making the representation less sensitive to noise or outliers. We illustrate these results on both synthetic and real data, showing how regularization improves standard clustering scores.

In order to improve the quality of the embedding, two main types of regularization have been proposed. The first artificially increases the degree of each node by a constant factor [Chaudhuri et al., 2012, Qin and[START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF], while the second adds a constant to all entries of the original adjacency matrix [START_REF] Amini | Pseudo-likelihood methods for community detection in large sparse networks[END_REF][START_REF] Joseph | Impact of regularization on spectral clustering[END_REF][START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF]. See Figure 3.2 for an illustration. In the practically interesting case where the original adjacency matrix A is sparse, the regularized adjacency matrix is dense but has a so-called sparse + low rank structure, enabling the computation of the spectral embedding on very large graphs as presented in Section 1.4. While [START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF] explains the effects of regularization through graph conductance and [START_REF] Joseph | Impact of regularization on spectral clustering[END_REF] through eigenvector perturbation on the Stochastic Block Model, there is no simple interpretation of the benefits of graph regularization. In this section, we show on a simple block model that the complete graph regularization forces the spectral embedding to separate the blocks in decreasing order of size, making the embedding less sensitive to noise or outliers in the data.

Indeed, [START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF] identified that, without regularization, the cuts corresponding to the first dimensions of the spectral embedding tend to separate small sets of nodes, so-called dangling sets, loosely connected to the rest of the graph. Our work shows more explicitly that regularization forces the spectral embedding to focus on the largest clusters. Moreover, our analysis involves some explicit characterization of the eigenvalues, allowing us to quantify the impact of the regularization parameter.

The rest of this section is organized as follows. Section 3.2.1 presents block models and an important preliminary result about their aggregation. Section 3.2.2 presents the main result, about the regularization of block models, while Section 3.2.5 extends this result to bipartite graphs.

Aggregation of Block Models

Let A ∈ R n×n + be the adjacency matrix of an undirected, weighted graph, that is a symmetric matrix such that A ij > 0 if and only if there is an edge between nodes i and j, with weight A ij . Assume that the n nodes of the graph can be partitioned into K blocks of respective sizes n 1 , . . . , n K so that any two nodes of the same block have the same neighborhood, i.e., the corresponding rows (or columns) of A are the same. Without any loss of generality, we assume that the matrix A has rank K. We refer to such a graph as a block model.

Let Z ∈ R n×K be the associated membership matrix, with Z ij = 1 if index i belongs to block j and 0 otherwise. We denote by W = Z Z ∈ R K×K the diagonal matrix of block sizes. Now define Ā = Z AZ ∈ R K×K . This is the adjacency matrix of the aggregate graph, where each block of the initial graph is replaced by a single node; two nodes in this graph are con-nected by an edge of weight equal to the total weight of edges between the corresponding blocks in the original graph. We denote by D = diag( Ā1 K ) the degree matrix and by L = D -Ā the Laplacian matrix of the aggregate graph.

The following result shows that the solution to the generalized eigenvalue problem (3.2) follows from that of the aggregate graph: Proposition 3. Let x be a solution to the generalized eigenvalue problem: Lx = λDx.

(3.3)

Then either Z x = 0 and λ = 1 or x = Zy where y is a solution to the generalized eigenvalue problem:

Ly = λ Dy. (3.4)
Proof. Consider the following reformulation of the generalized eigenvalue problem (3.3):

Ax = Dx(1 -λ). (3.5)
Since the rank of A is equal to K, there are n -K eigenvectors x associated with the eigenvalue λ = 1, each satisfying Z x = 0. By orthogonality, the other eigenvectors satisfy x = Zy for some vector y ∈ R K . We get:

AZy = DZy(1 -λ), so that Āy = Dy(1 -λ).
Thus y is a solution to the generalized eigenvalue problem (3.4).

Regularization of Block Models

Let A be the adjacency matrix of some undirected graph. We consider a regularized version of the graph where an edge of weight α is added between all pairs of nodes, for some constant α > 0. The corresponding adjacency matrix is given by:

A α = A + αJ,
where J = 1 n 1 n is the all-ones matrix of same dimension as A. We denote by D α = diag(A α 1 n ) the corresponding degree matrix and by L α = D α -A α the Laplacian matrix.

We first consider a simple block model where the graph consists of K disjoint cliques of respective sizes

n 1 > n 2 > • • • > n K nodes, with n K ≥ 1.
In this case, we have A = ZZ , where Z is the membership matrix.

The objective of this section is to demonstrate that, in this setting, the k-th dimension of the spectral embedding isolates the k -1 largest cliques from the rest of the graph, for any k ∈ {2, . . . , K} Lemma 1. Let λ 1 ≤ λ 2 ≤ . . . ≤ λ n be the eigenvalues associated with the generalized eigenvalue problem:

L α x = λD α x. (3.6) We have λ 1 = 0 < λ 2 ≤ . . . ≤ λ K < λ K+1 = . . . = λ n = 1.
Proof. Since the Laplacian matrix L α is positive semi-definite, all eigenvalues are non-negative [START_REF] Chung | Spectral graph theory[END_REF]. We know that the eigenvalue 0 has multiplicity 1 on observing that the regularized graph is connected. Now for any vector x,

x A α x = x Ax + αx Jx = ||Z x|| 2 + α(1 n x) 2 ≥ 0,
so that the matrix A α is positive semi-definite. In view of (3.5), this shows that λ ≤ 1 for any eigenvalue λ. The proof then follows from Proposition 3, on observing that the eigenvalue 1 has multiplicity n -K.

Lemma 2. Let x be a solution to the generalized eigenvalue problem (3.6) with λ ∈ (0, 1). There exists some s ∈ {+1, -1} such that for each node i in block j,

sign(x i ) = s ⇐⇒ n j ≥ α 1 -λ λ n.
Proof. In view of Proposition 3, we have x = Zy where y is a solution to the generalized eigenvalue problem of the aggregate graph, with adjacency matrix:

Āα = Z A α Z = Z (A + αJ)Z. Since A = ZZ and W = Z Z, we have Āα = W 2 + αZ JZ. Using the fact that Z1 K = 1 n , we get J = 1 n 1 n = ZJ K Z with J K = 1 K 1 K the all-ones matrix of dimension K × K, so that: Āα = W(I K + αJ K )W,
where I K is the identity matrix of dimension K × K. We deduce the degree matrix:

Dα = W(W + αnI K ),
and the Laplacian matrix:

Lα = Dα -Āα = αW(nI K -J K W).
The generalized eigenvalue problem associated with the aggregate graph is:

Lα y = λ Dα y.

We get:

α(nI K -J K W)y = λ(W + αnI K )y.
Observing that J K Wy ∝ 1 K , we conclude that:

(αn(1 -λ) -λW)y ∝ 1 K , (3.7) and since W = diag(n 1 , . . . , n K ), ∀j = 1, . . . , K, y j ∝ 1 λn j -α(1 -λ)n . (3.8)
The result then follows from the fact that x = Zy.

Lemma 3. The K smallest eigenvalues satisfy:

0 = λ 1 < µ 1 < λ 2 < µ 2 < • • • < λ K < µ K ,
where for all j = 1, . . . , K,

µ j = αn αn + n j .
Proof. We know from Lemma 1 that the K smallest eigenvalues are in [0, 1). Let x be a solution to the generalized eigenvalue problem (3.6) with λ ∈ (0, 1). We know that x = Zy where y is an eigenvector associated with the same eigenvalue λ for the aggregate graph. Since 1 K is an eigenvector for the eigenvalue 0, we have y Dα 1 K = 0. Using the fact that Dα = W(W + αnI K ), we get:

K ∑ j=1 n j (n j + αn)y j = 0.
We then deduce from (3.7) and (3.8) that λ ∈ {µ 1 , . . . , µ K } and

K ∑ j=1 n j (n j + αn) 1 λ/µ j -1 = 0.
This condition cannot be satisfied if λ < µ 1 or λ > µ K as the terms of the sum would be either all positive or all negative. Now let y be another eigenvector for the aggregate graph, with y Dα y = 0, for the eigenvalue λ ∈ (0, 1). By the same argument, we get:

K ∑ j=1
n j (n j + αn)y j y j = 0, and

K ∑ j=1 n j (n j + αn) 1 λ/µ j -1 1 λ /µ j -1 = 0.
with λ ∈ {µ 1 , . . . , µ K }. This condition cannot be satisfied if λ and λ are in the same interval (µ j , µ j+1 ) for some j as the terms in the sum would be all positive. There are K -1 eigenvalues in (0, 1) for K -1 such intervals, that is one eigenvalue per interval.

The main result of the section is the following, showing that the k -1 largest cliques of the original graph can be recovered from the spectral embedding of the regularized graph in dimension k.

Theorem 1. Let X be the spectral embedding of dimension k, as defined by (3.2), for some k in the set {2, . . . , K}. Then sign(X) gives the k -1 largest blocks of the graph.

Proof. Let x be the j-th column of the matrix X, for some j ∈ {2, . . . , k}. In view of Lemma 3, this is the eigenvector associated with eigenvalue λ j ∈ (µ j-1 , µ j ), so that

α 1 -λ j λ j n ∈ (n j-1 , n j ).
In view of Lemma 2, all entries of x corresponding to blocks of size n 1 , n 2 . . . , n j-1 have the same sign, the other having the opposite sign.

Illustration

We illustrate Theorem 1 with a toy graph consisting of 3 cliques of respective sizes 5, 3 and 2. We compute the embedding in dimension 2 and report the plots in Figure 3.3. As we can see, the embedding without regularization isolates the smallest clique from the two others. On the other hand, in the regularized embedding, the first axis isolates the biggest clique from the others and the second axis isolates the two biggest cliques from the small one.

Extensions

Theorem 1 can be extended in several ways. First, the assumption of distinct block sizes can easily be relaxed. If there are L distinct values of block sizes, say m 1 , . . . , m L blocks of sizes n 1 > . . . > n L , there are L distinct values for the thresholds µ j and thus L distinct values for the eigenvalues λ j in [0, 1), the multiplicity of the j-th smallest eigenvalue being equal to m j . The spectral embedding in dimension k still gives k -1 cliques of the largest sizes.

Secondly, the graph may have edges between blocks. Taking A = ZZ + εJ for instance, for some parameter ε ≥ 0, the results are exactly the same, with α replaced by + α. A key observation is that regularization really matters when ε → 0, in which case the initial graph becomes disconnected and, in the absence of regularization, the spectral embedding may isolate small connected components of the graph. In particular, the regularization makes the spectral embedding much less sensitive to noise, as will be demonstrated in the experiments.

Finally, degree correction can be added by varying the node degrees within blocks. Taking A = θZZ θ, for some arbitrary diagonal matrix θ with positive entries, similar results can be obtained under the regularization A α = A + αθ Jθ. Interestingly, the spectral embedding in dimension k then recovers the k -1 largest blocks in terms of normalized weight, the ratio of the total weight of the block to the number of nodes in the block. For example, a block with 5 nodes connected by edges of weight 10 has a normalized weight of (5 2 × 10)/5 = 50 while a block of 10 nodes connected by edges of weight 4 has a normalized weight of 40. Such case is likely to happen for communication networks. There can be a large group of connected individuals with small communication rate and a much smaller group but with a very high communication rate.

Regularization of Bipartite Graphs

Let B = R n 1 ×n 2 + be the biadjacency matrix of some bipartite graph with respectively n 1 , n 2 nodes in each part, i.e., B ij > 0 if and only if there is an edge between node i in the first part of the graph and node j in the second part of the graph, with weight B ij . This is an undirected graph of n 1 + n 2 nodes with adjacency matrix:

A = 0 B B 0
The spectral embedding of the graph (3.2) can be written in terms of the biadjacency matrix as follows:

BX 2 = D 1 X 1 (I -Λ) B X 1 = D 2 X 2 (I -Λ) (3.9)
where X 1 , X 2 are the embeddings of each part of the graph, with respective dimensions

n 1 × k and n 2 × k, D 1 = diag(B1 n 2 ) and D 2 = diag(B 1 n 1 ).
In particular, the spectral embedding of the graph follows from the generalized SVD of the biadjacency matrix B.

The complete regularization adds edges between all pairs of nodes, breaking the bipartite structure of the graph. Another approach consists in applying the regularization to the biadjacency matrix, i.e., in considering the regularized bipartite graph with biadjacency matrix:

B α = B + αJ,
where J = 1 n 1 1 n 2 is here the all-ones matrix of same dimension as B. The spectral embedding of the regularized graph is that associated with the adjacency matrix:

A α = 0 B α B α 0 (3.10)
As in Section 3.2.2, we consider a block model so that the biadjacency matrix B is block-diagonal with all-ones block matrices on the diagonal. Each part of the graph consists of K groups of nodes of respective sizes n 1,1 > . . . > n 1,K and n 2,1 > . . . > n 2,K , with nodes of block j in the first part connected only to nodes of block j in the second part, for all j = 1, . . . , K.

We consider the generalized eigenvalue problem (3.6) associated with the above matrix A α . In view of (3.9), this is equivalent to the generalized SVD of the regularized biadjacency matrix B α .

We have the following results:

Lemma 4. Let λ 1 ≤ λ 2 ≤ . . . ≤ λ n be the eigenvalues associated with the generalized eigenvalue problem (3.6). We have

λ 1 = 0 < λ 2 ≤ . . . ≤ λ K < λ K+1 = . . . = λ n-2K < . . . < λ n = 2.
Lemma 5. Let x be a solution to the generalized eigenvalue problem (3.6) with λ ∈ (0, 1). There exists s 1 , s 2 ∈ {+1, -1} such that for each node i in block j of part p ∈ {1, 2},

sign(x i ) = s p ⇐⇒ n 1,j n 2,j (n 1,j + αn 1 )(n 2,j + αn 2 ) ≥ 1 -λ.
Lemma 6. The K smallest eigenvalues satisfy:

0 = λ 1 < µ 1 < λ 2 < µ 2 < • • • < λ K < µ K ,
where for all j = 1, . . . , K,

µ j = 1 - n 1,j n 2,j (n 1,j + αn 1 )(n 2,j + αn 2 )
.

Theorem 2. Let X be the spectral embedding of dimension k, as defined by (3.2), for some k in the set {2, . . . , K}. Then sign(X) gives the k -1 largest blocks of each part of the graph.

Like Theorem 1, the assumption of decreasing block sizes can easily be relaxed. Assume that block pairs are indexed in decreasing order of µ j . Then the spectral embedding of dimension k gives the k -1 first block pairs for that order. It is interesting to notice that the order now depends on α:

• When α → 0 + , the block pairs j of highest value ( n 1 n 1,j + n 2 n 2,j ) -1 (equivalently, highest harmonic mean of proportions of nodes in each part of the graph) are isolated first.

• When α → +∞, the block pairs j of highest value n 1,j n 2,j n 1 n 2 (equivalently, the highest geometric mean of proportions of nodes in each part of the graph) are isolated first.

The results also extend to non-block diagonal biadjacency matrices B and degree-corrected models, as for Theorem 1.

The proof of Theorem 2 follows the same workflow as that of Theorem 1. Let Z 1 ∈ R n 1 ×K and Z 2 ∈ R n 2 ×K be the left and right membership matrices for the block matrix

B ∈ R n 1 ×n 2 . The aggregated matrix is B = Z 1 BZ 2 ∈ R K×K . The diagonal matrices of block sizes are W 1 = Z 1 Z 1 and W 2 = Z 2 Z 2 .
We have the equivalent of Proposition 3: Proposition 4. Let x 1 , x 2 be a solution to the generalized singular value problem:

Bx 2 = σD 1 x 1 B x 1 = σD 2 x 2
Then either Z 1 x 1 = Z 2 x 2 = 0 and σ = 0 or x 1 = Z 1 y 1 and x 2 = Z 2 y 2 where y 1 , y 2 is a solution to the generalized singular value problem:

By 2 = σ D1 y 1 , B y 1 = σ D2 y 2 .
Proof. Since the rank of B is equal to K, there are n -K pairs of singular vectors (x 1 , x 2 ) associated with the singular values 0, each satisfying Z 1 x 1 = 0 and Z 2 x 2 = 0. By orthogonality, the other pairs of singular vectors satisfy x 1 = Z 1 y 1 and x 2 = Z 2 y 2 for some vectors y 1 , y 2 ∈ R K . By replacing these in the original generalized singular value problem, we get that (y 1 , y 2 ) is a solution to the generalized singular value problem for the aggregate graph.

In the following, we focus on the block model described in Section 3.2.5, where B = Z 1 Z 2 .

Proof of Lemma 4. The generalized eigenvalue problem (3.6) associated with the regularized matrix A α is equivalent to the generalized SVD of the regularized biadjacency matrix B α :

B α x 2 = σD α,1 x 1 B α x 1 = σD α,2 x 2 , with σ = 1 -λ.
In view of Proposition 4, the singular value σ = 0 has multiplicity n -K, meaning that the eigenvalue λ = 1 has multiplicity n -K. Since the graph is connected, the eigenvalue 0 has multiplicity 1. The proof then follows from the observation that if (x 1 , x 2 ) is a pair of singular vectors for the singular value σ, then the vectors x = (x 1 , ±x 2 ) are eigenvectors for the eigenvalues 1σ, 1 + σ.

Proof of Lemma 5. By Proposition 4, we can focus on the generalized singular value problem for the aggregate graph:

Bα

y 2 = σ Dα,1 y 1 B α y 1 = σ Dα,2 y 2 , Since Bα = W 1 (I K + αJ K )W 2 , and Dα,1 = W 1 (W 2 + αn 1 I), Dα,2 = W 2 (W 1 + αn 2 I),
we have:

W 1 (I K + αJ K )W 2 y 2 = W 1 (W 2 + αn 1 I)y 1 σ, W 2 (I K + αJ K )W 1 y 1 = W 2 (W 1 + αn 2 I)y 2 σ.
Observing that J K W 1 y 1 ∝ 1 K and J K W 2 y 2 ∝ 1 K , we get:

(W 2 + αn 2 I K )y 1 σ -W 2 y 2 ∝ 1 K , (W 1 + αn 1 I K )y 2 σ -W 1 y 1 ∝ 1 K .
As two diagonal matrices commute, we obtain:

(W 1 + αn 1 I K )(W 2 + αn 2 I K )y 1 σ -W 1 W 2 y 1 = η 1 (W 1 + αn 1 I K ) + η 2 W 2 1 K , (W 1 + αn 1 I K )(W 2 + αn 2 I K )y 2 σ -W 1 W 2 y 2 = η 1 W 1 + η 2 (W 2 + αn 2 I K ) 1 K ,
for some constants η 1 , η 2 , and

         y 1,j = η 1 (n 1,j + αn 1 ) + η 2 n 2,j (n 1,j + αn 1 )(n 2,j + αn 2 )σ -n 1,j n 2,j
,

y 2,j = η 1 n 1,j + η 2 (n 2,j + αn 2 ) (n 1,j + αn 1 )(n 2,j + αn 2 )σ -n 1,j n 2,j
.

Letting s 1 = -sign(η 1 (n 1,j + αn 1 ) + η 2 n 2,j ) and s 2 = -sign(η 1 n 1,j + η 2 (n 2,j + αn 2 )), we get:

sign(y 1,j ) = s 1 ⇐⇒ sign(y 2,j ) = s 2 ⇐⇒ n 1,j n 2,j (n 1,j + αn 1 )(n 2,j + αn 2 ) ≥ σ = 1 -λ,
and the result follows from the fact that x 1 = Z 1 y 1 and x 2 = Z 2 y 2 .

Proof of Lemma 6. The proof is the same as that of Lemma 3, where the threshold values follow from Lemma 5:

µ j = 1 - n 1,j n 2,j (n 1,j + αn 1 )(n 2,j + αn 2 ) .
Proof of Theorem 2. Let x be the j-th column of the matrix X, for some j ∈ {2, . . . , k}. In view of Lemma 6, this is the eigenvector associated with eigenvalue λ j ∈ (µ j-1 , µ j ). In view of Lemma 4, all entries of x corresponding to blocks of size n 1,1 , n 1,2 . . . , n 1,j-1 have the same sign, the other having the opposite sign.

Interpretable Graph Embedding

Unlike some ad-hoc methods [START_REF] Funke | Statistical manifold embedding for directed graphs[END_REF][START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF][START_REF] Chen | Directed graph embedding[END_REF], vanilla spectral embedding is not applicable to bipartite or directed graphs. Moreover, the interpretation in terms of random walk is valid only for the full embedding, using all eigenvectors of the Laplacian. This is not feasible for large graphs (e.g., more than 10 000 nodes); in this case, the dimension of the embedding space must be much lower than the number of nodes and there is no clear interpretation of the geometry of the embedding.

In this section, we propose a graph embedding technique based on the generalized singular value decomposition (SVD) [START_REF] Ewerbring | Canonical correlations and generalized svd: applications and new algorithms[END_REF] of the adjacency matrix that can be interpreted in terms of graph structure. Specifically, the distance between nodes in the embedding space corresponds to the distance between their respective neighborhoods in the graph, for some appropriate metric. Moreover, this interpretation is valid in any dimension k, considering the neighborhoods in the best rank-k approximation of the graph rather than in the graph itself. The proposed embedding applies to any type of graphs, including bipartite graphs and directed graphs. It turns out to be equivalent to the spectral embedding of the associate coneighbor graph, after some appropriate scaling depending on the spectrum of the Laplacian matrix.

The rest of the section is organized as follows. We first present some useful matrix norm in Section 3.3.1 and introduce the co-neighbor graph in Section 3.3.2. In Section 3.3.3, we present the proposed embedding technique in the specific case of bipartite graphs, the general cases of directed and undirected graphs being considered in sections 3.3.4 and 3.3.5, respectively. The link with spectral embedding is explained in Section 3.3.6.

Matrix norm

For any positive vector α ∈ R n , we denote by •, • α the following weighted dot product of R n :

∀u, v ∈ R n , u, v α = n ∑ i=1 u i v i α i .
We refer to the corresponding metric as the α metric. We denote by || • || α the corresponding norm:

||u|| 2 α = u diag(α) -1 u = ||diag(α) -1 2 u|| 2 ,
where || • || is the standard Euclidian norm of R n .

For any positive vectors α ∈ R n , β ∈ R m , we define the weighted spectral norm of any matrix M of size n 1 × n 2 for the weight vectors α, β as:

||M|| α,β = sup u∈R n 2 :||u|| β =1 ||Mu|| α . Observe that ||M|| α,β = ||diag(α) -1 2 Mdiag(β) -1 2 || σ ,

Co-neighbor graph

Consider any non-negative matrix A of size n 1 × n 2 . This may be viewed as the biadjacency matrix of a weighted bipartite graph G with n 1 nodes in one part and n 2 nodes in the other.

Equivalently, each row of A may be interpreted as a data sample and each column as a feature. Thus there are n 1 data samples and n 2 features, each row of A corresponding to the feature vector of a data sample. We denote by d = A1 and f = A 1 the weight vectors of data samples and features, respectively, and by D = diag(d) and F = diag( f ) the corresponding diagonal matrices. If the matrix A is binary, the vectors d and f correspond to the degrees the nodes of each part of the bipartite graph G. We assume that the vectors d and f are positive (i.e., there is no all-zero row or column).

The proposed embedding is related to the spectral embedding of the co-neighbor graph. The co-neighbor graph is a weighted graph of n 1 nodes (the data samples) with adjacency matrix:

S = AF -1 A .
The weight between nodes i and j can be interpreted as the similarity between the corresponding data samples i and j:

S ij = n 2 ∑ l=1 A il A jl f l = d i d j w p i , p j β .
(3.11)

Note that the matrix S is symmetric positive semi-definite and has the same rank than A.

Besides, the degree of a node in the co-neighbor graph is the same as in the original graph S1 = A1 = d.

Graph embedding

The normalized weight vectors define probability distributions over the data samples and the features, respectively:

α = d w , β = f w ,
where w = d1 = f 1 = 1 A1 is the total weight of the nodes. We have:

β = α P, (3.12)
where P = D -1 A. Observe that each row of P gives the feature distribution of a data sample, as a probability distribution over {1, . . . , n 2 }.

Generalized SVD. Consider the generalized SVD of the matrix A associated with the diagonal matrices D and F:

AV = DUΣ, A U = FVΣ, with U DU = I, V FV = I. (3.13)
The matrices U = (u 1 , . . . , u r ) and V = (v 1 , . . . , v r ) are formed by the left and right generalized singular vectors of A associated with the generalized singular values

σ 1 ≥ σ 2 ≥ . . . ≥ σ r > 0,
where r is the rank of A, and Σ = diag(σ 1 , . . . , σ r ) is the diagonal matrix of positive generalized singular values. We have: k) be the respective restrictions of the matrices U, V to their first k columns, and let Σ (k) = diag(σ 1 , . . . , σ k ). Define the matrix:

A = DUΣV F. Now for any k ≤ r, let U (k) , V (
A (k) = DU (k) Σ (k) V (k) F.
This is an optimal rank-k approximation of the matrix A for the norm || • || α,β defined in section 3.3.1:

Proposition 5. We have:

min M:rank(M)=k ||A -M|| α,β = ||A -A (k) || α,β = σ k+1 w if k < r, 0 if k = r.
Proof. As the regular SVD provides an optimal low-rank approximation of a matrix, we have:

min M:rank(M)=k ||D -1 2 AF -1 2 -M|| = ||D -1 2 AF -1 2 -Ũ(k) Σ Ṽ(k) || = σ k+1 if k < r, 0 if k = r,
where Ũ(k) , Ṽ(k) are the restrictions of the matrices Ũ, Ṽ to their first k columns. Observing that M has rank k if and only if D

1 2 MF 1 2 has rank k, we get: min M:rank(M)=k ||D -1 2 (A -M)F -1 2 || = ||D -1 2 (A -A (k) )F -1 2 || = σ k+1 if k < r, 0 if k = r.
The proof follows from the definition of the norm || • || α,β .

This is illustrated by Figure 3.4 for some toy graph, where the width of each edge in the lowrank approximation graph corresponds to its weight (the original graph is unweighted). It turns out that the weight vectors d and f are preserved by the low-rank approximation, provided k is larger than the number of connected components of the graph: Proposition 6. Let K be the number of connected components of the graph G. We have:

1. σ 1 = . . . = σ K = 1 > σ K+1 , 2. if K = 1, then u 1 = 1/ √ w and v 1 = 1/ √ w, 3. α u k = 0 and β v k = 0 for all k > K, 3. Let u 1 = 1/ √ w and v 1 = 1/ √ w.
Since u 1 , v 1 are generalized singular vectors of A for the singular value 1, there exists some generalized singular decomposition U , V of A such that u 1 , v 1 are the respective first columns of U , V . Now the matrices D 1 2 U , F 1 2 V form a regular singular decomposition of the normalized adjacency matrix D -1 2 AF -1 2 . We deduce that U (K) = U (K) G and V (K) = V (K) H for some orthogonal matrices G, H of size K × K. Now for any k > K, we have u k DU (K) = 0 and v k FV (K) = 0 so that u k DU (K) = 0 and v k FV (K) = 0. In particular, u k D1 = 0 and v k F1 = 0.

4. Writing A (k) = ∑ k l=1 σ l Du l v l F and using the fact that v l F1 = 0 for all l > K, we get Full embedding. We first consider the full embedding, in dimension r (the rank of A). Define:

A (k) 1 = A1 = d for all k > K. Similarly, A (k) 1 = A 1 = f .
X = √ w(I -1α )UΣ.
Each column of X has its center of mass at the origin, for the weight vector α:

α X = 0.
The embedding is then given by the rows of X. We denote the corresponding vectors by x 1 , . . . , x n . Thus each data sample is embedded as a vector of dimension r.

Interestingly, the features can be embedded in the same space. Let:

Y = √ w(I -1β )V.
Each column of Y has its center of mass at the origin, for the weight vector β:

β Y = 0.
The co-embedding of features is given by the rows of Y. We denote the corresponding vectors by y 1 , . . . , y n 2 .

Proposition 7. We have:

X = PY.
Proof. Using (3.12), we get:

PY = √ wP(I -1β )V, = √ w(I -1α )PV, = √ w(I -1α )UΣ = X.
In view of Proposition 7, we have:

x i = n 2 ∑ j=1 p ij y j .
Thus each data sample i is embedded by some vector x i which is the weighted average of its features. This is illustrated by Figure 3.5, for a toy bipartite graph connecting movies to actors. See how each movie is at the barycenter of its actors. In particular, two data samples with the same feature vectors (up to some multiplicative constant) have the same embedding. Geometry of the embedding. As noted earlier, the rows of P can be interpreted as feature distributions. Specifically, the feature distribution of data sample i is given by:

p i = P e i .
In view of (3.12), the weighted average feature distribution is given by:

β = P α = n 1 ∑ i=1 α i p i .
The main result of the section is the following.

Theorem 3. We have for all i, j = 1, . . . , n 1 :

x i , x j = p i -β, p j -β β = p i , p j β -1.
Proof. Using the fact that:

PF -1 P = D -1 SD -1 = UΣ 2 U ,
we get:

XX = w(I -1α )UΣ 2 U (I -1α ) , = w(I -1α )PF -1 P (I -1α ) , = (P -1β )diag(β) -1 (P -1β ) ,
for which the proof follows.

Theorem 3, shows that the geometry of the embedding reflects the graph structure. In particular:

Vector norms. We have

||x i || 2 = ||p i -β|| 2 β = ||p i || 2 β -1.
The norm of data sample i in the embedding space is equal to the distance between its feature distribution p i and the average feature distribution β for the β metric. Data samples that have feature distributions close to the average are close to the origin in the embedding space; data samples that have unusual feature distributions (compared to the average) are far from the origin. Note that ||x i || 2 ≤ 1/ min j β j -1.

Distances. We have ||x

i -x j || 2 = ||p i -p j || 2 β .
The distance between two data samples in the embedding space is equal to the distance between their respective feature distributions for the β metric.

Angles. The cosine similarity between data samples i and j is the embedding space is given by:

s(i, j) = x i , x j ||x i ||||x j || = p i , p j β -1 ||p i -β|| β ||p j -β|| β .
In particular, we have s(i, j) = 1 if and only if p i = p j (same feature distributions) and s(i, j) < 0 whenever p i , p j β = 0 (disjoint feature distributions).

Illustration Consider the toy bipartite graph displayed in Figure 3.6. The movies, indexed by chronological order, are the feature nodes so that β = (3/8, 3/8, 2/8). The sample nodes are the characters. For example, p Greedo = (1, 0, 0) and p Vador = (1/3, 1/3, 1/3). So, in the full embedding, we have

x Greedo , x Vador = (1 × 1 3 )/( 3 8 ) -1, = - 1 9 .
The same derivation with Jabba, p Jabba = (1/2, 0, 1/2), leads to Embedding in low dimension. Now consider the embedding in dimension k, for any k ∈ {K + 1, . . . , r}. Let:

x Jabba , x Vador = 1 2 × 1 3 . 8 3 + 8 2 -1, = 1
X (k) = √ w(I -1α )U (k) Σ (k) .
Again, each column of X (k) has its center of mass at the origin:

α X (k) = 0.
We embed the data samples as vectors of dimension k, corresponding to the rows of X (k) . We denote the corresponding vectors by x

(k) 1 , . . . , x (k)
n . Similarly, we co-embed the features as vectors of dimension k, corresponding to the rows of

Y (k) = √ w(I -1β )V (k) .
We denote the corresponding vectors by y

(k) 1 , . . . , y (k)
m . We have:

X (k) = P (k) Y (k) with P (k) = D -1 A (k)
. The proof is similar to that of Proposition 7, using the following extension of (3.12):

β = α P (k) . (3.14)
Theorem 4. We have for all i, j = 1, . . . , n:

x

(k) i , x (k) j = p (k) i , p (k) j β -1.
Proof. The proof is similar, using (3.14) and the fact that:

P (k) F -1 P (k) = D -1 S (k) D -1 = U (k) Σ (k) 2 U (k) .
In view of Theorem 4, the geometry of the embedding in dimension k can be interpreted as the full embedding but for the best rank-k approximation of the graph rather than for the graph itself.

Co-neighbor graph. As introduced in Section 3.3.2, this embedding is related to the spectral embedding of the co-neighbor graph S.

In view of (3.13), we have:

SU = DUΣ 2 , U DU = I. (3.15)
This is the generalized eigenvalue problem of S associated with the diagonal matrix of weights D. Observe that the spectral embedding of the co-neighbor graph is based on the generalized spectral decomposition of the associate Laplacian matrix L = D -S. In view of (3.15), its eigenvalues (in non-decreasing order) satisfy:

λ 1 = . . . = λ K = 0 < λ K+1 = 1 -σ 2 K+1 ≤ . . . ≤ λ r = 1 -σ 2 r < λ r+1 = . . . = λ n = 1.
Thus this embedding corresponds to the spectral embedding of the normalized co-neighbor graph after scaling by 1λ 1 , . . . , 1λ n and shift so that the center of mass is at the origin. Observe that only eigenvectors associated with eigenvalues less than 1 are kept (the other vanish), so that the dimension of the embedding is at most r.

For most real graphs, the biadjacency matrix A is sparse but the similarity matrix S is too dense to be stored. Consider for instance the node Beyoncé in the Twitter graph. This artist has 15.5 million followers, hence the corresponding node has an in-degree of 15.5 million. In the coneighbor graphs, these followers form a clique of the same size, generating potentially 240 000 billion of nonzero entries in S. However, using techniques introduced in Section 1.4, there is no need to actually store the similarity matrix to compute the embedding.

Case of directed graphs

Consider a directed graph G of n nodes with adjacency matrix A. Observe that the matrix A is square but not symmetric in general. Let A + be the reduced adjacency matrix, obtained after removing all-zero rows and columns. The graph G can then be considered as a bipartite graph with biadjacency matrix A + . There are n + nodes in the first part (nodes of G having positive out-degrees) and n -nodes in the second part (nodes of G having positive in-degrees). We simply apply the embedding described in section 3.3.3 to this bipartite graph.

The geometry of the embedding can be interpreted using Theorems 3 and 4, with p i and p (k) i corresponding to the distributions of successors of node i in the graph G and in its best rank-k approximation, respectively. The embedding also corresponds to the spectral embedding of the co-neighbor graph, after appropriate scaling and shifting. The similarity S ij between two nodes i, j of G depends on their out-degrees and on the weighted dot-product of their distributions of successors p i , p j , according to (3.11).

Case of undirected graphs

Consider any undirected graph G of n nodes, with adjacency matrix A. Now A is both square and symmetric. Without loss of generality, we assume that all nodes have positive degrees.

Similarly, the graph G may be seen as a bipartite graph with biadjacency matrix A. There are n nodes in each part. We then directly apply the embedding described in section 3.3.3 to this bipartite graph.

Again, the geometry of the embedding can be interpreted using Theorems 3 and 4, with p i and p (k) i corresponding to the distributions of successors of node i in the graph G and in its best rank-k approximation, respectively. The embedding also corresponds to the spectral embedding of the co-neighbor graph, after appropriate scaling and shifting. This is not the spectral embedding of the graph G. In particular, two nodes with the same neighborhoods have the same representation in the embedding space, which is not the case for the spectral embedding of the graph G.

Link with spectral embedding

We have seen that the proposed embedding is related to the spectral embedding of the cocitation graph, after appropriate scaling and shifting. This is not a practical method, however, as the similarity matrix S is most often too dense to be stored, meaning that the co-neighbor graph cannot be built. In this section, we show that the proposed embedding is also related to the spectral embedding of some extended graph, which can be built in practice.

Let A be any non-negative matrix of size n 1 × n 2 , as considered in section 3.3.3. This is the biadjacency matrix of some bipartite graph G. Now consider the extended graph G, with adjacency matrix and degree matrix:

à = 0 A A 0 , D = D 0 0 F .
Note that the graph G has n 1 + n 2 nodes. In view of (3.13), we have:

à U V = D U V Σ, with U V D U V = I.
So the generalized singular values σ 1 , . . . , σ r of A are the generalized positive eigenvalues of Ã. The spectral embedding of the graph G is based on the spectral decomposition of the Laplacian matrix L = D -Ã. We get:

L U V = D U V Λ.
with Λ = I -Σ. So our embedding corresponds to the spectral embedding of the extended graph restricted to positive generalized eigenvalues, after scaling by 1λ 1 , . . . , 1λ r and shifting so that the center of mass is at the origin.

In the particular case where A is the adjacency matrix of some undirected graph G, the embedding differs from the spectral embedding of the graph G, which is based on the spectral decomposition of the Laplacian L = D -A (instead of L).

About embedding metrics

Assessing the quality of an embedding without a supervised downstream task is not trivial. A good embedding should fit the data in the sense that nodes which are considered close according to the topology of the graph should also be close in the embedding space. Such fit functions F include:

• Spectral: F = -∑ i,j A ij ||x i -x j || 2 .
• Node2Vec: F = ∑ (i,j)∈E

x i x j .

However, such functions have either trivial or ill-defined optimums. There are two main strategies to overcome this issue. Spectral methods rely on orthonormality constraints while iterative methods such as [START_REF] Grover | node2vec: Scalable Feature Learning for Networks[END_REF] or Tang et al. [2015b] rely on a penalization function D. For example, in node2vec:

D = -∑ i log ∑ j exp(x i x j ) .
This idea is similar to the one used for graph visualization where nodes apply repulsive forces on each other [START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF]. However, in both cases, the cost of evaluating D is quadratic in the number of nodes. Node2vec relies on negative sampling in order to maintain tractable computations but this comes at the cost of an extra source of randomness and supplementary hyper-parameters. The final quality function is simply Q = F -D.

Here, we propose a different one, as it is inspired by the modularity [Newman, 2006a] for graph clustering, we refer to it as cosine modularity:

Q = ∑ ij A ij w -γ d i f j w 2 1 + cos(x i , x j ) 2 , (3.16)
for some resolution parameter γ which is set to 1 by default. See Example 3.4.1 for an illustration.

Example. 3.4.1: Cosine modularity

[1]: from sknetwork.embedding import cosine_modularity from sknetwork.data import karate_club karate = karate_club(metadata=True) A = karate.adjacency X = karate.position round(cosine_modularity(A, X), 2)

[1]: 0.35

This quality function has a few interesting properties:

1. It maintains the standard Q = F -D structure.

2. As for the standard modularity, Q ∈ [-γ, 1].

3. Its asymptotic complexity is only O(m).

Proof. Let us denote s ij = 1 + cos(x i , x j ) 2 .

1. In equation (3.16), we identify:

F = ∑ ij A ij w s ij , D = ∑ ij γ d i f j w 2 s ij . 2. As s ij ∈ [0, 1], ∀(i, j), Q ≤ ∑ ij A ij w = 1, Q ≥ -∑ ij γ d i f j w 2 = -γ.
3. F is a sum over m elements. For Q, let us denote by x the embedding projected onto the unit-sphere: xi = x i /||x i || 2 . Then:

D = γ 2 1 + ∑ ij d i f j w 2 x i xj , = γ 2   1 + ∑ i d i xi w ∑ j f j xj w   .
The proof of the last point highlights one major limitation of this quality function. It is sufficient for x to satisfy x, 1 α = 0 or x, 1 β = 0 in order to yield the minimal penalization.

Another approach to quantify the quality of an embedding is self-supervision through linkprediction. The idea is to hide a certain fraction of the edges, embed the nodes, predict the most likely missing edges based on the embedding and compare this prediction to the actually missing edges. The likelihood of an edge is usually either given by x i , x j or cos(x i , x j ).

Still, evaluating the performance of a link prediction algorithm is also a challenging task. For example, the brute force computation of the edges sorted by likelihood has a complexity of O(n 2 log(n)) which is prohibitive for large graphs. The selection of the hidden edges and the choice of the metric are also non-trivial. For more information on this topic, we refer to the work of [START_REF] Yang | Evaluating link prediction methods[END_REF] as this lies beyond the scope of the present thesis.

Experiments

This section describes part of the experiments related to the embeddings. Sections 3.5.1 and 3.5.2 respectively introduce some metrics and the datasets. Section 3.5.3 justifies the choice of the normalized Laplacian. Section 3.5.4 outlines the influence of graph regularization for spectral clustering introduced in Section 3.2 while Section 3.5.5 studies the influence of the shift in the GSVD embedding. Finally, Section 3.5.6 highlights the importance of normalization for the spectral embedding. We defer general benchmarks to Chapters 4 and 5 respectively in order to focus on the different spectral embeddings in the present section.

Metrics

We consider a large set of metrics from the clustering literature. All metrics are upper-bounded by 1 and the higher the score the better. For supervised metrics, we denote by c the predicted clustering and by y the ground truth labels.

Homogeneity (H), Completeness (C) and V-measure score (V) [START_REF] Rosenberg | V-measure: A conditional entropy-based external cluster evaluation measure[END_REF]. Supervised metrics. A cluster is homogeneous if all its data points are members of a single class in the ground truth. A clustering is complete if all the members of a class in the ground truth belong to the same cluster in the prediction. Then, the V-measure score is the harmonic mean of homogeneity and completeness.

Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF]. Supervised metric. This is the corrected for chance version of the Rand Index (equation (3.17)) which is itself an accuracy on pairs of samples. A pair of samples is correctly classified if δ c i c j = δ y i y j . This means that either belong to the same cluster in the prediction and in the ground truth or they belong to different clusters in both partitions.

RI(c, y) = 1

n 2 ∑ i,j
1(δ c i c j = δ y i y j ).

(3.17)

Then, the ARI is

ARI = RI -E(RI) max(RI) -E(RI) ,
where the expected value is calculated with respect to a random cluster assignment.

Adjusted Mutual Information (AMI) [START_REF] Vinh | Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance[END_REF] Supervised metric. Adjusted for chance version of the mutual information.

Fowlkes-Mallows Index (FMI) [START_REF] Fowlkes | A method for comparing two hierarchical clusterings[END_REF]]. Supervised metric. Geometric mean between precision and recall on the edge classification task, as described for the ARI. [Newman, 2006a]. Unsupervised metric as in equation (4.3). Fraction of edges within clusters compared to that is some null model where edges are shuffled at random.

Modularity (M)

Normalized Standard Deviation (NSD)

Unsupervised metric. 1 minus normalized standard deviation in cluster size. This metric provides insight on the balance of the clustering.

Datasets

This section describes the datasets used in our experiments. Table 3.1 presents the main features of the graphs.

Stochastic Block-Model (SBM)

We generate 100 instances of the same stochastic block model [START_REF] Holland | Stochastic blockmodels: First steps[END_REF]. There are 100 blocks of size 20, with intra-block edge probability set to 0.5 for the first 50 blocks and 0.05 for the other blocks. The inter-block edge probability is set to 0.001. Other sets of parameters can be tested using the code available online. The ground-truth cluster of each node corresponds to its block.

20newsgroup (NG)

This dataset consists of around 18 000 newsgroups posts on 20 topics. This defines a weighted bipartite graph between documents and words. The label of each document corresponds to the topic.

Wikipedia for Schools (WS) [START_REF] Haruechaiyasak | Article recommendation based on a topic model for wikipedia selection for schools[END_REF]. This is the graph of hyperlinks between a subset of Wikipedia pages. The label of each page is its category (e.g., countries, mammals, physics).

Cora (CO) and CiteSeer (CS)

Citation networks between scientific publications. These are standard datasets for node classification [START_REF] Fey | Splinecnn: Fast geometric deep learning with continuous b-spline kernels[END_REF][START_REF] Huang | Adaptive sampling towards fast graph representation learning[END_REF][START_REF] Wijesinghe | Dfnets: Spectral cnns for graphs with feedback-looped filters[END_REF]. Labels correspond to social communities.

Wikipedia Vitals (WV) Graph of hyperlinks between a selection of

Influence of Laplacian normalization

In this experiment, we justify our choice to favor the normalized Laplacian with respect to the regular one. In order to do so, we compare their performance on a semi-supervised task. We extract 16 eigenvectors and apply a 1-nearest-neighbor classifier with 1% of the nodes, selected uniformly at random, as labeled seeds. We do not use regularization nor unit-sphere normalization. For this experiment, all graphs are treated as undirected. The process is repeated 10 times for each graph. We report the resulting macro F1-scores in Table 3.2. Recall that, for a multilabel classification, the macro F1-score is the average of F1-scores obtained for the prediction of each class in a one-against-all setting. As we can see, the normalized Laplacian outperforms the regular one by a significant margin.

Laplacian CO CS WS WV regular 0.14 ± 0.01 0.17 ± 0.01 0.06 ± 0.00 0.10 ± 00 normalized 0.44 ± 0.05 0.28 ± 0.03 0.29 ± 0.06 0.50 ± 0.02 Table 3.2: Macro F1-scores (mean ± standard deviation).

Influence of regularization

We now illustrate the impact of regularization on the quality of spectral embedding. In all experiments, we skip the first dimension of the spectral embedding as it is not informative (the corresponding eigenvector is the all-ones vector, up to some multiplicative constant).

Experimental setup

All graphs are embedded in dimension 20, with different regularization parameters. To compare the impact of this parameter across different datasets, we use a relative regularization parameter (w/n 2 )α, where w = 1 n A1 n is the total weight of the graph, as illustrated in Example 3.5.1.

Example. 3.5.1: Regularized spectral embedding

[1]: from sknetwork.embedding import Spectral spectral = Spectral(n_components=20, normalized_laplacian=True, regularization=0.01, relative_regularization=True)

We use the k-means algorithm to cluster the nodes in the embedding space. The parameter k is set to the ground-truth number of clusters. We use the Scikit-learn [Pedregosa et al., 2011a] implementation of the metrics, when available.

Results

We report the results in Table 3.3 for α ∈ {0, 0.1, 1, 10}. We see that the regularization generally improves performance, the optimal value of α depending on both the dataset and the score function. As suggested by Lemma 3, the optimal value of the regularization parameter should depend on the distribution of cluster sizes, on which we do not have any prior knowledge.

To test the impact of noise on the spectral embedding, we add isolated nodes with self loop to the graph and compare the clustering performance with and without regularization. The number of isolated nodes is given as a fraction of the initial number of nodes in the graph. Scores are computed only on the initial nodes. The results are reported in Table 3.4 for the Wikipedia for Schools dataset. We observe that, in the absence of regularization, the scores drop even with only 1% noise. The computed clustering is a trivial partition with all initial nodes in the same cluster. This means that the 20 first dimensions of the spectral embedding focus on the isolated nodes. On the other hand, the scores remain approximately constant in the regularized case, which suggests that regularization makes the embedding robust to this type of noise.

SBM

Influence of embedding shift

The embedding described in Section 3.3.3 is shifted so that its center of mass is at the origin for the weight vector α. In this experiment, we intend to assess whether this shift leads to different results in practice. Note that, this shift does not change the nearest neighbor ordering with respect to the Euclidean distance, however, it might change it with respect to the cosine similarity.

The setup is similar to the one of Section 3.5.3, except that embeddings are normalized onto the unit-sphere so that the nearest neighbor is defined with respect to the cosine similarity instead of the Euclidean distance. In addition to the macro F1-score, we report the mean norm of the nodes before the shift and the norm of the shift vector ∆ = (UΣ) α in Table 3.5. As the macro F1-score is not modified by the shift for any of the datasets, we only report it once. This can be surprising as the norm of the shift is not negligible with respect to the average norm of the node vectors. Thus, understanding the influence of this shift requires further investigation. We omit it in further experiments for simplicity.

Cosine similarity versus L2 norm

In this experiment, we highlight the importance of embedding normalization for spectral clustering. Recall that projecting the node vectors onto the unit sphere is equivalent to use the cosine similarity in the embedding space rather than the Euclidean distance.

We embed the nodes of the Cora dataset in dimension 16 and then apply k-means for k = 5. If the embedding is not normalized, this yields a V-score of 0.01 while the score is 0.21 for 0.35 0.12 1.00 1.00 macro-F1 (µ ± σ) 0.43 ± 0.06 0.26 ± 0.02 0.28 ± 0.01 0.52 ± 0.02 Table 3.5: Some measurements. µ: mean. σ: standard deviation. the normalized embedding. This is a direct consequence of the cluster sizes. Indeed, without normalization, the NSD of the clustering is 0.04, which means that most of the nodes belong to a unique giant cluster while the four other clusters have very few nodes. With the normalization, the NSD is 0.80, which means that the clusters sizes are very balanced.

α = 0 noise H C V
The aggregate graphs are displayed in Figure 3.7. Our conclusion is that cosine similarity should be preferred for spectral clustering. soft clustering. This section studies some of these variations. Especially, we show in Section 4.1.2 that a generalized version of the Louvain algorithm can be used to maximize different types of modularity function that we introduce in Section 4.1.1.

Modularity functions

In this section, we define the different modularity functions used for graph clustering and present some of their properties. Note that we do not discuss soft-modularity functions here i.e. modularity functions that take a soft partition as input. However, we refer the interested reader to the work of [START_REF] Hollocou | Modularity-based Sparse Soft Graph Clustering[END_REF].

Definitions

Consider a graph G = (V, E), and two probability distributions p and q over V × V. The distribution q is sometimes referred to as the null model of the graph. The general form of the modularity function, given a resolution parameter γ is

M = ∑ i,j p ij -γq ij δ c i c j . (4.1)
Or, aggregating the terms cluster by cluster:

M = ∑ k (P k -γQ k ) , (4.2)
where P k = ∑ (i,j)∈C k ×C k p ij is the probability of sampling a pair of nodes in cluster k according to distribution p. Similar definition for Q k .

The standard choice of p and q for undirected graphs are p ij ∝ A ij and q ij ∝ d i d j , which lead to the modularity as defined by Newman [2006a]:

M = ∑ i,j A ij w -γ d i d j w 2 δ c i c j . (4.3)
Another choice is p ij ∝ A ij and q ij ∝ 1, which is referred to as the constant Potts model [START_REF] Traag | Narrow scope for resolution-limit-free community detection[END_REF]:

M = ∑ i,j A ij w -γ 1 n 2 δ c i c j . (4.4)
If the graph is directed, two different approaches have been proposed. On the one hand, [START_REF] Dugué | Directed Louvain: maximizing modularity in directed networks[END_REF] suggests to use p ij ∝ A ij and q ij ∝ d i f j , which leads to:

M = ∑ i,j A ij w -γ d i f j w 2 δ c i c j . (4.5)
On the other hand, [START_REF] Lambiotte | Laplacian dynamics and multiscale modular structure in networks[END_REF] and [START_REF] Kim | Finding communities in directed networks[END_REF] suggest to use p ij = A ij π j /d j and q ij = π i π j , where π denotes the Pagerank probability distribution without damping factor. This leads to:

M = ∑ i,j A ij π j d j -γπ i π j δ c i c j .
The authors prove that this modularity is the same as equation ( 4.3) applied to the modified adjacency Ãij = 1 2

A ij π j d j + A ji π i d i .
In the case of bipartite graphs, [START_REF] Barber | Modularity and community detection in bipartite networks[END_REF] proposed to use p ij = p ji ∝ B ij and q ij ∝ d i f j if (i, j) ∈ V 1 × V 2 and q ij = 0 otherwise. Which leads to the bimodularity function:

M = ∑ i∈V 1 ,j∈V 2 B ij w -γ d i f j w 2 δ c i c j .
Proposition 8. Let B be the biadjacency matrix of a bipartite graph and let

A = 0 B 0 0 ∈ R (n 1 +n 2 )×(n 1 +n 2 )
be the adjacency matrix of this graph, seen as a directed one. Given a clustering c, the bimodularity defined in equation ( 4.1.1) coincides with the directed-modularity as defined in equation (4.5).

Proof. First, let us note that A and B have the same total weight w. Next, note that d A = (0, d B ) and f A = (0, f B ). Let Z = (Z 1 , Z 2 ) be the membership matrix. Then, starting from Equation 4.5:

M = ∑ i,j A ij w -γ d A i f A j w 2 δ c i c j , = Tr Z A w - γ w 2 d A ( f A ) Z , = Tr Z 1 B w - γ w 2 d B ( f B ) Z 2 .

Properties

A key property of the modularity functions, at the core of Louvain algorithm, is following the aggregation property described in Proposition 9.

Proposition 9. The modularity of a partition of the nodes of a graph G = (V, E) is equal to the modularity of the trivial partition where each node is assigned to its own cluster for the graph G = ( Ṽ, Ẽ) such that Ṽ and Ẽ are the set of nodes and edges obtained by merging all nodes in the same cluster in the original partition.

Proof. Let à = Z AZ be the adjacency matrix of the aggregated graph. We have

d = Ã1 = Z AZ1 = Z d, f = Ã 1 = Z A Z1 = Z f , w = 1 Ã1 = 1 Z AZ1 = 1 A1 = w.
Hence, pkk = (Z pZ) kk and qkk = (Z qZ) kk hold for all the probability distributions listed in Subsection 4.1.1:

M = ∑ i,j p ij -γq ij δ c i c j , = ∑ i,j p ij -γq ij Z i Z j , = Tr Z (p -γq) Z , = Tr I K Z pZ -γZ qZ I K , = ∑ k,k ( pkk -γ qkk ) δ kk .
Another property of the generalized modularity function defined in equation ( 4.1) is its connection to block models for graphs. Consider a model in which, given a clustering and ω in > ω out ≥ 0, the weight of an edge between i and

j is Âij = ω in , if δ c i c j = 1, ω out , otherwise.
We denote by µ(δ) the edge sampling probability for this model. In particular,

log(µ ij ) = log ω in ω out δ c i c j + log(ω out ), ∀(i, j),
up to an additive constant.

We want to find a clustering such that the resulting distribution is both close to the edge sampling distribution p and far from the null model q. To do so, we seek to minimize over δ the following loss:

L = KL (p, µ(δ)) -γKL (q, µ(δ)) , (4.6)
where KL denotes the Kullback-Leibler divergence. Theorem 5 is the main result of this chapter.

Theorem 5. The set of optimal clusterings for L is equal to the set of optimal clustering for the generalized modularity metric:

arg min δ L = arg max δ M. (4.7)

Generalized Louvain algorithm

This section describes the Generalized Louvain algorithm and discusses some implementation issues. Even though the Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] was introduced to maximize the modularity for undirected networks as in equation ( 4.3), the general idea can be applied to any type of objective function that satisfies two properties:

• Easy computation of variation in score resulting from the modification of the cluster assignment of a single node.

• Aggregation property as described in Proposition 9.

The pipeline for the Louvain algorithm consists in two nested loops which we call "Macro-Louvain" and "Micro-Louvain", respectively described in Algorithms 4.1 and 4.2.

Local updates

The efficiency of the Generalized Louvain Algorithm relies on the simple closed form solution of arg max

c ∆M(c i ← c)
in Algorithm 4.2. Note that it is not necessary to test all possible values of c but only the different labels among the neighbors of i.

From equation (4.2), we have:

∆M(c i ← c j ) = ∑ k (∆P k -γ∆Q k ) , = (∆P c i -γ∆Q c i ) + ∆P c j -γ∆Q c j ,
as i and j are the only two modified clusters. Actually, we only need to derive the variations for M corresponding to equation (4.5) to cover all cases listed in subsection 4.1.1.

For the case of equation (4.5), we have, on the one hand,

∆P c i = 1 w 2A ii -∑ c l =c i (A il + A li ) , ∆P c j = 1 w ∑ c l =c j (A il + A li ) .
For the null model, on the other hand, it is convenient to rewrite

Q k = 1 w 2 D k F k , where D k = ∑ c i =k d i is the total out-degree of cluster k and F k = ∑ c i =k f i its total in-degree. Then, ∆Q c i = 1 w 2 [ f i (d i -D c i ) + d i ( f i -F c i )] , ∆Q c j = 1 w 2 f i D c j + d i F c j .
The case of equation (4.3) is simply the particular case where A ij = A ji and d i = f i . Besides, it is sufficient to replace d and f by 1 to recover (4.4). Note that the values D k and F k can easily be stored and updated during the computations. Z ← Z Ẑ.

5:

A ← Z AZ. 6: until Convergence. 7: return c(Z).

Algorithm 4.2 Micro-Louvain

Require: Adjacency A.

1: Initialize c = [1, 2, . . . , n]. 2: repeat 3: for i = 1 ≤ • • • ≤ n do 4: c i = arg max c j
∆M(c i ← c j ).

5:

end for 6: until Convergence. 7: return Z(c).

Convergence

The modularity is non-decreasing under the local updates performed in Micro-Louvain and the aggregation property 9 ensures that it is not decreasing under each update performed by Macro-Louvain. As the modularity is upper-bounded, the Generalized Louvain Algorithm converges towards a local maximum of the modularity. In practice, the stopping criterion for Macro and Micro-Louvain can be a maximum number of passes, a tolerance parameter for the increase in modularity, or both.

Parallel computing for Bilouvain

In the case of bimodularity optimization, the local update for a node in V 1 only depends on the cluster assignment of the nodes in V 2 and vice versa. It is thus technically possible to perform all local updates for nodes of V 1 in parallel (the same goes for V 2 ). However, as Louvain is a greedy method, it is highly dependent on the order in which the cluster assignments are updated. Such an arbitrary order has no guarantee to result in better partitions than a random order.

Recall that, as for k-means initialization, a common practice to maximize modularity is to run several instances of Louvain with different node orderings and select the best partition afterwards. We illustrate this in Example 4.1.1. Two different orderings of the nodes yield distinct partitions with similar modularity. For these algorithms, we test different combinations of embedding dimension d and number of clusters k: d ∈ {8, 16, 32} and k ∈ {0.8 × K, K, 1.2 × K} where K is the number of clusters in the ground truth.

Propagation [START_REF] Raghavan | Near linear time algorithm to detect community structures in large-scale networks[END_REF]. This method iterates over the nodes and assigns them the most represented label among their neighbors until convergence. The label of each node is initialized with its index.

We report the best ARI score obtained by each algorithm in Table 4.1 and the corresponding runtimes in Table 4.2. All computations are performed in the laptop setup. We use our Scikitnetwork implementations for all algorithms.

Our first observation is the generally good performance of all Louvain methods, except for Potts, combined with some of the shortest running times. Note that, among these methods, DiLouvain performs best on the directed graph WV and Bilouvain performs best on the bipartite graph NG, even though it is outperformed by GSVD. On the other hand, the standard Laplacian Eigenmaps

The goal is to recover these labels without supervision. We use the V-measure score to compare two different pipelines.

K-means

The data is centered to have zero mean, then we apply a principal component analysis with whitening. In the L2 version, we apply k-means directly on the output of the PCA. In the cosine version, we first normalize each row so that it has an l 2 norm equal to 1.

Louvain We apply a principal component analysis without whitening and construct a nearestneighbor graph on which we apply DiLouvain, as the resulting graph is directed by construction.

We perform a grid-search optimization for the following hyper-parameters: dimension of the PCA n dim ∈ {8, 16, 32}, number of clusters for k-means k ∈ {8, 9, 10, 11, 12}, number of nearest neighbors to construct the graph n neigh ∈ {1, 3, 5} and resolution for DiLouvain γ ∈ {0.1, 0.5, 1.}. We report the best score for each pipeline in Table 4 Chapter 5

Semi-supervised node classification

This chapter is dedicated to the application of heat diffusion to node classification. Heat diffusion is a very intuitive and interpretable physical process which has inspired the heat kernel, a standard tool in graph mining with many applications, including semi-supervised classification of nodes. The idea is to compute the system state at some finite time t after free diffusion from a set of hot source nodes. In this work, we use a different diffusion based on the solution to a Dirichlet problem, with both hot and cold sources. These new boundary conditions lead to meaningful steady states, which can be exploited for node classification.This chapter covers the work presented in De Lara and Bonald [2020].

Introduction

Heat diffusion, describing by the evolution of temperature T in an isotropic material, is governed by the following heat equation, initially developed by Joseph Fourier:

∂T ∂t = α∆T, (5.1)
where α is the thermal conductivity of the material and ∆ is the Laplace operator. In steady state, this equation simplifies to ∆T = 0 and the function T is said to be harmonic. The Dirichlet problem consists in finding the equilibrium in the presence of boundary conditions, that is when the temperature T is fixed on the boundary of the region of interest. Figure 5.1 is a small tribute to these scientists. In this chapter, we show how to apply a discrete version of the Dirichlet problem to the semisupervised classification of the nodes of a graph: given labels known for some nodes of the graph, referred to as the seeds, how to infer the labels of the other nodes? The number of seeds is typically small compared to the total number of nodes (e.g., 1%), hence the name of semisupervised classification.

We propose to solve one Dirichlet problem per label by setting the temperatures of the seeds accordingly, in a one-against-all strategy, and to classify the nodes with respect to the deviation of temperature to the mean. We prove the consistency of our algorithm on a simple block model and its efficiency through experiments on real graphs.

The application of the heat equation to graph mining is well known in the literature. It is usually associated with kernel learning [START_REF] Kolla | Diffusion kernels on graphs and other discrete structures[END_REF] and referred to as heat kernel. It has been used to many different tasks like pattern matching [START_REF] Thanou | Learning heat diffusion graphs[END_REF], node ranking [START_REF] Ma | Mining social networks using heat diffusion processes for marketing candidates selection[END_REF][START_REF] Ma | Mining web graphs for recommendations[END_REF], graph embedding [START_REF] Donnat | Learning structural node embeddings via diffusion wavelets[END_REF], graph clustering [START_REF] Tremblay | Graph wavelets for multiscale community mining[END_REF], and semi-supervised classification [START_REF] Zhu | Semi-supervised learning with graphs[END_REF][START_REF] Merkurjev | A semi-supervised heat kernel pagerank mbo algorithm for data classification[END_REF][START_REF] Berberidis | Adadif: Adaptive diffusions for efficient semi-supervised learning over graphs[END_REF], Li et al., 2019]. In most cases, the learning task relies on the transient state of the diffusion process, using hot sources only. Our approach is different in that we solve one Dirichlet problem per label, using a one-against-all strategy, with both hot sources (the seeds of the considered label) and cold sources (the seeds of the other labels). The algorithm is parameterfree, unlike existing techniques based on the heat kernel, whose performance critically depends on some time parameter used to stop the diffusion. Our theoretical analysis also shows that temperature centering is critical, i.e., classification must rely on temperature deviations to the mean for each Dirichlet problem.

The rest of this chapter is organized as follows. In section 5.2, we introduce the Dirichlet problem on graphs. Section 5.3 describes our algorithm for node classification. The analysis showing the consistency of our algorithm on a simple block model is presented in section 5.4. Finally, Section 5.5 presents the experiments.

Dirichlet problem on graphs

In this section, we introduce the Dirichlet problem on graphs, show its interpretation in terms of random walk in the graph and characterize the solution. The difference with the heat kernal is also highlighted.

Heat equation

Consider a graph G with n nodes. We first assume that the graph is undirected and unweighted.

We assume that nodes are indexed from 1 to n and denote by A the corresponding adjacency matrix. This is a symmetric, binary matrix. Let d = A1 be the degree vector, which is assumed positive, and D = diag(d). The Laplacian matrix is defined by

L = D -A.
Now let S be some strict subset of {1, . . . , n} and assume that the temperature of each node i ∈ S is set at some fixed value T i . We are interested in the evolution of the temperatures of the other nodes. Heat exchanges occur through each edge of the graph proportionally to the temperature difference between the corresponding nodes. Then,

∀i / ∈ S, dT i dt = n ∑ j=1 A ij (T j -T i ), that is ∀i / ∈ S, dT i dt = -(LT) i ,
where T is the vector of temperatures. This is the heat equation in discrete space, where -L plays the role of the Laplace operator in (5.1). At equilibrium, T satisfies Laplace's equation:

∀i / ∈ S, (LT) i = 0.

(5.2)

We say that the vector T is harmonic. With the boundary conditions T i for all i ∈ S, this defines a Dirichlet problem in discrete space.

Random walk

Consider a random walk in the graph G with a probability of moving from node i to node j equal to A ij /d i . Let X 0 , X 1 , X 2 , . . . be the sequence of nodes visited by the random walk. This defines a Markov chain on {1, . . . , n} with transition matrix P = D -1 A.

Observing that L = D(I -P), Laplace's equation can be written equivalently ∀i / ∈ S, T i = (PT) i .

(5.3)

In other words, the temperature of each node i ∈ S at equilibrium is the average of the temperature of its neighbors. This implies in particular that the solution to the Dirichlet problem is unique, provided that the graph is connected. For the sake of completeness, we provide a proof of this result in the considered discrete case:

Proposition 11. If the graph is connected, there is at most one solution to the Dirichlet problem.

which is similar to iteration (5.3) but in the absence of boundary. In both (5.9) and (5.10), the number of iterations, say N, plays a key role, as the state T converges to the equilibrium with uniform temperatures when N → +∞.

Note that, in the case of the diffusion in discrete time, it is the weighted average temperature T = 1 w ∑ i d i T i that is preserved over time:

T(t + 1) = d w T(t + 1) = d w P T(t) = d w T(t) = T(t).

Extensions

All results extend to weighted graphs, with a positive weight assigned to each edge; this weight can then be interpreted as the thermal conductivity of the edge in the diffusion process. The results also apply to directed graphs. Indeed, a directed graph G of n nodes, with adjacency matrix A, can be considered as a bipartite graph of 2n nodes, with adjacency matrix:

0 A A 0
The diffusion can be applied to this bipartite graph, which is undirected. Observe that each node of the directed graph G is duplicated in the bipartite graph and is thus characterized by 2 temperatures, one as heat source (outgoing edges) and one as heat destination (incoming edges). It is not necessary for the directed graph to be strongly connected; only the associate bipartite graph needs to be connected.

Node classification algorithm

In this section, we introduce a node classification algorithm based on the Dirichlet problem.

The objective is to infer the labels of all nodes given the labels of a few nodes called the seeds.

Our algorithm is a simple modification of the popular method proposed by [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF]. Specifically, we propose to center temperatures before classification.

Binary classification

When there are only two different labels, the classification can be done by solving one Dirichlet problem. The idea is to use the seeds with the first label as hot sources, setting their temperature at 1, and the seeds with the second label as cold sources, setting their temperature at 0. The solution to this Dirichlet problem gives temperatures between 0 and 1.

A natural approach, proposed by [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF], consists in assigning label 1 to all nodes with temperature above 0.5 and label 2 to other nodes. The analysis of section 5.4 suggests that it is preferable to set the threshold to the mean temperature,

T = 1 n n ∑ i=1 T i .
Specifically, all nodes with temperature above T are assigned label 1, the other are assigned label 2. Equivalently, temperatures are centered before classification: after centering, nodes with positive temperature are assigned label 1, the others are assigned label 2.

(temperature 0). After centering the temperatures (so that the mean temperature of each diffusion is equal to 0), each node is assigned the label that maximizes its temperature. This algorithm, we refer to as Dirichlet classifier, is parameter-free. T S = 0.

3:

for i ∈ S do 4:
if y i = k then 5:

T S i = 1. 6: end if 7:
end for 8:

T (k) ← Dirichlet(S, T S ).

9: k) mean(T (k) ) 10: end for 11: for i ∈ S do 12:

∆ (k) ← T (
x i = arg max k=1,...,K (∆ (k) i ) 13: end for 14: return x, labels of nodes outside S

The key difference with the vanilla method lies in temperature centering (line 9 of the algorithm). Another variant proposed by [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF] consists in rescaling the temperature vector by the weight of the considered label in the seeds (see equation ( 9) in their paper).

Time complexity

The time complexity depends on the algorithm used to solve the Dirichlet problem. We here focus on the approximate solution by successive iterations of (5.6). Let m be the number of nonzero entries in A. Using the Compressed Sparse Row format for the adjacency matrix, each matrix-vector product has a complexity of O(m). The complexity of Algorithm 5.1 is then O(NKm), where N is the number of iterations. Note that the K Dirichlet problems are independent and can thus be computed in parallel.

Analysis

In this section, we prove the consistency of Algorithm 5.1 on a simple block model. In particular, we highlight the importance of temperature centering in the analysis.

Block model

Consider a graph of n nodes consisting of K blocks of respective sizes n 1 , . . . , n K , forming a partition of the set of nodes. There are s 1 , . . . , s K seeds in these blocks, which are respectively assigned labels 1, . . . , K. Intra-block edges have weight p and inter-block edges have weight q. We expect the algorithm to assign label k to all nodes of block k whenever p > q, for all k = 1, . . . , K.

Dirichlet problem

Consider the Dirichlet problem when the temperature of the s 1 seeds of block 1 is set to 1 and the temperature of the other seeds is set to 0. We have an explicit solution to this Dirichlet problem.

Lemma 7. Let T k be the temperature of non-seed nodes of block k at equilibrium. We have:

(s 1 (p -q) + nq)T 1 = s 1 (p -q) + n Tq, (s k (p -q) + nq)T k = n Tq k = 2, . . . , K,
where T is the average temperature, given by:

T = s 1 n n 1 (p -q) + nq s 1 (p -q) + nq / 1 - K ∑ k=1 (n k -s k )q s k (p -q) + nq .
Proof. In view of (5.2), we have:

(n 1 (p -q) + nq)T 1 = s 1 p + (n 1 -s 1 )pT 1 + ∑ j =1

(n js j )qT j ,

(n k (p -q) + nq)T k = s 1 q + (n k -s k )pT k + ∑ j =k
(n js j )qT j , k = 2, . . . , K.

We deduce:

(s 1 (p -q) + nq)T 1 = s 1 p + Vq, (s k (p -q) + nq)T k = s 1 q + Vq k = 2, . . . , K. with V = K ∑ j=1
(n js j )T j .

The proof then follows from the fact that

n T = s 1 + K ∑ j=1
(n js j )T j = s 1 + V.

Classification

We now state the main result of the chapter: the Dirichlet classifier is a consistent algorithm for the block model, in the sense that all nodes are correctly classified whenever p > q.

Theorem 6. If p > q, then x i = k for all non-seed nodes i of each block k, for any parameters n 1 , . . . , n K (block sizes) and s 1 , . . . , s K (numbers of seeds).

Proof. Let δ

(1) k = T k -T be the deviation of temperature of non-seed nodes of block k for the Dirichlet problem associated with label 1. In view of Lemma 7, we have:

(s 1 (p -q) + nq)δ (1) 1 = s 1 (p -q)(1 -T), (s k (p -q) + nq)δ (1) k = -s k (p -q) T k = 2, . . . , K,
For p > q, using the fact that T ∈ (0, 1), we get δ

(1)

1 > 0 and δ

(1) k

< 0 for all k = 2, . . . , K. By symmetry, for each label l = 1, . . . , K, δ (l) l > 0 and δ (l) k < 0 for all k = l. We deduce that for each block k, x i = arg max l δ (l) k = k for all non-seed nodes i of block k.

Observe that the temperature centering is critical for consistency. In the absence of centering, non-seed nodes of block 1 are correctly classified if and only if their temperature is the highest in the Dirichlet problem associated with label 1. In view of Lemma 7, this means that for all k = 2, . . . , K,

s 1 q n 1 (p -q) + nq s 1 (p -q) + nq + s 1 (p -q) 1 - K ∑ j=1 (n j -s j )q s j (p -q) + nq > s k q n k (p -q) + nq s k (p -q) + nq .
This condition might be violated even if p > q, depending on the parameters n 1 , . . . , n K and s 1 , . . . , s K . In the practically interesting case where s 1 << n 1 , . . . , s K << n K for instance (low fractions of seeds), the condition requires:

s 1 (n 1 (p -q) + nq) > s k (n k (p -q) + nq).
For blocks of same size, this means that only blocks with the largest number of seeds are correctly classified. The classifier is biased towards labels with a large number of seeds. This sensitivity of the vanilla algorithm to the label distribution of seeds will be confirmed in the experiments on real graphs.

Experiments

In this section, we show the impact of temperature centering on the quality of classification using both synthetic and real data. First, in Sections 5.5.1 and 5.5.2, we only focus on 3 algorithms: the vanilla algorithm (without temperature centering), the weighted version proposed by [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF] (also without temperature centering) and our algorithm (with temperature centering). Then, in Section 5.5.3, we provide a general benchmark of classification methods to assess the performance of our algorithm.

Synthetic data

We first use the stochastic block model (SBM) [Airoldi et al., 2008] to generate graphs with an underlying structure in clusters. This is the stochastic version of the block model used in the analysis. There are K blocks of respective sizes n 1 , . . . , n K . Nodes of the same block are connected with probability p while nodes in different blocks are connected probability q. We denote by s k the number of seeds in block k and by s the total number of seeds.

We first compare the performance of the algorithms on a binary classification task (K = 2) for a graph of n = 10 000 nodes with p = 10 -3 and q = 10 -4 , in two different settings:

• Seed asymmetry: Both blocks have the same size n 1 = n 2 = 5000 but different numbers of seeds, with s 1 /s 2 ∈ {1, 2, . . . , 10} and s 2 = 250 (5% of nodes in block 2).

• Block size asymmetry: The blocks have different sizes with ratio n 1 /n 2 ∈ {1, 2, . . . , 10} and seeds in proportion to these sizes, with a total of s = 1 000 seeds (10% of nodes).

For each configuration, the experiment is repeated 10 times. Randomness comes both from the generation of the graph and from the selection of the seeds. We report the F1-scores in Figure 5.3 (average ± standard deviation). Observe that the variability of the results is very low due to the relatively large size of the graph. As expected, the centered version is much more robust to both types of asymmetry. Besides, in case of asymmetry in the seeds, the weighted version of the algorithm tends to amplify the bias and leads to lower scores than the vanilla version. We show in Figure 5.4 the same type of results for K = 10 blocks and p = 5.10 -2 . For the block size asymmetry, the size of blocks 2, . . . , 10 is set to 1 000.

Real data

We use datasets from the NetSet1 and SNAP 2 collections (see These datasets can be categorized into 3 groups:

• Citations networks: Cora (CO) and CiteSeer (CS) are citation networks between scientific publications. These are standard datasets for node classification [START_REF] Fey | Splinecnn: Fast geometric deep learning with continuous b-spline kernels[END_REF][START_REF] Huang | Adaptive sampling towards fast graph representation learning[END_REF][START_REF] Wijesinghe | Dfnets: Spectral cnns for graphs with feedback-looped filters[END_REF].

• Wikipedia graphs: Wikipedia for schools (WS) [START_REF] Haruechaiyasak | Article recommendation based on a topic model for wikipedia selection for schools[END_REF], Wikipedia vitals (WV) and Wikilinks (WL) are graphs of hyperlinks between different selections of Wikipedia pages. In WS and WV, pages are labeled by category (People, History, Geography...). For WL, pages are labeled through clusters of words used in these articles. As these graphs are directed, we use the extension of the algorithm described in §5.2.5, with nodes considered as heat sources.

• Social networks: DBLP and Amazon are social networks with partial ground-truth communities [START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF]. As nodes are partially labeled and some nodes have several labels, the results for these datasets are presented separately, with specific experiments based on binary classification.

For the citation networks and the Wikipedia graphs, we compare the classification performance of the algorithms in terms of macro-F1 score and two seeding policies:

• Uniform sampling, where seeds are sampled uniformly at random.

• Degree sampling, where seeds are sampled in proportion to their degrees.

In both cases, the seeds represent 1% of the total number of nodes in the graph. The process is repeated 10 times for each configuration. We do not display the results for the weighted version of the algorithm as they are very close to those obtained with the vanilla algorithm.

We report the results in Tables 5.2 and 5.3 for uniform sampling and degree sampling, respectively. We see that centered version outperforms the vanilla one by a significant margin. Graph Convolutional Network (GCN) [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]. These neural networks perform nonlinear message passing between the nodes of the graph in order to minimize both a classification loss and an embedding loss. Here, we use the original network design with one hidden layer for which we test dimensions h ∈ {16, 32, 64}. The network is implemented with the DGL Python library [START_REF] Wang | Deep graph library: Towards efficient and scalable deep learning on graphs[END_REF] and trained for 500 epochs using Adam stochastic optimization [START_REF] Kingma | A method for stochastic optimization[END_REF]. In the absence of exogenous node features, we use a one-hot encoding of node degrees as features, as suggested by [START_REF] Xu | How powerful are graph neural networks?[END_REF], combined with binning in dimension n feat ∈ {16, 32, 64}.

For Dirichlet and the standard Heat Kernel based on the free diffusion, we test N ∈ {5, 15, 25}.

We report the best scores obtained by each algorithm in Table 5.5. We also report the corresponding median runtimes, computed in the server setup, in Table 5.6. Computations longer than 6 hours trigger a Time Out.

As we can see, Dirichlet and Pagerank perform very well on most graphs, however, recall that Pagerank has an extra hyper-parameter to tune with respect to Dirichlet. GSVD also obtains good results, especially on the bipartite and directed graphs which is coherent with the results of Chapter 4. On the other hand, despite being the fastest algorithm on all datasets (probably because of the small world property and the hard stopping criterion), Propagation yields some of the poorest results, especially on the bipartite and directed graphs. Finally, GCN gives globally very poor results while being slower than other algorithms by several orders of magnitude. Some explanation about this poor performance can be found in [START_REF] Hou | Measuring and improving the use of graph information in graph neural networks[END_REF], it is possible that the GCN overfits the node attributes instead of actually learning from the graph structure. Secondly, these representations are index dependent: up to indexing of its nodes, a same graph admits n! equivalent representations. In a classification task, the label of a graph is independent from the indices of its nodes, so the model used for prediction should be invariant to node ordering as well. Some neural network architectures have been specifically designed to tackle this issue for standard vector data such as clouds of points [Qi et al., 2017, Segol and[START_REF] Segol | On universal equivariant set networks[END_REF]. However, these networks are not able to learn directly from graph data.

Section 6.2 presents some classic techniques to overcome these difficulties. Besides, due to this specific type of input, graph classification is not included in Scikit-network. However, our paper [START_REF] Lara | A Simple Baseline Algorithm for Graph Classification[END_REF] has been implemented in another open-source graph library, Karate Club [START_REF] Rozemberczki | An api oriented open-source python framework for unsupervised learning on graphs[END_REF].

Related work

Graph classification methods can schematically be divided into three categories: graph kernels, sequential methods and embedding methods. In this section, we briefly present these different approaches, focusing on methods that only use the structure of the graph and no exogenous information, such as node features, to perform classification as we only want to compare the capacity of the algorithms to capture structural information.

Kernel methods Kernel methods [Nikolentzos et al., 2017a,b, 2018, Neumann et al., 2016] perform pairwise comparisons between the graphs of the dataset and apply a classifier, usually a support vector machine (SVM), on the similarity matrix. In order to maintain the number of comparisons tractable when the number of graphs is large, they often use Nyström algorithm [START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF] to compute a low rank approximation of the similarity matrix. The key is to construct an efficient kernel that can be applied to graphs of varying sizes and captures useful features for the downstream classification.

Sequential methods Some methods tackle the varying sizes of graphs by processing them as a sequence of nodes. Earliest models used random walk based representations [START_REF] Callut | Classification in graphs using discriminative random walks[END_REF][START_REF] Xu | Protein classification using random walk on graph[END_REF]. More recently, [START_REF] Jin | Learning graph-level representations with gated recurrent neural networks[END_REF] or [START_REF] You | Graphrnn: A deep generative model for graphs[END_REF] transform a graph into a sequence of fixed size vectors, corresponding to its nodes, which is fed to a recurrent neural network. The two main challenges in this approach are the design of the embedding function for the nodes and the order in which the embeddings are given to the recurrent neural network. The algorithm of Section 6.4 falls into this category.

Embedding methods Embedding methods [START_REF] Gomez | Dynamics based features for graph classification[END_REF][START_REF] Barnett | Feature-based classification of networks[END_REF][START_REF] Dutta | High order stochastic graphlet embedding for graph-based pattern recognition[END_REF][START_REF] Narayanan | graph2vec: Learning distributed representations of graphs[END_REF], derive a fixed number of features for each graph which is used as a vector representation for classification. Even though deriving a good set of features is often a difficult task, this approach has the benefit of being compatible with any standard classifier in a plug and play fashion (SVM, random forest, multilayer perceptron...). The model of Section 6.3 belongs to this class of methods as it relies on spectral features of the graph.

Spectral method

In this section, we propose a simple and fast algorithm based on the spectral decomposition of graph Laplacian to perform graph classification and get a first reference score for a dataset. We show in Section 6.5 that this method obtains competitive results compared to state-of-the-art algorithms. Note that, even though we present a very simple algorithm here, spectral methods have been incorporated into more complex ones with good experimental results [Chauhan et al., 2020, Verma and[START_REF] Verma | Hunt for the unique, stable, sparse and fast feature learning on graphs[END_REF].

Let G = (V, E) be an undirected and unweighted graph and A ∈ {0, 1} n×n its Boolean adjacency matrix with respect to an arbitrary indexing of the nodes. G is assumed to be connected, otherwise, we extract its largest connected component. Let D = diag(A1) be the matrix of node degrees, the normalized Laplacian of G is defined as

L = I -D -1/2 AD -1/2 . (6.1)
We use the k smallest positive eigenvalues of L in ascending order as input of the classifier: X = (σ 1 , . . . , σ k ).

If the graph has less than k nodes, we use right zero padding to get a vector of appropriate dimensions: X = (σ 1 , . . . , σ n-1 , 0, . . . , 0).

The normalized Laplacian matrix of a graph is a well-known object in spectral learning [Belkin andNiyogi, 2002a, Kamvar et al., 2003]. However, for node clustering or classification most of the attention is usually directed to its eigenvectors and not its spectrum. A major benefit of the ordered spectrum representation for graph classification is that it does not depend on the indexing of the nodes.

Some Laplacian eigenvalues properties

The eigenvalues of the normalized Laplacian matrix lie between 0 and 2. Such a property is very convenient for the downstream use of a standard classifier without heavy rescaling or preprocessing. The multiplicity of the eigenvalue 0 corresponds to the number of connected components in the graph, hence the omission of σ 0 in our representation as we only consider the largest connected component. Other values are also known to denote the presence of specific structures in the graph [START_REF] Chung | Spectral graph theory[END_REF].

For example, an eigenvalue equal to 2 denotes a bipartite structure.

Physical interpretations

In Bonald et al. [2018b], each eigenvalue of the Laplacian corresponds to the energy level of a stable configuration of the nodes in the embedding space. The lower the energy, the stabler the configuration. In [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF], these eigenvalues correspond to frequencies associated to a Fourier decomposition of any signal living on the vertices of the graph. Thus, the truncation of the Fourier decomposition acts as low-pass filter on the signal. Characterizing a graph by the smallest eigenvalues of its normalized Laplacian is thus comparable to characterizing a melody by its lowest fundamental frequencies.

Finally, there have been some attempts to connect spectral decomposition to graph isomorphism [START_REF] Van Dam | Which graphs are determined by their spectrum? Linear Algebra and its applications[END_REF]Haemers, 2003, Kolla et al., 2017], however, to the best of our knowledge, this is still an open problem.

The choice of the classifier is left to the discretion of the user. In our experiments, we chose a random forest classifier (RFC) which offers a good computational speed versus accuracy tradeoff.

An illustration of the model is proposed in figure 6.2. Figure 6.2: Spectral model. L refers to equation (6.1) and ĉ is the predicted class.

Variational RNN method

In this section, we propose a method to sequentially embed graph information in order to perform classification. By construction, this recurrent graph classifier overcomes the common difficulties listed in Section 6.1. Besides, we propose to use an additional node prediction block to help the model to capture the intrinsic structure of the graphs. The complete model is denoted variational recurrent graph classifier (VRGC).

We propose to use a sequential approach to embed graphs with a variable number of nodes and edges into a vector space of a chosen dimension. This latent representation is then used for classification. Node index invariance is approximated through specific pre-processing and aggregation.

Let G = (V, E) be an undirected and unweighted graph with V a set of nodes and E a set of edges. The graph G can be represented, modulo any permutation π over its nodes i ∈ {1, . . . , n}, by its boolean adjacency matrix A π such that A π ij = 1 if nodes indexed by i and j are connected in the graph and A π ij = 0 otherwise. We use this adjacency matrix as a raw representation of the graph.

The proposed VRGC model is composed of three main parts: node ordering and embedding, classification and regularization with variational auto-regression (VAR). See Figure 6.3 for an illustration. Node ordering and embedding Before being processed by the neural network, the adjacency matrix of a graph is transformed on-the-fly [START_REF] You | Graphrnn: A deep generative model for graphs[END_REF]. First, a node is selected at random and used as root for a breadth first search (BFS) over the graph, as illustrated in Figure 6.4. The rows and columns of the adjacency matrix are then reordered according to the sequence of nodes returned by the BFS. Next, each row i (corresponding to the i th node in the BFS ordering) is truncated to keep only the connections of node i with the min(i, d) nodes that preceded in the BFS. This way, each node is d-dimensional, and each truncated matrix is zero-padded in order to have dimensions (n G , d), with n G the size of the larger graph in the dataset. This can be seen as applying a sliding window of width d to the adjacency matrix. See Figure 6.5 for an illustration.

After node ordering and pre-embedding, each graph is processed as a sequence of ddimensional nodes by a gated recurrent unit (GRU) neural network [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF]. The GRU is a special RNN able to learn long term dependencies by solving vanishing gradient effect. The choice of GRU over Long Short Term Memory networks is arbitrary as they have equivalent long-term modeling power [START_REF] Chung | Empirical evaluation of gated recurrent neural networks on sequence modeling[END_REF]. In order to help the recurrent network training, we propose to add a simple fully connected network between pre-embedding Figure 6.5: Truncation procedure from the root node R.

and recurrent embedding. Therefore, the node will be presented to the GRU in the shape of continuous vectors instead of binary adjacency vectors.

Finally the GRU sequentially embeds each node i by using i -1 and information contained in a memory cell h i-1 that theoretically embeds all previously seen information. The embedded node sequence {h i } n G i=1 then feeds both the VAR and the classifier as discussed in subsequent sections. See Figure 6.6 for an illustration. Classification After the embedding step, we use an additional GRU dedicated to classification that takes {h i } n G i=1 as input. Its last memory cell, denoted hn G , feeds a softmax multilayer perceptron (MLP) which performs class prediction. Formally, let c be the class index, the classifier is trained by minimizing the cross-entropy loss between ground-truth and p(G, r) the softmax class membership probability vector for a given graph G that has been sorted by a BFS rooted with node r. We call this objective term L classi f . As discriminating patterns might be spread across the whole graph, the network is required to model long-term dependencies. By construction, GRUs have such ability. See Figure 6.7 for an illustration.

Regularization with variational auto-regression

As the structure of a graph is the concatenation of the interactions between all nodes and their respective neighbors, learning a good representation without using node attributes requires for the model to capture the structure of the graph while classifying. Accordingly, we add an auto-regression block to our model. Given a node to process, the network makes a prediction for the neighborhood of the next node. Multitask learning is a powerful leverage to learn rich representation in NLP [START_REF] Sanh | A hierarchical multi-task approach for learning embeddings from semantic tasks[END_REF]. In particular, such representation for sequence classification has already been used for sentiment analysis [START_REF] Latif | Variational autoencoders for learning latent representations of speech emotion[END_REF][START_REF] Xu | Variational autoencoder for semi-supervised text classification[END_REF]. We use a variational auto-encoder (VAE) [START_REF] Kingma | Auto-encoding variational bayes[END_REF] to learn a representation of each node i given h i-1 . The first layers of the encoder are shared with the classifier and corresponds to the graphs preprocessing (blue part in Figures 6.3 and 6.6). The subsequent encoder layers, the latent sampling and the decoder constitute the VAR. For each graph G with embedded nodes i ∈ {1, . . . , n G }, the fully connected variational auto-encoder takes the latent states {h i } i∈{1,...,n G } as input. Let {z i } i∈{1,...,n G } be the latent random variables for the following model p(G, z|h) = p(i = 1, z 1 ) n G ∏ i=2 p θ (i|z i )q φ (z i |h i-1 ).

In practice, p θ and q φ are modelled by neural networks parameterized by θ and φ, which require differentiable functions for training. However, p θ (i|z i ) models a binary adjacency vector representing the connections between node i and previously visited nodes j < i. Therefore, we use sigmoid continuous relaxation to train our model, and hard binary sampling at test time. We use a Gaussian variational posterior distribution. Training is done by maximizing the variational lower bound of the log-likelihood of the observation as in Kingma's VAE.

The VAE-like loss for VAR regularization is the following:

L pred = E p d (G) n G ∑ i=2 KL q φ (z i |h i-1 )||q(z i ) -E p d (G) n G ∑ i=2 E q φ (z i |h i-1 ) [log p θ (i|z i )] ,
which is a lower bound of the negative marginal log-likelihood E p d (G) [log p θ (G)]. p θ and q φ are the respective densities of i|z and z|h, whose distribution are parameterized by θ and φ, respectively. KL denotes the Kullback-Leibler divergence, p d is the empirical distribution of G and q(z i ) is the density of the prior distribution of latent variables {z i } n G i=2 . We chose the standard Gaussian prior for q(z i ).

The regularization part is illustrated in Figure 6.8.

In the end, the model is trained by minimizing the total loss L = L classi f + αL pred , where α is a hyper-parameter.

Aggregation of the results at test time

The node ordering step introduces randomness to our model. On the one hand, it helps learn more general graph representations during the training phase, but on the other hand, it might produce different outputs for the same graph during the testing phase, depending on the root of the BFS. In order to counter this side effect, we add the following aggregation step for the testing phase. Each graph is processed N times by the model with N different roots for BFS ordering. The N class membership probability vectors are

Experiments

This section is dedicated to our experiments on graph classification. First, we introduce the datasets in Section 6.5.1 and the results in sections 6.5.2, 6.5.3, 6.5.4 and 6.5.5.

Datasets

We evaluated our model against four standard datasets from biology: Mutag (MT), Enzymes (EZ), Proteins Full (PF) and National Cancer Institute (NCI1) [START_REF] Kersting | Benchmark data sets for graph kernels[END_REF]. 2016] (FB), Dynamic-Based Features [START_REF] Gomez | Dynamics based features for graph classification[END_REF] (DyF), Stochastic Graphlet Embedding [START_REF] Dutta | High order stochastic graphlet embedding for graph-based pattern recognition[END_REF] (SGE) and Family of Graph Spectral Distances [START_REF] Verma | Hunt for the unique, stable, sparse and fast feature learning on graphs[END_REF] (FGSD).

All values are directly taken from the aforementioned papers as they use a setup similar to ours. For algorithms presenting results with and without node features, we reported the results without node features. For those presenting results with several sets of hyper-parameters, we reported the results for the parameters that performed best on the largest number of datasets.

Results are reported in Table 6.3.

VRGC obtains state-of-the-art results on three out of these four datasets and the second best result on the fourth one. However, its standard deviation is usually higher than for other algorithms. The spectral method also gets competitive scores while it is much faster to train. As we can see, RFC provides the best results for all datasets except DD where MLP has an accuracy of 75.6 against 75.4. Our intuition to explain these good results is that the decision tree classifier, which is at the core of RFC, is an algorithm based on level thresholding. As explained in section 6.3, our embedding represents a sequence of energy levels, being above or below a certain level is thus likely to be meaningful for classification.

Influence of k for the spectral model

We experimented with different embedding dimensions for RFC: k ∈ {1, 5, 10, 25, 50}. The hyper-parameters are the same as in the other experiments. Results are reported in We see that even the first energy level is sufficient to obtain a non-trivial classification. k = 5 provides results competitive with the state of the art while k = 50 provides results relatively similar to k = avg(|V|). We did not experiment with larger values of k as it would mostly result into additional zero padding for most graphs. Note that, embedding all graphs for k = 50 took less than a minute in the laptop setting.

Node indexing invariance for VRGC

The model is designed to be independent from node ordering of the graph with respect to different BFS roots. Inputs representing the same graph (up to node ordering) should be close from one another in the latent embedding space. As the preprocessing is performed on each graph at each epoch, a same graph is processed many times by the model during training with different embeddings. This creates a natural regularization for the network.

We illustrate this in Figure 6.9. We embed five graphs from the dataset EZ with 20 different BFS each and plot their TSNE projections. As we can see, the projections corresponding to the same graphs form a heap in the low dimensional representation of the latent space. Chapter 7

Conclusions and perspectives

In this work, we have sketched a map of modern graph mining and related challenges.

First, we addressed the data representation problem. In particular, we explained how linear operators combined with advanced linear algebra techniques could help make various graph algorithms scale to industrial levels. The next step would be to study these techniques in the context of distributed data and computing power.

Then, we introduced Scikit-network, an open-source Python library for large scale graph mining inspired by Scikit-learn. In the future, we hope to keep integrating new features such as anomaly detection or link prediction and to continue improving the package.

We have provided a simple explanation for the well-known benefits of regularization on spectral embedding. Specifically, regularization forces the embedding to focus on the largest clusters, making the embedding more robust to noise. This result was obtained through the explicit characterization of the embedding for a simple block model and extended to bipartite graphs. An interesting perspective of this work is the extension to stochastic block models, using, for instance, the concentration results proved in [START_REF] Lei | Consistency of spectral clustering in stochastic block models[END_REF][START_REF] Le | Concentration and regularization of random graphs[END_REF]. Another area of interest is the impact of regularization on other downstream tasks such as link prediction. Finally, we would like to further explore the impact of the regularization parameter while exploiting the theoretical results presented in this section.

We have proposed a novel embedding based on the Generalized Singular Value Decomposition, which applies to undirected, directed, and bipartite graphs. We have explained how the distances in the embedding space could be easily interpreted in terms of neighborhood distribution for the best rank-k approximation of the graph. Efficiency of this embedding has been demonstrated on real datasets for both node clustering and node classification. For this new and standard spectral algorithm, we cannot stress enough the importance of normalization.

Projecting each node vector onto the unit sphere drastically improves performance for all the downstream tasks we have experimented with.

We showed how the classic modularity function for graph clustering extends to a large family of functions based on graph sampling. We then explained how optimizing such a function corresponds to fitting a simple average model to the data through the Kullback Leibler divergence. We illustrated how the Louvain heuristic can be used to perform greedy maximization for all these modularity functions. Future work could include the study of soft-clustering membership based on Louvain algorithms as graph embeddings.

We have proposed a novel approach to node classification based on heat diffusion. Specifically, we propose to center the temperatures of each solution to the Dirichlet problem before classification. We have proved the consistency of this algorithm on a simple block model and we have shown that it drastically improves classification performance on real datasets with respect to the vanilla version. In future work, we plan to extend this algorithm to soft classification, using the centered temperatures to get a confidence score for each node of the graph. Another interesting research perspective is to extend our proof of consistency of the algorithm to stochastic block models.

We explained the specific challenges of the graph classification problem. Then, we described and evaluated two different approaches to overcome them. However, we believe that further refinement of these algorithms should be more context dependent and include domain knowledge. For example, classifying molecular compounds is not the same as identifying malware programs.

Finally, we wish to study the impact simplex projection for membership matrices. In this work, we always relied on naive scaling to normalize either soft-clustering memberships or classification confidence. However, such scaling does not yield the closest simplex point to the membership, hence this new projection could lead to different results which are worth exploring. In particular, the simplex projection is known to yield sparser vectors than the naive scaling.

Figure 1 . 1 :

 11 Figure 1.1: Some toy graphs.

Figure 1 . 2 :

 12 Figure 1.2: Node embedding.

Figure 1 . 3 :

 13 Figure 1.3: Node clustering.

Figure 1 . 4 :

 14 Figure 1.4: Node classification.

Figure 1 . 5 :

 15 Figure 1.5: Node ranking.

Figure 1 . 6 :

 16 Figure 1.6: Computation times in seconds for partial decomposition of matrices.

Figure 1 . 7 :Figure 1

 171 Figure 1.7: Computation times in seconds for different values of the number of components, k.

Figure 2

 2 Figure 2.1: Scikit-network's logo.

Figure 2 . 2 :

 22 Figure 2.2: Visualization of a graph and a dendrogram as SVG images.

Figure 3 . 1 :

 31 Figure 3.1: Virtual mechanical system.

Figure 3 . 2 :

 32 Figure 3.2: Regularization of a toy graph.

  α = 1.

Figure 3

 3 Figure 3.3: Embeddings of a toy graph. The size of a marker is proportional to the size of the corresponding clique.

  where || • || σ denotes the spectral norm of R n 1 ×n 2 , given by the largest singular value. The Froebenius norm is denoted by || • || F .

  Figure 3.4: Low-rank approximation of a toy graph.

Figure 3

 3 Figure 3.5: Co-embedding of movies (blue) and actors (red).

  Figure 3.6: A toy bipartite graph.

Figure 3 . 7 :

 37 Figure 3.7: Aggregate graphs resulting from spectral clustering. Node sizes are proportional to cluster sizes.

Figure 4

 4 Figure 4.1: Movie-actor bipartite graph.

Figure 5 . 1 :

 51 Figure 5.1: Left to right: Joseph Fourier, Pierre-Simon de Laplace and Johann P.G.L. Dirichlet.

Algorithm 5. 1

 1 Dirichlet classifierRequire: Seed set S and associated labels y ∈ {1, . . . , K}.1: for k in {1, . . . , K} do 2:

Figure 5 . 3 :

 53 Figure 5.3: Binary classification performance on the SBM.

Figure 5 . 4 :

 54 Figure 5.4: Multi-label classification performance on the SBM.

Figure 5 . 5 :

 55 Figure 5.5: Impact of the fraction of seeds for label 1 on macro-F1 score (mean ± standard deviation).

Figure 6 . 1 :

 61 Figure 6.1: Trimethylamine molecule represented as a graph.

Figure 6

 6 Figure 6.3: Macroscopic representation of VRGC.

Figure 6

 6 Figure 6.4: Node indexing using breadth first search. R: root node.

Figure 6 . 6 :

 66 Figure 6.6: Node ordering and embedding.

Figure 6 . 7 :

 67 Figure 6.7: Classification block.

Figure 6 . 8 :

 68 Figure 6.8: Regularization with VAR plus final aggregation.

Figure 6

 6 Figure 6.9: TSNE projection of the latent state preceding classification.

Table 1

 1

	task	COO CSR CSC DOK
	Memory space (Mo)	6.8	3.8	3.8	46.8
	Node insertion (ms)	5.7	11	11	108
	Node removal (ms)	-	2.6	2.5	743
	Edge insertion (ms)	-	0.16 0.15	0.01
	Edge removal (ms)	-	0.7	0.7	0.01
	Out-neighbors access (µs)	301	0.6	-	1992
	In-neighbors access (µs)	331	-	0.6	3220

.2: Performance on Wikivitals.

  Table 1.4 presents some characteristics of these datasets.

	code	name	n		m	density
	MK	Kangaroo	2.10 1	2.10 2	10 -1
	JZ	Jazz musicians	2.10 2	5.10 3	10 -1
	Shf	Hamsterster	2.10 3	2.10 4	10 -3
	RC	Reactome	6.10 3	3.10 5	10 -2
	GW	Gowalla	2.10 5	2.10 6	10 -4
	SK	Skitter	2.10 6	2.10 7	10 -5
	OR	Orkut	3.10 6	2.10 8	10 -5
		(a) Undirected graphs.
	code	name	n 1	n 2	m	density
	SW	Southern women 1	2.10 1	2.10 1	9.10 1	10 -1
	UL	Unicode languages	2.10 2	6.10 1	1.10 3	10 -1
	SX	Sexual escorts	1.10 4	6.10 3	4.10 4	10 -3
	Mut	MovieLens user-tag	4.10 3	2.10 4	4.10 4	10 -3
	R2	Reuters-21578	2.10 4	4.10 4	1.10 6	10 -3
	M3	MovieLens 10M	7.10 4	1.10 4	1.10 7	10 -2
	RE	Reuters	8.10 5	3.10 5	6.10 7	10 -4
		(b) Bipartite graphs.

Table 1

 1

.4: Some characteristics of the datasets.

Table 2 .

 2 1: Overview of graph software features. : Available. : Partial overlap. : Not available.

	Path. Module relying on SciPy for graph traversals: shortest paths, breadth first search,
	depth first search, diameter estimation, etc.
	• Topology. Module for topological structures extraction such that connected components,
	bipartite structure, k-core, clustering coefficient etc.
	• Visualization. Module for visualizing graphs and dendrograms in SVG (Scalable Vector
	Graphics) format. Examples are displayed in Figure 2.2.
	These modules are only partially covered by existing graph softwares (see Table 2.1). Another
	interesting feature of Scikit-network is its ability to work directly on bipartite graphs, represented
	by their biadjacency matrix.
		Modules	scikit-network NetworkX iGraph graph-tool
		Data
		Clustering
		Hierarchy
		Embedding
		Ranking
		Classification
		Path
		Topology
		Visualization
					Henri Matisse
					Gustav Klimt	Pablo Picasso
	Edouard Manet	Henri Matisse	Pablo Picasso	Paul Cezanne Pierre-Auguste Renoir
	Paul Cezanne			Egon Schiele	Claude Monet
	Claude Monet				Edouard Manet
	Pierre-Auguste Renoir	Vincent van Gogh	Vincent van Gogh
	Edgar Degas				Edgar Degas
					Egon Schiele
					Gustav Klimt
	Peter Paul Rubens		Leonardo da Vinci
	Michel Angelo		Rembrandt	Michel Angelo Rembrandt
					Peter Paul Rubens
		Leonardo da Vinci

  in the server setup. The graph has 3 072 441 nodes and 117 184 899 edges. As we can see, Scikit-network is highly competitive.

		scikit-network NetworkX iGraph graph-tool
	Louvain	139	1 978
	PageRank	48	236	45
	HITS	109	80	144
	Spectral	534
	Table 2.2: Execution times (in seconds). : Not available. : Memory overflow.
	We also give in Table 2.3 the memory usage of each package when loading the graph. Thanks
	to the CSR format, scikit-network has a minimal footprint.
		scikit-network NetworkX iGraph graph-tool
	RAM usage	1 222	17 765	10 366
	Table 2.3: Memory usage (in MB). : Memory overflow.

Table 3

 3

	code name	n	m	# classes
	SBM Stochastic Block Model 2.10 3 5.10 3	100
	CO	Cora	3.10 3 1.10 4	7
	CS	CiteSeer	3.10 3 1.10 4	6
	NG	20newsgroup	1.10 4 2.10 6	19
	WS	Wikipedia for Schools	5.10 3 2.10 5	14
	WV	Wikipedia Vitals	1.10 4 1.10 6	11

Wikipedia pages 1 , labeled by category (People, History, Geography...). .1: Some characteristics of the datasets.

Table 3

 3

	α	H	C	V	ARI	AMI	FMI	M	NSD
	0	0.19	0.27	0.22	0.0	0.01	0.03	0.45	0.76
	0.1	0.33	0.35	0.34	0.0	0.01	0.01	0.52	0.91
	1	0.36	0.37	0.36	0.0	0.01	0.01	0.50	0.92
	10	0.28	0.34	0.30	0.0	0.00	0.02	0.36	0.78
					NG
	α	H	C	V	ARI	AMI	FMI	M	NSD
	0	0.40	0.70	0.51	0.19	0.50	0.34	0.21	0.55
	0.1	0.44	0.70	0.54	0.22	0.54	0.35	0.21	0.59
	1	0.46	0.67	0.54	0.20	0.54	0.33	0.20	0.60
	10	0.37	0.55	0.45	0.13	0.44	0.26	0.17	0.56
					WS
	α	H	C	V	ARI	AMI	FMI	M	NSD
	0	0.23	0.29	0.25	0.05	0.25	0.26	0.25	0.49
	0.1	0.26	0.29	0.28	0.10	0.27	0.26	0.29	0.61
	1	0.23	0.24	0.23	0.04	0.23	0.20	0.30	0.65
	10	0.19	0.22	0.20	0.00	0.19	0.20	0.23	0.53

.3: Impact of regularization on clustering performance.

  .3. See how Louvain outperforms k-means by a significant margin.

	pipeline k-means (L2) k-means (cosine) Louvain
	V	0.69	0.74	0.93
		Table 4.3: V-measure scores.

Table 5

 5

	.1). Some of them have

Table 5

 5

	algorithm	CO	CS		NG	WS	WV
	GSVD + KNN 0.47 ± 0.08 0.30 ± 0.03 0.77 ± 0.02 0.33 ± 0.04 0.55 ± 0.02
	Propagation	0.31 ± 0.09 0.28 ± 0.05 0.07 ± 0.00 0.04 ± 0.00 0.03 ± 0.01
	Pagerank	0.54 ± 0.09 0.38 ± 0.04 0.62 ± 0.01 0.24 ± 0.03 0.50 ± 0.02
	GCN	0.15 ± 0.02 0.20 ± 0.03			0.12 ± 0.02 0.24 ± 0.06
	Heat Kernel	0.50 ± 0.09 0.37 ± 0.05 0.02 ± 0.01 0.19 ± 0.03 0.39 ± 0.04
	Dirichlet	0.53 ± 0.09 0.39 ± 0.04 0.60 ± 0.06 0.22 ± 0.04 0.51 ± 0.03
	Table 5.5: Macro F1-score (mean ± standard deviation). : Time Out.
		algorithm	CO	CS	NG WS WV
		GSVD + KNN 0.10 0.11 1.89 0.18 0.62
		Propagation	0.00 0.00 0.25 0.00 0.02
		Pagerank	0.41 0.40 0.96 0.40 0.69
		GCN	3.69 5.31		28.4 110
		Heat Kernel	0.12 0.10 1.34 0.17 0.41
		Dirichlet	0.04 0.05 1.54 0.23 0.74

.6: Running time for best score (seconds). : Time Out.

Table 6

 6 Table 6.1: Basic characteristics of the datasets. Bias indicates the proportion of the largest class.

	.1

3 Influence of the classifier for the spectral model

  ± 2.0 29.0 ± 1.2 70.0 ± 1.3 62.9 ± 1.0 DyF 86.3 ± 1.3 26.6 ± 1.2 73.1 ± 0.4 66.6 ± 0.Spectral 88.4 ± 7.0 43.7 ± 6.3 73.6 ± 3.5 75.2 ± 2.1 VRGC 86.3 ± 8.6 48.4 ± 6.2 74.8 ± 3.0 80.7 ± 2.2 Table 6.3: Mean accuracy (%) and standard deviation. Besides RFC, we experimented with different standard classifiers combined to our spectral embedding. Namely: k-nearest neighbors classifier (kNN), 2-layers perceptron with Relu nonlinearity (MLP), support vector machine with one versus one classification (SVM) and ridge regression classifier (RRC). Results are reported in Table 6.4. Table 6.4: Mean accuracy (%) of some classifiers combined to the spectral model.

			MT	EZ	PF	NCI1
	EMD	86.1 ± 0.8 36.8 ± 0.8	-	72.7 ± 0.2
	PM		85.6 ± 0.6 28.2 ± 0.4	-	69.7 ± 0.1
	FB		84.7 3
	SGE		87.3	40.7	71.9	-
	FGSD	92.1	-	73.4	79.8
		MT	PTC	EZ	PF	DD	NCI1
	RFC	88.4	62.8	43.7	73.6	75.4	75.2
	1NN	86.8	59.3	37.3	65.6	69.6	68.3
	15NN	85.7	61.9	33.7	70.4	75.0	69.6
	MLP	86.3	60.5	31.8	71.6	75.6	62.3
	SVM	85.3	60.8	31.3	73.0	75.0	63.9
	RRC	84.2	59.6	26.7	71.5	75.0	62.2

6.5.

  Table 6.5. Table 6.5: Mean accuracy (%) of the spectral model for different dimensions.

	k	MT	PTC	EZ	PF	DD	NCI1
	1	76.2	56.1	23.8	64.0	57.2	58.2
	5	86.8	62.5	39.0	69.6	73.9	72.5
	10	86.8	61.4	42.8	71.7	75.5	75.5
	25	88.4	62.8	42.7	72.8	75.7	75.2
	50	88.4	62.8	43.7	73.6	75.1	75.2

https://docs.scipy.org/doc/scipy/reference/sparse.html

https://graphs.telecom-paristech.fr/Home_page.html

https://scikit-network.readthedocs.io/en/latest/

https://github.com/sknetwork-team/scikit-network

https://scikit-network.readthedocs.io/en/latest/

https://netset.telecom-paris.fr/

A (k) 1 = d and A (k) 1 = f for all k > K.Proof. We have, 1. D -1 SU = UΣ 2 , where S is the adjacency matrix of the co-neighbor graph. Thus the square generalized singular values are eigenvalues of the stochastic matrix D -1 S. In particular, we have σ 1 = 1, and the multiplicity of the generalized singular value 1 is the number of connected components of the co-neighbor graph. This is also the number of connected components of the bipartite graph G, because two nodes are connected in G if and only if they are connected in the co-neighbor graph.2. If K = 1, then the generalized singular value σ 1 has multiplicity 1. It is then sufficient to verify that the vectors u 1 = 1/ √ w and v 1 = 1/ √ w satisfy Au 1 = Dv 1 , A v 1 = Fu 1 and u 1 Du 1 = v 1 Fv 1 = 1.

https://en.wikipedia.org/wiki/Wikipedia:Vital_articles/Level/4

https://netset.telecom-paris.fr/

https://snap.stanford.edu/data/

Chapter 4

Node clustering

This chapter is dedicated to node clustering. As mentioned in Chapter 1, this task consists in grouping nodes into a certain number of coherent clusters.

Formally, we denote by {C 1 , . . . , C K } a partition of the set of nodes V into K clusters. The associated cluster index in denoted by c i.e. if node i belongs to cluster C k , then c i = k while δ is the coclustering variable i.e., δ c i c j = 1 if i and j are in the same cluster, 0 otherwise.

Finally, Z ∈ R n×K is the membership matrix such that Z ik = 1(c i = k). In particular,

The rest of this chapter is organized as follows. Section 4.1 presents a family of algorithms based on the Louvain heuristic while Section 4.2 is dedicated to experiments, making the connection with Section 3.3.

Note that, even though we consider graph clustering from a knowledge discovery point of view here, it is also a common tool from graph coarsening in graph neural networks [START_REF] Simonovsky | Dynamic edge-conditioned filters in convolutional neural networks on graphs[END_REF]Komodakis, 2017, Deng et al., 2019]. Graph coarsening is the graph-equivalent operation of pooling in convolutional neural networks.

Louvain algorithms

Initially introduced in [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], Louvain algorithm is a popular clustering method for graphs based on the greedy maximization of an objective function called modularity. Its excellent speed-performance trade-off has made this algorithm the common denominator for community detection of many graph softwares such as NetworkX, Gephi, IGraph, Neo4J... Even-though this algorithm was originally designed for hard node clustering of undirected graphs, it can naturally be extended to handle bipartite, directed or signed graphs as well as

Proof. Starting from Equation (4.6), we derive arg min δ L = arg min δ KL (p, µ(δ)) -γKL (q, µ(δ)) ,

As log(µ ij ) = log( ω in ω out )δ c i c j + log(ω out ) and log( ω in ω out ) > 0:

Theorem 5 is complementary to the results of [START_REF] Newman | Equivalence between modularity optimization and maximum likelihood methods for community detection[END_REF] which connect modularity maximization and block models through statistical inference. Now, the following property is related to the so-called resolution limit of modularity optimization.

Proposition 10. The optimal partition at resolution γ for a graph with L connected components with respective total weights w 1 , . . . , w L is the concatenation of optimal partitions for each connected component at resolution γ l = w l w γ.

Proof. We denote by V 1 , . . . , V L the set of nodes corresponding to the connected components.

M |V l ,γ l denotes the modularity of the partition of V l at resolution γ l and M * |V l ,γ l the optimal one. This upper bound is reached by the concatenation of optimal partitions of V 1 , . . . , V l , hence the result.

Finally, note that finding the partition that maximizes the modularity is an NP-complete problem [START_REF] Brandes | Maximizing modularity is hard[END_REF]. Still, many algorithms offer to find local optimums. Some methods rely on spectral division [Newman, 2006b, White and[START_REF] White | A spectral clustering approach to finding communities in graphs[END_REF], some on simulated annealing [Guimera et al., 2004, Newman and[START_REF] Newman | Finding and evaluating community structure in networks[END_REF], but the most commonly used is Louvain greedy optimization which is the subject of Section 4.

Illustration

Figure 4.1 illustrates two partitions of the same graph connecting movies to actors. The first partition is obtained by maximization of the modularity, seeing the graph as an undirected one, while the second is obtained by maximization of the bimodularity. As we can see, these partitions do not even have the same number of clusters even though some of them look similar.

General benchmark

In this section, we compare the performance of several clustering algorithms on labeled datasets. For this experiment, we compute the Adjusted Rand Index with respect to the ground truth labels for the following algorithms:

Louvain algorithms: We refer to them as Louvain for the maximization of (4.3), DiLouvain for (4.5), BiLouvain for (4.1.1) and Potts for (4.4). We try different resolutions values for each one of them: γ ∈ {0.1, 0.2, . . . , 1.5}. We omit (4.1.1) as the pagerank vector without damping is usually not defined on real directed graphs.

Spectral algorithms:

We experiment with the standard Laplacian Eigenmaps [Belkin and Niyogi, 2002b] and with the algorithm described in Section 3.3 which we refer to as GSVD. has poor ARI on all graphs with the longest running times. Finally, the Propagation algorithm is very fast but yields very poor results.

In Chapter 5, we present a similar benchmark in a semi-supervised setting.

Digits classification

In order to further illustrate the interest of Louvain clustering, we perform an experiment on the MNIST dataset subset of Scikit-learn. The data is composed of 1 797 images in 8 × 8 gray scale pixels format. Each image represents a digit between 0 and 9. Some samples are displayed in Proof. We first prove that the maximum and the minimum of the vector T are achieved on the boundary S. Let i be any node such that T i is maximum. If i / ∈ S, it follows from (5.3) that T j is maximum for all neighbors j of i. If no such node belongs to S, we apply again this argument until we reach a node in S. Such a node exists because the graph is connected. It achieves the maximum of the vector T. The proof is similar for the minimum. Now consider two solutions T, T to Laplace's equation. Then ∆ = T -T is a solution of Laplace's equation with the boundary condition ∆ i = 0 for all i ∈ S. We deduce that ∆ i = 0 for all i (because both the maximum and the minimum are equal to 0), that is T = T.

Now let P S

ij be the probability that the random walk first hits the set S in node j when starting from node i. Observe that P S is a stochastic matrix, with P S ij = δ ij (Kronecker delta) for all i ∈ S. By first-step analysis, we have:

(5.4)

The following result provides a simple interpretation of the solution to the Dirichlet problem in terms of random walk in the graph: the temperature of any node is the average of the temperatures of the nodes at the boundary, weighted by the probabilities of reaching each of these nodes first:

Proposition 12. The solution to the Dirichlet problem is

(5.5)

Proof. The vector T defined by (5.5) satisfies for all i ∈ S:

where we have used (5.4). Thus, T satisfies (5.3). The proof then follows from Proposition 11.

Solution to the Dirichlet problem

We now characterize the solution to the Dirichlet problem in discrete space. Without any loss of generality, we assume that nodes with unknown temperatures (i.e., not in S) are indexed from 1 to ns so that the vector of temperatures can be written

where X is the unknown vector of temperatures at equilibrium, of dimension ns. Writing the transition matrix in block form as

so that:

(5.7)

Note that the inverse of the matrix I -Q exists whenever the graph is connected, which implies that the matrix Q is sub-stochastic with spectral radius strictly less than 1 [START_REF] Chung | Spectral graph theory[END_REF].

The exact solution of (5.7) requires to solve a (potentially large) linear system. In practice, a very good approximation is provided by a few iterations of (5.6), the rate of convergence depending on the spectral radius of the matrix Q. The small-world property of real graphs suggests that a few iterations are enough for real graphs [START_REF] Watts | Collective dynamics of small-world networks[END_REF]]. This will be confirmed by the experiments.

Heat kernel

Now consider the heat diffusion in the absence of boundary conditions:

First, note that, in this case, the average temperature of the nodes T is preserved over time:

The solution of this differential equation is given by:

The matrix H(t) = e -Lt is referred to as the heat kernel. It can be expressed through the spectral decomposition of L, L = UΛU with Λ = diag(λ 1 , . . . , λ n ) the diagonal matrix of eigenvalues λ 1 = 0 ≤ λ 2 ≤ . . . ≤ λ n and U = (u 1 , . . . , u n ) an orthogonal matrix of eigenvectors:

If the graph is connected, then λ 2 > 0 and

Since u 1 ∝ 1, the temperature becomes uniform at equilibrium. The interest of the diffusion lies in the transient states, for some finite time t.

Note that the exact computation of the heat kernel requires the spectral decomposition of the Laplacian matrix, which is not feasible for large graphs. In practice, an approximation of the system state at some finite time t can be derived directly from the heat diffusion (5.8). For some sufficiently small time step δ > 0, the state at time t follows from t/δ iterations of the following updates: T ← T -δLT.

(5.9)

Another approach consists in applying the diffusion in discrete time: Figure 5.2: Binary classification of the Karate Club graph [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF] with 2 seeds (indicated with a black circle). Red nodes have label 1, blue nodes have label 2.

Note that the temperature of each node can be used to assess the confidence in the classification: the closer the temperature to the mean, the lower the confidence. This is illustrated by Figure 5.2 (the lighter the color, the lower the confidence). In this case, only one node is misclassified and has indeed a temperature close to the mean, as illustrated in Example 5.3.1.

Example. 5.3.1: Binary classification

Multi-class classification

In the presence of more than 2 labels, we use a one-against-all strategy: the seeds of each label alternately serve as hot sources (temperature 1) while all the other seeds serve as cold sources Finally, we assess the classification performance of the algorithms in the case of seed asymmetry. Specifically, we first sample s = 1% of the nodes uniformly at random and progressively increase the number of seeds for the dominant class of each dataset, say label 1.

The process is repeated 10 times for each configuration. Figure 5.5 shows the macro-F1 scores.

We see that the performance of the centered algorithm remains steady in the presence of seed asymmetry.

General benchmark

We evaluate our Dirichlet algorithm against the following ones from the literature based on the macro F1-score. For each dataset, we generate 10 different seed sets by selecting 1% of the nodes uniformly at random. Each algorithm is tested using all the seed sets and the results are averaged.

We add some datasets from Chapter 4 and experiment with the following algorithms:

GSVD + KNN This two step method first performs the spectral embedding of the nodes of the graph [von Luxburg, 2007, Belkin andNiyogi, 2002c], then applies a k-nearest neighbors classification [START_REF] Omohundro | Five balltree construction algorithms[END_REF][START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF] in the embedding space, using the labels of the seeds. Here, we use the GSVD embedding projected onto the unit-sphere as described in Chapter 3. We test this algorithm with a dimension of the embedding space n dim ∈ {8, 16, 32} and the number of nearest neighbors k ∈ {1, 3, 5}.

Propagation [START_REF] Raghavan | Near linear time algorithm to detect community structures in large-scale networks[END_REF]. This method iterates over the nodes and assigns them the most represented label among their neighbors until convergence. It has the benefit of being hyper-parameter free and does not require the one-versus-all strategy in the case of a multi-class problem.

Pagerank [START_REF] Lin | Semi-supervised classification of network data using very few labels[END_REF]. This method is similar to the Heat Kernel except that the diffusion is replaced by a Personalized PageRank [Page et al., 1999b], with a restart distribution taken uniform over the hot sources. The pagerank is approximated using an iteration similar to (5.3) for a number of iterations N ∈ {5, 15, 25} and a damping factor α ∈ {0.75, 0.85, 0.95}.

Chapter 6

Graph classification

This chapter is dedicated to graph classification. Given a collection of graphs, the goal is to infer the class of unlabeled ones based on a labeled training set. Standard applications are molecular compounds classification in biology [START_REF] Childs | Identification and classification of ncrna molecules using graph properties[END_REF][START_REF] Kudo | An application of boosting to graph classification[END_REF] or malware detection in cyber-security [START_REF] Hu | Large-scale malware indexing using function-call graphs[END_REF][START_REF] Gascon | Structural detection of android malware using embedded call graphs[END_REF][START_REF] Kinable | Malware classification based on call graph clustering[END_REF].

The rest of this chapter is organized as follows, Section 6.1 introduces the specificity of graph classification with respect to supervised classification on standard vector data, Section 6.2 presents some related work, Section 6.3 presents a first algorithm published in [de [START_REF] Lara | A Simple Baseline Algorithm for Graph Classification[END_REF] while Section 6.4 presents a second algorithm published in [Pineau and de Lara, 2019a]. Finally, Section 6.5 covers the experiments.

• De [START_REF] Lara | A Simple Baseline Algorithm for Graph Classification[END_REF]. A simple baseline algorithm for graph classification.

Relational Representation Learning, NeurIPS 2018 Workshop.

• Pineau, E., & De Lara, N. (2019). Variational recurrent neural networks for graph classification.

In Representation Learning on Graphs and Manifolds Workshop.

Challenges

Many natural or synthetic systems have a natural graph representation where entities are described through their mutual connections: chemical compounds (see Figure 6.1), social or biological networks, for example. Therefore, automatic mining of such structures is useful in a variety of applications. However, graph classification raises two main difficulties to leverage standard machine learning algorithms.

First, most of these algorithms take vectors of fixed size as inputs. In the case of graphs, usual representations such as edge list or adjacency matrix do not match this constraint. The size of the representations is graph dependent (number of edges in the first case, number of nodes squared in the second).

Step All graphs represent chemical compounds, nodes are molecular substructures (typically atoms) and edges represent connections between these substructures (chemical bound or spatial proximity). In MT, the compounds are either mutagenic or not mutagenic. EZ contains tertiary structures of proteins from the 6 Enzyme Commission top level classes; it is the only multiclass dataset of this section. PF is a subset of the Dobson and Doig dataset representing secondary structures of proteins being either enzyme or not enzyme. In NCI1, compounds either have an anti-cancer activity or do not.

General benchmark

Each dataset is divided into 10 folds such that the class proportions are preserved in each fold for all datasets. These folds are then used for cross-validation i.e, one fold serves as the testing set while the other ones compose the training set. Results are averaged over all testing sets. We built the folds using Scikit-learn Pedregosa et al. [2011b] StratifiedKFold function with the random seed fixed to 1 in order to get reproducible results.

Spectral method

For the general benchmark, the embedding dimension is set to the average number of nodes for each dataset. Another experiment illustrates the influence of k, see Table 6.5. We use a unique set of hyper-parameters for the classifier is used for all datasets. We used the random forest classifier from Scikit-learn with class_weights: balanced. The other non-default hyper parameters were selected by randomized cross validation over the different datasets: 500 decision trees with 1 sample minimum per leaf, a maximum depth of 100 and use of bootstrap.

VRGC This model is implemented in Pytorch [START_REF] Paszke | [END_REF] and trained with the Adam stochastic optimization method [START_REF] Kingma | A method for stochastic optimization[END_REF] in the server setup 1.3.2. The input size d n of the recurrent neural network is chosen for each dataset according to the algorithm described in [START_REF] You | Graphrnn: A deep generative model for graphs[END_REF], namely 11 for MT, 25 for EZ, 80 for PF and 11 for NCI1. α is set to 0.1. For training, batch size is set to 64, and the learning rate to 10 -3 , decreased by 0.3 at iterations 400 and 1000. We use the same hyper-parameters for every dataset. Table 6.2 details the architecture of the network.

We compare our results to those obtained by Earth Mover's Distance [Nikolentzos et al., 2017b] (EMD), Pyramid Match [Nikolentzos et al., 2017b] In the mid 2000s, the arrival of social networks has amplified this phenomenon, creating new use-cases for these algorithms. Relationships between entities can be of multiple types: user-user symmetric relationships for Facebook or LinkedIn, follower-followee asymmetric ones for Twitter or even user-content bipartite ones for Netflix or Amazon. They all come with their own challenges and the applications are numerous: centrality calculus for influence measurement, node clustering for knowledge discovery, node classification for recommendation or embedding for link prediction, to name a few.

In the meantime, the context in which graph algorithms are applied has rapidly become more constrained. On the one hand, the increasing size of the datasets with millions of entities, and sometimes billions of relationships, bounds the asymptotic complexity of the algorithms for industrial applications. On the other hand, as these algorithms affect our daily lives, there is a growing demand for explanability and fairness in the domain of artificial intelligence in general. Graph mining is no exception. For example, the European Union has published a set of ethics guidelines for trustworthy AI. This calls for further analysis of the current models and even new ones. This thesis provides specific answers via a novel analysis of not only standard, but also extensions, variants, and original graph algorithms. Scalability is taken into account every step of the way. Following what the Scikit-learn project does for standard machine learning, we deem important to make these algorithms available to as many people as possible and participate in graph mining popularization