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

$$ (7.2) \qquad (l - \Delta_{\theta})_{i,j}^{l+1} = ((l - \Delta_{\theta})^l (l - \Delta_{\theta}))_{i,j} $$

$$ (7.3) \qquad = \sum_{v \in \mathcal{V}} \sum_{\substack{p \in P_i^v \\ l(p)=l}} \left( \prod_{0 \le k < l} \frac{w_{(u_k, u_{k+1})}}{d_{u_k}} \right) \frac{w_{v,j}}{d_v} \psi_{p^*} \psi_{(v,j)^*} $$

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

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

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

$$ (8.1) \qquad \mathbf{P}_{\nu_i}(p) = \mathbf{1}_{\{i\}}(u_0)\mathbf{1}_{\{\Gamma\}}(u_l) \prod_{0 \le 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) $$

for a path $p = ((u_0, u_1), \dots, (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.

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

$$ (8.2) \qquad ((\mathbf{L}_{\theta} + \mathbf{Q})^{-1}\mathbf{Q}g)(i) = \mathbf{E}_{p \sim \nu_i} (\psi_{p_i^*}(g_j)). $$

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 3.2 with different values of $\sigma^2$ for each node).

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

$$ (8.3) \qquad ((\Delta_{\theta}^{q})_{\mathcal{V}}^{-1})_{i,j} = \left[ \sum_{l \ge 0} (l - (\Delta_{\theta}^{q})_{\mathcal{V}})^l \right] $$

$$ (8.4) \qquad = \sum_{p \in P_i^j} \left( \prod_{0 \le k < l(p)} \frac{w_{(u_k, u_{k+1})}}{d_{u_k} + q_{u_k}} \right) \psi_{p^*}, $$
Proof. The result is obvious for $l = 1$. We proceed by induction and assume the result true for some $l \in \mathbb{N}$. Then, recalling that $(l - \Delta_{\theta})_{v,j} = \frac{w_{v,j}}{d_v}$, we have:

$$ (7.2) \qquad (l - \Delta_{\theta})_{i,j}^{l+1} = ((l - \Delta_{\theta})^l (l - \Delta_{\theta}))_{i,j} $$

$$ (7.3) \qquad = \sum_{v \in \mathcal{V}} \sum_{\substack{p \in P_i^v \\ l(p)=l}} \left( \prod_{0 \le k < l} \frac{w_{(u_k, u_{k+1})}}{d_{u_k}} \right) \frac{w_{v,j}}{d_v} \psi_{p^*} \psi_{(v,j)^*} $$

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

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

**8. Proof of Proposition 2.2.** We will show an extension of Proposition 2.2 to heterogeneous values of $q$. That is, we associate to each node $i$ some non-negative weight $q_i \in \mathbb{R}$ (at least one of them needs to be positive), that we record in the diagonal matrix $\mathbf{Q}$ with $\mathbf{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

$$ (8.1) \qquad \mathbf{P}_{\nu_i}(p) = \mathbf{1}_{\{i\}}(u_0)\mathbf{1}_{\{\Gamma\}}(u_l) \prod_{0 \le 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) $$

for a path $p = ((u_0, u_1), \dots, (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.

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

$$ (8.2) \qquad ((\mathbf{L}_{\theta} + \mathbf{Q})^{-1}\mathbf{Q}g)(i) = \mathbf{E}_{p \sim \nu_i} (\psi_{p_{\Gamma}^*}(g_j)). $$

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 3.2 with different values of $\sigma^2$ for each node).

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

$$ (8.3) \qquad ((\Delta_{\theta}^{\mathcal{Q}})_{\mathcal{V}}^{-1})_{i,j} = \left[ \sum_{l \ge 0} (l - (\Delta_{\theta}^{\mathcal{Q}})_{\mathcal{V}})^l \right]_{i,j} $$

$$ (8.4) \qquad = \sum_{p \in P_i^j} \left( \prod_{0 \le k < l(p)} \frac{w_{(u_k, u_{k+1})}}{d_{u_k} + q_{u_k}} \right) \psi_{p^*}, $$
where Equation 8.4 follows from applying Proposition 7.1 to $\mathcal{G}^\Gamma$.

The result follows by right multiplication with $(\mathbf{D}+\mathbf{Q})^{-1}\mathbf{Q}$, where $((\mathbf{D}+\mathbf{Q})^{-1}\mathbf{Q})_{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$. ■

**9. Proof of Theorem 2.4.** 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}_M$ over $\mathcal{M}(\mathcal{G})$ such that

$$ (9.1) \qquad \mathbf{P}_{\mathcal{D}_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)} ((2 - 2 \cos(\theta_C))), $$

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

**Theorem 9.1.**

$$ (9.2) \qquad ((\mathbf{L}_{\theta} + \mathbf{Q})^{-1}\mathbf{Q}g)(i) = \mathbf{E}_{\mathcal{D}_M}(\tilde{f}(i, \phi, g)). $$

We now proceed with the proof.

*Determinantal Point Processes.* The first step in our proof consists in re-stating the definition of $\mathcal{D}_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_{|\mathcal{X}|}(\mathbf{C})$, whose eigenvalues must all lie in $[0, 1]$, known as the *marginal kernel* of the DPP, denoted $\text{DPP}(K)$. Precisely, we say that $X \sim \text{DPP}(K)$ if

$$ (9.3) \qquad \mathbf{P}_{\text{DPP}(K)}(A \subseteq X) = \det(K)_{A,A} $$

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}_M$.* Let us now describe the marginal kernel $K$ that we will associate to distribution $\mathcal{D}_M$. To this end, we consider a *twisted discrete differential* operator $\nabla : \mathbf{C}^\mathcal{V} \to \mathbf{C}^\mathcal{E}$, mapping complex values defined on the complex planes $\mathbf{C}_v$ (associated to the nodes of $\mathcal{G}$) to copies of the complex planes $\mathbf{C}_e$ associated to the *edges* of $\mathcal{G}$. Its expression relies on a splitting of the connection maps $\psi_e : \mathbf{C}_{s_e} \xrightarrow{\psi_{s_e,e}} \mathbf{C}_e \xrightarrow{\psi_{e,t_e}} \mathbf{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} = \text{id}_{\mathbf{C}_e, \mathbf{C}_{t_e}}$ and $\psi_{s_e,e}(z) = e^{i\theta_e} \cdot z$), and reads [35]:

$$ (9.4) \qquad (\nabla f)(e) = \sqrt{w_e} \psi_{t_e, e}(f(t_e)) - \sqrt{w_e} \psi_{s_e, e}(f(s_e)). $$

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

$$ (9.5) \qquad \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} $$
where Equation (8.4) follows from applying Proposition 7.1 to $\mathcal{G}^\Gamma$.
The result follows by right multiplication with $(\mathbf{D}+\mathbf{Q})^{-1}\mathbf{Q}$, where $((\mathbf{D}+\mathbf{Q})^{-1}\mathbf{Q})_{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$. ■

**9. Proof of Theorem 2.4.** 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}_M$ over $\mathcal{M}(\mathcal{G})$ such that
$$ (9.1) \qquad \mathbf{P}_{\mathcal{D}_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)} ((2 - 2 \cos(\theta_C))), $$

and the estimator $\tilde{f}(i, \phi, g) = \psi_{r_{\phi}(i) \to i}(g(r_{\phi}(i)))$ 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.

**Theorem 9.1.** Letting $\|f\|_{\mathbf{Q}}^2 = \sum_{i \in \mathcal{V}} q_i |f(i)|^2$, we have both:

$$ (9.2) \qquad ((\mathbf{L}_{\theta} + \mathbf{Q})^{-1}\mathbf{Q}g)(i) = \mathbf{E}_{\mathcal{D}_M}(\tilde{f}(i, \phi, g)), $$

$$ (9.3) \qquad \mathbf{E}_{\phi \sim \mathcal{V}} (\|\tilde{f}(\phi, g) - f_*\|_{\mathbf{Q}}^2) = \|g\|_{\mathbf{Q}}^2 - \|f_*\|_{\mathbf{Q}}^2. $$

We now proceed with the proof.

**Determinantal Point Processes.** The first step in our proof consists in re-stating the definition of $\mathcal{D}_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_{|\mathcal{X}|}(\mathbf{C})$, whose eigenvalues must all lie in $[0, 1]$, known as the *marginal kernel* of the DPP, denoted $\text{DPP}(K)$. Precisely, we say that $X \sim \text{DPP}(K)$ if

$$ (9.4) \qquad \mathbf{P}_{\text{DPP}(K)}(A \subseteq X) = \det(K)_{A,A} $$

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}_M$.** Let us now describe the marginal kernel $K$ that we will associate to distribution $\mathcal{D}_M$. To this end, we consider a *twisted discrete differential* operator $\nabla : \mathbf{C}^\mathcal{V} \to \mathbf{C}^\mathcal{E}$, mapping complex values defined on the complex planes $\mathbf{C}_v$ (associated to the nodes of $\mathcal{G}$) to copies of the complex planes $\mathbf{C}_e$ associated to the *edges* of $\mathcal{G}$. Its expression relies on a splitting of the connection maps $\psi_e : \mathbf{C}_{s_e} \xrightarrow{\psi_{s_e,e}} \mathbf{C}_e \xrightarrow{\psi_{e,t_e}} \mathbf{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} = \text{id}_{\mathbf{C}_e, \mathbf{C}_{t_e}}$ and $\psi_{s_e,e}(z) = e^{i\theta_e} \cdot z$), and reads [41]:

$$ (9.5) \qquad (\nabla f)(e) = \sqrt{w_e} \psi_{t_e, e}(f(t_e)) - \sqrt{w_e} \psi_{s_e, e}(f(s_e)). $$

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

$$ (9.6) \qquad \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} $$
We finally take $\mathcal{X} = \mathcal{V} \cup \mathcal{E}$ and define the marginal kernel $K$:

$$ (9.6) \qquad K = \nabla_Q (L_{\theta} + Q)^{-1} \nabla_Q^*, \text{ where } \nabla_Q = \begin{bmatrix} \nabla \\ \sqrt{Q} \end{bmatrix}. $$

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

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

$$ (9.7) \qquad L_{\theta} + Q = \nabla_Q^* \nabla_Q, $$

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}_{\text{DPP}(K)}(X) = \det(K)_{X,X}$ when $|X| = \text{rk}(K)$ the rank of $K$ ($\text{rk}(K) = |\mathcal{V}|$ in our case), and $\mathbf{P}_{\text{DPP}(K)}(X) = 0$ otherwise [29].

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

$$ (9.8) \qquad \det(K)_{\phi,\phi} = \frac{\det(\nabla_Q^* \nabla_Q)_{\phi,\phi}}{\det(L_{\theta} + Q)}. $$

The next step then consists in expliciting $\det(\nabla_Q^* \nabla_Q)_{\phi,\phi} = \det(\nabla_Q)_{\phi,:} \det(\nabla_Q^*)_{:,\phi}$ depending on $\phi$: we start by inspecting $\det(\nabla_Q^*)_{:,\phi}$ (similar arguments apply to $\det(\nabla_Q)_{\phi,:}$). In this case, $\phi$ indexes the columns of $\nabla_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_Q^*)_{:,c_{\phi}}$ can only be non-zero if $|c_{\phi}| = m$, so that $c_{\phi}$ must be either a unicycle or a rooted tree. We can then can compute $\det(\nabla_Q^* \nabla_Q)_{c_{\phi},c_{\phi}}$ from $\det(\nabla_Q^*)_{:,\phi}$ and $\det(\nabla_Q)_{\phi,:}$ explicitly in those cases.

**Lemma 9.3.**

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

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

$$ (9.10) \qquad \det(\Delta_Q^*)_{:,c_{\phi}} = \sum_{\sigma \in S_m} \prod_{1 \le i \le m} \operatorname{sgn}(\sigma)(\Delta_Q^*)_{i,\sigma(i)}, $$

where $\operatorname{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,
We finally take $\mathcal{X} = \mathcal{V} \cup \mathcal{E}$ and define the marginal kernel $K$:

$$ (9.7) \qquad K = \nabla_Q(\mathbf{L}_{\theta} + \mathbf{Q})^{-1}\nabla_Q^*, \quad \text{where } \nabla_Q = \begin{bmatrix} \nabla \\ \sqrt{\mathcal{Q}} \end{bmatrix}. $$

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

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

$$ (9.8) \qquad \mathbf{L}_{\theta} + \mathbf{Q} = \nabla_Q^* \nabla_Q, $$

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}_{\text{DPP}(K)}(X) = \det(K)_{X,X}$ when $|X| = \text{rk}(K)$ the rank of $K$ ($\text{rk}(K) = |\mathcal{V}|$ in our case), and $\mathbf{P}_{\text{DPP}(K)}(X) = 0$ otherwise [34].

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

$$ (9.9) \qquad \det(K)_{\phi,\phi} = \frac{\det(\nabla_Q^* \nabla_Q)_{\phi,\phi}}{\det(\mathbf{L}_{\theta} + \mathbf{Q})}. $$

The next step then consists in expliciting $\det(\nabla_Q^* \nabla_Q)_{\phi,\phi} = \det(\nabla_Q)_{\phi,:} \det(\nabla_Q^*)_{:,\phi}$ depending on $\phi$: we start by inspecting $\det(\nabla_Q^*)_{:,\phi}$ (similar arguments apply to $\det(\nabla_Q)_{\phi,:}$). In this case, $\phi$ indexes the columns of $\nabla_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_Q^*)_{:,c_\phi}$ can only be non-zero if $|c_\phi| = m$, so that $c_\phi$ must be either a unicycle or a rooted tree. We can then can compute $\det(\nabla_Q^* \nabla_Q)_{c_\phi,c_\phi}$ from $\det(\nabla_Q^*)_{:,\phi}$ and $\det(\nabla_Q)_{\phi,:}$ explicitly in those cases.

**Lemma 9.3.**

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

**Proof.** Suppose first that $c_\phi$ is a unicycle (this is the case treated in [41]). 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

$$ (9.11) \qquad \det(\Delta_Q^*)_{:,c_\phi} = \sum_{\sigma \in S_m} \prod_{1 \le i \le m} \operatorname{sgn}(\sigma)(\Delta_Q^*)_{i,\sigma(i)}, $$

where $\operatorname{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:

$$ (9.11) \quad \det(\Delta_Q^*, c_\phi) = (-1)^m \operatorname{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). $$

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

$$ (9.12) \quad \det(\nabla_Q^* \nabla_Q)_{c_\phi, c_\phi} = (2 - 2 \cos(\theta_C)) \prod_{e \in c_\phi} w_e, $$

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

$$ (9.13) \quad \det(\Delta_Q^*)_{c_\phi} = \sqrt{q_r} \times (-1)^m \operatorname{sgn}(\sigma_C) \prod_{e \in (c_\phi \cap \mathcal{V}) \setminus C} \sqrt{w_e} \psi_{s_e, e}, $$

which results in

$$ (9.14) \quad \det(\nabla_Q)_{c_\phi, c_\phi} = q_r \prod_{e \in c_\phi \cap \mathcal{E}} w_e. $$

As a corollary, we obtain

$$ (9.15) \quad \mathbf{P}_{\text{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)} ((2 - 2 \cos(\theta_C)))}{\det(\mathbf{L}_{\theta} + \mathbf{Q})}, $$

which is exactly $\mathbf{P}_{\mathcal{D}_M}(\phi)$.

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

**Proposition 9.4 (Cauchy-Binet formula).** For two matrices $A \in M_{m,n}(\mathbf{C})$ and $B \in M_{n,m}(\mathbf{C})$ with $m < n$, we have:

$$ (9.16) \quad \det(AB) = \sum_{T} \det(A)_{:,T} \det(B)_{T,:}, $$

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

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

$$ (9.17) \quad (\mathbf{L}_{\theta} + \mathbf{Q})^{-1}_{i,j} = (-1)^{i+j} \frac{\det(\mathbf{L}_{\theta} + \mathbf{Q})_{\mathcal{V}\setminus\{j\}, \mathcal{V}\setminus\{i\}}}{\det(\mathbf{L}_{\theta} + \mathbf{Q})}. $$

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

$$ (9.18) \quad \det(\mathbf{L}_{\theta} + \mathbf{Q})_{\mathcal{V}\setminus\{j\}, \mathcal{V}\setminus\{i\}} = \sum_{\substack{\phi \subseteq \mathcal{V} \cup \mathcal{E} \\ |\phi| = |\mathcal{V}|-1}} (-1)^{i+j} \det\left(\left[ (\nabla_Q^*)_{:,\phi} \delta_j \right] \det\left( \begin{bmatrix} (\nabla_Q)_{\phi,:} \\ \delta_i \end{bmatrix} \right)\right), $$
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:

$$ (9.12) \quad \det(\Delta_Q^*, c_\phi) = (-1)^m \operatorname{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). $$

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

$$ (9.13) \quad \det(\nabla_Q^* \nabla_Q)_{c_\phi, c_\phi} = (2 - 2 \cos(\theta_C)) \prod_{e \in c_\phi} w_e, $$

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

$$ (9.14) \quad \det(\Delta_Q^*)_{c_\phi} = \sqrt{q_r} \times (-1)^m \operatorname{sgn}(\sigma_C) \prod_{e \in (c_\phi \cap \mathcal{V}) \setminus C} \sqrt{w_e} \psi_{s_e, e}, $$

which results in

$$ (9.15) \quad \det(\nabla_Q)_{c_\phi, c_\phi} = q_r \prod_{e \in c_\phi \cap \mathcal{E}} w_e. $$

As a corollary, we obtain

$$ (9.16) \quad \mathbf{P}_{\text{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)} ((2 - 2 \cos(\theta_C)))}{\det(\mathbf{L}_{\theta} + \mathbf{Q})}, $$

which is exactly $\mathbf{P}_{\mathcal{D}_M}(\phi)$.

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

**Proposition 9.4 (Cauchy-Binet formula).** For two matrices $A \in M_{m,n}(\mathbf{C})$ and $B \in M_{n,m}(\mathbf{C})$ with $m < n$, we have:

$$ (9.17) \quad \det(AB) = \sum_{T} \det(A)_{:,T} \det(B)_{T,:}, $$

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

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

$$ (9.18) \quad (\mathbf{L}_{\theta} + \mathbf{Q})^{-1}_{i,j} = (-1)^{i+j} \frac{\det(\mathbf{L}_{\theta} + \mathbf{Q})_{\mathcal{V}\setminus\{j\}, \mathcal{V}\setminus\{i\}}}{\det(\mathbf{L}_{\theta} + \mathbf{Q})}. $$

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

$$ (9.19) \quad \det(\mathbf{L}_{\theta} + \mathbf{Q})_{\mathcal{V}\setminus\{j\}, \mathcal{V}\setminus\{i\}} = \sum_{\substack{\phi \subseteq \mathcal{V} \cup \mathcal{E} \\ |\phi| = |\mathcal{V}|-1}} (-1)^{i+j} \det\left(\left[ (\nabla_Q^*)_{:,\phi} \delta_j \right] \det\left( \begin{bmatrix} (\nabla_Q)_{\phi,:} \\ \delta_i \end{bmatrix} \right)\right), $$
where $\delta_i$ is the i-th vector in the usual basis of $\mathbf{R}^V$. An argument analogous to Lemma 9.3 then shows that the product of determinants

$$ (9.19) \qquad \det\left(\left[\left(\nabla_Q^*\right)_{:,\phi} \delta_j\right]\right) \det\left(\begin{bmatrix} \left(\nabla_Q\right)_{\phi,:} \\ \delta_i \end{bmatrix}\right) $$

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 \to 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 9.3. We thus obtain

$$ (9.20) \qquad (L_{\theta} + Q)^{-1}_{i,j} = \frac{1}{q_j} \sum_{\phi \in \mathcal{M}(G)} \mathbf{P}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (\phi) \mathbf{1}_{c_{\phi}(j)}(i) \psi_{j \to i}^{\phi} $$

$$ (9.21) \qquad = \frac{1}{q_j} \mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} \left( \mathbf{1}_{c_{\phi}(j)}(i) \psi_{j \to i}^{\phi} \right), $$

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

$$ (9.22) \qquad \left( (L_{\theta} + Q)^{-1} Q g \right) (i) = \langle \delta_i, (L_{\theta} + Q)^{-1} Q g \rangle $$

$$ (9.23) \qquad = \sum_{j \in \mathcal{V}} q_j (L_{\theta} + Q)^{-1}_{i,j} g(j) $$

$$ (9.24) \qquad = \sum_{j \in \mathcal{V}} \mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} \left( \psi_{j \to i}^{\phi} (g(j)) \mathbf{1}_{c_{\phi}(j)}(i) \right) $$

$$ (9.25) \qquad = \mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (\tilde{f}(i, \phi, g)). $$

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

**10. Complexity Analysis.** We will describe the complexity of different implementations of the MTSF-sampling algorithm. We recall the MTSF-sampling procedure in Algorithm 10.1. We consider here the generalization to weighted graphs with heterogeneous $q_i$'s, and use a generalized randomSuccessor(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.

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

$$ (10.1) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (T_{\phi}) \leq \mathrm{tr} \left( (L+Q)^{-1} (D+Q) \right), $$
1022 where $\delta_i$ is the $i$-th vector in the usual basis of $\mathbf{R}^V$. An argument analogous to Lemma 9.3

1023 then shows that the product of determinants

$$1024 \qquad (9.20) \qquad \det\left(\left[\left(\nabla_Q^*\right)_{:,\phi} \delta_j\right]\right) \det\left(\begin{bmatrix} \left(\nabla_Q\right)_{\phi,:} \\ \delta_i \end{bmatrix}\right)$$

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 \to 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 Lemma 9.3. We thus obtain

$$1025 \qquad (9.21) \qquad (L_\theta + Q)^{-1}_{i,j} = \frac{1}{q_j} \sum_{\phi \in \mathcal{M}(G)} \mathbf{P}_{\phi \sim \mathcal{D}_M}(\phi) \mathbf{1}_{c_\phi(j)}(i) \psi_{j \to i}$$

$$1026 \qquad (9.22) \qquad = \frac{1}{q_j} \mathbf{E}_{\phi \sim \mathcal{D}_M} \left( \mathbf{1}_{c_\phi(j)}(i) \psi_{j \to i}^{\phi} \right),$$

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

$$1027 \qquad (9.23) \qquad \left( (L_\theta + Q)^{-1} Q g \right) (i) = \langle \delta_i, (L_\theta + Q)^{-1} Q g \rangle$$

$$1028 \qquad (9.24) \qquad = \sum_{j \in \mathcal{V}} q_j (L_\theta + Q)^{-1}_{i,j} g(j)$$

$$1029 \qquad (9.25) \qquad = \sum_{j \in \mathcal{V}} \mathbf{E}_{\phi \sim \mathcal{D}_M} \left( \psi_{j \to i}^{\phi} (g(j)) \mathbf{1}_{c_\phi(j)}(i) \right)$$

$$1030 \qquad (9.26) \qquad = \mathbf{E}_{\phi \sim \mathcal{D}_M} (\tilde{f}(i, \phi, g)).$$

**Variance computation.** The expression of the variance relies on the observation that $\tilde{f}(\phi, g)$ is linear in $g$, so that we can write $\tilde{S}_\phi g = \tilde{f}(\phi, g)$ for some linear operator/matrix $\tilde{S}_\phi$, and otherwise consists in a simple computation. Note that we have $\mathbf{E}_{\phi \sim \mathcal{D}_M} (\tilde{S}_\phi g) = f_*$ and

$$1031 \qquad (9.27) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_M} (\tilde{S}_\phi) = (L_\theta + Q)^{-1} Q.$$

First, notice that $\tilde{S}_\phi^* Q \tilde{S}_\phi = Q$. We then have

$$1032 \qquad (9.28) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_M} \left( \|\tilde{S}(\phi, g) - f_*\|_{\mathbb{Q}}^2 \right) = \sum_{i \in \mathcal{V}} q_i \left( \mathbf{E} \left( \left| (\tilde{S}_\phi g)(i) \right|^2 \right) - \left| \mathbf{E} \left( (\tilde{S}_\phi g)(i) \right) \right|^2 \right)$$

$$1033 \qquad (9.29) \qquad = \left( \sum_{i \in \mathcal{V}} \langle g, \mathbf{E}_{\phi \sim \mathcal{D}_M} (\tilde{S}_\phi^*(q_i \delta_i \delta_i^*) \tilde{S}_\phi) g \rangle \right) - \|f_*\|_{\mathbb{Q}}^2$$

$$1034 \qquad (9.30) \qquad = \|g\|_{\mathbb{Q}}^2 - \|f_*\|_{\mathbb{Q}}^2,$$

where $\delta_i$ is usual basis vector of $\mathbf{C}^V$ with non-zero entry $\delta_i(i) = 1$, Equation (9.29) follows from the linearity of the expectation, and Equation (9.30) from the previous observation that $\tilde{S}_\phi^* Q \tilde{S}_\phi = Q$.
<caption>**Algorithm 10.1** MTSF sampling algorithm [22].</caption>

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

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

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

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

Let us briefly comment on the value of the trace in Equation 10.1. We have:

$$ (10.3) \qquad \text{tr} \left( (L+Q)^{-1}(D+Q) \right) = \sum_{i \in \mathcal{V}} (L+Q)_{i,i}^{-1} (q_i + d_i), $$

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

**Remark 10.2.** Algorithm 10.1 can be used in conjunction with Equation 9.2 to estimate $q(\tilde{L}_{\theta} + qI)^{-1}g$ (see Remark 4.1). In this context, we take $Q = qD$, and obtain the bound $\mathbf{E}_{\phi \sim \mathcal{D}_M}(T_{\phi}) \in \mathcal{O}\left((1 + \frac{1}{q}) |\mathcal{V}|\right)$.
Remark 9.5. We claimed in Section 3.2 that the quantity $s(q)$ in Equation (3.4) is given, in expectation, by the number of roots of MTSFs sampled according to $\mathcal{D}_M$. This fact hinges on two well-known properties of determinantal point processes (see, e.g., [44]):

1031. if $X$ follows a DPP with over $\mathcal{X}$ with kernel $K$, and if $S \subseteq \mathcal{X}$, then $X \cap S$ also follows a DPP, with kernel $K_S = K_{S,S}$;

1032. the cardinality $|X|$ of a DPP-sample is in general a random quantity, and its expectation is given by

$$ (9.31) \qquad \mathbf{E}_{X \sim \text{DPP}(K)}(|X|) = \text{tr}(K). $$

In particular, the kernel $K_V$ obtained by restricting the kernel $K$ of $\mathcal{D}_M$ (Equation (9.7)) describes a distribution over the nodes $V$, and reads

$$ (9.32) \qquad K_V = q(L_\theta + qI)^{-1} $$

or, more generally for heterogeneous $q_i$'s:

$$ (9.33) \qquad K_V = \sqrt{Q}(L_\theta + qI)^{-1} \sqrt{Q}. $$

The claim follows from Equation (9.31):

$$ (9.34) \qquad \mathbf{E}_{\phi \sim \mathcal{V}}(|\phi \cap \mathcal{V}|) = \text{tr} \left( \sqrt{Q}(L_\theta + qI)^{-1} \sqrt{Q} \right), $$

which translates to an efficient estimator for the effective degrees of freedom of the connection-aware smoothing problem.

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

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

**Proposition 10.1.**

$$ (10.1) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_M}(T_\phi) \le \text{tr} \left( (L+Q)^{-1}(D+Q) \right), $$

where we abuse notation and write $\phi \sim \mathcal{D}_M$ for $\phi$ sampled using Algorithm 10.1. This bound is easily obtained by considering the special case of the trivial connection with $\psi_e = \text{id}_{C_{se}, C_{te}}$ for all $e \in \vec{\mathcal{E}}$, in which case Algorithm 10.1 samples spanning forests of $\mathcal{G}$ from the distribution mentioned in Remark 2.5, and the expected number of steps is known to be [61]:

$$ (10.2) \qquad \text{tr} \left( (L+Q)^{-1}(D+Q) \right). $$
**10.2. Implementation-specific Complexity.** The time complexity of Algorithm 10.1 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 $c_v \in \mathbb{N}$ to each node $v \in \mathcal{V}$ encountered during the random walk (initialized to $c_v = 0$), based on a global counter $c \in \mathbb{N}$. This counter is incremented each time a new random walk begins (instruction at line 3), and we set $c_{u'} = 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 $c_v = 0$ for all nodes $v$ spanned by $C$ except $u'$. A cycle can then be detected in $O(1)$ time by checking the value of $c_{u'}$ at each step of the random walk: there is a cycle if $c_{u'} = c$, and no cycle if $c_{u'} < c$.

In a run of Algorithm 10.1, resetting the values $c_v$ when cycles are discarded requires going through at most $T_\phi$ nodes, which results in $O\left(\frac{|\mathcal{E}|}{q_{\min}}\right)$ overall expected *time complexity*.

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

*Multiple-counters detection.* The $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 $c$'s (indexed by a global $c$). We will keep track of couples of values $(id_v, val_v) \in \mathbb{N}^2$ for each $v$, initialized to $(0,0)$) and with updates based on two global counters $id, val \in \mathbb{N}^2$, and of values cap(id) associated to each id. $id_v$ should be thought of as the ID of a counter at $v$ (amongst multiple others), and $val_v$ as its value (we only need to store the value associated to the largest $id_v$). These counters are reset to 0 whenever $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 $id_{u'} = id$, $val_{u'} = val$, and increment val.

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

Using these counters, a cycle is detected at $u'$ if $val_{u'} \le cap(id_{u'})$.

The number 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.

**Remark 10.3.** The $O\left(\frac{|\mathcal{E}|}{q_{\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 \to 0$, 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.

**Remark 10.4.** The propagation maps used in Algorithm 2.2 can be computed during the
<caption>**Algorithm 10.1** MTSF sampling algorithm [22].</caption>

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

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

Let us briefly comment on the value of the trace in Equation (10.1). We have:

$$ (10.3) \qquad \operatorname{tr} (\mathbf{L} + \mathbf{Q})^{-1} (\mathbf{D} + \mathbf{Q}) = \sum_{i \in \mathcal{V}} (\mathbf{L} + \mathbf{Q})^{-1}_{i,i} (q_i + d_i), $$

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

**Remark 10.2.** Algorithm 10.1 can be used in conjunction with Equation (9.2) to estimate $q(\tilde{\mathcal{L}}_\theta + q\mathbf{l})^{-1}g$ (see Remark 4.1). In this context, we take $\mathbf{Q} = q\mathbf{D}$, and obtain the bound $\mathbf{E}_{\phi \sim \mathcal{D}_M}(T_\phi) \in \mathcal{O}\left((1 + \frac{1}{q}) |\mathcal{V}|\right)$.

**10.2. Implementation-specific Complexity.** The time-complexity of Algorithm 10.1 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 $c_v \in \mathbb{N}$ to each node $v \in \mathcal{V}$ encountered during the random walk (initialized to $c_v = 0$),
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.

**11. Variance Reduction.** We show that the estimators in Equations 2.10 and 2.12 are unbiased (Proposition 2.11).

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

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

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

**Proposition 11.1.** $E_{\phi \sim D_M}(\bar{f}(i, \phi, g)) = f_*(i)$.

*Proof.* We will express $\bar{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:

$$ (11.2) \qquad E_{\phi \sim D_M}(\tilde{f}(i, \phi, g) | \pi \sqsubseteq \mathcal{E} \phi) = \sum_{\phi \in \mathcal{M}(\mathcal{G})} \mathbf{P}_{D_M}(\phi | \pi \sqsubseteq \mathcal{E} \phi) \tilde{f}(i, \phi, g) $$

$$ (11.3) \qquad = \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) $$

$$ (11.4) \qquad = \bar{f}(i, \phi, g). \qquad \blacksquare $$

*Gradient step as control variates.* We show that $E_{\phi \sim D_M}(\hat{f}(\phi, g)) = f_*$. First recall that $E_{D_M}(\bar{f}) = f_*$ (and $f_* = q(L_\theta + qI)^{-1}g$) so that, by linearity of the expectation:

$$ (11.5) \qquad E_{\phi \sim D_M}(\hat{f}(\phi, g)) = f_* - \alpha \mathbf{P}(g - g) = f_* $$

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

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

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 3). 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.
based on a global counter $c \in \mathbb{N}$. This counter is incremented each time a new random walk begins (instruction at line 3), and we set $c_{u'} = 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 $c_v = 0$ for all nodes $v$ spanned by $C$ except $u'$. A cycle can then be detected in $O(1)$ time by checking the value of $c_{u'}$ at each step of the random walk: there is a cycle if $c_{u'} = c$, and no cycle if $c_{u'} < c$.

In a run of Algorithm 10.1, resetting the values $c_v$ when cycles are discarded requires going through at most $T_{\phi}$ nodes, which results in $O(\frac{|\mathcal{E}|}{q_{\min}})$ overall expected time complexity.

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

**Multiple-counters detection.** The $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 $c$'s (indexed by a global $c$). We will keep track of couples of values $(id_v, val_v) \in \mathbb{N}^2$ for each $v$, initialized to $(0,0)$) and with updates based on two global counters $id, val \in \mathbb{N}^2$, and of values $cap(id)$ associated to each id. $id_v$ should be thought of as the ID of a counter at $v$ (amongst multiple others), and $val_v$ as its value (we only need to store the value associated to the largest $id_v$). These counters are reset to 0 whenever $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 $id_{u'} = id$, $val_{u'} = val$, and increment val.

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

Using these counters, a cycle is detected at $u'$ if $val_{u'} \le cap(id_{u'})$.

The number 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.

**Remark 10.3.** The $O(\frac{|\mathcal{E}|}{q_{\min}})$ 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 \to 0$, DPP($K$) approaches a distribution over spanning forests of unicycles, with a Wilson-like sampling algorithm described in [39]. 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.

**Remark 10.4.** The propagation maps used in Algorithm 2.2 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.

**11. Variance Reduction.** We show that the estimators in Equations (2.12) and (2.14) are unbiased (Proposition 2.9).
<figure><img src="image_22.png" /><figcaption>Figure 11.1: Mean error as a function of α.</figcaption></figure>

**12. Supporting Arguments for Tikhonov Smoothing.** The goal of this Section is to argue that, when the connection is sufficiently consistent, eigenvectors of L and L_θ can be used similarly in graph-signal-processing applications, such as the smoothing experiment in Section 3.

*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 $L = U\Lambda U^*$ of L, and denote respectively by x the vector with entries $x_i = e^{-i\omega_i}$, and $D_x$ the diagonal matrix with $(D_x)_{i,i} = x_i$. In this situation, we have

$$ (12.1) \qquad L_{\theta} = D_{x}LD_{x}^{*} = (D_{x}U)\Lambda(D_{x}U)^{*}, $$

and the eigenvectors $v_i$ of $L_\theta$ are given by $D_x u_i$, with $u_i$ the eigenvectors of L.

*Noisy connection.* Equation 12.1 no longer holds for incoherent connections. We provide illustrations of the resulting eigenvectors of $L_\theta$ for connections of the form:

$$ (12.2) \qquad \theta_{i,j} = \omega_j - \omega_i + \eta\varepsilon_{i,j}, $$
Rao-Blackwell estimator. Let us first re-define $\bar{f}$ for heterogeneous values of $q$. The difference resides in the aggregation function $h_{\phi}$:

1123

1124

1125

$$h_{\phi}(r, g) = \frac{\sum_{j \in c_{\phi}(r)} \psi_{\phi}^{j} g(j)}{\sum_{r' \in c_{\phi}(r)} q_{r'}}$$

1126

1127

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

1128

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

*Proof.* We will express $\bar{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:

$$ (11.2) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (\tilde{f}(i, \phi, g) | \pi \sqsubseteq \mathcal{E} \phi) = \sum_{\phi \in \mathcal{M}(\mathcal{G})} \mathbf{P}_{\mathcal{D}_{\mathcal{M}}} (\phi | \pi \sqsubseteq \mathcal{E} \phi) \tilde{f}(i, \phi, g) $$

$$ (11.3) \qquad = \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) $$

$$ (11.4) \qquad = \bar{f}(i, \phi, g). \qquad \blacksquare $$

1129

1130

1131

The variance of this improved estimator is derived similarly to that of $\tilde{f}$. In particular, $\bar{f}$ is linear in $g$ and, denoting by $\bar{S}_{\phi}$ the operator such that $\bar{S}_{\phi}g = \bar{f}(\phi, g)$, which is such that $\bar{S}_{\phi}^{*}Q\bar{S}_{\phi} = Q\bar{S}_{\phi}$ and $\mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (\bar{S}_{\phi}) = (L_{\theta} + Q)^{-1}Q$, we obtain from the same argument that.

1132

**Proposition 11.2.** $\mathbf{E}_{\phi \sim \mathcal{D}_{\mathcal{M}}} (\|\bar{f}(\phi, g) - f_{*}\|_{Q}^{2}) = \langle g, Qf_{*} \rangle - \|f_{*}\|^{2}$.

1133

Finally, we show the following concentration bound.

1134

**Proposition 11.3.** Let $\varepsilon, \delta \in (0, 1)$, and consider sampling $m$ MTSFs $\{\phi_k\}_{k=1}^m$. Then, if

1135

$$ m \geq \frac{6}{\varepsilon^2} \log \left( \frac{|\mathcal{V}|}{\delta} \right), $$

1136

it holds for any signal $g \in \mathcal{C}^{\mathcal{V}}$ that

1137

$$ (11.6) \qquad \mathbf{P} \left( \left\| \left( \frac{1}{m} \sum_{k=1}^{m} \bar{f}(\phi_k, g) \right) - f_* \right\|_2 \leq \varepsilon \|g\|_2 \right) \geq 1 - \delta. $$

1138

1139

Let us begin by stating the Matrix-Bernstein concentration bound from [80].

1140

1141

1142

**Theorem 11.4 (1.4 from [80]).** Consider a set of $m$ self-adjoint matrices $\{X_k\}_{k=1}^m$ in $M_n(\mathbf{C})$, sampled independently. Assume that, almost surely and for each such matrix, both the conditions

1143

$$ (11.7) \qquad \mathbf{E}(X_k) = 0 $$
for different values of $\eta$, and $\epsilon_{i,j}$ uniformly distributed in $[-1, 1]$ (like in Section 3). Specifically, we represent $u_i$, $v_i$, and $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

$$ (12.3) \qquad \min_{r \in U(C)} \frac{\|u_i - rv_i\|_2}{n}, $$

and use the following graphs:

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

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

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 2.6). Recall that the eigenvectors of $L_\theta$ are only defined up to a global rotation, and this is apparent in our illustrations.

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

**Discussion.** We observe in Figure 12.1 that $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$): $D_x^* v_2$ somewhat recovers the shape of $u_2$, but $D_x^* v_3$ appears completely unrelated to $u_3$. Results for the SBM (Figure 12.2) 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$, $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 12.1), with greater error associated to higher values of $\eta$, and smooth decay. Results on the cycle-graph (Figure 12.3) 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 $L$ and $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 $L_\theta$ understood as an analytic perturbation of $L$ (this is similar to the setting mentioned in Remark 1.1). The situation is much more nuanced for strongly-incoherent connections, and likely requires more involved mathematical descriptions.

**13. Towards Extensions: an Alternative Construction.** We exhibit how the loop-erasure procedure in Algorithm 2.3 can be understood as stemming from the Feynman-Kac formula in Proposition 2.2. The aim is twofold: first, to better understand the relation between the
1144

$$
1145 \qquad (11.8) \qquad \|\mathbf{X}_k\| \le R
$$

$$
\|\mathbf{X}_k\| \le R
$$

1146 are satisfied, where $\|\cdot\|$ denotes the operator norm$^{29}$.

1147 Then, for all $\varepsilon \ge 0$,

$$
1148 \qquad (11.9) \qquad \mathbf{P} \left( \left\| \sum_{k=1}^{m} \mathbf{X}_k \right\| \ge \varepsilon \right) \le n e^{\frac{-\varepsilon^2}{\sigma^2 + R \frac{1}{\sqrt{s}}}}
$$

1149 where we define

$$
1150 \qquad (11.10) \qquad \sigma^2 = \left\| \sum_{k=1}^{m} \mathbf{E} (\mathbf{X}_k^2) \right\|.
$$

1151 Note that this bound is quite powerful and, in particular, that the matrices $X_k$ need not be
1152 sampled from the same distribution. Our setting is quite simple in comparison, as all the
1153 matrices we consider are sampled from the same distribution.

1154 Specifically, we consider the matrices $X_k = \frac{1}{m}(\bar{S}_{\phi_k} - K)$ associated to a set of $m$ multi-
1155 type spanning forests $\{\phi_k\}_{k=1}^m$, with $K = (L_\theta + Q)^{-1}Q$ denoting the “matrix” solution to the
1156 connection-aware Tikhonov smoothing problem. We assume that all $m$ multi-type spanning
1157 forests are sampled according to $\mathcal{D}_M$, and forget the subscripts in the following. Note that
1158 we have

$$
1159 \qquad (11.11) \qquad \mathbf{E}(\mathbf{X}_k) = 0
$$

1160 for all $k$.

1161 In order to determine $R$, we bound the spectrum of each $X_k$:

$$
1162 \qquad (11.12) \qquad \|\mathbf{X}_k\| \le \frac{1}{m} \left( \|\bar{\mathbf{S}}_{\phi_k}\| + \|\mathbf{K}\| \right).
$$

1163 Since we have both $\|\mathbf{K}\| = \frac{q}{q+\lambda_1}$ (where $\lambda_1$ is the smallest eigenvalue of $L_\theta$) and $\|\bar{\mathbf{S}}_{\phi_k}\| \le 1$ (by,  
1164 e.g., Gershgorin's circle theorem), we can take

$$
1165 \qquad (11.13) \qquad R = \frac{1}{m} \left( 1 + \frac{q}{q + \lambda_1} \right).
$$

1166 In order to simplify our bounds later-on, we can notice that

$$
1167 \qquad (11.14) \qquad R \le \frac{3}{m}.
$$

<footnote n="29"><sup>29</sup>Associated to the usual l<sub>2</sub>-norm on C<sup>n</sup>.</footnote>
Feynman-Kac formula in Proposition 2.2 and the MTSF-based estimation of Theorem 2.4 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 $G^\Gamma$ that differ by the orientation of at most one cycle (denoted $P_i^\Gamma(1)$): that is, $p \simeq q$ if $p = aCb$ and $q = aC^*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:

$$ (13.1) \qquad \sum_{p \in P_i^\Gamma(1)} \mathbf{P}_{\nu_i}(p) \psi_{p_i^*} = \sum_{\overline{p(C)} \in P_i^\Gamma(1) / \simeq} \left( \mathbf{P}_{\nu_i}(p(C)) \psi_{p(C)_i^*} + \mathbf{P}_{\nu_i}(p(C^*)) \psi_{p(C^*_i)^*} \right) $$

$$ (13.2) \qquad = \sum_{\overline{p(C)} \in P_i^\Gamma(1) / \simeq} \mathbf{P}_{\nu_i}(\overline{p(C)}) (\psi_{a^*} \circ (\psi_C + \psi_{C^*}) \circ \psi_{b^*}) $$

$$ (13.3) \qquad = \sum_{\overline{p(C)} \in P_i^\Gamma(1) / \simeq} 2\mathbf{P}_{\nu_i}(\overline{p_C}) \left( \left( \psi_{a^*} \circ \left( \frac{\psi_C + \psi_{C^*}}{2} \right) \right) \circ \psi_{b^*} \right), $$

where $p_C = aCb$, and $\mathbf{P}_{\nu_i}(\overline{p_C}) = \mathbf{P}_{\nu_i}(p_C)$ in Equation 13.2 (we use $\mathbf{P}_{\nu_i}(p_C) = \mathbf{P}_{\nu_i}(p_{C^*})$). Let us rephrase this observation: whereas propagation along a cycle $C$ in Proposition 2.2 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 13.3 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 2.6, 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 2.3 (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 $G$, resulting in unbiased MTSF-based estimators free from the weak-inconsistency Condition 2.3? Also, could similar observations be leveraged to generalize these estimators to higher-dimensional connections (see Remark 2.3)?
As for the value of $\sigma^2$, we obtain from the triangle inequality that

$$ (11.15) \qquad \sigma^2 = \frac{1}{m} \|K - K^2\| $$

$$ (11.16) \qquad \leq \frac{1}{m} \left( \frac{q}{q + \lambda_1} + \frac{q^2}{(q + \lambda_n)^2} \right). $$

1168 This expression can be made more interpretable by lower-bounding $\lambda_n$, and we find that:

$$ (11.17) \qquad \sigma^2 \leq \frac{1}{m} \left( \frac{q}{q + \lambda_1} + \frac{q^2}{(1 + q + \max_{v \in \mathcal{V}} d_v)^2} \right), $$

1170 which is proved by leveraging the bounds on $\lambda_n$ in Section 12 below.

1171 It follows from our previous bounds on $R$ and $\sigma^2$ that

$$ (11.18) \qquad \left( \sigma^2 + R \frac{\varepsilon}{3} \right) \leq \frac{1}{m} \left( \varepsilon + \frac{q}{q + \lambda_1} + \frac{q^2}{(1 + q + \max_{v \in \mathcal{V}} d_v)^2} \right). $$

1173 Finally, $\| \sum_{k=1}^m X_k \| \leq \varepsilon$ implies by definition that

$$ (11.19) \qquad \left\| \sum_{k=1}^{m} X_k g \right\|_{2} \leq \varepsilon \|g\|_{2} $$

1175 for any signal $g \in \mathbf{C}^\mathcal{V}$, so that putting all the previous results together results in:

$$ (11.20) \qquad \mathbf{P} \left( \left\| \left( \frac{1}{m} \sum_{k=1}^{m} \bar{S}_{\phi_k} \right) g - f_* \right\|_2 \leq \varepsilon \|g\|_2 \right) \geq 1 - \delta $$

1177 whenever we have

$$ (11.21) \qquad m \geq 2 \frac{\varepsilon + \frac{q}{q + \lambda_1} + \frac{q^2}{(1 + q + \max_{v \in \mathcal{V}} d_v)^2}}{\varepsilon^2} \log \left( \frac{|\mathcal{V}|}{\delta} \right). $$

1179 The bound from Proposition 11.3 follows when taking $\varepsilon < 1$.

1181 Gradient step as control variates. We show that $\mathbf{E}_{\phi \sim \mathcal{D}_M}(\hat{f}(\phi, g)) = f_*$. First recall that
1182 $\mathbf{E}_{\mathcal{D}_M}(\bar{f}) = f_*$ (and $f_* = q(L_\theta + qI)^{-1}g$) so that, by linearity of the expectation:

$$ (11.22) \qquad \mathbf{E}_{\phi \sim \mathcal{D}_M} (\hat{f}(\phi, g)) = f_* - \alpha \mathbf{P}(g - g) = f_* $$

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

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

1188 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 3). This contrasts with the (similar) variance reduction technique proposed in [63], developed for connection-free graphs, for which a good choice of step-size was less straightforward.
<figure><img src="image_23.png" /><figcaption>(a) $\eta_1 \approx 0.11$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) $\eta_2 \approx 0.75$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) 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).</figcaption></figure>

<caption>Figure 12.1: Results for the $\epsilon$-graph. Left: second eigenvector. Right: third eigenvector.</caption>
<figure><img src="image_19.png" /><figcaption>Figure 11.1: Mean error as a function of $\alpha$.</figcaption></figure>

**12. Bounds on $\lambda_n$.** We show that the largest eigenvalue $\lambda_n$ of $L_\theta$ is such that

$$ (12.1) \qquad 1 + \max_{v \in \mathcal{V}} d_v \le \lambda_n \le 2 \max_{v \in \mathcal{V}d_v}, $$

which is nothing but a straightforward generalization of the same bound for the largest eigenvalue of the combinatorial Laplacian.

The right-hand-side inequality is immediately obtained from Gershgorin circle's theorem. To obtain the left-hand-side inequality, let $v_{\max} \in \mathcal{V}$ denote one vertex of $\mathcal{G}$ with maximum degree $d_{v_{\max}} = \max_{v \in \mathcal{V}} d_v$, and consider the signal $f \in \mathbf{C}^\mathcal{V}$ defined by

$$ (12.2) \qquad f(i) = \begin{cases} d_{v_{\max}} & \text{if } i = v_{\max} \\ -e^{i\theta_{v,i}} & \text{if } i \text{ is adjacent to } v_{\max} \\ 0 & \text{otherwise.} \end{cases} $$

Clearly, $\frac{\langle f, L_\theta f \rangle}{\langle f, f \rangle} = 1 + \max_{v \in \mathcal{V}} d_v$, from which the lower-bound on $\lambda_n$ ensues.

**13. Supporting Arguments for Tikhonov Smoothing.** The goal of this Section is to argue that, when the connection is sufficiently consistent, eigenvectors of $L$ and $L_\theta$ can be
<figure><img src="image_26.png" /><figcaption>Figure 12.2: Second eigenvector for the SBM-graph. The two communities are depicted as triangles and squares respectively. Left: $\eta_1 \approx 0.11$. Right: $\eta_2 \approx 0.75$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 12.3: Normalized errors for the cycle-graph. Left: second eigenvector. Right: third eigenvector. Vertical lines represent $\eta_0$ (orange), $\eta_1$ (green) and $\eta_2$ (purple).</figcaption></figure>
used similarly in graph-signal-processing applications, such as the smoothing experiment in Section 3.

**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 $L = U\Lambda U^*$ of $L$, and denote respectively by $x$ the vector with entries $x_i = e^{i\omega_i}$, and $D_x$ the diagonal matrix with $(D_x)_{i,i} = x_i$. In this situation, we have

$$ (13.1) \qquad L_{\theta} = D_{x}LD_{x}^{*} = (D_{x}U)\Lambda(D_{x}U)^{*}, $$

and the eigenvectors $v_i$ of $L_{\theta}$ are given by $D_x^* u_i$, with $u_i$ the eigenvectors of $L$.

**Noisy connection.** Equation (13.1) no longer holds for incoherent connections. We provide illustrations of the resulting eigenvectors of $L_{\theta}$ for connections of the form:

$$ (13.2) \qquad \theta_{i,j} = \omega_j - \omega_i + \eta\epsilon_{i,j}, $$

for different values of $\eta$, and $\epsilon_{i,j}$ uniformly distributed in $[-1, 1]$ (like in Section 3). Specifically, we represent $u_i$, $v_i$, and $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

$$ (13.3) \qquad \min_{r \in U(C)} \frac{\|u_i - r(D_x^* v_i)\|_2}{n}, $$

and use the following graphs:

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

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

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 2.6). Recall that the eigenvectors of $L_{\theta}$ are only defined up to a global rotation, and this is apparent in our illustrations.

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

**Discussion.** We observe in Figure 13.1 that $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$): $D_x^* v_2$ somewhat recovers the shape of $u_2$, but $D_x^* v_3$ appears completely unrelated to $u_3$. Results for the SBM (Figure 13.2) 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$, $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 13.1), with greater error
<figure><img src="image_20.png" /><figcaption>(c) 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).</figcaption></figure>

<caption>Figure 13.1: Results for the $\varepsilon$-graph. Left: second eigenvector. Right: third eigenvector.</caption>

<footnote n="1243">1243 associated to higher values of $\eta$, and smooth decay. Results on the cycle-graph (Figure 13.3)
1244 are very different, with three different behaviors (similar for both eigenvectors), occurring at
1245 sharp thresholds: perfect correspondence with $u_i$ at $\eta = 0$, two small, low-error, plateau (one
1246 of them occurring as soon as $\eta > 0$, and containing $\eta = \eta_0$), and a higher, essentially constant,
1247 error otherwise. For the SBM (not shown), the error decayed linearly with $\eta$.</footnote>
<figure><img src="image_21.png" /><figcaption>Figure 13.2: Second eigenvector for the SBM-graph. The two communities are depicted as triangles and squares respectively. Left: $\eta_1 \approx 0.11$. Right: $\eta_2 \approx 0.75$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 13.3: 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).</figcaption></figure>

These illustrations provide weak evidence that eigenvectors of $L$ and $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 $L_\theta$ understood as an analytic perturbation of $L$ (this is similar to the setting mentioned in Remark 1.1). The situation is much more nuanced for strongly-incoherent connections, and likely requires more involved mathematical descriptions.

**14. Power iteration: $A_\theta$ v.s. $q(L_\theta + qI)^{-1}$.** We compare in the following the performance of our synchronization estimators (Section 4) with that of the power iteration applied to the matrix $A_\theta$ (as proposed in [75]).

*Setup.* We perform the experiment on three different graph/connection models, where the connection is generated according to the model of Equation (3.1) for different incoherence levels.

* With low-incoherence connections (weak-inconsistency condition satisfied), with $\eta = \frac{\pi}{2n}$, on two different random-graph models.

- On the $\varepsilon$-graph model of Section 3.2.

- On the DC-SBM 1 model of Section 3.2.
Note that, for these low-incoherence models, the spectral gaps of the matrices $L_\theta$ and $A_\theta$ should be similar to those of the matrices $L$ and $A$ (which we analyzed in Remark 4.1). Let us stress that low-incoherence connections of this type are encountered in practice, in problems such as ranking [16, 83] for instance.

• In higher-incoherence connections (no guarantee that the weak-inconsistency condition is satisfied), with $\eta = \frac{\pi}{100}$, on the same $\varepsilon$-graph model.

For each graph we estimate the solution $x \in U(\mathbf{C})^V$ to the angular-synchronization problem by performing $k$ power-method iterations, for $k \in \{1, 2, 4, 7, 13, 24, 45, 85, 160, 300\}$. Each power-method iteration consists in either a multiplication by $A_\theta$, or in an approximate smoothing by $\hat{f}$ with $m = 5$ MTSFs (resp., by $m = 5$ iterations of diagonally-preconditioned CG). We set $q = 10^{-4} \times \bar{d}$, and measure for each $k$ the resulting synchronization error (Equation (4.8)) and associated runtime. Each runtime measurement is averaged over 100 runs. The final results are averaged over 10 uniformly-selected initializations $f_0$ of the power method, themselves averaged over 5 graph-realizations for each model.

The results are depicted in Figure 14.1.

**Comments.** There are two main observations.

1. For low-inconsistence connections, the predonimant parameter is the expansivity of the graph. In particular, the spectral gap of the matrix $q(L_\theta + qI)^{-1}$ is much larger than that of $A_\theta$ for the $\varepsilon$-graph, and we obtain a significant gain in performance with both the $\hat{f}$-based and CG-based iterations (for $\hat{f}$, the synchronization error is better in 2 iterations than in 300 iterations for $A_\theta$). The situation is more nuanced for the DC-SBM 1 graph: in that case, the spectral gap of $A_\theta$ is quite small (due to the clear community structure exhibited for this graph model$^{30}$), so that we could expect a similar speed-up, but not actually small enough for the power-method based on $A_\theta$ to converge very slowly (as we observe in Figure 14.1)$^{31}$.

2. For higher-incoherence connections, the bottom eigenvalue $\lambda_1$ of $L_\theta$ (resp. the top eigenvalue $\alpha_1$ of $A_\theta$) gets larger (resp. smaller), which in this example appears to slow-down (resp. speed-up) the convergence of the power-method based on the matrix $q(L_\theta+qI)^{-1}$ (resp. $A_\theta$). This translates to a smaller improvement in the performance of our method against the power-method based on $A_\theta$. Nevertheless, it is hard to derive a general conclusion: the eigenvalues $\lambda_2$ and $\alpha_2$ are also affected by the incoherence of the connection, and the convergence rates of the power-methods depend on the *rates at which the ratios* $\frac{\lambda_2+q}{\lambda_1+q}$ *and* $\frac{\alpha_1}{\alpha_2}$ *evolve* (when the incoherence is low, this problem should be amenable to a pertubative analysis, as sketched in Section 13).

Let us note that the choice of $q$ (which has no bearing on the end-result of the iterations, but does impact the runtime of the inverse power methods) is still subject to the trade-off discussed in Section 4: $q$ should be small enough to allow for a fast convergence of the power iteration, but large enough so that each iteration can be performed efficiently. Here, we purposefully choose a small value of $q$ to showcase the gain in performance of our methods

<footnote>$^{30}$In particular, the parameters of this model ensure that we are well-above the theoretical community-reconstruction threshold for DC-SBMs (see, e.g. [30]).</footnote>

<footnote>$^{31}$For the DC-SBM 2 model, the community structure is even more clear-cut, and we do observe a similar speed-up.</footnote>
against the power-iteration based on $A_\theta$, but do not tailor this choice to the exact graph nor to the incoherence level.

**15. Illustration on MNIST images: an average digit.** We propose a variant of the illustration from Section 4.2 in which the incoherence of the connection is not due to the addition of (controlled) noise, but is inherent to the data itself. Our goal is to extract from a set of similar images an “average representative” of this dataset. Here, we consider a set of 5000 “1”s in the MNIST dataset [46] (a dataset of 28 × 28 handwritten digits), that all represent similar images but, due to different calligraphic styles, present some discrepancies: in particular, the images are not aligned properly, and two very similar digits may differ by a simple rotation.

For each pair of images $(I, J)$ in the dataset, we compute a rotation $r_{I,J}$ registering the two images (based on image moments), extract for each image $I$ its set of $s = 35$ nearest neighbors (with respect to the distance $\|I - r_{J,I}(J)\|_2$ in $\mathbb{R}^{28 \times 28}$), and build the corresponding nearest-neighbor graph and connection.

We then compute a global alignment via angular synchronization, using respectively the exact eigenvector of $L_\theta$ and a MTSF-based estimate (obtained from $\hat{f}$, with $k = 10$ power-method iterations, $m = 5$ MTSFs sampled at each step, and $q = 10^{-4} \times \bar{d}$), and average the 5000 aligned images. For comparison, we also include the average obtained without the rotations.

The results are depicted in Figure 15.1.

**16. Towards Extensions: an Alternative Construction.** We exhibit how the loop-erasure procedure in Algorithm 2.3 can be understood as stemming from the Feynman-Kac formula in Proposition 2.2. The aim is twofold: first, to better understand the relation between the Feynman-Kac formula in Proposition 2.2 and the MTSF-based estimation of Theorem 2.4 and, second, to propose a possible proof outline for generalizations of our work.

*Loop Erasure from Feynman-Kac.* Let us begin by defining an equivalence relation on paths in $P_i^\Gamma$ (recall that $P_i^\Gamma$ denotes the set of paths from $i$ to $\Gamma$ in the extended graph $\mathcal{G}^\Gamma$): two paths $p, q \in P_i^\Gamma$ are equivalent, which we denote by $p \simeq q$, if the only difference between the two is the orientation of (at most) one cycle; that is, $p \simeq q$ if $p = aCb$ and $q = aC^*b$, with $C$ a (possibly empty) cycle in $\mathcal{G}$. For convenience, we denote by $p(C)$ and $p(C^*)$ the two possible representatives of each of the induced equivalence classes, and by $\overline{p(C)}$ the equivalence class itself.

For any node $i \in \mathcal{V}$, we can write the expected value obtained at $i$ by performing parallel transport, from which we obtain:

$$ (16.1) \qquad \sum_{p \in P_i^\Gamma} \mathbf{P}_{\nu_i^\Gamma}(p) \psi_{p^*} = \sum_{\overline{p(C)} \in P_i^\Gamma / \sim} \left( \mathbf{P}_{\nu_i^\Gamma}(p(C)) \psi_{p(C)^*} + \mathbf{P}_{\nu_i^\Gamma}(p(C^*)) \psi_{p(C^*)} \right) $$

$$ (16.2) \qquad = \sum_{\overline{p(C)} \in P_i^\Gamma / \sim} \mathbf{P}_{\nu_i^\Gamma}(\overline{p(C)}) (\psi_{a^*} \circ (\psi_C + \psi_{C^*}) \circ \psi_{b^*}) $$

$$ (16.3) \qquad = \sum_{\overline{p(C)} \in P_i^\Gamma / \sim} 2 \mathbf{P}_{\nu_i^\Gamma}(\overline{p(C)}) \left( \left( \psi_{a^*} \circ \left( \frac{\psi_C + \psi_{C^*}}{2} \right) \right) \circ \psi_{b^*} \right), $$
where $p(C) = aCb$, and we denote by $\mathbf{P}_{\nu_i^\Gamma}(\overline{p(C)}) = \mathbf{P}_{\nu_i^\Gamma}(p(C))$ in Equation 16.2 the probability of sampling either of the two representatives $p(C)$ or $p(C^*)$ (note that $\mathbf{P}_{\nu_i^\Gamma}(p_C) = \mathbf{P}_{\nu_i^\Gamma}(p_{C^*})$). Let us re-phrase this observation in a more qualitative language: whereas each estimation in the Feynman-Kac estimator proceeds by parallel transport along a path $p(C)$, with the transport application depending on the orientation of the cycle $C$, we observe from Equation (16.3) that we can instead perform an orientation-agnostic estimation resulting in the same expectation, by applying the map $\frac{\psi_C + \psi_{C^*}}{2}$ instead of the connection map $\psi_C$; the factor 2 in Equation (16.3) in turn accounts for the possibility of sampling either $p(C)$ or $p(C^*)$.

Under the weak-inconsistency condition, this allows to probabilistically interpret the maps $\frac{\psi_C + \psi_{C^*}}{2}$ in Equation (16.3), which act as multiplications by the factors $\cos(\theta_C)$:

$$ (16.4) \qquad \left( \frac{\psi_C + \psi_{C^*}}{2} \right) (z) = \cos(\theta_C) \cdot z. $$

We can thus understand Equation (16.3) as a variant of the Feynman-Kac estimator involving a Bernoulli trial (with success probability $\cos(\theta_C)$), going as follows. Suppose that the path $p$ contains a unique cycle $C$, and takes the form $p = aCb(j, \Gamma)$, then:

1. if the trial succeeds, set $\psi_{(ab)^*}(g(j))$ as the estimation on node $i$.

2. If the trial fails, set 0 as the estimation.

Let us stress that, in expectation, this procedure still results in the smoothed signal $f_*$.

This argument extends to any number $k$ of cycles, and translates to a significant reduction of variance for the Feynman-Kac estimator$^{32}$: any path $p$ with $k$ cycles, sampled according to $\nu_i^\Gamma$, belongs to an equivalence class of $2^k$ paths, for each of which the Feynman-Kac estimator results in a different value at node $i$, while the modified estimator assigns the same value to all $2^k$ paths.

One can go a step further, and get rid of the weak-inconsistency condition. Indeed, for any cycle $C$ such that $\cos(\theta_C) < 0$, 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 the path $p$, thereby extending the estimation to any $U(\mathbf{C})$-connection.

The crucial observation follows, as one can identify:

* the modified Feynman-Kac estimator, and in particular the Bernoulli trial for each cycle;

* the cycle-acceptance trial of Algorithm 10.1.

More specifically, consider the first node $i$ in the queue of the MTSF-sampling algorithm, at which the first loop-erased random walk starts: this walk constructs a path $p$ and, at each encountered cycle $C$ in $\mathcal{G}$, may be interrupted due to cycle-acceptance, with probability $1 - \cos(\theta_C)$, in which case the ensuing estimation at node $i$ is $\tilde{f}(i) = 0$.

Similarly, consider the modified Feynman-Kac estimator, which performs successive Bernoulli trials with success probability $\cos(\theta_C)$, and notice that the result of the estimation at node $i$ is decided as soon as a trial fails, in which case it is going to be 0. Effectively, the resulting path

<footnote>$^{32}$Note that this requires to be careful in the definition of the successive cycles, as one cycle may be contained within another. Here, a convenient choice of definition is to define the set of cycles according to the following procedure: whenever a cycle is observed in the path $p$ constructed in the random walk, add it to the set of cycles and remove it from $p$.</footnote>
$p$ upon which the parallel transport will be applied describes a part of a multi-type spanning forest: the branch on a unicycle, and the associated cycle.

In both cases, if no cycle is accepted (resp. no Bernoulli trial fails), the walk ends upon reaching node $\Gamma$, and the estimation at node $i$ is described by the parallel transport from $\Gamma$ to $i$ along the path $p \subseteq \mathcal{E}$.

The question is then the following. Can such a construction be extended to all the nodes in $\mathcal{G}$, so as to recover the (unbiasedness of the) MTSF-based estimator $\tilde{f}$? This question is important for a few reasons.

* First and foremost, this construction could both allow to extend MTSF-based estimators to non-weakly-inconsistent connections, and provide a much more intuitive construction of multi-type spanning forests.

* Second, it is likely to be more flexible: instead of interpreting the factors $\frac{\psi_C+\psi_{C^*}}{2}$ probabilistically, one could, e.g., systematically apply the map $\frac{\psi_C+\psi_{C^*}}{2}$, which would result in connection-aware MTSF-based estimators with no unicycles.

**Remark 16.1.** A similar cycle-equivalence argument has also been used to prove the determinantality of the spanning-forest-of-unicycles measure of [25, 41] in [38] (in the weakly-inconsistent setting) and, more recently, to derive the exact expected number of steps for Algorithm 10.1 algorithm in [23].
<figure><img src="image_23.png" /><figcaption>Figure 14.1: Runtime-precision comparisons for the power iteration based on $A_\theta$ [75] v.s. the approximate (inverse) power methods from Section 4. The results obtained from the exact eigenvectors of $L_\theta$ and $A_\theta$ are represented as horizontal lines.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 15.1: Top: some “1”’s from the MNIST dataset. Bottom: average “1” in the dataset, as obtained by three different methods (left: exact synchronization; center: MTSF-based synchronization; right: no synchronization). The results of the synchronization are only given up to a global rotation, which is apparent in the recovered images.</figcaption></figure>