## 2 Materials and Methods
### Identifying Graph Topologies from Observations of Stationary Graph Signals
#### 2.1.3 Proximal Gradient Algorithm for Batch Topology Identification  
Exploiting the problem structure in (6), a batch proximal gradient (PG) algorithm is developed in this section to recover the network topology; see [11] for a comprehensive tutorial treatment on proximal methods. Based on this module, an online algorithm for tracking the (possibly dynamically-evolving) GSO from streaming graph signals is obtained in Section 2.2. PG methods have been popularized for \(\ell_{1}\)-norm regularized linear regression problems, through the class of iterative shrinkage-thresholding algorithms (ISTA); see e.g., [43; 44; 45]. The main advantage of ISTA over off-the-shelf interior point methods is its computational simplicity. We show this desirable feature can permeate naturally to the topology identification context of this paper, addressing the open algorithmic questions in ([6], Section IV). To make the graph learning problem amenable to this optimization method, we dualize the constraint \(|\mathbf{S}\mathbf{\hat{C}}_{\mathbf{y}}-\mathbf{\hat{C}}_{\mathbf{y}}\mathbf{ S}|_{F}^{2}\leq\epsilon\) in (6) and write the composite, non-smooth optimization  
\[\mathbf{S}^{\star}\in\operatorname*{argmin}_{\mathbf{S}\in\mathcal{S}}F( \mathbf{S}):=|\mathbf{S}|_{1}+\underbrace{\frac{\mu}{2}|\mathbf{S}\mathbf{\hat {C}}_{\mathbf{y}}-\mathbf{\hat{C}}_{\mathbf{y}}\mathbf{S}|_{F}^{2}}_{g( \mathbf{S})}. \tag{8}\]  
The quadratic function \(g(\cdot)\) is convex and \(M\)-smooth [i.e., \(\
abla g(\cdot)\) is \(M\)-Lipschitz continuous] and \(\mu>0\) is a tuning parameter. Notice that (8) and (6) are equivalent optimization problems, since for each \(\epsilon\) there exists a value of \(\mu\) such that the respective minimizers coincide. To derive the PG iterations, first notice that the gradient of \(g(\mathbf{S})\) in (8) has the simple form  
\[\
abla g(\mathbf{S})=\mu\big{[}(\mathbf{S}\mathbf{\hat{C}}_{\mathbf{y}}- \mathbf{\hat{C}}_{\mathbf{y}}\mathbf{S})\mathbf{\hat{C}}_{\mathbf{y}}-\mathbf{ \hat{C}}_{\mathbf{y}}(\mathbf{S}\mathbf{\hat{C}}_{\mathbf{y}}-\mathbf{\hat{C} }_{\mathbf{y}}\mathbf{S})\big{]}, \tag{9}\]  
which is Lipschitz continuous with constant \(M=4\mu\lambda_{\max}^{2}(\mathbf{\hat{C}}_{\mathbf{y}})\). With \(\alpha>0\) and \(\mathcal{S}\) a convex set, introduce the proximal operator of a function \(\alpha f(\cdot):\mathbb{R}^{N\times N}\to\mathbb{R}\) evaluated at matrix \(\mathbf{M}\in\mathbb{R}^{N\times N}\) as  
\[\mathbf{Z}(\mathbf{M})=\operatorname{prox}_{\alpha f,\mathcal{S}}(\mathbf{M}): =\operatorname*{argmin}_{\mathbf{X}\in\mathcal{S}}\left[f(\mathbf{X})+\frac{1} {2\alpha}|\mathbf{X}-\mathbf{M}|_{F}^{2}\right]. \tag{10}\]  
With these definitions, the PG updates with fixed step size \(\gamma<\frac{2}{M}\) to solve the batch problem (8) are given by (\(k=1,2,\ldots\) denote iterations)  
\[\mathbf{S}_{k+1}:=\operatorname{prox}_{\gamma|\cdot|_{1},\mathcal{S}}\left( \mathbf{S}_{k}-\gamma\
abla g(\mathbf{S}_{k})\right). \tag{11}\]  
It follows that GSO refinements are obtained through the composition of a gradient-descent step and a proximal operator. Instead of directly optimizing the composite cost in (8), the PG update rule in (11) can be interpreted as the result of minimizing a quadratic overestimator of \(F(\mathbf{S})\) at judiciously chosen points (here the current iterate \(\mathbf{S}_{k}\)); see [12] for a detailed justification. Evaluating the proximal operator efficiently is key to the success of PG methods. For our specific case of sparse graph learning with partial connectivity information, i.e., \(\mathcal{S}\) as defined in (7) and \(f(\mathbf{S})=\|\mathbf{S}\|_{1}\), the proximal operator \(\mathbf{Z}\) in (10) has entries given by  
\[Z_{ij}(M_{ij})=\left\{ \begin{array}{cc}0,&i=j\\ s_{ij},&(i,j)\in\Omega\\ \max(0,M_{ij}-\alpha),&\text{otherwise.}\end{array} \right.