Generalized Laplacian Regularized Framelet Graph Neural Networks

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},W)$ denote a weighted graph, where $\mathcal{V}=\{v_{1},v_{2},\cdots,v_{N}\}$ and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ represent the node set and the edge set, respectively. ${\bf X}\in\mathbb{R}^{N\times c}$ is the feature matrix for $\mathcal{G}$ with $\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{N}\}$ as its rows, and $W=[w_{i,j}]\in\mathbb{R}^{N\times N}$ is the weight matrix on edges with $w_{i,j}>0$ if $(v_{i},v_{j})\in\mathcal{E}$ and $w_{i,j}=0$ otherwise. For undirected graphs, we have $w_{i,j}=w_{j,i}$ which means that $W$ is a symmetric matrix. For directed graphs, it is likely that $w_{i,j}\neq w_{j,i}$ which means that $W$ may not be a symmetric matrix. In most cases, the weight matrix is the graph adjacency matrix, i.e., $w_{i,j}\in\{0,1\}$ with elements $w_{i,j}=1$ if $(v_{i},v_{j})\in\mathcal{E}$ and $w_{i,j}=0$ otherwise. Furthermore, let $\mathcal{N}_{i}=\{v_{j}:(v_{i},v_{j})\in\mathcal{E}\}$ denote the set of neighbours of node $v_{i}$ and $\mathbf{D}=\text{diag}(d_{1,1},...,d_{N,N})\in\mathbb{R}^{N\times N}$ denote the diagonal degree matrix with $d_{i,i}=\sum^{N}_{j=1}w_{i,j}$ for $i=1,...,N$. The normalized graph Laplacian is defined as $\widetilde{\mathbf{L}}=\mathbf{I}-\mathbf{D}^{-\frac{1}{2}}(\mathbf{W}+\mathbf{I})\mathbf{D}^{-\frac{1}{2}}$. Lastly, for any vector $\mathbf{x}=(x_{1},...,x_{c})\in\mathbb{R}^{c}$, we use $\|\mathbf{x}\|_{2}=(\sum^{c}_{i=1}x^{2}_{i})^{\frac{1}{2}}$ to denote its L₂-norm, and similarly for any matrix $\mathbf{M}=[m_{i,j}]$, $\|\mathbf{M}\|:=\|\mathbf{M}\|_{F}=(\sum_{i,j}m^{2}_{i,j})^{\frac{1}{2}}$ is used to denote its Frobenius norm.

Generally speaking, most GNN frameworks are designed under the homophily assumption in which the labels of nodes and neighbours in the graph are mostly identical. The recent work by Zhu et al (2020a) emphasises that the general topology GNN fails to obtain outstanding results on the graphs with different class labels and dissimilar features in their connected nodes, which we call heterophilic or low homophilic graphs. The definition of homophilic and heterophilic graphs are given by:

The first paper that explores the notation of graph $p$-Laplacian and utilizes it in the regularization problem in graph structure data can be traced back to the work (Zhou and Schölkopf, 2005), in which the regularization scheme was further explored to build a flexible GNN learning model in the study (Fu et al, 2022). In this paper, we refer to the notation similar to the one used in the paper (Fu et al, 2022) to define the graph $p$-Laplacian operator. We first define the graph gradient as follows:

Without confusion, we will simply denote $\nabla_{W}$ as $\nabla$ for convenience. For $[i,j]\notin\mathcal{E}$, $(\nabla\mathbf{F})([i,j]):=0$. The graph gradient of $\mathbf{F}$ at a vertex $i$, $\forall i\in\{1,...,N\}$, is defined as $\nabla\mathbf{F}(i):=((\nabla\mathbf{F})([i,1]);\dots;(\nabla\mathbf{F})([i,N]))$ and its Frobenius norm is given by $\|\nabla\mathbf{F}(i)\|_{2}:=(\sum^{N}_{j=1}(\nabla\mathbf{F})^{2}([i,j]))^{\frac{1}{2}}$ which measures the variation of $\mathbf{F}$ around node $i$. Note that we have two explanations for the notation $\nabla\mathbf{F}$: one as a graph gradient (over edges) and one as node gradient (over nodes). The meaning can be inferred from its context in the rest of the paper. We also provide the definition of graph divergence, analogous to the Laplacian operator in continuous setting, which is the divergence of the gradient of a continuous function.

Given the above definitions on graph gradient and divergence, we reach the definition of the graph $p$-Laplacian.

Note that in above definition, we use the fact that $\mathbf{u}$ acting on all nodes gives a vector in $\mathbb{R}^{N}$. Let $\psi_{p}(\mathbf{U})=(\psi_{p}(\mathbf{U}_{;,1}),...,\psi_{p}(\mathbf{U}_{:,N}))$ for $\mathbf{U}\in\mathbb{R}^{N\times N}$. One (Fu et al, 2022) can show that $p$-Laplacian can be decomposed as $\Delta_{p}=\psi_{p}(\mathbf{U})\mathbf{\Lambda}\psi_{p}(\mathbf{U})^{T}$ for some diagonal matrix $\boldsymbol{\Lambda}$.

Framelet is a type of wavelet frame. The study by Sweldens (1998) is the first to present a wavelet with a lifting scheme that provides a foundation in the research of wavelet transform on graphs. With the increase of computational power, Hammond et al (2011) proposed a framework for the wavelet transformation on graphs and employed Chebyshev polynomials to make approximations on wavelets. Dong (2017) further developed tight framelets on graphs. The new design has been applied for graph learning tasks (Zheng et al, 2021b) with great performance enhancement against the classic GNNs. The recent studies show that framelet can naturally decompose the graph signal and re-aggregate them effectively, achieving a significant result on graph noise reduction (Zhou et al, 2021a) with the use of a double term regularizer on the framelet coefficients. Combing with singular value decomposition (SVD), the framelets have been made applicable to directed graphs (Zou et al, 2022).

A simple method of building more versatile and stable framelet families was suggested by Yang et al (2022) which is known as Quasi-Framelets. In this study, we will introduce graph framelets using the same architecture described in the paper (Yang et al, 2022). We begin by defining the filtering functions for Quasi-Framelets:

Particularly $g_{0}$ aims to regulate the highest frequency while $g_{R}$ to regulate the lowest frequency, and the rest to regulate other frequencies in between.

Consider a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with its normalized graph Laplacian $\widetilde{\mathbf{L}}$. Let $\widetilde{\mathbf{L}}$ have the eigen-decomposition $\widetilde{\mathbf{L}}=\mathbf{U}\Lambda\mathbf{U}^{T}$ where $\mathbf{U}$ is the orthogonal spectral bases with its spectra $\Lambda=\text{diag}(\lambda_{1},\lambda_{2},...,\lambda_{N})$ in increasing order. For a given set of Quasi-Framelet functions $\mathcal{F}=\{g_{0}(\xi),g_{1}(\xi),...,g_{R}(\xi)\}$ defined on $[0,\pi]$ and a given level $J$ ($\geq 0$), define the following Quasi-Framelet signal transformation matrices

	$\displaystyle\mathcal{W}_{0,J}$	$\displaystyle=\mathbf{U}g_{0}(\frac{\boldsymbol{\Lambda}}{2^{m+J}})\cdots g_{0}(\frac{\boldsymbol{\Lambda}}{2^{m}})\mathbf{U}^{T},$		(5)
	$\displaystyle\mathcal{W}_{r,0}$	$\displaystyle=\mathbf{U}g_{r}(\frac{\boldsymbol{\Lambda}}{2^{m}})\mathbf{U}^{T},\;\;\text{for }r=1,...,R,$		(6)
	$\displaystyle\mathcal{W}_{r,\ell}$	$\displaystyle=\mathbf{U}g_{r}(\frac{\boldsymbol{\Lambda}}{2^{m+\ell}})g_{0}(\frac{\boldsymbol{\Lambda}}{2^{m+\ell-1}})\cdots g_{0}(\frac{\boldsymbol{\Lambda}}{2^{m}})\mathbf{U}^{T},$		(7)
		$\displaystyle\text{for }r=1,...,R,\ell=1,...,J.$

Note that in the above definition, $m$ is called the coarsest scale level which is the smallest $m$ satisfying $2^{-m}\lambda_{N}\leq\pi$. Denote by $\mathcal{W}=[\mathcal{W}_{0,J};\mathcal{W}_{1,0};...;\mathcal{W}_{R,0};$ $\mathcal{W}_{1,1};...,\mathcal{W}_{R,J}]$ as the stacked matrix. It can be proved that $\mathcal{W}^{T}\mathcal{W}=\mathbf{I}$, thus providing a signal decomposition and reconstruction process based on $\mathcal{W}$. We call this graph Framelet transformation.

In order to alleviate the computational cost imposed by eigendecomposition for the graph Laplacians, the framelet transformation matrices can be approximated by Chebyshev polynomials. Empirically, the implementation of the Chebyshev polynomials $\mathcal{T}^{n}_{r}(\xi)$ with a fixed degree $n$, $n=3$ is sufficient to approximate $g_{r}(\xi)$. In the sequel, we will simplify the notation by using $\mathcal{T}_{r}(\xi)$ instead of $\mathcal{T}^{n}_{r}(\xi)$. Then the Quasi-Framelet transformation matrices are defined in Eqs. (6) - (7) can be approximated by,

	$\displaystyle\mathcal{W}_{0,J}$	$\displaystyle\approx\mathcal{T}_{0}(\frac{1}{2^{m+J}}\widetilde{\mathbf{L}})\cdots\mathcal{T}_{0}(\frac{1}{2^{m}}\widetilde{\mathbf{L}}),$		(8)
	$\displaystyle\mathcal{W}_{r,0}$	$\displaystyle\approx\mathcal{T}_{r}(\frac{1}{2^{m}}\widetilde{\mathbf{L}}),\;\;\;\text{for }r=1,...,R,$		(9)
	$\displaystyle\mathcal{W}_{r,\ell}$	$\displaystyle\approx\mathcal{T}_{r}(\frac{1}{2^{m+\ell}}\widetilde{\mathbf{L}})\mathcal{T}_{0}(\frac{1}{2^{m+\ell-1}}\widetilde{\mathbf{L}})\cdots\mathcal{T}_{0}(\frac{1}{2^{m}}\widetilde{\mathbf{L}}),$		(10)
		$\displaystyle\;\;\;\;\text{for }r=1,...,R,\ell=1,...,J.$

Remark 1: The approximated transformation matrices defined in Eqs. (8)-(10) simply depend on the graph Laplacian. For directed graphs, we directly take in the Laplacian normalized by the out degrees in our experiments. We observe this strategy leads to improved performance in general.

For a graph signal $\mathbf{X}$, the framelet (graph) convolution similar to the spectral graph convolution that can be defined as

$\displaystyle\theta\star\mathbf{X}=\mathcal{W}^{T}(\textrm{diag}(\theta))(\mathcal{W}\mathbf{X}),$

(11)

where $\theta$ is the learnable network filter. We also call $\mathcal{W}\mathbf{X}$ the framelet coefficients of $\mathbf{X}$ in Fourier domain. The signal then will be filtered in its spectral-domain according to learnable filter $\textrm{diag}(\theta)$.

In reality, graph data normally have large, noisy and complex structures (Leskovec and Faloutsos, 2006; He and Kempe, 2015). This brings the challenge of choosing a suitable neural network for fitting the data. It is well known that one of the common computational issues for most classic GNNs is over-smoothing which can be quantified by Dirichlet energy that converges to zero as the number of layers increases. This observation leads to an investigation into the so-called implicit layers on regularized graph neural networks.

The first paper to consider GNNs layer as a regularized signal smoothing process is done by Zhu et al (2021), in which the classic GNN layers are interpreted as the solution of the regularized optimization problem, with certain approximation strategies to avoid matrix inversion in the closed-form solution. The regularized layer also can be linked to the implicit layer (Zhang et al, 2003), more specific examples are given by Zhu et al (2021). In general, given the node features, the output of the GNN layer can be written as the solution for the following optimization problem

$\displaystyle\mathbf{F}=\operatorname*{arg\,min}_{\mathbf{F}}\left\{\mu\|\mathbf{X}-\mathbf{F}\|^{2}_{F}+\frac{1}{2}\text{tr}(\mathbf{F}^{T}\widetilde{\mathbf{L}}\mathbf{F})\right\},$

(12)

where $\mathbf{F}$ is the graph signal. Eq. (12) defines the layer output as the solution for the (layer) optimization problem with a regularization term $\text{tr}(\mathbf{F}^{T}\widetilde{\mathbf{L}}\mathbf{F})\}$, which is able to enforce smoothness of a graph signal. The closed form solution is given by $\mathbf{F}=(\mathbf{I}+\frac{1}{2\mu}\widetilde{\mathbf{L}})^{-1}\mathbf{X}$. It is computationally inefficient because of matrix inversion. However, we can rearrange it to $(\mathbf{I}+\frac{1}{2\mu}\widetilde{\mathbf{L}})\mathbf{F}=\mathbf{X}$ and interpret it as a fixed point solution to $\mathbf{F}^{(t+1)}=-\frac{1}{2\mu}\widetilde{\mathbf{L}}\mathbf{F}^{(t)}+\mathbf{X}$. The latter is then the implicit layer in GNN at layer $t$. Such an iteration is motivated by the recurrent GNNs as shown in several recent works (Gu et al, 2020; Liu et al, 2021a; Park et al, 2021). Different from explicit GNNs, the output features $\mathbf{F}$ of a general implicit GNN are directly modeled as the solution of a well defined implicit function, e.g., the first order condition from an optimization problem. Denote the rows of $\mathbf{F}$ by $\mathbf{f}_{i}$ ($i=1,2,...,N$) in column shape. It is well known that the regularization term (12) is

$$\frac{1}{2}\text{tr}(\mathbf{F}^{T}\widetilde{\mathbf{L}}\mathbf{F})=\frac{1}{2}\sum_{(v_{i},v_{j})\in\mathcal{E}}\left\|\sqrt{\frac{w_{i,j}}{d_{j,j}}}\mathbf{f}_{j}-\sqrt{\frac{w_{i,j}}{d_{i,i}}}\mathbf{f}_{i}\right\|^{2},$$

which is the so-called graph Dirichlet energy (Zhou and Schölkopf, 2005). As proposed in the work (Zhou and Schölkopf, 2005), one can replace the graph Laplacian matrix in the above equation to the pre-defined graph $p$-Laplacian, then the $p$-Laplacian based energy denoted as $\mathcal{S}_{p}(\mathbf{F})$ can be defined as (for any $p\geq 1$):

$\displaystyle\mathcal{S}_{p}(\mathbf{F})=\frac{1}{2}\sum_{(v_{i},v_{j})\in\mathcal{E}}\left\|\sqrt{\frac{w_{i,j}}{d_{j,j}}}\mathbf{f}_{j}-\sqrt{\frac{w_{i,j}}{d_{i,i}}}\mathbf{f}_{i}\right\|^{p},$

(13)

where we adopt the definition of element-wise $p$-norm as in paper (Fu et al, 2022). Finally, the regularization problem becomes

$\displaystyle\mathbf{F}=\operatorname*{arg\,min}_{\mathbf{F}}\left\{\mu\|\mathbf{X}-\mathbf{F}\|^{2}_{F}+\mathcal{S}_{p}(\mathbf{F})\right\}.$

(14)

With strong generalizability of the regularization form in Eq. (14), it is natural to consider to deploy Eq. 14 to multi-scale graph neural networks (i.e., graph framelets) to explore whether the benefits of allocating $p$-Laplacian based regularizer still remains and how the changes of $p$ within included regularizer interacts with the feature propagation process via each individual filtered domains. In the next section, we will accordingly proposed our model and provide detailed discussion and analysis on it.

Instead of simply taking the convolution result as the framelet layer output (derived in Eq. (11)), we apply the $p$-Laplacian regularization on the framelet reconstructed signal by imposing the following optimization to define the new layer,

$\displaystyle\mathbf{F}=\operatorname*{arg\,min}_{\mathbf{F}}\mathcal{S}_{p}(\mathbf{F})+\mu\|\mathbf{F}-\mathcal{W}^{T}\textrm{diag}(\theta)\mathcal{W}\mathbf{X}\|^{2}_{F},$

(15)

where $\mu\in(0,\infty)$. In which we recall that the first term $\mathcal{S}_{p}(\mathbf{F})$ in the above equation is the $p$-Laplacian based energy defined in Eq. (13). $\mathcal{S}_{p}(\mathbf{F})$ measures the total signals’ variation throughout the graph-based format of $p$-Laplacian. The second term dictates that optimal $\mathbf{F}$ should not deviate excessively from the framelet reconstructed signal for the input signal $\mathbf{X}$. Each component of the network filter $\theta$ in the frequency domain is applied to the framelet coefficients $\mathcal{W}\mathbf{X}$. We call the model in Eq.(15) pL-UFG2.

One possible variant to model (15) is to apply the regularization individually on the reconstruction at each scale level, i.e., for all the $r,\ell$, define

	$\displaystyle\mathbf{F}_{r,\ell}$	$\displaystyle=\operatorname*{arg\,min}_{\mathbf{F}_{r,\ell}}\mathcal{S}_{p}(\mathbf{F}_{r,\ell})+\mu\|\mathbf{F}_{r,\ell}-\mathcal{W}^{T}_{r,\ell}\textrm{diag}(\theta_{r,\ell})\mathcal{W}_{r,\ell}\mathbf{X}\|^{2}_{F},$
	$\displaystyle\mathbf{F}$	$\displaystyle=\mathbf{F}_{0,J}+\sum^{R}_{r=1}\sum^{J}_{\ell=0}\mathbf{F}_{r,\ell}.$		(16)

We name the model with the propagation in the above equation as pL-UFG1. In our experiment, we note that pL-UFG1 performs better than pL-UFG2.

In pL-fUFG, we take a strategy of regularizing the framelet coefficients. Informed by earlier experience, we consider the following optimization problem for each framelet transformation in Fourier domain. For each framelet transformation $\mathcal{W}_{r,j}$ defined in Eqs. (9)-(10), define

$\displaystyle\mathbf{F}_{r,\ell}=\operatorname*{arg\,min}_{\mathbf{F}_{r,\ell}}\mathcal{S}_{p}(\mathbf{F}_{r,\ell})+\mu\|\mathbf{F}_{r,\ell}-\textrm{diag}(\theta_{r,\ell})\mathcal{W}_{r,\ell}\mathbf{X}\|^{2}_{F}.$

(17)

Then the final layer output is defined by the reconstruction as follows

$\displaystyle\mathbf{F}=\mathcal{W}^{T}_{0,J}\mathbf{F}_{0,J}+\sum^{R}_{r=1}\sum^{J}_{\ell=0}\mathcal{W}^{T}_{r,\ell}\mathbf{F}_{r,\ell}.$

(18)

Or we can only aggregate all the regularized and filtered framelet coefficients in the following way

$\displaystyle\mathbf{F}=\mathbf{F}_{0,J}+\sum^{R}_{r=1}\sum^{J}_{\ell=0}\mathbf{F}_{r,\ell}.$

(19)

With the form of both pL-UFG and pL-fUFG, in the next section we show an even more generalized regularization framework by incorporating $p$-Laplacian with graph framelets.

For convenience, we write the graph gradient for multiple channel signals $\mathbf{F}$ as

$$\nabla_{W}\mathbf{F}([i,j])=\sqrt{\frac{w_{i,j}}{d_{j,j}}}\mathbf{f}_{j}-\sqrt{\frac{w_{i,j}}{d_{i,i}}}\mathbf{f}_{i}.$$

or simply $\nabla\mathbf{F}$ if no confusion. For an undirected graph we have $\nabla_{W}\mathbf{F}([i,j])=-\nabla_{W}\mathbf{F}([j,i])$. With this notation, we can re-write the $p$-Laplacian regularizer in (13) (the element-wise $p$-norm) in the following.

		$\displaystyle\mathcal{S}_{p}(\mathbf{F})=\frac{1}{2}\sum_{(v_{i},v_{j})\in\mathcal{E}}\left\|\nabla_{W}\mathbf{F}([i,j])\right\|^{p}$		(20)
	$\displaystyle=$	$\displaystyle\frac{1}{2}\sum_{v_{i}\in\mathcal{V}}\left[\left(\sum_{v_{j}\sim v_{i}}\left\|\nabla_{W}\mathbf{F}([i,j])\right\|^{p}\right)^{\frac{1}{p}}\right]^{p}=\frac{1}{2}\sum_{v_{i}\in\mathcal{V}}\|\nabla_{W}\mathbf{F}(v_{i})\|_{p}^{p}$

where $v_{j}\sim v_{i}$ stands for the node $v_{j}$ that is connected to node $v_{i}$ and $\nabla_{W}\mathbf{F}(v_{i})=\left(\nabla_{W}\mathbf{F}([i,j])\right)_{v_{j}:(v_{i},v_{j})\in\mathcal{E}}$ is the node gradient vector for each node $v_{i}$ and $\|\cdot\|_{p}$ is the vector $p$-norm. In fact, $\|\nabla_{W}\mathbf{F}(v_{i})\|_{p}$ measures the variation of $\mathbf{F}$ in the neighbourhood of each node $v_{i}$. Next, we generalize the regularizer by considering any positive convex function $\phi$ as

$\displaystyle\mathcal{S}^{\phi}_{p}(\mathbf{F})=\frac{1}{2}\sum_{v_{i}\in\mathcal{V}}\phi(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p}).$

(21)

It is clear when $\phi(\xi)=\xi^{p}$, we recover the $p$-Laplacian regularizer. In image processing field, several penalty functions have been proposed in the literature. For example, $\phi(\xi)=\xi^{2}$ is known in the context of Tikhonov regularization (Li et al, 2022; Belkin et al, 2004; Assis et al, 2021). When $\phi(\xi)=\xi$ (i.e. $p=1$), it is the classic total variation regularization. When $\phi(\xi)=\sqrt{\xi^{2}+\epsilon^{2}}-\epsilon$, it is referred as the regularized total variation. An example work can be found in the work (Cheng and Wang, 2020). The case of $\phi(\xi)=r^{2}\log(1+\xi^{2}/r^{2})$ corresponds to the non linear diffusion (Oka and Yamada, 2023).

In pL-UFG and pL-fUFG, we use $\mathbf{Y}$ to denote $\mathcal{W}^{T}\textrm{diag}(\theta)\mathcal{W}\mathbf{X}$ and     $\textrm{diag}(\theta_{r,j})\mathcal{W}_{r,j}\mathbf{X}$ respectively, as then we propose the generalized $p$-Laplacian regularization models as

$\displaystyle\mathbf{F}=\operatorname*{arg\,min}_{\mathbf{F}}\mathcal{S}^{\phi}_{p}(\mathbf{F})+\mu\|\mathbf{F}-\mathbf{Y}\|^{2}_{F}.$

(22)

We derive an iterative algorithm for solving the generalized $p$-Laplacian regularization problem (22), presented as the following theorem.

Here $\boldsymbol{\alpha}^{(k)}$, $\boldsymbol{\beta}^{(k)}$ and $\mathbf{M}^{(k)}$ are calculated according to the current features $\mathbf{F}^{(k)}$. When $\phi(\xi)=\xi^{p}$, from Eq. (23) we can have the algorithm for solving pL-UFG and pL-fUFG.

Proof: Define $\mathcal{L}^{\phi}_{p}(\mathbf{F})$ as the objective function in Eq. (22), consider a node feature $\mathbf{f}_{i}$ in $\mathbf{F}$, then we have

		$\displaystyle\frac{\partial\mathcal{L}^{\phi}_{p}(\mathbf{F})}{\partial\mathbf{f}_{i}}=2\mu(\mathbf{f}_{i}\!-\!\mathbf{y}_{i})\!+\!\frac{1}{2}\sum_{v_{k}\in\mathcal{V}}\frac{\partial}{\partial\mathbf{f}_{i}}\phi(\|\nabla_{W}\mathbf{F}(v_{k})\|_{p})$
	$\displaystyle=$	$\displaystyle 2\mu(\mathbf{f}_{i}-\mathbf{y}_{i})+\frac{1}{2}\frac{\partial}{\partial\mathbf{f}_{i}}\phi(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})+\frac{1}{2}\sum_{v_{j}\sim v_{i}}\frac{\partial}{\partial\mathbf{f}_{i}}\phi(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})$
	$\displaystyle=$	$\displaystyle 2\mu(\mathbf{f}_{i}-\mathbf{y}_{i})+\frac{1}{2}\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})\frac{1}{p}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})^{1-p}\frac{\partial}{\partial\mathbf{f}_{i}}\!\left(\!\sum_{v_{j}\sim v_{i}}\|\nabla_{i,j}(w,\mathbf{f})\|^{p}\!\right)\!$
		$\displaystyle+\!\frac{1}{2p}\sum_{v_{j}\sim v_{i}}\frac{\partial}{\partial\mathbf{f}_{i}}\|\nabla_{j,i}(w,\mathbf{f})\|^{p}\cdot\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})^{1-p}$
	$\displaystyle=$	$\displaystyle 2\mu(\mathbf{f}_{i}-\mathbf{y}_{i})+\frac{1}{2p}\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})^{1-p}\cdot$
		$\displaystyle\quad\quad\left(\!\sum_{v_{j}\sim v_{i}}p\|\nabla_{W}\mathbf{F}([i,j])\|^{p-2}\sqrt{\frac{w_{i,j}}{d_{i,i}}}(-\nabla_{W}\mathbf{F}([i,j]))\!\right)$
		$\displaystyle+\frac{1}{2}\sum_{v_{j}\sim v_{i}}\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})^{1-p}\|\nabla_{W}\mathbf{F}([i,j])\|^{p-2}\sqrt{\frac{w_{i,j}}{d_{i,i}}}\nabla_{W}\mathbf{F}([j,i])$
	$\displaystyle=$	$\displaystyle 2\mu(\mathbf{f}_{i}\!-\!\mathbf{y}_{i})\!+\!\sum_{v_{j}\sim v_{i}}\!\frac{1}{2}\|\nabla_{W}\mathbf{F}([i,j])\|^{p-2}\!\!\sqrt{\frac{w_{i,j}}{d_{i,i}}}\nabla_{W}\mathbf{F}([j,i])\cdot$
		$\displaystyle\quad\quad\left[\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{i})\|^{p-1}_{p}}+\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{j})\|^{p-1}_{p}}\right]$

Setting $\frac{\partial\mathcal{L}^{\phi}_{p}(\mathbf{F})}{\partial\mathbf{f}_{i}}=0$ gives the following first order condition:

		$\displaystyle\bigg{(}\sum_{v_{j}\sim v_{i}}\frac{1}{2}\left[\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{i})\|^{p-1}_{p}}+\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{j})\|^{p-1}_{p}}\right]\cdot\|\nabla_{i,j}(w,\mathbf{f})\|^{p-2}\frac{w_{i,j}}{d_{i,i}}+2\mu\bigg{)}\mathbf{f}_{i}$
	$\displaystyle=$	$\displaystyle\sum_{v_{j}\sim v_{i}}\frac{1}{2}\left[\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{i})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{i})\|^{p-1}_{p}}+\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}(v_{j})\|_{p})}{\|\nabla_{W}\mathbf{F}(v_{j})\|^{p-1}_{p}}\right]\cdot\|\nabla_{W}\mathbf{F}([i,j])\|^{p-2}\frac{w_{i,j}}{\sqrt{d_{i,i}d_{j,j}}}\mathbf{f}_{j}+2\mu\mathbf{y}_{i}.$

This equation defines the message passing on each node $v_{i}$ from its neighbour nodes $v_{j}$. With the definition of $M_{i,j}$, $\alpha_{ii}$ and $\beta_{ii}$, it can be turned into the iterative formula in Eq. (23). This completes the proof.

In this section, we present some theoretical support on how the $p$-Laplacian regularizer interact with the framelets in the model. Specifically, we show that the proposed algorithm in Eq. (23) provides a sequence of ($p$-Laplacian-based) spectral graph convolutions and diagonal scaling of the node features over the framelet filtered graph signals.This is indicated by the following analysis.

First considering the iteration Eq. (23) with initial values $\mathbf{F}^{(0)}=\mathbf{Y}=\mathcal{W}^{T}\textrm{diag}(\theta)\mathcal{W}\mathbf{X}$ or $\textrm{diag}(\theta_{r,j})\mathcal{W}_{r,j}\mathbf{X}$ without loss of generality³³3In our practical algorithm we choose $\mathbf{F}^{(0)}=0$, this will gives $\mathbf{F}^{(1)}=\beta^{(0)}\mathbf{Y}$ with a diagonal matrix $\beta^{(0)}$., we have

		$\displaystyle\mathbf{F}^{(K)}=\boldsymbol{\alpha}^{(K\!-\!1)}\mathbf{D}^{-1/2}\mathbf{M}^{(K-1)}\mathbf{D}^{-1/2}\mathbf{F}^{(K\!-\!1)}\!+\!\boldsymbol{\beta}^{(K-1)}\mathbf{Y}$
	$\displaystyle=$	$\displaystyle\boldsymbol{\alpha}^{(K-1)}\widetilde{\mathbf{M}}^{(K-1)}\mathbf{F}^{(K-1)}+\boldsymbol{\beta}^{(K-1)}\mathbf{Y}$
	$\displaystyle=$	$\displaystyle\boldsymbol{\alpha}^{(K-1)}\widetilde{\mathbf{M}}^{(K-1)}\left(\boldsymbol{\alpha}^{(K-2)}\widetilde{\mathbf{M}}^{(K-2)}\mathbf{F}^{(K-2)}+\boldsymbol{\beta}^{(K-2)}\mathbf{Y}\right)+\boldsymbol{\beta}^{(K-1)}\mathbf{Y}$
	$\displaystyle=$	$\displaystyle\boldsymbol{\alpha}^{(K-1)}\widetilde{\mathbf{M}}^{(K-1)}\boldsymbol{\alpha}^{(K-2)}\widetilde{\mathbf{M}}^{(K-2)}\mathbf{F}^{(K-2)}+\boldsymbol{\alpha}^{(K-1)}\widetilde{\mathbf{M}}^{(K-1)}\boldsymbol{\beta}^{(K-2)}\mathbf{Y}+\boldsymbol{\beta}^{(K-1)}\mathbf{Y}$
	$\displaystyle=$	$\displaystyle\left(\prod_{k=0}^{K-1}\boldsymbol{\alpha}^{(k)}\widetilde{\mathbf{M}}^{(k)}\right)\mathbf{Y}+\boldsymbol{\beta}^{(K-1)}\mathbf{Y}+\sum_{k=0}^{K-1}\left(\prod_{l=K-k}^{K-1}\boldsymbol{\alpha}^{(l)}\widetilde{\mathbf{M}}^{(l)}\right)\boldsymbol{\beta}^{(K-k-1)}\mathbf{Y}.$

The result $\mathbf{F}^{(K)}$ depends on the key operations $\boldsymbol{\alpha}^{(k)}\widetilde{\mathbf{M}}^{(k)}$ for $k=0,1,...,K-1$, where

	$\displaystyle\widetilde{M}^{(k)}_{i,j}=\frac{w_{i,j}}{\sqrt{d_{i,i}d_{j,j}}}\left\|\nabla_{W}\mathbf{F}^{(k)}[(i,j)]\right\|^{p-2}\!\!\!\cdot\left[\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}^{(k)}(v_{i})\|_{p})}{\|\nabla_{W}\mathbf{F}^{(k)}(v_{i})\|^{p-1}_{p}}+\frac{\phi^{\prime}(\|\nabla_{W}\mathbf{F}^{(k)}(v_{j})\|_{p})}{\|\nabla_{W}\mathbf{F}^{(k)}(v_{j})\|^{p-1}_{p}}\right]$
	$\displaystyle=\frac{w_{i,j}\cdot w^{(k)}_{i,j}(p)}{\sqrt{d_{i,i}d_{j,j}}}\left\|\nabla_{W}\mathbf{F}^{(k)}[(i,j)]\right\|^{p-2},$		(25)
	$\displaystyle\alpha_{i,i}^{(k)}=1/\left(\sum_{v_{j}\sim v_{i}}{\widetilde{M}^{(k)}_{i,j}}+2\mu\right),$		(26)

where $\widetilde{M}_{ij}$ has absorbed $d_{i,i}$ and $d_{j,j}$ into $M_{ij}$ defined in Theorem 3.1, i.e., $\widetilde{M}_{ij}=M_{i,j}/\sqrt{d_{i,i}d_{j,j}}$.