DPP Determinantal Point Process GD Gradient Descent GI Graph Interpolation GSP Graph Signal Processing GTR Graph Tikhonov Regularization LA Linear Algebra LASSO Least Absolute Shrinkage and Selection Operator LERW Loop-erased Random Walks LS Least Squares RLA Randomized Linear Algebra RSF Random Spanning Forest RST Random Spanning Tree SDD Symmetric Diagonally Dominant SOTA

Keywords:

An extensive range of problems in machine learning deals with data structured over networks/graphs. The examples vary from drug discovery to traffic forecasting, social network analysis to recommender systems, or epidemic analysis to natural language processing. Along with the exploding size and number of data, a big chunk of the proposed algorithms has focused on analyzing, representing and leveraging the network structures in the last decades. In many of them, the graph Laplacian is the central object. They are notable matrix representations of graphs that relate to various aspects. For example, by analyzing the Laplacian algebraically, one can measure the connectivity, count the number of spanning trees or generate a visualization of a graph. Further examples can be found in the problems of node clustering, graph sparsification, signal processing on graphs and many more. However, the algebraic analysis of Laplacian can be frustratingly time-consuming if the graph consists of a large number of nodes. Despite many numerical tools specialized for the graph Laplacian, in certain cases, the computational time remains one of the main issues.

Random spanning forests (RSFs), a random process on graphs, is a probabilistic tool for analyzing graphs with strong links with the Laplacian. In a nutshell, RSFs provide meaningful random sketches (spanning forests) of the graph to analyze some properties of interest. These sketches are cheap to obtain by a variant of Wilson's algorithm. Moreover, by only looking at a few of them, one can deduce significant statistical and/or algebraic information on the overall graph. The existing applications of RSFs, including graph wavelet analysis, semi-supervised learning on graphs or centrality analysis, show that RSFs provide useful information while binding diverse aspects of graphs.

In this manuscript, we present ways of leveraging the rich connections between RSFs and the graph Laplacian to speed up the algebraic analysis. In doing so, we propose randomized algorithms to approximate several quantities of interest, such as the regularized inverse of the Laplacian or effective resistances. Interestingly, these quantities are required in a wide range of applications on graphs, including graph signal filtering, hyper-parameter tuning, graph-based optimization and sparsification. We illustrate these methods on real-life networks and compare them with state-ofthe-art methods.

vi

Resumé

Un large éventail de problèmes d'apprentissage automatique traite des données structurées sur des réseaux/graphes. Les exemples vont de la découverte de médicaments aux prévisions de trafic, de l'analyse des réseaux sociaux aux systèmes de recommandation, ou de l'analyse des épidémies au traitement du langage naturel. Parallèlement à l'explosion des données, une grande partie des algorithmes proposés se sont concentrés sur l'analyse, la représentation et l'exploitation des structures de réseau au cours des dernières décennies. Dans beaucoup d'entre eux, les laplaciens de graphes sont les objets centraux. Ce sont des représentations matricielles notables des graphes qui se rapportent à divers aspects. Par exemple, en analysant les laplaciens de manière algébrique, on peut mesurer la connectivité, compter le nombre d'arbres spanning ou générer une visualisation d'un graphe. D'autres exemples peuvent être trouvés dans les problèmes de regroupement de noeuds, de sparsification de graphes, de traitement du signal sur les graphes et bien d'autres encore. Cependant, l'analyse algébrique du Laplacien peut prendre un temps frustrant si le graphe est constitué d'un grand nombre de noeuds. Malgré de nombreux outils numériques spécialisés pour les laplaciens, dans certains cas, le temps de calcul reste l'un des principaux problèmes.

Les Forêts Couvrante Aléatoire (FCA), un processus aléatoire sur les graphes, est un outil probabiliste pour analyser les graphes avec des liens forts avec les Laplaciens. En bref, les RSF fournissent des esquisses aléatoires du graphe pour analyser les propriétés d'intérêt. Ces esquisses sont faciles à obtenir par l'algorithme de Wilson. De plus, en ne regardant que quelques-uns d'entre elles, on peut déduire des informations statistiques et/ou algébriques significatives sur le graphe global. Les applications existantes des FCA, notamment l'analyse des ondelettes des graphes, l'apprentissage semi-supervisé sur les graphes ou l'analyse de la centralité, montrent que les FCA fournissent des informations utiles tout en liant divers aspects des graphes.

Dans ce manuscrit, nous présentons des moyens d'exploiter les riches connexions entre les FCA et les laplaciens des graphes pour accélérer l'analyse algébrique. Ce faisant, nous proposons des algorithmes aléatoires pour approcher plusieurs quantités d'intérêt, telles que l'inverse régularisé du Laplacien ou les résistances effectives. Il est intéressant de noter que ces quantités sont requises dans un large éventail d'applications sur les graphes, notamment le filtrage de signaux sur les graphes, l'ajustement d'hyperparamètres, l'optimisation basée sur les graphes et la sparsification. Nous illustrons ces méthodes sur des réseaux réels et les comparons aux méthodes de pointe.

Acknowledgement

Firstly, I would like to thank my supervisors Pierre-Olivier Amblard, Simon Barthelmé and Nicolas Tremblay. Their brilliant guidance and eternal patience are what made this thesis possible. Throughout three years, their door was always open for my instant inquiries and questions. I will always be grateful that they introduced me to the beautiful theory of random spanning forests and gave me the opportunity to develop things that I am proud of. Thanks to their pedagogical skills and efforts, I substantially improved my skills not only in scientific methodology but also in academic writing and communication.

My next acknowledgements go to the members of the jury, Samuel Vaiter, Agnès Desolneux, Olivier Michel, Nicolas Courty and Alexandre Gaudillière. I thank my reviewers Samuel Vaiter and Agnès Desolneux for their constructive comments and useful suggestions. By implementing their remarks, I improved a lot the clarity of some of the significant points that I made in the thesis. I would like to thank all the jury members for carefully examining my thesis and for their excellent questions and remarks during the defence.

Next, I would like to express my gratitude to Alexandre Gaudillière and Clothilde Mélot for our fruitful conversations and discussions, and I am really happy to have known them. Our discussions with Alex on random spanning forests (RSF) always led me to interesting research questions and ideas and Clothilde's talks inspired me about the applications of RSFs on signal processing. It is thanks to their hospitality, I have always enjoyed my trips to Marseille very much. I thank the members of the GAIA team for creating a vibrant and nurturing research environment. I appreciated a lot our weekly meetings and seminars, where I have the opportunity to listen to various lines of work of current research and broaden my horizons with different perspectives.

Then I would like to thank all of my lab-mates who have made GIPSA Lab an amazing workplace. For this, I thank a lot to Lorenzo, Malik, Cosme, Cyril, Sara, Gaëlle, Khan, Victor, Hugo, Hung, Clément, Yannis, Léa, Hélène, Sanjana, Matthieu, Martin, Olivier, Alexandre (Bleuzé), Alexandre (Constantin), Fatih, Davoud, Arthur and all of my office mates Lhoussein, Julien, Morgane and Juliette. Special thanks ix to Hugo (Jaquard) for all our discussions on RSFs and various other subjects which I always learnt a lot. Many times during these three years, I got by with a little help from my friends as the wise words of the Beatles say. In fact, "a little help" is an underestimation of their support and encouragement. I would like to take this opportunity to appreciate their help and company during my PhD. A huge thank you to Clémentine and Roza a lot for being such cool dudes, making my time in Grenoble always exciting and amusing, and for all of their support and help from the very beginning. I would like to thank Charles and Iris for all the fun times we shared in gatherings, concerts, travelling etc. and for all the interesting discussions. Thanks to Katia for her funny jokes and harp concerts. Many thanks to Barış Can (a.k.a. Cançam) for his fun company and useful advice over these years and for being an terrific partner in Bughouse chess (a chess variant). Special thanks to my friends from Turkey, Onat, Mustafa, Mert (Ak), Canberk, Onur, Barış, Sıla, Eylem, Kamil, Hakan, Atınç, Dogukan, Süleyman, Işınsu, Aybüke, Sercan, Emre, Mert, İrfan, Ugur, Özgür, Akcan, Bugra, Nazlı, Aslı and many others that I am not able to mention due to lack of space, for our continuous contact throughout the years. Although there have been long distances between us, I never felt lonely thanks to them, even during the isolation times of COVID-19. I thank my friends from Nomad Book Club, Nina, Iris, Caitong, Guillaume and all the past and present members for filling my Thursday evenings with exhilarating discussions. I would like to thank a lot to Elif Vural for introducing to me the exciting field of graph signal processing, for all the interesting discussions, for our collaboration from which I learnt tremendously, and for her support and precious advice throughout the years.

Finally, I would like to express my deepest gratitude to my family. I cannot thank enough my parents Oya and İhsan for being excellent role models and always supporting me in the path I wanted to follow, my grandmother Nurten for being such a great inspiration, and my sister İpek for being a great friend that I always enjoy spending time with. I love them dearly and this thesis would not be possible without their constant support and encouragement.

Grenoble, 2023

x Introduction 1

" "Graphs are ubiquitous structures."

-Almost any paper on graph ML Since they first appeared for modeling the bridges of Königsberg, graphs, i.e.a set of vertices and edges, have been used to model and solve numerous problems in physics, chemistry, biology, computer science, social sciences and many more. They are so ubiquitous because they are natural models to use whenever the data represent mutual relations between individual elements. Today, processing data structured on graphs is one of the major research questions in signal processing and machine learning. Voluminous research focuses on leveraging the graph structure to accomplish advanced tasks, such as managing/analyzing Internet (a colossal computer network), understanding the human brain [START_REF] Farzad V Farahani | Application of graph theory for identifying connectivity patterns in human brain networks: a systematic review[END_REF], drug discovery [START_REF] Jiang | Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models[END_REF], analyzing social networks [START_REF] Otte | Social network analysis: a powerful strategy, also for the information sciences[END_REF] or large-scale traffic forecasting/planning [START_REF] Hasanzadeh | A graph signal processing approach for real-time traffic prediction in transportation networks[END_REF]. Further examples appear in natural language processing [START_REF] Nastase | A survey of graphs in natural language processing[END_REF], finance [START_REF] Vinıécius De | Learning undirected graphs in financial markets[END_REF], food networks [START_REF] Dunne | Network structure and robustness of marine food webs[END_REF], power grid networks [START_REF] Giuliano | The power grid as a complex network: a survey[END_REF] and many more. In many of these examples and others, graphs come with some features/signals, also called graph signals, on the vertex set. Here are some examples:

• Transportation networks [START_REF] Hasanzadeh | Piecewise stationary modeling of random processes over graphs with an application to traffic prediction[END_REF]: A prominent example is road network.

In a typical setup, each location corresponds to a node, and two nodes are connected by an edge whenever there is a road between them. The signals over such networks are usually some traffic features, such as the traffic speed or volume.

• Social networks [START_REF] Villalba-Diez | Geometric Deep Lean Learning: Evaluation Using a Twitter Social Network[END_REF]: In many social studies, along with the individual information over the subjects, the structure of their interactions is also of central importance. Social media networks such as Facebook and Twitter are probably the most-known examples. In those cases, every user is considered a node and two nodes are connected if there is an interaction between the corresponding users, such as friends on Facebook or followers on Twitter. Given this graph/network structure, any feature of a user, such as their likes/dislikes, can be analyzed as a graph signal. • Temperature networks [START_REF] Perraudin | Stationary signal processing on graphs[END_REF]: A temperature network is a network of temperature sensors over a region in which each sensor is a node. Here the edges can be established according to different criteria. The common practice is to build the nearest neighbour graph from the proximity of the sensors. The signal, in this case, is naturally the temperature measurements collected from the sensors.

• Neurological networks [START_REF] Huang | A graph signal processing view on functional brain imaging[END_REF]: Given numerous unanswered questions, the analysis of the brain is indeed an exciting and arduous task. The graphbased study, in fact, gives a natural way to analyze brain signals. In this type of analysis, we often consider the different regions of the brain as the nodes and the edges are established according to their structural or functional properties such as correlation in their activities. Analyzing neurological signals over such a network can decipher interesting information about the brain.

• Other biological networks [GDJSRLHVRT+21; CPFSC21; YB20]: The brain is not only the biological subject of graph-signal analysis. In the vast literature, we find many studies related to drug discovery, analysis of protein structures or gene interactions.

Given the large volume and diversity of such datasets, developing appropriate tools for processing and analyzing them has become quite central. In fact, research in machine learning and signal processing for graphs and graph signals has been emerging in the last decades [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Ricaud | Fourier could be a data scientist: From graph Fourier transform to signal processing on graphs[END_REF]. In turn, these tools have been used for solving real-life problems varying from fake-news detection in social media [START_REF] Monti | Fake news detection on social media using geometric deep learning[END_REF] to spread analysis for COVID-19 [START_REF] Panagopoulos | Transfer graph neural networks for pandemic forecasting[END_REF], from decoding brain signals [START_REF] Ouahidi | Pruning Graph Convolutional Networks to select meaningful graph frequencies for fMRI decoding[END_REF] to weather forecasting [START_REF] Keisler | Forecasting Global Weather with Graph Neural Networks[END_REF], and so on. In most of these tools and methods, the matrix representations of graphs, especially the graph Laplacian, are of significant importance.

Graph Laplacian. Many analysis tools extract information from the matrix represen-

A Path Graph Its Laplacian L   1 -1 0 -1 2 -1 0 -1 1  
tations of the graphs, particularly from the graph Laplacian. The graph Laplacian is a symmetric square matrix in the size of number of the vertices. It can be considered as the discrete Laplacian operator on graphs when applied to a function over the vertices:

(Lf ) i = j∈N (i) (f (i) -f (j)),
where N (i) denotes the nodes that share an edge with node i (formal definition can be found in Section 2.1.1). The graph Laplacian is a simple but useful matrix representation of graphs. The algebraic analysis relates them with many useful properties of graphs. For example, many tasks in connectivity analysis of the graph can be accomplished by observing the eigendecomposition of the graph Laplacian. Simply by looking at the multiplicity of zero eigenvalues, one can count the number of connected components i.e. disconnected sub-parts of the graph. Moreover, the second smallest eigenvalue, also called spectral gap, gives a measure of the connectivity of the graph, which is mainly used in applications of stability and robustness of dynamic networks. Similarly, the corresponding eigenvector, also known as the Fiedler vector, yields the sparsest cut partition of the graph i.e. the partition that leaves the least number of edges inter-parts. Spectral graph drawing is another problem that uses the eigendecomposition of graph Laplacian. In this problem, one seeks a drawing (a mapping from vertices to Euclidean coordinates) of a graph that minimizes the distance between the vertices that are close to each other in the graph. It turns out that certain eigenvectors of graph Laplacian is the analytic solutions to this problem. Closer to the main themes of this thesis, we find an interesting link between the eigenvalues of L to the enumeration of a certain type of sub-graphs. In particular, the celebrated Kirchoff's matrix-tree theorem states A spanning tree. Its edges are given in red, and the dashed lines indicate the edges of the graph.

that the product of the non-zero eigenvalues of L equals the number of all spanning trees, i.e. connected sub-graphs that contains no cycles, on a graph.

Apart from these examples, the graph Laplacian has an important role in analyzing electrical resistor networks. In these networks, we model each resistor by an edge where the edge weight depicts the conductance of the resistor, and each node corresponds to its connection points. These models are extensively employed in circuit theory in order to calculate the currents flowing through the resistors/edges or the voltages induced at the nodes when a fixed potential difference (voltage) is applied between two (or more) nodes. Surprisingly, these calculations boil down to solving linear systems involving the graph Laplacian (See Chapter 4).

In the last decades, even more, examples have emerged in signal processing and machine learning involving graphs.

Laplacian in graph signal processing. Graph Signal Processing (GSP) is the subfield of signal processing that deals with the signals defined over the vertices of the graph. In the last decade, many of the classical signal processing tools have been adapted in order to deal with such signals. Examples include filtering [SVF11; SNFOV13], sampling [TAB17; PTGV18], translation operation [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF], wavelet analysis [HVG11; NO12; ACGM20] or uncertainty principle [START_REF] Tsitsvero | Signals on graphs: Uncertainty principle and sampling[END_REF]. In many of these adaptations, the graph Fourier transform plays a significant role. By analogy, it allows to represent the graph signals in the graph frequency domain. In turn, one can define graph translation and filtering schemes by using an analogous of the celebrated convolution theorem [START_REF] Oppenheim | Discrete-time signal processing[END_REF]. However, due to the irregularity of the domain of the signal, the definition of the Fourier transform is not straight-forward. One popular definition in GSP builds an analogy on the fact that the eigenfunctions of the Laplacian operator are the basis functions in the Fourier transform. Developing on this, the authors in [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] suggest using the eigenvectors of the graph Laplacian as the Fourier basis and the eigenvalues as the graph frequencies. In this way, higher eigenvalues correspond to a higher frequency component in the Fourier analysis. By projecting the signal on these basis, one can calculate the graph frequency response of the signal to each component.

Laplacian in machine learning. An early use of the graph Laplacian occurs in the semi-supervised learning on graphs [START_REF] Zhu | Semi-supervised learning with graphs[END_REF][START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. In this problem, we are given a few labels over the vertices and the goal is to infer the labels of the others by using the given labels and the underlying graph. A baseline solution given by [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] suggest using the graph Laplacian to formulate this problem. Their main assumption is that the labeling function is fixed at a subset V ⊂ V and smooth over the rest U = V \ V , i.e.it does not vary too much through the edges. For K classes of labels, this formulation seeks a classification function F : V × {1, . . . , K} → R such that for each node i, argmax k F (i, k) yields the label of node i. Let us denote the known labels by Y ∈ R |V |×K , then, [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] suggest solving the following constrained problem:

F = argmin F i∈U j∼i (F (i) -F (j)) 2 s. t. ∀i ∈ V, k ∈ {1, . . . , K}, F (i, k) = Y i,k ,
In fact, the minimized term can be written in terms of the graph Laplacian as FL U F where F = [F (i, k)] i∈U,k∈{1,...,K} . Moreover, the closed form solution involves the following in the graph Laplacian:

F = (L U ) -1 L U |V Y.
Later on, this formulation is considered in an unconstrained setup and generalized in [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. In these forms as well, the solution takes the form of the matrix inverse of the graph Laplacian. More recent examples in machine learning that extensively use graph Laplacian include spectral clustering [Von07; TPGV16], graph embedding [START_REF] Xu | Understanding graph embedding methods and their applications[END_REF], sparsification [START_REF] Daniel | Graph sparsification by effective resistances[END_REF] and deep learning on graphs [START_REF] Wu | A comprehensive survey on graph neural networks[END_REF]. In all of these applications and many more, the analysis of L is of utmost importance.

Laplacian-based Numerical Algebra. As listed above, the number of applications involving the graph Laplacian is myriad. In many of these applications, the solution requires computing an (regularized or pseudo) inverse or eigendecomposition of the graph Laplacian. However, computing these quantities directly does not scale well with the size of the graph and even becomes impractical for very large graphs. In turn, many researchers have focused on developing numerical algebra tools specified for Laplacian. The most prominent ones are closely related to spectral graph theory [START_REF] Spielman | Spectral graph theory[END_REF] and solving Laplacian linear systems [START_REF] Nisheeth | Lx= b[END_REF]. These studies give a rich collection of algebraic tools that avoid the expensive direct computation and, instead, approximate the required solution within certain error ranges. In doing so, they leverage some common properties of the graph Laplacian, such as its algebraic properties, its strong links with graph theory, or the sparsity of real-life graphs. Compared to these works, we take in this thesis an alternative path by using randomization in order to develop Laplacian-specific algebraic tools.

Randomized Linear Algebra (RLA). In a nutshell, RLA is a set of tools that avoid the direct computation of an algebraic operation via randomization. The main advantage of these tools is that they can approximate very expensive operations such as matrix inversion, and singular value decomposition, very cheaply by only taking random samples from the input matrix. More details and examples can be found in [START_REF] Drineas | RandNLA: randomized numerical linear algebra[END_REF]. Closer to the themes of this thesis, [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] show that determinantal point processes (DPP), i.e. a random point process with many tractable properties, can be easily adapted as RLA algorithms in order to solve linear systems involving a large set of input matrices covering the graph Laplacian. Inspired by these examples and many more, we focus on a particular determinantal point process that allows us to develop randomized tools for Laplacian-based numerical algebra.

Random Spanning Forests (RSF).

A forest in a graph is a subgraph which contains no cycle in it. It is called spanning if it contains all the vertices of the graph. Finally,

A spanning forest with three trees. The edges of the forest are in red, and the graph is depicted by dashed lines.

the random spanning forests are the random processes where we choose spanning forest at random on a given graph. One can come up with infinitely many options for the distribution of such a process. In this thesis, we focus on the distribution introduced in [START_REF] Avena | Random forests and networks analysis[END_REF]. A particular case of this distribution coincides with the uniform spanning trees which is the process of uniformly selecting a spanning tree. They are well-known objects in probability theory, computer science and statistical physics [START_REF] Grimmett | The random-cluster model[END_REF]. In this distribution, the probability of having a fixed spanning forest is mainly regulated by the weights of its edges and the number of its connected components (trees). We restrict ourselves to this distribution due to two main reasons:

• The RSFs sampled from this distribution have rich theoretical connections with the graph Laplacian via certain DPPs,

• There exists an efficient algorithm to sample from this distribution, called Wilson's algorithm.

Thanks to these facts, we know that the samples of RSFs can be obtained cheaply, and their algebraic properties closely relate to the algebraic properties of the original graph. The main goal of this thesis is to leverage these facts in order to develop RLA algorithms for Laplacian-based numerical algebra.

The main contributions of this thesis can be listed as follows:

• We shed light on the rich connections of RSFs with determinantal point processes and the graph Laplacian; We go through the elegant theory and useful properties of RSFs as DPPs. Throughout our tour, while for some of the properties, we reproduce the existing proofs, for some others, we give mostly algebraic proofs that are more accessible for the researchers from different domains (See Theorems 2.4.4, 2.4.5).

• In turn, we develop and analyze RLA algorithms based on RSFs for approximating the solutions of a wide variety of problems involving the graph Laplacian. These problems are:

-Tikhonov regularization and interpolation for graph signals: The problem of denoising, i.e. removing noisy parts of a signal, and completing missing parts of a signal are long-standing problems in signal processing. We first give a unifying optimization framework for these problems for graph signals i.e. signals defined over vertices. Then, we give novel randomized algorithms that approximate the solution for this framework.

-Estimating the trace of regularized inverse of symmetric diagonally dominant (SDD) matrices: Trace (sum of diagonal entries of a matrix) is a central operation in linear algebra. However, computing the trace of a matrix is a challenging problem whenever the input matrix is not directly accessible e.g. multiplication of large matrices or inverse of a large matrix. We give novel algorithms based on RSFs for estimating the regularized inverse of SDD matrices (a class of matrices that contain the graph Laplacian). In turn, our algorithms are at least comparable, and they usually outperform state-of-the-art methods.

-Estimating effective resistances in electrical networks: Effective resistance is a metric that takes its root from electrical resistor networks 1 and they are used for measuring the similarity of a pair of nodes. They play a central role in many graph-related applications such as graph sparsification, clustering or graph learning. However, the computation of this useful quantity does not scale well with increasing sizes of graphs as one needs to invert a graph Laplacian. To avoid this expensive computation, we give RSF-based algorithms to estimate effective resistances. The proposed methods are comparable and even better than the state-of-the-art algorithms when the number of desired effective resistances is small.

-Estimating various graph filters: This aspect of our work makes interesting links between RSFs and graph filtering i.e. spectral signal filtering for graph signals [START_REF] Tremblay | Design of graph filters and filterbanks[END_REF]. Graph filtering is an essential tool for dealing with these type of signals and have been used in various applications. Yet again, the associated computations often require diagonalizing L, which is expensive. In this work, we give a set of graph filters that can be approximated via RSF-based estimators by avoiding the diagonalization of L. These filters either cover or can mimic by optimizing a few hyperparameters some of the frequently used filters, such as the ideal low-pass filter.

• We provide a bias-variance analysis of the proposed algorithms for stating their theoretical performance.

• We illustrate these algorithms in real-life applications and datasets while comparing these algorithms with the existing algorithms.

The main structure of the thesis is as follows; we start by giving necessary technical machinery in Chapter 2 on graph theory, randomized linear algebra and random spanning forests. Then, by Chapter 3, we introduce the RSF-based methods for solving graph Tikhonov regularization. We compare these methods with existing methods and show their uses in other graph-related problems while illustrating real-life applications. In Chapter 4, we extend the range of RSFs-based methods into a new set of problems, including estimating the trace of regularized inverse of SDD matrices, effective resistances and certain graph filters. Finally, we finish in Chapter 5 with a general conclusion, discussion on open questions and future works.

The peer-reviewed publications associated with this thesis can be found in the following: Also, some materials are introduced in my master thesis [START_REF] Yusuf Yigit Pilavci | Random spanning forests: theory and applications[END_REF].

Background 2 " Mathematics compares the most diverse phenomena and discovers the secret analogies that unite them.

-Joseph Fourier

In this chapter, we go through the necessary technical definitions and tools to be able to explain the main contributions of this thesis. We collect these diverse materials under four sections:

• Graph theory contains the basic definition of a graph, graph matrices (e.g. adjacency, Laplacian etc.), graph-related structures (e.g. walks, paths, cycles, trees, forests etc.), random processes on graphs.

• Determinantal Point Process (DPP) contains the basic definitions on DPP and their useful properties.

• Randomized Linear Algebra (RLA) presents the use of randomization for numerical linear algebra with emphasis on least-squares problems.

• Random Spanning Forest (RSF) introduces a random object on graphs, called random spanning forests, which has fascinating probabilistic properties.

As they are listed above, these concepts might seem unconnected. However they have elegant mathematical links that tie them to each other. One of the main goals of this chapter is to bring these links to spotlight.

Graph Theory

We dedicate this section to defining all graph-related objects, which will be repeatedly used in the rest.

Definition 2.1.1 (Graph). A graph G = (V, E, w) is a triple which consists of n nodes/vertices V = {v 1 , v 2 , . . . , v n }, m edges connecting the vertices E ⊆ V × V
and the weight function w : V × V → R ≥0 which maps every edge (i, j) ∈ E to a positive real value, and w(i, j) = 0 for any pair (i, j) ∈ E.

In many problems of physics, mathematics, biology, computer science and machine learning, the data are known up to pairwise relations of data points. In such cases, graphs/networks are natural choices to interpret, visualize or analyze the data. In doing so, certain distinctions over graphs facilitate the understanding of real-life networks. Let us list them in the following paragraphs:

A directed graph is a graph in which every edge (i, j) has an orientation from i to j and w(i, j) is not necessarily equal to w(j, i). If all edges in G verifies w(i, j) = w(j, i), then we call it an undirected graph as the orientation of the edges are no longer informative about the graph1 . The weight function w of a graph G maps every edge to either 1 or 0, i.e. ∀(i, j) ∈ E, w(i, j) ∈ {0, 1}, then G is called unweighted. We denote the neighbors of a node i as N (i) := {j : j ∈ V and (i, j) or (j, i) ∈ E}, i.e. the set of nodes which has an incident edge (j, i) or (i, j) in G. The sum Then, we pass its degree distribution given in (a) to the Kalisky model [START_REF] Kalisky | Tomography and stability of complex networks[END_REF] (a statistical model for generating a graph given a degree distribution) to generate another graph (c) that has the same degree distribution.

d i := j∈N (i) w(i, j) is called the degree of node i.
The degrees d 1 , d 2 , . . . , d n are limited but useful devices for characterizing graphs. They are often analyzed by statistical tools. In particular, the degree sequence d 1 , d 2 , . . . , d n of a graph is considered as a sample of a probabilistic model, also called a degree distribution. In this way, one can reduce large and complex real-life networks to some degree distributions with a few hyper-parameters. However, when it comes to detecting finer details of the graph, the degree distribution might fail.

We illustrate such a case with a toy example in Fig. 2.1. In a nutshell, we show that we can obtain two graphs with completely different characteristics whereas they have exactly the same degree sequence. In this case, it makes sense to take our focus on some of the substructures of a graph that contain descriptive patterns, such as connected components, paths, cycles or trees. We finish this section by giving the formal definitions of some of such substructures.

Definition 2.1.2 (Subgraph). A graph H = (V s , E s , w s ) is a subgraph of a graph G, if it verifies V s ⊆ V, E s ⊆ E and w s (i, j) = w(i, j) for all (i, j) pairs.

Graph Theory

Definition 2.1.8 (Spanning Tree). A tree, denoted by τ = (V τ , E τ , w), is a subgraph which contains no cycle in it. If V τ = V, then τ is called a spanning tree.

Taken from shorturl.at/cfkpu Definition 2.1.9 (Rooted Spanning Tree). A rooted tree is a directed tree in which every edge is oriented towards a special node called the root.

Definition 2.1.10 (Rooted Spanning Forest). A rooted forest, denoted by φ =

(V φ , E φ , w), in G is a set of disjoint rooted trees with V φ ⊆ V, E φ ⊆ E. Similarly, if V φ = V, φ is called a rooted spanning forest.
We illustrate some of these graph structures in Fig. 2.2

Graph Matrices

Matrix representations of graphs are helpful tools to examine probabilistic, algorithmic and combinatorial aspects of graphs. In this section, we will revisit the essential matrix representations of graphs.

The simplest matrix representation of a graph is the (weighted) adjacency matrix, defined as W := [w(i, j)] i,j ∈ R n×n . Another simple matrix, called the degree matrix, contains the degrees in a diagonal matrix D where ∀i ∈ V, D i,i := d i . Finally, we define the combinatorial graph Laplacian as L := D -W. One can find various ways to define the Laplacian operator for graphs, such as the normalized Laplacian I -D -σ WD σ-1 with σ ∈ {0, 1}. In this thesis, we refer to L as the graph Laplacian unless otherwise stated. We stress that there are various theoretical properties of L which connect diverse concepts such as determinantal point process [START_REF] Burton | Local characteristics, entropy and limit theorems for spanning trees and domino tilings via transfer-impedances[END_REF], combinatorics [START_REF] Chaiken | Matrix tree theorems[END_REF], computer science and graph signal processing [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. This thesis aims to give a unifying view of these concepts and leverage the connections for randomized linear algebra applications. Thus L is an object of central importance.

The edge incidence matrix B ∈ R m×n is a less common but still useful representation of graphs. It depicts the orientation of the edges as follows:

∀i, j, B k,i =          -w(i, j), e k = (i, j) ∈ E w(j, i), e k = (j, i) ∈ E

0, otherwise

In case of undirected graphs, we retain B in this form by encoding every undirected edge e k = (i, j) as B k,i = -w(i, j) and B k,j = w(i, j), for all i < j without loss of generality. The following identity relates B with the graph Laplacian L:

L = B B. (2.1)
This relation also proves the positive semi-definiteness of L.

Spectral properties, i.e. eigen/singular value decomposition, of these matrices, are studied within many fields, from signal processing [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] to clustering [START_REF] Von | A tutorial on spectral clustering[END_REF].

Similarly, they are of central importance in this thesis. Therefore, we revisit some of these properties and their uses in the rest.

We generically write the spectral decomposition of W as W = PΞR where P =

Spectrum of W [p 1 | . . . |p n ] and R = [r 1 | . . . |r n ]
contain the left and right eigenvectors p i ∈ C n and r i ∈ C n , and Ξ is the diagonal matrix containing the eigenvalues

ξ 1 ≤ ξ 2 ≤ . . . ≤ ξ n .
For an undirected graph, the adjacency matrix W is entry-wise nonnegative and symmetric, thus P = R. The spectrum of W is insightful for understanding the structure of the graph. For example, the multiplicity of the largest eigenvalue gives the number of connected components in the graph. Another example is that if the eigenvalues are symmetric with respect to 0, i.e. ∀i ∈ 1, . . . , n, ξ i = -ξ n-i , then the graph is bipartite.

As W is an entry-wise non-negative matrix, it is possible to examine its spectrum with the celebrated Perron-Frobenius theorem [Per07; FFFFM12]. Along with other properties, the Perron-Frobenius theorem associates the largest eigenvalue with a positive eigenvector. An influential application of this property in directed graphs takes place in the famous algorithm of Google, called PageRank [PBMW99] for sorting web pages by their relevance given their hyperlink graph i.e. whenever a page links to another, there is a directed edge between them. PageRank approximates different websites according to their ranking score, which is based on the principal eigenvector of the graph of the world wide web (WWW) (See Example 2.1.1).

A $25 Billion Algorithm a

a The approximate market value of Google in 2004 [START_REF] Bryan | The $25,000,000,000 eigenvector: The linear algebra behind Google[END_REF] Information retrieval i.e. retrieving the relevant information in a humongous amount of data is a long-standing problem in computer science. In the last two decades, research in this field has provided many efficient tools with the help of machine learning and computer science. Google's web search engine has been probably the most popular one. In the early stages of search engines, one main issue was efficiently ordering different web pages in search results according to their relevancy. In [START_REF] Page | The PageRank citation ranking: Bringing order to the web[END_REF], the founders of Google proposed the celebrated PageRank algorithm to rank web pages, which relies on probability and algebraic graph theory.

Fig. 2.3.:

Given a graph (in the left), PageRank is an algorithm to rank nodes according to their importance in the graph. In the right, we plot the graph by resizing the nodes according their PageRank scores. Bigger nodes have higher scores.

Consider the hyperlink graph of the world wide web (WWW) in which each website is a node, and there is a directed edge whenever there is a hyperlink from one to another. Given such a graph, one seeks to assign a ranking score r i to each node/web page i. The PageRank algorithm proposes to use the ranking scores r i 's which verify:

∀i ∈ V, r i = 1 d i j∈N (i) w(i, j)r j ≥ 0.
Intuitively, the ranking scores of its neighbours define the ranking score of a node. The web pages pointed by the pages with high scores have higher scores.

Besides, this formulation corresponds to solving the following eigenvector problem: r = D -1 Wr.

By the Perron-Frobenius theorem, we know the largest eigenvalue is 1 a and it is associated with an entry-wise positive eigenvector r. PageRank draws a simple iterative scheme, also called the Power method, to estimate r as follows:

r (k+1) = D -1 W r (k) ,
starting from arbitrary initialization r 0 . As k goes to infinity, r k converges to r. An illustration of this algorithm is in Fig. 2.3.

a This is due to the bounds on the spectral radius given by the Perron Frobenius theorem i.e. mini

n j=1 (D -1 W)i,j = 1 ≤ ρ(D -1 W) ≤ maxi n j=1 (D -1 W)i,j = 1 16 Chapter 2 Background
Similar to the adjacency matrix, the spectrum of the (combinatorial) graph Laplacian L also provides many interesting insights into the graph. We detail some of them for undirected graphs as they are of central importance for the main concepts of this thesis. Let us denote the eigen-decomposition by L = UΛU where U

= [u 1 | . . . |u n ] ∈
R n×n contains the eigenvectors u i 's in the columns and Λ is the diagonal matrix which consists of the eigenvalues λ 1 = 0 ≤ λ 2 ≤ . . . ≤ λ n in its diagonal entries. Note that L is a positive semi-definite matrix as L = B B, and its smallest eigenvalue is 0. The multiplicity of the smallest eigenvalue at 0 equals the number of connected components. Moreover, the second smallest eigenvalue λ 2 , also called algebraic connectivity or Fiedler value, is considered a measure of the connectivity of the graph. For example, if the graph has two connected components, then one also has λ 2 = 0, and it increases as one adds edges connecting the two components. Another significant property of L is due to the matrix-tree theorem [Kir47]:

Theorem 2.1.1 (Matrix-Tree Theorem). Given a graph Laplacian L of a connected, unweighted and undirected graph, let T be the set of all spanning trees. Then, one has:

∀i ∈ V, |T | = det L -i|-i
where L -i|-i is the submatrix of L obtained by deleting i-th row and column.

Proof. Among many possible proofs of the matrix-tree theorem, we prefer to provide an algebraic one which will be re-invoked in the further chapters. The Cauchy-Binet determinant formula yields:

det L -i|-i = S⊆E |S|=n-1 det B S|-i det(B ) S|-i , = S⊆E |S|=n-1 (det B S|-i ) 2 .
Now we invoke a technical result from [START_REF] Van Mieghem | Graph spectra for complex networks[END_REF] which is attributed to Poincaré. To do so, we consider the square submatrices of the edge incidence matrix B in size of (n -1) × (n -1). Notice that the rows and columns of such submatrices are respectively indexed by n -1 edges and nodes of the graph. Let us denote these edges and nodes by S ⊂ E and -i = V \ i. Theorem 2 in [START_REF] Van Mieghem | Graph spectra for complex networks[END_REF] states 2.1 Graph Theory that det B S|-i equals to ±1 if and only if S is a spanning tree of G. Otherwise it equals to 0. Plugging this result in the Cauchy-Binet formula finishes the proof:

det L -i|-i = S⊆E |S|=n-1 (det B S|-i ) 2 = S⊆E |S|=n-1 I(S is a spanning tree)
where I is the indicator function.

At last but not least, we detail the spectral decomposition of B. Let us denote its singular value decomposition by

B = VΣU where V = [v 1 | . . . |v n ] ∈ R m×m and U = [u 1 | . . . |u n ] ∈ R n×n
are the orthonormal basis and Σ ∈ R m×n is a rectangular diagonal matrix which contains the singular values σ 1 ≤ . . . ≤ σ n in its diagonal entries. Note that U is also the orthonormal basis of L and ∀i ∈ V, one has λ i = σ 2 i as an immediate result of the identity in (2.1). As mentioned in the proof of Theorem 2.1.1, Poincaré's theorem [START_REF] Van Mieghem | Graph spectra for complex networks[END_REF] makes interesting connections between the edge incidence matrix B and the spanning trees of the corresponding graph. These links are of central importance for this thesis. Thus, to conclude our visit on algebraic graph theory, we reproduce and extend this theorem for weighted graphs and spanning forests. where S r is the edge set of the tree that contains node r ∈ R and S = r∈R S r . Otherwise, this absolute determinant is equal to 0.

Proof. As long as S forms a spanning forest, it includes |R| disconnected trees due to the fixed number of edges in trees (i.e.a tree with t nodes consists of exactly t -1 edges). Similar to Theorem 2.1.2, the rest of the proof is also two-fold: In the first part, we consider the case that the given conditions hold.

In the second, otherwise is examined to finish the proof.

Assume that S forms a spanning forest with |R| trees and it does not include any path between the nodes in R = {r 1 , . . . , r |R| }. Then, each tree in such forest contains exactly one distinct node r ∈ R. This allows us to write a row-column permutation of square matrix B S|-R in the following form:

perm(B S|-R ) =          B (r 1 ) Sr 1 |-r 1 B (r 2 ) Sr 2 |-r 2 . . . B (r |R| ) Sr |R| |-r |R|         
where perm(B S|-R ) is a particular row-column permutation of B S|-R in which each tree's reduced edge incidence matrix B

(r)

Sr|-r appears in the diagonals. The absolute determinant of this block diagonal matrix writes:

| det perm(B S|-R ) | = r∈R | det B (r) Sr|-r | where | det B (r) Sr|-r | is equal to a,b∈Sr w(a, b) 1/2 due to Prop. 2.1.2. Noticing | det perm(B S|-R ) | = | det B S|-R | writes: | det B S|-R | = r∈R   a,b∈Sr w(a, b)   1/2 =   a,b∈S w(a, b)   1/2
and finishes the first part of the proof. The contrary of the constraints given in the theorem yields either S includes a cycle or there exists a path between the nodes in R. The reduced edge incidence matrix for a graph with this path writes:

B {S|-R} =              w(r, 1) 0 . . . 0 -w(1, 2) w(1, 2) 0 . . . 0 0 -w(2, 3) w(2, 3) 0 . . . 0 . . . . . . 0 . . . 0 -w(k, p) 0 X Y             
Regardless of matrices X and Y (and so, the rest of the graph), we can conclude that this matrix is singular (i.e. det B S|-R = 0) because it has a non-zero right null vector x T = 1 w(r,1) w(1,2) w(r,1) w(2,3) . . . w(r,1) w(k,p) 0 which gives x T B S|-R = 0 T . Notice that a direct edge from r to p is a particular case of this generic path which generates a zero row in B S|-R . Moreover, choosing r = p does not alter the given matrix form. This completes the second part and the proof.

Remark 1 (Switching to rooted trees/forests). Notice that in both theorems (and also for the whole thesis), we work on undirected graphs. Thus we omit the orientation within the trees and forests for the sake of simplicity. However these objects are still considered to be rooted trees/forests by setting the corresponding orientation. For example, in Theorem 2.1.1, the root of the spanning trees counted by the determinant is the deleted row/column index i i.e. |T i | = det L -i|-i . This is correct because of the following lemma. Lemma 2.1.4 (Unique Orientation). Let τ be a undirected tree and let us assign an orientation to each edge. Then, we define an orientation of a tree as the set of the directed edges in it. Over all possible orientations for τ , there exists only one unique orientation for an arbitrary node i ∈ V such that τ and this orientation together yield a tree rooted in node i.

Proof. Consider a simpler case where τ is a single path. In this case, the only orientation that makes node i the root is the orientation in which all edges are oriented towards i. Then, for more generic trees, we recall that in a tree, two nodes are connected only by a single path. As we have one particular orientation for each path, putting them together yields one unique orientation for τ to be rooted in i.

Thus, one has

|T i | = |T | = det L -i|-i .
In the case of forests and Theorem 2.1.3, the vertex set R is considered as the root set since they are ensured to be in separate connected components. Then the number of all spanning forests with the root set R reads:

|F R | = det L -R .

Random Walks on Graphs

A random walk W on a graph is a stochastic process that starts from a node s, and at every step it goes to another node with a probability w(i,j) d i . As a result, the probability of having a fixed walk ω in this process equals to:

P(W = ω|W 0 = s) = (i,j)∈ω w(i, j) d i , (2.2)
Such random walks naturally are discrete Markov chain.

Graph Theory

Definition 2.1.11 (Discrete Markov Chain). A discrete Markov chain M (i) is a random process over a discrete state space V that satisfies the Markovian property for all k ∈ N, i, j k ∈ V:

P(M (k) = i|M (k -1) = j k-1 , . . . , M (1) = j 1 ) = P(M (k) = i|M (k -1) = j k-1 ).
Markov chains are often described by their transition probabilities i.e. the probability of switching from one state to another. A transition matrix P = [P(M (k) = i|M (k -1) = j)] i,j∈V is the matrix that contains these probabilities. We associate a Markov chain with a graph by setting the vertices as the state space and the edges as the transitions with non-zero probabilities. One can generate random walks distributed according to Eq. (2.2) by simulating the Markov chain with the following transition probabilities:

P(M (k) = j|M (k -1) = i) = w(i, j) d i . (2.3)
To better understand Markov chains, we revisit some of their useful properties. We list them in the following with short descriptions:

• Periodicity: The period of a state x is defined as:

p(x) = gcd(n ∈ N + : P n (x, x) > 0).
where gcd returns the greatest common integer divisor of a given set. A state x in a Markov chain is called periodic if p(x) > 1. Otherwise, it is called aperiodic. If all states are aperiodic, the Markov chain is called aperiodic.

• Transient/Recurrent States: If a state i is recurrent, then the random walks starting i almost surely return to i while the number of steps goes to infinity. Otherwise, it is transient.

• Ergodicity: A state is ergodic if it is aperiodic and recurrent. If all states are ergodic, then the Markov chain is called ergodic.

In ergodic Markov chains, the probability distribution of M (k) over all states converges to a distribution as k goes to infinity [START_REF] Lai | Markov chains[END_REF]. This distribution is also called the stationary distribution and it verifies:

π P = π with i∈V π i = 1.
In other words, the stationary distribution of an ergodic Markov chain is the right eigenvector of its transition matrix. This result becomes intuitive when we look at the algebraic view of random walks. Let us consider a random walk starting from node i and write the vector δ i ∈ R |V| which is the Kronecker delta. After taking k steps, the probability distribution of the current state of the random walk reads:

r k = δ i P k .
The stationary distribution, in fact, is the limit of r k as k goes to infinity. By the Perron-Frobenius theorem, this limit converges to the right eigenvector of P.

Loop-Erased Random Walks

One can find many kinds of random walks by slightly deviating from the original definition in Eq. (2.2). However, all of them do not necessarily come with the Markovian property. Due to this, their distribution often is not mathematically tractable whenever we diverge from the original random walks. Thus any theoretical analysis usually cannot go further. Loop-erased random walks (Loop-erased Random Walks (LERW)s) are fascinating exceptions in this manner. They were first proposed by [START_REF] Francis | A SELF-AVOIDING RANDOM WALK[END_REF] in polymer physics. Since then, they have been "arguably the most tractable model among non-Gaussian models in statistical physics [START_REF] Kozma | The scaling limit of loop-erased random walk in three dimensions[END_REF]". More related to this thesis, they play a significant role in the sampling of Uniform Spanning Tree (UST) i.e. uniformly sampled spanning trees on a graph. [START_REF] Burton | Local characteristics, entropy and limit theorems for spanning trees and domino tilings via transfer-impedances[END_REF] shows that a path from node i to j in a UST has the same distribution as LERWs. In his seminal paper [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF], Wilson leveraged this fact to propose the current most efficient exact algorithm to generate uniform spanning trees. [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF] later gave the expected complexity of this algorithm by using the law of LERWs and graph matrices such as W or L. All of these connections lie at the very heart of this thesis. Therefore, we conclude our tour on graph theory by introducing some of the essential results related to LERWs.

Let us start with the formal definition of LERWs:

(v i , v i+1 ).
Then we define a new subset of indices I = {i 1 , i 2 , . . . , i k }'s of 1, . . . , l+ 1 that satisfies:

i 1 = 1 i j+1 = max{i : v i = v i j } + 1.
(2.4)

From the subset I, we form a loop-erased random walk LE(W ) ⊆ E as follows:

LE(W ) = {(v i 1 , v i 2 ), (v i 2 , v i 3 ), . . . , (v i k-1 , v i k )}.
As illustrated in 2.6, LERWs are random walks in which the cycles are deleted in the order that they show up (chronological order). In comparison to usual random walks, this description can seem complicated. However, its probability law surprisingly remains tractable in terms of the graph matrices.

Theorem 2.1.5 (Law of LERWs [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF]). A loop-erased random walk LE(W ) on G = (V, E, w) that is stopped at the subset ∆ ⊂ V has the following probability distribution:

P(LE(W ) = γ) = det L -∆∪s(γ) det L -∆ (i,j)∈γ w(i, j)
where γ is a fixed path and s(γ) denotes the nodes visited in γ.

Proof. See Lemma 1 in [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF].

This beautiful result ties not only diverse concepts such as probability theory and graph-related algebra but also becomes handy in the development and analysis of some of well-known graph algorithms in computer science. [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF] uses this result as a building block for the analysis of Wilson's algorithm as Wilson's algorithm repeatedly uses LERWs interrupted at some boundary. We will in Section 2.4.1 further detail this analysis along with the formal description of Wilson's algorithm.

Determinantal Point Processes

DPPs appear in most diverse phenomena. One of the earliest occurrences is in quantum physics; The particles called fermions obey Pauli's exclusion principle, which states that two fermions can not occupy the same quantum state. In [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF],

Macchi describes this behaviour of fermions by a DPP. Another early example is found in random matrix theory i.e. a field of mathematics that examines the properties of matrices whose entries are random variables [START_REF] Edelman | Random matrix theory[END_REF]. DPPs arise in the spectral analysis of the random matrices whose each entry is an independent random variable with complex normal distribution. The random eigenvalues of such matrices are distributed according to a DPP [START_REF] Johansson | Random matrices and determinantal processes[END_REF]. Another interesting example pops up in number theory. Consider a random sequence S = {X 1 , . . . , X N } of N numbers drawn between 0 and 9 uniformly. The locations whenever there is a descent in this sequence i.e. {i ∈ S|X i+1 ≤ X i }, is surprisingly a DPP. [KT+12] lists many interesting examples. A notable instance is associated with uniform spanning trees on graphs. The edges in a UST follow a DPP law which leads to tractable probabilities by using algebraic graph theory. We will detail this example and extend it for weighted trees in Section 2.4. Before doing so, we cover the necessary formal definitions and properties of DPPs.

Basics of DPPs

A point process X over a finite and discrete set Ω is a random subset of Ω. DPPs are point processes with a particular analytical form for their probability distribution.

Definition 2.2.1 (DPP). A point process X over a finite and discrete set Ω is called a determinantal point process if its inclusion probabilities i.e. the probability of any subset in X , verify for every fixed subset S ∈ Ω:

P(S ⊆ X ) = det K S
where K S denotes the submatrix of the semi-definite positive matrix K ∈ R |Ω|×|Ω| restricted to the rows and columns indexed by S. The matrix K is also called marginal kernel and its spectrum verifies 0 K I.

This compact description of the inclusion probabilities yield many tractable properties. For example, if we take two objects i, j ∈ Ω and take a look on their inclusion probability, we see:

P(i, j ∈ X ) = det K {i,j} = K i,i K i,j K j,i K j,j = K i,i K j,j -K i,j K j,i = P(i ∈ X )P(j ∈ X ) -K i,j K j,i
(2.5)

For simplicity, we assume K is a symmetric matrix. Then, one can easily see the negative association:

P(i, j ∈ X ) ≤ P(i ∈ X )P(j ∈ X ) P(i ∈ X |j ∈ X ) ≤ P(i ∈ X ).
Notice that having j ∈ X decreases the probability of also having i ∈ X . One can see this as a repulsive process where j repels i.

A bit technical but useful property follows the definition:

Lemma 2.2.1 (Restriction [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]). Given a DPP X on a ground set Ω with a marginal kernel K and a fixed subset α ⊆ Ω, X ∩ α is also drawn from a DPP with marginal kernel K α .

In other words, a restriction of a DPP is also a DPP.

The choice of the kernel matrix K can be various. The similarity matrices, in which each entry depicts the mutual similarity between the objects, are popular choices. In this case, as the similarity K i,j between two objects i and j grows, we end up with lower joint probabilities of i and j, which endorses the repulsion of two similar objects, in other words, diversity. However, notice that the choice of K such that all of its principal minors are probabilities is not obvious. Therefore, L-ensembles, a more natural subclass of DPPs, are often preferred:

Definition 2.2.2 (L-ensemble). An L-ensemble X is a point process associated with a matrix L ∈ R n×n which satisfies:

∀S ⊆ Ω, P(X = S) ∝ det L S ,
where L is a positive semi-definite matrix.

An L-ensemble is a DPP2 described by its atomic probabilities, the probability of sampling a fixed element or subset, which are easier to grasp for those outside this field. In addition, proportionality gives a higher degree of freedom by dropping some of the constraints of a proper matrix to design the DPP. Moreover, one can recover the corresponding K from a given L (See Thm. 2.2 in [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]):

K = L(L + I) -1 = I -(I + L) -1 . (2.6)
An immediate result of this identity is that K and L are spanned by the same eigenvectors. Let us denote the eigen-pairs of K by (ξ i , x i ) i∈Ω and the eigenvalues of L by 1 ≤ 2 ≤ . . . n . Notice that each x i is also the eigenvector of L with the eigenvalue:

∀i ∈ Ω, i = ξ i 1 -ξ i .
The sampling algorithm in [HKPV06; KT+12] uses this link to sample a DPP given L. This algorithm is exact but requires the eigendecomposition of L as input. Therefore, one needs to diagonalize a n × n matrix as a preprocessing, which yields a time complexity O(n 3 ). This cost is avoided for large n ( i.e. n = 10 5 ) by low-rank approximations or other approximate sampling algorithms providing a significant speed-up [START_REF] Derezinski | Exact sampling of determinantal point processes with sublinear time preprocessing[END_REF].

The tractable properties of DPPs are not limited to the inclusion and atomic probabilities. One can calculate the expected number of items in X :

Proposition 2.2.2 ([KT+12]

). The expectation and variance of |X | is as follows:

E[|X |] = tr(K), Var(|X |) = tr(K -K 2 ) (2.7)
Proof. The proof of expectation is one line:

E[|X |] = E   i∈Ω I(i ∈ X )   = i∈Ω P(i ∈ X ) = i∈Ω K i,i = tr(K).
(2.8)

One can derive the variance as follows:

Var(|X |) = E[|X | 2 ] -E[|X |] 2 = E[|X | 2 ] - i K i,i 2 = E      i∈Ω I(i ∈ X )   2    - i∈Ω K 2 i,i - i∈Ω j∈Ω\i K i,i K j,j = E   i∈Ω I(i ∈ X ) + i∈Ω j∈Ω\i I(i, j ∈ X )   - i∈Ω K 2 i,i - i∈Ω j∈Ω\i K i,i K j,j = i∈Ω P(i ∈ X ) + i∈Ω j∈Ω\i P(i, j ∈ X ) - i∈Ω K 2 i,i - i∈Ω j∈Ω\i K i,i K j,j = i∈Ω K i,i + i∈Ω j∈Ω\i (K i,i K j,j -K 2 i,j ) - i∈Ω K 2 i,i - i∈Ω j∈Ω\i K i,i K j,j =   i∈Ω K i,i   -   i∈Ω j∈Ω\i K 2 i,j   -   i∈Ω K 2 i,i   = tr(K -K 2 )
(2.9)

In terms of the eigenvalues of L, one has:

E [|X |] = i∈Ω i i + 1 with Var (|X |) = i∈Ω i ( i + 1) 2 .
Other tractable properties such as conditional or marginal probabilities are also available for DPPs. We refer to [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] for more details.

A fixed-sized DPP, also called k-DPP, is a determinantal point process conditioned over its size. Unlike DPPs, their probabilities are not interpretable sole determinants but ratios of determinants.

Definition 2.2.3 (Fixed-size DPPs). A fixed-size DPP is a special point process in which the number of items is set to a fixed number k. The atomic probabilities for a fixed set S with |S| = k read:

P(X = S) = det L S S ⊆Ω |S |=k det L S .
The fixed size DPPs give a rise naturally in certain applications. For example, consider the cellular network problem introduced in [START_REF] Bardenet | On a few statistical applications of determinantal point processes[END_REF]. The goal is to place antennas over a field so that every location receives the most reception that it can get. A natural solution to this problem is to model and sample antenna locations from a DPP. Imagine that we also have a strict constraint on the number of antennas that we can place due to the budget or environmental issues. In this case, k-DPPs come into play as we want to place exactly k antennas in the most diverse locations.

We Proof. See result 1 in [START_REF] Barthelmé | Asymptotic equivalence of fixed-size and varying-size determinantal point processes[END_REF].

The uniform spanning tree process on a graph is one instance in which projection DPPs show up. In a nutshell, given a graph with n nodes, one considers the random subset of edges that forms a uniform spanning tree. As aforementioned, this random process is precisely described by a DPP. Moreover, its kernel is a projection matrix with n -1 non-zero eigenvalues. In addition, a graph-theoretic verification for k = n -1 is that the number of edges needed to generate a tree that reaches all n vertices of a graph without creating a cycle is exactly n -1. These connections are some of the support pillars of this thesis. In Section 2.4, we detail these connections with formal descriptions.

Randomized Linear Algebra

In the last two decades, the volume of numerical data has exploded, and the need for efficient matrix algorithms has become more and more prominent. However, the classical and frequently used tools of linear algebra might not provide the required efficiency for the big data. For example, matrix inversion is essential in solving many machine learning problems. Inverting an n × n square matrix requires O(n 3 ) operations via the Gauss-Jordan elimination method. For sparse matrices, this cost is often avoided via Cholesky decomposition, which remains in O(n 3 ) in the worst case.

When it comes to dealing with billions of data points i.e. n = 10 9 , this computation might take days, thus be impractical. Along with the computational issues, there are also hardware limitations. For example, as the data grows exponentially, it becomes less and less feasible to fit the whole data to the RAMs (random access memories) of today's computers. Due to this, the computer needs to access the main memory multiple times, yielding runtime inefficiencies.

Randomized linear algebra (RLA) addresses both of these concerns by introducing randomness to linear algebra (LA). RLA proposes to imitate some of the main operations in LA with a randomized scheme, also called sketch-and-solve in [START_REF] Drineas | RandNLA: randomized numerical linear algebra[END_REF]; first, we sample a smaller portion of the data/the corresponding matrix; then, we apply the expensive operation on this portion which can be handled cheaply. We obtain by this way an estimate of the result.

This scheme has been applied to some of the main problems in LA, including solving dense/sparse linear systems, least-squares (also regularized), eigen or singular value decomposition [START_REF] Martinsson | Randomized numerical linear algebra: Foundations and algorithms[END_REF]. They have been able to outperform direct computations in solving some of these problems e.g. least-squares [START_REF] Drineas | Faster least squares approximation[END_REF] or low-rank approximation problems [START_REF] Drineas | Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix[END_REF]. In the scope of this thesis, we take look on a more recent RLA algorithm [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] for tackling regularized least-squares problems. This algorithm propose a DPP for creating sketches such that the corresponding estimates of the result are an unbiased estimate i.e. the expectation of the estimates equals the exact result. This algorithm can also be seen as a basis for the methods that will be introduced in Chapter 3. Therefore, we will revisit the least-squares problem and the DPP-based RLA algorithm in the following section.

30

Chapter 2 Background

Linear Least-Squares Problem

In linear least-squares problems, given the observations b ∈ R n and the predictors A ∈ R n×p , we seek the solution of:

x = argmin x∈R p ||Ax -b|| 2 .
(2.10)

Depending on n and p, the solution takes a different form. In our case, we restrict ourselves to the overdetermined case i.e. n > p where we usually have the solution which has a closed-form solution as follows:

x = (A A) -1 A b = A † b,
where A † denotes the Moore-Penrose inverse of A and we assume that A is a tall and full-rank matrix, thus (A A) is invertible.

A well-known application of this formulation is linear regression [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. In this case, each row of A contains a data point a i where a i ∈ R p and b consists of an observation per data point. The main goal is to learn the linear relation between the data points a i 's and the observations b i 's. To do so, we fit a hyperplane (a line in case of p = 1) that is parameterized by x and minimizes ||Ax-b|| 2 . Real life examples are limitless and each is not necessarily interesting. For solely entertainment purposes, in Example 2.3.1, we take a liberty of generating a fictional example.

Mission: Killing the Sun

Fig. 2.7.: A man shooting at the sun

According to the local news in 2016 [Hur], two men who were overwhelmed by the extremely high temperatures during summer in Adana/-Turkey started to take shots at the sun for the purpose of killing it. While the local police arrested them for the public shooting, the vital status of the sun is reported to be okay.

The mayor invited the public on their social media accounts to keep their calm in the face of these sweltering temperatures and advised that shooting at the sun would not be a feasible solution. We imagine an amusing line-fitting problem between the temperatures (data points) in Adana during the summer, and the number of the bullets fired at the sun (observation).

Our purpose here is to find the line that best explains the relation between the temperature and the number of bullets. Solving this problem under the least-squares formulation, one fits a line as shown in Fig 2 .8.

In many real life applications, the solution of least-squares suffers from ill-posed problems i.e. A is not full rank, and poor generalization i.e. overfitting. In the former case (A A) becomes singular thus cannot be inverted. The latter is a long standing issue in many machine learning models. In case of Least Squares (LS), the minimizer x might take a form which is extremely successful for minimizing ||Ax -b|| but completely useless to predict the observation for unseen data. This is because x can take extreme values to minimize ||Ax -b|| which is not necessarily robust depending on A and b. The regularization addresses both issues. The main idea is to add an additional term to the optimization formula to control x such that it is well-defined and bounded.

In case of L 2 regularization with a regularization matrix R ∈ R p×p , the problem becomes:

x = argmin x∈R p ||Ax -b|| 2 + x Rx (2.11)
and the solution still preserves a closed form as:

x = (A A + R) -1 A b, (2.12)
As interpretable and straightforward as it is, this solution comes with computational issues as p grows. Inverting a p × p matrix takes O(p 3 ) operations via direct computations of the Gauss-Jordan elimination method. Some of the intermediate calculations (e.g. pivoting) in Gauss-Jordan can be avoided via Cholesky decomposition [START_REF] Saad | Iterative methods for sparse linear systems[END_REF], yet the number of elementary operations remains in the cubic order of p. As the exact methods become inefficient whenever p is very large, state-of-the-art methods are the approximate ones in which we estimate the result much more cheaply with a small approximation error > 0. Depending on the algorithmic properties and the application, there are several schools; such as centralized/distributed [START_REF] Mou | A distributed algorithm for efficiently solving linear equations and its applications (special issue jcw)[END_REF] or deterministic/randomized methods [MT20] etc. We revisit in the following section some of these methods.

Approximate Methods for Solving LS Problem

Aligned with the main themes of this thesis, we will examine approximate methods in two groups; deterministic and randomized (or stochastic). Within this taxonomy, we classify the algorithms according to whether they contain a randomized element or not. In our tour over the most used algorithms, we reproduce some of them and compare their strength and weaknesses.

Deterministic Methods

Iterative solvers are the most used algorithms that fall under the deterministic methods for solving linear systems. In a nutshell, they use an iteration scheme that approximates the exact solution. They usually start with an arbitrary initial solution3 , and per iteration, the solution is updated such that it gets closer and closer to the exact solution. The iterative algorithms are fast to solve linear systems in form of Ax = b by using only matrix multiplications with A. In this way, they benefit from the sparsity if A is a sparse matrix i.e. number of non-zero entries is in the same order of magnitude as p. In addition, their approximation error = ||x -x|| is tractable in terms of eigen/singular values of A. Thanks to this fact, one can easily derive theoretical guarantees for convergence. In the following, we will shortly introduce the main ideas of these methods and reproduce some of the well-known algorithms and results.

Let us define the square matrix B := A A + R and y := A b ∈ R p and rewrite the solution in 2.12 as: x = B -1 y.

Iterative solvers vary in computing x depending on how they formulate the update f : R p → R p within the following iteration scheme:

x k+1 := f (x k ), k = 1, 2, . . . (2.13)
In our case, it is possible to draw an iteration for a given splitting B = M -N:

x k+1 = M -1 Nx k + M -1 y, x k+1 = x k + M -1 (y -Bx k ).
(2.14)

The basic iterative methods, e.g. Jacobi, Gauss-Seidel, successive over-relaxation (SOR) or Richardson, split B such that M can be inverted much easier than B. For example, the Jacobi iteration takes the following form:

x k+1 := D -1 (L + U)x k + D -1 y,
where we split Proof. We refer the reader to [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] for the detailed proof.

B = D -L -U,
This result indicates that not all splittings yield convergent iterations to the solution. In particular, one needs to know the spectral radius of the corresponding iteration matrix to guarantee the convergence. As ρ(G) is not easily accessible in practice, [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] suggests using ||G|| 2 > ρ(G) as an upper bound.

In fact, this result extends on the convergence rate i.e. a quantity to measure how fast the iterations converges to the solution. Formally, we define the convergence rate as follows:

µ := lim k→∞ max x 1 ∈R p ||x k -x || ||x 1 -x || 1 k .
For all the basic iterative methods above, this limit equals to ρ(G) [START_REF] Saad | Iterative methods for sparse linear systems[END_REF].

The basic iterative solvers are simple to implement and analyze. However, they are not applicable whenever the matrix B does not satisfy the constraints that are necessary for convergence e.g. diagonal dominance for Jacobi iterations. More robust methods with faster convergence guarantees are the Krylov subspace methods. These methods have a vast literature and practice [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. In the following, we summarize the main theory and reproduce one of the most used algorithms, called the conjugate gradient method.

Definition 2.3.1 (Krylov Subspace). A Krylov subspace K k (B, r) based on the matrix B and the vector r is a subspace of R n with a size m ≤ n that is defined as:

K k (B, r) = span(r, Br, B 2 r, . . . , B k-1 r) Lemma 2.3.2 ([Saa03]). The solution x = B -1 b belongs to K n (B, y).
Proof. The minimal polynomial of B, denoted by q(B), is the unique polynomial with the minimal degree such that q(B) = 0. Given the eigenvalues of B λ 1 ≤ λ 2 ≤ . . . ≤ λ n , one can rewrite this polynomial as:

q(t) = n i=1 (t -λ n ).
Rewriting q(t), one has:

q(t) = α 0 + α 1 t + α 2 t 2 + . . . α n t n
with the constant term α 0 = n i=1 λ i = det B. When B is non-singular, α 0 is ensured to be non-zero. Thus, one rewrite q(A) = 0 as:

-α 0 I = α 1 B + α 2 B 2 + . . . B n -α 0 I = B(α 1 I + α 2 B + . . . B n-1 ) B -1 = - 1 α 0 (α 1 I + α 2 B + . . . B n-1 )
Then, multiplying both sides with the vector b finishes the proof:

B -1 b = - 1 α 0 (α 1 I + α 2 B + . . . B n-1 )y = x ∈ K n (B, y)
There are several algorithms proposed that are based on Krylov subspaces. The main idea is to approximate x iteratively by following a set of candidate solutions x i 's that 2.3 Randomized Linear Algebra are also belong to K k (B, b). The best-known ones are Arnoldi's method, generalized minimized residual method and conjugate gradient method Let us complete this section by reproducing conjugate gradient algorithm. A generic gradient descent algorithm minimizes the following loss function

F (x) := 1 2 x Bx -x y.
Notice that x minimizes F (x). Given an arbitrary initial solution x 0 , an intuitive way to approach x is to take a step from x 0 in the opposite direction of the gradient of F (x). An iteration that implements this idea is as follows:

x k+1 = x k -α∇F (x k ) (2.15)
where the gradient ∇F (x) reads Bxy and α ∈ R adjusts the size of the step. On top of this, the conjugate gradient method suggests taking orthogonal directions with optimal step sizes such that there will be no need to take another step in the direction of previous steps. As a result, the updates form a Krylov subspace based on K k (B, b).

We summarize the algorithm in Alg. 1. The convergence rate of this algorithm is

Algorithm 1 Conjugate Gradient Algorithm 1: Inputs: B ∈ R n×n , b ∈ R n # Input matrix and vector x 1 ∈ R n # Initial point K ∈ N + , # Maximum number of iterations ∈ R + # Tolerance 2: Initialize: r 0 ← b -Bx 1 # Initialize the residual p 0 ← r 0 3: for k ← 1 to K do 4: α k ← r k r k r k Br k # Calculate the step size 5: x k+1 ← x k + α k p k # Update the iterate 6: r k+1 ← r k -α k Bp k # Update the residual 7: if ||r k+1 || 2 ≤ then # If
β k ← r k+1 r k+1 r k r k 11: p k+1 ← r k+1 + β k p k 12: end for 13: return x k+1 given as 1-2 √ κ i.e. at every iteration k, one has ||x k+1 -x || 2 ≤ 1 -2 √ κ k ||x 1 -x || 2
where κ = λmax(B) λ min (B is the condition number.

Iterative algorithms cover a significant part of approximate methods for solving least-squares problems. They provide fast convergence with a small approximation error, while the main computational cost is due to the matrix-vector products with B. Thus they scale linearly with the number of non-zero entries in B. Moreover, it is possible to improve their performance via preconditioning [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. On the other hand, we note that, for using CG, B is restricted to the real, symmetric and positive semi-definite matrices. In addition, there are some cases where the product with B is also not available. For example, when B is too big, its vector products might be costly to evaluate.

Another branch of approximate methods, so-called randomized methods, is better adapted for dealing with these issues by leveraging probability theory and Monte Carlo simulations. In the following section, we revisit some of the relevant ones with this thesis.

Randomized Algorithms

The earliest randomized algorithms to solve linear systems in the form of Bx = y are attributed to Ulam and von Neumann and later developed by [START_REF] George | Matrix inversion by a Monte Carlo method[END_REF]. These algorithms are based on random walks on a Markov chain whose initial distribution and transition matrix are closely related to B and y, respectively. However, these algorithms require strict constraints on the spectrum of B to guarantee convergence to the desired solution. In addition, how to design the Markov chain to minimize the approximation error is still an open question. We refer to [Ökt05] for more details.

A set of algorithms given by the RLA literature addresses a broader class of matrices in the setup of the least square problem. Generically, these algorithms use the sketch-and-solve scheme to give an estimate x that is an ( , δ)-approximation, i.e.:

||Ax -b|| 2 ≤ (1 + )||Ax -b|| 2 with probability 1 -δ, ( 2.16)
A generic description of such algorithms given by [START_REF] Drineas | RandNLA: randomized numerical linear algebra[END_REF] is as follows:

Consider the ordinary least square problem in Eq 2.10 (without regularization) where the solution is the pseudo-inverse A † b. Let us define a random variable J ∈ {1, 2, . . . , n} over the row indices of A with a probability mass function p : {1, 2, . . . , n} → (0, 1) that verifies n i=1 p(i) = 1. Then, a meta-algorithm to estimate x takes the following form:

• Get k samples j 1 , j 2 , . . . , j k of J according to p with replacement,

Randomized Linear Algebra

• Form the matrix à such that ∀i ∈ {1, . . . , k}, Ãi|: = 1 kp(j i ) A j i |: and the vector b such that ∀i ∈ {1, . . . , k}, bi = 1 kp(j i ) b j i

• Solve the following least-squares problem to get approximation of x :

x = argmin x∈R p || Ãx -b|| 2 = Ã † b.
(2.17)

As typically k n, computing x can be handled much more cheaply than the original problem. On the other hand, the performance of this algorithm strongly depends on the probability mass function. [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] analyzes this scheme with some of the conventional probability mass functions e.g. uniform or weighted by row norms. This analysis concludes that the algorithm described above produces ( , 0.1)approximations for some k n. This is to say that one can avoid the expensive inverse operation by sampling rows of A and b and operating on these rows which, in turn, yields a fairly small approximation error.

Closer to the main themes of this thesis, more recent RLA algorithms deploy determinantal point processes [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF]. These methods give rise to the elegant connections between the determinantal point processes and linear algebra. As a result, they benefit from certain theoretical guarantees such as unbiasedness.

[DM21] proposes two new DPP-based estimators for solving the ordinary and regularized least-squares problem. The main idea in these estimators is to sample rows by a fixed size DPP such that the solution to the Eq. (2.17) yields an unbiased estimation of x . In case of the ordinary LS, the suggested DPP is a fixed size L-ensemble with size of p and L = AA . Theorem 2 in [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] states that the corresponding DPP X verifies:

E[A -1 X b X ] = argmin x∈R p ||Ax -b|| 2 = A † b,
where A X is the submatrix of A restricted to the rows by X . The second estimator uses the same idea to solve the regularized least-squares. In this case, X takes form of a DPP (not fixed size) with L = 1 λ AA to sample row indices. Then, it verifies:

E[A † X b X ] = argmin x∈R p ||Ax -b|| 2 + λ||x|| 2 = (A A + λI) -1 Ab.
In the actual problem, λ is the hyper-parameter to adjust the balance between the confidence on the predictors/observations and the regularization. As λ increases, the solution will be more regularized and vice versa. Similarly, as λ increases, the expected size of X becomes smaller. This means that the estimate A † X b X relies less on the predictor A and the observations b. For further analysis, such as the exact error or algorithmic complexity, we refer the reader to [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF].

These estimators give a simple and elegant solution to least-squares problems when the direct or other approximate approaches are expensive. Within the contributions of this thesis, we will present in Chapter 3 a cheap method of reducing the expected error (or variance) of these unbiased estimators. Surprisingly, this error reduction method connects diverse concepts such as gradient descent optimization and variance reduction Monte Carlo simulations. In addition, we will also present some estimators for graph-regularized least-squares problems. Surprisingly, these estimators can be reconsidered as the DPP-based estimators explained above.

Random Spanning Forests

Similar to spanning trees, combinatorialists have been studying spanning forests to extract information about the graph. The earliest results on the enumeration, and also randomization, over spanning trees/forests dates back to Kirchhoff's matrix tree theorem [START_REF] Kirchhoff | Ueber die Auflösung der Gleichungen, auf welche man bei der Untersuchung der linearen Vertheilung galvanischer Ströme geführt wird[END_REF]. Later on, these results emerged in physics with exciting links to certain random models in statistical mechanics. A significant instance (for this thesis) is that uniform spanning trees/forests are particular cases of a celebrated model in statistical physics, called random-cluster model [START_REF] Grimmett | The random-cluster model[END_REF] which unifies the percolation theory and Ising-Potts' model i.e. a statistical model for ferromagnetic materials. In parallel, probabilists have extended the fruitful links between USTs (also Uniform Spanning Forest (USF)) and algebraic graph theory highlighting the connections with random walks, harmonic functions and electrical networks [START_REF] Benjamini | Special invited paper: uniform spanning forests[END_REF]. Then computer scientists have come into play and developed efficient algorithms for sampling USTs [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF]. As Wilson's algorithm [BLPS01; Ald90] remains the most efficient exact algorithm (it exactly generates a spanning tree from the uniform distribution) in the past three decades, more recent algorithms trade the exactness with the time complexity and propose approximate sampling algorithms [START_REF] Kelner | Faster generation of random spanning trees[END_REF].

With all these rich theoretical properties, USFs (along with USTs) have become apt tools to analyze complex graphs. In this thesis as well, we aim to leverage these properties for randomized linear algebra. Therefore, we dedicate the rest of this section to introduce random spanning forests (a more generic form of uniform spanning forests) and bring their deep connections with the graph Laplacian into the light.

Given a graph G = (V, E, w), consider all rooted spanning forests F over G. The cardinal of this set |F| may vary a lot. The exact expression reads |F| = det(L + I) by extending the matrix-tree theorem [CA02]. As a result, the number of spanning forests on a graph with n can grow up to n n+2 (in the case of a complete graph) depending on the structure. As this number is usually very large for general graphs, analyzing a graph by enumerating all spanning forests remains a complicated task.

A down-to-earth alternative is by randomization, defining the distribution of forests by associating probabilities on every spanning forest in F. Such distributions over forests have been thoroughly studied by [START_REF] Benjamini | Special invited paper: uniform spanning forests[END_REF] with the limits in case of infinite graphs i.e. n → ∞. In this thesis, we only study finite graphs and focus on the following random objects studied in [ACGM18]:

Definition 2.4.1 (RSF). A random (rooted) spanning forest is a random directed subgraph whose outcome space is F with the following parametric distribution:

P (Φ Q = φ) = 1 Z Q i∈ρ(φ) q i (i,j)∈φ w(i, j) (2.18)
where Q = {q 1 , . . . , q n } is a set of real parameters that verify min i q i > 0 and Z Q is the normalization constant over all forests. In case of q 1 = q 2 = . . . = q n = q, our notation for the corresponding RSF is Φ q .

A particular case occurs when some q i 's are set to 0 for a subset of the nodes V ⊂ V by leaving at least one node with a non-zero parameter. The roots of the forests generated by this configuration are always contained in the set V \ V . Another interesting case is as q i → ∞ for i ∈ V ⊆ V. In this case, V is ensured to be a subset of ρ(Φ Q ) while the distribution takes the following form:

P (Φ Q = φ) ∝ i∈ρ(φ)\V q i w(φ),
where w(φ) = (i,j)∈φ . The RSF Φ Q is a fascinating object with many interesting properties. As it is of central importance in this thesis, let us reproduce some of the significant results related to Φ Q .

Wilson's algorithm

One of the biggest advantages of Φ Q is the existence of a simple and efficient algorithm to sample it. This algorithm is due to an adaptation of Wilson's algo-rithm [START_REF] Avena | Random forests and networks analysis[END_REF]. Given a root node r, Wilson's algorithm originally generates samples of the random spanning tree T with the following distribution:

P(T = τ ) ∝ (i,j)∈τ w(i, j), τ ∈ T r . (2.19)
where T r is the set of all spanning trees rooted in r. This algorithm uses LERWs to generate the paths of a spanning tree as follows:

1. Initialize the boundary ∆ = {r} and the tree T = ∅, 2. Run a random walk starting from an arbitrary node in V \ ∆, 3. Interrupt the random walk whenever it hits ∆, 4. Erase the cycles in the order that they appeared to obtain a path γ, 5. Update T and ∆ by T ← T ∪ γ and ∆ ← ∆ ∪ s(γ) where s(γ) denotes the nodes included in the walk γ, 6. If ∆ \ V = ∅, terminate the algorithm and return the tree T . Otherwise go to step 2.

This algorithm keeps updating T with the loop erased walks in G until all vertices are covered. Therefore, when it terminates, the result in T is a spanning tree. Wilson's algorithm given in Alg. 2 in fact is an efficient implementation of this description.

Here we note that RandomSuccessor returns a random neighbor v of the given node u with the probability w(u,v) du .

Having these detailed descriptions at hand, in the rest, we will answer three important questions about Wilson's algorithm:

1. How does one arrive at the distribution in Eq. (2.19) from this algorithm? 2. What is the time complexity of the algorithm?

3. How can we adapt this algorithm to sample Φ Q ?

One answer, that is often cited, to the first question was given by Wilson in his seminal paper. This proof heavily uses a representation of the random process conducted by the algorithm, which is called the stack representation. The stack representation proposed by [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF], is a sophisticated tool to analyze Wilson's algorithm. In this analysis, one considers a procedure over the stack representation, called cycle popping algorithm. This algorithm also terminates with a spanning tree and the distribution of the trees is the same as the spanning tree distribution induced by Wilson's algorithm. [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF] derives this distribution as in Eq. (2.19). A

Algorithm 2 RandomTreeWithRoot [Wil96] 1: Inputs: G = (V, E, w) r ∈ V 2: Initialize: ∀i ∈ V, InTree[i] ← false InTree[r] ← true ∀i ∈ V, Next[i] ← -1 3: for i ← 1 to |V| do 4: u ← i 5: while not InTree[u] do 6: Next[u] ← RandomSuccessor(u, G) 7: u ← Next[u] 8: end while 9: u ← i 10: while not InTree[u] do 11: InTree[u] ← true 12: u ← Next[u] 13:
end while 14: end for 15: return Next more down-to-earth approach [Mar00] uses the probability law of LERW given in Theorem 2.1.5. In this case, one considers Wilson's algorithm as a series of LERWs interrupted at evolving boundaries through the algorithm. Then, one derives the law of the spanning trees by using the law of LERWS given in 2.1.5. In the following, we first go through this analysis to give a detailed answer to question 1 as it is more compatible with the main themes of this thesis. Later, we take a slight detour to formally introduce the stack representation and cycle popping algorithm (along with the links to Wilson's algorithm) because we will need them in the following chapters.

LERWs in Wilson's Algorithm

Given the root r, k = 0, T 0 = ∅ and s(T 0 ) = {r}, Wilson's algorithm can be considered in two steps:

1. Run a LERW that starts from an arbitrary V \ s(T k ) and is interrupted at s(T k ).

2. Let γ be the path generated by the LERW until interruption. Then, update

• T k+1 ← T k ∪ γ, • s(T k+1 ) ← s(T k ) ∪ s(γ), • k ← k + 1,
and go to step 1. In Wilson's algorithm, these two steps are repeated until s(T k ) = V or equivalently T k is a spanning tree. In other words, each branch in T k is a LERW interrupted at certain boundaries. In the following, we use this fact to derive the distribution for T k .

Theorem 2.4.1 (Correctness of Wilson's algorithm [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF]). Wilson's algorithm generates random spanning trees distributed according to the law of Eq. (2.19).

Proof. At an arbitrary step k of Wilson's algorithm, we sample a branch of a spanning tree by using a LERW interrupted at s(T k-1 ). The distribution of such random walks read:

P(T k = γ ∪ T k-1 |T k-1 ) = det L -s(T k-1 )∪s(γ) det L -s(T k-1 ) (i,j)∈γ w(i, j).
Assume that the algorithm terminates at step K. Then, the probability law of the resulting tree is:

P(T K = τ ) = K k=1 P(T k = γ ∪ T k-1 |T k-1 ).
After telescopic cancellations, one recovers the distribution:

P(T K = τ ) = K k=1 det L -s(T k-1 )∪s(γ k ) det L -s(T k-1 ) (i,j)∈γ k w(i, j) = 1 det L -r (i,j)∈τ w(i, j).
Noticing that det L -r is the normalization constant (same for all r ∈ V) for the random spanning trees finishes the proof.

Here the tractable law of LERWs allows us to show the correctness of Wilson's algorithm i.e. it samples from the correct distribution. Moreover, this analysis is easy to adapt for different types of probabilistic event, such as the probability of having a tree (not spanning) in a random spanning tree. We refer the reader to [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF] for more details.

2.4 Random Spanning Forests

Stack Representation and Cycle Popping Algorithm

Define a random infinite stack (a collection of items that are inserted and deleted in last-in-first-out order)

S (i) = [S (i) 1 , S (i) 2 . . .] per node i which verifies P(S (i) l = k) = w(i,k) d i
or alternatively S (i) j = RandomSuccessor(G, i) at every level j. Now we build a directed graph G S by looking at the top level i.e. l = 1 of all stacks. Formally,

G S = (V, E S , w) such that E S = {(i, j) : S (i) 1 = j}.
An immediate result is that G S has n -1 edges that form cycles and paths and if there is no cycle, then G S is a spanning tree. Given these facts, we define a cycle-popping algorithm that pops the stacks associated with the cycles until there is none left, and we are left with a random spanning tree. Wilson's algorithm in fact is an implementation of the described cycle-popping algorithm. The only difference in Wilson's algorithm is the fact that we erase the cycles in the order they appear whereas in the cycle popping algorithm described, the order can be arbitrary. The following lemma shows that the order of the cycles popped is irrelevant to the resulting spanning tree. Therefore, the result of these two algorithms is the same random spanning tree.

Lemma 2.4.2 (Invariance to the order of cycles [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF]). If the cycle-popping algorithm terminates, then the order of cycles popped does not alter the resultant spanning tree.

Proof. Assign a color l to each edge (i, S (i) l ) according to their level in the stack. Then we define a colored cycle c in G S as a cycle formed by colored edges. Note that two colored cycle might have the same set of edges, yet they are not the same as long as they contain a different set of colors on the edges.

Let c 1 , c 2 , . . . , c k , c be a sequence of colored cycles that are observed during the execution of the cycle popping algorithm. We prove the theorem by showing that the cycle c will be popped in case of popping any other feasible sequence.

To do so, we look at simpler case: Assume that c is a cycle that can be popped instead of c 1 . The question that we need to answer is:

If we start from c instead of c 1 , can we still pop c?

We answer this question by considering the following cases:

• Case 1: c and c 1,...,k shares no common vertex. In this case, the answer is yes because after popping c , one can still pop the sequence c 1,...,k and reach c.

• Case 2: c and c 1,...,k has some common vertices. Let c i be the first cycle that shares some vertices with c i.e. i = min{i ∈ N + : c i ∩ c = ∅}. If c and c i do not have the same edges, then they share at least one vertex v that has a different successor. However as the stack of v has not been popped by c 1 , . . . , c i-1 up until c i , its successor in c i must be the same with its successor in c . Applying this reasoning to all vertices in c ∩ c i shows that c and c i have the same edges. Moreover, the edges have the same set of colors, because the cycles c 1 , . . . , c i-1 share no vertex with c and c i . Therefore, the colors on the edges of c and c i remains untouched.

Then under the facts c = c i and c 1,...,i-1 ∩ c = ∅, we may follow the order

c = c i , c 1 , . . . , c i-1 , c i+1 , . . . , c k by reaching c.
With these two cases, we generalize that if a cycle c is encoded within the stacks, the cycle popping algorithm will eventually pop it regardless of the order that we follow to pop cycles. In case that the algorithm terminates, the resultant structure, a spanning tree, is invariant from the order of popped cycles.

By using this lemma, Wilson shows that the resulting tree structure by the cycle popping/his algorithm is distributed according to Eq. (2.19). As we already recovered this result in Theorem 2.4.1, we prefer not to repeat it in this manuscript.

Second Question. Now, we shift our focus to the second question; what is the time complexity of Wilson's algorithm? In this case, we refer by the time complexity to the expected number of calls of RandomSuccessor which is the most repeatedly called function in the algorithm. [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF] shows that this expectation is equal to the mean commute time ν, i.e. the expected number of steps of a random walk in a Markov chain until it returns to its initial state of the Markov chain used in Wilson's algorithm. The mean commute time is defined as:

ν := i∈V π i (t i,j + t j,i )
where π is the stationary distribution of the Markov chain and t i,j is the expected number of steps to reach j from i. In case of Wilson's algorithm, the Markov chain admits the transition probabilities P i,j = w(i,j)

d i except that it is absorbed at r. 2.4 Random Spanning Forests S = S (1) S (2) S (3) S (4) S (5)               2 1 4 3 ∅ 5 3 4 2 ∅ 2 4 5 1 ∅ 4 2 5 1 ∅ . . .

Fig. 2.9.: A toy example on the order invariance of the popped cycles. The infinite matrix S

contains the initial configuration of the infinite stacks in its columns. As we pop the cycles, the corresponding stacks pop the element at their first row. We follow different orders for popping the cycles but the popped cycles in the are the same ones yielding the same spanning tree as a result.

Then [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF] calculates the expected number of calls to RandomSuccessor in terms of P as follows:

E[# of calls to RandomSuccessor] = ν = tr(I -P -r ) -1 .
(2.20)

where P -r is the reduced matrix by deleting the rows and columns corresponding to r.

Third Question. Finally, adapting Wilson's algorithm for spanning forests is straightforward: • Add an additional node called Γ to G,

Γ 1 4 3 2 q 3 q 4 q 2 q 1 Γ 1 4 3 2 1 4 3 2 Wilson's algorithm Cut the edges to Γ
• Add an edge from every node to Γ with an edge weight w(i, Γ) = q i ,

• Run Wilson's algorithm on the extended graph by choosing Γ as the root.

Alg. 3 gives an implementation of this adaptation (also See Fig. 2.10). This adaptation samples a random spanning tree on an extended graph with an additional node Γ. The distribution of these trees reads:

P(T Γ = τ ) ∝ (i,j)∈τ w(i, j)
Notice that by cutting the edges in T Γ that are incident to Γ, one obtains a spanning forest. Moreover, the relation between these forests and trees is a bijection i.e. for each instance of T Γ , one reaches a unique forest by cutting the edges incident to Γ. Then these forests admit the same distribution as T Γ . By separating the edges linked to Γ, one recovers:

P(Φ q = φ) ∝ i∈ρ(φ) q i (i,j)∈φ w(i, j).
Finally, we obtain the time complexity of RandomForest by adapting Eq. (2.20):

E    # of calls to RandomSuccessor in RandomForest    = tr (Q + L) -1 (Q + D) ≤ n + 2m min i q i . (2.21)
The upper-bound we give here is to give an understanding of this unknown trace in terms of graph parameters. However, it is usually not tight. We refer the reader to App. A.1 for the derivation of the upper bound.

DPPs in Random Spanning Forests

In the scope of this thesis, the most interesting properties of RSFs are their links with DPPs. In particular, the edges and roots of Φ q are sampled from certain DPPs whose kernel matrices are tractable in terms of the edge incidence matrix or graph Laplacian. These results have their roots from the results in random spanning trees [START_REF] Burton | Local characteristics, entropy and limit theorems for spanning trees and domino tilings via transfer-impedances[END_REF], later summarized in [START_REF] Avena | Random forests and networks analysis[END_REF]. We give in this section some important ones.

Random Spanning Forests

Algorithm 3 RandomForest [START_REF] Avena | Random forests and networks analysis[END_REF] 1: Inputs:

G = (V, E, w) Q = {q 1 , . . . , q n } 2: Initialize: # Initially, the forest is empty ∀i ∈ V, InForest[i] ← false ∀i ∈ V, Next[i] ← -1 ∀i ∈ V, d[i] ← j∈N (i) w(i, j) # Degrees 3: for i ← 1 to |V| do 4: u ← i 5:
# Start a random walk to create a forest branch 6:

while not InForest[u] do # Stop if u is in the forest 7: if rand ≤ q i q i +d[u] then #If true, u becomes a root 8: InForest[u] ← true # Add u to the forest 9: Next[u] ← -1 # Set next
Q are DPPs). Given a graph G = (V, E, w), the roots ρ(Φ Q ) and edges S of Φ Q is an L-ensemble with: ∀V ⊆ V, S ⊆ E, P((ρ(Φ Q ), S) = (V, S)) ∝ det L V ∪S with L = Q 1/2 B Q 1/2 B = Q Q 1/2 B BQ 1/2 BB , where Q = diag(Q).
Proof. We prove this theorem algebraically by showing that for any fixed set V ⊂ V and S ⊂ E, the determinant det L V ∪S is proportional to the distribution of RSFs (at the right hand side):

det L V ∪S ∝ P(ρ(Φ Q ), S = V, S) = 1 Z Q                    r∈V q r (i,j)∈S
w(i, j), if S forms a spanning forest with the root set V , 0, otherwise.

Recall that the adaptation of Wilson's algorithm for forests samples a spanning tree over an extended graph with a node Γ. Let us start by writing the edge incidence matrix of this graph:

B ext = Q 1/2 -Q 1/2 1 B 0 ∈ R (m+n)×(n+1)
An immediate result of this definition links B ext with L:

L = (B ext :|-Γ ) B ext :|-Γ .
In addition, given a fixed set S of edges on the extended graph, Theorem 2.1.3 yields

| det B ext S |-Γ | =              (i,j)∈S w(i, j), if S forms a span- ning tree rooted in Γ, 0,
otherwise.

S consists of the edges incident to Γ and the edges from E. As aforementioned, when we delete the former group, we obtain a rooted spanning forest φ = 2.4 Random Spanning Forests (V, S, w) with roots V are the neighbors of Γ in S . Therefore, there is a bijection between the roots V and the edges incident to Γ in S . Relying on this fact we rewrite:

| det B ext S∪V |-Γ | =              r∈V √ q r (i,j)∈S
w(i, j), if S forms a spanning forest with the root set V , 0, otherwise.

Then we finish the proof by writing the determinant formula for matrix multiplications:

det L V ∪S = det(B ext V ∪S|-Γ ) 2 =              r∈V q r (i,j)∈S
w(i, j), if S forms a spanning forest with the root set V , 0, otherwise.

. One recovers the probability P((ρ(Φ Q ), S) = (V, S)) by normalizing the right hand side by Z Q .

This process of roots and edges is in fact a projection DPP based on L-ensemble L by Lemma 2.2.4. In order to see this, first notice that the roots-edges process is a fixed size DPP as |S| + ρ(Φ Q ) = n is always true. Moreover, one has rank L = n. 4 Then the marginal kernel of the projection DPP can be written as M = n k=1

x k x T k where x k is the k-th eigenvector of L associated to a non-zero eigenvalue.

Remark 2 (Restriction to the roots). Restricting a DPP to a subset of its outcome space also yields a DPP restricting its marginal kernel (See Lemma 2.2.1). We can calculate the marginal kernel of ρ(Φ Q ) by using this property.

Theorem 2.4.4. ρ(Φ Q ) is a DPP with marginal kernel K = (L + Q) -1 Q.
Proof. The proof have three main steps:

1. We first derive the marginal kernel of the root-edge process in terms of blocks of L.

2. Then we invoke Lemma 2.2.1 to derive the kernel only for the roots ρ(Φ Q ).

3. Finally, we show a DPP with the resulting kernel is equivalent to a DPP with the kernel K = (L + Q) -1 Q.

For the first step, let us diagonalize L = XDX where X ∈ R (n+m)×n contains the non-trivial eigenvectors and D ∈ R n×n is a diagonal matrix whose entries contain the non-zero eigenvalues. The marginal kernel of the root-edge process reads:

M = XX .
In order to write this matrix in terms of the blocks of L, we use the following trick:

lim α→∞ αL(αL + I) -1 = XX . a
In left hand side, we write the block matrix form of L within the matrix inverse:

lim α→∞ αL αQ + I αQ 1/2 B αBQ 1/2 I + αBB -1 = L lim α→∞ Q + α -1 I Q 1/2 B BQ 1/2 α -1 I + BB -1
Before going further calculations with the block matrix inversion formula, recall that we want to derive the kernel only for the root set. This means by Lemma 2.2.1 that we are interested in deriving only the first n × n block of M. Then let us restrict our calculations on the region of interests. This means that by the block matrix inversion, one has:

M {1,...,n} = Q Q 1/2 B lim α→∞ S -(BB + α -1 I) -1 BQ 1/2 S where S = (Q + α -1 I -Q 1/2 B (BB + α -1 I) -1 BQ 1/2 ) -1
is the Schur's complement. After several algebraic manipulations (See Appendix A.3), this limit converges to:

M {1,...,n} = lim α→∞ (Q -Q 1/2 B (BB + α -1 I) -1 BQ 1/2 )S = Q 1/2 (L + Q) -1 Q 1/2 .
(2.22) We finish the proof by showing the equivalence of the two kernels M {1,...,n} and K. As both (L + Q) and Q 1/2 are square matrices, one has for all subsets X ⊆ V:

det(Q 1/2 (L + Q) -1 Q 1/2 ) X = det(Q 1/2 ) X det(L + Q) -1 ) X det(Q 1/2 ) X = det(Q) X det(L + Q) -1 ) X = det(Q(L + Q) -1 ) X = det K X

Random Spanning Forests

a This is easy to check diagonalizing αL(αL + I) -1 as X(I + (αD) -1 ) -1 X . As α → ∞, the limit converges to XX .

The connection of K with RSFs goes even beyond. The following theorem indicates that the off-diagonal entries K i,j corresponds to certain probabilities in RSFs.

Theorem 2.4.5. For a given graph G = (V, E, w), in a realization of RSF Φ Q , the probability of node i to be rooted in node j reads:

P (r Φ Q (i) = j) = K i,j (2.23)
Proof. Let F ij be all spanning forests in which node i is rooted at j. Then, one writes the probability P (r Φ Q (i) = j) by using the definition of Φ Q :

P (r Φ Q (i) = j) = φ∈F ij l∈ρ(φ) q l a,b∈φ w(a, b) φ∈F l∈ρ(φ) q l a,b∈φ w(a, b) = q j φ∈F ij l∈ρ(φ)/j q l a,b∈φ w(a, b) φ∈F l∈ρ(φ) q l a,b∈φ w(a, b) = q j A ij Z (2.24)
where A ij and Z are scalars. In the rest, we leverage the revisited theorems to give close form expressions for A ij and Z.

Given the extended graph G = (V , E , w) and a subset S = S i ∪ S Γ with |S| = |V| -1, one can write the resultant absolute determinant for R = {i, Γ} by using Theorem 2.1.3:

| det B {S|-i,-Γ} | =   a,b∈S Γ w(a, b) a,b∈S i w(a, b)   1/2
whenever S forms a spanning forest on G with two trees in which node i and Γ are disconnected. Similarly, one obtains the same result by plugging node j instead of i. Then, the following holds:

| det B {S|-i,-Γ} det B {S|-j,-Γ} | =   a,b∈S Γ w(a, b) a,b∈S i =S j w(a, b)   52
Chapter 2 Background whenever S forms a spanning forest on G with two trees in which Γ is not connected to either i or j. This implies that i and j are contained the same tree. Notice that adding the edge (j, Γ) to such subset S reconstructs a spanning tree on Γ. Then, one can acquire a rooted spanning forest from this tree as done in the sampling procedure. In the resultant forest, the node j is set as the root of the tree that contains both i and j. In other words, it is an element of the forest set F ij . Moreover, the edge weight product of this forest is equal to

q j | det B {S|-i,-Γ} det B {S|-j,-Γ} |.
Then, the sum of the edge weight products of all forests in F ij equals to:

φ∈F ij l∈ρ(φ) q l a,b∈φ w(a, b) = q j S⊆E ,|S|=|V |-1 | det B S|-i,-Γ det B S|-j,-Γ | (2.25)
By using Lemma A.4.1, we exchange the absolute value operator with the sign of the expression. Then, by using the Cauchey-Binet, one has

φ∈F ij l∈ρ(φ) q l a,b∈φ w(a, b) = (-1) i+j q j S⊆E ,|S|=|V |-1 det B S|-i,-Γ det B S|-j,-Γ = (-1) i+j q j det L -i,-Γ|-j,-Γ (2.26)
where L is the graph Laplacian of G . This completes the computation of A ij in Eq. (2.24). In a similar fashion, one computes Z as det L -Γ|-Γ by doing the same calculations without i and j. Noticing L -Γ|-Γ = (L + Q) reads the famous Cramer's formula for inverse matrices:

P (r Φ Q (i) = j) = A ij Z = q j (-1) i+j det(L + Q) -i,|-j det(L + Q) = q j (L + Q) -1 i,j = K i,j

Random Partitions induced by Φ Q

A partition of the vertex set is any set

P := (V 1 , . . . , V K ) that satisfies ∀k, k ∈ {1, . . . , K}, V k ∩ V k = ∅ and ∪ K k=1 V k = V.
A straight-forward partition in a spanning forest φ is given by its connected components, or equivalently trees. We denote this partition by π(φ

) := (V 1 , V 2 , . . . , V |ρ(φ)| ).
Let us enumerate the trees of φ from 1 to |ρ(φ)|. We define the tree number of a node as the number of the tree that it belongs to and it is denoted by t(π(φ), i) where t : (π(F), V) → {1, . . . , |ρ(φ)|} maps every node to its tree number.

Deleting the edges between the parts Fig. 2.11.: Fixing the partition created by Φ Q ensures that the edges between parts are never sampled. Moreover, every tree sampled in such process has to be a spanning tree of a part. Therefore, given the partition, sampling Φ Q boils down to sampling independent spanning trees in each part.

Another surprising theoretical result on Φ Q emerges when we look at the conditional probabilities of roots over the random partition π(Φ Q ):

Proposition 2.4.6 (Roots in a fixed partition [START_REF] Avena | Random forests and networks analysis[END_REF]). For a fixed subset V ⊆ V and a fixed partition P, the conditional probability of roots for a fixed partition of random spanning forests is:

P(ρ(Φ Q ) = V |π(Φ Q ) = P) = u∈V q u v∈V t(P,v) q v .
(2.27)

Proof. Consider the subgraphs of G by deleting the edges in between the subsets given by the fixed partition (See Fig. 2.11). Notice that sampling Φ Q on a fixed partition P is equivalent to sampling random spanning trees independently on each such subgraph. Therefore we need to look at the probability distribution of the roots in these trees to derive the conditional probability. To do so, let us denote the random spanning tree restricting Φ Q to only spanning trees by T . The probability distribution of T reads:

P(T = τ ) ∝ (i,j)∈τ w(i, j)
Then, the probability of T being rooted in rooted in node r is:

P(T ∈ T r ) = q r τr∈Tr (i,j)∈τr w(i, j) r ∈V q r τ r ∈T r (i,j)∈τ r w(i, j) .
(2.28)

Recall that given a undirected spanning tree, for each node r, there exists exactly one orientation of edges which makes r root. Therefore, the sum 54 Chapter 2 Background τr∈Tr (i,j)∈τr w(i, j) is a constant for all r ∈ V. Then, after the cancellation in both numerator and denominator, one has:

P(T ∈ T r ) = q r r ∈V q r .
As this is the probability of each tree sampled in Φ Q conditioned over the partition P, recalling the independent nature of the trees finishes the proof:

P(ρ(Φ Q ) = V |π(Φ Q ) = P) = u∈V q u v∈V t(P,v) q v .
This last result on RSFs covers all the technical results that we need for the rest. In the next section, we take a break before proceeding to the main materials of the thesis and summarize some of the important results.

Conclusion

Overall, we have covered many results from diverse concepts in this section. We finalize this chapter by highlighting the ones that will be repeatedly used in the following chapters. Therefore, we kindly suggest the reader to keep these following couple of pages at reach as a cheat-sheet through the rest of the journey.

Graph Laplacian

Given a graph G = (V, E, w), the graph Laplacian L = D -W is a simple but One matrix to rule them all... resourceful object in graph theory, combinatorics, probability theory, computer science, signal processing and machine learning. Closer to the main themes of this thesis, in this section, we take a look on two significant results on certain determinants of L:

We first look at its link with the spanning trees that can be found on G. This link is the celebrated matrix-tree theorem which gives the number of all the spanning trees on a unweighted graph in a simple closed-form expression:

∀i ∈ V, |T | = det L -i|-i .

Conclusion

We later generalize it for weighted graphs:

∀i ∈ V, det L -i|-i = τ ∈T (i,j)∈Eτ w(i, j).
These results extend on rooted spanning forests on G. In Theorem 2.1.3, we adapt Poincaré's results for spanning forests that yields a useful device to count all spanning forests on G given the root set R:

∀R ⊆ V, det L -R = φ∈F R (i,j)∈E φ w(i, j).
Second, we give the links of L with a special type of random walk, called loop-erased random walks. These random walks are obtained by erasing cycles as they appear in the usual random walks. Due to the cycle erasure, they are not Markovian random processes i.e. the next vertex to land does not solely depend on the current vertex. Unlike many non-Markovian process, LERWs surprising have a tractable distribution, which can be written in terms of principle minors of L:

P(LE(W ) = γ) = det L -∆∪s(γ) det L -∆ (i,j)∈γ w(i, j)
where ∆ ⊂ V is the boundary where the LERW LE(W ) is stopped, γ is a fixed path and s(γ) denotes the nodes visited in γ.

Determinantal Point Processes and Randomized Linear Algebra

A determinantal point process X is a special type of point process over a set Ω that verifies:

P(S ⊆ X ) = det K S .
for some square matrix K. This definition also yields that:

E[|X |] = E   i∈Ω I(i ∈ X )   = i∈Ω P(i ∈ X ) = i∈Ω K i,i = tr(K). 56 Chapter 2 Background A = Sample rows         † b
Fig. 2.12.: An illustration of the sample-and-solve scheme. First we sample a small portion of the matrix A, then proceed to the expensive operation.

The key feature of DPPs is that they can encode negative association between different elements in Ω. For a symmetric K, one can examine this phenomena by looking at the joint probabilities in the case of two elements i, j ∈ Ω:

P(i, j ∈ X ) = P(i ∈ X )P(j ∈ X ) -K i,j K j,i ≤ P(i ∈ X )P(j ∈ X ).
As the calculation suggests, having node i in X reduce the chances of also having node j in X .

Although DPPs by definition seem a bit abstract, real-life examples and use cases of DPPs are numerous. In our case, we are interested in their use in randomized linear algebra. RLA is a set of numerical linear algebra tools that replaces expensive direct operations with randomized approximations. We consider solving least-squares problem (also regularized) with such approximation scheme. In this case, we are given some observations b ∈ R n and predictors A ∈ R n×p and we solve:

x = argmin x∈R p ||Ax -b|| 2 .
The exact solution x reads:

x = A † b.
In the regularized case, it is:

x = (A A + λR) -1 A b.
Generically these inversions requires O(np 2 ) computations with the direct methods. RLA proposes randomized methods to avoid this expensive computation. In a nutshell, these methods first sample a smaller portion i.e. some rows of A and (also b depending on the randomization) at random. Then, they do the expensive operations on the sampled portion which avoid the expensive global operations (Also see Fig. 2.12). Indeed, the type of the randomization i.e. the probability

Conclusion

distribution is of central importance in such a scheme. This is where DPPs come into play. Sampling the rows from certain DPPs yield unbiased estimators for both leastsquares and regularized least-squares solution. Closer to this, one main contribution of this thesis is a new estimator whose variance is significantly reduced w.r.t. the DPP-based estimators (More details in Section 2.2). In addition, we in Chapter 3 introduce elegant estimators for graph-regularized least-squares problems i.e. the solution is regularized to be smooth on a graph. Surprisingly these estimators are also DPP-based estimators as the random process (random spanning forests) that they use is a DPP.

Random Spanning Forests

The last but not least concept we present in this chapter is the random spanning forests i.e. process of randomly selecting spanning forest on a given graph. In fact, most of the theoretical results we rely on in this thesis is closely related with the RSFs. We denote a RSF by Φ Q and write its distribution as follows:

P (Φ Q = φ) = 1 Z Q i∈ρ(φ) q i (i,j)∈φ w(i, j)
where Q = {q 1 , . . . , q n } with min i q i > 0 and Z Q is the normalization constant over all forests. Let us highlight some of the significant properties and theorems of Φ Q :

• Wilson's algorithm for RSFs. Wilson's algorithm is originally developed for generating spanning trees at uniform. With a slight modification, it can be used for generating samples of Φ Q . We give this modified algorithm Alg. 3, we also note the time complexity of this algorithm in terms of number of random variables sampled:

E    # of calls to RandomSuccessor in RandomForest    = tr (Q + L) -1 (Q + D) ,
where Q = diag(Q). While it is much more efficient than naive sampling i.e. enumerating all forests and choosing one at random, it is known as the most efficient algorithm for generating spanning/trees.

• The root-edge process. The roots and edges in Φ Q are sampled from an L-ensemble with:

∀V ⊆ V, S ⊆ E, P((ρ(Φ Q ), S) = (V, S)) ∝ det L V ∪S with L = Q 1/2 B Q 1/2 B = Q Q 1/2 B BQ 1/2 BB
A consequence of this along with the restriction lemma (Lemma 2.2.1) is that the root set ρ(Φ Q ) is a DPP with the marginal kernel

K = (Q + L) -1 Q. • P(r Φ Q (i) = j) = K i,j
. This identity is another surprising result between spanning forests/trees and Laplacian-based algebra. Its proof is an algebraic proof mainly based on Poincaré's theorem for counting spanning forests with root constraints5 .

• Roots in fixed partitions. The final theorem we visit gives the distribution of the roots when we generate a sample of Φ Q conditioned over a fixed partition i.e. Φ Q |π(Φ Q ) = P. If the partition is fixed, sampling Φ Q boil downs to sampling independent trees within each subset of the partition. Then within each subset, the probability having a node i as the root becomes

q i j∈V t(P,i) q j
where V t(P,i) is the subset that contains i. By independence of the trees, one has:

P(ρ(Φ Q ) = V |π(Φ Q ) = P) = u∈V q u v∈V t(P,v) q v .
Note that for q 1 = q 2 = • • • = q n = q, this boils down to a uniform distribution.

With this last point, we reach to the end of our short(!) summary on the technical background for the next material.

Graph Tikhonov

Regularization and Interpolation with RSFs 3 " C'est par la logique qu'on démontre, c'est par l'intuition qu'on invente.

-Henri Poincaré

This chapter is dedicated to presenting the contributions of the thesis for solving the Graph Tikhonov Regularization (GTR) (a graph adaptation of regularized linear regression) and Graph Interpolation (GI) problems. These contributions are twofold; i/ the RSF-based Monte Carlo estimators for solving GTR and GI and ii/ their improved versions via variance reduction techniques. In order to do so, we first give formal definitions of these problems and show how they relate to each other. Then we leverage the elegant theory of RSFs to define the methods proposed that yield unbiased estimates for these problems. In parallel, we give several ways to improve the estimators. Under certain conditions, these techniques can substantially decrease their variance with a cheap additional computational cost. In addition, we find several other use cases for these estimators where one needs to avoid expensive computations within intermediate operations, such as semi-supervised learning on graphs, Newton's method and Alternating Direction Method of Multipliers (ADMM).

In parallel, we also provide empirical illustrations or comparisons of the techniques proposed with state-of-the-art methods.

Problem Definition

As in many signal processing applications, the measurements over graph signals are either noisy or incomplete i.e. missing over some subset of nodes. A big chunk of tools to overcome these issues proposes to smooth out the signal over the given graph. In the case of denoising, smoothing means that we eliminate specific components of the measurements that are highly varying over the graph by assuming that the original signal is smooth on the graph. Similarly, in an interpolation setup, we complete the missing parts of the signal so that it remains smooth on the underlying graph.

A unifying approach tackles both problems with a single optimization problem. This approach adapts the Tikhonov regularization for graph signals and draws the following optimization problem: where ||a|| Q = a Qa is the Mahalonobis distance and Q is a diagonal matrix where each Q i,i is a node-wise regularization parameter. We define the selection matrix M ∈ R n×p as a rectangular diagonal matrix where M i,i = I(i ∈ V ). This problem admits a closed form solution in the following form:

x = (Q + L) -1 QMy (3.2)
This generic linear regression formulation, in fact, has been revisited in various contexts, from solving heat equations in physics to semi-supervised learning on graphs in machine learning. Among many, two cases exemplify the primary uses of GTR and GI:

Graph signal filtering: Graph signal filtering means eliminating certain components

In this case, one considers

M = I ∈ R n×n and Q = qI ∈ R n×n
of the graph signal and retain the parts of interest. The common practice in linear filtering schemes is through the frequency analysis of the signal i.e. eliminating certain frequency components. A fundamental tool in this type of analysis is the

Problem Definition

Fourier transform that interprets the signal in terms of its frequency components. However, the adaptation of the Fourier transform for graph signals is not evident due to the irregularity of the signal domain. By analogy, there exist multiple definitions of Fourier transform [RBTGV19; TGB18]. The one that we are interested in, takes the eigenvalues Λ = {λ 1 , . . . , λ n } and the eigenvectors U = [u 1 | . . . |u n ] of the graph Laplacian as the graph frequencies and the graph Fourier basis, respectively. Then the graph Fourier transform of a signal is defined as:

ŷ := U y. (3.3)
The vector ŷ contains the Fourier coefficients of the signal corresponding to the different frequency components. One recovers the original signal via the inverse transform y = Uŷ. Given these definitions, we describe a filtering operation in three steps; i/ first compute the Fourier coefficients as ŷ = U y, ii/ then filter out these coefficients by computing ŷ = diag(g(Λ))ŷ where g : Λ → R n is the transfer function of the filter, iii/ finally compute the filtered signal y = Uŷ by going back to the original domain via the inverse transform. We can summarize these three steps in the following expression:

y = U diag(g(Λ))Uy.
Depending on the characteristic of the transfer function g, this filtering scheme eliminates some components and amplifies others. For example, for the ideal lowpass filters or k-bandlimited filters, which keep only the first k frequency components and eliminate the rest (of the high-frequency components), the transfer function becomes:

g k (λ) =    1 if λ < λ k , 0 otherwise.
Interestingly, the solution x = Ky where K = q(L + qI) -1 can be considered under this graph filtering scheme. This is because L and K share the same eigenvectors and one can write:

K = U diag(g(Λ))U with g q (λ) = q q + λ .
In other words, the solution x is a filtered graph signal with the transfer function g q (λ) = q q+λ . Moreover, g q (λ) is a low-pass filter as it diminishes for larger values of the graph frequencies. In Fig. 3.1, we illustrate this behavior on a graph built by temperature sensors over Bretagne/France. Here, the graph signal is the average temperature measurements collected from various locations in Bretagne. We apply the filter g q (λ) on this signal with three different values of q. We observe that 62 Chapter 3 Graph Tikhonov Regularization and Interpolation with RSFs for decreasing values of q, the filter eliminates more and more high frequency components. In turn, the output becomes smoother on the underlying graph.

Dirichlet boundary problem on graphs: In solving differential equations, the

In this case, for some subset

V ⊆ V one has M ∈ R n×|V | with Mi,i = I(i ∈ V ) and Q = limq→∞ qMM .
Dirichlet boundary condition specifies the values of the solution at a given boundary. As it is a generic constraint, it can be used with any differential system. However, we are interested in the case when it is imposed on Laplace's equation. Laplace's equation is a particular partial differential equation (PDE) in which we seek a multivariate function u(x 1 , . . . , x p ) that satisfies: where ∆ is the continuous Laplacian operator, δD is the boundary of the closed region D ⊂ R p and the function g gives the fixed boundary values. The heat diffusion, fluid flow, or electrostatics are only some of the physical phenomena in which this PDE or its variants appear [START_REF] Haeseler | Laplacians on infinite graphs: Dirichlet and Neumann boundary conditions[END_REF][START_REF] Martinetti | Designing the sound of a cut-off drum[END_REF]. By analogy, this problem translates to the discrete graph domain. In this case, the subset V ⊂ V is considered as the boundary and we seek a solution over U = V \ V of the following linear system:

∆u = p i=1 ∂ 2 u ∂x 2 i = 0, s.t. ∀x 1 , . . . , x p ∈ δD, u(x 1 , . . . , x p ) = g(x 1 , . . . , x p ).

Problem Definition

∀i ∈ U, j∈N (i) w(i, j)(x i -x j ) = 0, s.t. ∀i ∈ V, x i = g(i).
(3.4)

In turn, we have a closed-form solution for this problem:

x i =    ((L U |U ) -1 L U |V g) i , if i ∈ U g(i), otherwise, (3.5)
where g = [g(i)] i∈U . We illustrate this scheme in Fig. 3.2. In this illustration, we are given a graph with a lot of unlabeled nodes (unlabeled nodes in blue). We propagate the information contained in the labeled nodes by using the solution in Eq. 3.5. This solution is also called the harmonic function. A notable application of harmonic functions is the electrical networks i.e. networks of resistors (a circuit element that has a proportional relation between the potential difference applied to it and the current flowing through it). In these networks/graphs, each resistor is an edge, and each point where two resistors meet is a node. In this abstraction, the conductances of the resistors are denoted by the edge weights. A typical problem is the calculation of the currents flowing through the resistors whenever there is a potential difference between some nodes. This problem boils down to solving the problem in Eq. 3.4.

To show this, let us denote the potentials/voltages at each node i by x i and assume that the potentials in V are given by g(i)'s for all i ∈ V . By Ohm's law, the current flowing in edge (i, j) equals to w(i, j)(x i -x j ). Then, by Kirchoff's current law state that the current sum j∈N (i) w(i, j)(x i -x j ) at each node i must be equal to zero.

As a result, we recover the problem in Eq. (3.4).

State-of-the-Art

In these applications and in many more that boil down to this form, the efficient calculation of x is of utmost necessity. The main computational cost is inverting a matrix involving the graph Laplacian L. In the worst-case scenario, the naive computations go through the Gauss-Jordan elimination process that takes O(n 3 ) elementary operations. This cost is prohibitive as n increases. This shortcoming is solved by using either direct methods after useful matrix factorizations (e.g. Cholesky decomposition) or approximate methods (e.g. iterative solvers). The former approaches provide the exact solution by avoiding some cumbersome intermediate operations. They usually yield an efficient computation. However, in the worst case, the number of operations scales up to O(n 3 ). Therefore, state-of-the-art is the approximate methods. Among those, the strongest ones are the iterative solvers, notably the conjugate gradient descent algorithm with proper preconditioning. The time complexity of Conjugate Gradient Descent (CGD) scales linearly with the number of non-zero entries in the inverted matrix, which is the number of edges m in our case. More details on the taxonomy of the approximate methods are in Section 2.3.3. This chapter will present the RSF-based algorithms to approximate x. Appealing properties of these estimators are that:

• They illustrate a practical use of elegant links between RSFs and the graph Laplacian.

• They give unbiased estimates of x (i.e. their expectation equals to x),

• The time complexity of both estimators scale linearly with m,

3.2 State-of-the-Art
• The theoretical analysis is tractable.

On the other hand, the intrinsic weakness of these methods is their convergence rate to x. They are subject to the Monte Carlo convergence rate O(N -1/2 ) where N is the number of Monte Carlo samples. Therefore, they require many samples to reach high accuracy. As this issue remains the main limiting factor in all Monte Carlo estimators, the vast literature on Monte Carlo has focused on decreasing the expected error of the estimators. In order to do so, many algorithms have been proposed [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF]. The main line of the algorithms proposed suggests using an additional statistic to efficiently eliminate some of the variance of the estimator. We manage to adapt two of them for the RSF-based estimators. We detail these adaptation in Section 3.3.1 and 3.3.2.

Estimating x via RSFs

Consider the case V = V. Then, one recovers:

x = Ky with K = (Q + L) -1 Q.
Now, recall Theorem 2.4.4 which states that the matrix K is the kernel matrix of the DPP of the random roots of Φ Q . Moreover, it verifies the following identity:

P(r Φ Q (i) = j) = K i,j .
We present in this section two unbiased estimators of x that take their roots from this beautiful result. Moreover, they easily extend to the case where there are missing measurements i.e. V ⊂ V.

Definition 3.3.1 (x). Let y ∈ R n be the (complete) measurements over V. Define an estimator x = [x 1 , . . . , xn ] as follows:

∀i ∈ V, xi := y ir with i r = r Φ Q (i). (3.6)
The expectation and variance analysis is tractable and uncover the link between x i The estimator x is unbiased for x,

ii The weighted sum of the node-wise variance of x reads:

E[||x -x|| 2 Q ] = i∈V q i Var(x i ) = y (Q -K QK)y.
Proof. The proof of (i) is easy. Define a random matrix S :

= [I(r Φ Q (i) = j)] i,j .
An immediate result of this definition is:

E[ S] = K.
Noticing x = Sy finishes the first part of the proof:

E[x] = E[ S]y = Ky = x.
To prove (ii), we replace the variance with its definition:

i∈V q i Var(x i ) = i∈V q i E[x 2 i ] -E[x i ] 2 = i∈V q i   j∈V P(r Φ Q (i) = j)y 2 j -x2 i   . = i∈V q i   j∈V K i,j y 2 j -x2 i   .
Then define y (2) = [y 2 1 , . . . , y 2 n ] and q = [q 1 , . . . , q n ]. Rewriting the summation above, one has:

i∈V q i Var(x i ) = i∈V q i (Ky (2) ) i - i∈V q i (Ky) 2 i = q Ky (2) -yK QKy.

Estimating x via RSFs

67

Noticing that q is the eigenvector of (Q + L) -1 with the eigenvalue 1 finishes the proof:

i∈V q i Var(x i ) = y (Q -K QK)y.
Implementing x is straightforward; first sample a spanning forest via Wilson's algorithm, then within each tree of the forest, propagate the signal value of the root throughout the tree (See Fig. 3.3). In order to generate N samples of x, we need to sample N spanning forests via Wilson's algorithm. As this is the heaviest computation, one can calculate the time complexity of the proposed method in terms of the total number of steps taken by the random walks in Wilson's algorithm to sample N forests. Let us denote this number K N . By using Eq. (2.21), one has:

E[K N ] = N E[K 1 ] = N tr(Q + L) -1 (D + Q).
Notice that if y = c1, the variance of x is zero. Since K is a smoothing operator on the graph, an intuitive thinking is that there is nothing to smooth out in y. If there is some variation in y, the variance of x becomes non-zero. In fact it is maximized when y equals to the right eigenvector of K that corresponds to the smallest eigenvalue. In this case, the quadratic form y K QKy is minimized and the variance is maximized. An intuition for a such vector can be taken from the particular case Q = qI. As previously discussed, this is the case of graph signal filtering by using the Laplacian eigenvectors as the Fourier basis. Then the eigenvector discussed is u n the eigenvector of L with the largest eigenvalue λ n or the highest graph frequency. In other words, the vector which has the highest variation on the graph results in the highest variance in x (Also see Fig. 3.4).

It is possible to adapt x when we only have measurements over a subset V ⊂ V. In this case, M ∈ R n×p takes a form of a rectangular selection matrix. Let y V = My = [y , 0 ] ∈ R n Then it is possible to write x in the following form:

x = (Q + L) -1 Qy V .
Then, we can design x to estimate this inversion. In this configuration, we sample Φ Q by setting: In this example, we consider the case of graph signal filtering with q = 1. The variance of x then becomes y (I -K 2 )y. We plot the spectrum of the matrix (I -K 2 ) w.r.t. the graph Laplacian eigenvalues/graph frequencies. As a filter, (I -K 2 ) acts like a high-pass filter. Thus the quadratic formula y (I -K 2 )y will take higher values when y is a highly varying signal on the graph.

q i =    Q i,i , if i ∈ V,
and we calculate x over these forests. Notice that the roots of the forests generated by this configuration are always a subset of V . Therefore, the propagated values in x are always the known signal measurements. Indeed, x efficiently estimates x. However, there is still a big room for improvement without compromising the runtime efficiency. In our case, we find techniques that yield significant improvements in the literature on variance reduction techniques [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF]. Besides their cheap cost, these techniques give us a way to leverage our additional knowledge embedded in the nature of the RSFs ( e.g. certain statistics) to substantially decrease the expected error/variance of x. In the following, we delve into the application of these VR techniques, yielding two new estimators for x.

Variance Reduction via Conditional Monte Carlo

Let X be a random variable with an unknown expectation E[X] = µ X . The classical Monte Carlo algorithm would take N samples of X and compute the empirical mean:

X = 1 N N i=1 X i .
Now assume that there exists another random variable Y for which the conditional expectation E[X|Y ] is easy to calculate. The conditional Monte Carlo technique also

Estimating x via RSFs

Sample a forest

Average & Propagate Fig. 3.5.: An illustration of x. We follow the same first step with x. Then we propagate the (weighted) average signal within each tree. We assume

q 1 = • • • = q n = q.
called Rao-Blackwellization [Bla47; Rao45], allows us to use our knowledge on Y to reduce the error of the naive estimation on µ X . In doing so, it defines a new estimator as follows:

Z := 1 N N i=1 E[X i |Y ]. As E[E[X|Y ]] = E[X],
Z is an unbiased estimator. Moreover, it has a reduced variance due to the law of total variance:

Var(X) = E[Var(X|Y )] + Var(E[X|Y ]).
(3.7)

In our case, we find a simple way of deploying this idea for x:

Definition 3.3.2 (x). Given a fixed partition P, we define the conditional expectation x = [x 1 , . . . , xn ] ∈ R n as follows:

∀i ∈ V, xi := E[x i |π(Φ Q ) = P] = n j=1 y j P(r Φ Q (i) = j|π(Φ) = P) = n j=1
y j q j k∈V t(P,i)

q k = j∈V t(P,i)
q j y j k∈V t(P,i)

q k . (3.8)
In a practical sense, after sampling Φ Q , the new estimator x is simply the empirical weighted mean within each tree, thus it is easy to implement. Note that for q 1 = • • • = q n = q, it boils down to the simple averaging. This is illustrated in Fig. 3.5. After sampling a spanning forest, one only needs to compute local averages and propagate them within trees for calculating x. On the theory side, one can easily ensure that x is unbiased and has a lower variance than x by Eq. 3.7. Moreover, we have a closed form expression for the variance of the new estimator:

Proposition 3.3.2 (Variance of x). The weighted sum of the node-wise variance of x equals to:

E[||x -x|| 2 Q ] = y (QK -K QK)y.
Proof. We start by replacing the variance by its definition:

E[||x -x|| 2 Q ] = i∈V q i Var(x i ) = i∈V q i E[x 2 i ] -E[x i ] 2 , = i∈V q i E[x 2 i ] -||x|| 2 Q = i∈V q i E[x 2 i ] -y K QKy.
(3.9)

As in Prop. 3.3.1, we define a random matrix S that verifies x = Sy. Here S can be considered as a weighted averaging operator. Then, rewriting x = Sy, one has:

i∈V q i Var(x i ) = i∈V q i E[( Sy) 2 i ] -y K QKy, = (E[|| Sy|| 2 Q ]) -y K QKy.
(3.10)

Here we invoke the identity a SQ S = Q S:

i∈V q i Var(x i ) = E[y S Q Sy] -y K QKy, = y E[Q S]y -y K QKy,
= y (QK -K QK)y.

(3.11)

a This result might not be seem straight-forward but easy to check by calculating the results of the matrix products Similar to x, the variance of x is zero for a constant measurement vector y = c1 and it becomes non-zero when y contains some variation between nodes. However, x has a significant difference for highly varying signals. This is easy to observe in case of graph signal filtering which is illustrated in Fig. 3.6.

Overall, x gives theoretical guarantees for the improvement in estimating x while it does not compromise the run-time efficiency. One can calculate the local weighted averages by a single pass over the signal y, yielding O(n) time complexity, which is dominated by the cost of sampling a forest. In practice, we observe over various graphs that x brings substantial performance improvement. Nevertheless, we

I -K 2 K -K 2

Fig. 3.6.:

In this example, we again consider the case of graph signal filtering with q = 1.

The variance of x then becomes y (K -K 2 )y. We plot the spectrum of the matrix (K -K 2 ) w.r.t. the graph Laplacian eigenvalues/graph frequencies and compare it with the case of x. As a filter, (K -K 2 ) acts like a band-pass filter i.e. its transfer function only takes high values for a particular range of frequencies, called the pass-band, that are neither high, nor low. Thus the quadratic formula y (K -K 2 )y will take higher values whenever y is mostly composed by the frequency components that are in the pass-band. Also note that in comparison to x, x is always better as illustrated.

are capable of even going beyond this performance by applying another variance reduction technique. We detail this technique in the next section.

Variance Reduction via Control Variates

Consider again the problem of estimating µ X from the samples of the random variable X. The control variate (CV) method proposes to take another random variable Y , also called control variate, that has:

• a known or easy to calculate expectation µ Y ,

• non-zero correlation with the estimated quantity i.e. Cov(X, Y ) = 0.

Given this ingredient, we define another estimator as follows:

Z := X -α(Y -µ Y ),
where α ∈ R is a hyper-parameter. An immediate result is that Z is an unbiased estimator of µ X regardless of the value of α, as

E[Z] = E[X] -α(E[Y ] -µ Y ) = µ X .
However, the variance analysis is closely related to the value of α. Let us write the variance of Z in terms of the variance of X, Y and α:

Var(Z) = Var(X) + Var(αY ) -2 Cov(X, αY ), = Var(X) + α 2 Var(Y ) -2α Cov(X, Y ).
(3.12)

To ensure Var(Z) < Var(X), one has to verify that:

α 2 Var(Y ) -2α Cov(X, Y ) < 0 Assuming α > 0, one obtains: α < 2 Cov(X, Y ) Var(Y ) .
This bound ensures variance reduction. Yet, we can still go beyond this analysis and derive the optimum value for α. To do so, we take the derivative of Var(Z) with respect to α and set it to zero.

∂Var(Z) ∂α = 2α Var(Y ) -2 Cov(X, Y ) = 0.

Estimating x via RSFs

Solving this equation gives:

α = Cov(X, Y ) Var(Y ) . (3.13)
Plugging α , we get a reduced variance as follows:

Var(Z) = Var(X) - Cov(X, Y ) 2 Var(Y ) .
In many applications, α is not necessarily accessible as it is a ratio of two variances that are not often available. A common practice to overcome this issue is to estimate it from the samples of X and Y :

α = Cov(X, Y ) Var(Y ) ,
where Cov(X, Y ) and Var(Y ) are sample variance formulas. In doing so, one can either use extra samples of X and Y to compute α or use the same samples for computing both α and Z. The former option gives an unbiased estimator at the cost of having extra samples. In the latter, the resulting estimator might have a bias as α and Z are calculated from the same samples. However, in [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF], the author discusses that the bias diminishes faster than the Monte Carlo error w.r.t. the number of samples, i.e. O(N -1 ) compared to O(N -1/2 ). Therefore, as N grows, the bias becomes smaller compared to the Monte Carlo approximation error.

We find a simple adaptation of the control variate technique on the forest-based estimators x and x. This adaptation is mainly inspired by the gradient descent algorithm for solving Prob. 3.1.1. Let us briefly revisit this algorithm:

Definition 3.3.3 (Gradient Descent). Consider the following optimization problem:

x = argmin x∈R n F (x) s.t. F (x) = 1 2 x Kx -x y.
Notice that the solution of this problem is also x = Ky. The gradient descent algorithm suggests iteratively approximating the solution. It picks an arbitrary initial vector x 0 ∈ R n and updates it as follows: where ∇F (x k ) = (K -1 x k -y) and α ∈ R. This iteration updates the solution by taking a step in the opposite direction of the gradient of the loss function F (x) which is the steepest descent to the solution (See Fig. 3.7).

x k+1 := x k -α∇F (x k ), k = 1,
We present another unbiased estimator by proposing to use the gradient descent update as the control variate.

Definition 3.3.4. Given the forest estimate x, we apply a single step of gradient descent algorithm, which yields the following estimator:

z := x -α(K -1 x -y).
The unbiasedness of z is easy to see as

E[z] = E[x] -α(K -1 E[x] -y) = Ky.
The variance highly depends on the choice of α as discussed before. From Eq. (3.13), the optimal choice is:

α = tr(Cov(K -1 x, x)) tr(Cov(K -1 x)) = tr(K -1 Cov(x)) tr(K -2 Cov(x)) , 3.3 Estimating x via RSFs
where Cov(x) = [Cov(x i , xj )] i,j ∈ R n×n is the auto-covariance matrix of x. Computing α requires computing Cov(x) which is not available without expensive operations. On the other hand, one can calculate α from the samples of x much more cheaply. In numerical experiments, this choice of α often reduces the variance near optimal despite the aforementioned bias. Yet, one may still want to ensure the variance reduction regardless of how much it is reduced. When this is the case, we suggest using the classical arguments of gradient descent literature for the step size. By doing so, we find a cheap constant for α that guarantees variance reduction:

Proposition 3.3.3. Define the maximum degree of a graph as d max := max i∈V {d i }.

Fixing α = 2q q+2dmax guarantees that:

||z -x|| 2 ≤ ||x -x|| 2 . (3.14)
Proof. The inequality in (3.14) yields a step size such that the result gets closer to the solution at every step. Expanding the definition of z:

||x -x -α(K -1 x -y)|| 2 ≤ ||x -x|| 2 , ||(I -αK -1 )(x -x)|| 2 ≤ ||x -x|| 2 . (3.15)
This inequality is satisfied if the eigenvalues µ 1 , . . . , µ n of the matrix (I -αK -1 ) verifies:

|µ i | ≤ 1, ∀i ∈ {1, . . . , n},
In terms of the graph Laplacian eigenvalues, one has:

∀i ∈ {1, . . . , n}, |1 -α λ i q -α| ≤ 1, 0 ≤ α ≤ 2q λ i + q . 0 ≤ α ≤ 2q λ n + q .
Recalling λ n ≤ 2d max [START_REF] Van Mieghem | Graph spectra for complex networks[END_REF] finishes the proof.

Indeed, this choice of α is cheaper compared to α and introduces no biases. However, it might provide a much smaller variance reduction compared to optimal in certain cases. This is mostly because λ n ≤ 2d max is not necessarily a tight bound, and the constraint in (3.14) might be too excessive for taking a safe step, and thus it might slow down the algorithm. In order to clearly illustrate these differences, we consider the following case study:

Different Choices of α

In this example, we consider the signal denoising problem on two graphs generated by two random models, namely random k-regular and Barabasi-Albert. We are motivated to choose these models as they have opposite characteristics in terms of their degree distribution. While the random kregular model always generates graphs in which all degrees are fixed to k, the graphs generated by the Barabasi-Albert model typically have spread degree distributions. For simplicity, we take a random signal generated from N (0, I). We apply smoothing via the direct computation (to compute x) and the forest-based estimators x and z. In Fig. 3.8, we plot the mean squared error of these estimate with x for varying values of α. Min. error In these plots, we focus on three points that correspond to α (the optimal solution), α (the estimate from samples) and α = 2q q+2dmax (the constant from Prop. 3.3.3). In the case of a regular graph, both options for α give a very close performance to the optimal choice α . This is because λ n ≤ 2d max in a regular graph in this case gives a tight upper bound thus α = 2q q+2dmax is very close to the optimal choice. This drastically changes when we switch to a highly irregular graph. In this case, α continues to give a performance next to the optimal, whereas choosing α = 2q q+2dmax becomes too excessive for reducing the error.

Estimating x via RSFs

This covers our analysis of variance reduction techniques on the forest-based estimators. Indeed, the variance reduction literature is much broader than these two techniques. We observe two main issues in bringing other methods to play; the first and fundamental issue is that it is not often obvious how to adapt these methods efficiently for our particular case. The variance reduction methods often require some prior knowledge on the estimated quantity that is not necessarily available in our case. The second is that even if the necessary conditions for applying a VR method are satisfied, anticipating if it would substantially decrease the variance or not is another question. This requires a good understanding of the random variable and its distribution. In particular, we need to know how x (or other forest estimates) varies with different samples of Φ Q , which is not straightforward to grasp. In the case of the proposed methods, we are able to overcome these issues. Yet it may be still possible to go beyond with other methods as long as these issues are solved.

In the following sections, we compare the proposed methods with the State-Of-The-Art (SOTA) algorithms in a graph signal denoising setup. Then, we will take a look at possible use cases of the proposed estimators within certain machine learning and optimization applications. In particular, we will look into certain cases where the computation Ky is needed, and the expensive direct computation can be avoided by the forest-based estimators.

Empirical Comparisons with SOTA Algorithms

In this section, we present our empirical comparisons of the forest-based estimators x and z with the state-of-the-art algorithms, namely Chebyshev polynomial approximation (Pol) [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] and the conjugate gradient (CG) method (with and without preconditioning) [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] in approximating the exact solution x. The experimentation setup we use is heavily adapted from our journal paper [START_REF] Yusuf | Graph tikhonov regularization and interpolation via random spanning forests[END_REF].

Experimentation Setup

Datasets. We run our experiments over three benchmark datasets, namely Cora, Citeseer and Pubmed. The following with ∈ N (0, σ 2 ) and the underlying graph. By solving the problem in (3.1), we compute the denoised signal as follows:

x = Ky with K = q(qI + L) -1 (3.16)

We assume the original signal x is a k-bandlimited signal which verifies:

x = k i=1 α i u i
where the u i 's are the lowest k graph Fourier basis and the α i 's are the lowest k graph Fourier coefficients of the signal. In the experiments, we artificially generate such signals by generating α i 's at random. Then, we add Gaussian noise in order to obtain the noisy measurements y. In the whole experimentation, we set k to 50, 50 and 500 for Cora, Citeseer and Pubmed datasets, respectively. The noise variance σ is adjusted per signal such that the signal-to-noise ratio (SNR) of the noisy signal equals to 2.

Parameter tuning. Over all graphs and graph signal realizations, we set the parameter q to the value that gives the best denoising performance (that is, to the value of q minimizing ||x -x|| 2 ).

Performance Metrics. We compare all the algorithms by looking at two performance metrics. Given the solution by an algorithm x , we look at the approximation error ||x -x || 2 and the reconstruction error ||x -x || 2 . While the former measures the quality of the approximation to x, the latter gives the denoising performance of the algorithm.

Experimentation Procedure. The procedure we follow in the experiments is threefold:

• Generate a noisy signal y,

• Tune the hyperparameter q by a grid-search,

• Run the algorithms and measure the run-time, approximation and reconstruction error.

Estimating x via RSFs

We note that all the experiments are implemented in Julia programming language and run in a single thread of a laptop.

Results

As aforementioned, the algorithms in comparisons are: the RSF based estimators x and z, Chebyshev polynomial approximation and conjugate gradient1 method with and without preconditioning. In the preconditioned case, we use the Algebraic Multigrid (AMG) algorithm2 as the preconditioner. In all of these algorithms, there is an iteration parameter. This parameter typically adjust the trade-off between the run-time and the approximation error of the algorithm. These parameters are:

1. The number of forests for x and z 2. The polynomial's degree in Chebyshev approximation 3. The number of iterations for CG and preconditioned CG (PCG).

During the experiments, we vary these parameters in order to observe the performance of the algorithms over changing run-time. For each algorithm and each graph, it is swept through 17 logarithmically spaced values between 1 and 100: {1, 2, . . . , 62, 78, 100}, corresponding to the 17 plotted markers forming each curve.

In Fig. 3.9, we present our results over 3 benchmark graph datasets. From the approximation error results, one can see that the deterministic methods CG and Chebyshev polynomial approximation converges very quickly to the exact solution whereas the forest estimators are stuck with a Monte Carlo convergence rate (which is linear in log-scale). Therefore, the RSF-based methods cannot compete with the deterministic methods for obtaining a low approximation error. On the other hand, from the results on the reconstruction error, we observe that the forest-based estimators are comparable with the other methods. In other words, it is not necessary to have a good approximation in order to get a good reconstruction of the signal. This is because there is no point to compute an inaccurate quantity (in this case, it is x) with a high accuracy (See [START_REF] Bottou | Stochastic gradient descent tricks[END_REF]). In addition, we can see by both performance metrics (more clearly in reconstruction error), z is superior to x. The control variate estimator z achieves better approximation and reconstruction without taking too much additional time.

Pubmed Graph

Fig. 3.9.: Approximation and reconstruction error of the algorithms CG (blue and cyan), polynomial approximation (orange) and the forest estimators x (light green) and z (dark green) for solving graph signal denoising problem. The dark dashed line indicates the time taken by the direct solution (backslash operator in Julia

). In the reconstruction error plots, we also add the error of the initial measurements y in green dashed line and the error of the exact solution x in red line.

Several Use Cases of the RSF-based Estimators

The scope of Eq. (3.1) is not limited to the signal denoising problem. We already know that solving this optimization problem under certain constraints corresponds to solving the Dirichlet boundary problem and semi-supervised learning on graphs. On top of that, we encounter this formulation in solving some graph-related optimization problems. In these cases, solving these problems requires repeatedly solving Prob. (3.1.1) which leaves us with the expensive computation of x at every step. In the following, we show how forest-based estimators can be used in order to circumvent this issue. We conclude each section by giving some illustrations of the forest estimators.

Semi-supervised Learning on Graphs

We consider the node classification problem on graphs. Given a few label information A few labels given on a graph:

We seek the others: over vertices, the main goal is to correctly infer the labels for the other nodes by using the graph structure. Among many possible formulations and solutions, a well-known baseline algorithm is due to the semi-supervised learning algorithm proposed in [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] and later unified under Prob. 3.1.1 by [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. A formal definition of the problem and the algorithm is as follows: Problem 3.4.1 (Semi-supervised learning for node classification). Let us consider a node classification problem over K classes. We denote the vertices with the known labels by V ⊂ V. In a node classification problem, we typically have |V | n as we have a few labeled nodes and define a label encoding as follows:

∀i ∈ V, k ∈ {1, . . . , K}, Y i,k =    1 , if node i belongs to class k 0 , otherwise
Given the matrix Y ∈ R n×K , the main goal is to infer labels for the vertices in V \ V . In a very classical setup, we evaluate a function F : V × {1, . . . , K} → R (also called classification function) such that for each node i ∈ V \ V we assign its label by computing arg max k∈{1,...,K} F (i, k). A baseline method for computing all evaluations F = [F (i, k)] i,k boil down a particular case of Prob. 3.1.1:

F = argmin F∈R n×K ||Y -F|| 2 Q + tr(F LF). ( 3
F = (Q + L) -1 QY. (3.18)
Again, this leaves us with the expensive computation of Ky for computing each column of F.

Adapting the forest-based estimators for approximating all columns F is easy; one only needs to sample Φ Q by passing Q = diag(Q) and compute x (or z) per class k by passing y = Y :|k . As we can re-use the same forests to estimate each column F, we avoid the cumbersome of re-sampling.

Different configuration for Q yields a different type of regularization. For example, for Q = qD, the solution becomes F = q(qI + L rw ) -1 Y where one obtains a regularization with the random walk Laplacian L rw = D -1 L. On top of this, one can obtain a regularization with the normalized graph Laplacian L n = D -1/2 LD -1/2 by slight modifications on the formulation [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. All of these cases are still straightforward to adapt for the RSF estimators. We find a particular case of this formulation in which the adaptation of the forest estimators might not be obvious. This case is due to the celebrated algorithm for semi-supervised learning called label propagation [START_REF] Zhu | Semi-supervised learning with graphs[END_REF]. Let us detail this algorithm by showing how it is linked to the formulation in Eq. (3.17). Then, we show how to use the RSF-based estimators to approximate the solution given by this algorithm.

As described in [START_REF] Zhu | Semi-supervised learning with graphs[END_REF], label propagation consists of the following steps:

1. Set k ← 1. 2. Set F (k) ← D -1 WY and increment k, 3. Reassign F (k) V |: ← Y V |: ,
4. Go to step 2 unless F (k) satisfies a convergence criteria e.g. ||F (k) -F (k-1) || ≤ for some > 0.

This iterative algorithm propagates the known labels through the neighborhoods at every step. In doing so, it uses the random walk operator P = D -1 W to diffuse

Several Use Cases of the RSF-based Estimators

the information encoded in Y. As a result it converges to a closed form solution as follows:

Fi,k =      Y i,k if i ∈ V -(L U |U ) -1 L U |V Y V |: i,k otherwise , (3.19)
where U = V \ V . This solution, in fact, is a limited case of the solution of Eq. (3.17). We show this case in the following proposition:

Proposition 3.4.2 (Label Propagation). For the solution in Eq. (3.18), consider the configuration of the hyperparameters:

Q i,i =    q if i ∈ V 0 otherwise
Then, the limit lim q→∞ F converges to the solution given in Eq. (3.19).

Proof. Let I V ∈ R n×n be a matrix with (I V ) i,i = 1 whenever i ∈ V and zero elsewhere. For the given parameter setting, the matrix K takes a form of:

K = (Q + L) -1 Q = (qI V + L) -1 qI V = (I V + q -1 L) -1 I V .
We re-write this inverse in block matrix form (without loss of generality, the nodes in V are enumerated 1, . . . , |V |):

K = q -1 L V |V + I q -1 L V |U q -1 L U |V q -1 L U |U -1 I V .
The non-zero blocks of this matrix can be computed by block matrix inversion formula:

K = A 0 -(L U |U ) -1 L U |V A 0 . where A = (q -1 L V |V + I -q -1 L V |u (L u|u ) -1 L u|V ) -1 .
As α goes to infinity, one has

lim q→∞ K = lim q→∞ A 0 -(L U |U ) -1 L U |V A 0 = I 0 -(L U |U ) -1 L U |V I 0 .
This proposition in fact gives us a way to adapt the RSF estimators for approximating the solution in Eq. (3.19). We sample Φ Q with the given parameter setting. By definition, we only have non-zero q values over the vertices in V . As q → ∞, the loop-erased random walks in Wilson's algorithm are always interrupted whenever they reach a node in V . Therefore, all sampled forests are conditioned to be rooted in V . Moreover, x converges to x, as the roots in Φ Q given the partitions are no longer random as we have a predefined root set V .

We finish this section by illustrating the forest estimates in solving graph based SSL over some benchmark datasets, namely Cora, Citeseer and Pubmed. These datasets were previously introduced in Table 3.1

In the simulations, we consider a setup where k n vertices per class are labeled and the goal is to classify other vertices given the labeled ones and the underlying graph. In doing so, we use two solutions, namely label propagation (LP) [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] and generalized SSL (gSSL) [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. In both solutions, one needs to calculate the classification function F by computing x per each class under different parametrizations. As we previously discussed, one can estimate both cases by using the forest-based estimator x without taking expensive matrix inverses. In the following illustration, we compare the classification accuracy yielded by x with that of x. In doing so, we use the following procedure:

• Select k n vertices at random per class as the labeled nodes,

• Compute the classification functions F LP and F gSSL given LP and gSSL by using x and x per class.

• For each vertex i, assign arg max c Fi, c as its class and calculate the classification accuracy as the ratio of correctly predicted labels to the total number of predictions.

In Fig. 3.10, we plot the classification accuracy of LP, gSSL and their forest estimation while varying the number of labeled vertices per class k. We average the results over 50 realizations. The empirical results show that the forest estimator for gSSL gives close estimates for the exact solution in Citeseer and Cora whereas it fails to provide the same performance in Pubmed. This is probably due to insufficient number of labeled vertices or poor hyper-parameter choice for the forest estimators. Moreover, one can also deduce x need much less forest realizations to reach the exact solution of LP rather than the generalized SSL. However, sampling a forest for LP might take more time compared to the time needed for gSSL. For example, in the Pubmed graph, for k = 20, sampling a single forest for LP (resp. the generalized SSL) takes 6.3 × 10 -2 (resp. 1.4 × 10 -3 ) seconds averaged over 100 repetitions in a single threaded run time of a laptop. In each plot, the performance of the exact solution is given in the black line. For gSSL, we compute x over 50 (light green) and 500 (dark green) forest realization. In the case of LP, the number of realization becomes 5 (light green) and 50 (dark green). We also add the performance of a classifier that returns a label at random (random classifier) and returns the same label all the time (constant classifier).

86

Chapter 3 Graph Tikhonov Regularization and Interpolation with RSFs

RSF-based Quasi Newton's Method

The graph Laplacian regularization is not restricted to the least-squares problem. In many other problems on graphs, one minimizes the loss function:

L(x) = f y,q (x) + x Lx
to find the solution:

x = argmin x∈R n L(x).
In case of LS i.e. f y,q (x) = q||x -y|| 2 2 , this problem has a closed-form solution. However, this is not the case for many other formulations such as f y,q (x) = q||x-y|| 1 or cross-entropy loss function. In such cases, the optimization literature is capable of serving a large diversity of algorithms to approximate the solution. For convex and differentiable f y,q (x), which we are interested in this thesis, the popular approach is the gradient descent algorithms which are generically introduced in 2.3.3. As long as we have access to the gradient ∇f y,q (x), one can approach the solution with a convergence rate. If the loss function is also twice-differentiable, a faster approach called Newton's method (sometimes called Newton-Raphson method) brings the second derivatives into the game and derives the following iteration [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]:

x k+1 = x k -H -1 ∇L(x) (3.20)
where

H = ∂ 2 L(x) ∂x i ∂x j i,j
is the Hessian matrix which contains the partial second derivatives of the loss function. Let us assume that for i = j, fy,q(x) ∂x i ∂x j = 0. This implies that only prior correlation between the estimated parameters are induced by the graph. Also, by the convexity one has ∂ 2 fy,q(x) ∂x 2 i ≥ 0. Then the Hessian matrix boils down to a positive semi-definite matrix H = L + diag(g) where

g i = ∂ 2 fy,q(x) ∂x 2 i for all i ∈ V.
Newton's method yields much faster convergence than the gradient descent. Intuitively it takes a more direct route to the solution by using extra information on the curvature. However, this extra information comes with a cost. Mainly, computing the inverse of the Hessian matrix is expensive when n is very large. To avoid this expensive operation, H -1 is often replaced with its approximations. This branch of algorithms is called quasi-Newton methods. SR1 (Symmetric rankone) [START_REF] Andrew R Conn | Convergence of quasi-Newton matrices generated by the symmetric rank one update[END_REF], BFGS (Broyden-Fletcher-Goldfarb-Shanno) [START_REF] Fletcher | A new approach to variable metric algorithms[END_REF] and DFS (Davidon-Fletcher-Powell) [START_REF] William C Davidon | Variable metric method for minimization[END_REF] are some of the well-known algorithms that approximate H -1 deterministically. As these algorithms are cheap approximations, their convergence rate is naturally slower (super-linear) compared to the rate of Newton's 3.4 Several Use Cases of the RSF-based Estimators method (quadratic). However, in large dimensions, they may reach the solution more quickly.

Interestingly, the forest-based estimators also can be used for approximating the update in Newton's iteration. In particular, at every iteration k, we first sample RSFs with the parameter set Q = {g 1 , . . . , g n }. Then we evaluate x over these forests with the measurements y = ∇L(x k ) i g i i∈V

.

It is straightforward to see that this setup x gives an unbiased estimate of the update at every iteration. However, as the updates are stochastic, the convergence might not be guaranteed. Nevertheless, in practice, it is possible to guide the iterations towards the solution by adjusting the step size.

In this case, we consider the following iteration:

x k+1 = x k -β x,
where β ∈ R n . Coupling with line search approaches [START_REF] Boyd | Convex optimization[END_REF], it is possible to keep the iterations within safe ranges. We illustrate this technique on an image denoising application where we consider a Poisson prior for f y,q (x).

Removing Shot Noise from an Image

Fig. 3.11.: The original (left) and noisy image (right). The peak signal-tonoise ratio for the noisy image is 17.13

In image de-noising, it is important to know or assume the probabilistic noise model correctly. A big chunk of de-noising algorithms [START_REF] Luisier | Image denoising in mixed Poisson-Gaussian noise[END_REF] assumes that the noise is generated by a Gaussian or Poisson process. While the former (Gaussian) assumption can be easily motivated by the celebrated central limit theorem, the latter (Poisson) is due to a physical phenomenon that is more specific to the image processing domain. An imaging device is an array of sensors where each sensor measures the intensity of the light by counting the photons hitting the sensor. This count turns out to be random at a constant level of light. Due to this randomness, the measurements collected from the sensors are noisy. This type of noise is called shot or Poisson noise (See Fig. 3.11). Further derivations deduce that the Poisson distribution is a natural model for the random number of photons hitting each sensor, hence the name.

We use the graph Laplacian regularization in this illustration to remove the shot noise. Here we consider a 2D-grid graph in which each sensor is a node connected to the sensor on its left, right, top and bottom. Let us denote the random number of photons counted by the sensor i by Y i .

Then we assume:

P(Y i = y i |λ = exp(x i )) = exp(x i ) y i exp(-exp(x i )) y i ! ,
where λ is the parameter of the Poisson distribution. Given this prior, we adapt our fidelity term with the following log-likelihood function:

f y,q (x) = -q i∈V log P(Y i = y i |λ = exp(x i )) = i∈V -y i x i + exp(x i ) + log y i !.
Plugging this leaves us with the following optimization problem:

x = argmin x∈R n q   i∈V -y i x i + exp(x i )   + 1 2 x Lx.
As this problem does not have a closed-form solution and the loss function is convex and twice-differentiable, Newton's method is suitable for approximating the solution. Let us quickly write the gradient and the Hessian for this loss function:

∀i ∈ V, ∇L(x) i = (Lx) i + q(exp(x i ) -y i ) and H = L + q diag g,

Several Use Cases of the RSF-based Estimators

where ∀i ∈ V, g i = exp(x i ) ≥ 0. Then the forest-based estimation of the Newton's update boils down evaluating x under the parameter setting Q = {q exp(x 1 ), , . . . , q exp(x n )} and y = ∇L(x k ) i qg i i∈V . In Fig. 3.12, we give an illustration of this setup.

L-1 Regularization on Graphs

The graph-based regularization can be found in various forms. A generalization can be made by the p-Laplacian operator:

r G,p (x) =   (i,j)∈E w(i, j) p/2 |x i -x j | p   .
In case p = 2, one recovers r G,2 (x) = x Lx. In this section, we focus on another popular choice, p = 1. When applied, this type of regularization, sometimes called the edge LASSO (Least Absolute Shrinkage and Selection Operator), minimizes the variation of the solution over the graph in the LASSO sense. In turn, the solution usually forms a piece-wise constant signal over the vertices. Let us rewrite the optimization problem by combining it with this new regularization:

x = argmin x∈R n q 2 ||x -y|| 2 2 + ||Bx|| 1 . (3.21)
This formulation yields a convex optimization problem. However, it does not have a closed-form solution, unlike the case of p = 2, as the corresponding loss function is not differentiable everywhere. Yet, there are various algorithms to compute the minimizer x . The most straightforward approaches, such as ISTA, use the proximal gradient methods. More advanced methods such as ADMM and IRLS provide faster convergence to x by some extra operations. In the case of ADMM and IRLS, this operation boils down to inverting a matrix in the form of the regularized Laplacian at every iteration. This section shows how to adapt the forest estimators for approximating this expensive inversion. We restrict ourselves to the case of ADMM. However, one can find the same adaptation for IRLS in our journal paper [START_REF] Yusuf | Graph tikhonov regularization and interpolation via random spanning forests[END_REF].

90

Chapter 3 Graph Tikhonov Regularization and Interpolation with RSFs The ADMM algorithm for the problem in Eq. (3.21) takes the following three steps iteration scheme [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]:

x k+1 = argmin x∈R n q 2 ||x -y|| 2 2 + ρ 2 ||Bx -z k + u k || 2 2 z k+1 = argmin z∈R m ||z|| 1 + ρ 2 ||Bx k+1 -z + u k || 2 2 u k+1 = u k + (Bx k+1 -z k+1 ).
The parameters x k , y k and u k are arbitrarily initialized and are updated at every step, and ρ is a user-defined parameter. The most time-consuming part in this iteration is the first step, where we need to compute the inverse:

x k+1 = (qI + ρL) -1 (qy + ρB z k -ρB u k ).
Instead of directly computing this inverse at every step, one can save a significant amount of time by replacing them with forest-based estimates. Moreover, unlike in the case of Newton's method, we do not have to re-sample RSFs at every step. After generating N samples of Φ q/ρ , we can quickly re-evaluate the forest estimates at every step k by simply updating the input signal.

We give a proof of concept for this forest-based ADMM algorithm by a simple example in Fig. 3.13. In this illustration, we consider the image denoising problem with Gaussian noise. In order to denoise the given measurements in Fig. 3.13b, we solve Eq. (3.21) in which the underlying graph is 128 × 128 2-D grid. In doing so, we use the ADMM algorithm with two types of updates for calculating x k+1 at every iteration; these are the exact computation and the forest-based estimation via x. In qualitative results (Figs. 3.13c and 3.13d), we observe that the final result is almost the same for both types of updates. By monitoring the loss function through the iterations, we also observe that the forest-based updates, for a few realizations of forests (N = 10), closely follow the exact updates in minimizing the objective function.

Conclusion

Network structured data often come with some signal/information over the vertices, called graph signals. These kinds of signals are subject to noise or incompleteness. Many classical signal processing tools have been adapted for graph signals in dealing with these issues. A well-known technique gives a unified solution to these two problems by solving the optimization problem 3.1.1. This problem admits a linear closed-form solution with y. Moreover, it appears within several graph-related contexts, namely graph signal filtering and Dirichlet boundary problem. In the former, we see that a particular case of this solution, i.e. Q = qI and y ∈ R n , corresponds to filtering a graph signal with the low-pass transfer function g(λ) = q q+λ . In the latter, we consider discrete partial differential equations defined on graphs that model various physical phenomena. We show that the solution of these equations is again a particular case (Q = lim q→∞ I V ) of the solution of Prob. 3.1.1. Indeed this one-line solution is of central importance in various fields. However, directly computing it might not be straightforward as n grows. The time complexity for computing the inverse is O(n 3 ), which is prohibitive whenever n is around thousands or millions. Having state-of-the-art as the approximate methods in large n, our contributions in this chapter are in avoiding the expensive computations in estimating x. Let us finish by summarizing them:

• We propose an unbiased estimator x for x based on Φ Q ,

• The expected time complexity of this algorithm is found as tr(Q+L) -1 (Q+D) ≤ 2m q +n which remains comparable with SOTA algorithms in terms of the number operations at every sample (resp. every iteration for CGD),

• We show how to adapt two variance reduction techniques, namely conditional Monte Carlo and control variate, to significantly improve the performance of x (in approximating x), yielding two new estimators x and z,

• We show three different real-life use cases in which these estimators can accelerate the complicated intermediate steps.

• The proposed methods are easy to implement and they have a low memory footprint i.e. O(n). Moreover, they are easy to parallelize as we can start random walks from arbitrary nodes for sampling RSFs.

Given these fruitful results and diverse applications, in the next chapter, we shift our focus on the estimation of other useful quantities, including tr(K) or effective resistances. In turn, we find forest-based algorithms and significant improvements on the existing ones.

Conclusion

RSF Estimation of Some Important Graph Quantities

4

" All exact science is dominated by the idea of approximation.

-Bertrand Russell

The forest estimators for x are not the only outcome of our hike through random spanning forests. Estimating the trace tr(K), the effective resistances and other graph filters than q q+λ are three other objectives which we can achieve by using RSFs. In a higher abstraction, the first problem can be considered as the estimation of the trace of an inverted matrix where the inverse is expensive to take. Various applications require such an estimation. Closer to our work, the estimation of tr(K) is of central importance in the hyper-parameter tuning in Prob. 3.1.1 i.e. tuning of q. Similarly, effective resistances have a broad set of applications [SS11; SXGRS20; AALG17; WPKV14]. Initially studied in electrical networks, effective resistances define a measure over the vertices, which has many theoretical links with random walks and USTs. However, directly computing them might be computationally demanding as it requires access to the entries of the pseudo inverse L † . Finally, we show that the graph filters that can be obtained by RSFs are not limited to a single transfer function q q+λ . In fact, one can reach a wide family of filters via RSFs by slight modifications. We dedicate this chapter to presenting our novel RSF algorithms for estimating tr(K) and the effective resistances and we show how to leverage RSF estimators for other types of graph filters.

x = argmin x∈R n q||x -y|| 2 2 + x Lx.
Recall that q is a user-defined parameter. However, the good value of q in the sense of denoising performance cannot be explicitly known. Fortunately, there are several ways, such as AIC, BIC, cross-validation methods and SURE [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF], to find a near-optimal value for q. All of these algorithms either need to compute the diagonal entries K i,i 's, or the trace tr(K) where K = q(L + qI) -1 . However, in both cases, one needs to calculate K for different values of q. As aforementioned, the direct computation comes with a prohibitive cost for large n. This leaves us with approximate methods. The state-of-the-art algorithm for approximating the trace is Hutchinson's estimator [START_REF] Michael F Hutchinson | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF]:

Definition 4.1.1 (Hutchinson's estimator). Let p ∈ R n be a random vector with each element p i ∈ {-1, 1} and independently distributed as:

P(p i = ±1) = 1 2 .
Given N samples from p, Hutchinson's estimator for the trace of K is defined as follows:

h := 1 N N k=1 p (k) Kp (k) .

Estimation of the Trace of the Regularized Inverse of the Laplacian

It is easy to show that this estimator is unbiased with the permutation property of the trace and E[pp ] = I:

E[h] = 1 N N k=1 E[p Kp] = 1 N N k=1 tr(E[pp ]K) = tr(K).
Moreover, the variance is tractable:

Var(h) = 2 N tr(K 2 ) - n i=1 K 2 i,i = 2 N n i =j K 2 i,j .
In other words, the variance is equal to the squared sum of the off-diagonal entries in K.

In some other methods, such as Girard's estimator, the distribution of p differs, but they still provide an unbiased estimation. In fact, as long as the random vector p satisfies E[p] = 0 and E[pp ] = I, one obtains an unbiased estimator by the averaging of quadratic terms. However, Hutchinson's method is the estimator that provides the minimum variance [START_REF] Michael F Hutchinson | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF]. Even so, it is still subject to two significant bottlenecks:

• Computation of each quadratic term p (k) Kp (k) requires computing Kp (k) which requires inverting 1 q (L + qI),

• As it is a Monte Carlo estimator, it has the convergence rate of O(N -1/2 ), thus one might need many samples to reach an very accurate result.

Taking the inverse is not practical with the exact methods, such as Gaussian elimination. Therefore, the solution is often calculated via iterative methods such as gradient descent algorithms. In dealing with the second issue, one can find many generic variance reduction methods in the literature [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF]. However, it is not obvious how to adapt these methods for Hutchinson's estimator to obtain considerable performance gain.

Previously in [START_REF] Barthelme | Estimating the inverse trace using random forests on graphs[END_REF], the authors (my supervisors) proposed a simple estimator for tr(K) with strong theoretical properties:

s := |ρ(Φ q )|.
The expectation and variance of this estimator are already noted as:

E[s] = tr(K) and Var(s) = tr(K -K 2 ).

96

Chapter 4 RSF Estimation of Some Important Graph Quantities Moreover, they give an empirical comparison with Girard's estimator coupled with several iterative methods for computing Kp. The experimentation over various graphs shows that the forest estimator often gives the best or at least comparable performance with state-of-the-art. Inspired by these results, we give two efficient ways of reducing the variance of the RSF trace estimator. These ways are well-known generic variance reduction techniques in the Monte Carlo literature. However, their adaptation for a given Monte Carlo estimator is not evident or sometimes not possible. The main contribution of this section shows how to adapt these variance reduction techniques for the forest-based trace estimator [START_REF] Barthelme | Estimating the inverse trace using random forests on graphs[END_REF] yielding two novel estimators for tr(K).

Variance Reduction for Trace Estimation

We adapt two variance reduction techniques, namely the control variate method and stratified sampling, for reducing the variance of |ρ(Φ q )|.

Control Variate for Trace Estimation

In this section, we adapt the control variate technique introduced in Section 3.3.2. Let us briefly recall this technique and how it is applied to the previous estimator for x = Ky. Given the estimator x, we seek another random quantity with a known expectation. In our case, we chose this quantity as α(K -1 xy) and the resulting estimator is z = x -α(K -1 xy). Let us redefine z as a matrix-vector product z := Zy where Z = S -α(K -1 S -I) and S verifies x = Sy. We already know that

E[ Z] = K.
Then, a new unbiased estimator of tr(K) is defined by taking the trace of Z:

s := s -α(tr(K -1 S -I)) = s + α   n -s - 1 q n i=1 j∈N (i) w(i, j) Si,i I(r Φq (i) = j)   , (4.1)
where s = tr( S) = tr( S) = |ρ(Φ q )|. Here the sum over the neighbors of all the vertices accumulates the weighted diagonal entries of Si,i over the edges that are between different parts of the partition induced by Φ q . We illustrate these edges in The edges between trees in a spanning forest (in bold). One needs to track these edges to compute the control variate introduced in Eq. 4.1. For computing the control variate induced by S, one only needs the edges between trees that are incident to the roots. In this example, only edge that fits this description is (4, 6).

see that α = 2q q+2dmax guarantees variance reduction. It is also possible to adapt the formulas of α for estimating the optimal case α . Finally, the heuristic value α = q q+davg performs very well in practice across various graphs.

Under a proper selection of α, this technique can substantially decrease the variance. However, it also requires some additional computations. In particular, the additional computation after sampling a forest for computing s requires traversing the edges of the graph. By replacing S with S, this cost can be decreased at the expense of accuracy because for computing the control variate in Eq. (4.1), one only needs to accumulate some edge weights only around the roots instead of all vertices. In this case, one does not have to compute Si,i 's and only need to traverse the neighbours of the roots. However, this version, denoted by s, is worse than s in approximating tr(K) as S is less accurate S.

Stratified Sampling for Trace Estimation

Another variance reduction technique that is applicable to our case is stratified sampling. In a nutshell, stratified sampling is a divide & conquer strategy for improving Monte Carlo estimation. Let us go back to the generic example of Monte Carlo estimation where we approximate the expectation µ X from the samples of the random variable X. Now consider another random variable Y whose outcome space

Ω(Y ) is divided into K strata (disjoint parts) C 1 , C 2 , . . . , C K such that Ω(Y ) = ∪ K k=1 C k .
We assume that Y verifies for all C k :

• The probabilities P(Y ∈ C k ) are known or easy to compute,

• The conditional random variable X|Y ∈ C k is easy to sample.

98

Chapter 4 RSF Estimation of Some Important Graph Quantities

These constraints might seem a bit technical and hard to satisfy for all cases. However, whenever such a random variable Y is available, stratified sampling can decrease the variance significantly by defining a new estimator as follows:

Z := K k=1   1 N k N k j=1 X j |Y ∈ C k   P(Y ∈ C k )
N k is the number of samples taken from the conditional random variable X|Y ∈ C k .

The stratified estimator Z is also unbiased since:

E[Z] = K k=1   1 N k N k j=1 E[X j |Y ∈ C k ]   P(Y ∈ C k ) = K k=1 E[X|Y ∈ C k ]P(Y ∈ C k ) = E[X] = µ X .
Also, the variance of Z reads:

Var(Z) = K k=1 P(Y ∈ C k ) 2 1 N k Var(X|Y ∈ C k ).
This variance is guaranteed to be less than that of the original estimator for certain configurations of N k 's. For example, if one can collect samples of X|Y ∈ C k such that N k ∝ P(Y ∈ C k ) for all k, the variance of Z becomes:

Var(Z) = K k=1 1 N P(Y ∈ C k ) Var(X|Y ∈ C k )P(Y ∈ C k ) 2 = K k=1 1 N Var(X|Y ∈ C k )P(Y ∈ C k ) = 1 N E Y [Var(X|Y ∈ C k )] ≤ 1 N Var(X).
The last inequality is due to the law of total variance. In fact, the optimal configuration for variance reduction is:

∀k, N k ∝ Var(X|Y ∈ C k )P(Y ∈ C k )
However, this is not practical in most cases as the conditional variance Var(X|Y ∈ C k ) is often not known explicitly. In the case of forest-based trace estimation, we find a novel way to apply stratified sampling to reduce the variance. This adaptation yields a new estimator for tr(K) which often provides the best empirical performance.

S =

S (1) S (2) S (3) S (4) S (5) S (6) S (7) In order to explain our stratification scheme for trace estimation, we need to recall some aspects of the stack representation of Wilson's algorithm. Recall the infinite stack

          2 7 6 1 Γ Γ 1 Γ 7 2 5 4 4 2 3 7 1 5 1 4 1 . . .
S (i) = [S (i) 1 , S (i) 2 , .
. . ] defined for each node i. Every element in S (i) is an independent random variable with the probability law P(S (i)

j = k) = w(i,k) d i .
In the case of sampling RSFs, we add a new node Γ to this representation and slightly change the probability distribution as:

P(S (i) j = k) =    w(i,j) q+d i , k ∈ N (i), q q+d i , k = Γ.
Note that whenever a node has Γ at the top of its stack, it becomes a root. Now we define a new random set by only looking at the first layer of the stacks:

ρ i (Φ q ) := {v ∈ V : S (v) j = Γ and j ≤ i}.
In other words, the random set ρ i (Φ q ) contains the random roots in Φ q sampled until their i-th stack popped. See Fig. 4.2 for an example at i = 1. In Wilson's algorithm, this set can be seen as the roots sampled up until i-th visit by the random walks. Notice that as i → ∞, we recover the original roots ρ(Φ q ). For stratification, we look at the cardinality of this set at i = 1. As all the stacks are independently sampled at the beginning i.e. i = 1, this cardinality equals a sum of Bernoulli random variables:

|ρ 1 (Φ q )| = v∈V Ber q q + d i ,
where Ber(p) ∈ {0, 1} denotes a Bernoulli random variable that is equal to 1 with probability p. Our main idea here is to use |ρ 1 (Φ q )| as the additional random variable in stratification of the forest trace estimator s. Indeed, to do so, we need to show that this new random variable aligns with the constraints of stratified sampling.

The first constraint is to have probabilities that are easy to compute. In order to verify that this is the case with |ρ 1 (Φ q )|, let us look at its probability distribution. By definition, |ρ 1 (Φ q )| is a sum of independent Bernoulli variables where each does not necessarily have the same parameter. Such random variables, also known as Poisson-Binomial random variables, are one of the classical random variables in probability theory [START_REF] Yuan | On the number of successes in independent trials[END_REF]. As a result, their distribution has been well studied, and efficient direct or approximate numerical computations have already been proposed [START_REF] Hong | On computing the distribution function for the Poisson binomial distribution[END_REF].

In the case of large n, the normal approximations (motivated by the central limit theorem) are able to provide good performance at a very small cost. In this approach, we approximate the distribution of |ρ 1 (Φ q )| with a normal distribution N (µ, σ 2 )

where µ = E[|ρ 1 (Φ q )|] = v∈V q q+d i and σ 2 = Var(|ρ 1 (Φ q )|) = i∈V qd i
(q+d i ) 2 . This approximation highly simplifies any computation related to the probabilities that we need to calculate as the cumulative/probability mass functions are explicitly known. Indeed, one can go beyond this approximation via more sophisticated approaches. However, by doing so, we do not obtain much improvement in practice for trace estimation. This is probably because the Monte Carlo error is dominant by far.

The second constraint dictates that the conditional random variable |ρ(Φ q )| | |ρ 1 (Φ q )| ∈ C k must be easy to sample for each stratum C k . The sampling procedure per k can be considered two-fold:

• Generate a sample S of ρ 1 (Φ q ) that verifies |ρ 1 (Φ q )| ∈ C k ,
• Sample |ρ(Φ q )| given the sample generated at the first step.

The second step is easy by a simple modification of Wilson's algorithm for forests. As S are the roots that are sampled at the first sight, in Wilson's algorithm, we set S as the predefined set of roots, and we do not allow any other to be a root at the first visit of random walks. The forest generated this way is a sample of Φ q |ρ 1 (Φ q ) = S. For the first step i.e. generating S that verifies |S| ∈ C k , the naive approach is the rejection sampling where we keep sampling from ρ 1 (Φ q ) until the constraint |S| ∈ C k is satisfied. This approach may terminate very quickly if P(|ρ

1 (Φ q )| ∈ C k ) 0.
Indeed, this constraint is directly related to how the strata C 1 , . . . , C K is chosen. In our case, we approximate the distribution of |ρ 1 (Φ q )| by a normal distribution. Thus, we can design the strata C 1 , . . . , C K so that P(|ρ

1 (Φ q )| ∈ C k )
0 is verified for all k ∈ {1, . . . , K}. Both in theory and practice, we do not see a particular benefit of having strata with small probabilities for each stratum. Therefore, we choose to use rejection sampling in the application of trace estimation.

Estimation of the Trace of the Regularized Inverse of the Laplacian

Up to this point, we show that |ρ 1 (Φ q )| is a theoretically and practically suitable random quantity for stratification of s. Then let us define the stratified estimator for tr(K) as follows:

s st := K i=1 1 N i      N i j=1 |ρ 1 (φ (j) )|∈C i |ρ(φ (j) )|      P(|ρ 1 (Φ q )| ∈ C i ). (4.2)
This concludes the methods that we propose to reduce the variance of the forestbased trace estimators. In the next section, we give empirical comparisons of these two approaches (control variate and stratified sampling) with the state-of-the-art algorithms over various real-life graphs.

Empirical Results on Trace Estimation

This section is adapted from the results section of our paper [START_REF] Yusuf | Variance Reduction for Inverse Trace Estimation via Random Spanning Forests[END_REF]. In these results, we compare the forest-based trace estimators s, s, s and s st with Hutchinson's estimator combined with several linear solvers.

Firstly, notice that all methods in this comparison are Monte Carlo methods. Recall that the variance of a Monte Carlo estimator over N samples reads:

σ 2 N ≈ σ 2 1 N ,
where σ 1 is the standard deviation over a single sample. 1 We leverage this fact to compare the effective runtimes of all methods i.e. the time needed to reach a fixed relative error . First, we run all methods with N = 100. This gives us the average runtime for the computation per sample and the sample variance σ2 N . Then, we approximate σ1 = √ N σN for each method. By using this approximation, we solve = σ1 tr(K) √ k for = 0.002 to calculate the number of iterations k needed to reach error. Finally, we calculate the effective runtime per method by multiplying k by the average time for generating a single sample.

The linear solvers we use in Hutchinson's estimator are Algebraic Multigrid (AMG)2 , Conjugate Gradient (CG)3 , CG preconditioned with AMG and sparse Cholesky decomposition using CHOLMOD [START_REF] Chen | Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate[END_REF]. Here we use the block implementation of CG given in [START_REF] Dianne | The block conjugate gradient algorithm and related methods[END_REF]. For the control variate methods, we choose α = q q+davg for s and s. In stratified sampling, we divide the sample space into 5 strata C 1 , . . . , C 5 verifying P(|ρ (Φ q )| ∈ C k ) ≈ 0.2 for all k = 1, . . . , 5. Here, we empirically choose the number of strata so that P(|ρ (Φ q )| ∈ C k ) 0. For too many strata, this probability approaches zero and generating samples from each stratum take much more time. For too few, the effect of the stratification diminishes. We use the proportional allocation where we set N k = N P(|ρ (Φ q )| ∈ C k ) per stratum k. The graphs4 we use in this comparison can be listed as:

• K-random regular: A random regular graph with with n = 10 4 nodes and m = 10 5 edges generated by K-random regular model with the parameter k = 20,

• Barabasi-Albert: A random graph with n = 10 4 nodes and m = 99900 edges generated from a Barabasi-Albert model with the parameter k = 10.

• Collab-CM: A collaboration network of n = 21363 nodes and m = 91342 edges. Each node is an author in Arxiv on condense matter physics and each edge depicts their collaboration,

• Citation-HEP: A citation network with n = 34401 nodes and m = 420828 edges. Each node is a paper in Arxiv on high energy physics and there is an edge whenever a paper is cited by another,

• 3D-Grid: A 3 dimensional grid graph with n = 125000 nodes and m = 375000,

• Amazon: A real-life network in Amazon with n = 262111 nodes and m = 899792 edges. Each node corresponds to a product and an edge between two products indicates that both products are bought by the same clients.

Fig. 4.3 summarizes our comparisons on effective runtimes while we vary q so that tr(K)/n varies approximately between 0.01 and 0.65.

In relatively small and sparse graphs, such as Collab-CM, the direct method (solving via sparse Cholesky decomposition) gives the best performance which is closely followed by RSF methods. When it comes to larger and denser graphs, the RSF methods often provide the best performance or a comparable performance with the state-of-the-art algorithms.

Estimators for Effective Resistances

In the second part of this chapter, effective resistances, another significant set of graph quantities, will be the central object. They have important roles in many applications, yet they are not easy to access in the case of large graphs. In the following, after giving their formal definition and, certain properties, we focus on the problem of efficiently computing them. In doing so, we revisit existing algorithms and propose two forest-based estimators. We finish this section with the empirical comparison of these algorithms.

A Crash Course on Effective Resistances

A distance metric, e.g. the geodesic distance (the length of the shortest path between a pair of vertices), over the vertices of a graph describes the distance of each pair of vertices on the graph. Such metrics have a vast number of applications in network science. However, distance metrics may be defined in various ways [Sol11;

RR15], and there is no definitive answer as to which one is the best. Depending on the theory and the application that the metric is deployed (mostly the latter), network analysts may favour one over others. The effective resistances are a distance metric over vertices. Initially studied in electrical networks, the effective resistances have many theoretical links with random walks and USTs [START_REF] Lyons | Probability on Trees and Networks[END_REF]. Therefore, they contain significant information related to the probabilistic and combinatorial properties of the graph. Motivated by these facts and many others, they have been used in many network science applications, including graph clustering [START_REF] Levi Alev | Graph clustering using effective resistance[END_REF], sparsification [START_REF] Daniel | Graph sparsification by effective resistances[END_REF] and learning [START_REF] Rakshith S Srinivasa | Fast graph attention networks using effective resistance based graph sparsification[END_REF]. In this section, we detail their connections with USTs (and also RSFs) and give RSF-based methods for their numerical computation.

Let us start with the formal definition of effective resistances: Definition 4.2.1 (Effective Resistances). A resistor network with n nodes (connection points of resistors) can be depicted by a graph G = (V, E, w) where each edge is a resistor with a resistance value 1/w(i, j). The effective resistance R i,j is the reciprocal of the current induced when unit potential difference/voltage is applied between i and j. It can be seen as the resistance measured between node i and j. One can calculate R i,j in terms of the entries of L † :

R i,j = L † i,i + L † j,j -L † i,j -L † j,i .
The last algebraic identity takes its roots from the combination of Ohm's law and Kirchoff's current law. In order to elaborate, let us denote the voltage at node i by v i and the current at edge (i, j) by c i,j . Applying a unit voltage between i and j means v i = 1 and v j = 0. Then, by the Ohm's law, one has c x,y = w(x, y)(v x -v y ) for any pair (x, y) ∈ E. In a matrix form, this constraint can be summarized as:

CB u v = c. (4.3)
where B u is the unweighted edge incidence matrix and C ∈ R m×m is the diagonal matrix which contains the edge weights in its diagonal entries. In addition, Kirchoff's current law indicates that the total current flow f x = y∈N (x) c x,y at every node x = i, j equals to zero. The current flow f i at node i (or f j at node j) gives us the effective conductance, which is the reciprocal of effective resistance between i and j.

Let us write the flow vector as

f = 1 R i,j (δ i -δ j ) where δ i = [I(i = j)] j∈V ∈ R n . Then one has: f = B u c = B u CB u v = Lv 4.2 Estimators for Effective Resistances

105

Multiplying both sides with f L † , one has:

f L † f = f I - 1 n J v,
where J = 11 ∈ R n×n . Also by Kirchoff's current law f i = -f j as the total current going in and out from the circuit must be the same. Then, one has f 1 = 0 and

f v = 1 R i,j (v i -v j ) = 1 R i,j : 1 R 2 i,j (L † i,i + L † j,j -L † i,j -L † j,i ) = f v = 1 R i,j R i,j = L † i,i + L † j,j -L † i,j -L † j,i
Along with these algebraic descriptions, effective resistances also have probabilistic properties. One well-known property, closely related to our work, connects effective resistances with RSTs:

Theorem 4.2.1 (Edges in RSTs [START_REF] Nisheeth | Lx= b[END_REF]). Given a graph G = (V, E, w), the probability of having an edge (i, j) in the RST, T reads:

P((i, j) ∈ T ) = w(i, j)R i,j .
Proof. The RST T is a projection DPP over E that verifies for all S ⊂ E:

P(S ∈ T ) = det(B(B B) † B) S
Then writing the probability for a single edge e = (i, j) finishes the proof:

P(e ∈ T ) = det(B(B B) † B) e = det(BL † B) e = w(e)(L † i,i +L † j,j -L † i,j -L † j,i ) = w(e)R i,j .
Due to all of these properties and the others omitted in this thesis, effective resistances have been playing a central role in many applications of network science. On the other hand, direct computation is demanding as it requires computing the entries of L † which is expensive when the graph is large. Therefore, the state-of-the-art algorithms are approximate methods that avoid this computation at the cost of precision. In the following sections, we visit some of these algorithms and propose two forest-based estimators for effective resistances.

106

Chapter 4 RSF Estimation of Some Important Graph Quantities

Computing Effective Resistances

In the case of large graphs, the existing algorithms to compute R i,j 's are randomized approximate approaches. Each algorithm deploys different type of randomization in order to leverage probabilistic properties of effective resistances or to optimize a certain aspect of the algorithm. In the following, we list three algorithms and highlight their strengths and weaknesses:

• Estimation by Spanning Trees [START_REF] Hayashi | Efficient Algorithms for Spanning Tree Centrality[END_REF]: By Theorem 4.2.1, we already know the direct link between RSTs and effective resistances. In [START_REF] Hayashi | Efficient Algorithms for Spanning Tree Centrality[END_REF], the authors use this link to give an unbiased estimator for R i,j 's for all (i, j) ∈ E. In a nutshell, after sampling N spanning trees T 1 , . . . , T N by Wilson's algorithm, this estimator computes:

∀(i, j) ∈ E, R st i,j := 1 N N k=1 I((i, j) ∈ T k ) w(i, j) ,
where I is the indicator function. The unbiasedness is easy to check via Theorem 4.2.1. By Hoeffding's inequality, the authors give the number of spanning trees to sample as

N = log(2m/δ) 2 2
to guarantee |R i,j -R st i,j | ≤ for each (i, j) ∈ E with the probability at least 1 -δ. However, this estimator is restricted to estimate effective resistances only over the edges.

• Estimation by Random Projections [SS11]: This method suggests estimating R i,j for all pairs of (i, j) ∈ V by using efficient Laplacian solvers [START_REF] Daniel | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF] and Johnson-Lindenstrauss' lemma [START_REF] William | Extensions of Lipschitz mappings into a Hilbert space[END_REF]. The estimator takes three main steps:

1. Sample a random matrix A ∈ R k×m where k = O( log n / 2 ) for some error parameter > 0 and each

A i,j ∈ {-1/ √ k, -1/ √ k} is a normalized Rademacher random variable. 2. Compute Y = AB.
3. For each row of Y, denoted by y i , solve the linear system

Lz i = y i to compute Z = [z 1 | . . . |z k ].
4. Finally, one can calculate the estimate R rp i,j for any pair (i, j) ∈ V as:

R rp i,j := ||Z(δ i -δ j )|| 2 2 .
By the Johnson-Lindenstrauss lemma, this algorithm verifies:

∀(i, j) ∈ V, (1 -)R i,j ≤ R rp i,j ≤ (1 + )R i,j

Estimators for Effective Resistances

In [START_REF] Daniel | Graph sparsification by effective resistances[END_REF], the authors propose to approximate the inversion at step 3 with the algorithm proposed in [START_REF] Daniel | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF]. This algorithm computes the estimates zi 's such that each verifies ||z i -L † y i || L ≤ δ||L † y i || L for some δ > 0 at the time complexity O(m log c n log( 1 /δ)) for some constant c > 1. Then Z computed with the error parameter δ ≤ 3 2(1-)w min (1+ )n 3 wmax5 verifies:

∀(i, j) ∈ V, (1 -) 2 R i,j ≤ Rrp i,j ≤ (1 + ) 2 R i,j
where Rrp i,j := || Z(δ i -δ j )|| 2 2 . The authors develop on these facts and give the time complexity of computing Z, and so the whole algorithm, as Õ( m log r / 2 ) where r = wmax /w min6 .

• Local Algorithms [START_REF] Peng | Local Algorithms for Estimating Effective Resistance[END_REF]: Closer to our work, the authors in [START_REF] Peng | Local Algorithms for Estimating Effective Resistance[END_REF] leverage the links of effective resistances with random walks, Markov chains and spanning trees. In turn, they propose a family of randomized algorithms for estimating R i,j 's. Given a pair of nodes s and t, these algorithms estimate R s,t without needing the whole graph as the input. They operate locally around the nodes s and t (e.g. they use random walks that start from node s and stop at node t). For the detailed analysis of all of these algorithms, we refer the reader to [START_REF] Peng | Local Algorithms for Estimating Effective Resistance[END_REF]. In the experiments, we implement two of these algorithms (the most precise and the fastest) in order to compare them with the algorithms proposed in this manuscript.

This finishes our revisit on the SOTA algorithms for estimating R i,j 's. We note that all these algorithms are randomized, but they differ in their pros/cons. For example, the estimator R st i,j can approximate R i,j only if (i, j) ∈ E whereas one can reach an arbitrary R i,j with i, j ∈ V via random projections and some of the local algorithms7 . In the following sections, we give the forest-based methods to estimate R i,j 's.

Estimating ER via 2-Rooted Forests

The first algorithm we propose is inspired by the following fact: Lemma 4.2.2 (R i,j as a determinantal ratio [START_REF] Ravindra | Graphs and matrices[END_REF]). Given a graph G = (V, E, w), one has:

∀i, j ∈ V, R i,j = det L -i,-j|-i,-j det L -i|-i = φ∈F {i,j} w(φ) τ ∈T i w(τ )
where T i is the set of all spanning trees rooted in i and F {i,j} is the set of all spanning forests with exactly two disjoint trees rooted in i and j.

Proof. The second part of the identity is straight-forward by adapting Theorems 2.1.2 and 2.1.3. Then we prove this lemma by showing that R i,j equals the ratio of these two determinants. In order to do so, we go back to electrical networks. Assume we apply the unit voltage between the nodes i and j in order to obtain the effective resistance R i,j . Recall that, by Ohm's law, one has the relation c x,y = w(x, y)(v x -v y ) between the currents c x,y and voltages v x , v y at every edge (x, y). Moreover, by Kirchoff's law, the total current flow f x at every node x ∈ V\{i, j} is equal to 0. The current flow f i gives the effective conductance 1 R i,j . A matrix form [START_REF] Ravindra | Graphs and matrices[END_REF] (alternative to Eq. 4.3) to summarize these constraints is:

Av = δ i with A =     L -i,-j|: δ i δ j     ∈ R n×n
where δ i = [I(i = j)] j∈V ∈ R n . By solving this linear system, one computes the voltages induced at every node via Cramer's rule as:

∀k ∈ V \ {i, j}, v k = (-1) i+k det L -i,-j|-k,-j det L -i,-j|-i,-j .
Then rewriting the flow at node i gives us: j|-i,-j .

f i = k∈N (i) w(i, j)(v i -v k ), = k∈N (i) w(i, j) 1 - (-1) i+k det L -i,-j|-k,-j det L -i,-j|-i,-j , = det L -i|-i det L -i,-
Recalling f i = (R i,j ) -1 finishes the proof.

Estimators for Effective Resistances

Corollary 4.2.3 (2-Rooted Forests and ERs). For any pair i, j ∈ V, one has:

P(ρ(Φ q ) = {i, j}||ρ(Φ q )| = 2) ∝ R i,j .
Proof. The proof is easy by multiplying and dividing the determinantal ratio above by q 2 φ∈F 2 w(φ):

R i,j = φ∈F 2 w(φ) τ ∈T i w(τ ) q 2 φ∈F {i,j} w(φ) q 2 φ∈F 2 w(φ) = φ∈F 2 w(φ) τ ∈T i w(τ ) P(ρ(Φ q ) = {i, j}||ρ(Φ q )| = 2).
The proportionality follows by realizing that φ∈F 2 w(φ)

τ ∈T i w(τ ) is invariant w.r.t. i and j a . a Recall that ∀x, y ∈ V, τ ∈Tx w(τ ) = τ ∈Ty w(τ )
Let us define the RSFs conditioned over |ρ(Φ q )| as Φ (k) q := Φ q ||ρ(Φ q )| = k. and consider a random variable X i,j = I(ρ(Φ (2) q ) = {i, j}). By the corollary above, we have:

E[X i,j ] = P(ρ(Φ (2) q ) = {i, j}) ∝ R i,j .
Then, one can estimate R i,j with X i,j up to a constant. Let X = [X i,j ] i,j∈V be the matrix form of this estimator. By definition, only non-zero entries of X are X i,j = X j,i = 1 and the rest is equal to 0. In fact, this estimation can be improved by deploying the conditional Monte Carlo technique used in 3.3.1. Recall Prop. 2.4.6 which states that the distribution of roots is uniform Φ q when the partition is fixed.

Then an improved version of X with the conditional Monte Carlo method is:

∀i, j ∈ V, Xi,j :=      1 |V t(π(Φ (2) q ,i) ||V t(π(Φ (2) q ,j) | if r Φ (2) q (i) = r Φ (2) q (j) 0 otherwise . (4.4)
The estimator X admits the same expectation as X by the law of iterated expectation,

The matrix X for a given spanning forest with two trees.

and it gives a better estimation than X by the law of total variance. Alternatively, for a 2-rooted spanning forests with two parts V 1 and V 2 , one can write X a low-rank representation as:

X = 1 |V 1 ||V 2 | 1 V 1 1 V 2 + 1 V 2 1 V 1 , (4.5) 110
Chapter 4 RSF Estimation of Some Important Graph Quantities

where

1 V 1 = [I(i ∈ V 1 )] i and 1 V 2 = [I(i ∈ V 2 )] i .
On the other hand, Xi,j can estimate R i,j up to a graph related constant c, i.e. E[ X] = cR where,

c = τ ∈T i w(τ ) φ∈F 2 w(φ)
.

We already know that this constant is invariant in i and j. In fact, c here can be written in terms of the effective resistances as follows:

c = τ ∈T i w(τ ) φ∈F 2 w(φ) =   i,j∈V φ∈F i,j w(φ) τ ∈T i w(τ )   -1 =   i,j∈V R i,j   -1 .
Indeed, this quantity is solely related to the graph itself. In fact, it is the reciprocal of the Kirchoff's index, that is often used to measure the stability of a network [START_REF] Ellens | Graph measures and network robustness[END_REF]. However, it is not easy to compute directly. In order to circumvent this issue, we propose to use the following identity (Foster's theorem [START_REF] Ronald | The average impedance of an electrical network[END_REF]):

(i,j)∈E w(i, j)R i,j = n -1.
As this identity holds for all graphs, an unbiased estimator for c is:

C = (i,j)∈E w(i, j) Xi,j n -1
Then we define a ratio estimator as follows:

R 2rf i,j := Xi,j C .
This estimator, like many other ratio estimators, is biased in most cases.

Bias analysis. The bias can be approximated by Taylor expansion on E[R 2rf i,j ] (at second order of derivatives):

E[R 2rf i,j ] -R i,j ≈ - Cov( Xi,j , C) E[C] 2 + E[ Xi,j ] E[C] 3 Var(C).
This bias usually diminishes at the rate of O(N -1 ) as N grows8 . Moreover, one can correct the bias by estimating it from the samples by the sample covariance and variance formulas and subtracting it from the estimate. Many recent works [START_REF] Jose | Unbiased Monte Carlo computation of smooth functions of expectations via taylor expansions[END_REF][START_REF] Babu Moka | Unbiased estimation of the reciprocal mean for non-negative random variables[END_REF] propose various ways of improvement (including unbiased estimators), taking extra computational time.

Variance analysis. The variance of the estimator is the other main (usually the most important) element of the approximation error of the proposed estimator. Similar to the bias, it can be approximated by Taylor expansion as follows:

Var(R 2rf i,j ) ≈ Var( Xi,j ) E[C] 2 -2 E[ Xi,j ] E[C] 3 Cov( Xi,j , C) + E[ Xi,j ] 2 E[C] 4 Var(C) = 1 c 2 (Var( Xi,j ) -2R i,j Cov( Xi,j , C) + R 2 i,j Var(C)) (4.6)
In order to better understand this variance expression, we derive an upper-bound in terms of graph-related constants such as n, m, c or R i,j . In doing so, instead of an individual pair (i, j), we consider the sum of the variance over all pairs i,j∈V Var(R 2rf i,j ) to get a global sense. As implied in Eq. 4.6, i,j∈V Var(R 2rf i,j ) consists of three additive terms. For each, we find the following upper-bounds (not necessarily tight) for unweighted graphs:

• i,j∈V Var( Xi,j ) ≤ 2 n-1 -c 2 i,j R 2 i,j , • -2 i,j∈V R i,j Cov( Xi,j , C) ≤ 2c 2 i,j R 2 i,j , • i,j∈V R 2 i,j Var(C) ≤ m 2 (n-1) 4 -c 2 i,j∈V R 2 i,j .
The derivations can be found in App. A.5. Putting all of these together yields that the total variance is bounded:

i,j∈V Var(R 2rf i,j ) 1 c 2   2 n -1 + m 2 (n -1) 4 i,j∈V R 2 i,j   .
This bound indicates that for larger c =

1 i,j∈V R i,j
, the total variance tends to be smaller. This indicates that R 2rf tends to perform better when the effective resistances are small all over the graph.

Although R 2rf estimates all pairs over a single forest, in practice it yields a poor performance compared to other algorithms. In the following, we list several reasons and possible solutions to improve the performance by going beyond the work done in this thesis:

• Recycling forests with 1,3,4 trees. In order to sample X and C, one needs to sample spanning forests with exactly two trees for a fixed value of q. However, for an arbitrary value of q, sampling a two-rooted forest might be a rare event.

In order to find the proper value of q, [ACGM18] suggests using an annealing process. This process iteratively decreases the value of q until the desired number of roots is approximately reached. However, even plugging this value, one might still end up with many forests with 1,3,4 trees which cannot be used in X and C. The up-down sampler given in [START_REF] Luıés | Linking and cutting spanning trees[END_REF] may be a possible fix for this issue. Originally designed for sampling uniform spanning trees, the up-down sampler for RSFs takes a Markov chain whose states are the spanning forests of the graph (or the spanning trees on the extended graph with Γ), and stationary distribution is the distribution of Φ q . There are three types of transitions authorized between the states/forests; adding and deleting an edge, adding an edge and deleting a root and finally adding a root and deleting an edge. Setting the transition probabilities accordingly, one can guarantee the stationary distribution is the same with the distribution of Φ q . An MCMC (Markov Chain Monte Carlo) algorithm to sample from the distribution of Φ q starts from an arbitrary state in this Markov chain and run a random walk.

After taking a sufficient number of steps, the random walk is stopped and the current state is distributed similar to Φ q . An alternative use of this random walk may allow us to transform the spanning forests with 1,3,4 roots to the spanning forests 2 roots. Instead of an arbitrary state, we start with a spanning forest with 1,3,4 roots sampled via Wilson's algorithm and do the random walk until reaching a spanning forests with two trees. As we sample our initial state from the stationary distribution, the final step we reach by the random walks is also distributed according the stationary distribution. In this way, one recycles the unused samples in computing X and C. One of the main difficulties with this idea is that one might need to take too many steps until reaching a tworooted forest depending on the graph. Also, one needs to implement complex data structures to maintain transitions efficiently as proposed in [START_REF] Luıés | Linking and cutting spanning trees[END_REF].

• Variance Reduction. High variance is the other main issue with this estimator. One may attempt to solve this issue by variance reduction techniques with the help of rich machinery of statistics in RSFs. A particular example that we come up with is the adaptation of the control variate method. In this case, we use the following relation [START_REF] Gutman | Generalized inverse of the Laplacian matrix and some applications[END_REF] between L † and R:

L † = - 1 2 R - 1 n (RJ + JR) + 1 n 2 JRJ . (4.7)
Multiplying both sides by L, one has:

I - 1 n J = - 1 2 LR I - 1 n J .
By plugging R 2rf instead of R, one can show that the expectation E LR 2rf I -1 n J converges to a known matrix 2 1 n J -I . Thus one can use LR 2rf I -1 n J as 4.2 Estimators for Effective Resistances a control variate. Moreover, by the low-rank representation of X in Eq. (4.5), computing the matrix product LR 2rf boils down to cheap matrix vector products.

At this moment, unfortunately, R 2rf yields poor estimator. Besides its simplicity, its main advantage is that given a spanning forest, it provides estimates for all pairs of effective resistances. Nevertheless, its performance may be improved when the suggested fixes above are applied which we keep as the future work. In the next section, we give another estimator for effective resistances. In contrast to R 2sf , this one runs locally and gives an estimate for a single R i,j at each run.

Estimating ER via Local Forests

In the proof of Lemma 4.2.3, we compute the effective conductance I i,j = (R i,j ) -1 as follows:

I i,j = k∈N (i) w(i, k) 1 - (-1) i+k det L -i,-j|-k,-j det L -i,-j|-i,-j . ( 4.8)
A closer look on the determinant ratio gives us with Theorem 2.1.3 and Lemma A.4.1:

(-1) i+k det L -i,-j|-k,-j det L -i,-j|-i,-j = φ∈F k→i,j w(φ) φ∈F i,j w(φ) = q 2 φ∈F k→i,j w(φ) q 2 φ∈F i,j w(φ) = P(r Φq (k) = i) P(ρ(Φ q ) = {i, j}) = P(r Φq (k) = i|ρ(Φ q ) = {i, j}).
where F k→i,j is the set of forests with exactly two roots i and j and the node k = i, j is rooted in i. Then, plugging this probability in Eq. 4.8 gives:

I i,j = k∈N (i) w(i, k)(1 -P(r Φq (k) = i|ρ(Φ q ) = {i, j})) = k∈N (i) w(i, k)P(r Φq (k) = j|ρ(Φ q ) = {i, j}) 114 Chapter 4 RSF Estimation of Some Important Graph Quantities
Let us denote the conditional random forest Φ q |ρ(Φ q ) = {i, j} by Φ {i,j} . Then, an unbiased estimator for I i,j is:

Ĩi,j := k∈N (i) w(i, k)I(r Φ {i,j} (k) = j).
Based on this, we define a reciprocal estimator for R i,j over N samples as follows:

R lf i,j := N N k=1 Ĩ(k) i,j
.

The samples from Φ {i,j} can be generated by setting q i,j → ∞ and the rest to zero in Alg. 3. Moreover, one does not need to sample a complete spanning forest to compute Ĩi,j . Whenever the nodes in N are rooted, one can compute Ĩi,j without waiting all the nodes are rooted, hence the name. Therefore, in practice, we initiate the LERWs in Alg. 3 from N (i) and interrupt it whenever all the nodes in N (i) are rooted. We summarize this implementation for this estimator in Alg. 4.

Algorithm 4 EstimateERviaLocalForests

Inputs:

G = (V, E, w), {i, j} ∈ V,N Initialize: Ĩi,j ← 0 ∀k ∈ V \ {i, j}, q k ← 0 q i = q j → ∞ for i = 1 : N do
Run RandomForest on G by launching LERWs from N (i).

Interrupt the algorithm whenever all the nodes in N (i) are rooted, yielding a forest φ (not necessarily spanning).

Ĩi,j ← Ĩi,j + k∈N (i) w(i, k)I(r φ (k) = j) end for return N/ Ĩi,j .

Indeed, in this form, R lf is a biased estimator as:

∀i, j ∈ V, E[R lf i,j ] = E   N N k=1 Ĩ(k) i,j   = N E[ N k=1 Ĩ(k) i,j ]
.

The bias here can be approximated via Taylor expansion on E[R lf i,j ]:

E[R lf i,j ] -R i,j ≈ 1 E[ Ĩi,j ] 3
Var( Ĩi,j ) N

Estimators for Effective Resistances

By bounding Var( Ĩi,j

) = E[ Ĩ2 i,j ] -E[ Ĩi,j ] 2 ≤ d 2 i -I 2 i,j
, one has a bounded bias:

E[R lf i,j ] -R i,j 1 N I 3 i,j (d 2 i -I 2 i,j ).
Moreover, the variance can be approximately bounded as follows: We complete our theoretical analysis on R lf i,j by noting its variance. The Taylor expansion on the variance gives us:

Var(R lf i,j ) ≈ 1 N Var( Ĩi,j ) I 4 i,j ≤ 1 N d 2 i -I 2 i,j I 4 i,j
.

Notice that the bias diminishes at the rate of O(N -1 ) compared to O(N -1/2 ) for the standard deviation. Thus, this estimator converges to the true expectation R i,j for large N . Nevertheless, one can decrease or completely eliminate the bias by making some extra effort. A simple correction is to estimate the bias term from the samples via the sample mean and variance formulas and subtract it. If the time budget is sufficient enough, one can deploy more advanced methods [START_REF] Babu Moka | Unbiased estimation of the reciprocal mean for non-negative random variables[END_REF] that can provide an unbiased estimation. Yet, in certain cases, they might increase the variance of the estimator. In the experiments of the next section, we restrict ourselves to simple bias correction for the sake of simplicity. Yet we finish this section by stressing that there is still room for improvement.

Empirical Results on ER Estimation

In this section, we empirically compare the forest-based estimator R lf with the state-of-the-art algorithms in terms of approximation error and run-time. However, a fair comparison is not straight-forward as some algorithms are more favorable for estimating a small number of effective resistances while the others are dedicated to estimate effective resistances in a global sense (over all edges or all pairs). Therefore, we first divide the compared algorithms into two groups called local and global. The local algorithms run on a small portion of the graph and estimate efficiently R i,j for a given pair (i, j). In this group, we have algorithms 1 and 4 in [START_REF] Peng | Local Algorithms for Estimating Effective Resistance[END_REF], namely

EstEff-TranProb (TP) and EstEff-MC2 (MC2), and the estimation via local forests (LF) of Section 4.2.4. The global algorithms approximate the effective resistances over all edges or all pairs by operating over the whole graph. The run-time of these algorithms does not depend on the number of effective resistances demanded. Therefore, the number of estimated effective resistances grows, these algorithms become more favorable compared to the local algorithms. In this group, we put the spanning tree algorithm by [START_REF] Hayashi | Efficient Algorithms for Spanning Tree Centrality[END_REF] (ST) and the random projection algorithm by [START_REF] Daniel | Graph sparsification by effective resistances[END_REF] (RP).

Notice that all algorithms in this comparison are Monte Carlo algorithms and each has a parameter that adjusts the number of Monte Carlo trials. For global algorithms ST and RP, these parameters are respectively the number of sampled spanning trees N st and the dimension of the projection space k. We start our experimentation by estimating effective resistances via these algorithms over all edges. In doing so, we vary the parameters N st and k and measure the total run-time and the average relative error 1 m (i,j)∈E |R i,j -R i,j |/R i,j of the algorithms. This gives us the performance of the global algorithms in estimating a single effective resistance. For the local algorithms TP, MC2 and LF, we run them over 30 randomly selected edges (with replacement) in order to estimate the corresponding effective resistances. Finally, we compute their average relative error and the average run-time over the edges in order to deliver their performance over a single effective resistance. Their parameters are tuned for each edge as follows; we first run LF for 40 local forests, then we tune the parameter l in TP and in MC2 in order to reach the relative error given by LF. This allows us to compare the run-times of local algorithms for the approximately same relative error. In order to compare them with the global algorithms, we calculate a number called threshold number N th . This is the maximum number of effective resistances that can be estimated by a local algorithm before reaching the run-time of the fastest global algorithm at the same relative error. By looking at this metric, one can answer the following question; For how many effective resistances, the local algorithms are more useful than the global ones? As a result, the larger this metric for a local algorithm the better the performance.

In Fig. 4.4, we present our results over 4 benchmark datasets which are listed as: First of all, we observe that over all graphs, ST provides the best performance within global algorithms in estimating effective resistances over all edges. Note that both methods inherently benefit from the sparsity of the graph. In case of RP, the sparse approximate Cholesky decomposition exploits the sparsity to facilitate its intermediate operations. For ST, Wilson's algorithm generates spanning trees by taking much less time in sparse graphs. The difference we observe in these 4.2 Estimators for Effective Resistances 10 3 10 2 10 1 10 0 10 1 10 2 10 3 Runtime (s) simulations is probably due to the low variance of ST for estimating large effective resistances. One can observe this on the extreme case. In this case, we estimate an effective resistance R i,j where i and j is connected only by a single edge. By definition, we find the true value of R i,j as 1. Also, the estimation by ST always equals to 1 as this edge has to appear in all spanning trees. Therefore, it has zero variance. On the other hand, this is not necessarily the case for RP as the randomization we use here is independent from the graph. Therefore, the estimation given by RP still might have variance around the true value of R i,j . We see that there are many edges that fit this description on these graphs and thus we observe a significant difference between ST and RP. However, one needs to keep in mind that ST is not defined for the pairs (i, j) ∈ E whereas RP can estimate those effective resistances.

Secondly, all local algorithms are usually more efficient compared to the global ones for estimating a single effective resistances. However, in the global sense, they cannot provide the same efficiency. Thus, we report their threshold numbers N th in order to say up to how many edges they are still useful. By looking at the estimation over a single edge and the threshold numbers, we see that the forest-based algorithm LF often gives the best performance. Compared to other local algorithms TP and MC2, it takes much less time for the same relative error over a single estimation, and it gives a threshold number that is at least 4-5 times that of the closest local algorithm. As a result, LF gives a way to estimate approximately 0.5 -2% of the effective resistance over all edges without being slower than the global algorithms.

Graph Filtering with RSFs

Graph signal filtering is an essential tool in GSP. However, the necessary computations, such as computing the graph Fourier basis U by diagonalizing the graph Laplacian L = UΛU , are not scalable with the size of the graph. The RSF estimators x, x and z avoid the heavy computations for approximating a particular graph filtering where the transfer function is:

g q (λ) = q q + λ (4.9)
where λ is the graph frequency parameter. Indeed filtering with g q (λ) is a very useful operation as it returns a smooth signal on the graph. Yet depending on the application, different transfer functions might be desired. For example, the ideal low-pass filters are more favorable if the noise is known to be contained in the high 4.3 Graph Filtering with RSFs frequency components. In these filters, the transfer function takes the following form:

g k (λ) =    1, if λ ≤ λ k 0, otherwise (4.10)
where λ k is the k-th smallest eigenvalue of L.

When applied, this filter completely removes the components associated with the frequencies λ > λ k from the signal. If the noise is known to be contained in those frequency components, one can perfectly reconstruct the original signal. The step function behaviour of g k (λ) is desired in many applications. In our case, we may try to approximate g k (λ) by g q (λ) with the forest estimators by optimizing q such that g q (λ) gives the greatest decrease at λ k . Even so, g q (λ) is a very poor approximation.

In this section, we focus on the graph filtering other than g q (λ) that can be approximated by the forest-based estimators. In doing so, we first give a forest-based method to solve the Laplacian systems in the form of Lx = y which corresponds to low-pass filtering. Then, we switch our focus to another setup where we find a family of graph filters estimated by RSFs. In this setup, we duplicate the graph and add edges between the nodes of the original and copy graph. For certain ways of creating such graphs, the RSF methods yield a set of graph filters different from g q (λ). Instead of a single hyper-parameter as in g q (λ), these filters come with four parameters, yielding wider options for filtering. We illustrate them in a case study where we approximate some of the typical graph filters g k (λ) or, g q,p (λ) = q q + λ p , (4.11)

where p > 1 is an integer.

Solving Lx = y with RSFs

Linear systems of the form Lx = y where L is a graph Laplacian are an important class of systems of linear equations. Apart from their direct use in electrical networks, they also appear in the spring networks and maximum flow problems [START_REF] Nisheeth | Lx= b[END_REF]. It is possible to see the solution x = L † y as a graph filter where the transfer function is:

g † (λ) =    0 λ = 0 1 λ otherwise . 120
Chapter 4 RSF Estimation of Some Important Graph Quantities Yet again, directly computing x (O(n 3 )) is not practical for very large graphs. In order to avoid this, one can use classical methods for solving linear systems, such as iterative methods. The state-of-the-art algorithms solve the Laplacian systems in quasi-linear time in the number of edges, more specifically O(m log c n log 1 ) where c is a constant and is the error parameter [START_REF] Daniel | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF]. In this section, we show how RSFs can be leveraged for approximating L † y. However, the resulting method at its current state is not competitive with SOTA algorithms.

The forest estimator we find for L † y in fact, takes its root from the following relation between L † and the effective resistance matrix R:

L † = - 1 2 R - 1 n (RJ + JR) + 1 n 2 (JRJ) , (4.12)
where J = 11 ∈ R n×n . Then, whenever one has an estimate R of R, the product L † y can be approximated by matrix-vector products with R. Interestingly, the forest estimator R 2rf may provide such products very cheaply. By definition R 2rf is a symmetric matrix with two non-zero blocks in its off-diagonals. Then, the product x 2rf := (R 2rf y) with an arbitrary vector y reads:

∀i ∈ V, x 2rf i =      j∈V y j I(r Φq (i) =r Φq (j)) C|V t(Φq ,r 1 ) ||V t(Φq ,r 2 ) | if ρ(Φ q ) = {r 1 , r 2 } 0 otherwise
where r 1 and r 2 are the random roots of Φ q . In other words, given a forest φ with two roots, it computes a scaled average of y within the tree that i does not belong to, and propagates it through the tree of i (See Fig. 4.5). Assuming 1 y = 0 without loss of generality9 and rearranging on Eq. (4.12) allows us to define the estimator for L † y as follows:

l := - 1 2 I - 1 n J x 2rf .
In practice, this estimator computes x 2rf , removes its mean and scales the result by -1/2. This ensures that the estimates given by l are zero-mean vectors similar to L † y. As defined above, the random vector x 2rf takes a two-tree spanning forest and computes the averages within each tree. Unlike the previous estimator x, it swaps the scaled averages between trees and then propagates. At first sight, the swapping operation might seem counter-intuitive as we do not have such behaviour in the estimated product L † y. In fact, this situation can be explained by examining the overall process. First, notice that for 1 y = 0, the scaled averages computed in two trees ȳ1 and ȳ2 verify ȳ1 = -ȳ 2 . Then, swapping the scaled averages across the trees corresponds to changing only the signs of the values that are propagated. Later on, in calculating l, we multiply x 2rf by -1 2 , which reverts the sign change that is induced by the swapping.

The intrinsic weakness of this estimator is due to the high variance of R 2rf . Due to this fact, we do not expect a good performance in estimating L † y. Nevertheless, one may improve this estimator by improving the performance of R 2rf by taking our suggestions in Section 4.2.3.

Filters on the Duplicated Graph

Let us start by describing how we duplicate the graph. Given the adjacency matrix W of a graph, the duplicated graph G d admits the following adjacency matrix:

W d = W W b W b W ,
where W b ∈ R n×n contains non-negative edge weights. In other words, in order to

A B C D A' B' C' D'
We copy the graph and add edges between the nodes of the original and copy graph form the duplicated graph G d , we duplicate G and create edges between the nodes of the original and copy graph where the weights of these edges are contained in W b . Indeed the options for W b are limitless. Here we restrict ourselves to (entry-wise) non-negative matrices that are matrix functions evaluated at L i.e. W b = h(L) := Uh(Λ)U for some function h : R → R. In turn, we ensure that every block of L d is diagonalizable by U. This constraint seems a bit technical, but it allows us to keep U as the Fourier basis in the analysis of the RSF filters proposed in the following.

Let us delve into the filtering properties on G d . In order to do so, we first write the corresponding graph Laplacian L d . By observing h(L)1 = h(0), we can simplify L d to the following block matrix form:

L d = D d -W d = h(0)I + D 0 0 h(0)I + D - W h(L) h(L) W = h(0)I + L -h(L) -h(L) h(0)I + L .

122

Chapter 4 RSF Estimation of Some Important Graph Quantities

Given this graph, let us take a look at the filtering induced by the Tikhonov regularization, which computes the smooth signal as:

x = q(L + qI) -1 y.

First of all, G d contains 2n nodes thus 2n measurements y 1 , . . . , y 2n . Therefore, we need to clarify the definition of the measurements vector. A natural choice is to duplicate the original measurements. In order to explore more graph filters on G d , we define a parametrized duplication as y d = αy βy . Similarly, instead of using the same value of q in both the original and copy graph, we give the flexibility of choosing different values q 1 and q 2 so that we can explore the full potential of the new set of filters. Under this parametrization, the corresponding kernel K d becomes:

K d = (L d + Q d ) -1 Q d = (h(0) + q 1 )I + L -h(L) -h(L) (h(0) + q 2 )I + L -1 q 1 I 0 0 q 2 I = q 1 M 2 S q 2 h(L)S q 1 h(L)S q 2 M 1 S = K 1 K 2 K 3 K 4 (4.13)
where M 1,2 = (h(0) + q 1,2 )I + L and S = (M 1 M 2 -h 2 (L)) -1 . Then, the product with y d gives:

K d y d = (αK 1 + βK 2 )y (αK 3 + βK 4 )y .
Notice that all K 1,2,3,4 commutes with L thanks to our initial constraint. Then, we can derive the transfer function of the filtering implemented by (αK 1 + βK 2 )y as follows:

f θ (λ) = αq 1 (λ + h(0) + q 2 ) + βq 2 (h(λ)) (λ + h(0) + q 1 )(λ + h(0) + q 2 ) -h(λ) 2 , (4.14)
where we collect the hyper-parameters into θ = (q 1 , q 2 , α, β). This transfer function with four hyper-parameters gives us a big family of filters, and all of these filters can be approximated by the unbiased estimators x, x and z. On top of the hyperparameters, a significant design choice is the selection of the matrix function h. Unlike θ, we are not completely free on the choice of h as it needs to verify that h(L) is a non-negative matrix. Some simple examples are:

The functions that map the graph Laplacian to non-negative matrices have been analyzed in several works by [BCN05; BH08; MW79]. In fact, these works consider a more general class of matrices, called M -matrices, instead of the Laplacian. The M -matrices are defined as the real matrices which contain only non-positive values in the off-diagonal entries. As a result, all graph Laplacians are also M -matrices.

Then the following theorem gives us a necessary condition for the function h to build h(L) = W b .

Theorem 4.3.1 (Necessary condition for h [START_REF] Rb Bapat | On functions that preserve M-matrices and inverse M-matrices[END_REF]). Let h : R + ∪ {0} → R be analytic function. The function h maps M -matrices to non-negative matrices only if h verifies for all i = 1, 2, . . . and for x ∈ (0, ∞):

(-1) i h (i) (x) ≥ 0,
in other words, h is completely monotonic.

Proof. We refer the reader to the proof of Theorem 2.1 in [START_REF] Rb Bapat | On functions that preserve M-matrices and inverse M-matrices[END_REF].

Thanks to this necessary condition, we can check the feasibility of a candidate function for building G d .

Overall, we define here a filtering scheme on an altered version of the graph. The way we alter the graph allows us to keep the same graph Fourier basis with the filtering schemes on the original graph. Further analysis give us a wide set of transfer functions that can be obtained from this new filtering scheme. Most importantly, all of these filters can be approximated by the RSF-based estimators. This fact opens many possible use cases for the forest estimators beyond approximating g q (λ).

The following case study is dedicated to illustrate these filters. In this illustration, we empirically show that there exist feasible parameter set θ that allows us to approximate the ideal low-pass filter g k (λ) or the (p-q)-rational filter g p,q (λ) = q q+λ p better than g q (λ) .

Low-pass Forest Filters

Filtering a signal by low-pass filters without diagonalizing L is usually desired in GSP [START_REF] Ramakrishna | A user guide to low-pass graph signal processing and its applications: Tools and applications[END_REF]. In approximating such filters, one can deploy approximate approaches such as Chebyshev polynomials [START_REF] David I Shuman | Chebyshev polynomial approximation for distributed signal processing[END_REF] or graph ARMA filters [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF] yielding state-of-the-art performance. Here we look at a dif-ferent perspective on approximating low-pass filters. Since the signals filtered by g q and f θ can be cheaply estimated by RSF-based estimators without any bias, in this example, we aim to explore the potential of approximating two particular low-pass filters g k (λ) and g p,q (λ).

First, let us fix some notation. In a more general sense, we consider a target transfer function g (λ) and try to approximate g (λ) with a candidate function c Θ (λ) with a parameter set Θ. As a result, we look for filter parameters Θ such that:

Θ = argmin Θ n i=1 (g (λ i ) -c Θ (λ i )) 2 + µ||c Θ (0) -1|| 2 2 ,
over the graph frequencies λ i 's a and µ||c Θ (0) -1|| 2 2 is the regularization term that enforces that the estimated filter equals to 1 at λ = 0 similar to the target filters. Under this generic setup, we consider g k (λ) or g p,q (λ) as the target function, and g q (λ) or f θ (λ) as the candidate function. In both options of c Θ (λ), here we have a constrained optimization problem as we need to ensure q > 0 for g q (λ) or q 1,2 > 0 for f θ (λ). On the other hand, this type of approximation for the target functions is quite appealing as we only need to solve an optimization problem with a few parameters. In the case of g q (λ) being the candidate function, one can find a global minimizer solving the corresponding single variable polynomial under the constraint q > 0. Indeed a similar analysis should be done for f θ (λ), which comes in a much more complicated form. In this case study, we prefer to avoid this type of analysis and empirically show that f θ (λ) may provide a better approximation for the target functions even if the parameter set θ is a local minimizer. In the following, we illustrate this idea on a 3D-Grid graph with 1000 nodes and 3000 edges. We approximate two target filters g k (λ) and g p,q (λ) via g q and f Θ . For f Θ , we choose h(L) = d max I -L. Then we optimize the parameters of the candidate filters g q and f Θ for estimating the filters g k=200 (λ) and g p=2,q=10 (λ). In the first figure, we find the parameter q = 2.14 for g q and Θ = (0.61, 1.70, 8.02, 18.16) for f Θ . As you can see f Θ yields much more approximate filter when these parameters plugged compared to that of g q . The similar results are also observed in the case of approximating g k=200 (λ).

The filter with four parameter f Θ where Θ = (0.25, 2.93, 10.80, 13.88) have more flexibility compared to g q=1.63 , thus it yields a better approximation.

a In practice, this can be replaced by the estimation of the spectrum as it is expensive to calculate all the eigenvalues

Conclusion

In this section, we present the forest-based algorithms to approximate three important graph quantities, namely the trace of the regularized inverse of the graph Laplacian, the effective resistances and various graph filters. Firstly, the proposed methods for the trace estimation problem in Section 4.1.1 are the improved versions of the forest estimator introduced in [START_REF] Barthelme | Estimating the inverse trace using random forests on graphs[END_REF]. Although these improvements take their root from the existing variance reduction methods [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF] (control variate method and stratified sampling), it is not evident how to apply them on an arbitrary estimator as they require additional statistics. Our contribution here is to provide such statistics so that these variance reduction techniques can be adapted to improve estimation performance. In turn, we observe that the proposed methods reach the SOTA performance over various graph datasets. In estimating effective resistances, we propose two novel algorithms, which are local and global estimators.

Our analysis shows that the global algorithm we propose has interesting properties but has poor estimation performance. On the other hand, the local algorithm yields simple implementation and comparable performance with the SOTA algorithms as observed in Section 4.1.2. Finally, we show that the graph filters that can be approximated via RSFs are not limited to q q+λ . In doing so, we first show that L † y is another filtering operation that can be estimated via RSFs, in particular 2-rooted random forests. Then, we slightly modify the graph by extending the filters that can be obtained RSFs. In this modification, we duplicate the graph by keeping the original edges and connecting the nodes in the original and copy graphs. In turn, this gives a family of filters parametrized by four parameters. Finally, we illustrate this filter in approximating certain graph filters such as g k and g q,p .

" The song is ended but the melody lingers on...

-Irving Berlin

The graph Laplacian is at the heart of many applications involving graph-structured data. Given their extensive use, many numerical algebra tools (deterministic or randomized) have been adapted for dealing with the graph Laplacian. In this thesis, we give a brand new set of techniques for randomized Laplacian-based algebra. These techniques cover a wide variety of problems, including graph Tikhonov regularization, interpolation, node classification, graph 1 -regularization, trace estimation of the regularized inverse of L, estimation of the effective resistances and graph signal filtering. Although these problems seem entirely unrelated, their solutions are closely linked with random spanning forests. Some already known and others uncovered through the thesis, these links have been the supporting pillars of the techniques proposed through this thesis. They take their root in the connections between various concepts, namely graph theory, determinantal point processes and randomized linear algebra. Therefore, we started by giving the necessary background on these concepts and also random spanning forests by emphasizing their connections with RSFs. Then, we presented our randomized methods for the problems in Laplacian-based numerical algebra. For each method, we provided a bias-variance analysis, and we compared them with the existing algorithms over various real datasets whenever possible.

Limitations, Open Questions and Perspectives

In this section, we will discuss each proposed method of this thesis by pointing out its strengths and weaknesses while highlighting associated open questions and future directions.

In Chapter 3, we present a unifying framework in Prob. (3.1.1) for Tikhonov regularization and interpolation on graph signals. While the former draws a regularized unconstrained optimization problem, the latter dictates the value of the signal at certain nodes. With the help of the node-wise hyperparametrization (defining q i per node i), we show that both cases can be understood as an instance of Prob. (3.1.1). Besides the wide variety of applications, from graph signal filtering to electrical resistor networks, the explicit solution to the unified problem can be found in a closed-form expression. One of our contributions in this chapter is the forest-based estimator x for approximating this solution x for all cases. In a nutshell, this esti-Given a graph and signal... the unbiased estimator x operates as follows:

mator samples rooted spanning forests and propagates the signal value at the root within each tree. As simple as it is, the bias-variance analysis shows it has no bias and a tractable variance in terms of the input signal and the Laplacian. Moreover, it is possible to apply variance reduction techniques to reduce its expected error significantly. This brings us to the other contributions of Chapter 3. Thanks to the rich theory of RSFs and existing algorithms, we manage to adapt two variance reduction techniques, namely conditional expectation and the control variate method. We obtain new estimators with no bias and reduced variance in both cases.

On the other hand, the existing algorithms, namely preconditioned conjugate gradient descent and Chebyshev polynomial approximation, to compute x are deterministic algorithms with a high convergence rate. Specific to inverting Laplacian systems, one can find more sophisticated preconditioning methods for iterative methods based on spanning trees or graph sparsification [START_REF] Nisheeth | Lx= b[END_REF]. In turn, they speed up the convergence of the algorithm, yielding nearly linear time algorithm with m. In contrast, the proposed methods are subjected to the Monte Carlo convergence rate O(N -1/2 ) i.e. in order to divide the squared error ||x -x|| 2 2 by 2, one needs to double the sample size. Over various graph datasets in Section 3.3.3, we observe that this rate remains slow compared to the convergence of the deterministic methods. Yet again, in many applications, the exact solution x itself is an inaccurate quantity. For example, in a signal denoising setup, the noise-free signal differs from x by some error. Keeping this in mind, we also measure the denoising performance of our estimators and see that they are comparable with state-of-the-art algorithms.

Indeed, the slower convergence rate is the principal issue of forest-based estimators. On top of this, certain compiler and hardware optimizations give the deterministic methods an advantage, while these are not yet available for the RSF estimators. The deterministic methods access the graph Laplacian L only by the matrix-vector products. These operations are well-optimized in the compilation by considering run-time performance and memory accesses. In the case of forest-based estimators, one needs to run loop-erased-random-walks over the graph, which access the entries of L in arbitrary indices. This brings an additional cost to pull the random parts of L into the computer's cache memory. Given the architecture of today's computers, this cost is inevitable. However, one may still look for speeding up in Wilson's algorithm itself. In fact, this has been an open research question during the last decades, and many have attempted to propose a faster algorithm than Wilson's for sampling uniform spanning trees. Many attempts in this field yielded approximate algorithms i.e. samples a spanning tree from approximately uniform distribution [KM09; MST14; Sch18; ALGVV21]. While some use the links of USTs with the graph Laplacian, some others choose MCMC methods to approximate the distribution of UST with certain theoretical guarantees. Although these methods are theoretically appealing, they require implementing either an algebraic tool or a specific data structure which makes the implementation complicated. On the other hand, Wilson's algorithm is very easy to implement as it is presented in Alg. 2. The most time-consuming step in Wilson's algorithm is usually the first LERW that only terminates at a single node which is the given root (Γ for the forests). As no other path is already sampled, the first LERW might take a lot of time to visit a lot of nodes before its termination. On the other hand, the last LERWs terminate very quickly as they may stop at many nodes. Recently, [START_REF] Nunes | A transient equivalence between Aldous-Broder and Wilson's algorithms and a two-stage framework for generating uniform spanning trees[END_REF] states that the algorithm by Aldous [START_REF] David | The random walk construction of uniform spanning trees and uniform labelled trees[END_REF] and Broder [START_REF] Broder | Generating random spanning trees[END_REF], which was the fastest algorithm before Wilson's algorithm, has the exact opposite behaviour. Based on this fact, the authors in [START_REF] Nunes | A transient equivalence between Aldous-Broder and Wilson's algorithms and a two-stage framework for generating uniform spanning trees[END_REF] suggest a hybrid framework that deploys Aldous-Broder's algorithm at the initial steps and launches Wilson's algorithm later on to complete the spanning tree. They show that for certain graphs, this procedure still samples a UST faster than both algorithms. However, at its current stage, it is not generalized to general graphs. Overall, Wilson's algorithm remains the current best option for our estimators.

In Chapter 4, we focus on the RSF estimation of other graph quantities, which are the trace of the regularized inverse of Laplacian, effective resistances and various graph filters. Each quantity appears as the solution to a problem involving the graph Laplacian. Firstly, the trace tr(q(L + qI) -1 ) is often required in the hyper-parameter selection algorithms in order to tune the hyper-parameter q in GTR. The existing algorithms usually estimate it via Monte Carlo simulations, such as Hutchinson's estimator [START_REF] Michael F Hutchinson | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF] and the number of roots in RSFs [START_REF] Barthelme | Estimating the inverse trace using random forests on graphs[END_REF]. We work on the RSF-based algorithm, the number of roots s, to improve it with variance reduction methods. In doing so, we adapt two variance reduction algorithms, namely control variate and stratified sampling, for reducing the variance of s yielding the best, or at least comparable, results with the state-of-the-art algorithms. In adapting control variate method, we follow a very similar procedure while we transform the estimator x (resp. x) into z (resp. z). In turn, one only needs to track the edges connecting two nodes belong to different trees for computing the control variate. To implement this, we only add a few lines in Wilson's algorithm, and we observe in empirical tests that this is a cheap price to pay given the variance reduction provided. However, our implementation gives a naive way of tracking edges at the boundaries. An open question is whether there exist faster algorithms or, more probably, data structures that allow us to track these edges without traversing all the edges. A possible direction might be to implement link-cut trees as done in [START_REF] Luıés | Linking and cutting spanning trees[END_REF] in order to track the edges of interest meanwhile generating the spanning forest. The second variance reduction method we adapt is stratified sampling. We divide the outcome space of s = |ρ(Φ q )| by using a random quantity |ρ 1 (Φ q )| (the number of roots sampled at the first visit of LERWs in Wilson's algorithm) while collecting samples. This quantity has a tractable probability distribution (Poisson-Binomial) and it is easy to sample from. Therefore, we choose to use the stratification variable to decrease the variance of s.

130 Chapter 5 Conclusion S = S (1) S (2) S (3) S (4) S (5) S (6) S (7)           2 7 6 1 Γ Γ 1 Γ 7 2 5 4 4 2 3 7 1 5 1 4 1 . . .
In turn, the new estimator competes and even outperforms the existing algorithms. Indeed, one may achieve even better performances by combining these two methods by going beyond this work.

Another interpretation of the random set ρ 1 (Φ q ) comes from the stack representation of Wilson's algorithm (See also Fig. 5.1). We consider infinite stacks at each node. At every layer k of a stack S (i) k associated with node i, one has a random neighbor j with a probability w(i,j) q+d i and with probability of q q+d i , the node i becomes a root. Then, we form a graph by looking at the top element of each layer. Typically, such a graph contain cycles in it. A cycle-popping algorithm removes the edges of the cycles from the corresponding stacks by updating their top elements. Surprisingly, Wilson's algorithm (for forests) is an implementation of a cycle popping algorithm that runs until there is no cycle remaining. The resulting structure is a spanning forest. When we look at the root set of the spanning forest, we see that it evolves through this process until termination. While some roots are sampled in the first layer of the stacks, some others become a root in the continuation of the cycle popping procedure. In fact, ρ 1 (Φ q ) is the roots sampled at the first layer, and this definition can be generalized to ρ k (Φ q ) which are the roots sampled until k-th layer. Interestingly, we see that for k = 1, ρ 1 (Φ q ) is a DPP with a diagonal marginal kernel K 1 = q(qI + D) -1 . Similarly as k goes to infinity, one recovers the original root set 131 ρ(Φ q ) with the marginal kernel K k→∞ = K = q(qI + L) -1 . Then, the natural question is that what happens in between? We wonder whether ρ k (Φ q ) is still a DPP for some k > 1, or the only cases where the DPPs appear is k = 1 and k → ∞. If the former is true, what are their marginal kernels? These questions remain unanswered in this thesis, but one might start to investigate by attempting to write the distributions of these objects by using the distributions of the underlying LERWs.

We also propose two RSF-based estimators for approximating effective resistances without explicitly computing L † . Both rely on the close connections between two rooted spanning forests and effective resistances. In one, we compute the effective resistances in a global sense (over all pairs of nodes), while the other runs in a much smaller portion of the graph to estimate a single effective resistance. The bias-variance analysis shows that both estimators are biased. Nevertheless, the bias decays much faster than the variance as the sample size increases, leaving the variance as the main source of the approximation error. In initial tests, we observe the global estimator we propose has poor performance. However, it can provide an estimate over all pairs of effective resistances for a single sample of two rooted spanning forests. In its variance analysis, we see that the variance is inversely proportional with c 2 where c = ( i,j∈V R i,j ) -1 . For large graphs, c 2 is typically a small number yielding a large variance in our estimator. Yet again, one might find improvements via recycling unused forests that are not with two roots or control variate suggested in Section 4.2.3. Although it gives poor estimation, it gives cheap matrix-vector products. In turn, one can estimate L † y for an arbitrary vector y at very cheap cost as described in Section 4.3.1.

The local estimator on the other hand gives a good estimation performance and is comparable with state-of-the-art algorithms. It usually outperforms other existing local algorithms in the maximum number of effective resistances that can be estimated before reaching the best global algorithm. In this algorithm, we set node i and j as roots of estimating R i,j , and we launch LERWs at the neighbours of i. We stop the algorithm when all neighbours of i are rooted. In this way, the algorithm produces an estimate without generating a complete spanning forest of the graph. As aforementioned, the first LERW again take the most of the run-time. An interesting modification might be an early stop strategy in which we stop the first LERW after reaching a certain number of steps, which introduces an additional bias in our estimation. Indeed, an open question is if this version of the algorithm yields a substantial run-time gain without significantly increasing the expected error of the estimator. One might start to analyze the bias by deriving the distribution of the interrupted LERWs, which might not necessarily appear in determinantal ratios as in the distributions of the original laws.

Later, we give a set of graph filters that RSFs can approximate. In this set, we have three types of graph filters whose transfer functions are g q (λ) (the graph filter in GTR), g † (λ) (the filter by the pseudo-inverse of Laplacian) and f Θ (λ) (the filter obtained by duplicating the graph). As mentioned before, g † (λ) or solving Laplacian linear systems is at the reach of the RSF-based estimators thanks to our global forestbased estimator for effective resistances. The implementation is simple and efficient, but the high variance of the global estimator is also the biggest issue for estimating g † (λ). Therefore, we suggest using this estimator after applying the improvements in Section 4.2.3. This leaves us with the parametric filters g q (λ) and f θ (λ). By duplicating the graph and solving a version of graph Tikhonov regularization on the duplicated graph, we obtain a new set filter f θ (λ). Indeed, this filter is also within reach of the RSF-based estimators. Moreover, it is parametrized by four parameters instead of a single one q, which provides a much wider range of filters. We illustrate this extended range in approximating popular graph filters such as the ideal low-pass filters or g q,p>1 (λ) = q q+λ p and observe that the duplicated graph filter f θ (λ) may provide better approximations than g q (λ) with the help of more parameters. Given this increased capability, one might ask; What happens if we copy the graph more than twice? In this case, one would obtain more degree of freedom in the filters while the parameters to select increase by 2 per added layer. One can increase the versatility of

A B C D A' B' C' D' A" B" C" D"
the filters by choosing the wiring function h(λ) differently between different layers. All of these lead us to a complicated function analysis to discover the full potential of the filters in duplicated graphs, which we leave as future work.

As a final remark, we would like to stress that all the methods proposed on Laplacianbased numerical algebra can be adapted for symmetric diagonally dominant (SDD) matrices, which is a larger class of matrices than the graph Laplacian. The trick by [START_REF] Keith D Gremban | Combinatorial preconditioners for sparse, symmetric, diagonally dominant linear systems[END_REF] allows us to transform the linear systems with SDD matrices into larger linear systems with graph Laplacian. Using this connection, one can replace L with an SDD matrix in GTR, and yet the solution can still be approximated via the RSF estimators. A similar adaptation is already given in [START_REF] Barthelme | Estimating the inverse trace using random forests on graphs[END_REF] for the trace estimation problem. Given this extension, one might hope to extend larger classes of matrices than SDD matrices. In this vein, a promising work [Ökt05] gives a random walk-based solution to solve linear systems in the form of Ax = b where A is only constrained to be an entry-wise real matrix. Here, an interesting perspective would be replacing LERWs, and so RSFs, in the place of simple random walks for solving Ax = b. So far, this adaptation is not obvious to our eyes, but we find investigating it promising direction to extend the scope of our methods further.

Epilogue

At first sight, this thesis's main achievements can be considered the RSF-based algorithms for several graph-based problems. Especially, the forest-based algorithms for trace estimation and effective resistances show ensuring performance when compared to existing algorithms. Similarly, the estimators for GTR and interpolation are shown to be useful when one does not need to compute the solution at very high accuracy. These methods are very easy to implement and flexible for parallelization or distributed implementations. In the bigger picture, in this thesis, we reconsider RSFs as tools to sample a small portion of the graph and solve the problem at hand on this small portion to avoid heavy computations. Compared to the existing set of tools for SDD matrices [Vis+13; ST04] and DPP-based approaches [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF], the RSF-based tools bring a brand new and elegant perspective to the literature.

Up to our work, a similar perspective is usually presented with simple random walks by emphasizing the connections between the Markov chains and the graph-related algebra [START_REF] Ben-Hamou | Estimating graph parameters with random walks[END_REF]. In this context, for example, the works by [Wu16; WLSWC12] provide a random-walk insight to several practical machine learning algorithms solving GTR, interpolation or clustering. Here, taking a step forward from random walks to loop-erased random walks carries us to the prosperous realm of DPPs and RSFs. In turn, we give practical algorithms to solve very central problems. In addition, we observe that the direction we take has inspired recent research on randomized graph-related algebra. For example, in their recent work [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF], the authors look at the problems of graph sparsification and solving Laplacian-based linear systems with multi-type spanning forests, which are certain types of DPPs and closely related to RSFs. In the future, we believe that RSFs and other types of DPPs will have more impact on applications of SDD-based linear algebra and hopefully extend to the other classes of matrices.

Resumé Substantiel

  1 -1 0 -1 2 -1 0 -1 1  
représentations matricielles des graphes, en particulier du Laplacien de graphe. Le laplacien est une matrice carrée symétrique de la taille du nombre de sommets. Il

Bibliography peut être considéré comme l'opérateur laplacien discret sur les graphes lorsqu'il est appliqué à une fonction sur les sommets :

(Lf ) i = j∈N (i) (f (i) -f (j)),
où N (i) désigne les noeuds qui partagent un bord avec le noeud i (la définition formelle se trouve section 2.1.1). Le laplacien est une représentation matricielle simple mais utile des graphes. L'analyse algébrique le relie à de nombreuses propriétés utiles des graphes. Par exemple, de nombreuses tâches dans l'analyse de la connectivité du graphe peuvent être accomplies en observant la décomposition de la matrice du Laplacien du graphe. de l'arbre matriciel de Kirchoff stipule que le produit des valeurs propres non nulles de L est égal au nombre de tous les arbres couvrants, i.e. sous-graphes connectés qui ne contiennent pas de cycles, sur un graphe.

En dehors de ces exemples, le laplacien joue un rôle important dans l'analyse des réseaux de résistances électriques. Dans ces réseaux, nous modélisons chaque résistance par une arête où le poids de l'arête représente la conductance de la résistance, et chaque noeud correspond à ses points de connexion. Ces modèles sont largement utilisés dans la théorie des circuits afin de calculer les courants circulant à travers les résistances/arêtes ou les tensions induites aux noeuds lorsqu'une différence de potentiel fixe (tension) est appliquée entre deux noeuds (ou plus). Étonnamment, ces calculs se résument à la résolution de systèmes linéaires impliquant le laplacien du graphe (voir le chapitre 4).

Au cours des dernières décennies, des exemples encore plus nombreux sont apparus dans le traitement du signal et l'apprentissage automatique impliquant des graphes.

Laplacien dans le traitement du signal des graphes. Le GSP est le sous-domaine du traitement du signal qui traite des signaux définis sur les sommets du graphe. Au cours de la dernière décennie, de nombreux outils classiques de traitement du signal ont été adaptés afin de traiter de tels signaux. Citons par exemple le filtrage [SVF11; SNFOV13], l'échantillonnage [TAB17; PTGV18], l'opération de translation [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF], l'analyse en ondelettes [HVG11; NO12; ACGM20] ou le principe d'incertitude [START_REF] Tsitsvero | Signals on graphs: Uncertainty principle and sampling[END_REF]. Dans beaucoup de ces adaptations, la transformation de Fourier de graphe joue un rôle important. Par analogie, elle nous permet de représenter les signaux du graphe dans le domaine de fréquence du graphe. À son tour, on peut définir des schémas de traduction et de filtrage des graphes en utilisant un analogue du célèbre théorème de convolution [START_REF] Oppenheim | Discrete-time signal processing[END_REF]. Cependant, en raison de l'irrégularité du domaine du signal, la définition de la transformée de Fourier n'est pas unique. Une définition populaire de la PSG repose sur le fait que les fonctions propres de l'opérateur Laplacien sont les fonctions de base de la transformée de Fourier. En s'appuyant sur ce principe, les auteurs de l'article [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] suggèrent d'utiliser les vecteurs propres du laplacien du graphe comme base de Fourier et les valeurs propres comme fréquences du graphe. De cette façon, des valeurs propres plus élevées correspondent à une composante de fréquence plus élevée dans l'analyse de Fourier. En projetant le signal sur ces bases, on peut calculer la réponse en fréquence du graphique du signal à chaque composante.

Laplacien en l'apprentissage automatique. Une des premières utilisations du laplacien se produit dans l'apprentissage semi-supervisé sur des graphes [START_REF] Zhu | Semi-supervised learning with graphs[END_REF][START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. Dans ce problème, on nous donne quelques étiquettes sur les sommets et le but est de déduire les étiquettes des autres sommets en utilisant les étiquettes données et le graphe sous-jacent. Une solution de base donnée par [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] suggère d'utiliser le Laplacien du graphe pour formuler ce problème. Leur principale hypothèse est que la fonction d'étiquetage f est fixée à un sous-ensemble V ⊆ V et lisse sur le reste U = V \ V , i.e. f ne varie pas trop à travers les bords. Pour deux classes d'étiquettes (pour simplifier), cette formulation cherche une fonction de classification f qui n'est connue que sur un petit ensemble de sommets V ⊆ V.

Bibliography

Désignons les étiquettes connues par y ∈ R |V | , puis, [START_REF] Zhu | Semi-supervised learning with graphs[END_REF] suggérons de résoudre le problème contraint suivant :

f = argmin f i∈U j∼i (f (i) -f (j)) 2 s. t. ∀i ∈ V, f (i) = y i ,
En fait, le terme minimisé peut être écrit en termes de laplacien du graphe comme f L U f où f = [f (i)] i∈U . De plus, la solution sous forme fermée implique ce qui suit dans le Laplacien du graphe :

f = (L U ) -1 L U |V f .
Plus tard, cette formulation est considérée dans un cadre non contraint et généralisée dans [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. Dans ces formes également, la solution prend la forme de l'inverse matriciel du laplacien du graphe. Parmi les exemples plus récents en apprentissage automatique qui utilisent largement le laplacien, citons le clustering spectral [Von07; TPGV16], l'intégration de graphes [START_REF] Xu | Understanding graph embedding methods and their applications[END_REF], la sparsification [START_REF] Daniel | Graph sparsification by effective resistances[END_REF] sus aléatoires dans lesquels nous choisissons une forêt étendue au hasard sur un graphe donné. On peut imaginer une infinité d'options pour la distribution d'un tel processus. Dans cette thèse, nous nous concentrons sur la distribution introduite dans [START_REF] Avena | Random forests and networks analysis[END_REF]. Dans cette distribution, la probabilité d'avoir une forêt couvrante est principalement régulée par les poids de ses arêtes et le nombre de ses composants connectés (arbres). Nous nous limitons à cette distribution pour deux raisons principales :

• Les FCAs échantillonnés à partir de cette distribution ont des connexions théoriques riches avec le Laplacien des graphes via certains PPDs,

• Il existe un algorithme efficace pour échantillonner à partir de cette distribution, appelé algorithme de Wilson.

Grâce à ces faits, nous savons que les échantillons de FCAs peuvent être obtenus de manière peu coûteuse, et que leurs propriétés algébriques sont étroitement liées aux propriétés algébriques du graphe original. L'objectif principal de cette thèse est d'exploiter ces faits afin de développer des algorithmes ALA pour l'algèbre numérique basée sur le Laplacien.

Les principales contributions de cette thèse peuvent être énumérées comme suit :

• Nous mettons en lumière les riches connexions des FCAs avec les processus ponctuels déterminant aux et le Laplacien des graphes ; Nous passons en Bibliography revue l'élégante théorie et les propriétés utiles des FCA en tant que PPD. Tout au long de notre tour, alors que pour certaines des propriétés, nous reproduisons les preuves existantes, pour d'autres, nous donnons surtout des preuves algébriques qui sont plus accessibles pour les chercheurs de différents domaines (Voir Théorèmes 2.4.4, 2.4.5).

• En retour, nous développons et analysons des algorithmes ALA basés sur les FCAs pour approximer les solutions d'une grande variété de problèmes impliquant le Laplacien des graphes. Ces problèmes sont :

-Régularisation et interpolation de Tikhonov pour les signaux de graphes : Les problèmes de débruitage, d'élimination des parties bruyantes d'un signal et de complétion des parties manquantes d'un signal sont des problèmes de longue date en traitement du signal. Nous donnons d'abord un cadre d'optimisation unifié pour ces problèmes pour les signaux de graphes i.e. signaux définis sur des sommets. Ensuite, nous donnons de nouveaux algorithmes aléatoires qui approchent la solution de ce cadre.

-Estimation de la trace de l'inverse régularisé des matrices symétriques diagonalement dominantes (SDD) : La trace (somme des entrées diagonales d'une matrice) est une opération centrale en algèbre linéaire. Cependant, le calcul de la trace d'une matrice est un problème difficile lorsque la matrice d'entrée n'est pas directement accessible eg multiplication de grandes matrices ou inverse d'une grande matrice. Nous donnons de nouveaux algorithmes basés sur les FCA pour estimer l'inverse régularisé des matrices SDD (une classe de matrices qui contiennent le Laplacien du graphe). En retour, nos algorithmes sont au moins comparables, et ils surpassent généralement les méthodes de l'état de l'art.

-Estimation des résistances effectives dans les réseaux électriques : La résistance effective est une métrique qui prend sa source dans les réseaux de résistances électriques 1 et ils sont utilisés pour mesurer la similarité d'une paire de noeuds. Ils jouent un rôle central dans de nombreuses applications liées aux graphes, telles que la sparsification, le regroupement ou l'apprentissage des graphes. Cependant, le calcul de cette quantité utile ne s'adapte pas bien à l'augmentation de la taille des graphes, car il faut inverser un laplacien de graphe. Pour éviter ce calcul coûteux, nous donnons des algorithmes basés sur le FCA pour estimer les résistances effectives. Les méthodes proposées sont comparables et même meilleures que les algorithmes de l'état de l'art lorsque le nombre de résistances effectives souhaitées est faible.

-Estimation de divers filtres graphiques: Cet aspect de notre travail établit des liens intéressants entre les FCAs et le filtrage sur graphes i.e. spectral signal filtering for graph signals [START_REF] Tremblay | Design of graph filters and filterbanks[END_REF]. Le filtrage sur graphes est un outil essentiel pour traiter ce type de signaux et a été utilisé dans diverses applications. Pourtant, les calculs associés nécessitent souvent de diagonaliser L, ce qui est coûteux. Dans ce travail, nous donnons un ensemble de filtres graphiques qui peuvent être approximés via des estimateurs basés sur le FCAs en évitant la diagonalisation de L. Ces filtres couvrent ou peuvent imiter en optimisant quelques hyperparamètres certains des filtres fréquemment utilisés, tels que le filtre passe-bas idéal.

• Nous fournissons une analyse biais-variance des algorithmes proposés afin d'énoncer leurs performances théoriques.

• Nous illustrons ces algorithmes dans des applications et des ensembles de données réels tout en comparant ces algorithmes avec les algorithmes existants.

La structure principale de la thèse est la suivante ; nous commençons par donner les outils techniques nécessaires dans le Chapitre 2 sur la théorie des graphes, l'algèbre linéaire aléatoire et les forêts couvrantes aléatoire. Ensuite, au Chapitre 3, nous introduisons les méthodes basées sur les FCAs pour résoudre la régularisation de Tikhonov sur les graphes. Nous comparons ces méthodes aux méthodes existantes et montrons leurs utilisations dans d'autres problèmes liés aux graphes tout en illustrant des applications réelles. Dans le Chapitre 4, nous étendons la gamme des méthodes basées sur les FCAs à un nouvel ensemble de problèmes, notamment l'estimation de la trace de l'inverse régularisée des matrices SDD, les résistances effectives et certains filtres de graphes. Enfin, nous terminons dans le chapitre 5 avec une conclusion générale, une discussion sur les questions ouvertes et les travaux futurs.

Les publications évaluées par les pairs associées à cette thèse peuvent être trouvées dans ce qui suit :

2.3

Given a graph (in the left), PageRank is an algorithm to rank nodes according to their importance in the graph. In the right, we plot the graph by resizing the nodes according their PageRank scores. Bigger nodes have higher scores. . . . . . . . . . . . . . . . . . . . . . . . . .

2.4

A generic path between nodes in R (in blue) . . . . . . . . . . . . . . .

2.5

Given a spanning tree, there is a unique orientation per node that makes the corresponding node the root of the tree. . . . . . . . . . . .

2.6

Loop erasure procedure. We take a random walk starting from the node at the center (in red), apply the loop erasure defined in Definition 2.1.12 2.7 A man shooting at the sun . . . . . . . . . . . . . . . . . . . . . . . . .

2.8

Line fitting between the bullets fired at the sun and temperature in Adana/Turkey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.9

A toy example on the order invariance of the popped cycles. The infinite matrix S contains the initial configuration of the infinite stacks in its columns. As we pop the cycles, the corresponding stacks pop the element at their first row. We follow different orders for popping the cycles but the popped cycles in the are the same ones yielding the same spanning tree as a result. . . . . . . . . . . . . . . . . . . . . . . . . .

2.10

We extend the original graph with an additional node Γ. Then we launch Wilson's algorithm to sample a spanning tree rooted in Γ on the extended graph. Cutting the edges incident to Γ gives a spanning forest distributed by the law of Φ Q . . . . . . . . . . . . . . . . . . . . . . .

2.11

Fixing the partition created by Φ Q ensures that the edges between parts are never sampled. Moreover, every tree sampled in such process has to be a spanning tree of a part. Therefore, given the partition, sampling Φ Q boils down to sampling independent spanning trees in each part. . 2.12 An illustration of the sample-and-solve scheme. First we sample a small portion of the matrix A, then proceed to the expensive operation. . . .

3.1

An Illustration of the filtering properties of g q (λ) for different q values. In this example, we take the temperature dataset collected in [START_REF] Perraudin | Stationary signal processing on graphs[END_REF].

We build a k = 5 nearest neighbour graph from the locations of the temperature sensors over Bretagne/France. The signal is the average temperatures measured at each sensor and depicted by the colours (blue is colder). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2

We complete the missing signal by the solution in Eq. (3.4). . . . . . .

3.4

In this example, we consider the case of graph signal filtering with q = 1. The variance of x then becomes y (I -K 2 )y. We plot the spectrum of the matrix (I -K 2 ) w.r.t. the graph Laplacian eigenvalues/graph frequencies. As a filter, (I -K 2 ) acts like a high-pass filter. Thus the quadratic formula y (I-K 2 )y will take higher values when y is a highly varying signal on the graph. . . . . . . . . . . . . . . . . . . . . . . . .

3.5

An illustration of x. We follow the same first step with x. Then we propagate the (weighted) average signal within each tree. We assume q 1 = • • • = q n = q. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

160

List of Figures

3.6

In this example, we again consider the case of graph signal filtering with q = 1. The variance of x then becomes y (K -K 2 )y. We plot the spectrum of the matrix (K -K 2 ) w.r.t. the graph Laplacian eigenvalues/graph frequencies and compare it with the case of x. As a filter, (K -K 2 ) acts like a band-pass filter i.e. its transfer function only takes high values for a particular range of frequencies, called the pass-band, that are neither high, nor low. Thus the quadratic formula y (K -K 2 )y will take higher values whenever y is mostly composed by the frequency components that are in the pass-band. Also note that in comparison to x, x is always better as illustrated. . . . . . . . . . . . . . . . . . . . .

3.7

We plot the loss function F (x) for x = [x 1 , x 2 ] ∈ R 2 and mark the path followed by the gradient descent algorithm to reach the point that minimizes F (x) (q = 1.0 and α = 0.001). . . . . . . . . . . . . . . . . .

3.9

Approximation and reconstruction error of the algorithms CG (blue and cyan), polynomial approximation (orange) and the forest estimators x (light green) and z (dark green) for solving graph signal denoising problem. The dark dashed line indicates the time taken by the direct solution (backslash operator in Julia). In the reconstruction error plots, we also add the error of the initial measurements y in green dashed line and the error of the exact solution x in red line. . . . . . . . . . .

3.10

The accuracy vs number of labeled vertices per class. We give the classification accuracy of gSSL (the first column) and LP (second column) and their forest estimates over the datasets Cora, Citeseer and Pubmed.

In each plot, the performance of the exact solution is given in the black line. For gSSL, we compute x over 50 (light green) and 500 (dark green) forest realization. In the case of LP, the number of realization becomes 5 (light green) and 50 (dark green). We also add the performance of a classifier that returns a label at random (random classifier) and returns the same label all the time (constant classifier). . . . . . . . . . . . . .

List of Figures

constructs trees. Without loss of generality, each newly added node is connected to the one added in the previous step. We denote the edge set at k-step as E k . To show that the statement is true for an inductive step, at first, we assume that it holds for k-th step. Under this assumption, proving that it also holds for k + 1-th step confirms that it holds for the induction step and finishes the first part. By applying the block matrix determinant formula, we have:

| det(B E k+1 |-1 )| = w(k + 1, k + 2)| det B E k |-1 | (A.2) Given | det B E k |-1 | = (i,j)∈E k w(i, j) 1/2
, one easily obtains:

| det B E k+1 |-1 | = w(k + 1, k + 2)   (i,j)∈E k w(i, j)   1/2 =   (i,j)∈E k+1 w(i, j)   1/2
which proves the statement for the induction step and completes the first part.

Then rearranging the summation yields:

R 2 i,j .
Upper bound oni,j∈V R i,j Cov( Xi,j , C). Notice that both Xi,j and C are nonnegative random variables. Then the following inequality holds: Upper bound on i,j∈V R 2 i,j Var(C). We again plug the definition of the variance and it yields:

i,j∈V R 2 i,j Var(C) = i,j∈V R 2 i,j E[C 2 ] -c 2 = i,j∈V R 2 i,j   1 (n -1) 2 E[

( a )

 a Traffic data collected by a network of sensors in Grenoble [WMOKBB] (b) Twitter friendship network of the authors of the tweets that the user @PilavciYigit liked until 29/06/2022. (c) A partial map of Internet (d) Temperature sensor network over Bretagne/France [PV17]. The red colour depicts higher average temperatures.

Fig

  Fig. 1.1.: Real-life examples of graphs

  Fig. 2.1.: One can obtain completely different graphs for the same degree distribution. First, we generate a graph (b) from Barabasi-Albert model [AB02] (a statistical model for generating a random graph) with n = 1000 and the model parameter k = 1.Then, we pass its degree distribution given in (a) to the Kalisky model[START_REF] Kalisky | Tomography and stability of complex networks[END_REF] (a statistical model for generating a graph given a degree distribution) to generate another graph (c) that has the same degree distribution.

2

 2 Proposition 2.1.2 (Poincaré for weighted graphs). Consider a weighted graph G = (V, E, w) with an edge incidence matrix B. If S ⊆ E with |S| = |V| -1 is the edges of a spanning tree, then, | det B S|-m | = (i,j)∈S w(i, j) 1/2 for all m ∈ V. Otherwise, | det B S|-m | = 0. Proof. See App. A.2. Developing on this theorem gives the following extension on spanning forests: Theorem 2.1.3 (Poincaré for forests). Given a graph G = (V, E, w) with an edge incidence matrix B, consider a vertex subset R ⊆ V and an edge subset S ⊆ E with |S| = |V| -|R|. If S is an edge subset such that • S forms a spanning forest with |R| trees, • S does not form any path between the nodes of R (for rooted forests, R is the root set.), 18 Chapter Background then, | det B S|-R | =

Fig. 2 .

 2 Fig. 2.4.: A generic path between nodes in R (in blue)

Fig. 2

 2 Fig. 2.5.:Given a spanning tree, there is a unique orientation per node that makes the corresponding node the root of the tree.

  Fig. 2.8.: Line fitting between the bullets fired at the sun and temperature in Adana/Turkey

  the residual is less than the tolerance, stop the loop

Fig. 2 .

 2 Fig. 2.10.:We extend the original graph with an additional node Γ. Then we launch Wilson's algorithm to sample a spanning tree rooted in Γ on the extended graph. Cutting the edges incident to Γ gives a spanning forest distributed by the law of Φ Q

Problem 3.1. 1 .

 1 Given a graph G = (V, E, w) and measurements y = [y 1 , . . . , y p ] ∈ R p over the vertices V ⊆ V with |V | = p, we attempt to recover the original signal by solving: x = argmin x∈R n ||My -x|| 2 Q + x Lx. (3.1)

( a )

 a Fig. 3.1.:An Illustration of the filtering properties of g q (λ) for different q values. In this example, we take the temperature dataset collected in[START_REF] Perraudin | Stationary signal processing on graphs[END_REF]. We build a k = 5 nearest neighbour graph from the locations of the temperature sensors over Bretagne/France. The signal is the average temperatures measured at each sensor and depicted by the colours (blue is colder).

  Fig. 3.2.:We complete the missing signal by the solution in Eq. (3.4).

(Chapter 3

 3 Expectation and Variance of x). Two main properties of x are that: 66 Graph Tikhonov Regularization and Interpolation with RSFs

  Fig. 3.3.: An illustration of x. Given the graph and the signals (as colours on the nodes), wefirst sample a rooted spanning forest via Wilson's algorithm. Then we propagate the signal at the roots through the corresponding trees.

  Fig. 3.4.:In this example, we consider the case of graph signal filtering with q = 1. The variance of x then becomes y (I -K 2 )y. We plot the spectrum of the matrix (I -K 2 ) w.r.t. the graph Laplacian eigenvalues/graph frequencies. As a filter, (I -K 2 ) acts like a high-pass filter. Thus the quadratic formula y (I -K 2 )y will take higher values when y is a highly varying signal on the graph.

3. 3

 3 Estimating x via RSFs

  Fig. 3.7.:We plot the loss function F (x) for x = [x 1 , x 2 ] ∈ R 2 and mark the path followed by the gradient descent algorithm to reach the point that minimizes F (x) (q = 1.0 and α = 0.001).

Fig. 3 .

 3 Fig. 3.8.: Empirical mean squared error of x (orange horizontal line) and z (blue parabola) w.r.t α on two graphs generated by random models. On the left is a random regular graph with n = 1000 and m = 10000. On the right is a Barabasi-Albert model with parameter k = 10 resulting in n = 1000 and m = 9900. The green (resp. red) vertical dashed line shows α = 2q q+2dmax(resp. the estimated α from the samples). The blue dot represents the best possible variance reduction obtained for α = α . The number of Monte Carlo samples is set to N = 10, and the error results are averaged over 200 realizations. The signal is a random vector generated from N (0, I).

Fig. 3 .

 3 Fig. 3.10.:The accuracy vs number of labeled vertices per class. We give the classification accuracy of gSSL (the first column) and LP (second column) and their forest estimates over the datasets Cora, Citeseer and Pubmed. In each plot, the performance of the exact solution is given in the black line. For gSSL, we compute x over 50 (light green) and 500 (dark green) forest realization. In the case of LP, the number of realization becomes 5 (light green) and 50 (dark green). We also add the performance of a classifier that returns a label at random (random classifier) and returns the same label all the time (constant classifier).

Fig. 3 .

 3 Fig. 3.12.:The images denoised by Newton's method where the updates calculated byx (left) and x (middle). On the right, we plot the loss function through the iterations of Newton's method. We assume a 2D-grid graph where n = 128 × 128. The hyperparameter is set to 0.1, and α is selected at every iteration by approximate line search algorithm. For the forest-based updates, we sample N = 40 forests. The peak signal-to-noise ratios on the final images are found as 24.84 for both images.

( a )

 a Fig. 3.13.:We illustrate ADMM algorithm for image denoising application. We are given a noisy measurements y = x + where the original image x is unknown and the noise is Gaussian i.e. ∼ N (0, σ = 0.2). We denoise the image by solving Eq (3.21) and the regularization parameter q is found as 3.3 by a grid search. In (c) and (d), we give the solutions by ADMM the algorithm where the updates are calculated via the exact and forest-based solutions. In both cases, ADMM is terminated after 500 iterations and ρ = 0.2. The forest updates are averaged over 10 forests realizations per iteration. Finally, in (e), we plot the objective function across the iterations.

Contents 4 . 1

 41 Estimation of the Trace of the Regularized Inverse of the Laplacian 95 4.1.1 Variance Reduction for Trace Estimation . . . . . . . . . 97 4.1.2 Empirical Results on Trace Estimation . . . . . . . . . . . 102 4.2 Estimators for Effective Resistances . . . . . . . . . . . . . . . . 104 4.2.1 A Crash Course on Effective Resistances . . . . . . . . . 104 4.2.2 Computing Effective Resistances . . . . . . . . . . . . . . 107 4.2.3 Estimating ER via 2-Rooted Forests . . . . . . . . . . . . 108 94 4.2.4 Estimating ER via Local Forests . . . . . . . . . . . . . . 114 4.2.5 Empirical Results on ER Estimation . . . . . . . . . . . . 116 4.3 Graph Filtering with RSFs . . . . . . . . . . . . . . . . . . . . . . 119 4.3.1 Solving Lx = y with RSFs . . . . . . . . . . . . . . . . . 120 4.3.2 Filters on the Duplicated Graph . . . . . . . . . . . . . . 122 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.1 Estimation of the Trace of the Regularized Inverse of the Laplacian Reconsider the simplified version of the denoising problem 3.1.1:

Fig. 4. 1 .

 1 Fig. 4.1.:The edges between trees in a spanning forest (in bold). One needs to track these edges to compute the control variate introduced in Eq. 4.1. For computing the control variate induced by S, one only needs the edges between trees that are incident to the roots. In this example, only edge that fits this description is (4, 6).

  Fig. 4.2.: An example of the stack representation (left) and the graph induced by the first layer (right). The blue nodes indicate the roots sampled at the first layer.

  Fig. 4.3.: Effective Runtime vs tr(K)/n.

•

  Citeseer: A citation network with n = 2110 nodes and m = 3668 and edges, • Cora: A citation network with n = 2485 nodes and m = 5069 edges, • Pubmed: A citation network with n = 19717 nodes and m = 44324 edges, • Collab-CM: A collaboration network with n = 21363 and m = 91342.

  Fig. 4.4.: Relative error vs. runtime of the global (ST and RP) and local (LF, TP, MC2)algorithms for estimating effective resistances. For the global algorithms, we run them by varying their iteration parameter and measure their average relative error over all edges and total run-time. For the local algorithms, we measure their average relative error and run-time for a single edge which is given in their plots as the leftest marker. We compute the time taken by each local algorithms for estimating k effective resistances by multiplying their run-time for a single estimation by k. By varying k between 1 and m, we generate the points in the plots of the local algorithms.

  Fig. 4.5.: An illustration of the averaging operation in x 2rf. After averaged in one tree, the signal is propagated in the other tree.

  Fig. 4.6.:We approximate two target filters g k (λ) and g p,q (λ) via g q and f Θ . For f Θ , we choose h(L) = d max I -L. Then we optimize the parameters of the candidate filters g q and f Θ for estimating the filters g k=200 (λ) and g p=2,q=10 (λ). In the first figure, we find the parameter q = 2.14 for g q and Θ = (0.61, 1.70, 8.02, 18.16) for f Θ . As you can see f Θ yields much more approximate filter when these parameters plugged compared to that of g q . The similar results are also observed in the case of approximating g k=200 (λ).The filter with four parameter f Θ where Θ = (0.25, 2.93, 10.80, 13.88) have more flexibility compared to g q=1.63 , thus it yields a better approximation.

  Fig. 5.1.: An example of the stack representation (left) and the induced graphs from the first and the second layer. Blue nodes are the current root set.

( a )

 a Fig. 5.1.: Exemples de graphes dans la vie réelle

3. 3

 3 An illustration of x. Given the graph and the signals (as colours on the nodes), we first sample a rooted spanning forest via Wilson's algorithm. Then we propagate the signal at the roots through the corresponding trees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 8

 8 Empirical mean squared error of x (orange horizontal line) and z (blue parabola) w.r.t α on two graphs generated by random models. On the left is a random regular graph with n = 1000 and m = 10000. On the right is a Barabasi-Albert model with parameter k = 10 resulting in n = 1000 and m = 9900. The green (resp. red) vertical dashed line shows α = 2q q+2dmax (resp. the estimated α from the samples). The blue dot represents the best possible variance reduction obtained for α = α . The number of Monte Carlo samples is set to N = 10, and the error results are averaged over 200 realizations. The signal is a random vector generated from N (0, I). . . . . . . . . . . . . . . . . . . . . . .

  Fig. A.1.: A weighted graph and its reduced edge incidence matrix

-

  R i,j Cov( Xi,j , C) =i,j∈V R i,j (E[C Xi,j ] -E[C]E[ Xi,j ]) ≤ i,j∈V R i,j E[C]E[ Xi,j ] = c 2 i,j∈V R 2 i,j .

  (x,y)∈E (a,b)∈E Xx,y , Xa,b ] -c 2 the term within the expectation, one obtains: (x,y)∈E (i,j)∈E Xx,y , Xa,b = (x,y)∈E (a,b)∈E I(r Φq (a) = r Φq (b), r Φq (x) = r Φq (y)|ρ(Φ q )| = 2) |V t(Φq,a) ||V t(Φq,b) ||V t(Φq,x) ||V t(Φq,y) | .

•

  Yusuf Yigit Pilavcı et al. "Variance Reduction for Inverse Trace Estimation via Random Spanning Forests". In: GRETSI 2022 -XXVIIIème Colloque Francophone de Traitement du Signal et des Images. Nancy, France, Sept. 2022 • Yusuf Yigit Pilavcı et al. "Variance reduction in stochastic methods for largescale regularised least-squares problems". In: arXiv preprint arXiv:2110.07894

	(2021)
	• Yusuf Yigit Pilavcı et al. "Graph tikhonov regularization and interpolation via
	random spanning forests". In: IEEE transactions on Signal and Information
	Processing over Networks 7 (2021), pp. 359-374
	• Yusuf Yigit Pilavcı et al. "Smoothing graph signals via random spanning forests".
	In: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and
	Signal Processing (ICASSP). IEEE. 2020, pp. 5630-5634

  Deterministic Methods . . . . . . . . . . . . . . . . . . . 33 2.3.4 Randomized Algorithms . . . . . . . . . . . . . . . . . . 37 2.4 Random Spanning Forests . . . . . . . . . . . . . . . . . . . . . . 39 2.4.1 Wilson's algorithm . . . . . . . . . . . . . . . . . . . . . 40 2.4.2 DPPs in Random Spanning Forests . . . . . . . . . . . . 47 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.5.1 Graph Laplacian . . . . . . . . . . . . . . . . . . . . . . . 55 2.5.2 Determinantal Point Processes and Randomized Linear Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.5.3 Random Spanning Forests . . . . . . . . . . . . . . . . . 58

	2.3.3
	2.3.1 Linear Least-Squares Problem . . . . . . . . . . . . . . . 31
	2.3.2 Approximate Methods for Solving LS Problem . . . . . . 33

Contents 2.1 Graph Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1 Graph Matrices . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.2 Random Walks on Graphs . . . . . . . . . . . . . . . . . 21 2.2 Determinantal Point Processes . . . . . . . . . . . . . . . . . . . 25 2.2.1 Basics of DPPs . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Randomized Linear Algebra . . . . . . . . . . . . . . . . . . . . . 30

Definition 2.1.12 (

  LERW). Given a random walk W = (e 1 , e 2 , . . . , e l ), let us denote the nodes visited in W by V = (v 1 , v 2 , . . . , v l+1 ) such that ∀e i ∈ W, e i =

	(a) Before loop-erasure	(b) After loop-erasure
	Fig. 2.6.: Loop erasure procedure. We take a random walk starting from the node at the center
	(in red), apply the loop erasure defined in Definition 2.1.12

  D is a diagonal matrix which contains the diagonal entries of B and L and U contains respectively the lower and upper triangular parts of B. As D is easy to invert, one only needs to access the matrix B via matrix-vector products per iteration. Other basic iterative methods slightly differs in how they decompose the matrix B. As a result, they vary in the computational load per iteration and convergences properties. However, one can still derive general results on the convergence:

	Theorem 2.3.1 (Theorem 4.1 in [Saa03]). Define the iteration matrix G = M -1 N
	for any splitting in form of B = M -N. Assume G is non-singular, thus there exists
	a unique solution x . Then the iteration in Eq. (2.14) converges to the solution
	x if and only if the spectral radius of G i.e. max i |λ i (G)| where λ i (G) is the i-th
	eigenvalues of G, denoted by ρ(G) is strictly less than 1.

  (The roots and edges of Φ

		Theorem 2.4.3
	12:	u ← Next[u]
	13:	end if
	14:	end while
	15:
	18:	InForest[u] ← true
	19:	u ← Next[u]
	20:	end while
	21: end for
	22: return Next

of u to null 10: else # If false, continue the random walk 11: Next[u] ← RandomSuccessor(u, G) u ← i # Go back to the initial node 16: # Add the newly created branch to the forest 17: while not InForest[u] do

  RSF-based Quasi Newton's Method . . . . . . . . . . . . 87 3.4.3 L-1 Regularization on Graphs . . . . . . . . . . . . . . . 90 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

	3.4.2
	3.3.1 Variance Reduction via Conditional Monte Carlo . . . . . 69
	3.3.2 Variance Reduction via Control Variates . . . . . . . . . . 73
	3.3.3 Empirical Comparisons with SOTA Algorithms . . . . . . 78
	3.4 Several Use Cases of the RSF-based Estimators . . . . . . . . . . 82
	3.4.1 Semi-supervised Learning on Graphs . . . . . . . . . . . 82
	60

Contents 3.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2 State-of-the-Art . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3 Estimating x via RSFs . . . . . . . . . . . . . . . . . . . . . . . . 66

  table summarizes their properties: The first two datasets Cora and Citeseer are not connected graphs. Therefore, in those graphs, we only keep the largest connected component and eliminate the rest.

		Tab. 3.1.: Benchmark Graph datasets
		Dataset #Nodes #Edges #Classes
		Citeseer	2110	3668	6
		Cora	2485	5069	7
		Pubmed 19717	44324	3
		Graph signals. In this experimentation, we consider the graph signal denoising
		problem over various graphs, in which we are given the noisy measurements y = x+
	78	Chapter 3 Graph Tikhonov Regularization and Interpolation with RSFs

  The solution to this optimization scheme is a smooth function over the graph, which still preserves the original labelling information encoded in Y. In other words, Fgives a diffused version of the given labels in Y through the neighborhoods of the graph. As in Prob. 3.1.1, one has the closed form solution as follows:

.17)

82

Chapter 3 Graph Tikhonov Regularization and Interpolation with RSFs

  Depuis leur première apparition pour modeliser les ponts de Königsberg, les graphes, c'est-à-dire un ensemble de sommets et d'arêtes, ont été les objets de modélisation et de résolution de nombreux problèmes en physique, chimie, biologie, informatique, sciences sociales et bien d'autres encore. Ils sont si ubiquitaires car ils sont les modèles naturels à utiliser lorsque les données représentent des relations mutuelles entre des éléments individuels. Aujourd'hui, le traitement des données structurées en graphes est l'une des principales questions de recherche en traitement du signal et en apprentissage automatique. De nombreuses recherches se concentrent sur l'exploitation de la structure des graphes pour accomplir des tâches avancées, comme la gestion et l'analyse d'Internet (un réseau informatique colossal), la compréhension du cerveau humain[START_REF] Farzad V Farahani | Application of graph theory for identifying connectivity patterns in human brain networks: a systematic review[END_REF], la découverte de médicaments[START_REF] Jiang | Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models[END_REF], Dans une configuration typique, chaque emplacement correspond à un noeud, et deux noeuds sont reliés par une arête lorsqu'une route les relie. Les signaux sur ces réseaux sont généralement des caractéristiques du trafic, telles que la vitesse ou le volume du trafic. • Réseaux sociaux [VMS21] : Dans de nombreuses études sociales, outre les informations individuelles sur les sujets, la structure de leurs interactions est également d'une importance capitale. Les réseaux de médias sociaux tels que Facebook et Twitter sont probablement les exemples les plus connus. Dans ces cas, chaque utilisateur est considéré comme un noeud et deux noeuds sont connectés s'il existe une interaction entre les utilisateurs correspondants, comme des amis sur Facebook ou des followers sur Twitter. Compte tenu de cette structure de graphe/réseau, toute caractéristique d'un utilisateur, telle que ses goûts, peut être analysée comme un signal de graphe.

l'analyse des réseaux sociaux

[START_REF] Otte | Social network analysis: a powerful strategy, also for the information sciences[END_REF]

ou la prévision/planification du trafic à grande échelle

[START_REF] Hasanzadeh | A graph signal processing approach for real-time traffic prediction in transportation networks[END_REF]

. D'autres exemples apparaissent dans le traitement du langage naturel

[START_REF] Nastase | A survey of graphs in natural language processing[END_REF]

, la finance

[START_REF] Vinıécius De | Learning undirected graphs in financial markets[END_REF]

, les réseaux alimentaires

[START_REF] Dunne | Network structure and robustness of marine food webs[END_REF]

, les réseaux électriques

[START_REF] Giuliano | The power grid as a complex network: a survey[END_REF]

et bien d'autres encore. Dans bon nombre de ces exemples et d'autres, les graphes de la vie réelle sont accompagnés de certaines caractéristiques/signaux, également appelés signaux sur graphe, sur l'ensemble des sommets. En voici quelques exemples :

• Réseaux de transport

[START_REF] Hasanzadeh | Piecewise stationary modeling of random processes over graphs with an application to traffic prediction[END_REF]

: Les réseaux routiers en sont un exemple frappant.

  • Réseaux de température [PV17] : Un réseau de température est un réseau de capteurs de température sur une région dans laquelle chaque capteur est un noeud. Ici, les arêtes peuvent être établies selon différents critères. La pratique courante est de construire le graphe des plus proches voisins à partir des proximités des capteurs. Le signal, dans ce cas, est naturellement les mesures de température recueillies par les capteurs. • Réseaux neurologiques [HBMBRV17] : Compte tenu des nombreuses questions restées sans réponse, l'analyse du cerveau est un sujet passionnant. L'étude basée sur les graphes, en fait, donne une façon naturelle d'analyser les signaux du cerveau. Dans ce type d'analyse, on considère souvent les différentes régions du cerveau comme des noeuds et les arêtes sont établies en fonction de leur proximité structurelle ou de propriétés fonctionnelles telles que la corrélation dans leurs activités. L'analyse des signaux neurologiques sur un tel réseau permet de déchiffrer des informations intéressantes sur le cerveau.

		Compte tenu du grand volume et de la diversité de ces ensembles de données,
		le développement d'outils appropriés pour les traiter et les analyser est devenu
		problématique. En fait, la recherche en apprentissage automatique et en traite-
		ment du signal pour les graphes et les signaux de graphes a émergé au cours des
		dernières décennies [OFKMV18; RBTGV19]. À leur tour, ces outils ont été util-
		isés pour résoudre des problèmes de la vie réelle allant de la détection de fausses
		nouvelles dans les médias sociaux [MFEMB19] à l'analyse de la propagation pour
		COVID-19 [PNV21], du décodage des signaux du cerveau [OTLFPG22] aux prévi-
		sions météorologiques [Kei22], etc. Dans la plupart de ces outils et méthodes,
		les représentations matricielles des graphes, notamment le laplacien des graphes,
		revêtent une importance significative.
	Un graphe path	Laplacien de Graphe De nombreux outils d'analyse extraient des informations des
	Son Laplacian L

• Autres réseaux biologiques [GDJSRLHVRT+21; CPFSC21; YB20] : Le cerveau n'est pas seulement le sujet biologique de l'analyse signal sur graphes. Dans la littérature, nous trouvons de nombreuses études liées à la découverte de médicaments, à l'analyse des structures protéiques ou des interactions génétiques.

  En regardant simplement la multiplicité des valeurs propres nulles, on peut compter le nombre de composantes connectées. De plus, la deuxième plus petite valeur propre, également appelée écart spectral, donne une mesure de la connectivité du graphe, qui est principalement utilisée dans les applications de stabilité et de robustesse des réseaux dynamiques. De même, le vecteur propre correspondant, également connu sous le nom de vecteur de Fiedler, donne la plus petite partition du graphe, celle qui laisse le plus petit nombre d'arêtes entre les parties. Le dessin spectral de graphes est un autre problème qui utilise la décomposition des vecteurs propres du laplacien. Dans ce problème, on cherche un dessin (une correspondance des sommets aux coordonnées euclidiennes) d'un graphe qui minimise la distance entre les sommets qui sont proches les uns des autres dans le graphe. Il s'avère que certains vecteurs propres du Laplacien sont les solutions analytiques de ce problème. Plus proche des thèmes principaux de cette thèse, nous trouvons un lien intéressant entre les valeurs propres de L et l'énumération d'un certain type de sous-graphes. En particulier, le célèbre théorème

	Un arbre cou-
	vrant. Ses arêtes
	sont indiquées
	en rouge, et les
	lignes	pointil-
	lées indiquent
	les arêtes du
	graphe.

  et l'apprentissage profond sur les graphes[START_REF] Wu | A comprehensive survey on graph neural networks[END_REF]. Dans toutes ces applications et bien d'autres, l'analyse de L est de la plus haute importance.Algèbre Linéaire Aléatoire (ALA). En bref, l'ALA est un ensemble d'outils qui évitent le calcul direct d'une opération algébrique par le biais de la randomisation. Le principal avantage de ces outils est qu'ils permettent d'approximer des opérations très coûteuses, telles que l'inversion de matrice et la décomposition de la valeur singulière, de manière très économique en prenant uniquement des échantillons aléatoires de la matrice d'entrée. Plus de détails et d'exemples peuvent être trouvés dans[START_REF] Drineas | RandNLA: randomized numerical linear algebra[END_REF]. Plus proche des thèmes de cette thèse,[START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] montrent que les processus ponctuel déterminantaux (PPD), i.e. un processus ponctuel aléatoires avec de nombreuses propriétés tractables, peuvent être facilement adaptés comme algorithmes ALA afin de résoudre des systèmes linéaires impliquant un large ensemble de matrices d'entrée couvrant le Laplacien du graphe. Inspirés par ces exemples et bien d'autres, nous nous concentrons sur un processus ponctuel déterminantal particulier qui nous permet de développer des outils randomisés pour l'algèbre numérique basée sur le Laplacien. Une forêt dans un graphe est un sousgraphe qui ne contient aucun cycle. Elle est dite couvrante si elle contient tous les sommets du graphe. Enfin, les forêts couvrantes aléatoires sont les proces-

	Les Forêts Couvrantes Aléatoire (FCA) Une forêt
	d'envergure avec
	trois arbres. Les
	arêtes de la forêt
	sont en rouge,
	et le graphe
	est représenté
	par des lignes
	pointillées.

•

  Yusuf Yigit Pilavcı et al. "Variance Reduction for Inverse Trace Estimation via Random Spanning Forests". In: GRETSI 2022 -XXVIIIème Colloque Francophone de Traitement du Signal et des Images. Nancy, France, Sept. 2022 Bibliography • Yusuf Yigit Pilavcı et al. "Variance reduction in stochastic methods for largescale regularised least-squares problems". In: arXiv preprint arXiv:2110.07894 (2021) • Yusuf Yigit Pilavcı et al. "Graph tikhonov regularization and interpolation via random spanning forests". In: IEEE transactions on Signal and Information Processing over Networks 7 (2021), pp. 359-374 • Yusuf Yigit Pilavcı et al. "Smoothing graph signals via random spanning forests". In: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE. 2020, pp. 5630-5634 De plus, certains matériaux sont présentés dans ma thèse de maîtrise [PIL19]. Real-life examples of graphs . . . . . . . . . . . . . . . . . . . . . . . . 2.1 One can obtain completely different graphs for the same degree distribution. First, we generate a graph (b) from Barabasi-Albert model [AB02] (a statistical model for generating a random graph) with n = 1000 and the model parameter k = 1. Then, we pass its degree distribution given in (a) to the Kalisky model [KCbH04] (a statistical model for generating a graph given a degree distribution) to generate another graph (c) that has the same degree distribution. . . . . . . . . . . . . . . . . . . . . . 2.2 Given a graph in (a), this figure illustrates a subgraph (b), a walk (c), a cycle (d), a path (e), a spanning tree (f), a rooted spanning tree (g) and a rooted spanning forest (h). The blue nodes in (g) and (h) indicates the root nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

		List of Figures
		1.1
	158	Bibliography

  Φq (i) = r Φq (j), |ρ(Φ q )| = 2) |V t(Φq,i) | 2 |V t(Φq,j) | 2Notice that, the sum i,j∈VI(r Φq (i) =r Φq (j),|ρ(Φq)|=2) |V t(Φq ,i) | 2 |V t(Φq ,j) | 2is only non-zero over the spanning forests with exactly two trees. Moreover, given an arbitrary forest φ with two roots v 1 , v 2 , it always sums up to:i,j∈V I(r φ (i) = r φ (j)) |V t(φ,v 1 ) | 2 |V t(φ,v 2 ) | 2 = 2 |V t(φ,v 1 ) ||V t(φ,v 2 )| . Recall that |V t(φ,v 1 ) | + |V t(φ,v 2 )| = nfor any graphs and both terms are positive integers. Thus, the minimum of |V t(φ,v 1 ) ||V t(φ,v 2 ) | equals to n -1. As a result, we obtain the following upper-bound:

	Var( Xi,j ) = E	 	I(r   -	c 2 R 2 i,j
	i,j∈V	i,j∈V			i,j∈V
	i,j∈V	Var( Xi,j ) ≤	2 n -1	-c 2	i,j∈V

Graph representation of an electrical resistor networks takes every edge as a resistor with a unit conductance. The conductance equals the weight of the edge for weighted graphs.

Undirected graphs are often defined as the graphs in which all edges are bidirectional. Our definition here is equivalent to this and more generalized by using the edge weights

2.1 Graph Theory

Chapter 2 Background

2.1 Graph Theory

2.2 Determinantal Point Processes

The contrary is not necessarily true. See[START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF] 2.2 Determinantal Point Processes

2.3 Randomized Linear Algebra

Except for the initialization, they do not necessarily include any sort of randomization.2.3 Randomized Linear Algebra

2.4 Random Spanning Forests

This is, because firstly rank L = rank B ext :|-Γ and rank B ext :|-Γ ≥ n due to its n × n diagonal matrix component Q. Also B ext :|-Γ is a n × (m + n) long matrix, one has rank L ≤ n which implies the equality rank L = n.

It is worth to note that to best of our knowledge, this is a new proof compared to the existing probabilistic and combinatorial proofs [ACGM18; CK78; CA02]2.5 Conclusion

Julia implementation https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/

Julia implementation https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl

3.3 Estimating x via RSFs

3.4 Several Use Cases of the RSF-based Estimators

4.1 Estimation of the Trace of the Regularized Inverse of the Laplacian

This is also true for the stratified sampling with the proportional allocation i.e. Ni ∝ P(|ρ (Φq)| ∈ Ci) for all i

The implementation can be found in https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl

The implementation can be found in https://docs.juliahub.com/KrylovMethods 102 Chapter

RSF Estimation of Some Important Graph Quantities

The real-life data sets can be found in https://snap.stanford.edu/data/ 4.1 Estimation of the Trace of the Regularized Inverse of the Laplacian

wmax and wmin are maximum and minimum edge weights in the graph.

The soft big-O notation Õ hides the term log c n

In[START_REF] Peng | Local Algorithms for Estimating Effective Resistance[END_REF], the authors suggest using the local algorithms for a small number of vertex pairs 108 Chapter 4 RSF Estimation of Some Important Graph Quantities

Here, we compute the sample mean for X and C separately from i.i.d. samples and compute the ratio.4.2 Estimators for Effective Resistances

Note that L † (y + µ1) = L † y for any µ ∈ R. We choose to assume y is zero-mean.4.3 Graph Filtering with RSFs

Chapter 4 RSF Estimation of Some Important Graph Quantities

4.4 Conclusion

Bibliography

Algèbre numérique basée sur le Laplacien. Comme indiqué ci-dessus, le nombre d'applications impliquant le Laplacien des graphes est grand. Dans beaucoup de ces applications, la solution nécessite le calcul d'une décomposition inverse (régularisée ou pseudo) ou d'une décomposition de l'indice du Laplacien. Cependant, le calcul direct de ces quantités ne s'adapte pas bien à la taille du graphe et devient même impossible pour les très grands graphes. En conséquence, de nombreux chercheurs se sont attachés à développer des outils d'algèbre numérique spécifiques au laplacien. Les plus importants sont étroitement liés à la théorie spectrale des graphes[START_REF] Spielman | Spectral graph theory[END_REF] et à la résolution de systèmes linéaires laplaciens[START_REF] Nisheeth | Lx= b[END_REF]. Ces études donnent une riche collection d'outils algébriques qui évitent le calcul direct coûteux et, à la place, approchent la solution requise dans certaines plages d'erreur. Ce faisant, elles tirent parti de certaines propriétés communes du laplacien des graphes, telles que leurs propriétés algébriques, leurs liens étroits avec la théorie des graphes ou la rareté des graphes réels. Par rapport à ces travaux, nous prenons dans cette thèse une voie alternative en utilisant la randomisation afin de développer des outils algébriques spécifiques au Laplacien. Nous approchons la solution requise dans certains de ces problèmes en exploitant les riches connexions théoriques entre les Laplaciens et un processus aléatoire sur les graphes.

La représentation graphique d'un réseau de résistances électriques prend chaque bord comme une résistance avec une conductance unitaire. La conductance est égale au poids de l'arête pour les graphes pondérés.156Bibliography

In fact, this is only possible on a graph with two nodes and one edge. Otherwise, there must be at least one edge that connects nodes with the same root.A.5 Upper bounds on the terms in Eq. (4.6)

3.11

The original (left) and noisy image (right). The peak signal-to-noise ratio for the noisy image is 17.13 . . . . . . . . . . . . . . . . . . . . .

3.12

The images denoised by Newton's method where the updates calculated by x (left) and x (middle). On the right, we plot the loss function through the iterations of Newton's method. We assume a 2D-grid graph where n = 128×128. The hyperparameter is set to 0.1, and α is selected at every iteration by approximate line search algorithm. For the forestbased updates, we sample N = 40 forests. The peak signal-to-noise ratios on the final images are found as 24.84 for both images. . . . . .

3.13

We illustrate ADMM algorithm for image denoising application. We are given a noisy measurements y = x + where the original image x is unknown and the noise is Gaussian i.e. ∼ N (0, σ = 0.2). We denoise the image by solving Eq (3.21) and the regularization parameter q is found as 3.3 by a grid search. In (c) and (d), we give the solutions by ADMM the algorithm where the updates are calculated via the exact and forest-based solutions. In both cases, ADMM is terminated after 500 iterations and ρ = 0.2. The forest updates are averaged over 10 forests realizations per iteration. Finally, in (e), we plot the objective function across the iterations. . . . . . . . . . . . . . . . . . . . . . . .

4.1

The edges between trees in a spanning forest (in bold). One needs to track these edges to compute the control variate introduced in Eq. 4.1.

For computing the control variate induced by S, one only needs the edges between trees that are incident to the roots. ) algorithms for estimating effective resistances. For the global algorithms, we run them by varying their iteration parameter and measure their average relative error over all edges and total run-time. For the local algorithms, we measure their average relative error and runtime for a single edge which is given in their plots as the leftest marker.

We compute the time taken by each local algorithms for estimating k effective resistances by multiplying their run-time for a single estimation by k. By varying k between 1 and m, we generate the points in the plots of the local algorithms. . . . . . . . . . . . . . . . . . . . . . . . .

List of Figures

4.5

An illustration of the averaging operation in x 2rf . After averaged in one tree, the signal is propagated in the other tree. . . . . . . . . . . .

4.6

We approximate two target filters g k (λ) and g p,q (λ) via g q and f Θ . For f Θ , we choose h(L) = d max I -L. Then we optimize the parameters of the candidate filters g q and f Θ for estimating the filters g k=200 (λ) and g p=2,q=10 (λ). In the first figure, we find the parameter q = 2.14 for g q and Θ = (0.61, 1.70, 8.02, 18.16) for f Θ . As you can see f Θ yields much more approximate filter when these parameters plugged compared to that of g q . The similar results are also observed in the case of approximating g k=200 (λ). Here we give an upper-bound on the trace tr((Q + L) -1 (Q + D)) which gives the expected number of steps taken by Wilson's algorithm. First, we write this trace explicitly as follows:

Then, by Theorem 2.4.4, we know that

By plugging this upper bound on K i,i 's, one obtains:

A.2 Proof of Proposition 2.1.2

The proof of this theorem is two-fold:

• In the first part, we show

, ∀m ∈ V in case that S forms a spanning tree.

• In the second, we prove that

The first part is proven by induction in which we will show that the statement holds for a base case and an inductive step. We consider two nodes connected with an edge e as the base case. The statement is true for the base case since

where m = 1 without loss of generality. This leaves us to examine the inductive step. At each inductive step, we add a node and connect it to one of the others (See Fig. A.1). Note that this inductive step can build any possible tree and it only 1 2 3 4

Fig. A.2.: A cycle and its reduced edge incidence matrix

In the second part, we examine the case that S includes a cycle. Imagine that we build a circuit in a similar fashion. Then, whenever we construct a cycle, the corresponding reduced edge incidence matrix contains a singular principal block (See Fig . A.2). Thus, the resultant determinant

A.3 Intermediate steps in the proof of Theorem 2.4.4

Let us rewrite as S = (α

Then, the matrix M reads:

In order to prove the equality in Eq. (2.22), we re-arrange the terms in the matrix A by using some matrix inverse identities. The first identity we apply is

where X is a invertible square matrix. By plugging X = α -1 I and Y = B, we write A as follows:

Noticing I -L(α -1 I + L) -1 = (I + αL) -1 even more simplifies:

Let us go back to the original limit. Another matrix trick allows to write:

Plugging A as found above gives:

A.3 Intermediate steps in the proof of Theorem 2.4.4

167

The final manipulation we do is due to the identity (I + XY) -1 = I -X(I + YX) -1 Y.

We re-write the final matrix inverse by setting X = Q 1/2 and Y = (α -1 I + L) -1 Q 1/2 . By doing so, we finish this proof as follows:

A.4 Proof of Lemma A.4.1

Lemma A.4.1. Given the extended graph G = (V , E , w), for all edge subsets S ⊆ E with |S| = |V | -2 = n -1, the following inequality holds:

where i < j assumed without loss of generality. Then, with a certain column permutation, one can find:

where P i is the associated permutation matrix. A similar can be done for B {S|-j,-Γ} as perm(B {S|-j,-Γ} ) = P j B {S|-j,-Γ} = [b i X]. Note that the determinants of these permutation matrices read det P i = (-1) j-1 and det P j = (-1) i for i < j. Then, one can rewrite the initial determinant product as:

≤ 0 implies that the sign of the overall product is (-1) i+j and finishes the proof:

X T X is a semi-positive definite matrix, thus det X T X ≥ 0. By definition, the

Due to this, b T j b i and all entries of both b T j X and X T b i are non-positive. Moreover, X T X writes a M-matrix whose inverse is entry-wise non-negative. Summing these up gives b The condition on the numerator is satisfied only if Φ q has exactly two roots and the pairs (a, b) and (x, y) must be in different trees. Then for an arbitrary forest φ with exactly two trees rooted in v 1 and v 2 , this sum boils down to: Assuming the condition is verified by all edge pairs 1 , one reaches to the maximum value of this sum as follows:

(x,y)∈E (i,j)∈E

Applying this upper bound on the expectation term in Eq. (A.4) concludes this section.