## 2 Materials and Methods  
### Identifying Graph Topologies from Observations of Stationary Graph Signals  
We start with batch (offline) estimation of a time-invariant network topology as considered in [6]. This context is useful to better appreciate the alternative formulation in Section 2.1.1 and the proposed algorithm in Section 2.1.3. Together, these two elements will be instrumental in transitioning to the online setting studied in Section 2.2. To state the problem, recall the symmetric GSO \(\mathbf{S}\in\mathbb{R}^{N\times N}\) associated with the undirected graph \(\mathcal{G}\). Upon defining the vector of coefficients \(\mathbf{h}:=[h_{0},\ldots,h_{L-1}]^{T}\in\mathbb{R}^{L}\) and the polynomial graph filter \(\mathbf{H}:=\sum_{l=0}^{L-1}h_{l}\mathbf{S}^{l}\in\mathbb{R}^{N\times N}\)[3; 36], the Cayley-Hamilton theorem asserts that the model in (1) can be equivalently reparameterized as  
\[\mathbf{y}=\left(\sum_{l=0}^{L-1}h_{l}\mathbf{S}^{l}\right)\mathbf{x}=\mathbf{ H}\mathbf{x}, \tag{2}\]  
for some particular \(\mathbf{h}\) and filter order \(L\leq N\). Since \(\mathbf{S}\) is a local (one-hop) diffusion operator, \(L\) specifies the (multi-hop) range of vertex interactions when generating \(\mathbf{y}\) from \(\mathbf{x}\). It should now be clear that we are constraining ourselves to observations of signals generated via graph filtering. As we explain next, this assumption leads to a tight coupling between \(\mathbf{S}\) and the second-order statistics of \(\mathbf{y}\). The covariance matrix of \(\mathbf{y}=\mathbf{H}\mathbf{x}\) is given by [recall (2) and \(\mathbf{C}_{\mathbf{x}}=\mathbf{I}_{N}\)]  
\[\mathbf{C}_{\mathbf{y}}:=\mathbb{E}\left[\mathbf{y}\mathbf{y}^{T}\right]= \mathbb{E}\left[\mathbf{H}\mathbf{x}(\mathbf{H}\mathbf{x})^{T}\right]=\mathbf{ H}\mathbb{E}\left[\mathbf{x}\mathbf{x}^{T}\right]\mathbf{H}=\mathbf{H}^{2}. \tag{3}\]We relied on the symmetry of \(\mathbf{H}\) to obtain the third equality, as \(\mathbf{H}\) is a polynomial in the symmetric GSO \(\mathbf{S}\). Using the spectral decomposition of \(\mathbf{S}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{T}\) to express the filter as \(\mathbf{H}=\sum_{l=0}^{L-1}h_{l}(\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{T})^{l }=\mathbf{V}(\sum_{l=0}^{L-1}h_{l}\mathbf{\Lambda}^{l})\mathbf{V}^{T}\), we can readily diagonalize the covariance matrix in (3) as  
\[\mathbf{C}_{\mathbf{y}}=\mathbf{V}\left(\sum_{l=0}^{L-1}h_{l}\mathbf{\Lambda}^ {l}\right)^{2}\mathbf{V}^{T}. \tag{4}\]  
Such a covariance expression is the requirement for a graph signal to be stationary in \(\mathbf{S}\) ([7], Def. 2.b). In other words, if \(\mathbf{y}\) is graph stationary (equivalently if \(\mathbf{C}_{\mathbf{x}}=\mathbf{I}_{N}\)), (4) shows that the eigenvectors of the shift \(\mathbf{S}\), the filter \(\mathbf{H}\), and the covariance \(\mathbf{C}_{\mathbf{y}}\) are _all the same_. Thus given a batch of observations \(\{\mathbf{y}_{t}\}_{t=1}^{T}\) comprising independent realizations of (4), the approach in [6] advocates: (i) forming the _sample covariance_\(\hat{\mathbf{C}}_{\mathbf{y}}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{y}_{t} \mathbf{y}_{t}^{T}\) and extracting its eigenvectors \(\hat{\mathbf{V}}\) as spectral templates of \(\mathcal{G}\); then (ii) recover \(\mathbf{S}\) that is optimal in some sense by estimating its eigenvalues \(\mathbf{\Lambda}=\text{diag}(\lambda_{1},\ldots,\lambda_{N})\). Namely, one solves the inverse problem  
\[\min_{\mathbf{\Lambda},\mathbf{S}\in\mathcal{S}}\,f(\mathbf{S}),\quad\text{ subject to }\,d(\mathbf{S},\hat{\mathbf{V}}\mathbf{\Lambda}\hat{\mathbf{V}}^{T})\leq\epsilon, \tag{5}\]  
which is convex for appropriate choices of the function \(f(\mathbf{S})\), the constraint set \(\mathcal{S}\) and the matrix distance \(d(\cdot,\cdot):\mathbb{R}^{N\times N}\times\mathbb{R}^{N\times N}\mapsto \mathbb{R}_{+}\). The tuning parameter \(\epsilon>0\) accounts for the (finite) sample size-induced errors in estimating \(\mathbf{V}\), and could be set to zero if the eigenvectors were perfectly known.