literature [@li2018adaptive; @luo2019robust].

Each edge weight is then computed as $$a_{i,j} = \exp\left\{-(\mathbf{f}^i_l - \mathbf{f}^j_l)^{\top} \mathbf{M}_l (\mathbf{f}^i_l - \mathbf{f}^j_l)\right\},
\label{eq:gaussian_kernel}$$ from which we can calculate the corresponding adjacency matrix and graph Laplacian.

As the graph Laplacian can be computed from edge weights by definition, we convert the graph learning problem to an optimization problem over $\mathbf{M}_l$. We optimize $\mathbf{M}_l$ by minimizing Graph Laplacian Regularizer in Eq. , which enforces the graph to capture the underlying structure of features. Let $s_{i,j} = (x_i-x_j)^2$, we have the following problem formulation for graph learning over all the layers of the network, with the objective as the sum of GLR over all the layers: $$\begin{split}
&\begin{split}
    \min_{\mathbf{M}_l}
     & \sum_l \sum_{\{i,j\}}\exp\left\{-(\mathbf{f}^i_l-\mathbf{f}^j_l)^{\top} \mathbf{M}_l (\mathbf{f}^i_l-\mathbf{f}^j_l) \right\} \, s_{i,j}
     \end{split} \\
& \text{s.t.} \quad \,\mathbf{M}_l \succeq 0, \, \forall l.
\label{eq:optimization}
\end{split}$$

Since $\mathbf{M}_l$ is PSD and symmetric, Eq. is computationally expensive to solve in general. Instead, we decompose $\mathbf{M}_l$ into $$\mathbf{M}_l = \mathbf{R}_l \mathbf{R}^{\top}_l,
    \label{eq:learnable_r}$$ where $\mathbf{R}_l \in \mathbb{R}^{K \times S}$, and $S \le K$ is a hyper-parameter to control the complexity of computing distance between high-dimensional features. When $K$ is large and $S << K$, the number of parameters for $\mathbf{M}_l$ is significantly reduced while assuming 