<figure><img src="image_1.png" /></figure>

HAL
open science

# Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs *

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelme, Nicolas Tremblay

► To cite this version:

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelme, Nicolas Tremblay. Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs *. 2024. hal-04524778

HAL Id: hal-04524778

https://hal.science/hal-04524778

Preprint submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelme, Nicolas Tremblay

► To cite this version:

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelme, Nicolas Tremblay. Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs. SIAM Journal on Mathematics of Data Science, 2025. hal-04524778v2

HAL Id: hal-04524778

https://hal.science/hal-04524778v2

Submitted on 14 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs*

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelmé, Nicolas Tremblay
CNRS, Univ. Grenoble Alpes, Grenoble INP, GIPSA-lab, France

**Abstract.** Random diffusions are a popular tool in Monte-Carlo estimations, with well established algorithms such as Walk-on-Spheres (WoS) going back several decades. In this work, we introduce diffusion estimators for the problems of angular synchronization and smoothing on graphs, in the presence of a rotation associated to each edge. Unlike classical WoS algorithms, these estimators allow for *global* estimations by propagating along the branches of *multi-type spanning forests*, and we show that they can outperform standard numerical-linear-algebra solvers in challenging instances, depending on the topology and density of the graph.

**Key words.** RandNLA, graphs, random spanning forests, connection, graph signal processing, synchronization

**AMS subject classifications.** 68W20, 68Q87, 60G35, 62M99

**Introduction.** Data processing over graphs is of interest in many fields of research, such as discrete geometry, signal processing or graph machine learning for instance. A common setting is that of nodes supporting some kind of data, typically numerical values, whereas the edges describe the topology of the underlying space. In a number of situations, the edges can carry more precise geometric information regarding *how* two data points defined on adjacent nodes should be compared, and thus processed. In geometrical parlance, this information is known as a *connection* (an idea first formalized in [13]).

We focus here on a connection describing rotations of angle $\theta_{i,j}$ associated to each edge $(i, j)$, a setting which has found a growing number of applications over the last decade, including *angular synchronization* [63, 71], signal processing over directed graphs [24, 73], discrete geometry processing [60], or even direction of arrival estimation [48, 54]. Solving these problems requires computing the solution to large numerical-linear-algebra (NLA) problems, whose formulation relies on the *connection Laplacian* $L_{\theta}$ (the formal definition of $L_{\theta}$ is postponed to Section 1.2), a linear operator encoding both the topology of the graph *and* the rotations $\theta_{i,j}$.

Computing an exact solution to these problems (*e.g.*, the smoothing problem we will consider) requires solving a linear system, with generic time complexity in $O(n^3)$) remains prohibitive beyond moderately large graphs ($n \approx 10^4$), so that approximations computed from iterative Krylov-subspaces-based algorithms are typically used instead [56, 57]. The main drawback of these methods is that their convergence can depend in a complex manner on the spectrum of the matrix at hand [37] (linear-system solvers for instance require good-quality preconditioners to ensure fast convergence).

Randomized numerical linear algebra (RandNLA) is a successful and modern alternative [21, 45]: Monte-Carlo estimators allow flexible schemes of computation (*e.g.*, parallelized or distributed) and can exhibit both advantageous complexities and state-of-the-art practical performances, in spite of their slow convergence rates in $O(\sigma^2/m)$ (with $\sigma^2$ the variance of the estimator and $m$ the number of Monte-Carlo samples). Methods specialized in graphs, based

<footnote>* This work is an extension of the conference proceedings [31].</footnote>
# Random Multi-Type Spanning Forests for Synchronization on Sparse Graphs*

Hugo Jaquard, Pierre-Olivier Amblard, Simon Barthelmé, Nicolas Tremblay
CNRS, Univ. Grenoble Alpes, Grenoble INP, GIPSA-lab, France

**Abstract.** Random diffusions are a popular tool in Monte-Carlo estimations, with well established algorithms such as Walk-on-Spheres (WoS) going back several decades. In this work, we introduce diffusion estimators for the problems of angular synchronization and smoothing on graphs, in the presence of a rotation associated to each edge. Unlike classical WoS algorithms that are point-wise estimators, our diffusion estimators allow for *global* estimations by propagating along the branches of random spanning subgraphs called *multi-type spanning forests*. Building upon efficient samplers based on variants of Wilson's algorithm, we show that our estimators outperform standard numerical-linear-algebra solvers in challenging instances, depending on the topology and density of the graph.

**Key words.** RandNLA, graphs, random spanning forests, connection, graph signal processing, synchronization

**AMS subject classifications.** 68W20, 68Q87, 60G35, 62M99

**Introduction.** Data processing over graphs is of interest in many fields of research, such as discrete geometry, signal processing or graph machine learning. A common setting is that of nodes supporting some kind of data, typically numerical values, whereas the edges describe the topology of the underlying space. In a number of situations, the edges can carry more precise geometric information regarding *how* two data points defined on adjacent nodes should be compared, and thus processed. In geometrical parlance, this information is known as a *connection* (an idea first formalized in [13]).

We focus here on a connection describing rotations of angle $\theta_{i,j}$ associated to each edge $(i, j)$, a setting which has found a growing number of applications over the last decade, including *angular synchronization* [75, 84], signal processing over directed graphs [26, 87], discrete geometry processing [71], or even direction of arrival estimation [55, 65]. Solving these problems requires computing the solution of large numerical-linear-algebra (NLA) problems, whose formulation relies on the *connection Laplacian* $L_{\theta}$ (the formal definition of $L_{\theta}$ is postponed to Section 1.2), a linear operator encoding both the topology of the graph and the rotations $\theta_{i,j}$.

Computing an exact solution to many of these problems is prohibitive beyond moderately large graphs ($n \approx 10^4$)<fnref n="1" />, so that approximations via iterative Krylov-subspace-based algorithms (such as the popular Conjugate Gradient algorithm) are typically used instead [67, 68]. The main drawback of these methods is their spectrum-dependent convergence rate [43], that often requires hard-to-find, good-quality preconditioners to ensure fast convergence.

Randomized numerical linear algebra (RandNLA) is a successful and modern alternative [21, 52]: Monte-Carlo estimators allow flexible schemes of computation (e.g., parallelized or distributed) and can exhibit both advantageous complexities and state-of-the-art practical performances, in spite of their slow convergence rates in $O(\sigma/\sqrt{m})$ (with $\sigma^2$ the variance of the estimator and $m$ the number of Monte-Carlo samples). Methods specialized in graphs, based

<footnote>*This work is an extension of the conference proceedings [36].</footnote>

<footnote n="1">¹For instance, Tikhonov smoothing requires solving a linear system, with generic time complexity in $O(n^3)$.</footnote>
on random walks and decompositions of the graph, have also appeared (e.g. [39, 51, 65]).

In this paper, we leverage novel connection-aware random decompositions of the graph, and propose RandNLA estimators for two connection-Laplacian-based problems:

* A graph Tikhonov smoothing problem in the presence of a connection.

* The angular synchronization problem on graphs.

**Tikhonov smoothing.** Connection-aware graph Tikhonov smoothing amounts to solving:

$$ (0.1) \qquad \underset{f \in \mathbf{C}^n}{\operatorname{argmin}} \quad q \|f - g\|_2^2 + \frac{1}{2} \sum_{i,j} |f(j) - e^{i\theta_{i,j}} f(i)|^2, $$

where $g \in \mathbf{C}^n$, $q \in \mathbf{R}_+^*$ is an a priori known regularization parameter, and the sum is over a subset of ordered pairs of indices $(i, j)$ describing the edges of the graph. This problem appears in different contexts, such as vector field extension in discrete geometry processing [60], or directed graph signal processing [24]. Here, the regularization term penalizes functions that are not locally coherent with respect to the connection, and can be expressed using the connection Laplacian $L_\theta$. The solution to this first problem takes the form $q(L_\theta + qI)^{-1}g$, and can be leveraged in the second problem we consider: angular synchronization.

**Angular Synchronization.** The objective is to recover a set of $n$ unknown angles $\omega = (\omega)_i \subseteq [0, 2\pi)^n$ from measured pairwise offset measurements $\{\theta_{i,j}\}_{i,j} [63]$:

$$ (0.2) \qquad \theta_{i,j} = \omega_j - \omega_i + \varepsilon_{i,j} \quad \operatorname{mod} 2\pi, $$

where $\varepsilon_{i,j}$ represents some unknown degradation of the measurement. This task appears in many structured signal processing problems, where it is often a key component in state-of-the-art recovery methods, such as perceived luminance reconstruction [70], ptychography [23], ranking [16], clock synchronization [25] or phase reconstruction [1], and also appears in the statistical physics literature (e.g. [14]).

In practice, we may only observe a subset of all such measurements $\theta_{i,j}$, and the problem is naturally formulated on a graph $\mathcal{G}$ with $n$ nodes whose set of edges $\mathcal{E}$ is indexed by the number of measurements, and where edge $(i, j)$ (resp. $(j, i)$) carries the offset $\theta_{i,j}$ (resp. $(\theta_{j,i} = -\theta_{i,j})$). If noiseless measurements $\theta_{i,j} = \omega_j - \omega_i$ are available, exact recovery can be trivially performed up to a global phase shift, by propagating values according to the offset measurements along a spanning tree of $\mathcal{G}$ (a procedure explained in Figures 0.1a and 0.1b).

The problem becomes more involved when considering imperfect measurements $\theta_{i,j}$, with non-zero noise $\varepsilon_{i,j}$. In general, exact angular synchronization can no longer be performed in this noisy regime, as long as even one incoherent cycle is present in the graph (see Figure 0.1c).

A common workaround, first introduced in [70], consists in quantifying the incoherence of an angular assignment $s \in [0, 2\pi)^n$ as

$$ (0.3) \qquad \mathcal{I}(s) = \sum_{\{i,j\} \in \mathcal{E}} (2 - 2 \cos((s_j - s_i) - \theta_{i,j})), $$

before minimizing this incoherence over all possible assignments. While this problem is non-convex and NP-hard in general [72] (see also [10] for a discussion), different techniques allowing to recover a solution have been proposed [10, 27, 50, 63, 71], be it approximately via relaxations or for specific noise regimes or topologies. Most of those techniques also rely on the introduction of the connection Laplacian $L_\theta$ (see Section 1.2).
on random walks and decompositions of the graph, have also appeared (e.g. [45, 61, 77]).

In this paper, we leverage novel connection-aware random decompositions of the graph, and propose RandNLA estimators for two connection-Laplacian-based problems:

* A graph Tikhonov smoothing problem in the presence of a connection.

* The angular synchronization problem on graphs.

**Tikhonov smoothing.** Connection-aware graph Tikhonov smoothing amounts to solving:

$$ (0.1) \qquad f_* = \underset{f \in \mathbf{C}^n}{\operatorname{argmin}} q \|f - g\|_2^2 + \frac{1}{2} \sum_{i,j} |f(j) - e^{i\theta_{i,j}} f(i)|^2, $$

where $g \in \mathbf{C}^n$ is a complex graph signal one wishes to smooth (e.g., denoise), $q \in \mathbf{R}_+^*$ is an a priori known regularization parameter, and the sum is over a subset of ordered pairs of indices $(i, j)$ describing the edges of the graph. This problem appears in different contexts, such as vector field extension in discrete geometry processing [71], or directed graph signal processing [26]. Here, the regularization term (the sum) penalizes functions that are not locally coherent with respect to the connection. The solution to Eq. (0.1) can be expressed using the connection Laplacian $L_\theta$: $f_* = q(L_\theta + qI)^{-1}g$. This inverse operation has a cubic complexity in $n$ and is prohibitive for large graphs, requiring the use of efficient estimators. Connection-aware Tikhonov smoothing also appears as an intermediate step in angular synchronization (the second problem we consider) as well as in iterative solvers for yet other problems$^2$.

**Angular Synchronization.** The objective is to recover a set of $n$ unknown angles $\omega = (\omega)_i \subseteq [0, 2\pi)^n$ from measured pairwise offset measurements $\{\theta_{i,j}\}_{i,j} [75]$:

$$ (0.2) \qquad \theta_{i,j} = \omega_j - \omega_i + \varepsilon_{i,j} \mod 2\pi, $$

where $\varepsilon_{i,j}$ represents some unknown degradation of the measurement. This task appears in many structured signal processing problems, where it is often a key component in state-of-the-art recovery methods, such as perceived luminance reconstruction [83], ptychography [24], ranking [16], clock synchronization [28] or phase reconstruction [1], and also appears in the statistical physics literature (e.g. [14]).

In practice, we may only observe a subset of all such measurements $\theta_{i,j}$, and the problem is naturally formulated on a graph $\mathcal{G}$ with $n$ nodes whose set of edges $\mathcal{E}$ is indexed by the number of measurements, and where edge $(i, j)$ (resp. $(j, i)$) carries the offset $\theta_{i,j}$ (resp. $(\theta_{j,i} = -\theta_{i,j})$). If noiseless measurements $\theta_{i,j} = \omega_j - \omega_i$ are available, exact recovery can be trivially performed up to a global phase shift, by propagating values according to the offset measurements along a spanning tree of $\mathcal{G}$ (a procedure explained in Figures 0.1a and 0.1b).

The problem becomes more involved when considering imperfect measurements $\theta_{i,j}$, with non-zero noise $\varepsilon_{i,j}$. In general, exact angular synchronization can no longer be performed in this noisy regime, as long as even one incoherent cycle is present in the graph (see Figure 0.1c).

<footnote>$^2$This is for instance the case for the so-called “edge-lasso” over graphs [72] (which only appears in the literature on graphs without a connection, but is relevant on graphs with a connection as well), which corresponds to an $l_1$-regularization problem over graph signals. The solution to this problem can for instance be obtained from iteratively-reweighted least-squares (see, e.g., [61]), or from an ADMM strategy (e.g., [27]), or other proximal-operator-based algorithms [59]; all of these strategies may involve solving graph Tikhonov regularization problems.</footnote>
<figure><img src="image_3.png" /><figcaption>(a) Propagation from the large blue node (bottom left) to the large blue square (top right) along the path in blue. A rotation is applied when traversing each edge.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Synchronization by propagation from the large blue node to all other nodes in the graph, along the branches of the blue spanning tree (only possible if there is no incoherence).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Propagation from the large blue node to the large blue rectangle (bottom left), along two different paths (blue and purple). Here, the rotations are incoherent and the result of the propagation depends on the specific path.</figcaption></figure>

**Figure 0.1:** Propagations on the 4 × 2 grid graph, along different paths. The rotation angles $\theta_{i,j}$'s associated to each edge are represented as circular rotating arrows (only drawn along the paths we consider). Propagations always start from the large blue node, and propagated angles are represented as straight arrows (top left of each node). Left 0.1a: propagation along an arbitrary path. Center 0.1b: exact synchronization achieved by propagation along a spanning tree, in the absence of noise. Right 0.1c: synchronization impossible due to noise and incoherence along cycles.

**Our contributions.** We propose novel RandNLA estimators for the smoothing and angular synchronization problems. We rely on propagations along branches of *Multi-Type Spanning Forests* (MTSF): spanning subsets of both edges and nodes of a graph, whose connected components are either rooted trees or unicycles (connected subsets of edges containing *exactly one* cycle). Our main theorem concerns our estimator for the graph Tikhonov smoothing problem, and can be roughly described as follows (see Theorem 2.4 for a formal statement).

**Theorem 0.1.** For MTSFs sampled from the probability distribution in Equation (2.4), propagating the values from the roots of the trees to the other nodes yields an unbiased estimation of the solution of the connection-aware Tikhonov smoothing problem of Equation (0.1).

See Figure 0.2 for an illustration. We include a teasier runtime comparison with a standard deterministic solver in Figure 0.3, which shows that our estimators are *less sensitive* to the density of the graph, and can provide significant speed-ups for equivalent precision. In practice, we obtain a *fast and scalable* algorithm. MTSF-sampling is achieved via a Wilson-like random-walk-based algorithm with runtime *linear* in the number of edges of the graph. From the theoretical side, our arguments rely on the theory of *Determinantal Point Processes* (DPPs) [29, 44], and generalize the combinatorial analyses of [35, 51].

A preliminary version of this work already appeared in [31], where we presented our smoothing estimator along with two variance-reduction techniques, and an application to a ranking problem for a simple synthetic data model. We improve and extend this work in a
<figure><img src="image_2.png" /><figcaption>(a) Propagation from the large blue node (bottom right) to the large blue square (top left) along the path in blue. A rotation is applied when traversing each edge.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Synchronization by propagation from the large blue node to all other nodes in the graph, along the branches of the blue spanning tree (only possible if there is no incoherence).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Propagation from the large blue node to the large blue rectangle (bottom left), along two different paths (blue and purple). Here, the rotations are incoherent and the result of the propagation depends on the specific path.</figcaption></figure>

Figure 0.1: Propagations on the 4 × 2 grid graph, along different paths. The rotation angles $\theta_{i,j}$'s associated to each edge are represented as circular rotating arrows (only drawn along the paths we consider). Propagations always start from the large blue node, and propagated angles are represented as straight arrows (top left of each node). Left 0.1a: propagation along an arbitrary path. Center 0.1b: exact synchronization achieved by propagation along a spanning tree, in the absence of noise. Right 0.1c: synchronization impossible due to noise and incoherence along cycles.

A common workaround, first introduced in [83], consists in quantifying the incoherence of an angular assignment $s \in [0, 2\pi)^n$ as

$$ (0.3) \qquad \mathcal{I}(s) = \sum_{\{i,j\} \in \mathcal{E}} (2 - 2 \cos((s_j - s_i) - \theta_{i,j})) $$

before minimizing this incoherence over all possible assignments. While this problem is non-convex and NP-hard in general [10, 86], different techniques allowing to recover a solution have been proposed [10, 32, 60, 75, 84], be it approximately via relaxations or for specific noise regimes or topologies. Most of these methods rely on the connection Laplacian $L_{\theta}$.

**Our contributions.** We propose novel RandNLA estimators for the smoothing and angular synchronization problems. We rely on propagations along branches of *Multi-Type Spanning Forests* (MTSF): spanning subsets of both edges and nodes of a graph, whose connected components are either rooted trees or unicycles (connected subsets of edges containing *exactly one cycle*). Our main theorem concerns our estimator for the graph Tikhonov smoothing problem, and can be roughly described as follows (see Theorem 2.4 for a formal statement).

**Theorem 0.1.** For MTSFs sampled from the correct distribution (the one of Eq. (2.4)), propagating within each tree the value of its root to all of its other nodes yields an unbiased estimation of the solution of the connection-aware Tikhonov smoothing problem of Eq. (0.1).

See Figure 0.2 for an illustration. We include a teaser runtime comparison with a standard deterministic solver in Figure 0.3, which shows that our estimators are less sensitive to the density of the graph, and can provide significant speed-ups for equivalent precision.
<figure><img src="image_6.png" /><figcaption>(a) A MTSF in purple, with two unicycles (top left, bottom) and one tree (right). The large purple node is the root of the tree.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Estimation: we propagate the values from the root (blue circle) to its tree via the blue arrows. The estimation is 0 on the unicycles (represented by arrowless dots at the top left of the nodes).</figcaption></figure>

Figure 0.2: Theorem 0.1 illustrated on the 4 × 4 grid graph. Left: a MTSF. Right: estimation by propagations along the branches of the MTSF. The estimation in $C$ is represented using angled arrows (top left of the nodes).

number of different manners, some of which we list below.

**Theoretical results.** In addition to Theorem 2.4, we derive a connection-aware Feynman-Kac formula (Proposition 2.2), resulting in a local (node-wise) random-walk based estimator for the smoothing problem, similar to walk-on-spheres algorithms [49, 58]. In comparison, our MTSF-based estimators allow for *global* estimation on all the nodes at once. We further relate these two estimators, which may pave the way to generalizations (see Section 13 of the Supplementary Material). We also improve one of our variance-reduction techniques to better handle heterogeneous degree distributions, with significant improvements.

**Methodological and experimental contributions.** We compare our estimators with standard Krylov-subspaces methods on synthetic data, for both the angular synchronization and smoothing problems. For the angular synchronization problem, we leverage our smoothing estimator as an iterative step in existing approaches [10, 63, 71]. We analyze and contrast the behavior of these methods on different graph topologies which, up to our knowledge, has never been studied in the literature before. Our results show that MTSF-based estimations can outperform standard deterministic methods whenever the graph is not *very sparse* (the gain getting bigger as the density increases). See Figure 0.3. The code used for these experiments is publicly available<fnref n="1" />.

**Related Work.** Our main result can be understood as a specialization of the approach in [19,20], based on DPPs, a class of distributions exhibiting negatively-correlated samples [29, 38], and resulting in unbiased estimators for least-squares problems. In contrast to this general work, the specific structure of MTSFs allows efficient sampling and practical algorithms.

<footnote n="1">¹https://gricad gitlab.univ-grenoble-alpes.fr / gaia / synchromtsf</footnote>
<figure><img src="image_5.png" /><figcaption>(a) A MTSF in purple, with two unicycles (top left, bottom) and one tree (right). The large purple node is the root of the tree.

(b) Estimation. Tree: the value of the root (blue circle) is propagated to the rest of the tree via the blue arrows. Unicycles: the estimation is 0 (arrowless dots at the top left of the nodes).</figcaption></figure>

Figure 0.2: Theorem 0.1 illustrated on the 4 × 4 grid graph. Left: a MTSF. Right: estimation by propagations along the branches of the trees of the MTSF. The estimation in **C** is represented using angled arrows (top left of the nodes).

In practice, we obtain a *fast and scalable* algorithm. MTSF-sampling is achieved via a Wilson-like random-walk-based algorithm with runtime *linear* in the number of edges of the graph. On the theoretical side, our arguments rely on the theory of *Determinantal Point Processes* (DPPs) [34, 51], and generalize the combinatorial analyses of [41, 61].

A preliminary version of this work already appeared in [36], where we presented our smoothing estimator along with two variance-reduction techniques, and an application to a ranking problem for a simple synthetic data model. We improve and extend this work in a number of different manners, some of which we list below.

**Theoretical results.** In addition to Theorem 2.4, we derive a connection-aware Feynman-Kac formula (Proposition 2.2), resulting in a local (node-wise) random-walk based estimator for the smoothing problem, similar to walk-on-spheres algorithms [56, 69]. In comparison, our MTSF-based estimators allow for *global* estimation on all the nodes at once. We further relate these two estimators, which may pave the way to generalizations (see Section 16 of the Supplementary Material). We also improve one of our variance-reduction techniques to better handle heterogeneous degree distributions, with significant improvements.

**Methodological and experimental contributions.** We compare our estimators with standard Krylov-subspaces methods on synthetic data, for both the angular synchronization and smoothing problems. For the angular synchronization problem, we leverage our smoothing estimator as an iterative step in existing approaches [10, 75, 84]. We analyze and contrast the behavior of these methods on different graph topologies which, up to our knowledge, has never been studied in the literature before. Our results show that MTSF-based estimations can outperform standard deterministic methods whenever the graph is not *very sparse* (the gain getting bigger as the density increases). See Figure 0.3. The code used for these
<figure><img src="image_8.png" /><figcaption>Figure 0.3: Runtime-precision comparisons for the graph Tikhonov smoothing problem, comparing two MTSF-based estimators (MTSF and MTSF+GS) against a diagonally-preconditioned conjugate-gradient-descent ($D^{-1}$+CG) on degree-corrected stochastic block model graphs (DC-SBM, see Section 3 for a definition) with 10000 nodes. Results are averaged over 10 graphs generated from DC-SBM models with two different parametrizations, resulting in different average degrees $\bar{d}$: on the left, a DC-SBM with average degree $\bar{d} \approx 30$ and, on the right, a DC-SBM with average degree $\bar{d} \approx 420$. x-axis: average runtime. y-axis: average reconstruction error. Each data-point corresponds to a number $m \in \{1, 2, 3, 5, 8, 13, 22, 36, 60, 100\}$ of MTSFs (or of CG iterations) used in the estimation. The vertical (resp. horizontal) line is the computation-time (resp. error) of an exact Cholesky solver (see Section 3 for details).</figcaption></figure>

There exists similar applications of random spanning forests distributions in RandNLA for problems on connection-free graphs [3,4], such as the trace estimation of (regularized) inverse Laplacians [6] or, the Tikhonov regularization and interpolation of graphs signals [51]. Multi-Type Spanning Forests have also been used to build spectral sparsifiers for the (regularized) connection Laplacian in [22], resulting in randomized preconditioners. Our Feynman-Kac formula is inspired from a similar result in [34] for continuous-time random walks, but is specially tailored to the graph Tikhonov smoothing problem; discrete-time walks also exhibit intriguing links with propagations on MTSFs (see Section 13 of the Supplementary Material).

**Organisation of the Paper.** In Section 1, we introduce preliminary background on graphs and connections used in the rest of the paper. In Section 2, we present our Feynman-Kac formula (Section 2.1) and our MTSF-based estimators (Section 2.2), analyze an efficient sampling algorithm for MTSFs (Section 2.3), and describe our variance-reduction techniques (Section 2.4). We analyze the numerical behavior of our estimators in Section 3, and compare them with standard deterministic algorithms (Sections 3.1 and 3.2). In Section 4, we propose an iterative scheme leveraging our smoothing estimators to solve the angular synchronization problem (Section 4.1), illustrate its application to a denoising problem inspired from *cryogenic electron microscopy* (Section 4.2), and compare this scheme with standard deterministic methods (Section 4.3). Proofs are deferred to the Supplementary Material.
<figure><img src="image_6.png" /><figcaption>Figure 0.3: Runtime-precision comparisons for the graph Tikhonov smoothing problem, comparing two MTSF-based estimators (MTSF and MTSF+GS) against a diagonally-preconditioned conjugate-gradient-descent ($D^{-1}$+CG) on degree-corrected stochastic block model graphs (DC-SBM, see Section 3 for a definition) with 10000 nodes. Results are averaged over 10 graphs generated from DC-SBM models with two different parametrizations, resulting in different average degrees $\bar{d}$: on the left, a DC-SBM with average degree $\bar{d} \approx 30$ and, on the right, a DC-SBM with average degree $\bar{d} \approx 420$. x-axis: average runtime. y-axis: average reconstruction error. Each data-point corresponds to a number $m \in \{1, 2, 3, 5, 8, 13, 22, 36, 60, 100\}$ of MTSFs (or of CG iterations) used in the estimation. The vertical (resp. horizontal) line is the computation-time (resp. error) of an exact Cholesky solver (see Section 3 for details).</figcaption></figure>

experiments is publicly available<fnref n="3" />.

**Related Work.** Our main result can be understood as a specialization of the approach in [19,20], based on DPPs, a class of distributions exhibiting negatively-correlated samples [34, 44], and resulting in unbiased estimators for least-squares problems. {However, and in contrast to generic DPPs that are in general expensive to sample from, the specific structure of MTSFs allows for very efficient sampling algorithms and, in turn, practical graph RandNLA algorithms that can compete with deterministic state-of-the-art solvers.

There exists similar applications of random spanning forests distributions in RandNLA for problems on connection-free graphs [3, 4], such as the trace estimation of (regularized) inverse Laplacians [6] or, the Tikhonov regularization and interpolation of graphs signals [61]. Multi-Type Spanning Forests have also been used to build spectral sparsifiers for the (regularized) connection Laplacian in [22], resulting in randomized preconditioners. Our Feynman-Kac formula is inspired from a similar result in [40] for continuous-time random walks, but is specially tailored to the graph Tikhonov smoothing problem; discrete-time walks also exhibit intriguing links with propagations on MTSFs (see Section 16 of the Supplementary Material).

**Organisation of the Paper.** In Section 1, we introduce preliminary background on graphs and connections used in the rest of the paper. In Section 2, we present our Feynman-Kac formula (Section 2.1) and our MTSF-based estimators (Section 2.2), analyze an efficient sampling algorithm for MTSFs (Section 2.3), and describe our variance-reduction techniques (Sec-

<footnote n="3">³https://gricad gitlab.univ-grenoble-alpes.fr / gaia / synchromtsf</footnote>
**1. Background.** A graph $\mathcal{G}$ is defined as a set on nodes $\mathcal{V}$ interconnected by a set of edges $\mathcal{E}$. The edges of $\mathcal{E}$ are non-oriented edges. However, for the purpose of this paper, we will also need to consider their oriented counterparts, denoted by $\vec{\mathcal{E}}$. The size of $\vec{\mathcal{E}}$ is twice the size of $\mathcal{E}$: each edge $e \in \mathcal{E}$ is associated to two oriented edges in $\vec{\mathcal{E}}$. Choosing arbitrarily an orientation for $e$ and writing $s_e$ and $t_e$ its source and target, we have $e = (s_e, t_e) \in \vec{\mathcal{E}}$, and its reversely-oriented edge $e^* = (t_e, s_e) \in \vec{\mathcal{E}}$. We state our results for unweighted graphs, the weighted case is included in the Supplementary Material.

**1.1. Combinatorial Laplacian.** An elementary instance of graph-supported data is that of real values attached to the nodes of the graph, usually formalized as a vector $f \in \mathbf{R}^{\mathcal{V}}$, and called a *graph signal* [62]. The regularity of such a signal can be quantified using the *combinatorial Laplacian L*: $\mathbf{R}^{\mathcal{V}} \to \mathbf{R}^{\mathcal{V}}$, a symmetric semi-definite positive operator, conveniently expressed as $L = D - A$, where $D$ is the diagonal degree matrix ($D_{i,i} = d_i$ the degree of node $i$) and $A$ the adjacency matrix of the graph (with $A_{(i,j)} = 1$ if $(i,j) \in \vec{\mathcal{E}}$ and $0$ otherwise) [15]. The quadratic form associated to $L$ acts as:

$$ (1.1) \qquad \langle f, Lf \rangle = \frac{1}{2} \sum_{e \in \vec{\mathcal{E}}} |f(t_e) - f(s_e)|^2, $$

which associates a high value to functions with important local variations over the edges of the graph, and serves as a basis for the notion of *frequency* for graph signals [62, 67]. These operators are also related to combinatorial properties of the graph $\mathcal{G}$ [8], for instance:

(P1) ker $L$ is always non-empty, with dimension the number of connected components of $\mathcal{G}$. This observation forms one of the bases for spectral clustering on graphs [68].

(P2) $D^{-1}A$ is the matrix of the natural *random walk* on $\mathcal{G}$, which transitions from a node to one of its neighbors with uniform probability at each step, with $(D^{-1}A)_{i,j}^l$ the probability that a path of length $l$ starting at $i$ ends up in node $j$.

**1.2. Connection Laplacian: Definition and Basic Properties.** *Unitary connections* are a practical way to add additional information to the graph’s structure<fnref n="2" />. We describe in the following how we can capture additional rotations associated to edges in the form of a Laplacian-like operator. The main idea is to represent rotations as *unitary complex numbers*, which will naturally lead to consider *complex-valued* functions $f \in \mathbf{C}^{\mathcal{V}}$ when studying variational properties on $\mathcal{G}$. The entries of $f$ should be understood as belonging to different copies $\mathbf{C}_v$ of the complex plane $\mathbf{C}$ associated to each node $v \in \mathcal{V}$, see Figure 1.1³.

We will associate to each directed edge in $\vec{\mathcal{E}}$ a map representing the transformation along that edge, which is known as a *connection* [35]. A *unitary connection* $\Psi$ on a graph $\mathcal{G}$ is a collection of unitary linear maps $(\psi_e)_{e \in \vec{\mathcal{E}}}$, with each map $\psi_e: \mathbf{C}_{s_e} \to \mathbf{C}_{t_e}$ acting as multiplication by a unitary complex number $e^{i\theta_e}$, so that $\psi_e(z) = e^{i\theta_e} \cdot z$, with $\iota$ the complex imaginary unit (we sometimes abuse notation and write $\psi_e = e^{i\theta_e}$ when there is no risk of confusion). In addition, we require that $\psi_{e^*} = \psi_e^*$ (i.e. $\psi_{(t_e, s_e)} = \psi_{(s_e, t_e)}^*$), so that the

<footnote n="2">²By describing explicit geometrical transformations between the signal values of the nodes of $\mathcal{G}$.</footnote>

<footnote n="3">³This association of a copy of $\mathbf{C}$ is analogous to a fiber bundle over a manifold (e.g. its tangent bundle), and is known as a *complex line bundle* over $\mathcal{G}$.</footnote>
tion 2.4). We analyze the numerical behavior of our estimators in Section 3, and compare them with standard deterministic algorithms (Sections 3.1 and 3.2). In Section 4, we propose an iterative scheme leveraging our smoothing estimators to solve the angular synchronization problem (Section 4.1), illustrate its application to a denoising problem inspired from cryogenic electron microscopy (Section 4.2), and compare this scheme with standard deterministic methods (Section 4.3). Proofs are deferred to the Supplementary Material.

**1. Background.** A graph $\mathcal{G}$ is defined as a set on nodes $\mathcal{V}$ interconnected by a set of edges $\mathcal{E}$. The edges of $\mathcal{E}$ are non-oriented edges. However, for the purpose of this paper, we will also need to consider their oriented counterparts, denoted by $\vec{\mathcal{E}}$. The size of $\vec{\mathcal{E}}$ is twice the size of $\mathcal{E}$: each edge $e \in \mathcal{E}$ is associated to two oriented edges in $\vec{\mathcal{E}}$. Choosing arbitrarily an orientation for $e$ and writing $s_e$ and $t_e$ its source and target, we have $e = (s_e, t_e) \in \vec{\mathcal{E}}$, and its reversely-oriented edge $e^* = (t_e, s_e) \in \vec{\mathcal{E}}$. We state our results for unweighted graphs, but all our results generalize to the weighted case, as detailed in the Supplementary Material.

**1.1. Combinatorial Laplacian.** An elementary instance of graph-supported data is that of real values attached to the nodes of the graph, usually formalized as a vector $f \in \mathbf{R}^V$, and called a *graph signal* [74]. The regularity of such a signal can be quantified using the *combinatorial Laplacian* $\mathbf{L}: \mathbf{R}^V \rightarrow \mathbf{R}^V$, a symmetric semi-definite positive operator, conveniently expressed as $\mathbf{L} = \mathbf{D} - \mathbf{A}$, where $\mathbf{D}$ is the diagonal degree matrix ($\mathbf{D}_{i,i} = d_i$ the degree of node $i$) and $\mathbf{A}$ the adjacency matrix of the graph (with $\mathbf{A}_{(i,j)} = 1$ if $(i,j) \in \vec{\mathcal{E}}$ and 0 otherwise) [15]. The quadratic form associated to $\mathbf{L}$ acts as:

$$ (1.1) \qquad \langle f, \mathbf{L}f \rangle = \frac{1}{2} \sum_{e \in \vec{\mathcal{E}}} |f(t_e) - f(s_e)|^2, $$

which associates a high value to functions with important local variations over the edges of the graph, and serves as a basis for the notion of *frequency* for graph signals [74, 79]. These operators are also related to combinatorial properties of the graph $\mathcal{G}$ [8], for instance:

(P1) $\ker \mathbf{L}$ is always non-empty, with dimension the number of connected components of $\mathcal{G}$.

(P2) $\mathbf{D}^{-1}\mathbf{A}$ is the matrix of the natural *random walk* on $\mathcal{G}$, which transitions from a node to one of its neighbors with uniform probability at each step, with $((\mathbf{D}^{-1}\mathbf{A})^l)_{i,j}$ the probability that a path of length $l$ starting at $i$ ends up in node $j$.

**1.2. Connection Laplacian: Definition and Basic Properties.** *Unitary connections* are a practical way to add additional information to the graph’s structure<fnref n="4" />. We describe in the following how to capture additional rotations associated to edges in the form of a Laplacian-like operator. The main idea is to represent rotations as *unitary complex numbers*, which will naturally lead to consider *complex-valued* functions $f \in \mathbf{C}^V$ when studying variational properties on $\mathcal{G}$. The entries of $f$ should be understood as belonging to different copies $\mathbf{C}_v$ of the complex plane $\mathbf{C}$ associated to each node $v \in \mathcal{V}$, see Figure 1.1⁵.

<footnote n="4">⁴By describing explicit geometrical transformations between the signal values of the nodes of $\mathcal{G}$.</footnote>

<footnote n="5">⁵This association of a copy of $\mathbf{C}$ is analogous to a fiber bundle over a manifold (e.g. its tangent bundle), and is known as a *complex line bundle* over $\mathcal{G}$.</footnote>
transformation associated to an edge traversed in one direction or the other differs only by
conjugation. This notion is illustrated in Figure 1.1.

<figure><img src="image_9.png" /><figcaption>Figure 1.1: A connection $\Psi$ on the triangle graph $K_3$ (in the center). The connections maps $\psi_e$ are associated with angles $\theta_{(v_1,v_2)} = \theta_{(v_2,v_3)} = \theta_{(v_3,v_1)} = \frac{\pi}{4}$, and we represent their action on some $z_1 \in C_{v_1}$. A cyclic path $C = ((v_1, v_2), (v_2, v_3), (v_3, v_1))$ is depicted using directed arrows, we denote by $\psi_C = \psi_{(v_3,v_1)} \circ \psi_{(v_2,v_3)} \circ \psi_{(v_1,v_2)}$ the composition of the connection maps along this cycle (a rotation of $\theta_C = \frac{3\pi}{4}$). Top right: $z_1 \in C_{v_1}$ (bold arrow) and $z'_1 = \psi_C(z_1) \in C_{v_1}$ (blue dotted arrow). Top left, bottom: $z_2 = \psi_{(v_1,v_2)}(z_2) \in C_{v_2}$ and $z_3 = \psi_{(v_2,v_3)}(z_2) \in C_{v_3}$ (dotted arrows).</figcaption></figure>

**Remark 1.1.** In applications, such as angular synchronization, the $\theta_e$'s often represent a priori known rotations, and are frequently modulated by a scale parameter $\gamma \in \mathbf{R}_+$, resulting in connections such that $\psi_e(z) = e^{i\gamma\theta_e}z$. While we will forego this additional parameter in our theoretical results, it is easily recovered by replacing each $\theta_e$ with $\gamma\theta_e$.

*The connection Laplacian.* To a graph endowed with a unitary connection, we associate the connection Laplacian $L_{\theta} : \mathbf{C}^{\mathcal{V}} \rightarrow \mathbf{C}^{\mathcal{V}}$, defined by $L_{\theta} = D - A_{\theta}$, with $A_{\theta}$ a connection- aware adjacency matrix such that $(A_{\theta})_{i,j} = e^{-i\theta_{(i,j)}}$ if $\{i,j\} \in \mathcal{E}$, and $A_{i,j} = 0$ otherwise. This operator generalizes the Laplacian to non-trivial connections, and is associated to the
We will associate to each directed edge in $\vec{\mathcal{E}}$ a map representing the transformation along that edge, which is known as a *connection* [41]. A *unitary connection* $\Psi$ on a graph $\mathcal{G}$ is a collection of unitary linear maps $(\psi_e)_{e \in \vec{\mathcal{E}}}$, with each map $\psi_e : \mathbf{C}_{s_e} \to \mathbf{C}_{t_e}$ acting as multiplication by a unitary complex number $e^{i\theta_e}$, so that $\psi_e(z) = e^{i\theta_e} \cdot z$, with $\iota$ the complex imaginary unit (we sometimes abuse notation and write $\psi_e = e^{i\theta_e}$ when there is no risk of confusion). In addition, we require that $\psi_{e^*} = \psi_e^*$ (i.e. $\psi_{(t_e, s_e)} = \psi_{(s_e, t_e)}^*$), so that the transformation associated to an edge traversed in one direction or the other differs only by conjugation. This notion is illustrated in Figure 1.1.

<figure><img src="image_7.png" /><figcaption>Figure 1.1: A connection $\Psi$ on the triangle graph $K_3$ (in the center). The connection maps $\psi_e$ are associated with angles $\theta_{(v_1,v_2)} = \theta_{(v_2,v_3)} = \theta_{(v_3,v_1)} = \frac{\pi}{4}$, and we represent their action on some $z_1 \in \mathbf{C}_{v_1}$. A cyclic path $C = ((v_1, v_2), (v_2, v_3), (v_3, v_1))$ is depicted using directed arrows, we denote by $\psi_C = \psi_{(v_3,v_1)} \circ \psi_{(v_2,v_3)} \circ \psi_{(v_1,v_2)}$ the composition of the connection maps along this cycle (a rotation of $\theta_C = \frac{3\pi}{4}$). Top right: $z_1 \in \mathbf{C}_{v_1}$ (bold arrow) and $z'_1 = \psi_C(z_1) \in \mathbf{C}_{v_1}$ (blue dotted arrow). Top left, bottom: $z_2 = \psi_{(v_1,v_2)}(z_2) \in \mathbf{C}_{v_2}$ and $z_3 = \psi_{(v_2,v_3)}(z_2) \in \mathbf{C}_{v_3}$.</figcaption></figure>

**Remark 1.1.** In applications, such as angular synchronization, the $\theta_e$'s often represent a priori known rotations, and are frequently modulated by a scale parameter $\gamma \in \mathbf{R}_+$, resulting in connections such that $\psi_e(z) = e^{i\gamma\theta_e}z$. While we will forego this additional parameter in our theoretical results, it is easily recovered by replacing each $\theta_e$ with $\gamma\theta_e$.

**The connection Laplacian.** To a graph endowed with a unitary connection, we associate the *connection Laplacian* $L_\theta : \mathbf{C}^\mathcal{V} \to \mathbf{C}^\mathcal{V}$, defined by $L_\theta = D - A_\theta$, with $A_\theta$ a connection-aware adjacency matrix such that $(A_\theta)_{i,j} = e^{-i\theta_{(i,j)}}$ if $\{i, j\} \in \mathcal{E}$, and $A_{i,j} = 0$ otherwise. This
quadratic form:

$$ (1.2) \qquad \langle f, \mathbf{L}_{\theta} f \rangle = \frac{1}{2} \sum_{e \in \mathcal{E}} |f(t_e) - e^{t\theta_e} f(s_e)|^2. $$

For the trivial connection, with $\psi_e = \mathrm{id}_{C_{s_e}, C_{t_e}}$ (i.e. $\theta_e = 0$) for all edges $e \in \mathcal{E}$, we recover $\mathbf{L}_{\theta} = \mathbf{L}$ and the expression in Equation (1.1). Unlike this specific case, $\langle f, \mathbf{L}_{\theta} f \rangle$ also penalizes functions that are incoherent with respect to the connection, including constant functions.

Let us now answer a natural question: when is the kernel $\mathrm{ker} \mathbf{L}_{\theta}$ non-empty, and what are the functions $f \in \mathbf{C}^{\mathcal{V}}$ that are not penalized by $\langle f, \mathbf{L}_{\theta} f \rangle$?

(P1') Angular synchronization and ker $\mathbf{L}_{\theta}$. When the entries $f_i$ of the vector $f$ are unitary complex numbers $f_i = e^{t s_i}$, we recover from Equation (1.2) the expression in Equation (0.3). In this case, $|f(j) - e^{t\theta_{(i,j)}} f(i)|^2 = 2 - 2\Re((e^{t s_j})^* e^{t\theta_{(i,j)}} e^{t s_i})$, which, thanks to conjugation-invariance of the real part, results in:

$$ (1.3) \qquad \langle f, \mathbf{L}_{\theta} f \rangle = \sum_{\{i,j\} \in E} 2 - 2 \cos((s_j - s_i) - \theta_{(i,j)}). $$

Even though this equality only holds for $f \in U(\mathbf{C})^{\mathcal{V}}$, note that $\langle f, \mathbf{L}_{\theta} f \rangle = 0$ if and only if $f_i = e^{t\theta_{(j,i)}} f_j$ for all $(i, j) \in \mathcal{E}$, that is, if $s_i = s_j + \theta_{(j,i)}$. As a consequence, $\mathbf{L}_{\theta}$ is not only semi-definite positive (as seen from Equation (1.2)), but *positive definite*, unless there is an *exact solution* $x \in U(\mathbf{C})^{\mathcal{V}}$ to the associated angular synchronization problem, in which case the kernel $\mathrm{ker} \mathbf{L}_{\theta}$ is generated by $x$.

This behavior of the smallest eigenpair of $\mathbf{L}_{\theta}$ is in clear contrast to that of $\mathbf{L}$ (which is *always* zero, see property (P1)), and we explain in Section 2 how a generalization of property (P2) pertaining to *random propagations* can be used to solve connection-Laplacian-based problems.

## 2. Random Estimators for Tikhonov Smoothing under a Connection.
The goal of the next Section is to smooth a complex signal on a graph, by way of random propagations. More precisely, we fix a graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ endowed with a connection $\Psi$, and we aim at smoothing a signal $g \in \mathbf{C}^{\mathcal{V}}$ by solving:

$$ (2.1) \qquad \underset{f \in \mathbf{C}^{\mathcal{V}}}{\operatorname{argmin}} \quad q \|f - g\|_2^2 + \langle f, \mathbf{L}_{\theta} f \rangle, $$

with $q \in \mathbf{R}_+^*$. This is nothing but a re-writing of Equation (0.1), and $\langle f, \mathbf{L}_{\theta} f \rangle$ penalizes functions that are not coherent with respect to the connection $\Psi$. Denoting by $f_*$ the optimal solution to this problem, we have.

**Proposition 2.1.** *The solution to Problem (2.1) can be expressed as:*

$$ (2.2) \qquad f_* = q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1}g. $$

The proof of Proposition 2.1, based on straightforward **CR**-calculus, is detailed in Section 8 of the Supplementary Material.

Let us now describe how to recover the solution of this problem via local propagation over random paths on $\mathcal{G}$.
operator generalizes the Laplacian to non-trivial connections. Its quadratic form reads:

$$ (1.2) \qquad \langle f, \mathbf{L}_{\theta} f \rangle = \frac{1}{2} \sum_{e \in \vec{\mathcal{E}}} |f(t_e) - e^{t\theta_e} f(s_e)|^2. $$

For the trivial connection, with $\psi_e = \mathrm{id}_{C_{s_e}, C_{t_e}}$ (i.e. $\theta_e = 0$) for all edges $e \in \vec{\mathcal{E}}$, we recover $\mathbf{L}_{\theta} = \mathbf{L}$ and the expression in Eq. (1.1). Unlike this specific case, $\langle f, \mathbf{L}_{\theta} f \rangle$ also penalizes functions that are incoherent with respect to the connection, including constant functions.

Let us now answer a natural question: when is the kernel $\ker \mathbf{L}_{\theta}$ non-empty, and what are the functions $f \in \mathbf{C}^{\mathcal{V}}$ that are not penalized by $\langle f, \mathbf{L}_{\theta} f \rangle$?

(P1') Angular synchronization and ker $\mathbf{L}_{\theta}$. When the entries $f_i$ of the vector $f$ are unitary complex numbers $f_i = e^{i\theta_i}$, we recover from Eq. (1.2) the expression in Eq. (0.3). In fact, one obtains $|f(j) - e^{i\theta_{(i,j)}} f(i)|^2 = 2 - 2\mathfrak{Re}((e^{i\theta_j})^* e^{i\theta_{(i,j)}} e^{i\theta_i})$, which yields:

$$ (1.3) \qquad \langle f, \mathbf{L}_{\theta} f \rangle = \sum_{\{i,j\} \in \mathcal{E}} 2 - 2 \cos((s_j - s_i) - \theta_{(i,j)}). $$

Even though this equality only holds for $f \in U(\mathbf{C})^{\mathcal{V}}$, note that $\langle f, \mathbf{L}_{\theta} f \rangle = 0$ if and only if $f_i = e^{i\theta_{(j,i)}} f_j$ for all $(i, j) \in \vec{\mathcal{E}}$, that is, if $s_i = s_j + \theta_{(j,i)}$. As a consequence, $\mathbf{L}_{\theta}$ is not only semi-definite positive (as seen from Eq. (1.2)), but positive definite, unless there is an exact solution $x \in U(\mathbf{C})^{\mathcal{V}}$ to the associated angular synchronization problem, in which case the kernel $\ker \mathbf{L}_{\theta}$ is generated by $x$.

This behavior of the smallest eigenpair of $\mathbf{L}_{\theta}$ is in clear contrast to that of $\mathbf{L}$ (which is always zero, see property (P1)), and we explain in Section 2 how a generalization of property (P2) pertaining to random propagations can be used to solve connection-Laplacian-based problems.

## 2. Random Estimators for Tikhonov Smoothing under a Connection.
This section shows how to smooth a complex signal on a graph via random propagations. Let $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ be a graph endowed with a connection $\Psi$. We aim at smoothing a signal $g \in \mathbf{C}^{\mathcal{V}}$ by solving:

$$ (2.1) \qquad f_* = \underset{f \in \mathbf{C}^{\mathcal{V}}}{\operatorname{argmin}} q \|f - g\|_2^2 + \langle f, \mathbf{L}_{\theta} f \rangle, $$

with $q \in \mathbf{R}_+^*$. This is nothing but a re-writing of Eq. (0.1). Note that $\langle f, \mathbf{L}_{\theta} f \rangle$ penalizes functions that are not coherent with respect to the connection $\Psi$.

**Proposition 2.1.** *The solution to Problem (2.1) can be expressed as:*

$$ (2.2) \qquad f_* = q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1}g. $$

The proof, based on standard **CR**-calculus, is in Section 8 of the Supplementary Material. Let us now describe how to estimate $f_*$ via local propagation over random paths on $\mathcal{G}$.

### 2.1. Local Estimation: a Feynman-Kac Formula.
Feynman-Kac formulas express solutions of variational problems as the expectations of stochastic processes: we develop here an instance tailored to the Tikhonov smoothing problem of Eq. (2.1). We start with a few
**2.1. Local Estimation: a Feynman-Kac Formula.** Feynman-Kac formulas express solutions of variational problems as the expectations of stochastic processes: we develop here an instance tailored to the Tikhonov smoothing problem of Equation (2.1). We start with a few definitions. A path $p$ in $\mathcal{G}$ is an ordered sequence of oriented edges in $\mathcal{E}$. We consider the two operations of concatenation $pq$ of two paths $p$ and $q$ and of orientation-reversal $p^*$ of a path $p$, and denote by $P_i^j$ the set of paths from $i$ to $j$ in $\mathcal{G}$. Connection maps extend to paths inductively as $\psi_{pq} = \psi_q \circ \psi_p$, and we will for instance frequently encounter the map $\psi_{p^*} : C_j \to C_i$, acting by composition of the (inverse) rotations encountered along path $p \in P_i^j$ from $i$ to $j$.

The stochastic process we consider is a random walk with a non-zero probability of being interrupted at any node. It is conveniently defined on an extended graph $\mathcal{G}^\Gamma$, with nodes $\mathcal{V} \cup \{\Gamma\}$ and edges $\mathcal{E} \cup \bigcup_{v \in \mathcal{V}} \{\{v, \Gamma\}\}$ (so that all the nodes of $\mathcal{V}$ are connected to $\Gamma$), where $\Gamma$ is an additional boundary node. We then define the random walk $(u_t)_{t \ge 0}$, which starts at some fixed node $u_0 = i \in \mathcal{V}$ and transitions to some other node at time $t$ as follows:

* $u_{t+1} = \Gamma$ with probability $\frac{q}{d_{u_t+q}}$,

* $u_{t+1} = v \in \mathcal{V}$ with probability $\frac{A_{u_t, v}}{d_{u_t+v}}$.

The process ends upon reaching node $\Gamma$. In other words, this is the natural random walk on $\mathcal{G}^\Gamma$, with a stopping criterion at the boundary $\Gamma$. Writing $u_l = \Gamma$ for some time $l$, this process results in a random path $((u_0, u_1), \dots, (u_{l-1}, u_l))$. We denote by $P_i^\Gamma$ the set of all possible paths obtainable via this process, and by $\nu_i$ the resulting probability measure on $P_i^\Gamma$.

We can then show that propagating along paths sampled according to $\nu_i$ converges in expectation to the solution of Problem (2.1) on node $i$. More precisely:

**Proposition 2.2.** *Denoting by $j = u_{l-1}$ the last node reached before $p$ reaches $\Gamma$, we have:*

$$ (2.3) \qquad f_*(i) = \mathbf{E}_{p \sim \nu_i} (\psi_{p^* \Gamma}(g_j)), $$

where $p = p_\Gamma e_\Gamma$ with $e_\Gamma$ the last edge in $p$.

In other words, an unbiased estimate of $f_*(i)$ is obtained by drawing a random path from $i$ to $\Gamma$, and retropropagating the value of $g$ at the node just before $\Gamma$ ($j$ in Proposition 2.2), taking into account the rotations along the path. See Algorithm 2.1. This formula is a discrete-time analog of a Theorem in [34], and we prove a more general version in the Supplementary Material (Section 8, Proposition 8.1), building upon a connection-aware version of property (P2).

<caption>**Algorithm 2.1** Feynman-Kac estimator (Proposition 2.2).</caption>

1: $f_i \leftarrow 0$
2: **Repeat** $m$ times
3: &nbsp;&nbsp;&nbsp;&nbsp;Sample a path $p \in P_i^\Gamma$ from $\nu_i$ &nbsp;&nbsp;&nbsp;&nbsp;▷ By running an interrupted random walk on $\mathcal{G}$
4: &nbsp;&nbsp;&nbsp;&nbsp;$f_i \leftarrow f_i + \psi_{p^* \Gamma}(g_j)$ &nbsp;&nbsp;&nbsp;&nbsp;▷ $j$ the last node before interruption of $p$
5: **Output** $\frac{1}{m} f_i$

Algorithm 2.1 is simple but may be computationally expensive: estimation on one node requires sampling $m$ paths of unbounded length, which needs to be repeated for each of the $|\mathcal{V}|$ nodes in $\mathcal{G}$ (a locality issue inherent to Walk-on-Spheres-type algorithms [49]). We address these limitations in the remainder of this Section and propose novel propagation-based
<caption>**Algorithm 2.1** Feynman-Kac estimator (Proposition 2.2).</caption>

1: $f_i \leftarrow 0$
2: **Repeat** $m$ times
3: Sample a path $p \in P_i^\Gamma$ from $\nu_i$ ▷ By running an interrupted random walk on $\mathcal{G}$
4: $f_i \leftarrow f_i + \psi_{p_\Gamma^*}(g_j)$ ▷ $j$ the last node before interruption of $p$
5: **Output** $\frac{1}{m}f_i$ ▷ An unbiased estimator of $f_*(i)$

definitions. A *path* $p$ in $\mathcal{G}$ is an ordered sequence of oriented edges in $\vec{\mathcal{E}}$. We consider the two operations of concatenation $pq$ of two paths $p$ and $q$ and of orientation-reversal $p^*$ of a path $p$, and denote by $P_i^j$ the set of paths from $i$ to $j$ in $\mathcal{G}$. Connection maps extend to paths inductively as $\psi_{pq} = \psi_q \circ \psi_p$, and we will for instance frequently encounter the map $\psi_{p^*} : C_j \to C_i$, acting by composition of the (inverse) rotations encountered along path $p \in P_i^j$ from $i$ to $j$.

The stochastic process we consider is a random walk with a non-zero probability of being interrupted at any node. It is conveniently defined on an extended graph $\mathcal{G}^\Gamma$, with nodes $\mathcal{V} \cup \{\Gamma\}$ and edges $\mathcal{E} \cup \bigcup_{v \in \mathcal{V}} \{\{v, \Gamma\}\}$ (so that all the nodes of $\mathcal{V}$ are connected to $\Gamma$), where $\Gamma$ is an additional *boundary* node. We then define the random walk $(u_t)_{t \ge 0}$, which starts at some fixed node $u_0 = i \in \mathcal{V}$ and transitions to some other node at time $t$ as follows:

* $u_{t+1} = \Gamma$ with probability $\frac{q}{d_{u_t+q}}$ (recall that $d_{u_t}$ stands for the degree of node $u_t$),

* $u_{t+1} = v \in \mathcal{V}$ with probability $\frac{A_{u_t, v}}{d_{u_t+q}}$.

The process ends upon reaching node $\Gamma$. In other words, this is the natural random walk on $\mathcal{G}^\Gamma$, with a stopping criterion at the boundary $\Gamma$. Writing $u_l = \Gamma$ for some time $l$, this process results in a random path $((u_0, u_1), \dots, (u_{l-1}, u_l))$. We denote by $P_i^\Gamma$ the set of all possible paths obtainable via this process, and by $\nu_i$ the resulting probability measure on $P_i^\Gamma$.

We can then show that propagating along paths sampled according to $\nu_i$ converges in expectation to the solution of Problem (2.1) on node $i$. More precisely:

**Proposition 2.2.** Denoting by $j = u_{l-1}$ the last node reached before $p$ reaches $\Gamma$, we have:

$$ (2.3) \qquad f_*(i) = \mathbf{E}_{p \sim \nu_i} (\psi_{p_\Gamma^*}(g_j)), $$

where $p = p_\Gamma e_\Gamma$ with $e_\Gamma$ the last edge in $p$.

In other words, an unbiased estimate of $f_*(i)$ is obtained by drawing a random path from $i$ to $\Gamma$, and retropropagating the value of $g$ at the node just before $\Gamma$ ($j$ in Proposition 2.2), taking into account the rotations along the path. See Alg. 2.1. This formula is a discrete-time analog of a Theorem in [40], and we generalize it to weighted graphs in the Supplementary Material (Section 8, Proposition 8.1), building upon a connection-aware version of property (P2).

Alg. 2.1 is simple but may be computationally expensive: estimation on one node requires sampling $m$ paths, which needs to be repeated *for each* of the $|\mathcal{V}|$ nodes in $\mathcal{G}$ (a locality issue inherent to all Walk-on-Spheres-type algorithms [56]). We address these limitations in the remainder of this section and propose novel propagation-based estimators of $f_*$, allowing to update the Monte-Carlo estimates on *all the nodes at once*, by i) sampling and ii) propagating the signal, on specific subgraphs: the so-called MTSFs.
estimators of $f_*$, allowing to update the Monte-Carlo estimates on all the nodes at once, by sampling specific substructures of the graph, which is achieved with an expected number of steps of a random walk linear in the number of edges $|\mathcal{E}|$.

**Remark 2.3.** The restriction to complex signals on $\mathcal{G}$ and unitary connections (encoding 2D rotations) is mostly artificial, and Proposition 2.2 generalizes to transformations in e.g. $O(\mathbf{R}^d)$. This setting is encountered for instance for 3D molecule alignment in cryogenic electron microscopy (cryo-EM) [64], or in structure-from-motion problems [2], where the $\psi_e$'s act as 3D rotations (i.e. $\psi_e \in SO(\mathbf{R}^3)$).

**2.2. Multi-Type Spanning Forests: a first global estimator.** To build our global estimators, we rely on decompositions of the graph into Multi-Type Spanning Forests (MTSFs) [22]. A MTSF $\phi \subseteq \mathcal{V} \cup \mathcal{E}$ decomposes $\mathcal{G}$ into disjoint components that are either rooted trees, or unicycles. More precisely, a MTSF $\phi$ must have fixed cardinality $|\phi| = |\mathcal{V}|$, and is divided into (maximal) components $c_\phi \subseteq \phi$, which must be of one amongst two types:

* A rooted tree $c_\phi \subseteq \mathcal{V} \cup \mathcal{E}$, such that $c_\phi \cap \mathcal{E}$ is a connected cycle-free subset of edges, and $c_\phi \cap \mathcal{V}$ is reduced to a single node, connected to $c_\phi \cap \mathcal{E}$.

* A unicycle $c_\phi \subseteq \mathcal{E}$, which is a connected subset of edges containing exactly one cycle.

This structure is illustrated in Figure 2.1. We denote by $\mathcal{M}(\mathcal{G})$ the set of all MTSFs over the graph $\mathcal{G}$.

Note that in the absence of unicycles, a MTSF is a *rooted spanning forest* [3, 4], whereas a tree-free MTSF is a *spanning forest of unicycles* [35].

We will consider the distribution $\mathcal{D}_M$ over $\mathcal{M}(\mathcal{G})$ defined by:

$$ (2.4) \qquad \mathbf{P}_{\mathcal{D}_M}(\phi) \propto q^{|\phi \cap \mathcal{V}|} \prod_{C \in \mathcal{C}(\phi)} (2 - 2 \cos(\theta_C)), $$

where $\mathcal{C}(\phi)$ is the set of cycles belonging to the unicycles of $\phi$, and $\theta_C$ is the argument of the unitary complex number associated connection map $\psi_C$ obtained from a path traversing $C$ one time (as in Figure 1.1). Note that, while $\psi_C$ depends on the orientation of this path, $\cos(\theta_C)$ is insensitive to this orientation. Moreover, a cycle has a non-zero probability of being sampled if and only if it is *inconsistent* (if $C$ is perfectly consistent then $\theta_C = 0$ and $2 - 2 \cos(\theta_C) = 0$).

Distribution $\mathcal{D}_M$ is a rooted variant of the distribution considered in [22], and is in fact a DPP$^4$. We can now state our main theorem.

**Theorem 2.4.** Let $\phi$ be a MTSF sampled according to $\mathcal{D}_M$. Denote by a $\phi \xrightarrow{\phi} b$ the unique path from a to b in some tree of a MTSF $\phi$ and, if $v \in \mathcal{V}$ belongs to a rooted tree of $\phi$, by $r_\phi(v)$ the root of this tree. Consider the estimator

$$ (2.5) \qquad \tilde{f}(i, \phi, g) = \psi_{r_{\phi}(i)} \xrightarrow{\phi} i (g(r_{\phi}(i))) $$

propagating the value from $r_\phi(i)$ to $i$ if $i \in \mathcal{V}$ belongs to a rooted tree, and $\tilde{f}(i, \phi, g) = 0$ if $i$ lies in a unicycle. Then, we have:

$$ (2.6) \qquad f_*(i) = \mathbf{E}_{\mathcal{D}_M}(\tilde{f}(i, \phi, g)). $$

<footnote>$^4$So is the distribution of [22], but their non-combinatorial argument does not carry over to our rooted process. See the Section 9 of the Supplementary Material for a proof in this case</footnote>
**Remark 2.3.** The restriction to complex signals on $\mathcal{G}$ and unitary connections (encoding 2D rotations) is mostly artificial, and Proposition 2.2 generalizes to transformations in e.g. $O(\mathbf{R}^d)$. This setting is encountered for instance for 3D molecule alignment in cryogenic electron microscopy (cryo-EM) [76], or in structure-from-motion problems [2], where the $\psi_e$'s act as 3D rotations (i.e. $\psi_e \in SO(\mathbf{R}^3)$).

**2.2. Multi-Type Spanning Forests: a first global estimator.** To build our global estimators, we rely on decompositions of the graph into Multi-Type Spanning Forests (MTSFs) [22]. A MTSF $\phi \subseteq \mathcal{V} \cup \mathcal{E}$ decomposes $\mathcal{G}$ into disjoint components that are either rooted trees, or unicycles. More precisely, a MTSF $\phi$ must have fixed cardinality $|\phi| = |\mathcal{V}|$, and is divided into (maximal) components $c_\phi \subseteq \phi$, which must be either:

* a rooted tree $c_\phi \subseteq \mathcal{V} \cup \mathcal{E}$, such that $c_\phi \cap \mathcal{E}$ is a connected cycle-free subset of edges, and $c_\phi \cap \mathcal{V}$ is reduced to a single node, called root, connected to $c_\phi \cap \mathcal{E}$;

* a unicycle $c_\phi \subseteq \mathcal{E}$, which is a connected subset of edges containing exactly one cycle. See Figure 2.1 for an illustration. Note that in the absence of unicycles, a MTSF is a *rooted* spanning forest [3, 4], whereas a tree-free MTSF is a *spanning forest of unicycles* [41].

Denote by $\mathcal{M}(\mathcal{G})$ the set of all MTSFs over $\mathcal{G}$. Consider the distribution $\mathcal{D}_M$ over $\mathcal{M}(\mathcal{G})$:

$$ (2.4) \qquad \mathbf{P}_{\mathcal{D}_M}(\phi) \propto q^{|\phi \cap \mathcal{V}|} \prod_{C \in \mathcal{C}(\phi)} (2 - 2 \cos(\theta_C)), $$

where $\mathcal{C}(\phi)$ is the set of cycles belonging to the unicycles of $\phi$, and $\theta_C$ is the argument of the unitary complex number associated to the connection map $\psi_C$ obtained from a path traversing $C$ one time (as in Fig. 1.1). Note that, while $\psi_C$ depends on the path's orientation, $\cos(\theta_C)$ is insensitive to this orientation. Also, a cycle has a non-zero probability of being sampled if and only if it is *inconsistent* (if $C$ is perfectly consistent: $\theta_C = 0$ and $2 - 2 \cos(\theta_C) = 0$). $\mathcal{D}_M$ is a rooted variant of the distribution considered in [22], and is a DPP$^6$ over the set $\mathcal{V} \cup \mathcal{E}$. We can now state our main theorem.

**Theorem 2.4.** Let $\phi$ be a MTSF sampled according to $\mathcal{D}_M$. Denote by a $\phi \to b$ the unique path from a to b in some tree of a MTSF $\phi$ and, if $v \in \mathcal{V}$ belongs to a tree of $\phi$, by $r_\phi(v)$ the root of this tree. Consider the estimator

$$ (2.5) \qquad \tilde{f}(i, \phi, g) \begin{cases} = \psi_{r_\phi(i) \xrightarrow{\phi} i} (g(r_\phi(i))) & \text{if } i \text{ belongs to a tree} \\ = 0 & \text{if } i \text{ lies in a unicycle.} \end{cases} $$

When $i$ belongs to a tree, this estimator propagates the value of the signal $g$ from the root $r_\phi(i)$ to $i$, via the connections. When $i$ belongs to a unicycle, it is simply zero. Then, one has:

$$ (2.6) \qquad f_*(i) = \mathbf{E}_{\mathcal{D}_M}(\tilde{f}(i, \phi, g)), $$

i.e., the estimator is unbiased. Moreover, the variance of $\tilde{f}$ can be characterized as$^7$

$$ (2.7) \qquad \mathbf{E}_{\mathcal{D}_M} (\|\tilde{f}(\phi, g) - f_*\|^2) = \langle g, (I - K^2)g \rangle, $$

<footnote>$^6$So is the distribution of [22], but their non-combinatorial argument does not carry over to our rooted process. See the Section 9 of the Supplementary Material for a proof in this case *proof of what?*</footnote>

<footnote>$^7$Here, the variance of a complex-valued random variable $\tilde{z}$ is defined as $\text{var}(\tilde{z}) = \mathbf{E}(|\tilde{z}|^2) - |\mathbf{E}(\tilde{z})|^2$, and we compute in Eq. (2.7) the expected squared error $\mathbf{E}_{\mathcal{D}_M} (\|\tilde{f}(\phi, g) - f_*\|^2) = \sum_{i \in \mathcal{V}} \text{var}(\tilde{f}(i, \phi, g))$.</footnote>
The proof of Theorem 2.4 relies on the reformulation of $\mathcal{D}_M$ as a DPP, and an extension of arguments coming from both [35] and [51]. It is a special case of Theorem 9.1, proved in Section 9 of the Supplementary Material.

In a less technical phrasing, propagating for each tree in a MTSF $\phi$ the value of $g$ at its root $r$ to its other nodes $i$ (along path $r \xrightarrow{\phi} i$) results in an unbiased estimator of $f_*$. The estimation on nodes belonging to unicycles is then 0. See Algorithm 2.2 and Figure 2.1.

<figure><img src="image_10.png" /><figcaption>(a) A MTSF $\phi$ in purple, with two unicycles (top left, bottom) and one tree (right) with its root $r$ (large purple node).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Estimation: we set $\tilde{f}(i, \phi, g) = \psi_{\phi \to i}(g_r)$ for nodes in the tree (blue arrows), and $\tilde{f}(j, \phi, g) = 0$ on unicycles (teal crosses).</figcaption></figure>

<caption>Figure 2.1: Theorem 2.4 illustrated on the 4 × 4 grid graph, on one MTSF.</caption>

**Algorithm 2.2** MTSF-based estimator (Theorem 2.4).
1: $f_i \leftarrow 0 \quad \forall i \in \mathcal{V}$
2: **Repeat** $m$ times
3: Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_M$
4: **for** $i \in \mathcal{V}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $i$ belongs to a rooted tree of $\phi$ **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$f_i \leftarrow f_i + \psi_{r_\phi(i) \xrightarrow{\phi} i}(g(r_\phi(i)))$
7: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
8: **end for**
9: **Output** $\frac{1}{m}(f_1, \dots, f_{|\mathcal{V}|})$

**Remark 2.5.** If the connection $\Psi$ is trivial (i.e., $\psi_e = \text{id}_{C_{s_e}, C_{t_e}}$ for all edges), we recover the result from [51]. There are no inconsistent cycles in this case, and $\mathcal{D}_M$ reduces to another determinantal distribution $\mathcal{D}_F$ over rooted spanning forests $\phi \in \mathcal{V} \cup \mathcal{E} \ [3, 4]$.

Unlike Algorithm 2.1, Algorithm 2.2 allows to update the estimates on all the nodes for
<figure><img src="image_8.png" /></figure>

(a) A MTSF $\phi$ in purple, with two unicycles (top left, bottom) and one tree (right) with its root $r$ (large purple node).

<figure><img src="image_9.png" /></figure>

(b) Estimation: we set $\tilde{f}(i, \phi, g) = \psi_{\xrightarrow{\phi} i}(g_r)$ for nodes in the tree (blue arrows), and $\tilde{f}(j, \phi, g) = 0$ on unicycles (teal crosses).

<caption>Figure 2.1: Theorem 2.4 illustrated on the 4 × 4 grid graph, on one MTSF.</caption>

where $K = q(L_{\theta} + qI)^{-1}$ denotes the connection-aware Tikhonov smoothing operator.

The proof of Theorem 2.4 relies on the reformulation of $\mathcal{D}_M$ as a DPP, and an extension of arguments coming from both [41] and [61]. It is a special case of Theorem 9.1, which includes weighted graphs, and is proved in Section 9 of the Supplementary Material.

In a less technical phrasing, propagating for each tree in a MTSF $\phi$ the value of $g$ at its root $r$ to its other nodes $i$ (along path $r \xrightarrow{\phi} i$) results in an unbiased estimator of $f_*$. The estimation on nodes belonging to unicycles is then 0. See Alg. 2.2 and Figure 2.1.

<caption>**Algorithm 2.2** MTSF-based estimator (Theorem 2.4).</caption>

1: $f_i \leftarrow 0 \quad \forall i \in \mathcal{V}$
2: **Repeat** $m$ times
3: Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_M$ ▷ Alg. 2.3 can be used
4: **for** $i \in \mathcal{V}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $i$ belongs to a rooted tree of $\phi$ **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$f_i \leftarrow f_i + \psi_{r_{\phi}(i) \xrightarrow{\phi} i}(g(r_{\phi}(i)))$
7: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
8: **end for**
9: **Output** $\frac{1}{m}(f_1, \dots, f_{|\mathcal{V}|})$

**Remark 2.5.** If the connection $\Psi$ is trivial (i.e., $\psi_e = \text{id}_{C_{s_e}, C_{t_e}}$ for all edges), we recover the result from [61]. There are no inconsistent cycles in this case, and $\mathcal{D}_M$ reduces to another determinantal distribution $\mathcal{D}_F$ over rooted spanning forests $\phi \in \mathcal{V} \cup \mathcal{E} [3, 4]$.

Unlike Alg. 2.1, Alg. 2.2 allows to update the estimates on all the nodes for each sampled MTSF, in $O(|\mathcal{V}|)$ time. On the other hand, it becomes necessary to sample a MTSF $\phi$
each sampled MTSF, in $O(|V|)$ time. On the other hand, it becomes necessary to sample a MTSF $\phi$ according to $\mathcal{D}_M$. A first naive strategy is to use sampling algorithms designed for general DPPs, relying on the eigendecomposition of the DPP's kernel matrix (a Hermitian matrix of size $ (|V| + |E|) \times (|V| + |E|) $), which is much too expensive. In the next section, we recall the much more efficient random-walk-based algorithm proposed in [22], and analyze its expected computational cost.

**2.3. Sampling Multi-Type Spanning Forests.** Efficient sampling from $\mathcal{D}_M$ is possible thanks to loop-erased random walks, traditionally used in Wilson's algorithm for uniformly sampling spanning trees of a graph [69], under the following additional condition$^5$.

**Condition 2.6 (Weak-Inconsistency).** For all cycles $C$ in $G$, $\cos(\theta_C) \ge 0$.

We stress that this condition must hold for all the cycles of the graph, and not just for the cycles $C(\phi)$ in some MTSF $\phi \in \mathcal{M}(G)$. When Condition 2.6 is satisfied, $1 - \cos(\theta_C) \le 1$ defines a probability measure over the oriented cycles in $G$, which can be leveraged to efficiently sample MTSF $\phi$, as we will see in the following.

**Remark 2.7.** Even though efficient sampling from $\mathcal{D}_M$ can only be performed under Condition 2.6, Theorem 2.4 holds beyond this technical condition. Furthermore, in some applicable settings (e.g. directed graph signal processing [24], ranking [16, 71] or perceived luminance reconstruction [70]), this condition can be forced by tweaking the scale parameter $\gamma$.

**MTSF sampling algorithm.** The algorithm of [22], recalled here as Algorithm 2.3, samples a MTSF $\phi$ according to $\mathcal{D}_M$ by simulating multiple random walks on $G^\Gamma$, and constructs $\phi$ iteratively from these random paths. The precise way in which a path is turned into a part of $\phi$ relies on a loop-erasure procedure. The full sampling scheme is detailed in Algorithm 2.3, performing random walks that can be stopped in a number of ways: by being interrupted (i.e. reaching $\Gamma$), by building a cycle $\theta_C$ (kept with probability $1 - \cos(\theta_C)$), or by reaching a node already spanned by $\phi^6$. We suppose that the nodes of $G$ are arbitrarily ordered in a queue, and say that a node $i$ is spanned by $\phi$ if $i \in \phi$ or if $i$ is an endpoint of some edge $e \in \phi$. We also use a function randomSuccessor($u$) that, at each call, randomly outputs either $\Gamma$ with probability $\frac{A_{u,v}}{d_{u+q}}$, or some node $v$ with probability $\frac{A_{u,v}}{d_{u+q}}$.

One then proves the following by applying the arguments of [22] (which describes how to sample unrooted MTSFs), and keeping track of the roots.

**Proposition 2.8.** Suppose that Condition 2.6 holds. Then, Algorithm 2.3 outputs a MTSF $\phi$ distributed according to $\mathcal{D}_M$.

Let us now discuss the cost of Algorithm 2.3.

<footnote>$^5$Condition 2.6 has already been identified as a technical sampling condition for Wilson-like algorithms (in [33] and [22]). It also appears in a different guise in [34], under the name of trace-positivity, where it allows to simplify some technical statements.</footnote>

<footnote>$^6$Some implementation details of Algorithm 2.3 are not made completely explicit here, such as the precise way in which we detect cycles (instruction at line 11), or how to track the cycle-acceptance probabilities (i.e. $\theta_C$). Please refer to Section 10 of the Supplementary Material and our Julia implementation for more details</footnote>
according to $\mathcal{D}_M$. A first naive strategy is to use algorithms designed to sample generic DPPs. They however rely on the eigendecomposition of a kernel matrix of size $(|V|+|E|) \times (|V|+|E|)$, which is prohibitive. In the next section, we show how to take advantage of a very efficient random-walk-based algorithm (from [22]), and analyze its expected computational cost.

**2.3. Sampling Multi-Type Spanning Forests.** Efficient sampling from $\mathcal{D}_M$ is possible via loop-erased random walks, traditionally used in Wilson's algorithm to sample uniform spanning trees [82]. The simplest scenario arises when the following condition is verified$^8$.

**Condition 2.6 (Weak-Inconsistency).** For all cycles $C$ in $G$, $\cos(\theta_C) \ge 0$.

We stress that this condition constrains all the cycles of the graph, and not just the cycles $C(\phi)$ in some MTSF $\phi \in \mathcal{M}(G)$. If Condition 2.6 is satisfied, $1 - \cos(\theta_C) \le 1$ defines a probability measure over the oriented cycles in $G$, which can be leveraged to efficiently sample an MTSF $\phi$, as we detail in the following. If this condition does not hold, one can still efficiently sample MTSFs, by relying on importance sampling, as explained later in this section.

**MTSF sampling algorithm.** The algorithm of [22], recalled here as Alg. 2.3, samples a MTSF $\phi$ according to $\mathcal{D}_M$ by simulating multiple random walks on $G^\Gamma$, and constructs $\phi$ iteratively from these random paths. The precise way in which a path is turned into a part of $\phi$ relies on a loop-erasure procedure. The full sampling scheme is detailed in Alg. 2.3, performing random walks that can be stopped in a number of ways: by being interrupted (i.e. reaching $\Gamma$), by building a cycle $\theta_C$ (kept with probability $1 - \cos(\theta_C)$), or by reaching a node already spanned by $\phi^9$. We suppose that the nodes of $G$ are arbitrarily ordered in a queue, and say that a node $i$ is spanned by $\phi$ if $i \in \phi$ or if $i$ is an endpoint of some edge $e \in \phi$. We also use a function `randomSuccessor(u)` that, at each call, randomly outputs either $\Gamma$ with probability $\frac{q}{d_{u+q}}$, or some node $v$ with probability $\frac{A_{u,v}}{d_{u+q}}$.

One then proves the following by applying the arguments of [22] (which describes how to sample unrooted MTSFs), and keeping track of the roots.

**Proposition 2.7.** Suppose that Condition 2.6 holds. Then, Alg. 2.3 outputs a MTSF $\phi$ distributed according to $\mathcal{D}_M$.

Let us now discuss the cost of Alg. 2.3.

**Proposition 2.8.** The expected number of calls to `randomSuccessor` (i.e., the expected number of random walk steps), in Alg. 2.3, is upper-bounded by

$$ (2.8) \qquad \mathrm{tr} \left( (L + qI)^{-1} (D + qI) \right). $$

Also, Alg. 2.3 can be implemented with an expected time complexity upper-bounded by $O\left(\frac{|\mathcal{E}|}{q}\right)$.

<footnote>$^8$Condition 2.6 has already been identified as a technical sampling condition for Wilson-like algorithms (in [39] and [22]). It also appears in a different guise in [40], under the name of trace-positivity, where it allows to simplify some technical statements.</footnote>

<footnote>$^9$Some implementation details of Alg. 2.3 are not made completely explicit here, such as the precise way in which we detect cycles (instruction at line 11), or how to track the cycle-acceptance probabilities (i.e. $\theta_C$). Please refer to Section 10 of the Supplementary Material and our Julia implementation for more details</footnote>
<caption>**Algorithm 2.3** MTSF sampling algorithm [22].</caption>

1: $\phi \leftarrow \emptyset$
2: **while** $\phi$ not spanning **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Let $i \in \mathcal{V}$ be the first node in the queue not spanned by $\phi$
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u \leftarrow i$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p \leftarrow \epsilon$
6: &nbsp;&nbsp;&nbsp;&nbsp;▷ $u$ is the current node of the random walk
7: &nbsp;&nbsp;&nbsp;&nbsp;▷ $\epsilon$ the empty path
8: &nbsp;&nbsp;&nbsp;&nbsp;**while** ($u \neq \Gamma$) and ($p$ does not intersect $\phi$ or contain a cycle) **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u' \leftarrow \text{random\_successor}(u)$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Move to next node
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $u' \neq \Gamma$ **then**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$e \leftarrow (u, u')$, $p \leftarrow pe$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Add $e$ to the path $p$
14: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
15: &nbsp;&nbsp;&nbsp;&nbsp;**if** $p$ contains a cycle $C$ **then**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Remove this cycle from $p$ with probability $\cos(\theta_C)$
17: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
18: &nbsp;&nbsp;&nbsp;&nbsp;$$u_{\text{old}} \leftarrow u, u \leftarrow u'$$
19: &nbsp;&nbsp;&nbsp;&nbsp;▷ $u_{\text{old}}$ the previous node
20: &nbsp;&nbsp;&nbsp;&nbsp;**end while**
21: &nbsp;&nbsp;&nbsp;&nbsp;**if** $u = \Gamma$ **then**
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\phi \leftarrow \phi \cup p \cup u_{\text{old}}$
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Add the sampled path $p$ and the root $u_{\text{old}}$ to $\phi^7$
24: &nbsp;&nbsp;&nbsp;&nbsp;**else**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\phi \leftarrow \phi \cup p$
26: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
27: &nbsp;&nbsp;&nbsp;&nbsp;**end while**
28: &nbsp;&nbsp;&nbsp;&nbsp;**Output** $\phi$

**Proposition 2.9.** The expected number of random walk steps of Algorithm 2.3 is bounded from above by

$$ (2.7) \qquad \text{tr} \left( (L + qI)^{-1} (D + qI) \right). $$

Furthermore, Algorithm 2.3 can be implemented with $O\left(\frac{|\mathcal{E}|}{q}\right)$ expected time complexity.

See Section 10 for the proof. The bound is obtained by noting that Algorithm 2.3 necessitates fewer steps than the Wilson-like algorithm in [51], used to sample random spanning forests [3,4], with expected number of steps given by (2.7). Unlike the algorithm in [51] though, Algorithm 2.3 requires both tracking of the angular offsets accumulated along a cycle, and detecting said cycle, at additional computational cost. We discuss two possible implementations in Section 10 of the Supplementary Material, one of them resulting in the $O\left(\frac{|\mathcal{E}|}{q}\right)$ expected time complexity mentioned in Proposition 2.9. This translates to $O\left(\frac{|\mathcal{E}|}{q} + |\mathcal{V}|\right)$ expected runtime for Algorithm 2.2.

**Remark 2.10.** In case Condition 2.6 is not satisfied, one can still estimate $f_*$ up to a multiplicative constant ($af_*$ for some $a$), using the importance sampling strategy from [22].

<footnote n="7">⁷We abuse notations and denote here by $p$ the set of non-oriented edges in the path.</footnote>
<caption>**Algorithm 2.3** MTSF sampling algorithm [22].</caption>

1: $\phi \leftarrow \emptyset$
2: **while** $\phi$ not spanning **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Let $i \in \mathcal{V}$ be the first node in the queue not spanned by $\phi$
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u \leftarrow i$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ $u$ is the current node of the random walk
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p \leftarrow \epsilon$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ $\epsilon$ the empty path
6: &nbsp;&nbsp;&nbsp;&nbsp;**while** ($u \neq \Gamma$) and ($p$ does not intersect $\phi$ or contain a cycle) **do**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u' \leftarrow \text{random\_successor}(u)$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Move to next node
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $u' \neq \Gamma$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$e \leftarrow (u, u')$, $p \leftarrow pe$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Add $e$ to the path $p$
10: **end if**
11: **if** $p$ contains a cycle $C$ **then**
12: &nbsp;&nbsp;&nbsp;&nbsp;Remove this cycle from $p$ with probability $\cos(\theta_C)$
13: **end if**
14: $u_{\text{old}} \leftarrow u$, $u \leftarrow u'$
15: **end while**
16: **if** $u = \Gamma$ **then**
17: &nbsp;&nbsp;&nbsp;&nbsp;$\phi \leftarrow \phi \cup p \cup u_{\text{old}}$
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Add the sampled path $p$ and the root $u_{\text{old}}$ to $\phi^{10}$
18: **else**
19: &nbsp;&nbsp;&nbsp;&nbsp;$\phi \leftarrow \phi \cup p$
20: **end if**
21: **end while**
22: **Output** $\phi$

See Section 10 of the Supplementary Material for the proof. The bound is obtained by noting that Alg. 2.3 necessitates fewer steps than the Wilson-like algorithm in [61], used to sample random spanning forests [3,4], with expected number of steps exactly given by (2.8). Unlike the algorithm in [61] though, Alg. 2.3 requires both tracking of the angular offsets accumulated along a cycle, and detecting said cycle, at additional computational cost. We discuss two possible implementations in Section 10 of the Supplementary Material, one of them resulting in the $O(|\mathcal{E}|/q)$ expected time complexity mentioned in Proposition 2.8. This translates to $O(|\mathcal{E}|/q + |\mathcal{V}|)$ expected runtime for Alg. 2.2.

When Condition 2.6 does not hold: Importance Sampling. In case Condition 2.6 is not satisfied, one can still estimate $f_*$ using an importance sampling strategy<fnref n="11" />. Specifically, one can threshold the incoherence $2 - 2 \cos(\theta_C)$ and sample from

$$ (2.9) \qquad \mathbf{P}_{IS}(\phi) \propto q^{|\phi \cap \mathcal{V}|} \prod_{C \in \mathcal{C}(\phi)} \min(2, 2 - 2 \cos(\theta_C)) $$

<footnote n="10">¹⁰We abuse notations and denote here by $p$ the set of non-oriented edges in the path.</footnote>

<footnote n="11">¹¹The importance weights we consider have been proposed in [22], in the context of spectral sparsification.</footnote>
Specifically, one can threshold the incoherence $2 - 2 \cos(\theta_C)$ and sample from

$$ (2.8) \qquad \mathbf{P}_{IS}(\phi) \propto q^{\left|\phi \cap V\right|} \prod_{C \in \mathcal{C}(\phi)} \min(2, 2 - 2 \cos(\theta_C)) $$

using a straightforward variant of Algorithm 2.3. The estimation then uses the importance weights

$$ (2.9) \qquad w(\phi) = \prod_{C \in \mathcal{C}(\phi)} \max(1, 1 - \cos(\theta_C)). $$

Both the Feynman-Kac-based Algorithm 2.1 and the MTSF-based Algorithm 2.2 can be roughly understood as performing a random walk on the graph before stopping at some root node, and then propagating the value from this root to the starting point of the random walk. One main difference is that MTSFs allow to update the estimated values on all nodes jointly.

## 2.4. Variance Reduction for the MTSF-based estimator.
A paramount property of Monte-Carlo estimators is not only their unbiased behavior in expectation, but also their variance: they are significantly improved when used in conjunction with efficient variance reduction techniques. We propose two such improvements over $\tilde{f}$, based on the classical approaches of Rao-Blackwellization [7,55] and control variates [36] respectively, generalizing the variance reduction techniques introduced in [51,52].

Rao-Blackwellization. Rao-Blackwellization leverages the two laws of total expectation and variance, roughly stating that conditioning an estimator using another statistic extracted from the same sample still results in an unbiased estimator, with lower variance (see Section 11 of the Supplementary Material for more details in our case). Here, our technique consists in conditioning on the set of edges in $\phi$. This yields the following estimator:

$$ (2.10) \qquad \bar{f}(i, \phi, g) = \psi_{r_{\phi}(i) \to i} (h_{\phi}(r_{\phi}(i), g)) \quad \text{if } i \text{ belongs to a rooted tree,} $$

where $h_{\phi}$ is defined as

$$ (2.11) \qquad h_{\phi}(r, g) = \frac{\sum_{j \in c_{\phi}(r)} \psi_{r \to j} g(j)}{|c_{\phi}(r)|}, $$

and $c_{\phi}(r)$ is the set of nodes spanned by the tree containing $r$. If $i$ belongs to a unicycle, we once again set $\bar{f}(i, \phi, g) = 0$.

This results in Algorithm 2.4, which amounts to:

* Computing at the root $r$ of each tree the average $h_{\phi}(r, g)$ of the values of $g$ over the tree (obtained by propagating from each node $i$ in the tree to $r$).

* Propagating this average back to the other nodes of the tree.

As compared with Algorithm 2.2, this procedure can be implemented at little additional cost<fnref n="8" />.

<footnote n="8">⁸Note also the similarity to the Belief Propagation algorithm over trees [47].</footnote>
340 using a variant of Alg. 2.3, where a cycle is systematically kept, at line 12, if $\cos(\theta_C) \le 0$. 341 The estimation then uses the importance weights

$$ (2.10) \qquad w(\phi) = \prod_{C \in \mathcal{C}(\phi)} \max(1, 1 - \cos(\theta_C)), $$

343 which allows to estimate $cf_*$ for some unknown constant $c > 0$. Note that in the angular-synchronization pipeline we consider in Section 4, this is sufficient as we only need an estimation of $f_*$ up to a multiplicative constant. For smoothing purposes, one can further consider a 344 self-normalized importance-sampling estimator [58]: for $m$ MTSFs $\{\phi_k\}_{k=1}^m$ sampled according 345 to the distribution of Eq. (2.9), we set

$$ (2.11) \qquad \tilde{f}_{IS} = \frac{1}{\sum_{k=1}^{m} w(\phi_k)} \sum_{k=1}^{m} \tilde{f}(\phi_k, g), $$

349 which tends to $f_*$ as $m \to \infty$. This result holds regardless of the incoherence of the connection, 350 but the estimator is only guaranteed to be unbiased in the limit.

351 Both the Feynman-Kac-based Alg. 2.1 and the MTSF-based Alg. 2.2 can be roughly 352 described as performing a random walk on $\mathcal{G}$ before stopping at some root node, and then 353 retro-propagating the value from this root to the starting point of the random walk. The main 354 difference is that MTSFs allow to update the estimation on all nodes jointly. We will now see 355 how our estimator $\tilde{f}$ can be improved by implementing two variance reduction techniques.

356 **2.4. Variance Reduction for the MTSF-based estimator.** A paramount property of 357 Monte-Carlo estimators is not only their unbiased behavior in expectation, but also their 358 variance: they are significantly improved when used in conjunction with efficient variance 359 reduction techniques. We propose two such improvements over $\tilde{f}$, based on the classical ap- 360 proaches of *Rao-Blackwellization* [7,66] and *control variates* [42] respectively, generalizing the 361 variance-reduction techniques introduced in [61,63].

362 *Rao-Blackwellization.* Rao-Blackwellization leverages the two laws of total expectation and 363 variance, roughly stating that conditioning an estimator using another statistic extracted from 364 the same sample still results in an unbiased estimator, with lower variance (see Section 11 365 of the Supplementary Material for more details in our case). Here, our technique consists in 366 conditioning on the set of edges in $\phi$. This yields the following estimator:

$$ (2.12) \qquad \bar{f}(i, \phi, g) = \psi_{r_{\phi}(i) \to i} (h_{\phi}(r_{\phi}(i), g)) \quad \text{if } i \text{ belongs to a rooted tree,} $$

368 where $h_\phi$ is defined as

$$ (2.13) \qquad h_{\phi}(r, g) = \frac{\sum_{j \in c_{\phi}(r)} \psi_{r \to r} g(j)}{|c_{\phi}(r)|}, $$

370 and $c_\phi(r)$ is the set of nodes spanned by the tree containing $r$. If $i$ belongs to a unicycle, we 371 once again set $\bar{f}(i, \phi, g) = 0$. This results in Alg. 2.4, which amounts to:

372 • Computing at the root $r$ of each tree the average $h_\phi(r, g)$ of the values of $g$ over the 373 tree (obtained by propagating from each node $i$ in the tree to $r$).
<caption>**Algorithm 2.4** Rao-Blackwellized MTSF-based estimator (Theorem 2.4).</caption>

1: $f_i = 0 \forall i \in \mathcal{V}$, $h_r = 0 \forall r \in \mathcal{V}$
2: **Repeat** $m$ times
3: Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_M$
4: **for** $r \in \phi \cap \mathcal{V}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j \in c_\phi(r)$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$h_r = h_r + \psi_{j \to r}^\phi g(j)$
7: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
9: **for** $i \in \mathcal{V}$ belonging to a rooted tree of $\phi$ **do**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$f_i = f_i + \psi_{r_\phi(i) \to i}^\phi \left( \frac{h_{r_\phi(i)}}{|c_\phi(r_\phi(i))|} \right)$
11: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
12: **Output** $\frac{1}{m}(f_1, \dots, f_{|\mathcal{V}|})$

$\triangleright$ Via Algorithm 2.3

$\triangleright$ Propagate and average at the root

$\triangleright$ Propagate back

*Control variates.* Second is the introduction of control variates: an addition of another random quantity to $\bar{f}$, designed to have zero mean (so that the expectation remains unchanged), but resulting in an estimator with lower variance when designed properly. We propose to use a single *gradient-descent* step with parametrized step-size $\alpha$:

$$ (2.12) \qquad \hat{f}(\phi, g) = \bar{f}(\phi, g) - \alpha \mathbf{P}(q^{-1}(L_{\theta} + qI)\bar{f}(\phi, g) - g), $$

where $\mathbf{P} = (q^{-1}\mathbf{D} + I)^{-1}$ is a diagonal preconditioner for the system $q^{-1}(L_{\theta} + qI)f = g$, and $\alpha \in \mathbf{R}_+^*$. A good choice of step-size $\alpha$ is crucial in order to obtain a significant reduction of variance. We simply take $\alpha = 1$ in the following (see Section 11 of the Supplementary Material for empirical results backing this choice).

We prove in Section 11 of the Supplementary Material that generalizations of both $\bar{f}$ and $\hat{f}$ are unbiased estimators of $f_*$.

**Proposition 2.11.** $\mathbf{E}_{\mathcal{D}_M}(\bar{f}(i, \phi, g)) = \mathbf{E}_{\mathcal{D}_M}(\hat{f}(i, \phi, g)) = f_*(i)$

Both variance-reduction techniques are easy to implement, and do not incur large additional computational costs: The Rao-Blackwellization $\bar{f}$ has $\mathcal{O}(|\mathcal{V}|)$ additional cost, and the gradient-descent step in $\hat{f}$ entails a single matrix-vector multiplication, in $\mathcal{O}(|\mathcal{E}|)$ time.

**3. Numerical Results under Weak-Inconsistency.** We now analyze the behavior of the estimator proposed in Theorem 2.4, along with the improved versions discussed in Section 2.4, and compare their performance with (deterministic) conjugate-gradient-based solvers [56], on graphs with different topologies<fnref n="9" />. We perform experiments on the following connection model.

*Connection Model.* For a given graph $\mathcal{G}$ with $n$ nodes, we associate to each of its nodes $v$ an angle $\omega_v$ chosen uniformly in $[0, 2\pi)$, and we set

$$ (3.1) \qquad \theta_e = \omega_{t_e} - \omega_{s_e} + \eta\epsilon_e $$

<footnote n="9">⁹In this Section and in Section 4, we perform all our measurements using a single thread on a laptop with an intel i7-1185G7 processor.</footnote>
<caption>**Algorithm 2.4** Rao-Blackwellized MTSF-based estimator (Theorem 2.4).</caption>

1: $f_i = 0 \forall i \in \mathcal{V}$, $h_r = 0 \forall r \in \mathcal{V}$
2: **Repeat** $m$ times
3: Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_M$ ▷ Via Alg. 2.3
4: **for** $r \in \phi \cap \mathcal{V}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j \in c_\phi(r)$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$h_r = h_r + \psi_{j \to r}^\phi g(j)$ ▷ Propagate and average at the root
7: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
8: **end for**
9: **for** $i \in \mathcal{V}$ belonging to a rooted tree of $\phi$ **do**
10: &nbsp;&nbsp;&nbsp;&nbsp;$f_i = f_i + \psi_{r_\phi(i) \to i}^\phi (h_{r_\phi(i)})$ ▷ Propagate back
11: **end for**
12: **Output** $\frac{1}{m}(f_1, \dots, f_{|\mathcal{V}|})$

<footnote n="374">374 • Propagating this average back to the other nodes of the tree. As compared with Alg. 2.2, this procedure can be implemented at little additional cost¹².</footnote>

<footnote n="375">375 **Control variates.** Second is the introduction of control variates: an addition of another random quantity to $\bar{f}$, designed to have zero mean (so that the expectation remains unchanged), but resulting in an estimator with lower variance when designed properly. We propose to use a single *gradient-descent* step with parametrized step-size $\alpha$:</footnote>

<footnote n="376">376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
(2.14) &nbsp;&nbsp;&nbsp;&nbsp; $\hat{f}(\phi, g) = \bar{f}(\phi, g) - \alpha \mathbf{P}(q^{-1}(\mathbf{L}_{\theta} + q\mathbf{l})\bar{f}(\phi, g) - g)$, where $\mathbf{P} = (q^{-1}\mathbf{D} + \mathbf{l})^{-1}$ is a diagonal preconditioner for the system $q^{-1}(\mathbf{L}_{\theta} + q\mathbf{l})f = g$, and $\alpha \in \mathbf{R}_+^*$. A good choice of step-size $\alpha$ is crucial in order to obtain a significant reduction of variance. We simply take $\alpha = 1$ in the following (see Section 11 of the Supplementary Material for empirical results backing this choice).
We prove in Section 11 of the Supplementary Material that generalizations of both $\bar{f}$ and $\hat{f}$ are unbiased estimators of $f_*$, along with some results on the concentration of $\bar{f}$.
**Proposition 2.9.** *The variance-reduced estimators are unbiased:*</footnote>

<footnote n="387">387
388
389
390
391
392
393
(2.15) &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{E}_{\mathcal{D}_M}(\bar{f}(i, \phi, g)) = \mathbf{E}_{\mathcal{D}_M}(\hat{f}(i, \phi, g)) = f_*(i)$. Further, the variance of the Rao-Blackwellized estimator is given by</footnote>

<footnote n="387">387
388
389
390
391
392
393
(2.16) &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{E}_{\phi \sim \mathcal{D}_M} (\|\bar{f}(\phi, g) - f_*\|^2) = \langle g, (K - K^2)g \rangle$, where $K = q(\mathbf{L}_\theta + q\mathbf{l})^{-1}$, and we have the following finite-sample concentration bound. Let $\varepsilon, \delta \in (0, 1)$, and consider sampling $m$ MTSFs $\{\phi_k\}_{k=1}^m$. Then, provided that¹³</footnote>

<footnote n="398">398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
(2.17) &nbsp;&nbsp;&nbsp;&nbsp; $m \ge \frac{6}{\varepsilon^2} \log \left( \frac{|\mathcal{V}|}{\delta} \right)$, where $K = q(\mathbf{L}_\theta + q\mathbf{l})^{-1}$, and we have the following finite-sample concentration bound. Let $\varepsilon, \delta \in (0, 1)$, and consider sampling $m$ MTSFs $\{\phi_k\}_{k=1}^m$. Then, provided that¹³</footnote>

<footnote n="404">404
405
406
407
408
409
410
411
412
413
(2.17) &nbsp;&nbsp;&nbsp;&nbsp; $m \ge \frac{6}{\varepsilon^2} \log \left( \frac{|\mathcal{V}|}{\delta} \right)$,</footnote>

<footnote n="12">¹²Note also the similarity to the Belief Propagation algorithm over trees [54].</footnote>

<footnote n="13">¹³The factor 6 in Eq. (11.5) is not optimal: it is obtained from a more explicit bound appearing in the proof.</footnote>
for all edges e, with $\epsilon_e$ a perturbation uniformly distributed in $[-1, 1]$, and $\eta \in \mathbf{R}_+^*$ a scaling constant. We set $\eta = \frac{\pi}{2n}$ to ensure that Condition 2.6 is satisfied (so that Algorithm 2.3 in fact samples a MTSF according to $\mathcal{D}_M$).

**3.1. A First Runtime Experiment on Erdös-Rényi Graphs.** We first analyze the behavior of our estimator $\tilde{f}$ with respect to two parameters: the choice of regularization parameter $q$ and the mean degree $\bar{d}$ of the graph (controlling the graph's density). In our first experiment, we let $\bar{d}$ take values in $\{50, 100, 150, 200\}$, and $q$ take values equal to $q'\bar{d}$, with $q' \in \{10^{-3}, 10^{-1}, 1\}$. The complexity bound of Equation (2.7) becomes linear in $|\mathcal{V}|$ when using such a parametrization.

*Setup.* We generate Erdös-Rényi random graphs $\mathcal{G} \sim \text{ER}(n, p)$ of size $n = 10000$ for varying $p \in [0, 1]$, so that each edge $e$ independently appears with probability $p$. To control the density, we set $p = \frac{\bar{d}}{n-1}$ so that the expected mean degree of these random graphs is $\bar{d}$. For each such graph, we generate a random signal $f \in \mathbf{C}^V$ with independent complex Gaussian entries $f_v \sim \mathcal{N}_C(0, 1)$.

We then measure the running time of sampling one MTSF $\phi$ and applying the estimator $\tilde{f}$. As a reference, we also measure in the same manner the runtime of the matrix-vector multiplication $L_{\theta}f$, where $L_{\theta}$ is implemented as a sparse Hermitian matrix in CSC format. Note that this operation is the most expensive part of an iteration of the Conjugate-Gradient algorithm, and serves as a simple baseline to compare computation costs.

The results, depicted in Figure 3.1, are the average over 10000 time measurements, themselves averaged over 10 realizations of the graph $G$.

<figure><img src="image_12.png" /><figcaption>Figure 3.1: Runtime when varying $\bar{d}$, for different values of $q$.</figcaption></figure>
it holds that

$$ (2.18) \qquad \mathbf{P} \left( \forall g \in \mathcal{C}^{\mathcal{V}}, \left\| \left( \frac{1}{m} \sum_{k=1}^{m} \bar{f}(\phi_k, g) \right) - f_* \right\|_2 \le \varepsilon \|g\|_2 \right) \ge 1 - \delta. $$

Both variance-reduction techniques are easy to implement, and do not incur large additional computational costs: the Rao-Blackwellization $\bar{f}$ comes at $\mathcal{O}(|\mathcal{V}|)$ additional cost, and the gradient-descent step in $\hat{f}$ entails a single matrix-vector multiplication, in $\mathcal{O}(|\mathcal{E}|)$ time. Further, Eq. (11.6) shows that the number of forests required to reach a fixed precision grows logarithmically with the size of the graph, and that it is insensitive to the choice of $q$ and to the incoherence of the graph (which both affect the conditioning of the system, and the performance of conjugate-gradient solvers$^{14}$). Note though that these parameters do impact the sampling cost for MTSFs.

**3. Numerical Results under Weak-Inconsistency.** We now analyze the behavior of the estimator proposed in Theorem 2.4, along with the improved versions discussed in Section 2.4, and compare their performance with (deterministic) conjugate-gradient-based solvers [67], on graphs with different topologies$^{15}$. We perform experiments on the following connection model.

*Connection Model.* For a given graph $\mathcal{G}$ with $n$ nodes, we associate to each node $v$ an angle $\omega_v$ chosen uniformly at random in $[0, 2\pi)$, and we set

$$ (3.1) \qquad \theta_e = \omega_{t_e} - \omega_{s_e} + \eta \varepsilon_e $$

for all edges $e$, with $\varepsilon_e$ a perturbation uniformly distributed in $[-1, 1]$, and $\eta \in \mathbf{R}_+^*$ a scaling constant. We set $\eta = \frac{\pi}{2n}$ to ensure that Condition 2.6 is satisfied (so that Alg. 2.3 in fact samples a MTSF according to $\mathcal{D}_M$).

**3.1. A First Runtime Experiment on Erdös-Rényi Graphs.** We first analyze the behavior of our estimator $\tilde{f}$ with respect to two parameters: the choice of regularization parameter $q$ and the mean degree $\bar{d}$ of the graph (controlling the graph's density). In our first experiment, we let $\bar{d}$ take values in $\{50, 100, 150, 200\}$, and $q$ take values equal to $q'\bar{d}$, with $q' \in \{10^{-3}, 10^{-1}, 1\}$. The complexity bound of Eq. (2.8) becomes linear in $|\mathcal{V}|$ when using such a parametrization.

*Setup.* We generate Erdös-Rényi random graphs $\mathcal{G} \sim \text{ER}(n, p)$ of size $n = 10^4$ for varying $p \in [0, 1]$, so that each edge $e$ independently appears with probability $p$. To control the density, we set $p = \frac{\bar{d}}{n-1}$ so that the expected mean degree of these random graphs is $\bar{d}$. For each such graph, we generate a random signal $f \in \mathbf{C}^{\mathcal{V}}$ with independent complex Gaussian entries $f_v \sim \mathcal{N}_C(0, 1)$.

We then measure the running time of sampling one MTSF $\phi$ and applying the estimator $\tilde{f}$. As a reference, we also measure in the same manner the runtime of the matrix-vector multiplication $L_{\theta}f$, where $L_{\theta}$ is implemented as a sparse Hermitian matrix in CSC format. Note that this operation is the most expensive part of an iteration of the Conjugate-Gradient algorithm, and serves as a simple baseline to compare computation costs.

<footnote>$^{14}$See the *Discussion* paragraph of Section 3.2 for a reminder on the expected performance of CG</footnote>

<footnote>$^{15}$We perform all our measurements using a single thread on a laptop with an intel i7-1185G7 processor.</footnote>
**Analysis.** Two trends emerge in Figure 3.1. First, as $q'$ (hence $q$) increases, computing $\tilde{f}$ becomes less expensive, since the probability of the random walk stopping at some root node increases. Second, in contrast to the matrix-vector product $L_{\theta}f$, our estimator is less sensitive to the density of the graph: computing $\tilde{f}$ is more expensive than computing $L_{\theta}f$ for $\bar{d} = 50$, and systematically faster for $\bar{d} = 200$. This reflects the $O(|V|)$ expected complexity for $\tilde{f}$ (obtained with the parametrization $q = q'\bar{d}$).

**Remark 3.1.** Results are similar when varying the size $n$ of the graph (not shown).

## 3.2. Runtime-Precision Trade-offs in the context of Complex Graph Signal Denoising.

We compare the performance of our improved estimators $\bar{f}$ and $\hat{f}$ with conjugate-gradient methods [56]. The objective consider is to recover a signal $f_{\top} \in \mathbf{C}^{\mathcal{V}}$ given a noisy degradation $g = f_{\top} + \varepsilon$.

In the specific instance where $f_{\top} \in \mathbf{R}^{\mathcal{V}}$ and the connection on $\mathcal{G}$ is trivial ($\psi_e = \mathrm{id}_{\mathcal{C}_{se}}, \mathcal{C}_{te}$ for all $e \in \mathcal{E}$), a common assumption in the graph signal processing literature is that $f_{\top}$ is a smooth signal, that is, a linear combination of the first (low-frequency) eigenvectors of $L$ (associated to the smallest eigenvalues) [43, 53].

We here similarly assume that $f_{\top} \in \mathbf{C}^{\mathcal{V}}$ is B-bandlimited (i.e. $f_{\top} \in \mathrm{ran}(u_1, \dots, u_B)^{10}$, with $u_i$ the $i$-th eigenvector of $L_{\theta}$), so that solving the Tikhonov smoothing Problem (2.1):

$$ (3.2) \qquad \underset{f \in \mathbf{C}^{\mathcal{V}}}{\operatorname{argmin}} \quad q \|f - g\|_2^2 + \langle f, L_{\theta} f \rangle, $$

should allow to faithfully recover $f_{\top}$ from $g$, by penalizing high-frequency components in the optimal solution $f_* = q(L_{\theta} + qI)^{-1}g$. We refer the reader to Section 12 of the Supplementary Material for additional supporting arguments in this connection-aware setting.

Recall that this matrix inversion can be computed using a Cholesky decomposition (for an exact solution), or a conjugate-gradient-based iteration (for a high quality approximation). We compare our estimators with these methods, and consider two cost functions:

* The reconstruction error $e_r(f) = \|f - f_{\top}\|_2/n$, measuring the quality of the denoising.

* The approximation error $e_a(f) = \|f - f_*\|_2/n$, measuring the quality of the approximation of $f_*$.

Our estimators and the conjugate-gradient algorithms are respectively parametrized by the number of MTSFs used and the number of gradient steps, that we both denote by $m$.

**Setup.** For each graph, we generate a B-bandlimited signal $f_{\top}$, each $u_i$ being weighted by a random complex Gaussian value $a_i \sim \mathcal{N}_C(0, \sigma^2)$, with $\sigma^2$ such that the SNR is equal to 2. We then degrade $f_{\top}$ with some additive Gaussian noise $\varepsilon \sim \mathcal{N}_C(0, 1)$ (independently on each entry), and determine the optimal parameter $q_*$ for which $e_r(f_*) = e_r(q_*(L_{\theta} + q_*I)^{-1}g)$ is minimized$^{11}$, before measuring the errors and running time associated to different values of $m$ (taken amongst 10 logarithmically-spaced values from 1 to 100). The iterative algorithms are initialized at $g$.

We compare the following estimation strategies: the estimators $\bar{f}$ and $\hat{f}$, the conjugate-gradient descent with no preconditioner, with a simple diagonal preconditioner $P = (q^{-1}D +

<footnote>$^{10}$The linear subspace generated by the vectors $u_1, \dots, u_B$.</footnote>

<footnote>$^{11}$We perform our search in $(0, 30)$.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 3.1: Runtime when varying $\bar{d}$, for different values of $q$.</figcaption></figure>

The results, depicted in Figure 3.1, are the average over $10^4$ time measurements, themselves averaged over 10 realizations of the graph $G$.

**Analysis.** Two trends emerge in Figure 3.1. First, as $q'$ (hence $q$) increases, computing $\tilde{f}$ becomes less expensive, since the probability of the random walk stopping at some root node increases. Second, in contrast to the matrix-vector product $L_{\theta}f$, our estimator is less sensitive to the density of the graph: computing $\tilde{f}$ is more expensive than computing $L_{\theta}f$ for $\bar{d} = 50_2$ and systematically faster for $\bar{d} = 200$. This reflects the $O(|V|)$ expected complexity for $\tilde{f}$ (obtained with the parametrization $q = q'\bar{d})^{16}$.

## 3.2. Runtime-Precision Trade-offs in the context of Complex Graph Signal Denoising.

We compare the performance of our improved estimators $\tilde{f}$ and $\hat{f}$ with conjugate-gradient methods [67]. The objective is to recover a signal $f_{\top} \in \mathbf{C}^V$ given a noisy degradation $g = f_{\top} + \varepsilon$.

In the specific instance where $f_{\top} \in \mathbf{R}^V$ and the connection on $\mathcal{G}$ is trivial ($\psi_e = \text{id}_{C_{s_e}, C_{t_e}}$ for all $e \in \mathcal{E}$), a common assumption in the graph signal processing literature is that $f_{\top}$ is a smooth signal, that is, a linear combination of the first (low-frequency) eigenvectors of $L$ (associated to the smallest eigenvalues) [50, 64].

We here similarly assume that $f_{\top} \in \mathbf{C}^V$ is B-bandlimited (i.e. $f_{\top} \in \text{span}(u_1, ..., u_B)$), with $u_i$ the $i$-th eigenvector of $L_{\theta})$, so that solving the Tikhonov smoothing Problem (2.1):

$$ (3.2) \qquad \underset{f \in \mathbf{C}^V}{\operatorname{argmin}} \quad q \|f - g\|_2^2 + \langle f, L_{\theta} f \rangle, $$

should allow to faithfully recover $f_{\top}$ from $g$, by penalizing high-frequency components in the optimal solution $f_* = q(L_{\theta} + qI)^{-1}g$. We refer the reader to Section 13 of the Supplementary Material for additional supporting arguments in this connection-aware setting.

Recall that this matrix inversion can be computed using a Cholesky decomposition (for an exact solution), or a conjugate-gradient-based iteration (for a high quality approximation). We compare our estimators with these methods, and consider two cost functions:

* The reconstruction error $e_r(f) = \|f - f_{\top}\|_2/n$, measuring the quality of the denoising.

* The approximation error $e_a(f) = \|f - f_*\|_2/n$, measuring the quality of the approximation of $f_*$.

<footnote n="16">¹⁶Results are similar when varying the size $n$ of the graph (not shown).</footnote>
$l)^{-1}$, and with a CROUT ILU preconditioner<fnref n="12" />. Computing this high-quality ILU preconditioner is too expensive to be competitive with the other methods, and we only include these results on one type of graph, for illustration purposes.

Each runtime measurement is averaged over 100 runs. The results are averaged over 5 realizations of the noise $\varepsilon$ and 10 samples for each random graph model. We plot in Figure 3.2 the mean results.

*Graphs used.* We use the following graphs:

* An $\varepsilon$-graph obtained by sampling $|V| = 10000$ i.i.d. points $x_i$'s in $[0, 1]^3$, with an edge between $x_i$ and $x_j$ whenever $\|x_i - x_j\|_2 < 0.1$.

* A graph generated from a *Stochastic Block Model* (SBM), with $|V| = 10000$ nodes, each belonging to one of two communities $C_1$ and $C_2$ of size 5000 each. In this model, an edge is drawn randomly between two nodes $i, j$ with probability equal to $\frac{c_{k,l}}{n}$, where $k, l$ denote respectively the community label of $i$ and $j$. Here, we take $c_{1,1} = c_{2,2} = 36$ and $c_{1,2} = c_{2,1} = 4$, resulting in average degree $\bar{d} = 40$.

* A graph generated from a related *Degree-Corrected Stochastic Block Model* (DC-SBM 1), with two communities of size 5000 and an edge between nodes $i \in C_k$ and $j \in C_l$ present with probability proportional to $p_i p_j \frac{c_{k,l}}{n}$. $p_i$ is a randomly sampled positive real value, representing the intrinsic connectivity of node $i$, with $E(p_i) = 1$ and finite second moment [32]. Here, we take $c_{1,1} = c_{2,2} = 36$, $c_{1,2} = 4$, and the $p_i$'s are drawn from a normalized mixture of Gaussian distributions<fnref n="13" />, resulting in a graph with mean degree $\bar{d}$ close to 40. The objective of adding this model is to illustrate, at constant density, how the degree distribution affects the results.

* Another DC-SBM-graph (DC-SBM 2), with higher density (typically with an average degree more than 10 times that of the previous DC-SBM model). We take two communities of size 5000, $c_{1,1} = c_{2,2} = 480$, $c_{1,2} = 20$, and the $p_i$'s are drawn from another mixture of Gaussian distributions<fnref n="14" />. The objective of this second DC-SBM model is to illustrate how density affects the results.

* We also illustrate the results on a real-world graph: a relationship graph for internet Autonomous Systems (AS), recorded by the Center for Applied Internet Data Analysis (CAIDA) and provided in the SNAP datasets [41]<fnref n="15" />.

In the event of a randomly generated graph containing isolated nodes, we remove them so that the graph is connected (to make the interpretations simpler). We summarize in Table 3.1 the information concerning the graphs generated, and the associated optimal $q_*$'s (averaged over all realizations of the noise $\varepsilon$). All those graphs are endowed with a synthetic connection as specified in Equation (3.1). The parameters used for the (DC) SBM models ensures that the resulting graphs have a strong community structure.

<footnote n="12">¹²We use the implementations from the libraries: https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl, https://github.com/JuliaLinearAlgebra/Preconditioners.jl and https://github.com/haampie/IncompleteLU.jl. For the ILU preconditioner, we fix the drop threshold to 0.1.</footnote>

<footnote n="13">¹³Specifically, the connectivity parameters are independently sampled from a mixture of $N(50, 20)$, $N(500, 100)$ and $N(10000, 100)$, with weights of 0.59, 0.4 and 0.01 respectively, and then normalized.</footnote>

<footnote n="14">¹⁴Here from the mixture of $N(50, 20)$, $N(1000, 50)$, $N(5000, 100)$ and $N(10000, 100)$, with weights of 0.45, 0.1, 0.44 and 0.01.</footnote>

<footnote n="15">¹⁵We use the library available at https://github.com/JuliaGraphs/SNAPDatasets.jl.</footnote>
Our estimators and the conjugate-gradient algorithms are respectively parametrized by the number of MTSFs used and the number of gradient steps, that we both denote by $m$.

**Setup.** For each graph, we generate a $B$-bandlimited signal $f_T$, each $u_i$ being weighted by a random complex Gaussian value $a_i \sim N_C(0, 1)$. We then degrade $f_T$ with some additive Gaussian noise $\varepsilon \sim N_C(0, \sigma^2)$ (independently on each entry), with $\sigma^2$ such that the SNR is equal to 2, and determine the optimal parameter $q_*$ for which $e_r(f_*) = e_r(q_*(L_\theta + q_*I)^{-1}g)$ is minimized$^{17}$, before measuring the errors and running time associated to different values of $m$ (taken amongst 10 logarithmically-spaced values from 1 to 100). The iterative algorithms are initialized at $g$.

We compare the following estimation strategies: the estimators $\bar{f}$ and $\hat{f}$, the conjugate-gradient descent with no preconditioner, with a simple diagonal preconditioner $P = (q^{-1}D + I)^{-1}$, and with a CROUT ILU preconditioner$^{18}$. Computing this high-quality ILU preconditioner is too expensive to be competitive with the other methods, and we only include these results on one type of graph, for illustration purposes.

Each runtime measurement is averaged over 100 runs. The results are averaged over 5 realizations of the noise $\varepsilon$ and 10 samples for each random graph model. We plot in Figure 3.2 the mean results.

**Graphs used.** We use the following graphs:

* An $\varepsilon$-graph obtained by sampling $|V| = 10^4$ i.i.d. points $x_i$'s in $[0, 1]^3$, with an edge between $x_i$ and $x_j$ whenever $\|x_i - x_j\|_2 < 0.1$.

* A graph generated from a *Stochastic Block Model* (SBM), with $|V| = 10^4$ nodes, each belonging to one of two communities $C_1$ and $C_2$ of size 5000 each. In this model, an edge is drawn randomly between two nodes $i, j$ with probability$^{19}$ $\frac{c_{k,l}}{n}$, where $k, l$ denote respectively the community label of $i$ and $j$. Here, we take $c_{1,1} = c_{2,2} = 36$ and $c_{1,2} = c_{2,1} = 4$, resulting in an average degree $\bar{d} = 40$.

* A graph generated from a related *Degree-Corrected Stochastic Block Model* (DC-SBM 1), with two communities of size 5000 and an edge between nodes $i \in C_k$ and $j \in C_l$ present with probability $p_i p_j^{c_{k,l}/n}$. $p_i$ is a randomly sampled positive real value, representing the intrinsic *connectivity* of node $i$, with $E(p_i) = 1$ and finite second moment [37]. Here, we take $c_{1,1} = c_{2,2} = 36$, $c_{1,2} = 4$, and the $p_i$'s are drawn from a normalized mixture of Gaussian distributions$^{20}$, resulting in a graph with mean degree $\bar{d}$ close to 40. The objective of adding this model is to illustrate, at constant density, how the degree distribution affects the results.

* Another DC-SBM-graph (DC-SBM 2), with higher density (typically with an average degree more than 10 times that of the previous DC-SBM model). We take two communities of size 5000, $c_{1,1} = c_{2,2} = 480$, $c_{1,2} = 20$, and the $p_i$'s are drawn from another

<footnote>$^{17}$We perform our search in $(0, 30)$.</footnote>

<footnote>$^{18}$We use the implementations from the libraries: https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl, https://github.com/JuliaLinearAlgebra/Preconditioners.jl and https://github.com/haampie/IncompleteLU.jl. For the ILU preconditioner, we fix the drop threshold to 0.1.</footnote>

<footnote>$^{19}$It is usually assumed that $n > c_{k,l}$, so that we always have $\frac{c_{k,l}}{n} \in [0, 1]$. Most actual implementations use the probabilities $\min(\frac{c_{k,l}}{n}, 1)$ in case this assumption is not satisfied.</footnote>

<footnote>$^{20}$Specifically, the connectivity parameters are independently sampled from a mixture of $N(50, 20)$, $N(500, 100)$ and $N(10000, 100)$, with weights of 0.59, 0.4 and 0.01 respectively, and then normalized.</footnote>
<table><caption>Table 3.1: Experimental parameters associated to each graph. Values are averaged over all 10 realizations for random graph models.</caption><thead><tr><th>Graph</th><th>|V|</th><th>d&#773;</th><th>d<sub>min</sub></th><th>d<sub>max</sub></th><th>q<sub>*</sub></th></tr></thead><tbody><tr><td>&epsilon;-graph</td><td>10000</td><td>37.3</td><td>5.2</td><td>67.5</td><td>6.518</td></tr><tr><td>SBM</td><td>10000</td><td>39.87</td><td>17.7</td><td>66.6</td><td>21.04</td></tr><tr><td>DC-SBM 1</td><td>9833.7</td><td>33.4</td><td>1</td><td>906.1</td><td>2.634</td></tr><tr><td>DC-SBM 2</td><td>9932.8</td><td>418.8</td><td>1</td><td>1513.7</td><td>4.31</td></tr><tr><td>AS CAIDA</td><td>26475</td><td>4.032</td><td>1</td><td>2628</td><td>0.4808</td></tr></tbody></table>

For the ε-graph and the AS-graph, we arbitrarily set $B = 5$ when generating the bandlimited signal $f_T$. For (DC) SBMs, we take $B = 2$. It is known that the eigenvectors of the combinatorial Laplacian $L$ encode the community structure of these graphs (here, and taking $B = 2$ here results in signals coherent with the community structure (see Section 12 of the Supplementary Material for more on the first few eigenvectors of $L_θ$).

**Discussion.** Let us first comment on the approximation error for the ε-graph (Figure 3.2a, left). Both of our estimators converge linearly in log-log-scale with a mild slope, as expected for Monte-Carlo estimators. Any of the other three conjugate-gradient-based algorithms achieves a better approximation error within 3 steps than our estimators in 100 steps (as expected for Monte-Carlo convergence rates), and our methods are not competitive in this setting. However, in this denoising setting, the quality of the denoising is reflected in the reconstruction error (depicted in all the plots of Figure 3.2). Overall, we observe that:

(1) $\hat{f}$ consistently performs better than $f$, but how much better depends on the graph (see e.g. SBM and DC-SBM 1), and improvements are smaller for graphs with more heterogeneous degree distributions (DC-SBM graphs).

(2) Compared with MTSF-based estimators, the CG algorithms perform worse on graphs with heterogeneous degree distributions.

(3) CG solvers perform worse than our methods when the density increases.

CG's somewhat poor performance on DC-SBM 1 (bullet point (2)) is partly explained by the conditioning of the system considered: a bound on the condition number $\kappa$ is given by

$$ (3.3) \qquad \kappa = \frac{q_* + \lambda_n}{q_* + \lambda_1} \le \frac{q_* + 2 \max_{v \in V} d_v}{q_* + \lambda_1}, $$

where $\lambda_1 \le \dots \le \lambda_n$ are the eigenvalues of $L_\theta$, and where $\lambda_1$ is a measure of the quality of the optimal angular assignment (the lower $\lambda_1$ is, the better the quality of the synchronization), which is small for our synthetic connection (good synchronization due to low incoherence).

CG iterations are slower for high-density graphs ($O(|\mathcal{E}|)$ time per iteration), and $\bar{f}$ and $\hat{f}$ reach near-optimal reconstruction around 10 times faster than (preconditionned) CG on DC-SBM 2 (with expected complexities in $O(\frac{|\mathcal{E}|}{q} + |V|)$ and $O((1 + \frac{1}{q}) |\mathcal{E}| + |V|)$ respectively).

Our observations so far show that our estimators are competitive with standard methods,
<table><caption>Table 3.1: Experimental parameters associated to each graph. Values are averaged over all 10 realizations for random graph models.</caption><thead><tr><th>Graph</th><th>|V|</th><th>d&#773;</th><th>d<sub>min</sub></th><th>d<sub>max</sub></th><th>q<sub>*</sub></th></tr></thead><tbody><tr><td>&epsilon;-graph</td><td>10000</td><td>37.3</td><td>5.2</td><td>67.5</td><td>6.518</td></tr><tr><td>SBM</td><td>10000</td><td>39.87</td><td>17.7</td><td>66.6</td><td>21.04</td></tr><tr><td>DC-SBM 1</td><td>9833.7</td><td>33.4</td><td>1</td><td>906.1</td><td>2.634</td></tr><tr><td>DC-SBM 2</td><td>9932.8</td><td>418.8</td><td>1</td><td>1513.7</td><td>4.31</td></tr><tr><td>AS CAIDA</td><td>26475</td><td>4.032</td><td>1</td><td>2628</td><td>0.4808</td></tr></tbody></table>

mixture of Gaussian distributions<sup>21</sup>. The objective of this second DC-SBM model is to illustrate how density affects the results.

* We also illustrate the results on a real-world graph: a relationship graph for internet Autonomous Systems (AS), recorded by the Center for Applied Internet Data Analysis (CAIDA) and provided in the SNAP datasets [48]<sup>22</sup>.

In the event of a randomly generated graph containing isolated nodes, we remove them so that the graph is connected (to make the interpretations simpler). We summarize in Table 3.1 the information concerning the graphs generated, and the associated optimal $q_*$'s (averaged over all realizations of the noise $\epsilon$). All those graphs are endowed with a synthetic connection as specified in Eq. (3.1). The parameters used for the (DC) SBM models ensures that the resulting graphs have a strong community structure.

For the $\epsilon$-graph and the AS-graph, we arbitrarily set $B = 5$ when generating the bandlimited signal $f_T$. For (DC) SBMs, we take $B = 2$. It is known that the eigenvectors of the combinatorial Laplacian L encode the community structure of these graphs (here, and taking $B = 2$ here results in signals coherent with the community structure (see Section 13 of the Supplementary Material for more on the first few eigenvectors of $L_\theta$).

**Discussion.** Let us first comment on the approximation error for the $\epsilon$-graph (Figure 3.2a, left). Both of our estimators converge linearly in log-log-scale with a mild slope, as expected for Monte-Carlo estimators. Any of the other three conjugate-gradient-based algorithms achieves a better approximation error within 3 steps than our estimators in 100 steps (as expected for Monte-Carlo convergence rates), and our methods are not competitive in this setting. However, in this denoising setting, the quality of the denoising is reflected in the reconstruction error (depicted in all the plots of Figure 3.2). Overall, we observe that:

(1) $\hat{f}$ consistently performs better than $f$, but how much better depends on the graph (see e.g. SBM and DC-SBM 1), and improvements are smaller for graphs with more heterogeneous degree distributions (DC-SBM graphs);

(2) compared with MTSF-based estimators, the CG algorithms perform worse on graphs with heterogeneous degree distributions;

<footnote n="21"><sup>21</sup>Here from the mixture of $N(50, 20)$, $N(1000, 50)$, $N(5000, 100)$ and $N(10000, 100)$, with weights of 0.45, 0.1, 0.44 and 0.01.</footnote>

<footnote n="22"><sup>22</sup>We use the library available at https://github.com/JuliaGraphs/SNAPDatasets.jl.</footnote>
<figure><img src="image_13.png" /><figcaption>(a) ε-graph. Left: Approximation error. Right: Reconstruction error.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>(b) SBM Reconstruction error.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) DC-SBM 1 Reconstruction error.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) DC-SBM 2 Reconstruction error.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) AS CAIDA Reconstruction error.</figcaption></figure>

Figure 3.2: Runtime-precision trade-offs for Tikhonov smoothing with optimal *q*. Each data-point corresponds to a value of *m*. The vertical blue line records the average runtime of the Cholesky-based solver, the horizontal red line the average reconstruction error of the exact solution to the Tikhonov problem, and the horizontal yellow line the average reconstruction error of the noisy signal. Vertical bars represent standard deviations.
<figure><img src="image_11.png" /><figcaption>(a) ε-graph. Left: Approximation error. Right: Reconstruction error.</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>(b) SBM Reconstruction error.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) DC-SBM 1 Reconstruction error.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) DC-SBM 2 Reconstruction error.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) AS CAIDA Reconstruction error.</figcaption></figure>

Figure 3.2: Runtime-precision trade-offs for Tikhonov smoothing with optimal $q$. Each data-point corresponds to a value of $m$. The curve labelled MTSF corresponds to the estimator $\bar{f}$, and the curve labelled MTSF+GS to the estimator $\hat{f}$. The vertical blue line records the average runtime of the Cholesky-based solver, the horizontal red line the average reconstruction error of the exact solution to the Tikhonov problem, and the horizontal yellow line the average reconstruction error of the noisy signal. Vertical bars represent standard deviations.
especially for low-accuracy solutions, and in the following situations:

* On graphs resulting in poorly conditioned systems (with high-degree nodes and for low values of $q_*$).

* On graphs with higher density (to which MTSF-based methods less sensitive).

Further profiling also shows that the most expensive operations during the computation of $\bar{f}$ (and $\hat{f}$) are the instantiation of the data structures we use, and the successive memory accesses to the different $\theta_e$. These issues will likely benefit from further democratization of randomized computational schemes, with more efficient compiler optimizations and architectures specialized for these types of computations. Note also that the connection model used here is only representative of a subset of applications (e.g. [71], [16] or [24], see also Remark 1.1).

**4. Randomized Angular Synchronization.** We describe in this Section how our novel estimators can be successfully applied to angular synchronization. Recall from Section 1.2 that, to perform angular synchronization, we aim to minimize the incoherence

$$ (4.1) \qquad \underset{f \in U(C)^n}{\operatorname{argmin}} \langle f, L_{\theta} f \rangle, $$

which is NP-hard in general [72]. We go over existing approaches before introducing in Section 4.1 a randomized scheme based on a spectral relaxation of Problem (4.1).

*Spectral relaxation and other existing approaches.* Problem (4.1) is often relaxed to the following form [63, 71]:

$$ (4.2) \qquad \underset{\|f\|_2^2=n}{\operatorname{argmin}} \langle f, L_{\theta} f \rangle, $$

the solution of which is given by the eigenvector associated to the smallest eigenvalue of $L_\theta$. This solution differs from the exact solution of Problem (4.1) (even though guarantees on its quality exist, e.g. [23]), and can be computed using either an inverse power method, a Rayleigh quotient iteration, or a Lanczos iteration [57]. Most theoretical studies focus on Erdös-Rényi graphs, and suggest that spectral relaxations work best in small noise regimes (as compared with other algorithms, regarding the resulting error) [10, 40, 50].

Other existing methods include (see [18] for a similar discussion):

* Semi-definite relaxations [63], providing flexible theoretical tools and performance similar to spectral relaxations, but impractical past mid-sized instances ($n \approx 10^4$).

* The generalized power method proposed in [10], provably reaching the optimal solution of Problem (4.1) in the presence of (low) Gaussian noise, on complete graphs.

* Graph Neural Networks, with state-of-the-art performance for high noise [27, 30].

* Message-passing algorithms. For Gaussian noise, Approximate Message Passing has been conjectured to be statistically optimal among polynomial time algorithms [50], even at higher noise levels, but is limited to very dense graphs. Cycle-Edge Message-Passing allows exact recovery under a theoretical corruption model (with linear rate), but is more computationally-expensive than spectral relaxations [40].

* Descent techniques such as [46], with exact recovery under a (different) corruption model (see also [42]).
(3) CG solvers perform worse than our methods when the density increases (MTSF-based methods reach optimal reconstruction error 10 times faster on DC-SBM 2).

The relatively poor performance of CG on DC-SBM graphs is not unexpected. CG is widely known to be quite sensitive to ill-conditioning [67], which is why problem-specific preconditioners are often necessary. For instance, one can show that for a linear system of the form $Ax = b$ with solution $x_*$ (where A is self-adjoint and positive definite), the time-complexity required by CG to reach an approximation $\tilde{x}$ of $x_*$ such that $\frac{\|\tilde{x}-x_*\|_A^2}{\|x_*\|_A^2} \le \varepsilon$ is bounded in $O(t_A\sqrt{\kappa\log(\frac{1}{\varepsilon})})$, where $\|x\|_A^2 = \langle x, Ax \rangle$, $t_A$ is the cost of a matrix-vector product involving A, and $\kappa$ is the condition number of A [81]. In our case, $\kappa$ is bounded by<fnref n="23" />

$$ (3.3) \qquad \frac{q_* + 1 + d_{\max}}{q_* + \lambda_1} \le \kappa = \frac{q_* + \lambda_n}{q_* + \lambda_1} \le \frac{q_* + 2d_{\max}}{q_* + \lambda_1}, $$

where $\lambda_1 \le \dots \le \lambda_n$ are the eigenvalues of $L_\theta$ and $d_{\max}$ is the maximum degree. Note that $\lambda_1$ is a measure of the quality of the optimal angular assignment (the lower $\lambda_1$ is, the better the quality of the synchronization), which is small for our synthetic connection (good synchronization due to low incoherence). In the trivial-connection case, graphs with a highly-heterogenous degree distribution, as is the case here, tend to have poor conditioning; this follows for example from a perturbation argument (see, e.g., Thm.1 in [85]).

In addition to being sensitive to conditioning, CG is sensitive to graph density: each iteration has cost $t_A$ in $O(|\mathcal{E}|)$. On the other hand, we find in practice that our algorithm is much less sensitive to conditioning, and less sensitive to density. Showing this rigorously is certainly possible but would involves a lengthy exhaustion argument, depending on the comparisons between $q_*$, $\lambda_1$ and $d_{\max}$. Note that the concentration bound in Proposition 2.9 does not feature the condition number, which suggests that poorly-conditioned graphs do not require more Monte-Carlo replicates. In addition, setting $q \approx \bar{d}$ in the runtime bound of Proposition 2.8 gives $O(|V|)$ expected runtime for sampling MTSFs, which suggests lower sensitivity to graph density.

In general, we therefore expect MTSF-based solutions to perform better than CG on poorly-conditioned or high-density graphs.

Further profiling also shows that the most expensive operations during the computation of $\bar{f}$ (and $\hat{f}$) are the instantiation of the data structures we use, and the successive memory accesses to the different $\theta_e$. These issues will likely benefit from further democratization of randomized computational schemes [57], with more efficient optimizations and architectures specialized for these types of computations. Note also that the connection model used here is only representative of a subset of applications (e.g. [84], [16] or [26], see also Remark 1.1).

*Choice of q.* In practice, the value of $q_*$ is unknown, and the choice of the regularization parameter $q$ is a classical model-selection issue, that goes beyond the specific case of graphs, and can be approached from different perspectives.

* It may be inferred from an application-specific perspective (this is the case in the vector-field extension problem of [71]).

<footnote n="23">²³These bounds hold because $1+d_{\max} \le \lambda_n \le 2d_{\max}$, as we show in the Supplementary Material (Section 12).</footnote>
Note that message-passing algorithms and descent techniques (along with others we did not mention as well) can be applied to synchronization problems over more general groups, but are often not competitive on benchmarks for angular synchronization [27].

We will now discuss how to use our randomized estimators for eigenvector-computation.

**4.1. Proposed Approach.** One way to solve the spectral relaxation of Equation (4.2) is to perform an inverse power iteration, by setting $f_0 \in \mathbf{C}^V$ and iterating:

$$ (4.3) \qquad f_{r+1} = \frac{\mathbf{L}_{\theta}^{-1} f_r}{\|\mathbf{L}_{\theta}^{-1} f_r\|_2} $$

until convergence. Note that for this iteration to be well-defined, $\mathbf{L}_{\theta}$ needs to be invertible, i.e., the angular synchronization problems needs to be non-trivial.

Our approach stems from the observation that $\mathbf{L}_{\theta}$ and $q^{-1}(\mathbf{L}_{\theta} + q\mathbf{I})$ share the same eigenvectors. The power method applied to this regularized matrix consists in computing iterations of the form:

$$ (4.4) \qquad f_{r+1} = \frac{q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1} f_r}{\|q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1} f_r\|_2}, $$

where an estimation of $q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1} f_r$ can be obtained by our estimators $\bar{f}$ and $\hat{f}$ (or other deterministic algorithms). The convergence of this iteration is geometric with ratio

$$ (4.5) \qquad \frac{\mu_1}{\mu_2} = \frac{\lambda_2 + q}{\lambda_1 + q}, $$

where $\mu_1 \le \mu_2$ and $\lambda_1 \le \lambda_2$ denote respectively the two smallest eigenvectors of $q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1}$ and $\mathbf{L}_{\theta}$ [57]. The choice of $q$ is then a trade-off between:

* Low values of $q$ that induce smaller $\frac{\lambda_2+q}{\lambda_1+q}$ ratio, favoring faster convergence of the power iteration.

* Large values of $q$ that enable a faster sampling time of MTSFs (in $O(\|\varepsilon\|)$).

**Remark 4.1.** One could instead maximize $\langle f, A_{\theta}f \rangle$ in Equation (4.1). This formulation is common in theoretical works, that mostly focus on (often dense) Erdös-Rényi graphs, and applying the power method to $A_{\theta}$ is efficient on these graphs. More involved estimators are necessary in case the ratio $\alpha_1/\alpha_2$ of the top eigenvalues of $A_{\theta}$ is large. For trivial connections, such cases include: regular grids (with $\alpha_1/\alpha_2$ growing worse with the dimension), graphs with homogeneous spatial correlations (e.g., $\varepsilon$-graphs, nearest-neighbors-graphs), or graphs with homogeneous degree distributions and bottlenecks (poor expansion) such as in SBMs (for regular graphs, this is the Cheeger inequality). Our experiments (not shown) suggest that the same behavior arises for non-trivial connections. Further, Cramér-Rao bounds for angular synchronization suggest that angular synchronization is difficult on these graphs [11].

It is also possible to consider a normalized Laplacian $\tilde{\mathbf{L}}_{\theta} = \mathbf{D}^{-\frac{1}{2}}\mathbf{L}_{\theta}\mathbf{D}^{-\frac{1}{2}}$ in the spectral relaxation (4.2), which was proposed in [17, 71] and comes with similar guarantees [5]. Our approach is flexible enough to be used in this case. Indeed, our MTSF-based estimators extend
• From a statistical-signal-processing perspective, it may be estimated by inspecting different criteria, such as (for instance) Akaike's information criterion, the Bayesian information Criterion, leave-one-out cross validation, or Stein's unbiased risk estimator (in the case of graph Tikhonov regularization for trivial connections, this last criterion is discussed in detail in [61]).

The methods mentioned in this second point require the computation of the *effective degrees of freedom* of the model [31], for different values of $q$. Here, this quantity takes the form

$$ (3.4) \qquad s(q) = \operatorname{tr} (q(\mathbf{L}_{\theta} + q\mathbf{1})^{-1}). $$

In the trivial-connection setting, it is known that an unbiased estimator of $s(q)$ is obtained by counting the number of roots of MTSFs: $s(q) = \mathbf{E}_{\phi \sim D_M} (|\phi \cap V|)$ [6]; in conjunction with variance-reduction techniques, this method reaches state-of-the-art performance [62]. An important corollary of our work is the extension of this strategy to any connection, as we detail in the Supplementary Material (Remark 9.5).

**4. Randomized Angular Synchronization.** We describe in this section how our novel estimators can be successfully applied to angular synchronization. Recall from Section 1.2 that, to perform angular synchronization, we aim to minimize the incoherence

$$ (4.1) \qquad \underset{f \in U(\mathbf{C})^n}{\operatorname{argmin}} \langle f, \mathbf{L}_{\theta} f \rangle, $$

which is NP-hard in general [86]. We go over existing approaches before introducing in Section 4.1 a randomized scheme based on a spectral relaxation of Problem (4.1).

*Spectral relaxation and other existing approaches.* Problem (4.1) is often relaxed to the following form [75, 84]:

$$ (4.2) \qquad \underset{\|f\|_2^2=n}{\operatorname{argmin}} \langle f, \mathbf{L}_{\theta} f \rangle, $$

the solution of which is given by the eigenvector associated to the *smallest* eigenvalue of $\mathbf{L}_{\theta}$. This solution differs from the exact solution of Problem (4.1) (even though guarantees on its quality exist, e.g. [24]), and can be computed using either an inverse power method, a Rayleigh quotient iteration, or a Lanczos iteration [68]. Most theoretical studies focus on Erdös-Rényi graphs, and suggest that spectral relaxations work best in small noise regimes (as compared with other algorithms, regarding the resulting error) [10, 47, 60]. Other existing methods include (see [18] for a similar discussion):

• Semi-definite relaxations [75], providing flexible theoretical tools and performance similar to spectral relaxations, but impractical past mid-sized instances ($n \approx 10^4$).

• The generalized power method proposed in [10], provably reaching the optimal solution of Problem (4.1) in the presence of (low) Gaussian noise, on complete graphs.

• Graph Neural Networks, with state-of-the-art performance for high noise [32, 35].

• Message-passing algorithms. For Gaussian noise, Approximate Message Passing has been conjectured to be statistically optimal among polynomial time algorithms [60], even at higher noise levels, but is limited to very dense graphs. Cycle-Edge Message-Passing allows exact recovery under a theoretical corruption model (with linear rate), but is more computationally-expensive than spectral relaxations [47].
to quantities such as $f_0 = (\mathbf{L}_\theta + q\mathbf{D})^{-1}(q\mathbf{D})g'$ (see Section 9 of the Supplementary Material), which allows to estimate $q(\mathbf{L}_\theta + q\mathbf{I})^{-1}g$ by setting $g' = \mathbf{D}^{-\frac{1}{2}}g$ and leveraging the decomposition

$$ (4.6) \qquad q(\tilde{\mathbf{L}}_\theta + q\mathbf{I})^{-1}g = \mathbf{D}^{\frac{1}{2}}((\mathbf{L}_\theta + q\mathbf{D})^{-1}(q\mathbf{D}))g'. $$

## 4.2. Illustration: Intra-Class Denoising.
Before delving into precise computation time comparisons with state-of-the-art methods, we first illustrate our proposed method on a toy denoising problem, inspired from an application in cryo-EM [59, 74].

We consider $n$ copies $(I_i)_{i=1}^n$ of some image $I_*$ that have been rotated and degraded:

$$ (4.7) \qquad I_i = r_i(I_*) + \varepsilon, $$

with $r_i$ a rotation and $\varepsilon \sim \mathcal{N}(0, \sigma^2\mathbf{I})$ some Gaussian noise. The rotations $r_i$ are unknown, and the goal is to recover image $I_*$. This is a simplified version of the intra-class denoising problem in cryo-EM, where each (2D) image corresponds to a noisy projection of a (3D) molecule observed in an unknown orientation.

In the absence of rotations $r_i$, a solution consists in averaging the $I_i$'s, trivially recovering $I_*$ as $n \to \infty$. This strategy fails in our setting (due to the rotations), and we need to estimate the $r_i$'s before denoising. To do that, we first estimate angles $\theta_{i,j}$ between a subset of pair of images $\mathcal{E}$ (representing the edges of a graph) using image moments [26]$^{16}$, and perform angular synchronization to estimate the rotations $r_i$'s. We then rotate and average the images $I_i$'s.

We take $I_*$ the 256 × 256 Shepp-Logan phantom [61], $n = 1000$, $q = 10^{-3}$, $k = 20$ iterations of the power method, and randomly choose the underlying graph $\mathcal{G} \sim ER(n, p)$ with $p = \frac{5}{n}$. We uniformly sample rotations $r_i$ with angles in $(-\frac{\pi}{2}, \frac{\pi}{2})^{17}$. We display the recovered images in Figure 4.1 (solving the spectral relaxation (4.2) using both our method, with $m = 10$ MTSFs for the estimator $\bar{f}$ at each step, and an exact solver), for different noise levels $\sigma^2 \in \{1, 3, 5, 10\}$. We also plot the image reconstruction error $e_I(x) = \|I_* - x\|_2$ for $\sigma^2 = 5$ with varying values of $m$. Here, we have no guarantee that the connection resulting from the image-moment-based estimation is weakly-incoherent, and do not perform importance sampling. Instead, we always accept a cycle in Algorithm 2.3 if $\cos(\theta_C) \le 0$, and apply estimator $\bar{f}$ on the resulting MTSF.

Even though the connection may not be weakly-incoherent, the MTSF-based synchronization allows to obtain results of a quality very similar to exact synchronization on this problem, even at higher noise levels.

We will now investigate trade-offs between the quality of the synchronization and the runtime of our methods depending on the solver used to compute the power-method iteration in Equation 4.4, and the topology of the graph.

## 4.3. Numerical Evaluation.
We compare the performance of the method from Section 4.1 on synthetic graph data using four different iterative solvers: $\bar{f}$, $\hat{f}$ and a conjugate-gradient descent with and without diagonal preconditioning.

<footnote>$^{16}$Specifically, we estimate the orientation of each of the two images using the eigenvectors of the matrix of its covariant moments, and take $\theta_{i,j}$ the angle between these two sets of eigenvectors</footnote>

<footnote>$^{17}$We restrict the possible rotations so that we can use the simple image-moment-based registration, other methods could be used instead.</footnote>
• Descent techniques such as [53], with exact recovery under a (different) corruption model (see also [49]).

598

599

600

601

602

We will now discuss how to use our randomized estimators for eigenvector-computation.

603

604

605

**4.1. Proposed Approach.** One way to solve the spectral relaxation of Eq. (4.2) is to perform an inverse power iteration, by setting $f_0 \in \mathbf{C}^\mathcal{V}$ and iterating:

$$ (4.3) \qquad f_{r+1} = \frac{\mathbf{L}_{\theta}^{-1} f_r}{\|\mathbf{L}_{\theta}^{-1} f_r\|_2} $$

607

608

609

until convergence. Note that for this iteration to be well-defined, $L_\theta$ needs to be invertible, i.e., the angular-synchronization problem needs to be non-trivial. Depending on the topology of the graph $L_\theta$ may be poorly conditioned, yielding a very difficult robust estimation of $L_\theta^{-1} f_r$.

610

611

612

613

614

615

A classical workaround is to regularize the matrix $L_\theta$. Noting that the matrices $L_\theta^{-1}$ and its regularized version $(L_\theta+q\mathbf{l})^{-1}$ (with a better condition number) share the same eigenvectors, it is a classical result that the power method applied to a matrix of the form $(L_\theta+q\mathbf{l})^{-1}$ converges to the eigenvalue of $L_\theta$ closest to $-q$ [68]. In our case $q$ is positive, so that the power method applied to $(L_\theta+q\mathbf{l})$ converges to the desired spectral solution of angular synchronization, and the iteration reads

616

$$ (4.4) \qquad f_{r+1} = \frac{(L_{\theta} + q\mathbf{l})^{-1} f_r}{\|(L_{\theta} + q\mathbf{l})^{-1} f_r\|_2}, = \frac{q(L_{\theta} + q\mathbf{l})^{-1} f_r}{\|q(L_{\theta} + q\mathbf{l})^{-1} f_r\|_2}. $$

617

The convergence of this power method iteration is geometric with ratio

618

$$ (4.5) \qquad \frac{\mu_1}{\mu_2} = \frac{\lambda_2 + q}{\lambda_1 + q}, $$

619

620

where $\mu_1 \le \mu_2$ (resp. $\lambda_1 \le \lambda_2$) denote respectively the two largest (resp. smallest) eigenvalues of $q(L_\theta + q\mathbf{l})^{-1}$ (resp. $L_\theta$) [68]. The choice of $q$ is then a classical trade-off between:

621

622

• low values of $q$ that induce a smaller $\frac{\lambda_2+q}{\lambda_1+q}$ ratio, favoring faster convergence of the power iteration;

623

• large values of $q$ that result in a better conditioning of the system

624

625

626

627

628

629

630

631

632

**Remark 4.1.** One could instead maximize $\langle f, A_\theta f \rangle$ in Eq. (4.1). This formulation is common in theoretical works, that mostly focus on (often dense) Erdös-Rényi graphs, and applying the power method to $A_\theta$ is efficient on these graphs. However, it is often recommended to work with the smallest eigenvalue of $L_\theta$ rather than the largest of $A$ when the ratio $\alpha_1/\alpha_2$ of the top eigenvalues of $A_\theta$ is close to 1. For trivial connections, such cases include: regular grids (with $\alpha_1/\alpha_2$ growing worse with the dimension), graphs with homogeneous spatial correlations (e.g., $\varepsilon$-graphs, nearest-neighbors-graphs), or graphs with homogeneous degree distributions and bottlenecks (poor expansion) such as in SBMs (for regular graphs, this is a classical Cheeger inequality). Our experiments in Section 14 of the Supplementary Material suggest that the
<figure><img src="image_19.png" /><figcaption>Figure 4.1: Examples of recovery from noisy images, with different noise levels, using both the exact solution to the spectral relaxation (4.2) and the approach from Section 4.1.</figcaption></figure>
same observations can be made for non-trivial connections<sup>24</sup>. It is also possible to consider a normalized Laplacian $\tilde{L}_{\theta} = D^{-\frac{1}{2}}L_{\theta}D^{-\frac{1}{2}}$ in the spectral relaxation (4.2), which was proposed in [17, 84] and comes with similar guarantees [5].

Our approach simply consists in estimating each of the successive smoothings of Eq. (4.4) with MTSF-based estimators (we in fact recognize the term in $q(L_{\theta} + qI)^{-1}f_r$ as precisely what we can estimate with our MTSF-based estimators). In the remainder of this section, we illustrate the performance of our strategy on a toy problem, and compare it for different graph topologies to a smoothing performed by CG at each step of the power iteration. Further, we also illustrate in Section 14 of the Supplementary Material that, as compared to the power iteration applied to the matrix $A_{\theta}$, our strategy can result in very significant speed-ups on graphs for which the top spectral gap of $A_{\theta}$ is small (as discussed in Remark 4.1).

Note that our approach is flexible enough to compute the bottom eigenvector of the normalized connection Laplacian $\tilde{L}_{\theta}$ as well (Remark 4.1). Indeed, our MTSF-based estimators extend to quantities of the form $f_{\circ} = (L_{\theta} + qD)^{-1}(qD)g'$ (see Section 9 of the Supplementary Material), which allows to estimate $q(\tilde{L}_{\theta} + qI)^{-1}g$ by setting $g' = D^{-\frac{1}{2}}g$ and noting that:

$$ (4.6) \qquad q(\tilde{L}_{\theta} + qI)^{-1}g = D^{\frac{1}{2}}((L_{\theta} + qD)^{-1}(qD)) g'. $$

## 4.2. Illustration: Intra-Class Denoising.
Before delving into precise computation time comparisons with state-of-the-art methods, we first illustrate our proposed method on a toy denoising problem, inspired from an application in cryo-EM [70, 88].

We consider $n$ copies $(I_i)_{i=1}^n$ of some image $I_*$ that have been rotated and degraded:

$$ (4.7) \qquad I_i = r_i(I_*) + \varepsilon, $$

with $r_i$ a rotation and $\varepsilon \sim N(0, \sigma^2I)$ some Gaussian noise. The rotations $r_i$ are unknown, and the goal is to recover image $I_*$. This is a simplified version of the intra-class denoising problem in cryo-EM, where each (2D) image corresponds to a noisy projection of a (3D) molecule observed in an unknown orientation.

In the absence of rotations $r_i$, a solution consists in averaging the $I_i$'s, trivially recovering $I_*$ as $n \to \infty$. This strategy fails in our setting (due to the rotations), and we need to estimate the $r_i$'s before denoising. To do that, we first estimate angles $\theta_{i,j}$ between a subset of pair of images $\mathcal{E}$ (representing the edges of a graph) using image moments [29]<sup>25</sup>, and perform angular synchronization to estimate the rotations $r_i$'s. We then rotate and average the images $I_i$'s.

We take $I_*$ the 256×256 Shepp-Logan phantom [73], $n = 1000$, $q = 10^{-3}$, $k = 20$ iterations of the power method, and randomly choose the underlying graph $\mathcal{G} \sim ER(n, p)$ with $p = \frac{5}{n}$. We uniformly sample rotations $r_i$ with angles in $(-\frac{\pi}{2}, \frac{\pi}{2})^{26}$. We display the recovered images in Figure 4.1 (solving the spectral relaxation (4.2) using both our method, with $m = 10$ MTSFs for the estimator $\bar{f}$ at each step, and an exact solver), for different noise levels $\sigma^2 \in \{1, 3, 5, 10\}$.

<footnote n="24"><sup>24</sup>Further, Cramér-Rao bounds for angular synchronization suggest that angular synchronization is (statistically) difficult on these graphs [11].</footnote>

<footnote n="25"><sup>25</sup>Specifically, we estimate the orientation of each of the two images using the eigenvectors of the matrix of its covariant moments, and take $\theta_{i,j}$ the angle between these two sets of eigenvectors</footnote>

<footnote n="26"><sup>26</sup>We restrict the possible rotations so that we can use the simple image-moment-based registration, other methods could be used instead.</footnote>
**Setup.** We work with the SBM and DC-SBM graph models discussed in Section 3.2. We endow each graph $G$ with a random connection generated according to the model described in Section 3, for both coherent ($\eta = \frac{\pi}{2n}$) and incoherent connections ($\eta = \frac{\pi}{10}$), and aim to recover the vector $x$ with $x_i = e^{i\omega_i}$. Starting from $f_0$ taken uniformly in $U(\mathbf{C})^n$, we perform $k$ iterations of the power iteration of Equation (4.4), estimating the solution with $m$ MTSFs (resp. conjugate-gradient iterations) for each of our methods. We use $m = 3$ for coherent connections, and $m = 10$ for incoherent ones. For incoherent connections, we do use the importance sampling strategy from Remark 2.10$^{18}$. We average the synchronization errors

$$ (4.8) \qquad e_s(f) = \min_{r \in U(\mathbf{C})} \frac{\|f - rx\|_2}{n} $$

over 20 executions, and measure the mean runtimes over 100 runs. We set $q = 10^{-2} \times \bar{d}$ and take measurements for different values of $k$ (10 logarithmically-spaced values in $\{1, ..., 100\}$). The results are averaged over 10 realizations of the random graphs. See Figure 4.2.

We also take measurements for a Lanczos-iteration-based computation$^{19}$ (for the matrix $L_\theta$), and for a naive synchronization algorithm going as follows. First, fix a root node $v \in V$, then:

* sample a spanning tree of $G$ uniformly (a UST, using Wilson's algorithm [69]),

* propagate the value from $v$ to the other nodes (taking into account the offsets).

Note that this procedure does not depend on $k$ or $m$. Finally, we also evaluate a similar strategy propagating along a maximum spanning tree (MST), with respect to the edge-weights $w_{i,j} = |\cos(\theta_{i,j})|$ (a strategy inspired from surface reconstruction techniques such as [28]).

**Comments.** Results vary greatly depending on the graph and connection. For coherent connections ($\eta = \frac{\pi}{2n}$), we observe the following:

(1) Krylov-subspaces-based methods are sensitive to: all methods perform similarly on the SBM, but the CG-based power-methods and Lanczos iteration converge more slowly on the (less well-conditioned) DC-SBM 1 graph.

(2) Higher density results in significant slow-downs for Krylov-subspaces methods, and MTSF-based iterations are roughly 10 times faster than the preconditioned-CG iteration for equivalent precision on the DC-SBM 1.

For these coherent connections, the tree-based methods also achieve good synchronization quality, and are cheap to compute. This is no longer the case for incoherent connections (DC-SBM 1 with $\eta = \frac{\pi}{10}$): in this case, our estimators require more MTSFs (resp. CG iterations) to converge ($m = 10$), and achieve a synchronization error around 10 times smaller than tree-based methods.

**Remark 4.2.** We also experimented with two other estimation strategies (not shown): the spectral relaxation using the normalized Laplacian $\tilde{L}_\theta$, and a generalized-power-method-like algorithm, replacing the global normalization in Equation (4.4) by a component-wise normalization, so that each $f_k$ belongs to $U(\mathbf{C})^n$ (like in [10]). In our experiments, we did not observe any qualitative differences with the method described in Section 4.1.

<footnote>$^{18}$Here, there is no issue with estimating $af_*$ instead of $f_*$ due to the re-normalization in Equation (4.4).</footnote>

<footnote>$^{19}$We use the variant implemented in https://jutho.github.io/KrylovKit.jl/stable/man/eig/#KrylovKit.eigsolve.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 4.1: Examples of recovery from noisy images (left), with different noise levels, using both the exact solution to the spectral relaxation (4.2) and the approach from Section 4.1.</figcaption></figure>

We also plot the image reconstruction error $e_I(x) = \|I_* - x\|_2$ for $\sigma^2 = 5$ with varying values of $m$. Here, we have no guarantee that the connection resulting from the image-moment-based estimation is weakly-incoherent, and *perform importance sampling* according to the importance distribution of Eq. (2.9)$^{27}$, based on the estimator $\bar{f}$.

Even though the connection may not be weakly-incoherent, the MTSF-based synchronization allows to obtain results of a quality very similar to exact synchronization on this problem, even at higher noise levels.

We propose in Section 15 of the Supplementary Material a similar application on the MNIST dataset, where we do not control the incoherence of the connection, nor the noise.

We will now investigate trade-offs between the quality of the synchronization and the runtime of our methods depending on the solver used to compute the power-method iteration in Eq. (4.4), and the topology of the graph.

### 4.3. Numerical Evaluation.
We compare the performance of the method from Section 4.1 on synthetic graph data using four different iterative solvers: $\bar{f}$, $\hat{f}$ and a conjugate-gradient descent with and without diagonal preconditioning.

*Setup.* We work with the SBM and DC-SBM graph models discussed in Section 3.2. We endow each graph with a random connection generated according to the model described in Section 3, for both weakly-incoherent ($\eta = \pi/2n$) and incoherent connections ($\eta = \pi/10$), and aim to recover the vector $x$ with $x_i = e^{i\omega_i}$. Starting from $f_0$ taken uniformly in $U(\mathbf{C})^n$, we perform $k$ iterations of the power iteration of Eq. (4.4), estimating the solution with $m$ MTSFs (resp. conjugate-gradient iterations) for each method. We use $m = 3$ for weakly-incoherent connections, and $m = 10$ for incoherent ones. For incoherent connections, we use the importance-sampling distribution from Eq. (2.9). We average the synchronization errors

$$e_s(f) = \min_{r \in U(\mathbf{C})} \frac{\|f - rx\|_2}{n} \qquad (4.8)$$

<footnote>$^{27}$Here, there is no issue with estimating $af_*$ instead of $f_*$, due to the re-normalization in Eq. (4.4).</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 4.2: Runtime-precision trade-offs for angular synchronization for $m = 10$. Each data point corresponds to a value of $k$. The vertical yellow line records the average runtime of the Lanczos-based computation, the horizontal blue line the synchronization error for the Lanczos-based computation. The horizontal red (resp. orange) line is the mean synchronization error for the UST-based (resp. MST-based) synchronization, and the horizontal teal line the error for the random initialization.</figcaption></figure>

For spectral relaxations, our results suggest that MTSF-based solutions perform no-worse than deterministic solutions, and get comparatively better as the density of the graph increases. These techniques are mostly suitable for approximations, especially since previous studies (e.g. [50]) suggest that spectral-relaxation-based techniques mainly offer good-quality solutions at low noise-level, but are otherwise sub-optimal. One drawback of our strategy is that it requires to set a value for $q$: further refinements may include adapting our methodology to Rayleigh-quotient iterations [66], such that $q$ no longer needs to be fixed, and that better convergence rates may be obtained.

**5. Conclusion and Perspectives.** We proposed MTSF-based estimators for the problems of smoothing (in presence of a connection) and angular synchronization, reaching competitive performance on both problems, especially in low-precision regimes, and even *without* imple-
<figure><img src="image_18.png" /><figcaption>Figure 4.2: Runtime-precision trade-offs for angular synchronization. Each data point corresponds to a value of *k*. The vertical yellow line records the average runtime of the Lanczos-based computation, the horizontal blue line the synchronization error for the Lanczos-based computation. The horizontal red (resp. orange) line is the mean synchronization error for the UST-based (resp. MST-based) synchronization, and the horizontal teal line the error for the random initialization.</figcaption></figure>

over 20 executions, and measure the mean runtimes over 100 runs. We set $q = 10^{-2} \times \bar{d}$ and take measurements for different values of *k* (10 logarithmically-spaced values in {1, ..., 100}). The results are averaged over 10 realizations of the random graphs. See Figure 4.2.

We also take measurements for a Lanczos-iteration-based computation<sup>28</sup> (for the matrix $L_\theta$), and for a naive synchronization algorithm going as follows. First, fix a root node $v \in V$, then:

* sample a spanning tree of *G* uniformly (a UST, using Wilson's algorithm [82]),

* propagate the value from *v* to the other nodes (taking into account the offsets).

Note that this procedure does *not* depend on *k* or *m*. Finally, we also evaluate a similar strategy propagating along a *maximum spanning tree* (MST), with respect to the edge-weights $w_{i,j} = |\cos(\theta_{i,j})|$ (a strategy inspired from surface reconstruction techniques such as [33]).

<footnote n="28"><sup>28</sup>We use the one implemented in https://jutho.github.io/KrylovKit.jl/stable/man/eig/#KrylovKit.eigsolve.</footnote>
menting parallelization. Our estimators perform propagation along the branches of MTSFs and, as compared to standard deterministic methods, have another noteworthy advantage: they are insensitive to both the density of the graph *and* the conditioning of the system considered (which is typically bad for graphs with broad degree distributions). Such techniques can apply to a variety of other problems: one result we did not include is the extension of the smoothing estimator to the interpolation problem, which may find applications to *e.g.* synchronization in presence of anchor nodes.

This opens a number of research directions. From the theoretical side: could we lift the weak-inconsistency condition, following the development in Section 13 of the Supplementary Material? Could similar estimators be developed for $O(\mathbf{R}^d)$ synchronization? Instead of using a fixed number $m$ of MTSFs for the estimation (or CG iterations) of each of $k$ iterations of the power method, how should those $km$ MTSFs be used to achieve maximal precision?

From a more applied perspective, many applications of angular synchronization only consider the *exact* solution of the spectral relaxation, but how much loss in precision (and gain in speed) is actually acceptable? Our estimators can be seamlessly implemented on distributed systems (with communication complexity of the order of the number of steps in the sampling algorithm): could this be leveraged in practical applications? Finally, can random graph compositions such as MTSFs be useful in other iterative inference algorithms, for instance message-passing algorithms [40, 47]?

## REFERENCES

[1] B. ALEXEEV, A. S. BANDEIRA, M. FICKUS, AND D. G. MIXON, *Phase Retrieval with Polarization*, SIAM J. Imaging Sci., 7 (2014), pp. 35–66.

[2] M. ARIE-NACHIMSON, S. Z. KOVALSKY, I. KEMELMACHER-SHLIZERMAN, A. SINGER, AND R. BASRI, *Global motion estimation from point matches*, in 2012 Second international conference on 3D imaging, modeling, processing, visualization & transmission, IEEE, 2012, pp. 81–88.

[3] L. AVENA AND A. GAUDILLIÈRE, *Random spanning forests, Markov matrix spectra and well distributed points*, arXiv preprint arXiv:1310.1723, (2013).

[4] ——, *Two applications of random spanning forests*, Journal of Theoretical Probability, 31 (2018), pp. 1975–2004.

[5] A. S. BANDEIRA, A. SINGER, AND D. A. SPIELMAN, *A Cheeger inequality for the graph connection Laplacian*, SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 1611–1630.

[6] S. BARTHELME, N. TREMBLAY, A. GAUDILLIERE, L. AVENA, AND P.-O. AMBLARD, *Estimating the inverse trace using random forests on graphs*, in GRETSI 2019 - XXVIIème Colloque francophone de traitement du signal et des images, Lille, France, Aug. 2019.

[7] D. BLACKWELL, *Conditional expectation and unbiased sequential estimation*, The Annals of Mathematical Statistics, (1947), pp. 105–110.

[8] B. BOLLOBÁS, *Modern Graph Theory*, vol. 184, Springer Science & Business Media, 1998.

[9] P. BOUBOULIS, *Wirtinger's calculus in general Hilbert spaces*, arXiv preprint arXiv:1005.5170, (2010).

[10] N. BOUMAL, *Nonconvex Phase Synchronization*, SIAM J. Optim., 26 (2016), pp. 2355–2377.

[11] N. BOUMAL, A. SINGER, P.-A. ABSIL, AND V. D. BLONDEL, *Cramér-Rao bounds for synchronization of rotations*, Information and Inference: A Journal of the IMA, 3 (2014), pp. 1–39.

[12] D. BRYDGES, J. FRÖHLICH, AND E. SEILER, *On the construction of quantized gauge fields. I. General results*, Annals of Physics, 121 (1979), pp. 227–284.

[13] É. CARTAN, *Les groupes d'holonomie des espaces généralisés*, Acta mathematica, 48 (1926), pp. 1–42.

[14] S. CHEN, G. HUANG, G. PICCIOLI, AND L. ZDEBOROVÁ, *Planted X Y model: Thermodynamics and inference*, Physical Review E, 106 (2022), p. 054115.

[15] F. R. CHUNG, *Spectral graph theory*, vol. 92, American Mathematical Soc., 1997.
*Comments.* Results vary greatly depending on the graph and connection. For weakly-
incoherent connections ($\eta = \frac{\pi}{2n}$), we observe the following:

(1) Krylov-subspaces-based methods are sensitive to the conditioning : all methods perform similarly on the SBM, but not on the (less well-conditioned) DC-SBM graphs.

(2) Higher density results in significant slow-downs for Krylov-subspaces methods, and MTSF-based iterations have a much better performance in this case(DC-SBM 2).

(3) For weakly-inconsistent connections, the tree-based methods also achieve good synchronization quality, but this is no longer the case for incoherent connections.

If the first two factors combine (DC-SBM 2), our methods offer a “10 times” speed-up.

Remark 4.2. We also experimented with two other estimation strategies (not shown): the spectral relaxation using the normalized Laplacian $\overline{L}_{\theta}$, and a generalized-power-method-like algorithm, replacing the global normalization in Eq. (4.4) by a component-wise normalization, so that each $f_k$ belongs to $U(\mathbf{C})^n$ (like in [10]). In our experiments, we did not observe any qualitative differences with the method described in Section 4.1.

For spectral relaxations, our results suggest that MTSF-based solutions perform no-worse
than deterministic solutions, and get comparatively better as the density of the graph in-
creases. These techniques are mostly suitable for approximations, especially since previous
studies (e.g. [60]) suggest that spectral-relaxation-based techniques mainly offer good-quality
solutions at low noise-level, but are otherwise sub-optimal. One drawback of our strategy is
that it requires to set a value for q: further refinements may include adapting our methodology
to Rayleigh-quotient iterations [78], such that q no longer needs to be fixed, and that better
convergence rates may be obtained.

**5. Conclusion and Perspectives.** We propose MTSF-based estimators for connection-
aware smoothing and angular synchronization, reaching competitive performance on both
problems, especially in low-precision regimes, and even *without* implementing parallelization.
Our estimators perform propagation along the branches of MTSFs and, compared to standard
deterministic methods, have another noteworthy advantage: they are less sensitive to both
the density of the graph *and* the conditioning of the system considered (which is typically bad
for graphs with broad degree distributions). Our techniques can apply to a variety of other
problems: one result we did not include is the extension of the smoothing estimator to the
interpolation problem (*i.e.*, extending a signal that is known only on a few nodes to all the
graph), which may find applications to *e.g.* synchronization in presence of anchor nodes.

Our work opens a number of research directions. On the theoretical side: i) can the weak-
inconsistency condition (Condition 2.6) be lifted without resorting to importance sampling
strategies? We include a few preliminary developments in this direction in Section 16 of the
Supplementary Material. ii) Can similar estimators be developed for $O(\mathbf{R}^d)$ synchronization?
iii) Instead of using a fixed number $m$ of MTSFs for the estimation (or CG iterations) of each
of $k$ iterations of the power method, how should those $km$ MTSFs be used to achieve
maximal precision? Alternatively, we could re-use the set of $m$ MTSFs for each of the $k$
iterations of the power method.

From a more applied perspective, many applications of angular synchronization only con-
sider the *exact* solution of the spectral relaxation, but how much loss in precision (and gain in
[16] M. CUCURINGU, *Sync-Rank: Robust Ranking, Constrained Ranking and Rank Aggregation via Eigenvector and SDP Synchronization*, IEEE Trans. Netw. Sci. Eng., 3 (2016), pp. 58-79.

[17] M. CUCURINGU, Y. LIPMAN, AND A. SINGER, *Sensor network localization by eigenvector synchronization over the Euclidean group*, ACM Transactions on Sensor Networks (TOSN), 8 (2012), pp. 1-42.

[18] M. CUCURINGU AND H. TYAGI, *An extension of the angular synchronization problem to the heterogeneous setting*, arXiv preprint arXiv:2012.14932, (2020).

[19] M. DEREZINSKI AND M. W. MAHONEY, *Determinantal point processes in randomized numerical linear algebra*, Notices of the American Mathematical Society, 68 (2021), pp. 34-45.

[20] M. DEREZINSKI AND M. K. WARMUTH, *Reverse Iterative Volume Sampling for Linear Regression*, J. Mach. Learn. Res., 19 (2018), pp. 23:1-23:39.

[21] P. DRINEAS AND M. W. MAHONEY, *RandNLA: randomized numerical linear algebra*, Commun. ACM, 59 (2016), pp. 80-90.

[22] M. FANUEL AND R. BARDENET, *Sparsification of the regularized magnetic Laplacian with multi-type spanning forests*, arXiv preprint arXiv:2208.14797, (2022).

[23] F. FILBIR, F. KRAHMER, AND O. MELNYK, *On recovery guarantees for angular synchronization*, Journal of Fourier Analysis and Applications, 27 (2021), p. 31.

[24] S. FURUTANI, T. SHIBAHARA, M. AKIYAMA, K. HATO, AND M. AIDA, *Graph Signal Processing for Directed Graphs Based on the Hermitian Laplacian*, in Machine Learning and Knowledge Discovery in Databases - European Conference, ECML PKDD 2019, Würzburg, Germany, September 16-20, 2019, Proceedings, Part I, U. Brefeld, É. Fromont, A. Hotho, A. J. Knobbe, M. H. Maathuis, and C. Robardet, eds., vol. 11906 of Lecture Notes in Computer Science, Springer, 2019, pp. 447-463.

[25] A. GIRIDHAR AND P. R. KUMAR, *Distributed Clock Synchronization over Wireless Networks: Algorithms and Analysis*, in 45th IEEE Conference on Decision and Control, CDC 2006, San Diego, CA, USA, December 13-15, 2006, IEEE, 2006, pp. 4915-4920.

[26] R. C. GONZALEZ, *Digital image processing*, Pearson education india, 2009.

[27] Y. HE, G. REINERT, D. WIPF, AND M. CUCURINGU, *Robust Angular Synchronization via Directed Graph Neural Networks*, arXiv preprint arXiv:2310.05842, (2023).

[28] H. HOPPE, T. DEROSE, T. DUCHAMP, J. MCDONALD, AND W. STUETZLE, *Surface Reconstruction from Unorganized Points*, in Proceedings of the 19th annual conference on computer graphics and interactive techniques, 1992, pp. 71-78.

[29] J. B. HOUGH, M. KRISHNAPUR, Y. PERES, AND B. VIRÁG, *Determinantal Processes and Independence*, Probability Surveys, 3 (2006), pp. 206-229.

[30] N. JANCO AND T. BENDORY, *Unrolled algorithms for group synchronization*, IEEE Open Journal of Signal Processing, (2023).

[31] H. JAQUARD, M. FANUEL, P.-O. AMBLARD, R. BARDENET, S. BARTHELMÉ, AND N. TREMBLAY, *Smoothing Complex-Valued Signals on Graphs with Monte-Carlo*, in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1-5.

[32] B. KARRER AND M. E. NEWMAN, *Stochastic blockmodels and community structure in networks*, Physical review E, 83 (2011), p. 016107.

[33] A. KASSEL AND R. KENYON, *Random curves on surfaces induced from the Laplacian determinant*, The Annals of Probability, 45 (2017), pp. 932-964.

[34] A. KASSEL AND T. LÉVY, *Covariant symanzik identities*, Probability and Mathematical Physics, 2 (2021), pp. 419-475.

[35] R. KENYON, *SPANNING FORESTS AND THE VECTOR BUNDLE LAPLACIAN*, The Annals of Probability, 39 (2011), pp. 1983-2017.

[36] D. P. KROESE, T. TAIMRE, AND Z. I. BOTEV, *Handbook of monte carlo methods*, John Wiley & Sons, 2013.

[37] A. B. KUIJLAARS, *Convergence analysis of Krylov subspace iterations with methods from potential theory*, SIAM review, 48 (2006), pp. 3-40.

[38] A. KULESZA AND B. TASKAR, *Determinantal Point Processes for Machine Learning*, Found. Trends Mach. Learn., 5 (2012), pp. 123-286.

[39] R. KYNG AND Z. SONG, *A matrix chernoff bound for strongly rayleigh distributions and spectral sparsifiers from a few random spanning trees*, in 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), IEEE, 2018, pp. 373-384.
744 speed) is actually acceptable? Our estimators can be seamlessly implemented on distributed
745 systems (with communication complexity of the order of the number of steps in the sam-
746 pling algorithm): can this be leveraged in applications? Can random decompositions such as
747 MTSFs be useful in other inference algorithms, e.g., message-passing algorithms [47,54]?

748

REFERENCES

[1] B. ALEXEEV, A. S. BANDEIRA, M. FICKUS, AND D. G. MIXON, *Phase Retrieval with Polarization*, SIAM J. Imaging Sci., 7 (2014), pp. 35–66.

[2] M. ARIE-NACHIMSON, S. Z. KOVALSKY, I. KEMELMACHER-SHLIZERMAN, A. SINGER, AND R. BASRI, *Global motion estimation from point matches*, in 3DIMPVT, IEEE, 2012.

[3] L. AVENA AND A. GAUDILLIÈRE, *Random spanning forests*, Markov matrix spectra and well distributed points, arXiv preprint arXiv:1310.1723, (2013).

[4] ——, *Two applications of random spanning forests*, Journal of Theoretical Probability, (2018).

[5] A. S. BANDEIRA, A. SINGER, AND D. A. SPIELMAN, *A Cheeger inequality for the graph connection Laplacian*, SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 1611–1630.

[6] S. BARTHELME, N. TREMBLAY, A. GAUDILLIERE, L. AVENA, AND P.-O. AMBLARD, *Estimating the inverse trace using random forests on graphs*, in GRETSI 2019, Lille, France, Aug. 2019.

[7] D. BLACKWELL, *Conditional expectation and unbiased sequential estimation*, The Ann. of Mathematical Statistics, (1947).

[8] B. BOLLOBÁS, *Modern Graph Theory*, vol. 184, Springer Science & Business Media, 1998.

[9] P. BOUBOLIS, *Wirtinger's calculus in general Hilbert spaces*, arXiv preprint arXiv:1005.5170, (2010).

[10] N. BOUMAL, *Nonconvex Phase Synchronization*, SIAM J. Optim., 26 (2016), pp. 2355–2377.

[11] N. BOUMAL, A. SINGER, P.-A. ABSIL, AND V. D. BLONDEL, *Cramér-Rao bounds for synchronization of rotations*, Information and Inference: A Journal of the IMA, 3 (2014), pp. 1–39.

[12] D. BRYDGES, J. FRÖHLICH, AND E. SEILER, *On the construction of quantized gauge fields. I. General results*, Annals of Physics, 121 (1979), pp. 227–284.

[13] É. CARTAN, *Les groupes d'holonomie des espaces généralisés*, Acta mathematica, 48 (1926), pp. 1–42.

[14] S. CHEN, G. HUANG, G. PICCIOLI, AND L. ZDEBOROVÁ, *Planted X Y model: Thermodynamics and inference*, Physical Review E, 106 (2022), p. 054115.

[15] F. R. CHUNG, *Spectral graph theory*, vol. 92, American Mathematical Soc., 1997.

[16] M. CUCURINGU, *Sync-Rank: Robust Ranking, Constrained Ranking and Rank Aggregation via Eigenvector and SDP Synchronization*, IEEE Trans. Netw. Sci. Eng., 3 (2016), pp. 58–79.

[17] M. CUCURINGU, Y. LIPMAN, AND A. SINGER, *Sensor network localization by eigenvector synchronization over the Euclidean group*, ACM Transactions on Sensor Networks (TOSN), 8 (2012), pp. 1–42.

[18] M. CUCURINGU AND H. TYAGI, *An extension of the angular synchronization problem to the heterogeneous setting*, arXiv preprint arXiv:2012.14932, (2020).

[19] M. DEREZINSKI AND M. W. MAHONEY, *Determinantal point processes in randomized numerical linear algebra*, Notices of the American Mathematical Society, 68 (2021), pp. 34–45.

[20] M. DEREZINSKI AND M. K. WARMUTH, *Reverse Iterative Volume Sampling for Linear Regression*, J. Mach. Learn. Res., 19 (2018), pp. 23:1–23:39.

[21] P. DRINEAS AND M. W. MAHONEY, *RandNLA: randomized numerical linear algebra*, Com. ACM, (2016).

[22] M. FANUEL AND R. BARDENET, *Sparsification of the regularized magnetic Laplacian with multi-type spanning forests*, arXiv preprint arXiv:2208.14797, (2022).

[23] M. FANUEL AND R. BARDENET, *On the Number of Steps of CyclePopping in Weakly Inconsistent U(1)-Connection Graphs*. 2024.

[24] F. FILBIR, F. KRAHMER, AND O. MELNYK, *On recovery guarantees for angular synchronization*, Journal of Fourier Analysis and Applications, 27 (2021), p. 31.

[25] R. FORMAN, *Determinants of Laplacians on graphs*, Topology, 32 (1993), pp. 35–46.

[26] S. FURUTANI, T. SHIBAHARA, M. AKIYAMA, K. HATO, AND M. AIDA, *Graph Signal Processing for Directed Graphs Based on the Hermitian Laplacian*, in ECML PKDD, 2019.

[27] D. GABAY AND B. MERCIER, *A dual algorithm for the solution of nonlinear variational problems via*
[40] G. LERMAN AND Y. SHI, *Robust Group Synchronization via Cycle-Edge Message Passing*, Foundations of Computational Mathematics, 22 (2022), pp. 1665-1741.

[41] J. LESKOVEC AND A. KREVL, *SNAP Datasets: Stanford Large Network Dataset Collection*. http://snap.stanford.edu/data, June 2014.

[42] H. LIU, X. LI, AND A. M.-C. SO, *ReSync: Riemannian Subgradient-based Robust Rotation Synchronization*, arXiv preprint arXiv:2305.15136, (2023).

[43] P. LORENZO, S. BARBAROSSA, AND P. BANELLI, *Sampling and recovery of graph signals*, in Cooperative and Graph Signal Processing, Elsevier, 2018, pp. 261-282.

[44] O. MACCHI, *The coincidence approach to stochastic point processes*, Advances in Applied Probability, 7 (1975), pp. 83-122.

[45] P. MARTINSSON AND J. A. TROPP, *Randomized numerical linear algebra: Foundations and algorithms*, Acta Numer., 29 (2020), pp. 403-572.

[46] T. MAUNU AND G. LERMAN, *Depth descent synchronization in SO (d)*, International journal of computer vision, 131 (2023), pp. 968-986.

[47] M. MEZARD AND A. MONTANARI, *Information, physics, and computation*, Oxford University Press, 2009.

[48] L. A. S. MOREIRA, A. L. L. RAMOS, M. L. R. DE CAMPOS, J. A. APOLINÁRIO, AND F. G. SERRENHO, *A Graph Signal Processing Approach to Direction of Arrival Estimation*, in 27th European Signal Processing Conference, EUSIPCO 2019, A Coruña, Spain, September 2-6, 2019, IEEE, 2019, pp. 1-5.

[49] M. E. MULLER, *Some continuous Monte Carlo methods for the Dirichlet problem*, The Annals of Mathematical Statistics, (1956), pp. 569-589.

[50] A. PERRY, A. S. WEIN, A. S. BANDEIRA, AND A. MOITRA, *Message-Passing Algorithms for Synchronization Problems over Compact Groups*, Communications on Pure and Applied Mathematics, 71 (2018), pp. 2275-2322.

[51] Y. Y. PILAVCI, P. AMBLARD, S. BARTHELMÉ, AND N. TREMBLAY, *Graph Tikhonov Regularization and Interpolation Via Random Spanning Forests*, IEEE Trans. Signal Inf. Process. over Networks, 7 (2021), pp. 359-374.

[52] Y. Y. PILAVCI, P.-O. AMBLARD, S. BARTHELMÉ, AND N. TREMBLAY, *Variance Reduction in Stochastic Methods for Large-Scale Regularized Least-Squares Problems*, in 2022 30th European Signal Processing Conference (EUSIPCO), IEEE, 2022, pp. 1771-1775.

[53] G. PUY, N. TREMBLAY, R. GRIBONVAL, AND P. VANDERGHEYNST, *Random sampling of bandlimited signals on graphs*, Applied and Computational Harmonic Analysis, 44 (2018), pp. 446-475.

[54] T. C. RAIA, M. G. P. THOMAS, F. G. SERRENHO, AND J. A. APOLINÁRIO JR, *Gsp-based doa estimation for a multimission radar*, Proceedings of the XXXVIII Simpósio Brasileiro de Telecomunicações e Processamento de Sinais (SBrT2020), Florianópolis, Brazil, (2020), pp. 22-25.

[55] C. R. RAO, *Information and the accuracy attainable in the estimation of statistical parameters*, in Breakthroughs in Statistics: Foundations and basic theory, Springer, 1992, pp. 235-247.

[56] Y. SAAD, *Iterative methods for sparse linear systems*, SIAM, 2003.

[57] ——, *Numerical methods for large eigenvalue problems: revised edition*, SIAM, 2011.

[58] R. SAWHNEY AND K. CRANE, *Monte Carlo geometry processing: A grid-free approach to PDE-based methods on volumetric domains*, ACM Transactions on Graphics, 39 (2020).

[59] S. H. SCHERES, *RELION: implementation of a Bayesian approach to cryo-EM structure determination*, Journal of structural biology, 180 (2012), pp. 519-530.

[60] SHARP, NICHOLAS AND SOLIMAN, YOUSUF AND CRANE, KEENAN, *The vector heat method*, ACM Transactions on Graphics (TOG), 38 (2019), pp. 1-19.

[61] L. A. SHEPP AND B. F. LOGAN, *The Fourier reconstruction of a head section*, IEEE Transactions on nuclear science, 21 (1974), pp. 21-43.

[62] D. I. SHUMAN, S. K. NARANG, P. FROSSARD, A. ORTEGA, AND P. VANDERGHEYNST, *The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains*, IEEE signal processing magazine, 30 (2013), pp. 83-98.

[63] A. SINGER, *Angular synchronization by eigenvectors and semidefinite programming*, Applied and computational harmonic analysis, 30 (2011), pp. 20-36.

[64] A. SINGER AND H.-T. WU, *Vector diffusion maps and the connection Laplacian*, Communications on pure and applied mathematics, 65 (2012), pp. 1067-1144.

[65] D. A. SPIELMAN AND N. SRIVASTAVA, *Graph sparsification by effective resistances*, in Proceedings of the
finite element approximation, Computers & mathematics with applications, 2 (1976), pp. 17–40.

[28] A. GIRIDHAR AND P. R. KUMAR, *Distributed Clock Synchronization over Wireless Networks: Algorithms and Analysis*, in CDC, IEEE, 2006.

[29] R. C. GONZALEZ, *Digital image processing*, Pearson education india, 2009.

[30] L. GULIKERS, M. LELARGE, AND L. MASSOULIÉ, *Non-backtracking spectrum of degree-corrected stochastic block models*, arXiv preprint arXiv:1609.02487, (2016).

[31] T. HASTIE, R. TIBSHIRANI, J. H. FRIEDMAN, AND J. H. FRIEDMAN, *The elements of statistical learning: data mining, inference, and prediction*, vol. 2, Springer, 2009.

[32] Y. HE, G. REINERT, D. WIPF, AND M. CUCURINGU, *Robust Angular Synchronization via Directed Graph Neural Networks*, arXiv preprint arXiv:2310.05842, (2023).

[33] H. HOPPE, T. DEROSÉ, T. DUCHAMP, J. MCDONALD, AND W. STUETZLE, *Surface Reconstruction from Unorganized Points*, in SIGGRAPH, 1992.

[34] J. B. HOUGH, M. KRISHNAPUR, Y. PERES, AND B. VIRÁG, *Determinantal Processes and Independence*, Probability Surveys, 3 (2006), pp. 206–229.

[35] N. JANCO AND T. BENDORY, *Unrolled algorithms for group synchronization*, IEEE Open Journal of Signal Processing, (2023).

[36] H. JAQUARD, M. FANUEL, P.-O. AMBLARD, R. BARDENET, S. BARTHELMÉ, AND N. TREMBLAY, *Smoothing Complex-Valued Signals on Graphs with Monte-Carlo*, in ICASSP, 2023.

[37] B. KARRER AND M. E. NEWMAN, *Stochastic blockmodels and community structure in networks*, Physical review E, 83 (2011), p. 016107.

[38] A. KASSEL, *Learning about critical phenomena from scribbles and sandpiles*, ESAIM: Proceedings and Surveys, 51 (2015), pp. 60–73.

[39] A. KASSEL AND R. KENYON, *Random curves on surfaces induced from the Laplacian determinant*, The Annals of Probability, 45 (2017), pp. 932–964.

[40] A. KASSEL AND T. LÉVY, *Covariant symanzik identities*, Probability and Mathematical Physics, (2021).

[41] R. KENYON, *Spanning forests and the vector bundle japlacian*, The Ann. of Probability, (2011).

[42] D. P. KROESE, T. TAIMRE, AND Z. I. BOTEV, *Handbook of monte carlo methods*, 2013.

[43] A. B. KUIJLAARS, *Convergence analysis of Krylov subspace iterations with methods from potential theory*, SIAM review, 48 (2006), pp. 3–40.

[44] A. KULESZA AND B. TASKAR, *Determinantal Point Processes for Machine Learning*, Found. Trends Mach. Learn., 5 (2012), pp. 123–286.

[45] R. KYNG AND Z. SONG, *A matrix chernoff bound for strongly rayleigh distributions and spectral sparsifiers from a few random spanning trees*, in FOCS, IEEE, 2018.

[46] Y. LECUN, L. BOTTOU, Y. BENGIO, AND P. HAFFNER, *Gradient-based learning applied to document recognition*, Proceedings of the IEEE, 86 (1998), pp. 2278–2324.

[47] G. LERMAN AND Y. SHI, *Robust Group Synchronization via Cycle-Edge Message Passing*, Foundations of Computational Mathematics, 22 (2022), pp. 1665–1741.

[48] J. LESKOVEC AND A. KREVL, *SNAP Datasets: Stanford Large Network Dataset Collection*, 2014.

[49] H. LIU, X. LI, AND A. M.-C. SO, *ReSync: Riemannian Subgradient-based Robust Rotation Synchronization*, arXiv preprint arXiv:2305.15136, (2023).

[50] P. LORENZO, S. BARBAROSSA, AND P. BANELLI, *Sampling and recovery of graph signals*, in Cooperative and Graph Signal Processing, Elsevier, 2018, pp. 261–282.

[51] O. MACCHI, *The coincidence approach to stochastic point processes*, Adv. in Applied Probability, (1975).

[52] P. MARTINSSON AND J. A. TROPP, *Randomized numerical linear algebra: Foundations and algorithms*, Acta Numer., 29 (2020), pp. 403–572.

[53] T. MAUNU AND G. LERMAN, *Depth descent synchronization in SO(d)*, Int. journ. of comp. vision, (2023).

[54] M. MEZARD AND A. MONTANARI, *Information, physics, and computation*, Oxford University Press, 2009.

[55] L. A. S. MOREIRA, A. L. L. RAMOS, M. L. R. DE CAMPOS, J. A. APOLINÁRIO, AND F. G. SERRENHO, *A Graph Signal Processing Approach to Direction of Arrival Estimation*, in EUSIPCO 2019.

[56] M. E. MULLER, *Some continuous Monte Carlo methods for the Dirichlet problem*, The Annals of Mathematical Statistics, (1956), pp. 569–589.

[57] R. MURRAY, M. W. MAHONEY, N. B. ERICHSON, M. MELNICHENKO, O. A. MALIK, L. GRIGORI, P. LUSZCZEK, M. DEREZIŃSKI, M. E. LOPES, ET AL., *Randomized numerical linear algebra: A perspective on the field with an eye to software*, arXiv:2302.11474, (2023).
fortieth annual ACM symposium on Theory of computing, 2008, pp. 563–568.

[66] L. N. TREFETHEN AND D. BAU, *Numerical linear algebra*, SIAM, 2022.

[67] N. TREMBLAY, P. GONÇALVES, AND P. BORGNAT, *Design of graph filters and filterbanks*, in Cooperative and Graph Signal Processing, Elsevier, 2018, pp. 299–324.

[68] U. VON LUXBURG, *A tutorial on spectral clustering*, Stat. Comput., 17 (2007), pp. 395–416.

[69] D. B. WILSON, *Generating Random Spanning Trees More Quickly than the Cover Time*, in Proceedings of the Twenty-Eighth Annual ACM Symposium on the Theory of Computing, Philadelphia, Pennsylvania, USA, May 22-24, 1996, G. L. Miller, ed., ACM, 1996, pp. 296–303.

[70] S. X. YU, *Angular embedding: From jarring intensity differences to perceived luminance*, in 2009 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2009), 20-25 June 2009, Miami, Florida, USA, IEEE Computer Society, 2009, pp. 2302–2309.

[71] ——, *Angular Embedding: A Robust Quadratic Criterion*, IEEE Trans. Pattern Anal. Mach. Intell., 34 (2012), pp. 158–173.

[72] S. ZHANG AND Y. HUANG, *Complex quadratic optimization and semidefinite programming*, SIAM Journal on Optimization, 16 (2006), pp. 871–890.

[73] X. ZHANG, Y. HE, N. BRUGNONE, M. PERLMUTTER, AND M. J. HIRN, *MagNet: A Neural Network for Directed Graphs*, in Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, M. Ranzato, A. Beygelzimer, Y. N. Dauphin, P. Liang, and J. W. Vaughan, eds., 2021, pp. 27003–27015.

[74] Z. ZHAO AND A. SINGER, *Rotationally invariant image representation for viewing direction classification in cryo-EM*, Journal of structural biology, 186 (2014), pp. 153–166.
[58] A. B. OWEN, *Monte carlo theory, methods and examples*, 2013.

[59] N. PARIKH, S. BOYD, ET AL., *Proximal algorithms*, Foundations and trends® in Optimization, (2014).

[60] A. PERRY, A. S. WEIN, A. S. BANDEIRA, AND A. MOITRA, *Message-Passing Algorithms for Synchro-nization Problems over Compact Groups*, Com. on Pure and Applied Mathematics, (2018).

[61] Y. Y. PILAVCI, P. AMBLARD, S. BARTHELMÉ, AND N. TREMBLAY, *Graph Tikhonov Regularization and Interpolation Via Random Spanning Forests*, IEEE Trans. Signal Inf. Process. over Networks.

[62] Y. Y. PILAVCI, P.-O. AMBLARD, S. BARTHELME, AND N. TREMBLAY, *Variance reduction for inverse trace estimation via random spanning forests*, arXiv preprint arXiv:2206.07421, (2022).

[63] Y. Y. PILAVCI, P.-O. AMBLARD, S. BARTHELMÉ, AND N. TREMBLAY, *Variance Reduction in Stochastic Methods for Large-Scale Regularized Least-Squares Problems*, in EUSIPCO, IEEE, 2022.

[64] G. PUY, N. TREMBLAY, R. GRIBONVAL, AND P. VANDERGHEYNST, *Random sampling of bandlimited signals on graphs*, Applied and Computational Harmonic Analysis, 44 (2018), pp. 446–475.

[65] T. C. RAIA, M. G. P. THOMAS, F. G. SERRENHO, AND J. A. APOLINÁRIO JR, *Gsp-based doa estimation for a multimission radar*, SBrT2020, (2020).

[66] C. R. RAO, *Information and the accuracy attainable in the estimation of statistical parameters*, in Break-throughs in Statistics: Foundations and basic theory, Springer, 1992, pp. 235–247.

[67] Y. SAAD, *Iterative methods for sparse linear systems*, SIAM, 2003.

[68] ——, *Numerical methods for large eigenvalue problems: revised edition*, SIAM, 2011.

[69] R. SAWHNEY AND K. CRANE, *Monte Carlo geometry processing: A grid-free approach to PDE-based methods on volumetric domains*, ACM Transactions on Graphics, 39 (2020).

[70] S. H. SCHERES, *RELION: implementation of a Bayesian approach to cryo-EM structure determination*, Journal of structural biology, 180 (2012), pp. 519–530.

[71] SHARP, NICHOLAS AND SOLIMAN, YOUSUF AND CRANE, KEENAN, *The vector heat method*, TOG, (2019).

[72] J. SHARPNACK, A. SINGH, AND A. RINALDO, *Sparsistency of the edge lasso over graphs*, in Artificial Intelligence and Statistics, PMLR, 2012, pp. 1028–1036.

[73] L. A. SHEPP AND B. F. LOGAN, *The Fourier reconstruction of a head section*, IEEE Transactions on nuclear science, 21 (1974), pp. 21–43.

[74] D. I. SHUMAN, S. K. NARANG, P. FROSSARD, A. ORTEGA, AND P. VANDERGHEYNST, *The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains*, IEEE signal processing magazine, 30 (2013), pp. 83–98.

[75] A. SINGER, *Angular synchronization by eigenvectors and semidefinite programming*, Applied and computational harmonic analysis, 30 (2011), pp. 20–36.

[76] A. SINGER AND H.-T. WU, *Vector diffusion maps and the connection Laplacian*, Communications on pure and applied mathematics, 65 (2012), pp. 1067–1144.

[77] D. A. SPIELMAN AND N. SRIVASTAVA, *Graph sparsification by effective resistances*, in Proceedings of the fortieth annual ACM symposium on Theory of computing, 2008, pp. 563–568.

[78] L. N. TREFETHEN AND D. BAU, *Numerical linear algebra*, SIAM, 2022.

[79] N. TREMBLAY, P. GONÇALVES, AND P. BORGNAT, *Design of graph filters and filterbanks*, in Cooperative and Graph Signal Processing, Elsevier, 2018, pp. 299–324.

[80] J. A. TROPP, *User-friendly tail bounds for sums of random matrices*, Found. of comput. maths., (2012).

[81] N. K. VISHNOI ET AL., $Lx=b$, Foundations and Trends® in Theoretical Computer Science, (2013).

[82] D. B. WILSON, *Generating Random Spanning Trees More Quickly than the Cover Time*, in STOC, 1996.

[83] S. X. YU, *Angular embedding: From jarring intensity differences to perceived luminance*, in CVPR, 2009.

[84] ——, *Angular Embedding: A Robust Quadratic Criterion*, Trans. Pattern Anal. Mach. Intell., (2012).

[85] C. ZHAN, G. CHEN, AND L. F. YEUNG, *On the distributions of laplacian eigenvalues versus node degrees in complex networks*, Physica A: Statistical Mechanics and its Applications, 389 (2010), pp. 1779–1788.

[86] S. ZHANG AND Y. HUANG, *Complex quadratic optimization and semidefinite programming*, SIAM Journal on Optimization, 16 (2006), pp. 871–890.

[87] X. ZHANG, Y. HE, N. BRUGNONE, M. PERLMUTTER, AND M. J. HIRN, *MagNet: A Neural Network for Directed Graphs*, in NeurIPS, 2021.

[88] Z. ZHAO AND A. SINGER, *Rotationally invariant image representation for viewing direction classification in cryo-EM*, Journal of structural biology, 186 (2014), pp. 153–166.
# Supplementary Material

We work with a weighted graph $\mathcal{G}$ in all the following proofs, with $w_e \in \mathbf{R}_+$ the weight of edge $e$. We still denote by $\mathbf{L}_\theta$ the resulting (weighted) connection Laplacian, with off-diagonal entries $(\mathbf{L}_\theta)_{i,j} = -w_{i,j}e^{i\theta_{(j,i)}}$ if $\{i, j\} \in \mathcal{E}$ ($(\mathbf{L}_\theta)_{i,j} = w_{(i,j)} = 0$ otherwise), and diagonal $(\mathbf{L}_\theta)_{i,i} = d_i = \sum_j w_{i,j}$ the weighted degree of node $i$.

**6. Proof of Proposition 2.1.** Denote by $C: \mathbf{C}^\mathcal{V} \to \mathbf{R}$ the cost function

$$ (6.1) \qquad C(f) = q\|f - g\|_2^2 + \langle f, \mathbf{L}_\theta f \rangle. $$

We use a standard argument from **CR**-calculus, and seek the zeros of the Fréchet Wirtinger derivatives of $C$ (see e.g. [9]). Let us first compute these derivatives: for $f, h \in \mathbf{C}^\mathcal{V}$, we have

$$ (6.2) \qquad C(f+h) = q\langle(f+h) - g, (f+h) - g\rangle + \langle(f+h), \mathbf{L}_\theta(f+h)\rangle, $$

$$ (6.3) \qquad = C(f) + D_f^W C(h) + D_f^{W*} C(h) + (q\langle h, h\rangle + \langle h, \mathbf{L}_\theta h\rangle), $$

where $(h \mapsto q\langle h, h\rangle + \langle h, \mathbf{L}_\theta h\rangle) \in o(h)$, and the associated (conjugate) Fréchet Wirtinger derivatives of $C$ (at $f$) $D_f^W C$ and $D_f^{W*} C: \mathbf{C}^\mathcal{V} \to \mathbf{C}$ are given by:

$$ (6.4) \qquad D_f^W C(h) = \langle q(f-g) + \mathbf{L}_\theta, h \rangle $$

$$ (6.5) \qquad D_f^{W*} C(h) = \langle h, q(f-g) + \mathbf{L}_\theta \rangle. $$

The result follows from the fact that $D_f^W C = 0$ (resp. $D_f^{W*} C = 0$) identically if and only if $f = q(\mathbf{L}_\theta + q\mathbf{I})^{-1}g$.

**7. A Generalization of Property (P2).** Define $\Delta_\theta = \mathbf{D}^{-1}\mathbf{L}_\theta$ the weighted connection-aware random walk Laplacian $\Delta_\theta$. Proposition 7.1 generalizes property (P2).

**Proposition 7.1.** For all $l \ge 1$, we have

$$ (7.1) \qquad (l - \Delta_{\theta})_{i,j}^{l} = \sum_{\substack{p \in P_i^j \\ l(p)=l}} \left( \prod_{0 \le k < l} \frac{w_{(u_k, u_{k+1})}}{d_{u_k}} \right) \psi_{p^*}, $$

with $l(p)$ the length of a path $p = ((u_0, u_1), (u_1, u_2), \dots, (u_{l(p)-1}, u_{l(p)}))$.

For each path $p$, the product in the r.h.s. of Equation (7.1) is the probability of observing $p$ when performing a random walk from $i$ to $j$ on $\mathcal{G}$. For the trivial connection this expression is equal to $(\mathbf{D}^{-1}\mathbf{A})_{i,j}^l$, the probability of a path going from $i$ to $j$ in $l$ steps, but Proposition 7.1 more broadly captures *random propagations* along the corresponding sampled paths.

Observations similar to Proposition 7.1 have been laid out in a few works, such as [34] for continuous time random walk, and in a number of other situations, where the connection often stems from theoretical physics considerations [12].
# Supplementary Material

We work with a weighted graph $\mathcal{G}$ in all the following proofs, with $w_e \in \mathbf{R}_+$ the weight of edge $e$. We still denote by $\mathbf{L}_{\theta}$ the resulting (weighted) connection Laplacian, with off-diagonal entries $(\mathbf{L}_{\theta})_{i,j} = -w_{i,j}e^{i\theta_{(j,i)}}$ if $\{i, j\} \in \mathcal{E}$ ($(\mathbf{L}_{\theta})_{i,j} = w_{(i,j)} = 0$ otherwise), and diagonal $(\mathbf{L}_{\theta})_{i,i} = d_i = \sum_j w_{i,j}$ the weighted degree of node $i$.

**6. Proof of Proposition 2.1.** Denote by $C: \mathbf{C}^{\mathcal{V}} \to \mathbf{R}$ the cost function

$$ (6.1) \qquad C(f) = q\|f - g\|_2^2 + \langle f, \mathbf{L}_{\theta}f \rangle. $$

We use a standard argument from **CR**-calculus, and seek the zeros of the Fréchet Wirtinger derivatives of $C$ (see e.g. [9]). Let us first compute these derivatives: for $f, h \in \mathbf{C}^{\mathcal{V}}$, we have

$$ (6.2) \qquad C(f+h) = q\langle(f+h) - g, (f+h) - g\rangle + \langle(f+h), \mathbf{L}_{\theta}(f+h)\rangle, $$

$$ (6.3) \qquad = C(f) + D_f^W C(h) + D_f^{W*} C(h) + (q\langle h, h \rangle + \langle h, \mathbf{L}_{\theta}h \rangle), $$

where $(h \mapsto q\langle h, h \rangle + \langle h, \mathbf{L}_{\theta}h \rangle) \in o(h)$, and the associated (conjugate) Fréchet-Wirtinger derivatives of $C$ (at $f$) $D_f^W C$ and $D_f^{W*} C: \mathbf{C}^{\mathcal{V}} \to \mathbf{C}$ are given by:

$$ (6.4) \qquad D_f^W C(h) = \langle q(f-g) + \mathbf{L}_{\theta}, h \rangle $$

$$ (6.5) \qquad D_f^{W*} C(h) = \langle h, q(f-g) + \mathbf{L}_{\theta} \rangle. $$

The result follows from the fact that $D_f^W C = 0$ (resp. $D_f^{W*} C = 0$) identically if and only if $f = q(\mathbf{L}_{\theta} + q\mathbf{I})^{-1}g$.

**7. A Generalization of Property (P2).** Define $\Delta_{\theta} = \mathbf{D}^{-1}\mathbf{L}_{\theta}$ the weighted connection-aware random walk Laplacian $\Delta_{\theta}$. Proposition 7.1 generalizes property (P2).

**Proposition 7.1.** For all $l \ge 1$, we have

$$ (7.1) \qquad (l - \Delta_{\theta})_{i,j}^l = \sum_{\substack{p \in P_i^j \\ l(p)=l}} \left( \prod_{0 \le k < l} \frac{w_{(u_k, u_{k+1})}}{d_{u_k}} \right) \psi_{p^*}, $$

with $l(p)$ the length of a path $p = ((u_0, u_1), (u_1, u_2), \dots, (u_{l(p)-1}, u_{l(p)}))$.

For each path $p$, the product in the r.h.s. of Equation (7.1) is the probability of observing $p$ when performing a random walk from $i$ to $j$ on $\mathcal{G}$. For the trivial connection this expression is equal to $(\mathbf{D}^{-1}\mathbf{A})_{i,j}^l$, the probability of a path going from $i$ to $j$ in $l$ steps, but Proposition 7.1 more broadly captures *random propagations* along the corresponding sampled paths.

Observations similar to Proposition 7.1 have been laid out in a few works, such as [40] for continuous time random walk, and in a number of other situations, where the connection often stems from theoretical physics considerations [12].