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.

RandNLA, graphs, random spanning forests, connection, graph signal processing, synchronization

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 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 $\mathsf{L}_\theta$ (the formal definition of $\mathsf{L}_\theta$ is postponed to Section <ref>), 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 $\mathcal{O}(n^3)$) remains prohibitive beyond moderately large graphs ($n \simeq 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\left( \frac{\sigma}{m} \right)$ (with $\sigma^2$ the variance of the estimator and $m$ the number of Monte-Carlo samples). Methods specialized in graphs, based on random walks and decompositions of the graph, have also appeared  (e.g. [65, 39, 51]).

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:

\begin{equation}
\argmin_{f \in \C^n} \ q \Vert f - g \Vert_2^2 + \frac{1}{2} \sum_{i,j} \vert f(j) - e^{\iota \theta_{i,j}} f(i) \vert^2,
\label{eq:tikhonov_intro}
\end{equation}

where $g \in \C^n$, $q \in \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 $\mathsf{L}_\theta$. The solution to this first problem takes the form $q(\mathsf{L}_\theta + q\mathsf{I})^{-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] :

\begin{equation}
\theta_{i,j} = \omega_j - \omega_i + \eps_{i,j} \mod 2\pi,
\end{equation}

where $\eps_{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 <ref> and <ref>).

The problem becomes more involved when considering imperfect measurements $\theta_{i,j}$, with non-zero noise $\eps_{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 <ref>).

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.

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).

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.

Propagations on the $4 \times 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 <ref>: propagation along an arbitrary path. Center <ref>: exact synchronization achieved by propagation along a spanning tree, in the absence of noise. Right <ref>: synchronization impossible due to noise and incoherence along cycles.

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

\begin{equation}
\mathcal{I}(s) = \sum_{\{i,j\} \in \mathcal{E}} \left( 2 - 2 \cos{\left( (s_j - s_i) - \theta_{i,j} \right)} \right),
\label{eq:cosine_loss}
\end{equation}

before minimizing this incoherence over all possible assignments. While this problem is non-convex and $\mathrm{NP}$-hard in general [72] (see also [10] for a discussion), different techniques allowing to recover a solution have been proposed [71, 63, 10, 50, 27], 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 $\mathsf{L}_\theta$ (see Section <ref>).

§.§ 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 <ref> for a formal statement).

For MTSFs sampled from the probability distribution in Equation (<ref>), 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 (<ref>).

See Figure <ref> for an illustration. We include a teaser runtime comparison with a standard deterministic solver in Figure <ref>, which shows that our estimators are less sensitive to the density of the graph, and can provide significant speed-ups for equivalent precision.

A MTSF in purple, with two unicycles (top left, bottom) and one tree (right). The large purple node is the root of the tree.

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).

Theorem <ref> illustrated on the $4 \times 4$ grid graph. Left: a MTSF. Right: estimation by propagations along the branches of the MTSF. The estimation in $\mathbf{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. From the theoretical side, our arguments rely on the theory of Determinantal Point Processes (DPPs) [29, 44], and generalize the combinatorial analyses of [51, 35].

DC-SBM with mean degree $\overline{d} \simeq 30$.

DC-SBM with mean degree $\overline{d} \simeq 420$.

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 ($\mathsf{D}^{-1}$+CG) on degree-corrected stochastic block model graphs (DC-SBM, see Section <ref> 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 $\overline{d}$: on the left, a DC-SBM with average degree $\overline{d} \simeq 30$ and, on the right, a DC-SBM with average degree $\overline{d} \simeq 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 <ref> for details).

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 number of different manners, some of which we list below.

Theoretical results
In addition to Theorem <ref>, we derive a connection-aware Feynman-Kac formula (Proposition <ref>), 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 <ref> 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 [71, 63, 10].
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 <ref>.
The code used for these experiments is publically available[<https://gricad-gitlab.univ-grenoble-alpes.fr/gaia/synchromtsf>].

§.§ Related Work

Our main result can be understood as a specialization of the approach in [20, 19], 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.

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 <ref> of the Supplementary Material).

§.§ Organisation of the Paper

In Section <ref>, we introduce preliminary background on graphs and connections used in the rest of the paper. In Section <ref>, we present our Feynman-Kac formula (Section <ref>) and our MTSF-based estimators (Section <ref>), analyze an efficient sampling algorithm for MTSFs (Section <ref>), and describe our variance-reduction techniques (Section <ref>). We analyze the numerical behavior of our estimators in Section <ref>, and compare them with standard deterministic algorithms (Sections <ref> and <ref>). In Section <ref>, we propose an iterative scheme leveraging our smoothing estimators to solve the angular synchronization problem (Section <ref>), illustrate its application to a denoising problem inspired from cryogenic electron microscopy (Section <ref>), and compare this scheme with standard deterministic methods (Section <ref>). Proofs are deferred to the Supplementary Material.

§ 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 $\overrightarrow{\mathcal{E}}$. The size of $\overrightarrow{\mathcal{E}}$ is twice the size of $\mathcal{E}$:  each edge $e \in \mathcal{E}$ is associated to two oriented edges in $\overrightarrow{\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 \overrightarrow{\mathcal{E}}$, and its reversely-oriented edge $e^* = (t_e,s_e) \in \overrightarrow{\mathcal{E}}$. We state our results for unweighted graphs, the weighted case is included in the Supplementary Material.

§.§ 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~\R^\mathcal{V}$, and called a graph signal [62]. The regularity of such a signal can be quantified using the combinatorial Laplacian $\mathsf{L}~:~\R^\mathcal{V}~\rightarrow~\R^\mathcal{V}$, a symmetric semi-definite positive operator, conveniently expressed as $\mathsf{L} = \mathsf{D} - \mathsf{A}$, where $\mathsf{D}$ is the diagonal degree matrix ($\mathsf{D}_{i,i} = d_i$ the degree of node $i$) and $\mathsf{A}$ the adjacency matrix of the graph (with $\mathsf{A}_{(i,j)} = 1$ if $(i,j) \in \overrightarrow{\mathcal{E}}$ and $0$ otherwise) [15]. The quadratic form associated to $\mathsf{L}$ acts as:

\begin{equation}
\langle f, \mathsf{L} f \rangle = \frac{1}{2} \sum_{e \in \overrightarrow{\mathcal{E}}} \vert f(t_e) - f(s_e) \vert^2,
\label{eq:action_laplacian}
\end{equation}

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:

* $\ker{\mathsf{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].

* $\mathsf{D}^{-1} \mathsf{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 $(\mathsf{D}^{-1} \mathsf{A})^l_{i,j}$ the probability that a path of length $l$ starting at $i$ ends up in node $j$.

§.§ Connection Laplacian: Definition and Basic Properties

Unitary connections are a practical way to add additional information to the graph's structure[By describing explicit geometrical transformations between the signal values of the nodes of $\mathcal{G}$.]. 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 \C^\mathcal{V}$ when studying variational properties on $\mathcal{G}$. The entries of $f$ should be understood as belonging to different copies $\C_v$ of the complex plane $\C$ associated to each node $v \in \mathcal{V}$, see Figure <ref>[This association of a copy of $\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}$.].

We will associate to each directed edge in $\overrightarrow{\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 \overrightarrow{\mathcal{E}}}$, with each map $\psi_e~:~\C_{s_e}\rightarrow\C_{t_e}$ acting as multiplication by a unitary complex number $e^{\iota \theta_e}$, so that $\psi_e(z) = e^{\iota \theta_e} \cdot z$, with $\iota$ the complex imaginary unit (we sometimes abuse notation and write $\psi_e = e^{\iota \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 <ref>.

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).

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 \R_+$, resulting in connections such that $\psi_e (z) = e^{\iota \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 $\mathsf{L}_\theta:~\C^\mathcal{V}\rightarrow~\C^\mathcal{V}$, defined  by $\mathsf{L}_\theta = \mathsf{D} - \mathsf{A}_\theta$, with $\mathsf{A}_\theta$ a connection-aware adjacency matrix such that $(\mathsf{A}_\theta)_{i,j} = e^{-\iota \theta_{(i,j)}}$ if $\{i,j\} \in \mathcal{E}$, and $\mathsf{A}_{i,j} = 0$ otherwise. This operator generalizes the Laplacian to non-trivial connections, and is associated to the quadratic form:

\begin{equation}
\langle f, \mathsf{L}_\theta f \rangle = \frac{1}{2} \sum_{e \in \overrightarrow{\mathcal{E}}} \vert f(t_e) - e^{\iota \theta_{e}} f(s_e) \vert^2.
\label{eq:connection_laplacian_action}
\end{equation}

For the trivial connection, with $\psi_e = \id_{\C_{s_e},\C_{t_e}}$ (i.e. $\theta_e = 0$) for all edges $e \in \overrightarrow{\mathcal{E}}$, we recover $\mathsf{L}_\theta = \mathsf{L}$ and the expression in Equation (<ref>). Unlike this specific case, $\langle f, \mathsf{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 \mathsf{L}_\theta$ non-empty, and what are the functions $f \in \C^\mathcal{V}$ that are not penalized by $\langle f , \mathsf{L}_\theta f \rangle$?

($\mathrm{P1'}$) Angular synchronization and $\ker{L_\theta}$
When the entries $f_i$ of the vector $f$ are unitary complex numbers $f_i = e^{\iota s_i}$, we recover from Equation (<ref>) the expression in Equation (<ref>).
In this case, $\vert f(j) - e^{\iota \theta_{(i,j)}} f(i) \vert^2 = 2 - 2 \mathfrak{Re}((e^{\iota s_{j}})^* e^{\iota \theta_{(i,j)}} e^{\iota s_{i}})$, which, thanks to conjugation-invariance of the real part, results in:

\begin{equation}
\langle f, \mathsf{L}_\theta f \rangle = \sum_{\{i,j\} \in E} 2 - 2 \cos((s_j - s_i) - \theta_{(i,j)}).
\end{equation}

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

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

§ 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 \C^\mathcal{V}$ by solving:

\begin{equation}
\argmin_{f \in \C^\mathcal{V}} \ q \Vert f - g \Vert_2^2 + \langle f, \mathsf{L}_\theta f \rangle,
\label{eq:tikhonov}
\end{equation}

with $q \in \R_+^*$. This is nothing but a re-writing of Equation (<ref>), and $\langle f , \mathsf{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.

The solution to Problem (<ref>) can be expressed as:

\begin{equation}
f_* = q(\mathsf{L}_\theta + q \mathsf{I})^{-1} g.
\label{eq:solution_tikhonov}
\end{equation}

The proof of Proposition <ref>, based on straightforward $\mathbf{CR}$-calculus, is detailed in Section <ref> 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}$.

§.§ 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 (<ref>). We start with a few definitions.
A path $p$ in $\mathcal{G}$ is an ordered sequence of oriented edges in $\overrightarrow{\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} \rightarrow \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 \geq 0}$, which starts at some fixed node $u_0 = i \in 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{\mathsf{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),...,(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 (<ref>) on node $i$. More precisely:

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

\begin{equation}
f_*(i) = \E_{p \sim \nu_i}\left( \psi_{p^*_\Gamma}(g_j) \right),
\label{eq:feynmac_kac}
\end{equation}

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 <ref>), taking into account the rotations along the path. See Algorithm <ref>. This formula is a discrete-time analog of a Theorem in [34], and we prove a more general version in the Supplementary Material (Section <ref>, Proposition <ref>), building upon a connection-aware version of property <ref>.
Feynman-Kac estimator (Proposition <ref>).
[1]
$f_i \leftarrow 0$
Repeat $m$ times
Sample a path $p \in P^\Gamma_i$ from $\nu_i$ By running an interrupted random walk on $\mathcal{G}$
$f_i \leftarrow f_i + \psi_{p^*_\Gamma}(g_j)$ $j$ the last node before interruption of $p$
Output $\frac{1}{m} f_i$

Algorithm <ref> 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 $\vert \mathcal{V} \vert$ 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 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 $\vert \mathcal{E} \vert$.

The restriction to complex signals on $\mathcal{G}$ and unitary connections (encoding 2D rotations) is mostly artificial, and Proposition <ref> generalizes to transformations in e.g. $O(\R^d)$. This setting is encountered for instance for 3D molecule alignement 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(\R^3))$.

§.§ 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 $\vert \phi \vert = \vert \mathcal{V} \vert$, 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 <ref>.
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}_\mathcal{M}$ over $\mathcal{M}(\mathcal{G})$ defined by:

\begin{equation}
\prb_{\mathcal{D}_\mathcal{M}}(\phi) \propto q^{\vert \phi \cap \mathcal{V} \vert} \prod_{C \in \mathcal{C}(\phi)} \left( 2 - 2 \cos(\theta_C) \right),
\label{eq:prob_dm}
\end{equation}

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 <ref>). 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}_\mathcal{M}$ is a rooted variant of the distribution considered in [22], and is in fact a DPP[So is the distribution of [22], but their non-combinatorial argument does not carry over to our rooted process. See the Section <ref> of the Supplementary Material for a proof in this case].
We can now state our main theorem.

Let $\phi$ be a MTSF sampled according to $\mathcal{D}_\mathcal{M}$. Denote by $a \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
\begin{equation}
\tilde{f}(i,\phi,g) = \psi_{r_\phi(i) \xrightarrow[]{\phi} i} \left( g(r_\phi(i)) \right)
\end{equation}

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:

\begin{equation}
f_*(i) = \mathbf{E}_{\mathcal{D}_\mathcal{M}}(\tilde{f}(i,\phi,g)).
\label{eq:fermionic}
\end{equation}

The proof of Theorem <ref> relies on the reformulation of $\mathcal{D}_\mathcal{M}$ as a DPP, and an extension of arguments coming from both [35] and [51]. It is a special case of Theorem <ref>, proved in Section <ref> 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 <ref> and Figure <ref>.

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

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

Theorem <ref> illustrated on the $4 \times 4$ grid graph, on one MTSF.

MTSF-based estimator (Theorem <ref>).
[1]
$f_i \leftarrow 0 \ \forall i \in \mathcal{V}$
Repeat $m$ times
Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_\mathcal{M}$ Algorithm <ref> can be used
$i \in \mathcal{V}$
$i$ belongs to a rooted tree of $\phi$
$f_i \leftarrow f_i + \psi_{r_\phi(i) \xrightarrow[]{\phi} i} \left( g(r_\phi(i)) \right)$
Output $\frac{1}{m} (f_1,...,f_{\vert \mathcal{V} \vert})$

If the connection $\Psi$ is trivial (i.e., $\psi_e = \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}_\mathcal{M}$ reduces to another determinantal distribution $\mathcal{D}_\mathcal{F}$ over rooted spanning forests $\phi \in \mathcal{V} \cup \mathcal{E}$ [3, 4].

Unlike Algorithm <ref>, Algorithm <ref>  allows to update the estimates on all the nodes for each sampled MTSF, in $\mathcal{O}(\vert \mathcal{V} \vert)$ time. On the other hand, it becomes necessary to sample a MTSF $\phi$ according to $\mathcal{D}_\mathcal{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 $(\vert V \vert + \vert E \vert) \times (\vert V \vert + \vert E \vert)$), 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.

§.§ Sampling Multi-Type Spanning Forests

Efficient sampling from $\mathcal{D}_\mathcal{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[Condition <ref> 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.].
For all cycles $C$ in $\mathcal{G}$, $\cos(\theta_C) \geq 0$.

We stress that this condition must hold for all the cycles of the graph, and not just for the cycles $\mathcal{C}(\phi)$ in some MTSF $\phi \in \mathcal{M}(G)$.
When Condition <ref> is satisfied, $1 -\cos(\theta_C) \leq 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.

Even though efficient sampling from $\mathcal{D}_\mathcal{M}$ can only be performed under Condition <ref>, Theorem <ref> holds beyond this technical condition. Furthermore, in some applicative settings (e.g. directed graph signal processing [24], ranking [71, 16] 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 <ref>, samples a MTSF $\phi$ according to $\mathcal{D}_\mathcal{M}$ by simulating multiple random walks on $\mathcal{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 <ref>, 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$[Some implementation details of Algorithm <ref> 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 <ref> of the Supplementary Material and our Julia implementation for more details]. We suppose that the nodes of $\mathcal{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 $\mathrm{random\_successor}(u)$ that, at each call, randomly outputs either $\Gamma$ with probability $\frac{q}{d_u + q}$, or some node $v$ with probability  $\frac{\mathsf{A}_{u,v}}{d_u + q}$.

MTSF sampling algorithm [22].
[1]
$\phi \leftarrow \emptyset$
$\phi$ not spanning
Let $i \in \mathcal{V}$ be the first node in the queue not spanned by $\phi$
$u \leftarrow i$ $u$ is the current node of the random walk
$p \leftarrow \epsilon$
$\epsilon$ the empty path
($u \neq \Gamma$) and ($p$ does not intersect $\phi$ or contain a cycle)
$u' \leftarrow \mathrm{random\_successor}(u)$ Move to next node
$u' \neq \Gamma$
$e \leftarrow (u,u')$, $p \leftarrow p e$ Add $e$ to the path $p$
${p}$ contains a cycle $C$
Remove this cycle from $p$ with probability $\cos(\theta_C)$
$u_\mathrm{old} \leftarrow u$, $u \leftarrow u'$ $u_\mathrm{old}$ the previous node
$u = \Gamma$
$\phi \leftarrow \phi \cup p \cup u_\mathrm{old}$ Add the sampled path $p$ and the root $u_\mathrm{old}$ to $\phi$
$\phi \leftarrow \phi \cup p$
Output $\phi$

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

Suppose that Condition <ref> holds. Then, Algorithm <ref> outputs a MTSF $\phi$ distributed according to $\mathcal{D}_\mathcal{M}$.

Let us now discuss the cost of Algorithm <ref>.

The expected number of random walk steps of Algorithm <ref> is bounded from above by

\begin{equation}
\mathrm{tr}\left( (\mathsf{L} + q \mathsf{I})^{-1} (\mathsf{D} + q \mathsf{I}) \right).
\label{eq:complexity}
\end{equation}

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

See Section <ref> for the proof. The bound is obtained by noting that Algorithm <ref> 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 $\eqref{eq:complexity}$.
Unlike the algorithm in [51] though, Algorithm <ref> 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 <ref> of the Supplementary Material, one of them resulting in the $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q} \right)$ expected time complexity mentioned in Proposition <ref>. This translates to $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q} + \vert \mathcal{V} \vert \right)$ expected runtime for  Algorithm <ref>.We abuse notations and denote here by $p$ the set of non-oriented edges in the path.

In case Condition <ref> is not satisfied, one can still estimate $f_*$ up to a multiplicative constant ($a f_*$ for some $a$), using the importance sampling strategy from [22]. Specifically, one can threshold the incoherence $2 - 2\cos(\theta_C)$ and sample from

\begin{equation}
\prb_{IS}(\phi) \propto q^{\vert \phi \cap \mathcal{V} \vert} \prod_{C \in \mathcal{C}(\phi)} \min(2,2 - 2 \cos(\theta_C))
\label{eq:importance}
\end{equation}

using a straightforward variant of Algorithm <ref>. The estimation then uses the importance weights

\begin{equation}
w(\phi) = \prod_{C \in \mathcal{C}(\phi)} \max(1,1-\cos(\theta_C)).
\end{equation}

Both the Feynman-Kac-based Algorithm <ref> and the MTSF-based Algorithm <ref> 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.

§.§ 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 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 <ref> 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:

\begin{equation}
\overline{f}(i,\phi,g) = \psi_{r_\phi(i) \xrightarrow[]{\phi} i} \left( h_\phi(r_\phi(i),g) \right) \text{ if $i$ belongs to a rooted tree,}
\label{eq:rao}
\end{equation}

where $h_\phi$ is defined as

\begin{equation}
h_\phi(r,g) = \frac{\sum_{j \in c_\phi(r)} \psi_{j \xrightarrow[]{\phi} r} g(j)}{\vert c_\phi(r) \vert},
\end{equation}
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 $\overline{f}(i,\phi,g) = 0$.

This results in Algorithm <ref>, 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 <ref>, this procedure can be implemented at little additional cost[Note also the similarity to the Belief Propagation algorithm over trees [47].].

Rao-Blackwellized MTSF-based estimator (Theorem <ref>).
[1]
$f_i = 0 \ \forall i \in \mathcal{V}$, $h_r = 0 \ \forall r \in \mathcal{V}$
Repeat $m$ times
Sample $\phi \in \mathcal{M}(\mathcal{G})$ from $\mathcal{D}_\mathcal{M}$ Via Algorithm <ref>
$r \in \phi \cap \mathcal{V}$
$j \in c_\phi(r)$
$h_r = h_r + \psi_{j \xrightarrow[]{\phi} r} g(j)$ Propagate and average at the root
$i \in \mathcal{V}$ belonging to a rooted tree of $\phi$
$f_i = f_i + \psi_{r_\phi(i) \xrightarrow[]{\phi} i} \left( \frac{h_{r_\phi(i)}}{\vert c_\phi(r_\phi(i)) \vert} \right)$ Propagate back
Output $\frac{1}{m} (f_1,...,f_{\vert \mathcal{V} \vert})$

Control variates
Second is the introduction of control variates: an addition of another random quantity to $\overline{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$:

\begin{equation}
\hat{f}(\phi,g) = \overline{f}(\phi,g) - \alpha \mathsf{P} (q^{-1} (\mathsf{L}_\theta + q \mathsf{I}) \overline{f}(\phi,g) - g),
\label{eq:gradient_step}
\end{equation}

where $\mathsf{P} = \left( q^{-1} \mathsf{D} + \mathsf{I} \right)^{-1}$ is a diagonal preconditioner for the system $q^{-1} (\mathsf{L}_\theta + q \mathsf{I}) f = g$, and $\alpha \in \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 <ref> of the Supplementary Material for empirical results backing this choice).

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

$\mathbf{E}_{\mathcal{D}_\mathcal{M}}(\overline{f}(i,\phi,g)) = \mathbf{E}_{\mathcal{D}_\mathcal{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 $\overline{f}$ has $\mathcal{O}(\vert \mathcal{V} \vert)$ additional cost, and the gradient-descent step in $\hat{f}$ entails a single matrix-vector multiplication, in $\mathcal{O}(\vert \mathcal{E} \vert )$ time.

§ NUMERICAL RESULTS UNDER WEAK-INCONSISTENCY

We now analyze the behavior of the estimator proposed in Theorem <ref>, along with the improved versions discussed in Section <ref>, and compare their performance with (deterministic)  conjugate-gradient-based solvers [56], on graphs with different topologies[In this Section and in Section <ref>, we perform all our measurements using a single thread on a laptop with an intel i7-1185G7 processor.]. 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

\begin{equation}
\theta_e = \omega_{t_e} - \omega_{s_e} + \eta \eps_e
\label{eq:degradation}
\end{equation}

for all edges $e$, with $\eps_e$ a perturbation uniformly distributed in $[-1,1]$, and $\eta \in \R_+^*$ a scaling constant. We set  $\eta = \frac{\pi}{2n}$ to ensure that Condition <ref> is satisfied (so that  Algorithm <ref> in fact samples a MTSF according to $\mathcal{D}_\mathcal{M}$).

§.§ 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 $\overline{d}$ of the graph (controlling the graph's density). In our first experiment, we let $\overline{d}$ take values in $\{50,100,150,200\}$, and $q$ take values equal to $q' \overline{d}$, with $q' \in \{10^{-3},10^{-1},1\}$. The complexity bound of Equation (<ref>) becomes linear in $\vert \mathcal{V} \vert$ when using such a parametrization.

We generate Erdös-Rényi random graphs $\mathcal{G} \sim \mathrm{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{\overline{d}}{n-1}$ so that the expected mean degree of these random graphs is $\overline{d}$.
For each such graph, we generate a random signal $f \in \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 $\mathsf{L}_\theta f$, where $\mathsf{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 <ref>, are the average over $10000$ time measurements, themselves averaged over $10$ realizations of the graph $G$.

Runtime when varying $\overline{d}$, for different values of $q$.

Two trends emerge in Figure <ref>. 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 $\mathsf{L}_\theta f$, our estimator is less sensitive to the density of the graph: computing $\tilde{f}$ is more expensive than computing $\mathsf{L}_\theta f$ for $\overline{d}=50$, and systematically faster for $\overline{d} = 200$. This reflects the $\mathcal{O}(\vert \mathcal{V} \vert)$ expected complexity for $\tilde{f}$ (obtained with the parametrization $q = q' \overline{d}$).

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

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

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

In the specific instance where $f_\top \in \R^\mathcal{V}$ and the connection on $\mathcal{G}$ is trivial ($\psi_e = \mathrm{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 $\mathsf{L}$ (associated to the smallest eigenvalues) [53, 43].

We here similarly assume that $f_\top \in \C^\mathcal{V}$ is $B$-bandlimited (i.e. $f_\top \in \mathrm{ran}(u_1,...,u_B)$[The linear subspace generated by the vectors $u_1$,...,$u_B$.], with $u_i$ the $i$-th eigenvector of $\mathsf{L}_\theta$), so that solving the Tikhonov smoothing Problem (<ref>):

\begin{equation}
\argmin_{f \in \C^\mathcal{V}} \ q \Vert f - g \Vert_2^2 + \langle f, \mathsf{L}_\theta f \rangle,
\end{equation}

should allow to faithfully recover $f_\top$ from $g$, by penalizing high-frequency components in the optimal solution $f_* = q (\mathsf{L}_\theta + q \mathsf{I})^{-1} g$. We refer the reader to Section <ref> 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) = \Vert f - f_\top \Vert_2 / n$, measuring the quality of the denoising.

* The approximation error $e_a(f) = \Vert f - f_* \Vert_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$.

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 $\eps \sim \mathcal{N}_\C(0,1)$ (independently on each entry), and determine the optimal parameter $q_*$ for which $e_r(f_*) = e_r( q_* (\mathsf{L}_\theta + q_* \mathsf{I})^{-1} g)$ is minimized[We perform our search in $(0,30)$.], 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 $\overline{f}$ and $\hat{f}$, the conjugate-gradient descent with no preconditioner, with a simple diagonal preconditioner $\mathsf{P} = ( q^{-1} \mathsf{D} + \mathsf{I})^{-1}$, and with a CROUT ILU preconditioner[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$.]. 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 $\eps$ and $10$ samples for each random graph model. We plot in Figure <ref> the mean results.

Graphs used
We use the following graphs:

* An $\eps$-graph obtained by sampling $\vert \mathcal{V} \vert = 10000$ i.i.d. points $x_i$'s in $[0,1]^3$, with an edge between $x_i$ and $x_j$ whenever $\Vert x_i - x_j \Vert_2 < 0.1$.

* A graph generated from a Stochastic Block Model (SBM), with $\vert \mathcal{V} \vert = 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 $\overline{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 $\mathbf{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[Specifically, the connectivity parameters are independently sampled from a mixture of $\mathcal{N}(50,20)$, $\mathcal{N}(500,100)$ and $\mathcal{N}(10000,100)$, with weights of $0.59$, $0.4$ and $0.01$ respectively, and then normalized.], resulting in a graph with mean degree $\overline{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[Here from the mixture of $\mathcal{N}(50,20)$, $\mathcal{N}(1000,50)$, $\mathcal{N}(5000,100)$ and $\mathcal{N}(10000,100)$, with weights of $0.45$, $0.1$, $0.44$ and $0.01$.]. 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][We use the library available at <https://github.com/JuliaGraphs/SNAPDatasets.jl>.].

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 <ref> the information concerning the graphs generated, and the associated optimal $q_*$'s (averaged over all realizations of the noise $\eps)$. All those graphs are endowed with a synthetic connection as specified in Equation (<ref>). The parameters used for the (DC) SBM models ensures that the resulting graphs have a strong community structure.

Graph     $\vert \mathcal{V} \vert$     $\overline{d}$     $d_\mathrm{min}$     $d_\mathrm{max}$     $q_*$

$\eps$-graph     $10000$     $37.3$     $5.2$     $67.5$     $6.518$

SBM     $10000$     $39.87$     $17.7$     $66.6$     $21.04$

DC-SBM $1$     $9833.7$     $33.4$     $1$     $906.1$     $2.634$

DC-SBM $2$     $9932.8$     $418.8$     $1$     $1513.7$     $4.31$

AS CAIDA     $26475$     $4.032$     $1$     $2628$     $0.4808$

Experimental parameters associated to each graph. Values are averaged over all $10$ realizations for random graph models.

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

$\eps$-graph. Left: Approximation error. Right: Reconstruction error.

SBM Reconstruction error.

DC-SBM $1$ Reconstruction error.

DC-SBM $2$ Reconstruction error.

AS CAIDA Reconstruction error.

Runtime-precision trade-offs for Tikhonov smoothing with optimal $q$. Each datapoint 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.

Let us first comment on the approximation error for the $\eps$-graph (Figure <ref>, 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 <ref>). Overall, we observe that:

* $\hat{f}$ consistently performs better than $\overline{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).

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

* CG solvers perform worse than our methods when the density increases.

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

\begin{equation}
\kappa = \frac{q_* + \lambda_n}{q_* + \lambda_1} \leq \frac{q_* + 2 \max_{v \in V} d_v }{q_* + \lambda_1},
\end{equation}

where $\lambda_1 \leq ... \leq \lambda_n$ are the eigenvalues of $\mathsf{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 ($\mathcal{O}(\vert \mathcal{E} \vert)$ time per iteration), and $\overline{f}$ and $\hat{f}$ reach near-optimal reconstruction around $10$ times faster than (preconditionned) CG on DC-SBM $2$ (with expected complexities in $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q} + \vert \mathcal{V} \vert \right)$ and $\mathcal{O}\left( \left(1 + \frac{ 1} {q} \right) \vert \mathcal{E} \vert + \vert \mathcal{V} \vert \right)$ respectively).

Our observations so far show that our estimators are competitive with standard methods, 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 $\overline{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 <ref>).

§ RANDOMIZED ANGULAR SYNCHRONIZATION

We describe in this Section how our novel estimators can be successfully applied to angular synchronization. Recall from Section <ref> that, to perform angular synchronization, we aim to minimize the incoherence

\begin{equation}
\argmin_{f \in U(\C)^n} \ \langle f , \mathsf{L}_\theta f \rangle,
\label{eq:synchro}
\end{equation}

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

Spectral relaxation and other existing approaches
Problem (<ref>) is often relaxed to the following form [71, 63]:

\begin{equation}
\argmin_{\Vert f \Vert_2^2 = n} \ \langle f , \mathsf{L}_\theta f \rangle,
\label{eq:sr}
\end{equation}

the solution of which is given by the eigenvector associated to the smallest eigenvalue of $\mathsf{L}_\theta$. This solution differs from the exact solution of Problem (<ref>) (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, 50, 40].

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 \simeq 10^4$).

* The generalized power method proposed in [10], provably reaching the optimal solution of Problem (<ref>) 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]).

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.

§.§ Proposed Approach

One way to solve the spectral relaxation of Equation (<ref>) is to perform an inverse power iteration, by setting $f_0 \in \C^\mathcal{V}$ and iterating:

\begin{equation}
f_{r+1} = \frac{\mathsf{L}_\theta^{-1} f_r}{\Vert \mathsf{L}_\theta^{-1} f_r \Vert_2}
\label{eq:ipm}
\end{equation}

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

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

\begin{equation}
f_{r+1} = \frac{q (\mathsf{L}_\theta + q \mathsf{I})^{-1} f_r}{\Vert q (\mathsf{L}_\theta + q \mathsf{I})^{-1} f_r \Vert_2},
\label{eq:power_iteration}
\end{equation}

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

\begin{equation}
\frac{\mu_1}{\mu_2} = \frac{\lambda_2 + q}{\lambda_1 + q},
\end{equation}

where $\mu_1 \leq \mu_2$ and $\lambda_1 \leq \lambda_2$ denote respectively the two smallest eigenvectors of $q (\mathsf{L}_\theta + q \mathsf{I})^{-1}$ and $\mathsf{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 $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q} \right)$).

One could instead maximize $\langle f , \mathsf{A}_\theta f \rangle$ in Equation (<ref>). This formulation is common in theoretical works, that mostly focus on (often dense) Erdös-Rényi graphs, and applying the power method to $\mathsf{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 $\mathsf{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., $\eps$-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 $\widetilde{\mathsf{L}_\theta} = \mathsf{D}^{-\frac{1}{2}} \mathsf{L}_\theta \mathsf{D}^{-\frac{1}{2}}$ in the spectral relaxation (<ref>), which was proposed in [71, 17] and comes with similar guarantees [5].
Our approach is flexible enough to be used in this case. Indeed, our MTSF-based estimators extend to quantities such as $f_\circ = (\mathsf{L}_\theta + q\mathsf{D})^{-1} (qD) g'$ (see Section <ref> of the Supplementary Material), which allows to estimate $q (\widetilde{\mathsf{L}_\theta} + q \mathsf{I})^{-1} g$ by setting $g' = \mathsf{D}^{-\frac{1}{2}} g$ and leveraging the decomposition

\begin{equation}
q (\widetilde{\mathsf{L}_\theta} + q \mathsf{I})^{-1} g = \mathsf{D}^{\frac{1}{2}} ((\mathsf{L}_\theta + q\mathsf{D})^{-1} (q\mathsf{D})) g'.
\end{equation}

§.§ 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:

\begin{equation}
I_i = r_i(I_*) + \eps,
\label{eq:image_degradation}
\end{equation}

with $r_i$ a rotation and $\eps \sim \mathcal{N}(0,\sigma^2 \mathsf{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 \rightarrow \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][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], 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 \times 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})$[We restrict the possible rotations so that we can use the simple image-moment-based registration, other methods could be used instead.]. We display the recovered images in Figure <ref> (solving the spectral relaxation (<ref>) using both our method, with $m = 10$ MTSFs for the estimator $\overline{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) = \Vert I_* - x \Vert_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 <ref> if $\cos(\theta_C) \leq 0$, and apply estimator $\overline{f}$ on the resulting MTSF.

$r(I_*)+\eps, \sigma^2 = 1$

$\sigma^2 = 1$, exact

$\sigma^2 = 1$, MTSF

$r(I_*)+\eps, \sigma^2 = 5$

$\sigma^2 = 5$, exact

$\sigma^2 = 5$, MTSF

Error for $\sigma^2 = 5$, $m \in \{1,...,10\}$,

Examples of recovery from noisy images, with different noise levels, using both the exact solution to the spectral relaxation (<ref>) and the approach from Section <ref>.

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 <ref>, and the topology of the graph.

§.§ Numerical Evaluation

We compare the performance of the method from Section <ref> on synthetic graph data using four different iterative solvers: $\overline{f}$, $\hat{f}$ and a conjugate-gradient descent with and without diagonal preconditioning.

We work with the SBM and DC-SBM graph models discussed in Section <ref>. We endow each graph $G$ with a random connection
generated according to the model described in Section <ref>, for both coherent ($\eta = \frac{\pi}{2 n}$) and incoherent connections ($\eta =  \frac{\pi}{10}$), and aim to recover the vector $x$ with  $x_i = e^{\iota \omega_i}$.
Starting from $f_0$ taken uniformly in $U(\C)^n$, we perform $k$ iterations of the power iteration of Equation (<ref>), 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 <ref>[Here, there is no issue with estimating $a f_*$ instead of $f_*$ due to the re-normalization in Equation (<ref>).]. We average the synchronization errors

\begin{equation}
e_s(f) = \min_{r \in U(\C)} \ \frac{\Vert f - rx \Vert_2}{n}
\end{equation}

over $20$ executions, and measure the mean runtimes over $100$ runs. We set $q = 10^{-2} \times \overline{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 <ref>.

We also take measurements for a Lanczos-iteration-based computation[We use the variant implemented in <https://jutho.github.io/KrylovKit.jl/stable/man/eig/#KrylovKit.eigsolve>.] (for the matrix $\mathsf{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} = \vert \cos(\theta_{i,j}) \vert$ (a strategy inspired from surface reconstruction techniques such as [28]).

SBM, $\eta = \frac{\pi}{2n}$, $m = 3$.

DC-SBM $1$, $\eta = \frac{\pi}{2n}$, $m = 3$.

DC-SBM $1$, $\eta = \frac{\pi}{10}$, $m = 10$.

DC-SBM $2$, $\eta = \frac{\pi}{2n}$, $m = 3$.

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.

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

* 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.

* 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.

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

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.

§ 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 implementing 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 <ref> of the Supplementary Material? Could similar estimators be developed for $O(\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 the 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 decompositions such as MTSFs be useful in other iterative inference algorithms, for instance message-passing algorithms [47, 40]?

[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,

[4]
height 2pt depth -1.6pt width 23pt, 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.

[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,

[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 –

[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.

[40]
G. Lerman and Y. Shi, Robust Group Gynchronization 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. Pilavcı, 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), Florianopolis, 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]
height 2pt depth -1.6pt width 23pt, 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 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]
height 2pt depth -1.6pt width 23pt, 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.

Supplementary Material

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

§ PROOF OF PROPOSITION <REF>

Denote by $C: \C^\mathcal{V} \rightarrow \R$ the cost function

\begin{equation}
C(f) = q \Vert f - g \Vert_2^2 + \langle f , \mathsf{L}_\theta f \rangle.
\end{equation}

We use a standard argument from $\C \R$-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 \C^\mathcal{V}$, we have

\begin{align}
C(f + h) & = q \langle (f + h) - g, (f + h) - g \rangle + \langle (f+h), \mathsf{L}_\theta (f+h) \rangle, \\
& = C(f) + D_f^{W}C(h) + D_f^{W^*}C(h) + \left( q \langle h,h \rangle + \langle h, \mathsf{L}_\theta h \rangle \right),
\end{align}

where $(h \mapsto q \langle h,h \rangle + \langle h, \mathsf{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: \C^\mathcal{V} \rightarrow \C$ are given by:

\begin{align}
D_f^W C(h) = \langle q(f-g) + \mathsf{L}_\theta, h \rangle \\
D_f^{W^*} C(h) = \langle h, q (f-g) + \mathsf{L}_\theta \rangle.
\end{align}

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 (\mathsf{L}_\theta + q \mathsf{I})^{-1} g$.

§ A GENERALIZATION OF PROPERTY <REF>

Define $\Delta_\theta = \mathsf{D}^{-1} \mathsf{L}_\theta$ the weighted connection-aware random walk Laplacian $\Delta_\theta$. Proposition <ref> generalizes property <ref>.

For all $l \geq 1$, we have

\begin{equation}
(\mathsf{I} - \Delta_\theta)^l_{i,j} = \sum_{\substack{{p \in P_i^j} \\ {l(p) = l}}} \left( \prod_{0 \leq k < l} \frac{w_{(u_k,u_{k+1})}}{d_{u_k}} \right) \psi_{p^*},
\label{eq:random_walk_laplacian_connection_bis}
\end{equation}

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

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

Observations similar to Proposition <ref> 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].

The result is obvious for $l = 1$. We proceed by induction and assume the result true for some $l \in \N$. Then, recalling that $(\mathsf{I} - \Delta_\theta)_{v,j} = \frac{w_{v,j}}{d_v}$, we have:

\begin{align}
(\mathsf{I} - \Delta_\theta)^{l+1}_{i,j} & = \left( (\mathsf{I} - \Delta_\theta)^{l} (\mathsf{I} - \Delta_\theta) \right)_{i,j} \\
& = \sum_{v \in \mathcal{V}} \sum_{\substack{{p \in P_i^v} \\ {l(p) = l}}} \left( \prod_{0 \leq k < l} \frac{w_{(u_k,u_{k+1})}}{d_{u_k}} \right) \frac{w_{v,j}}{d_v} \psi_{p^*}  \psi_{(v,j)^*} \\
& = \sum_{\substack{{p \in P_i^j} \\ {l(p) = l + 1}}} \left( \prod_{0 \leq k < l + 1} \frac{w_{(u_k,u_{k+1})}}{d_{u_k}} \right) \psi_{p^*}.
\end{align}

For $l = 0$, we have $(\mathsf{I} - \Delta_\theta)^l = \mathsf{I}$.

§ PROOF OF PROPOSITION <REF>

We will show an extension of Proposition <ref> to heterogeneous values of $q$. That is, we associate to each node $i$ some non-negative weight $q_i \in \R$ (at least one of them needs to be positive), that we record in the diagonal matrix $\mathsf{Q}$ with $\mathsf{Q}_{i,i} = q_i$. These $q_i$'s and can be interpreted as edge weights in the extended graph $\mathcal{G}^\Gamma$, and we extend the definition of the measure $\nu_i$ over $P_i^\Gamma$ to account for these weights, so that

\begin{equation}
\prb_{\nu_i}(p) = \mathbf{1}_{\{i\}}(u_0) \mathbf{1}_{\{ \Gamma \}}(u_l) \prod_{0 \leq k < l} \left( \mathbf{1}_{\mathcal{V}\setminus \{\Gamma\}(u_k) } \frac{w_{(u_k,u_{k+1})}}{d_{u_k} + q_{u_k}} \right)
\end{equation}

for a path $p = ((u_0,u_1),...,(u_{l-1},u_l))$ of length $l$, and $\mathbf{1}_S(x)$ the indicator that $x \in S$.
We show the following Feynman-Kac formula, on weighted graphs and for heterogeneous $q_i$'s.

For $j$ the last node reached before absorption of a path $p$, we have

\begin{equation}
\left((\mathsf{L}_\theta + \mathsf{Q})^{-1} \mathsf{Q} g \right)(i) = \mathbf{E}_{p \sim \nu_i} \left(\psi_{p^*_\Gamma(g_j)}\right).
\end{equation}

Note that the l.h.s. is a more general smoothing operation than the solution to the Tikhonov problem, that can be leveraged in the presence of heteroscedastic noise (e.g. for the problem considered in Section <ref> with different values of $\sigma^2$ for each node).

Consider the connection-aware random-walk Laplacian $\Delta_\theta^\mathsf{Q}$ associated to the extended graph $\mathcal{G}^\Gamma$. Its restriction to the rows and columns indexed by $\mathcal{V}$ reads $(\Delta_\theta^q )_\mathcal{V} = (\mathsf{D} + \mathsf{Q})^{-1} (\mathsf{L}_\theta + \mathsf{Q})$. It is clearly invertible, and all its eigenvalues lie in $(0,2)$ (for instance by the Gershgorin circle theorem), hence the eigenvalues of $(\mathsf{I} - (\Delta^\mathsf{Q}_\theta)_\mathcal{V})$ are in $(-1,1)$, ensuring the convergence of the right-hand-side expression in:

\begin{align}
\label{eq:proof_feynmac_kac_1}
\left((\Delta_\theta^q )_\mathcal{V}^{-1}\right)_{i,j} & = \left[ \sum_{l \geq 0} \left(\mathsf{I} - (\Delta^q_\theta)_\mathcal{V} \right)^l \right] \\
& = \sum_{p \in P_i^j} \left( \prod_{0 \leq k < l(p)} \frac{w_{(u_k,u_{k+1})}}{d_{u_{k}} + q_{u_k}} \right) \psi_{p^*}, \label{eq:fk_point}
\end{align}

where Equation <ref> follows from applying Proposition <ref> to $\mathcal{G}^\Gamma$.

The result follows by right multiplication with $(\mathsf{D} + \mathsf{Q})^{-1} \mathsf{Q}$, where $\left( (\mathsf{D} + \mathsf{Q})^{-1} \mathsf{Q} \right)_{j,j} = \frac{q_j}{d_{u_{l-1}} + q_j}$ accounts for the probability of the last transition from $u_{l-1} = j$ to $\Gamma$.

§ PROOF OF THEOREM <REF>

Before presenting the proof of our main Theorem, we state it for weighted graphs with heterogeneous $q_i$'s.

Statement in the presence of weights
Consider the distribution $\mathcal{D}_\mathcal{M}$ over $\mathcal{M}(\mathcal{G})$ such that

\begin{equation}
\prb_{\mathcal{D}_\mathcal{M}}(\phi) \propto \prod_{r \in \phi \cap \mathcal{V}} q_r \prod_{e \in \phi \cap \mathcal{E}} w_e \prod_{C \in \mathcal{C}(\phi)} \left((2 - 2 \cos(\theta_C)\right),
\end{equation}

and the estimator $\tilde{f}(i,\phi,g) = \psi_{r_\phi(i) \xrightarrow[]{\phi} i} \left( g(r_\phi(i)) \right)$ propagating the value from $r_\phi(i)$ to $i$ if $i$ belongs to a rooted tree, and $\tilde{f}(i,\phi,g) = 0$ if $i$ lies in a unicycle.
Then, we have:

\begin{equation}
\left( (\mathsf{L}_\theta + \mathsf{Q})^{-1} \mathsf{Q} g \right) (i) = \mathbf{E}_{\mathcal{D}_\mathcal{M}}(\tilde{f}(i,\phi,g)).
\label{eq:heterogeneous_fermionic}
\end{equation}

We now proceed with the proof.

Determinantal Point Processes
The first step in our proof consists in re-stating the definition of $\mathcal{D}_\mathcal{M}$ as a Determinantal Point Process (DPP), which will provide us with powerful tools for reasoning about MTSFs. A (discrete) DPP over a finite set $\mathcal{X}$ associates a probability to each subset $X \subseteq \mathcal{X}$, and is defined by its marginal probabilities. It is parametrized by an Hermitian matrix $K \in M_{\vert \mathcal{X} \vert} (\C)$, whose eigenvalues must all lie in $[0,1]$, known as the marginal kernel of the DPP, denoted $\mathrm{DPP}(K)$. Precisely, we say that $X \sim \mathrm{DPP}(K)$ if

\begin{equation}
\prb_{\mathrm{DPP}(K)}(A \subseteq X) = \det(K)_{A,A}
\end{equation}

for all $A \subseteq \mathcal{X}$, where $\det(K)_{R,C}$ is the minor of $K$ restricted to the rows and columns indexed respectively by $R$ and $C$. We write $\det(K)_{:,C}$ (resp. $\det(K)_{R,:}$) in case $R = \mathcal{X}$ (resp. $C = \mathcal{X}$).

Marginal kernel of $\mathcal{D}_\mathcal{M}$
Let us now describe the marginal kernel $K$ that we will associate to distribution $\mathcal{D}_\mathcal{M}$. To this end, we consider a twisted discrete differential operator $\nabla: \C^\mathcal{V} \rightarrow \C^\mathcal{E}$, mapping complex values defined on the complex planes $\C_v$ (associated to the nodes of $\mathcal{G}$) to copies of the complex planes $\C_e$ associated to the edges of $\mathcal{G}$. Its expression relies on a splitting of the connection maps $\psi_e: \C_{s_e} \xrightarrow[]{\psi_{s_e,e}} \C_e \xrightarrow[]{\psi_{e,t_e}} \C_{t_e}$  such that $\psi_e = \psi_{e,t_e} \circ \psi_{s_e, e}$ (this is always possible, e.g. $\psi_{e,t_e} = \id_{\C_e,\C_{t_e}}$ and $\psi_{s_e,e}(z) = e^{\iota \theta_e} \cdot z$), and reads [35]:

\begin{equation}
(\nabla f)(e) = \sqrt{w_e} \psi_{t_e,e}(f(t_e)) - \sqrt{w_e} \psi_{s_e,e}(f(s_e)).
\end{equation}

Expliciting the entries of the matrix of $\nabla$, we have:

\begin{equation}
\nabla_{e,v} = \begin{cases}
- \sqrt{w_e} \psi_{s_e,e} & \text{if $v = s_e$} \\
\sqrt{w_e} \psi_{e,t_e} & \text{if $v = t_e$} \\
0 & \text{otherwise.}
\end{cases}
\end{equation}

We finally take $\mathcal{X} = \mathcal{V} \cup \mathcal{E}$ and define the marginal kernel $K$:

\begin{equation}
K = \nabla_\mathsf{Q} (\mathsf{L}_\theta + \mathsf{Q})^{-1} \nabla_Q^*, \text{ where } \nabla_\mathsf{Q} = \begin{bmatrix}
\nabla \\ \sqrt{\mathsf{Q}}
\end{bmatrix}.
\end{equation}

For trivial connections, $\nabla$ is the edge-vertex incidence matrix of $\mathcal{G}$, and we have $\mathsf{L} = \nabla^t \nabla$. In a similar manner, one shows that $\mathsf{L}_\theta = \nabla^* \nabla$.

We will show that $\mathcal{D}_\mathcal{M}$ is a DPP with kernel $K$. First, notice that

\begin{equation}
\mathsf{L}_\theta + \mathsf{Q} = \nabla_\mathsf{Q}^* \nabla_\mathsf{Q},
\end{equation}

so that $K$ is a projection operator (i.e. $K^2 = K$), and all its eigenvalues are in $\{0,1\}$. DPPs associated to such kernels are known as projection DPPs, and can be defined without resorting to marginal probabilities, with $\mathbf{P}_{\mathrm{DPP}(K)}(X) = \det(K)_{X,X}$ when $\vert X \vert = \mathrm{rk}(K)$ the rank of $K$ ($\mathrm{rk}(K) = \vert \mathcal{ V } \vert$ in our case), and $\mathbf{P}_{\mathrm{DPP}(K)}(X) = 0$ otherwise [29].

This will allow to show that the samples of $\mathrm{DPP}(K)$ are MTSFs, distributed according to $\mathcal{D}_\mathcal{M}$ (we generalize here an argument from [35]). Let $\phi \subseteq \mathcal{V} \cup \mathcal{E}$ with $\vert \phi \vert = \vert \mathcal{V} \vert$, and remark that:

\begin{equation}
\det(K)_{\phi,\phi} = \frac{\det(\nabla_\mathsf{Q}^* \nabla_\mathsf{Q})_{\phi,\phi}}{\det(\mathsf{L}_\theta + \mathsf{Q})}.
\end{equation}

The next step then consists in expliciting $\det(\nabla_\mathsf{Q}^* \nabla_\mathsf{Q})_{\phi,\phi} = \det(\nabla_\mathsf{Q})_{\phi,:} \det(\nabla_\mathsf{Q}^*)_{:,\phi}$ depending on $\phi$: we start by inspecting $\det(\nabla_\mathsf{Q}^*)_{:,\phi}$ (similar arguments apply to $\det(\nabla_\mathsf{Q})_{\phi,:} $). In this case, $\phi$ indexes the columns of $\nabla_\mathsf{Q}$, with two columns linearly independent if they do not belong to the same component $c_\phi$. If $c_\phi$ spans $m$ nodes, $\det(\nabla_\mathsf{Q}^*)_{:,c_\phi}$ can only be non-zero if $\vert c_\phi \vert = m$, so that $c_\phi$ must be either a unicycle or a rooted tree. We can then can compute $\det(\nabla_\mathsf{Q}^* \nabla_\mathsf{Q})_{c_\phi,c_\phi}$ from $\det(\nabla_\mathsf{Q}^*)_{:,\phi}$ and $\det(\nabla_\mathsf{Q})_{\phi,:} $ explicitly in those cases.

\begin{equation}
\det(\nabla_\mathsf{Q}^* \nabla_\mathsf{Q})_{c_\phi,c_\phi} = \begin{cases}
(2 - 2\cos(\theta_C)) \prod_{e \in c_\phi} w_e & \text{if $c_\phi$ is a unicycle with cycle $C$} \\
q_r \prod_{e \in c_\phi \cap \mathcal{E}} w_e & \text{if $c_\phi$ is a tree rooted in $r$}
\end{cases}
\end{equation}

Suppose first that $c_\phi$ is a unicycle (this is the case treated in [35]). We fix an orientation of the edges of $c_\phi$ such that all edges are oriented towards the cycle $C$, and that the edges belonging to $C$ form a directed cycle (there are only two such orientations, we choose one arbitrarily). The only non-zero permutations in the determinant

\begin{equation}
\det(\Delta_\mathsf{Q}^*)_{:,c_\phi} = \sum_{\sigma \in \mathcal{S}_m} \prod_{1 \leq i \leq m} \mathrm{sgn}(\sigma) (\Delta_\mathsf{Q}^*)_{i,\sigma(i)},
\end{equation}

where $\mathrm{sgn}(\sigma)$ denotes the signature of permutation (bijection from nodes to edges) $\sigma$,
correspond to these orientations, and map vertices $s_e$ to edges $e$ ($\sigma(s_e) = e$). These two bijections, denoted $\sigma_C$ and $\sigma_{C^*}$, differ only on nodes adjacent to edges in $C$, and correspond to the two possible orientations of the cycle. We then obtain:

\begin{equation}
\det(\Delta_\mathsf{Q}^*)_{:,c_\phi} = (-1)^m \mathrm{sgn}(\sigma_C) \prod_{e \in c_\phi \setminus C} \sqrt{w_e} \psi_{s_e,e} \left( \prod_{e' \in C} \sqrt{w_{e'}} \psi_{s_{e'},e'} - \prod_{e' \in C} \sqrt{w_{e'}} \psi_{t_{e'},e'} \right).
\end{equation}

Performing a similar computation for $\det(\Delta_\mathsf{Q})_{c_\phi,:}$ and multiplying the two resulting expressions then allows to show that (this computation also appears in [35]):

\begin{equation}
\det(\nabla_\mathsf{Q}^* \nabla_\mathsf{Q})_{c_\phi,c_\phi} = (2 - 2\cos(\theta_C)) \prod_{e \in c_\phi} w_e,
\end{equation}

where we used $\psi_e = \psi_{e,t_e} \circ \psi_{s_e, e}$ and $\psi_C + \psi_C^* = 2 \cos(\theta_C)$.

Similarly, for $c_\phi$ a rooted tree with root $r \in c_\phi \cap \mathcal{V}$, we have

\begin{equation}
\det(\Delta_\mathsf{Q}^*)_{:,c_\phi} = \sqrt{q_r} \times (-1)^m \mathrm{sgn}(\sigma_C) \prod_{e \in (c_\phi \cap \mathcal{e}) \setminus C} \sqrt{w_e} \psi_{s_e,e},
\end{equation}

which results in

\begin{equation}
\det(\nabla_\mathsf{Q})_{c_\phi,c_\phi} = q_r \prod_{e \in c_\phi \cap \mathcal{E}} w_e.
\end{equation}

As a corollary, we obtain

\begin{equation}
\prb_{\mathrm{DPP}(K)}(\phi) = \frac{\prod_{r \in \phi \cap \mathcal{V}} q_r \prod_{e \in \phi \cap \mathcal{E}} w_e \prod_{C \in \mathcal{C}(\phi)} \left((2 - 2 \cos(\theta_C)\right)}{\det(\mathsf{L}_\theta + \mathsf{Q})},
\end{equation}

which is exactly $\prb_{\mathcal{D}_\mathcal{M}}(\phi)$.

Unbiased estimator
It remains to show that $\tilde{f}(i,\phi,g)$ is an unbiased estimator of the desired quantity $\left( (\mathsf{L}_\theta + \mathsf{Q})^{-1} \mathsf{Q} g \right) (i)$. Our argument is inspired from that of [51], and we will require another determinantal tool.

For two matrices $A \in M_{m,n}(\C)$ and $B \in M_{n,m}(\C)$ with $m < n$, we have:

\begin{equation}
\det(AB) = \sum_{T} \det(A)_{:,T} \det(B)_{T,:},
\end{equation}

where $T$ ranges over all subsets of $\{1,...,n\}$ of size $m$.

We begin by rewriting $(\mathsf{L}_\theta + \mathsf{Q})^{-1}_{i,j}$ as an expectation, starting from Cramer's rule

\begin{equation}
(\mathsf{L}_\theta + \mathsf{Q})^{-1}_{i,j}  = (-1)^{i+j} \frac{\det(\mathsf{L}_\theta + \mathsf{Q})_{\mathcal{V} \setminus \{j\}, \mathcal{V} \setminus \{i\}}}{\det(\mathsf{L}_\theta + \mathsf{Q})}.
\end{equation}

We can rewrite the numerator using the Cauchy-Binet formula:

\begin{align}
\det(\mathsf{L}_\theta + \mathsf{Q})_{\mathcal{V} \setminus \{j\}, \mathcal{V} \setminus \{i\}} & = \sum_{\substack{{\phi \subseteq \mathcal{V} \cup \mathcal E} \\ {\vert \phi \vert = \vert \mathcal{V} \vert - 1}}} (-1)^{i+j} \det\left( \begin{bmatrix} (\nabla_\mathsf{Q}^*)_{:,\phi} \ \delta_j \end{bmatrix} \right) \det\left( \begin{bmatrix} (\nabla_\mathsf{Q})_{\phi,:} \\ \delta_i \end{bmatrix} \right),
\end{align}

where $\delta_i$ is the $i$-th vector in the usual basis of $\R^\mathcal{V}$. An argument analogous to Lemma <ref> then shows that the product of determinants

\begin{equation}
\det\left( \begin{bmatrix} (\nabla_\mathsf{Q}^*)_{:,\phi} \ \delta_j \end{bmatrix} \right) \det\left( \begin{bmatrix} (\nabla_\mathsf{Q})_{\phi,:} \\ \delta_i \end{bmatrix} \right)
\end{equation}

can only be non-zero if $\phi$ contains a tree $T_i^j \subseteq \mathcal{E}$ spanning both $i$ and $j$ (with no root), with contribution $\psi_{j \xrightarrow[]{\phi} i} \prod_{e \in T_i^j} w_e$ to the product, and if all the other components are rooted trees or unicycles, with associated contributions described in <ref>.
We thus obtain

\begin{align}
(\mathsf{L}_\theta + \mathsf{Q})^{-1}_{i,j} & = \frac{1}{q_j} \sum_{\phi \in \mathcal{M}(\mathcal{G})} \prb_{\phi \sim \mathcal{D}_\mathcal{M}}(\phi) \mathbf{1}_{c_\phi(j)}(i) \psi_{j \xrightarrow[]{\phi} i} \\
& = \frac{1}{q_j} \mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}} \left(\mathbf{1}_{c_\phi(j)}(i) \psi_{j \xrightarrow[]{\phi} i} \right),
\end{align}

where $\mathbf{1}_{c_\phi(j)}(i)$ is the indicator that $i$ belongs to the set of nodes $c_\phi(j)$ spanned by the tree rooted in $r$. Finally, we have

\begin{align}
\left( (\mathsf{L}_\theta + \mathsf{Q})^{-1} \mathsf{Q} g \right) (i) & = \langle \delta_i, (\mathsf{L}_\theta + \mathsf{Q})^{-1} \mathsf{Q} g \rangle \\
& = \sum_{j \in \mathcal{V}} q_j (\mathsf{L}_\theta + \mathsf{Q})^{-1}_{i,j} g(j) \\
& = \sum_{j \in \mathcal{V}} \mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}\left( \psi_{j \xrightarrow[]{\phi} i} \left( g(j) \right) \mathbf{1}_{c_\phi(j)}(i) \right) \\
& = \mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(\tilde{f}(i,\phi,g)).
\end{align}

§ COMPLEXITY ANALYSIS

We will describe the complexity of different implementations of the MTSF-sampling algorithm.
We recall the MTSF-sampling procedure in Algorithm <ref>. We consider here the generalization to weighted graphs with heterogeneous $q_i$'s, and use a generalized $\mathrm{random\_successor}(u)$ that outputs either $\Gamma$ with probability $\frac{q_u}{d_u + q_u}$, or some node $v$ with probability  $\frac{w_{u,v}}{d_u + q_u}$. This generalized algorithms provably samples MTSFs from the distribution $\mathcal{D}_\mathcal{M}$, with the argument from [22] still applying.

MTSF sampling algorithm [22].
[1]
$\phi \leftarrow \emptyset$
$\phi$ not spanning
Let $i \in \mathcal{V}$ be the first node in the queue not spanned by $\phi$
$u \leftarrow i$ $u$ is the current node of the random walk
$p \leftarrow \epsilon$
$\epsilon$ the empty path
($u \neq \Gamma$) and ($p$ does not intersect $\phi$ or contain a cycle)
$u' \leftarrow \mathrm{random\_successor}(u)$ Move to next node
$u' \neq \Gamma$
$e \leftarrow (u,u')$, $p \leftarrow p e$ Add $e$ to the path $p$
${p}$ contains a cycle $C$
Remove this cycle from $p$ with probability $\cos(\theta_C)$
$u_\mathrm{old} \leftarrow u$, $u \leftarrow u'$ $u_\mathrm{old}$ the previous node
$u = \Gamma$
$\phi \leftarrow \phi \cup p \cup u_\mathrm{old}$ Add the sampled path $p$ and the root $u_\mathrm{old}$ to $\phi$
$\phi \leftarrow \phi \cup p$
Output $\phi$

§.§ Number of Steps

Let us re-state the upper bound on the expected number of steps of the random walk in the sampling Algorithm <ref>. Denote by $T_\phi$ the number of random neighbors sampled to build $\phi \in \mathcal{M}(\mathcal{G})$ during an execution of Algorithm <ref>. Then, we have:

\begin{equation}
\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(T_\phi) \leq \mathrm{tr}\left( (\mathsf{L} + \mathsf{Q})^{-1} (\mathsf{D} + \mathsf{Q})  \right),
\label{eq:bound_mtsf_bis}
\end{equation}

where we abuse notation and write $\phi \sim \mathcal{D}_\mathcal{M}$ for $\phi$ sampled using Algorithm <ref>. This bound is easily obtained by considering the special case of the trivial connection with $\psi_e = \mathrm{id}_{\C_{s_e},\C_{t_e}}$ for all $e \in \overrightarrow{\mathcal{E}}$, in which case Algorithm <ref> samples spanning forests of $\mathcal{G}$ from the distribution mentioned in Remark <ref>, and the expected number of steps is known to be [51]:

\begin{equation}
\mathrm{tr}\left( (\mathsf{L} + \mathsf{Q})^{-1} (\mathsf{D} + \mathsf{Q})  \right).
\end{equation}

This is also the slowest scenario: in this case, $\theta_C = 0$ for all cycles, and $p$ can never contain a cycle in step $\mathrm{6}$ of Algorithm <ref>. Hence,  the only way to exit the $\mathrm{\mathbf{while}}$ loop is to reach $\Gamma$ or to intersect $\phi$, and we obtain the bound of Proposition <ref> as a consequence.

Let us briefly comment on the value of the trace in Equation <ref>. We have:

\begin{equation}
\mathrm{tr} \left( \mathsf{L} + \mathsf{Q})^{-1} (\mathsf{D} + \mathsf{Q}) \right) = \sum_{i \in \mathcal{V}} (\mathsf{L} + \mathsf{Q})^{-1}_{i,i} (q_i + d_i),
\end{equation}

where it is known that $(\mathsf{L} + \mathsf{Q})^{-1}_{i,i}$ equals the probability that $i \in \mathsf{f} \cap \mathcal{V}$ for $\mathsf{f} \sim \mathcal{D}_\mathcal{F}$, and hence lies in $[0,1]$. It follows that $\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(T_\phi)$ is in $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q_\mathrm{min}} \right)$.

Algorithm <ref> can be used in conjunction with Equation <ref> to estimate $q(\widetilde{\mathsf{L}_\theta} + q \mathsf{I})^{-1} g$ (see Remark <ref>). In this context, we take $\mathsf{Q} = q \mathsf{D}$, and obtain the bound $\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(T_\phi) \in \mathcal{O}\left( \left( 1 + \frac{1}{q} \right) \vert \mathcal{V} \vert \right)$.

§.§ Implementation-specific Complexity

The time complexity of Algorithm <ref> does not only depend on the number of steps taken by the random walk, but also on the cost of the cycle-detection (step 10). We propose two counter-based solutions.

One-counter detection
The first strategy consists in dynamically assigning a numerical value $\mathrm{c}_v \in \N$ to each node $v \in \mathcal{V}$ encountered during the random walk (initialized to $\mathrm{c}_v = 0$), based on a global counter $\mathrm{c} \in \N$. This counter is incremented each time a new random walk begins (instruction at line 3), and we set $\mathrm{c}_{u'} = \mathrm{c}$ whenever the random walk reaches $u'$ on the instruction at line 7 (if it not already spanned by MTSF). If a cycle $C$ is created in node $u'$ and discarded, we reset $\mathrm{c}_v = 0$ for all nodes $v$ spanned by $C$ except $u'$. A cycle can then be detected in $\mathcal{O}(1)$ time by checking the value of $\mathrm{c}_{u'}$ at each step of the random walk: there is a cycle if $\mathrm{c}_{u'} = \mathrm{c}$, and no cycle if $\mathrm{c}_{u'} < \mathrm{c}$.

In a run of Algorithm <ref>, resetting the values $\mathrm{c}_v$ when cycles are discarded requires going through at most $T_\phi$ nodes, which results in $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q_\mathrm{min}} \right)$ overall expected time complexity.

In practice, we found it faster to use an implementation of the following strategy.

Multiple-counters detection
The $\mathrm{c}_v$'s need to be reset in the one-counter cycle-detection strategy because, whenever a cycle is discarded, nodes spanned by this cycle should no longer be remembered as spanned by the path $p$. This expensive resetting step can be bypassed by considering multiple counters, in addition to the $\mathrm{c}$'s (indexed by a global $\mathrm{c}$). We will keep track of couples of values $(\mathrm{id}_v,\mathrm{val}_v) \in \N^2$ for each $v$, initialized to $(0,0)$) and with updates based on two global counters $\mathrm{id}, \mathrm{val} \in \N^2$, and of values $\mathrm{cap}(\mathrm{id})$ associated to each $\mathrm{id}$. $\mathrm{id}_v$ should be thought of as the ID of a counter at $v$ (amongst multiple others), and $\mathrm{val}_v$ as its value (we only need to store the value associated to the largest $\mathrm{id}_v$). These counters are reset to $0$ whenever $\mathrm{c}$ is incremented (and a new random walk is initiated), and are otherwise updated by applying the following rules for all nodes $u'$ reached by the random walk (if not already spanned by the MTSF).

* If no cycle is created and discarded at $u'$, set $\mathrm{id}_{u'} = \mathrm{id}$, $\mathrm{val}_{u'} = \mathrm{val}$, and increment $\mathrm{val}$.

* If a cycle is created at node $u'$ and discarded, store the value $\mathrm{cap}(\mathrm{id}) = \mathrm{val}$ (the maximum value for the $\mathrm{id}^{th}$ counter), set $\mathrm{cap}(k) = 0$ for all $k > \mathrm{id}_{u'}$, $\mathrm{val}_{u'} = \mathrm{val}$, and increment $\mathrm{id}$.

Using these counters, a cycle is detected at $u'$ if $\mathrm{val}_{u'} \leq \mathrm{cap}(\mathrm{id}_{u'})$.

The number $\mathrm{id}$ of counters used is unbounded, and we did not derive an expected time complexity for this multiple-counters strategy, but consistently obtained better performance in our measurements when using it.

The $\mathcal{O}\left( \frac{\vert \mathcal{E} \vert}{q_\mathrm{min}} \right)$ complexity can seem daunting for very small values of $q_i$'s, but one should remember that this is only an upper bound, not accounting for the presence of unicycles in the sample. As the $q_i \rightarrow 0$, $\mathrm{DPP}(K)$ approaches a distribution over spanning forests of unicycles, with a Wilson-like sampling algorithm described in [33]. Our implementation strategy can also be applied in this case and translates to the bounds on the number of steps of the random walks they derive (see also [22]).

Finally, we point out another implementation detail.

The propagation maps used  in Algorithm <ref> can be computed during the sampling process and, one actually only needs knowledge of: the root of each tree, whether or not a node $i$ belongs to a rooted tree and, if it does, the propagation map from the root $r_\phi(i)$ to $i$. The actual MTSF $\phi$ is never used in the estimation.

§ VARIANCE REDUCTION

We show that the estimators in Equations <ref> and <ref> are unbiased (Proposition <ref>).

Rao-Blackwell estimator
Let us first re-define $\overline{f}$ for heterogeneous values of $q$. The difference resides in the aggregation function $h_\phi$:

\begin{equation}
h_\phi(r,g) = \frac{\sum_{j \in c_\phi(r)} \psi_{j \xrightarrow[]{\phi} r} g(j) } {\sum_{r' \in c_\phi(r)} q_{r'}},
\end{equation}

where we recall that $c_\phi(r)$ denotes the set of nodes spanned by the tree containing $r$. We still take $\overline{f}(i,\phi,g) = \psi_{r_\phi(i) \xrightarrow[]{\psi} i} (h_\phi(r_\phi(i),g))$, and show that:

$\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}} (\overline{f}(i,\phi,g)) = f_*(i)$.

We will express $\overline{f}$ as a conditional expectation, the result will follow from the law of total expectation. Here, we choose a connected subset of edges $\pi \subseteq \mathcal{E}$, and condition on $\phi \cap \mathcal{E}$ containing $\pi$ as one of its (maximal) components, which we denote by $\pi \sqsubseteq_\mathcal{E} \phi$. For $i$ belonging to the connected component spanned by $\pi$, we have:

\begin{align}
\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(\tilde{f}(i,\phi,g) \ \vert \  \pi \sqsubseteq_\mathcal{E} \phi) & = \sum_{\phi \in \mathcal{M}(\mathcal{G})} \prb_{\mathcal{D}_\mathcal{M}}(\phi \ \vert \ \pi \sqsubseteq_\mathcal{E} \phi) \ \tilde{f}(i,\phi,g) \\
& = \sum_{\substack{ {\phi \in \mathcal{M}(\mathcal{G})} \\ {\pi \sqsubseteq_\mathcal{E} \phi}} } \frac{q_r \mathbf{1}_{\phi \cap c_\phi(i)}(r)}{\sum_{r' \in c_\phi(i)} q_{r'}} \ \tilde{f}(i,\phi,g) \\
& = \overline{f}(i,\phi,g).
\end{align}

Gradient step as control variates
We show that $\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}}(\hat{f}(\phi,g)) = f_*$. First recall that $\mathbf{E}_{\mathcal{D}_\mathcal{M}}(\overline{f})~=~f_*$ (and $f_* = q (\mathsf{L}_\theta + q \mathsf{I})^{-1} g$) so that, by linearity of the expectation:

\begin{equation}
\mathbf{E}_{\phi \sim \mathcal{D}_\mathcal{M}} (\hat{f}(\phi,g)) = f_* - \alpha \mathsf{P} (g - g) = f_*
\end{equation}

where, for heterogeneous $q_i$'s, we set $\mathsf{P} = (\mathsf{Q}^{-1} \mathsf{D} + \mathsf{I})^{-1}$.

This shows that $\hat{f}$ is unbiased. Let us now discuss the choice of $\alpha$. We plot in Figure <ref> the mean errors $\Vert \overline{f} - f_* \Vert$ and $\Vert \hat{f} - f_* \Vert$ over $5$ trials, for different values values of $\alpha$. We always take $m = 20$ MTSFs, and set $q = 10^{-2} \times \overline{d}$. We use the random graph models from Section <ref>.

$\eps$-graph mean error.

SBM mean error.

DC-SBM $1$ mean error.

AS CAIDA mean error.

Mean error as a function of $\alpha$.

The extent of the error reduction depends on the graph, but the optimal choice of $\alpha$ seems to always be close to $1$. This behavior can be observed across order-of-magnitude variations of $q$ (not shown), and setting $\alpha = 1$ resulted in meaningful variance reduction for our other experiments (e.g. Section <ref>). This contrasts with the (similar) variance reduction technique proposed in [52], developed for connection-free graphs, for which a good choice of step-size was less straightforward.

§ SUPPORTING ARGUMENTS FOR TIKHONOV SMOOTHING

The goal of this Section is to argue that, when the connection is sufficiently consistent, eigenvectors of $\mathsf{L}$ and $\mathsf{L}_\theta$ can be used similarly in graph-signal-processing applications, such as the smoothing experiment in Section <ref>.

Coherent connection
Let us first consider a perfectly coherent connection, such that $\theta_{i,j} = \omega_j - \omega_i$ for all edges $\{i,j\}$. Consider the eigendecomposition $\mathsf{L} = \mathsf{U} \mathsf{\Lambda} \mathsf{U}^*$ of $\mathsf{L}$, and denote respectively by $x$ the vector with entries $x_i = e^{- \iota \omega_i}$, and $\mathsf{D}_x$ the diagonal matrix with $(\mathsf{D}_x)_{i,i} = x_i$. In this situation, we have

\begin{equation}
\mathsf{L}_\theta = \mathsf{D}_x \mathsf{L} \mathsf{D}_x^* = (\mathsf{D}_x \mathsf{U}) \mathsf{\Lambda} (\mathsf{D}_x \mathsf{U})^*,
\label{eq:eigenvectors}
\end{equation}

and the eigenvectors $v_i$ of $\mathsf{L}_\theta$ are given by $\mathsf{D}_x u_i$, with $u_i$ the eigenvectors of $\mathsf{L}$.

Noisy connection
Equation <ref> no longer holds for incoherent connections. We provide illustrations of the resulting eigenvectors of $\mathsf{L}_\theta$ for connections of the form:

\begin{equation}
\theta_{i,j} = \omega_j - \omega_i + \eta \eps_{i,j},
\end{equation}

for different values of $\eta$, and $\eps_{i,j}$ uniformly distributed in $[-1,1]$ (like in Section <ref>). Specifically, we represent $u_i$, $v_i$, and $\mathsf{D}_x v_i$ (which should approximate $u_i$). We also measure the difference between eigenvectors $u_i$ and $v_i$ as a function of $\eta$ using the normalized error

\begin{equation}
\min_{r \in U(\C)} \ \frac{\Vert u_i - r v_i \Vert_2}{n},
\end{equation}
and use the following graphs:

* An $\eps$-graph built from $n = 100$ points uniformly sampled in $[0,1]^3$, with $\eps = 0.3$ (Figure <ref>).

* A SBM with two communities on size $50$ ($n = 100$), with $c_{1,1} = c_{2,2} = 19$ and $c_{1,2} = 1$ (Figure <ref>).

We only use a single realization of these graphs, equipped with a single realization of the connection. The normalized errors are computed across a linear range of values of $\eta$ in $[0,1]$. For each graph, we plot the real and imaginary parts of the entries of the corresponding eigenvectors, for $\eta_1 \simeq 0.11$ and $\eta_2 \simeq 0.75$. The corresponding errors are also highlighted, along with the value $\eta_0 = \frac{\pi}{2n}$ (for which we know that the connection always satisfies the weak-inconsistency Condition <ref>).
Recall that the eigenvectors of $\mathsf{L}_\theta$ are only defined up to a global rotation, and this is apparent in our illustrations.

$\eta_1 \simeq 0.11$

$\eta_2 \simeq 0.75$

Normalized error as a function of $\eta$. The colored vertical lines represent the values of $\eta_0$ (orange), $\eta_1$ (green) and $\eta_2$ (purple).

Results for the $\eps$-graph. Left: second eigenvector. Right: third eigenvector.

Second eigenvector for the SBM-graph. The two communities are depicted as triangles and squares respectively. Left: $\eta_1 \simeq 0.11$. Right: $\eta_2 \simeq 0.75$.

Finally, we provide similar illustrations for a cyclic graph on size $n = 100$, in a more controlled setup: the values of $\eps_{i,j}$ are no longer randomly sampled, and we instead set $\eps_e = 1$ for all $e$ along a coherent orientation. This way, $\eta_0$ corresponds exactly to the weak-inconsistency threshold in Condition <ref>. We plot the corresponding normalized errors in Figure <ref>.

Normalized errors for the cycle-graph. Left: second eigenvector. Right: third eigenvector. Vertical lines represent $\eta_0$ (orange), $\eta_1$ (green) and $\eta_2$ respectively (purple).

We observe in Figure <ref> that $\mathsf{D}_x^* v_i$ provides a good approximation of $u_i$ for low levels of noise ($\eta = \eta_1$), but this is no longer the case at higher noise-levels ($\eta = \eta_2$): $\mathsf{D}_x^* v_2$ somewhat recovers the shape of $u_2$, but $\mathsf{D}_x ^* v_3$ appears completely unrelated to $u_3$. Results for the SBM (Figure <ref>) suggest that major structural properties of the graph (here, the community structure) are still captured even for high noise: for both $\eta = \eta_1$ and $\eta = \eta_2$, $\mathsf{D}_x^* v_2$ clearly partition the graph into the same communities as $u_i$.

These behaviors are reflected in the normalized-error plots (Figure <ref>), with greater error associated to higher values of $\eta$, and smooth decay. Results on the cycle-graph (Figure <ref>) are very different, with three different behaviors (similar for both eigenvectors), occurring at sharp thresholds: perfect correspondence with $u_i$ at $\eta = 0$, two small, low-error, plateau (one of them occurring as soon as $\eta > 0$, and containing $\eta = \eta_0$), and a higher, essentially constant, error otherwise. For the SBM (not shown), the error decayed linearly with $\eta$.

These illustrations provide weak evidence that eigenvectors of $\mathsf{L}$ and $\mathsf{L}_\theta$ behave similarly in low-noise regimes. In this spirit, simple instances (e.g. $\omega_v = 0$ for all $v$ but non-trivial noise) may be amenable to perturbative analyses, with $\mathsf{L}_\theta$ understood as an analytic perturbation of $\mathsf{L}$ (this is similar to the setting mentioned in Remark <ref>).
The situation is much more nuanced for strongly-incoherent connections, and likely requires more involved mathematical descriptions.

§ TOWARDS EXTENSIONS: AN ALTERNATIVE CONSTRUCTION

We exhibit how the loop-erasure procedure in Algorithm <ref> can be understood as stemming from the Feynman-Kac formula in Proposition <ref>.
The aim is twofold: first, to better understand the relation between the Feynman-Kac formula in Proposition <ref> and the MTSF-based estimaton of Theorem <ref> and, second, to propose a possible proof outline for generalizations of our work. We are not aware of a similar observation appearing in the literature and, as such, the argument may also be of independent interest. For notational brevity, we consider here the simple case of paths sampled from $\nu_i$ containing a single cycle, the reasoning is straightforwardly extended to paths with multiple cycles at the price of heavier notations.

Loop-erasure from Feynman-Kac
Define the equivalence relation $p \simeq q$ on $P_i^\Gamma$ for paths in the extended graph $\mathcal{G}^\Gamma$ that differ by the orientation of at most one cycle (denoted $P_i^\Gamma(1)$): that is, $p \simeq q$ if $ p = a C b$ and $q = a C^* b$, with $C$ a (possibly empty) cycle. For convenience, we will denote the two elements of each of the induced equivalence classes as $p(C)$ and $p(C^*)$, with $\overline{p(C)}$ denoting the equivalence class itself. For node $i \in \mathcal{V}$, we then have:

\begin{align}
\label{eq:obs}
\sum_{p \in P_i^\Gamma(1)} \prb_{\nu_i}(p) \psi_{p_\Gamma^*} & = \sum_{\overline{p(C)} \in \faktor{P_i^\Gamma(1)}{\simeq}} \left( \prb_{\nu_i}(p(C)) \psi_{p(C)^*_\Gamma} + \prb_{\nu_i}(p(C^*)) \psi_{p(C^*)^*_\Gamma} \right)  \\
& = \sum_{\overline{p(C)} \in \faktor{P_i^\Gamma(1)} {\simeq}} \prb_{\nu_i}(\overline{p(C)}) \left( \psi_{a^*} \circ (\psi_C + \psi_{C^*}) \circ \psi_{b^*} \right)\label{eq:point} \\
& = \sum_{\overline{p(C)} \in \faktor{P_i^\Gamma(1)} {\simeq}} 2 \prb_{\nu_i}(\overline{p_C}) \left( \left( \psi_{a^*} \circ \left(\frac{\psi_C + \psi_{C^*}}{2} \right) \right) \circ \psi_{b^*} \right), \label{eq:point2}
\end{align}

where $p_C = aCb$, and $\prb_{\nu_i}(\overline{p_C}) = \prb_{\nu_i}(p_C)$ in Equation <ref> (we use $\prb_{\nu_i}(p_C)~=~\prb_{\nu_i}(p_{C^*})$).

Let us rephrase this observation: whereas propagation along a cycle $C$ in Proposition <ref> depends on its orientation, in expectation this dependency disappears when propagating not by $\psi_C$ but by $\frac{\psi_C + \psi_{C^*}}{2}$, the factor $2$ in Equation <ref> in turn accounting for the possibility of sampling either $p_C$ or $p_{C^*}$. Loop-erasures appear when probabilistically interpreting the maps $\frac{\psi_C + \psi_{C^*}}{2}$, acting by multiplication with $\cos(\theta_C)$, which, under Condition <ref>, defines a Bernoulli trial for cycle acceptance.

Extending this argument to general paths (with more than one cycle) provides an alternative proof that propagating along the first branch sampled in Algorithm <ref> (starting from node $i$) results in an unbiased estimator of $f_*(i)$. As a corollary of the same analysis, if $\cos(\theta_C) < 0$ for some cycle $C$ encountered in path $p$, one can perform a Bernoulli trial with success probability $-\cos(\theta_C) = (-1) \times \cos(\theta_C)$, and multiply the resulting estimate $-1$ each time such a cycle is constructed in path $p$, extending the estimation to any connection.

Now, the question is: can such a construction be extended to all nodes in $\mathcal{G}$, resulting in unbiased MTSF-based estimators free from the weak-inconsistency Condition <ref>?
Also, could similar observations be leveraged to generalize these estimators to higher-dimensional connections (see Remark <ref>)?