# Graph Denoising with Framelet Regularizer

Bingxin Zhou, Ruikun Li, Xuebin Zheng, Yu Guang Wang, and Junbin Gao B.Zhou is
with The University of Sydney Business School, The University of Sydney, NSW,
Australia and Institute of Natural Sciences, Shanghai Jiao Tong University,
Shanghai, China.
email<EMAIL_ADDRESS>R.Li, X.Zheng and J.Gao are with The
University of Sydney Business School, The University of Sydney, NSW,
Australia. Y.G.Wang is with Institute of Natural Sciences and School of
Mathematical Sciences, Key Laboratory of Scientific and Engineering Computing
of Ministry of Education, and Center for Mathematics of Artificial
Intelligence Institute, Shanghai Jiao Tong University, Shanghai, China;
Shanghai Artificial Intelligence Laboratory, Shanghai, China; and School of
Mathematics and Statistics, The University of New South Wales, Australia.
Manuscript received October 30, 2021.

###### Abstract

As graph data collected from the real world is merely noise-free, a practical
representation of graphs should be robust to noise. Existing research usually
focuses on feature smoothing but leaves the geometric structure untouched.
Furthermore, most work takes $\mathbb{L}_{2}$-norm that pursues a global
smoothness, which limits the expressivity of graph neural networks. This paper
tailors regularizers for graph data in terms of both feature and structure
noises, where the objective function is efficiently solved with the
alternating direction method of multipliers (ADMM). The proposed scheme allows
to take multiple layers without the concern of over-smoothing, and it
guarantees convergence to the optimal solutions. Empirical study proves that
our model achieves significantly better performance compared with popular
graph convolutions even when the graph is heavily contaminated.

###### Index Terms:

Graph Convolutional Networks, Constrained Optimization, Framelet Transforms,
Graph Denoising.

##  1 Introduction

Data quality is of crucial importance for a reliable modeling. When non-
sampling errors become too large to neglect, the validity of the subsequent
inference for the estimation result remains of limited confidence. Such
concern has been explored throughout for grid data and panel data [1, 2, 3].
As a special type of data with geometric structure and feature, graph data
suffers from the same concern. For instance, when diagnose COVID-19, panicky
suspect patients (nodes) might overstate their symptoms (feature), while some
patients choose to undercover their track (edges). In the early stage of the
outbreak, there could be a systematic error or unstandardized procedures that
falsely classify ’low/medium/high risky regions’, in which case the edge
connection becomes less informative.

Graph Neural Networks (GNNs) have gained a predominant position in graph
representation learning due to its promising empirical performance on various
applications of machine learning. Recent studies in [4, 5] showed that many
graph convolution layers, such as GCN [6] and GAT [7], smooth the graph signal
by intrinsically performing a ${\mathbb{L}}_{2}$-based graph Laplacian
denoising procedure. While the neighborhood-based local smoothness is desired
for graph communities, such methods usually enforce global smoothness, which
makes inter-cluster vertices indistinguishable, and the learned graph
representation suffers from the well-known issue of _over-smoothing_.
Consequently, the model performance drops drastically before it learns the
complicated relationship of graph geometry. In order to solve this problem,
most researchers either consider spectral methods [8, 9] to include high-pass
information, or introduce bias [10, 11] that balances between the fitness and
smoothness of the models.

While the graph feature smoothing has caught increasing attention, literature
on structure denoising is far from sufficient. Researchers usually focus on
the situation of defending graph adversarial attacks, however, such a scenario
assumes a minor manipulation on the graph by hostile attackers, which does not
universally exist in any case and does not threats all entities of a graph. On
the other hand, a fundamental assumption on any GNNs, especially those popular
message-passing networks, is that relevant vertices (or neighbors) share
similar properties. When the precise connection is contaminated or
indeterminant, the model estimation becomes unreliable.

The lack of exploration on graph smoothing and denoising naturally brings our
first research question: How to measure the smoothness of both feature signal
and structure space in a noisy graph? Instead of over-simplify or even neglect
the graph geometry, the design should be tailor-made for a graph that reflects
its key properties such as connectivity. The measure should also be aware of
pitfalls that exist in training GNNs, such as over-smoothing. In consideration
of the practicability, it’s crucial to ask: How to design a robust mechanism
that is adaptive to any graph convolution layers? In this way, our method can
be sliced to other graph convolutions to enhance their performance.

To answer the questions, we explore regularizers on graphs and propose a
framelet-based Double-Term smoother, namely DoT, that detects and erases the
noises in graph representation. We keep the key properties of graphs by
applying the graph-norm on measuring the noise level. The punishment on each
node is weighted by node degree, so that highly active nodes (with more
neighbors) tolerant a smaller risk to be polluted. For feature noises, we
replace the predominant ${\mathbb{L}}_{2}$ norm with ${\mathbb{L}}_{1}$ norm
so that the sparsity in the transformed domain is desired while the smoothness
is not overly pursued. The geometry noise, on the other hand, is measured by
${\mathbb{L}}_{2,1}$ norm where we concentrate errors on merely connected
nodes so that they are less likely to provoke un-neglectable estimation
errors. We validate our model both theoretically and empirically. We detail
the construction of our optimization target as well as its update rules. Its
convergence analysis is attached to support its performance. We also compare
three ablation models with popular convolution methods and several tasks,
where our proposed method outperforms significantly.

The rest of the paper is arranged as follows. Section 2 discusses relevant
research in literature. Section 3 introduces popular graph convolution designs
from the view of signal smoothing then Section 4 presents graph-related
regularizers. Section 5 formulate the objective function for framelet-based
DoT regularizer. The update rule and convergence analysis are provided under
the protocol of the ADMM algorithm, following which three ablation models are
listed. In Section 6, several real-world examples indicate the model’s
recovery ability on different levels and types of contamination. We further
conclude the work in Section 7.

$\blacktriangledown$latentrepresentationfeaturestructurefidelitycl.feature$\Join$cl.structurePredictorlabel
Figure 1: Conceptual architecture of the proposed denoising scheme. For an
input graph representation (either raw or latent), DoT slices its feature and
structure information to do recursive cleaning. The fidelity function is
extracted from the input noise signal to control the power of the cleaning
scheme. The purified representation is then propagated to the next layer for
label prediction (or other operations).

##  2 Related Work

###  2.1 Graph Neural Networks for representation learning

GNNs have shown great success in dealing with irregular graph-structured data
that traditional deep learning methods such as CNNs fail to manage [12, 13,
14, 15]. The main factor that contributes to their success is that GNNs are
able to learn the structure pattern of the graph while CNNs can only handle
regular grid-like structures.

A broad spectrum of GNN approaches has been proposed for defining the graph
convolutional operation, where the two main groups are spatial-based and
spectral-based methods. Spatial-based methods is dominantly researched and
practiced due to its intuitive characteristics [6, 7, 16, 11, 17]. As an
extension of convolutional neural networks (CNNs), spatial-based linear GNNs
pass messages along the graph Laplacian, which naturally performs local
smoothness on second-order graph differences [5, 18] and can be interpreted as
low-pass filters of graph signals [19]. Spectral methods that convert the raw
signal or features in the vertex domain into the frequency domain, on the
other hand, was paid less attention to [20, 21, 22, 23, 8, 9]. In fact,
spectral-based methods have already been proved to have a solid mathematical
foundation in graph signal processing [24]. Versatile Fourier [22, 6, 21],
wavelet transforms[8, 25] and framelets[9] have also shown their capabilities
in graph representation learning. In addition, with fast transforms being
available in computing strategy, a big concern related to efficiency could be
well resolved.

###  2.2 Graph Signal Smoothing and Denoising

Common to both categories of graph convolutions is the signal smoothing of
graph filters. Spatial methods usually implement polynomial [6, 26], Lanczos
[27] or ARMA [28, 29] filters. From the perspective of graph signal denoising,
an aggregation layer approximates smooth graph signals with pairwise quadratic
local variation minimization [30, 31]. While linear approximation achieves
fast optimization over graph denoising, the abandoned details at each layer
limit the model’s fitting power over complex nonlinear functions, which
contradicts the requirement of an effective deep neural network. In contrast,
stacking multiple network layers aggregate complicated relationships of a
large receptive field, but taking a weighted moving average over too many
neighbors could lose sharp changes of frequency response.

Instead, spectral methods, due to their close connection to signal processing,
possess detail-preserved global smoothness. Conventional spectral transforms
has been widely used for processing non-stationary signals, such as images
compression and restoration [32, 33], speech enhancement [34, 35] and so on.
Compared to the Fourier transform of the time-frequency domain, wavelets
transforms engage between the time and scale domain of signals, which
facilitates a multi-scale or multi-resolution view of the input. The signal
coefficients exhibit sparse distribution while noise coefficients spread
uniformly with a small amplitude. This property allows differentiate signal
and noise with threshold-based approaches [36, 37] or empirical risk
minimization approaches [38, 39].

The same philosophy of sparsifying signal representation is implanted to
smoothing spectral graph convolutions. Compared to spatial-based methods where
the sparsity is usually restricted in vertex domain [40], spectral methods
measure the graph sparsity in transformed domain, e.g., in Fourier [41] or
Framelet domain [42]. Usually ${\mathbb{L}}_{1}$ regularization is suggested
as the sparsity measurement.

ADMM [43] or Split-Bregman [44] are conventional solutions of recursive
optimization. While the authors of [45, 46] show that the the feed-forward
propagation of neural networks can be interpreted as optimization iterations,
solving the optimal graph representation over diverse prior design becomes
possible, such as sparsity over input signals [40] or first-order difference
of input [47].

##  3 Graphs Convolutions and ${\mathbb{L}}_{2}$ Smoothers

This section discusses the equivalent expression of graph convolution from the
perspective of graph signal denoising. We consider an undirected graph
${\mathcal{G}}=({\mathbb{V}},\mathbb{E},{\bm{X}})$ with $N=|{\mathbb{V}}|$
nodes. The edge connection is described by an adjacency matrix
${\bm{A}}\in\mathbb{R}^{N\times N}$ and the $d$-dimensional node feature is
stored in ${\bm{X}}\in\mathbb{R}^{N\times d}$. Here we neglect the nonlinear
activation function $\sigma(\cdot)$ for simplicity.

###  3.1 Spatial Graph Convolution

A typical way of denoising graph signal is via graph smoothing, where the
feature smoothness can be forced by a modified _Dirichlet energy function_.
While the smoothed representation should be close to the input representation,
one leverages a _fidelity term_ to guarantee this similarity, and this
similarity is usually measured by a convex distance measure, such as the
Euclidean distance (${\mathbb{L}}_{2}$ norm).

###### Theorem 1

Consider a graph smoothing problem on the input representation
${\bm{X}}{\bm{W}}$

$\min_{{\bm{U}}}{\rm
tr}({\bm{U}}^{\top}\tilde{{\bm{L}}}{\bm{U}})+\|{\bm{U}}-{\bm{X}}{\bm{W}}\|^{2}_{2},$
(1)

where the smoothness of graph is measured by the Dirichlet energy that ${\rm
tr}({\bm{U}}^{\top}\tilde{{\bm{L}}}{\bm{U}})=\frac{1}{2}\sum_{i,j=1}^{N}{\bm{A}}_{ij}(\frac{{\bm{X}}_{i}}{\sqrt{{\bm{D}}_{ii}}}-\frac{{\bm{X}}_{j}}{\sqrt{{\bm{D}}_{jj}}})^{2}$.
Here the
$\tilde{{\bm{L}}}:={\bm{I}}-{\bm{D}}^{-\frac{1}{2}}{\bm{L}}{\bm{D}}^{-\frac{1}{2}}={\bm{I}}-\tilde{{\bm{A}}}$
denotes the normalized graph Laplacian with respect to the adjacency matrix
${\bm{A}}$ or the degree matrix of the graph ${\bm{D}}$. The first-order
approximation to this problem’s optimal solution is
${\bm{U}}:=\tilde{{\bm{A}}}{\bm{X}}{\bm{W}}$.

###### Proof 1

To obtain the optimal solution to the objective function, we calculate the
partial derivatives

$\displaystyle\frac{\partial{\rm
tr}({\bm{U}}^{\top}\tilde{{\bm{L}}}{\bm{U}})+\|{\bm{U}}-{\bm{X}}{\bm{W}}\|^{2}_{2}}{\partial{\bm{U}}}=0,$
$\displaystyle\Rightarrow\;\;$ $\displaystyle
2\tilde{{\bm{L}}}{\bm{U}}+2({\bm{U}}-{\bm{X}}{\bm{W}})=0,$
$\displaystyle\Rightarrow\;\;$
$\displaystyle({\bm{I}}+\tilde{{\bm{L}}}){\bm{U}}={\bm{X}}{\bm{W}},$
$\displaystyle\Rightarrow\;\;$
$\displaystyle{\bm{U}}=({\bm{I}}+\tilde{{\bm{L}}})^{-1}{\bm{X}}{\bm{W}}.$

The first-order Taylor expansion gives that
$({\bm{I}}+\tilde{{\bm{L}}})^{-1}\approx{\bm{I}}-\tilde{{\bm{L}}}$. By
Definition, ${\bm{L}}={\bm{D}}-{\bm{A}}$ and
$\tilde{{\bm{L}}}={\bm{I}}-\tilde{{\bm{A}}}$, which suggests the first-order
approximation of the optimal result
${\bm{U}}\approx\tilde{{\bm{A}}}{\bm{X}}{\bm{W}}$.

This approximation in graph representation learning is known as graph
convolutional networks [6], which is a popular message-passing scheme that
averages one-hop neighborhood representation with normalized adjacency matrix
$\tilde{{\bm{A}}}$. The ${\bm{W}}$ is a learnable weight matrix that embeds
the high-dimensional raw input ${\bm{X}}$ to a lower dimension. There are many
other aggregation strategies by different spatial-based GNNs under the
message-passing [48] framework, such as GAT [7], APPNP [10] and SGC [49].
Though, most of them implicitly achieves this optimization objective. Authors
of [19] explored the equivalence of GCN’s propagation process to low-pass
filters that smooths the graph signal by brutally filtering out all the
detailed information. Other recent research [5, 40] identified that many other
message passing schemes, although designing different propagation rules, share
the same construction logic of the objective function in terms of denoising
graph signals under the smoothness and fitness constraints [30], which is
tricky to balance.

###  3.2 Spectral Graph Convolution

Compared to spatial counterparts, _Spectral graph convolution_ s transform
graph signals into frequency domain for further operation, i.e., graph
smoothing.

###### Theorem 2

For a given transform $\bm{\Phi}$, a spectral-based graph convolution smooths
the input signal ${\bm{X}}$ by optimizing

$\min_{{\bm{U}}}{\rm tr}({\bm{C}}^{\top}{\rm
diag}({\bm{\theta}}){\bm{C}})+\|{\bm{U}}-{\bm{X}}\|^{2}_{2},$ (2)

where ${\bm{C}}:=\bm{\Phi}{\bm{U}}$ denotes the transformed coefficients with
respect to the clean signal ${\bm{U}}$ in frequency domain. The optimal
solution to (2) is ${\bm{\theta}}\star{\bm{X}}=\bm{\Phi}^{-1}{\rm
diag}({\bm{\theta}})\bm{\Phi}{\bm{X}}$.

###### Proof 2

We first rewrite the equation by replacing ${\bm{C}}$ with
$\bm{\Phi}{\bm{U}}$, which gives

$\displaystyle\min_{{\bm{U}}}{\rm tr}({\bm{U}}^{\top}\bm{\Phi}^{\top}{\rm
diag}({\bm{\theta}})\bm{\Phi}{\bm{U}})+\|{\bm{U}}-{\bm{X}}\|^{2}_{2}$

Note that the graph Laplacian
${\bm{L}}={\bm{V}}^{\top}\Lambda{\bm{V}}=\bm{\Phi}^{\top}{\rm
diag}({\bm{\theta}})\bm{\Phi}$ where
$\\{({\bm{v}}_{i},\lambda_{i}):{\bm{v}}_{i}\in{\bm{V}},\lambda_{i}\in\Lambda\\}$
are eigenpairs with respect to ${\bm{L}}$. We then have

${\rm tr}({\bm{U}}^{\top}\bm{\Phi}^{\top}{\rm
diag}({\bm{\theta}})\bm{\Phi}{\bm{U}})={\rm
tr}({\bm{U}}^{\top}{\bm{L}}{\bm{U}})+\|{\bm{U}}-{\bm{X}}\|^{2}_{2}.$

Similar to above, the partial derivatives of this function brings the optimal
solution that ${\bm{U}}={\bm{L}}{\bm{X}}=\bm{\Phi}^{\top}{\rm
diag}({\bm{\theta}})\bm{\Phi}{\bm{X}}$, which is identical to the convolution
function (3) below.

The optimal solution to the above smoothing problem is identical to spectral
graph convolution. Formally, a layer of spectral graph convolution reads

${\bm{\theta}}\star{\bm{X}}=\bm{\Phi}^{-1}{\rm
diag}({\bm{\theta}})\bm{\Phi}{\bm{X}},$ (3)

where the basis $\bm{\Phi}$ is closely related to the eigenvectors of graph
Laplacian and is determined by the type of transform, such as Fourier [20],
Wavelet [50, 25] and Framelet [9] transforms. Compared to spatial
convolutions, this genre allows graph signal processing [51] on multiple
channels, which pursues both local and global smoothness.

As the coefficient are graph signals projected to different channels by low-
pass or high-pass filters, it naturally separates the local and global
smoothness procedure, which makes the over-smoothing and loss-of-expressivity
issue less of a concern. That has been said, the recovered graph signal still
to some extend suffers from the global over-smoothing issue from
${\mathbb{L}}_{2}$ regularizer, as it suppresses regional fluctuations. To
solve this issue, we refer to a similar denoising idea adopted in signal
processing. Meanwhile, as neither of the above functions pay explicit
attention to structure information, we formulate connection-related noise
measurements to fill this gap.

##  4 The Noise of Graph Signals

This section formulates the regularization schemes related to graph smoothing
and subspace extraction to answer the first research question. We start from
the conventional analysis model for graph signal denoising and then introduce
regularizers that filter the feature and structure noises.

###  4.1 Analysis-Based Graph Regularization

We assume the observed graph signal ${\bm{x}}={\bm{u}}+\epsilon$ constitutes a
graph function ${\bm{u}}:{\mathbb{V}}\mapsto\mathbb{R}$ and Gaussian white
noise $\epsilon$. A denoising task recovers the true graph signal ${\bm{u}}$
from the observation ${\bm{x}}$. We restrict a sparse representation of graphs
in the transformed domain and follow the _analysis based model_ [52] to
measure the sparsity priors (that is, the coefficient matrix in the
transformed domain is assumed sparse). The unconstrained minimization is set
to

$\displaystyle\min_{{\bm{u}}}\|{\mathcal{D}}{\bm{u}}\|_{1}+H({\bm{u}}),$

where ${\mathcal{D}}\in\mathbb{R}^{m\times N}$ is a linear transform generated
from discrete transforms, such as discrete Fourier transforms or Framelet
transforms. Its sparsity is constrained by the ${\mathbb{L}}_{1}$-norm. The
$H(\cdot)$ is a smooth convex function that measures the data fidelity, a
common choice of which is the empirical risk minimization
$H({\bm{u}})=\|{\bm{B}}{\bm{u}}-{\bm{x}}\|_{2}$ with some linear
transformation ${\bm{B}}$ and observations ${\bm{x}}$.

###  4.2 Penalty Functions Design

As both node and edge noise are assumed in observed data, we design explicit
regularizers individually to reveal the noiseless set of graph signal
${\bm{U}}$ and adjacency matrix ${\bm{Z}}$.

####  4.2.1 Vertex Feature with Spectral Sparsity

We start from node feature denoising with ${\mathbb{L}}_{1}$ wavelet
regularization. The sparsity and similarity of the signal approximation is
measured by

$\min_{{\bm{U}}}\|{\mathcal{D}}{\bm{U}}\|_{1,G}+\frac{1}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G},$
(4)

For a graph signal ${\bm{u}}$, instead of the Euclidean
${\mathbb{L}}_{p}$-norm, we define a graph ${\mathbb{L}}_{p}$-norm
$\|{\bm{u}}\|_{p,G}:=\left(\sum_{i}|{\bm{u}}[i]|^{p}\cdot
d[i]\right)^{\frac{1}{p}}$, where $d[i]$ indicates the degree of the $i$th
node with respect to ${\bm{u}}$. This work considers $p=1,2$. One can
interpret the graph norm as a weighted-sum version of regularizer, where
errors associated to high-degree nodes result in heavier penalties.

####  4.2.2 Edge Connection with Self-Expressiveness

One fundamental assumption in graph convolution design is that connected nodes
are more likely to share common characteristics, and they can be useful to
distinguish one group of nodes from another. A similar idea is adopted in
sparse subspace clustering, where inner-group information flow is described by
self-expressiveness. We thus leverage this idea and assume that every node of
a graph can be written as a sparse linear combination of its neighbor nodes.
For contaminant data with a noise ${\bm{E}}\in\mathbb{R}^{N\times d}$, we
restrict

$\displaystyle\min_{{\bm{Z}},{\bm{E}}}$
$\displaystyle\|{\bm{Z}}\|_{1}+\|{\bm{E}}\|_{2,1,G}$ $\displaystyle s.t.$
$\displaystyle~{}~{}{\bm{U}}={\bm{Z}}{\bm{U}}+{\bm{E}},{\rm
diag}({\bm{Z}})=0.$ (5)

Here the sum error
$\|{\bm{E}}\|_{2,1,G}=\sum_{i=1}^{N}D_{i}\sqrt{\sum_{j=1}^{D_{i}}|E_{i,j}|^{2}}$
is weighted by $D_{i}$, the degree of node $i$. This convex optimization finds
an effective sparse representation of node connectivity
${\bm{Z}}\in\mathbb{R}^{N\times N}$, which can be served as a noiseless graph
adjacency matrix. We regularize the error term ${\bm{E}}$ by a
${\mathbb{L}}_{2,1,G}$-norm to concentrates errors on features of small-degree
nodes.

###  4.3 Discrete Framelet Transform

As discussed earlier, a graph signal measures its noise level by the sparsity
of the transformed graph coefficients. This work takes fast undecimated graph
Framelet transforms [42, 9] for multi-channel fast denoising. Formally, the
Wavelet Frame, or Framelet, is defined by a filter bank
${\bm{\eta}}:=\\{a;b^{(1)},\dots,b^{(K)}\\}$ (with $K$ the number of high-pass
filters) and the eigen-pairs $\\{(\lambda_{j},{\bm{u}}_{j})\\}_{j=1}^{N}$ of
its graph Laplacian ${\bm{L}}$. The $a$ and $b^{(k)}$ are called the _low-
pass_ and _high-pass filters_ of the Framelet transforms, which preserves the
approximation and detail information of the graph signal, respectively. At
_scale level_ $l=1,\ldots,L$, the low-pass and the $k$th ($k=1,\dots,K$) high
pass undecimated Framelet basis at node $p$ are defined as

$\displaystyle\bm{\varphi}_{l,p}(v)$
$\displaystyle=\sum_{\ell=1}^{N}\hat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\overline{{\bm{u}}_{\ell}(p)}{\bm{u}}_{\ell}(v);$
$\displaystyle\bm{\psi}_{l,p}^{k}(v)$
$\displaystyle=\sum_{\ell=1}^{N}\widehat{b^{(k)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\overline{{\bm{u}}_{\ell}(p)}{\bm{u}}_{\ell}(v),$

where $\\{\alpha;\beta^{(1)},\dots,\beta^{(K)}\\}$ are the associated scaling
functions defined by ${\bm{\eta}}$. The corresponding Framelet coefficients
for node $p$ at scale $l$ of a given signal ${\bm{x}}$ are the projections
$\langle\bm{\varphi}_{l,p},{\bm{x}}\rangle$ and
$\langle\bm{\psi}_{l,p}^{k},{\bm{x}}\rangle$. Note that we formulate the
transforms with Haar-type filters [42] with dilation factor $2^{l}$ to allow
efficient transforms.

The computational cost of decomposition (and reconstruction) algorithm of
Framelet transforms can be heavy. We thus consider approximating the scaling
functions $\alpha$ and $\\{\beta^{(1)},\dots,\beta^{(K)}\\}$ by Chebyshev
polynomials up to $m$-order, respectively, denoted by ${\mathcal{T}}^{m}_{0}$
and $\\{{\mathcal{T}}^{m}_{k}\\}_{k=1}^{K}$. For a given level $L$, we define
the full set of Framelet coefficients

$\displaystyle{\mathcal{W}}_{k,1}{\bm{x}}$
$\displaystyle={\mathcal{T}}^{m}_{k}(2^{-H}{\bm{L}}){\bm{x}}$

at $l=1$. When $l=2,\dots,L$, the coefficients

$\displaystyle{\mathcal{W}}_{k,l}{\bm{x}}$
$\displaystyle={{\mathcal{T}}^{m}_{k}}(2^{-H-l}{\bm{L}}){{\mathcal{T}}^{m}_{0}}(2^{-H-l+1}{\bm{L}}){\cdots{\mathcal{T}}^{m}_{0}}(2^{-H}{\bm{L}}){\bm{x}},$

where the dilation scale $H$ satisfies $\lambda_{\max}\leq 2^{H}\pi$. In this
definition, the finest scale is $1/2^{H+L}$ that guarantees
$\lambda_{\ell}/2^{H+L-l}\in(0,\pi)$ for $\ell=1,2,...,N.$

##  5 A Layer of Graph Feature and Structure Denoising

This section constructs the graph denoising scheme concerning the second
research question. The signal sparsity is measured under the undecimated
Framelet transform system that transforms the graph signal to low-pass and
high-passes Framelet coefficients with the $m$-degree Chebyshev polynomial of
${\bm{L}}$ by a set of $N\times N$ orthonormal decomposition operator
$\bm{{\mathcal{W}}}$. The sequence length $(KL+1)$ is determined by the number
of high-pass filters $K$ and scale level $L$.

###  5.1 ADMM Undecimated Framelet Denoising Model

We now present the denoising layer design. Here we interpret the forward
propagation of a GNN layer as an optimization iteration. The input ${\bm{X}}$
is a noisy representation of the raw signal or its hidden representation. Our
target is to recover ${\bm{U}}$ and ${\bm{Z}}$, a clean set of graph feature
and structure representations that is free from pollution. The objective
function and the corresponding update rule is detailed below.

####  5.1.1 Objective Function Design

The objective function that penalizes both node and edge noises reads

$\displaystyle\min_{{\bm{U}},{\bm{Z}},{\bm{E}}}$
$\displaystyle\|{\bm{\nu}}\bm{{\mathcal{W}}}{\bm{U}}\|_{1,G}+\|{\bm{Z}}\|_{1}+\lambda_{1}\|{\bm{E}}\|_{2,1,G}+\frac{\lambda_{2}}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G},$
$\displaystyle s.t.$
$\displaystyle~{}~{}\operatorname{diag}({\bm{Z}})=0,\;{\bm{U}}={\bm{Z}}{\bm{U}}+{\bm{E}},\;{\bm{Z}}{\bm{1}}={\bm{1}}.$
(6)

Here
$\|{\bm{\nu}}\bm{{\mathcal{W}}}{\bm{U}}\|_{1,G}:=\sum_{(k,l)\in{\mathbb{B}}}\nu_{k,l}\|{\mathcal{W}}_{k,l}{\bm{U}}\|_{1,G}$
promotes the sparsity of the Framelet coefficients of the $k$th high-pass
element on the $l$th scale level, where
${\mathbb{B}}:=\\{(0,L)\\}\bigcup\\{(k,l),1\leq k\leq K,1\leq l\leq L\\}$. The
${\bm{\nu}}$ is a set of tuning parameters with respect to tight decomposition
operator ${\mathcal{W}}_{k,l}$ in high pass. The
${\bm{U}}\in\mathbb{R}^{N\times d}$ is the target (noiseless) signal
approximation of $N$ nodes, ${\bm{L}}$ is the associated graph Laplacian, and
${\bm{X}}$ is the noisy input representation.

The above function can be rewritten with a variable splitting strategy. Define
${\bm{Y}}={\bm{Z}}-\text{diag}({\bm{Z}})$,
${\bm{Q}}:=\bm{{\mathcal{W}}}{\bm{U}}$, i.e.,
${\bm{Q}}_{k,l}:={\mathcal{W}}_{k,l}{\bm{U}}\in\mathbb{R}^{N\times d}$, gives

$\displaystyle\min_{{\bm{U}},{\bm{Z}},{\bm{E}}}$
$\displaystyle\|{\bm{\nu}}{\bm{Q}}\|_{1,G}+\|{\bm{Z}}\|_{1}+\lambda_{1}\|{\bm{E}}\|_{2,1,G}+\frac{\lambda_{2}}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G},$
$\displaystyle s.t.$
$\displaystyle~{}~{}{\bm{U}}={\bm{Y}}{\bm{U}}+{\bm{E}},\;{\bm{Y}}{\bm{1}}={\bm{1}},\;{\bm{Y}}={\bm{Z}}-\text{diag}({\bm{Z}}).$
(7)

The associated augmented Lagrangian [53] reads

$\displaystyle{\mathcal{L}}({\bm{U}},{\bm{Z}},{\bm{E}},{\bm{Q}},{\bm{Y}};\Lambda_{1},\Lambda_{2},\Lambda_{3},\Lambda_{4},\mu_{1},\mu_{2},\mu_{3},\mu_{4})$
$\displaystyle=$
$\displaystyle\sum_{k,l}\nu_{k,l}\|{\bm{Q}}_{k,l}\|_{1,G}+\|{\bm{Z}}\|_{1}+\lambda_{1}\|{\bm{E}}\|_{2,1,G}+\frac{\lambda_{2}}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G}$
$\displaystyle+\frac{\mu_{1}}{2}\|{\bm{U}}-{\bm{Y}}{\bm{U}}-{\bm{E}}\|_{2}^{2}+\sum_{k,l}\frac{\mu_{2}}{2}\|{\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}\|_{2}^{2}$
$\displaystyle+\frac{\mu_{3}}{2}\|{\bm{Y}}{\bm{1}}-{\bm{1}}\|_{2}^{2}+\frac{\mu_{4}}{2}\|{\bm{Y}}-{\bm{Z}}+\text{diag}({\bm{Z}})\|^{2}_{2}$
$\displaystyle+\text{tr}\left(\Lambda_{1}^{\top}({\bm{U}}-{\bm{Y}}{\bm{U}}-{\bm{E}})\right)+\sum_{k,l}\text{tr}\left(\Lambda_{2;k,l}^{\top}({\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}})\right)$
$\displaystyle+\Lambda_{3}^{\top}({\bm{Y}}{\bm{1}}-{\bm{1}})+\text{tr}(\Lambda^{\top}_{4}({\bm{Y}}-{\bm{Z}}+\text{diag}({\bm{Z}}))).$
(8)

We call $\Lambda_{1},\Lambda_{2;k,l}\in\mathbb{R}^{N\times
d};\Lambda_{3}\in\mathbb{R}^{N};\Lambda_{4}\in\mathbb{R}^{N\times N}$ the
Lagrangian multipliers. The $\Lambda_{i}$ ($i=1,2,3,4$) can be considered as
the running sum of errors associated with its constraint, and the target is to
choose the optimal $\Lambda_{i}$ that minimize the residual of the $i$th
constraint. The augmented Lagrangian penalty parameters
$\mu_{1},\mu_{2;k,l},\mu_{3},\mu_{4}$ are positive by definition. They stand
for upper bounds over all the adaptive penalty parameters. In practice we
value an identical $\mu_{2}$ to $\mu_{2;k,l}$ for faster matrix inversion in
updating ${\bm{U}}$. We will omit the subscripts for the rest of paper.

###### Remark 1

Conventional splitting methods select $\mu_{i}$s in advance. This paper
leverages adaptive tuning methods that only require a proper initialization on
$\mu^{(0)}_{i}$s.

###### Remark 2

For simplicity, we can define the low-pass ${\bm{\nu}}_{0,L}=0$ and high-
passse ${\bm{\nu}}_{k,l}=4^{-l-1}\nu_{0}$ by initializing $\nu_{0}$ to avoid
the overwhelming fine-tuning work in ${\bm{\nu}}$.

###### Remark 3

We keep the $\text{diag}({\bm{Z}}^{(t)})$ term here for theoretical
completeness. For implementation we remove the term
$\text{diag}({\bm{Z}}^{(t)})$ as it equals zero.

###### Remark 4

With ${\widetilde{\bm{L}_{0}}}:={\bm{I}}-{\bm{Y}}$, we can rewrite
${\bm{U}}-{\bm{Y}}{\bm{U}}-{\bm{E}}={\widetilde{\bm{L}_{0}}}{\bm{U}}-{\bm{E}}$,
and ${\widetilde{\bm{L}_{0}}}$ denotes the normalized Laplacian of noiseless
adjacency matrix.

####  5.1.2 ADMM Update Scheme

The objective function is optimized with the alternating direction method of
multipliers (ADMM) [43]. There involves six (sets of) parameters to update.
This work updates the primal variables
${\bm{U}},{\bm{Z}},{\bm{E}},{\bm{Y}},{\bm{Q}}$ iteratively, but it is possible
to design a parallel update scheme, see [54]. At the $t+1$th iteration, we
update iteratively

  1. 1.

${\bm{U}}^{(t+1)}$

$\displaystyle=$
$\displaystyle\left(\lambda_{2}{\bm{D}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\widetilde{\bm{L}_{0}}}^{(t)}+\mu_{2}^{(t)}{\bm{I}}\right)^{-1}$
$\displaystyle\left(\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\bm{E}}^{(t)}+\lambda_{2}{\bm{D}}{\bm{X}}-{\widetilde{\bm{L}_{0}}}^{(t)\top}\Lambda_{1}^{(t)}\right.$
$\displaystyle\left.+\sum_{k,l}(\mu_{2}^{(t)}{\mathcal{W}}_{k,l}^{\top}{\bm{Q}}^{(t)}_{k,l}+{\mathcal{W}}_{k,l}^{\top}\Lambda_{2;k,l})\right)$
(9)

where the noiseless graph Laplacian ${\widetilde{\bm{L}_{0}}}^{(t)}$ is
updated with ${\bm{Y}}^{(t)}$. The matrix inversion reaches a precise solution
by the linear solver [55] or the Cholesky Factorization [56].

  2. 2.

${\bm{Z}}^{(t+1)}={\bm{R}}-\text{diag}({\bm{R}})$

where
$\displaystyle{\bm{R}}:=\mathcal{T}_{1/\mu^{(t)}_{4}}\left({\bm{Y}}^{(t)}+{\mu^{(t)}_{4}}^{-1}\Lambda^{(t)}_{4}\right).$
(10)

The
$\mathcal{T}_{\eta}(x)=\text{sign}(x)\max\\{|x|-\eta,0\\}=\text{ReLU}(x-\eta)-\text{ReLU}(-x-\eta)$
is a _soft-threshold operator_.

  3. 3.

${\bm{E}}_{i}^{(t+1)}$

$\displaystyle=\mathcal{T}_{1/\mu_{1}^{(t)}}^{i}\left({\widetilde{\bm{L}_{0}}}^{(t+1)}{\bm{U}}^{(t+1)}+{\Lambda_{1}^{(t)}}/{\mu_{1}^{(t)}}\right).$
(11)

The above update rule works on the $i$th row of ${\bm{E}}$. The
$\mathcal{T}_{\eta}^{i}(x)=({x_{i}}/{\|x_{i}\|_{2}})\max\\{\|x_{i}\|_{2}-\eta,0\\}$
denotes a row-wise soft thresholding for group
${\mathbb{L}}_{2}$-regularization.

  4. 4.

${\bm{Y}}^{(t+1)}$

$\displaystyle=\left(\mu^{(t)}_{1}({\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)}){\bm{U}}^{(t+1)\top}+\mu^{(t)}_{3}{\bm{1}}{\bm{1}}^{\top}+\mu^{(t)}_{4}{\bm{Z}}^{(t+1)}\right.$
$\displaystyle\left.+\Lambda^{(t)}_{1}{\bm{U}}^{(t+1)\top}-\Lambda^{(t)}_{3}{\bm{1}}^{\top}-\Lambda^{(t)}_{4}\right)\left({\mu^{(t)}_{4}}^{-1}{\bm{I}}\right.$
(12)
$\displaystyle\left.-{\mu^{(t+1)}_{4}}^{-1}\widetilde{{\bm{U}}}^{(t+1)}\left[\mu^{(t+1)}_{4}{\bm{I}}+\widetilde{{\bm{U}}}^{(t+1)\top}\widetilde{{\bm{U}}}^{(t+1)}\right]^{-1}\widetilde{{\bm{U}}}^{(t+1)\top}\right).$

The
$\widetilde{{\bm{U}}}^{(t+1)}=[\sqrt{\mu^{(t+1)}_{1}}{\bm{U}}^{(t+1)},\sqrt{\mu^{(t)}_{3}}{\bm{1}}]\in\mathbb{R}^{N\times(d+1)}$.
The ${\bm{1}}$ is a $N$ dimensional one vector. We assume $N\gg d$, i.e., the
graph size is much larger than the node feature dimension.

  5. 5.

${\bm{Q}}_{k,l}^{(t+1)}[i,:]$

$\displaystyle=\mathcal{T}_{\nu_{k,l}d_{i}/\mu_{2}^{(t)}}\left({\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}[i,:]-{\mu^{(t)}_{2}}^{-1}\Lambda_{2;k,l}^{(t)}[i,:]\right).$
(13)

The above update rule acts row-wisely where ${\bm{Q}}_{k,l}^{(t+1)}[i,:]$
denotes the $i$th row of ${\bm{Q}}_{k,l}^{(t+1)}$.

  6. 6.

Update $\Lambda^{(t+1)}$s and $\mu^{(t+1)}$s:

$\displaystyle\Lambda_{1}^{(t+1)}$
$\displaystyle=\Lambda_{1}^{(t)}+\mu_{1}^{(t)}({\bm{U}}^{(t+1)}-{\bm{Y}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)});$
$\displaystyle\Lambda_{2;k,l}^{(t+1)}$
$\displaystyle=\Lambda_{2;k,l}^{(t)}+\mu_{2}^{(t)}({\bm{Q}}_{k,l}^{(t+1)}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)});$
$\displaystyle\Lambda_{3}^{(t+1)}$
$\displaystyle=\Lambda_{3}^{(t)}+\mu_{3}^{(t)}({\bm{Y}}^{(t+1)}{\bm{1}}-{\bm{1}});$
$\displaystyle\Lambda_{4}^{(t+1)}$
$\displaystyle=\Lambda_{4}^{(t)}+\mu_{4}^{(t)}\left({\bm{Y}}^{(t+1)}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)})\right);$
$\displaystyle\mu_{i}^{(t+1)}$
$\displaystyle=\min\left(\rho\mu_{i}^{(t)},\mu_{i,\max}\right),\quad
i=1,2,3,4.$ (14)

The pseudo-code is summarized in Algorithm 1.

Input : iteration $t$, scale level $L$, high-pass number $K$

Output : ${\bm{U}}^{(t)},{\bm{Z}}^{(t)}$

1 Initialization: ${\mathbb{B}}=\\{(0,L))\\}\bigcup\\{(k,l),1\leq k\leq
K,1\leq l\leq L\\}$

2 Initialization:
${\bm{Y}}^{(0)},{\bm{Z}}^{(0)},{\bm{U}}^{(0)},{\bm{Q}}^{(0)},{\bm{E}}^{(0)};\Lambda^{(0)}s,\mu^{(0)}s$

3 for  _$k\leftarrow 1$ to $K$_ do

4 for  _$(k,l)\in{\mathbb{B}}$_ do

5 Update
${\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Y}}^{(t)},{\bm{Q}}_{k,l}^{(t)}$
by (1)-(11);

6 Update
$\Lambda_{1}^{(t)},\Lambda_{2;k,l}^{(t)},\Lambda_{3}^{(t)},\Lambda_{4}^{(t)}$
by (14);

7 Update $\mu_{1}^{(t)},\mu_{2}^{(t)},\mu_{3}^{(t)},\mu_{4}^{(t)}$ by (14)

8 end for

9

10 end for

Update loss function ${\mathcal{L}}$ by (5.1.1).

Algorithm 1 DoT for graph denoising

###  5.2 Convergence of the main algorithm

The original optimization problem (5.1.1), in our proposed scheme, is
considered as a component of neural network. In the last section, we leverage
the ADMM scheme to solve the augmented Lagrangian problem (5.1.1) where the
optimal ${\bm{Z}}$ and ${\bm{U}}$ are both outputs of network layers. The
$t$th iteration of the ADMM update scheme, in particular, acts as the $t$th
layer of neural network that updates from
$\\{{\bm{U}}^{(t-1)},{\bm{Z}}^{(t-1)}\\}$ to
$\\{{\bm{U}}^{(t)},{\bm{Z}}^{(t)}\\}$. This section provides the theorem that
guarantees the convergence of ADMM algorithm for (5.1.1).

###### Theorem 3

Let
$\\{\Gamma_{t}=({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)},\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},$
$\Lambda^{(t)}_{3},\Lambda^{(t)}_{4})\\}^{\infty}_{t=1}$ be the sequence
generated by the ADMM scheme (1)-(14) presented in subsection 5.1.2. Under the
assumption that ${\bm{U}}^{(t)}$ is bounded, the sequence $\\{\Gamma_{t}\\}$
satisfies the following properties:

  1. 1.

The generated sequence $\\{\Gamma_{t}\\}_{t=1}^{\infty}$ is bounded.

  2. 2.

The sequence $\\{\Gamma_{t}\\}_{t=1}^{\infty}$ has at least one accumulation
point
$\Gamma_{*}=({\bm{U}}_{*},{\bm{Z}}_{*},{\bm{E}}_{*},{\bm{Q}}_{*},{\bm{Y}}_{*},\Lambda_{1*},\Lambda_{2*},\Lambda_{3*},\Lambda_{4*})$
that satisfies first-order optimality KKT conditions:

$\displaystyle{\bm{U}}_{*}={\bm{Y}}_{*}{\bm{U}}_{*}+{\bm{E}}_{*};$
$\displaystyle\mathcal{Q}_{*}=\mathcal{W}_{*}{\bm{U}}_{*};$
$\displaystyle{\bm{Y}}_{*}\mathbf{1}=\mathbf{1};$
$\displaystyle{\bm{Y}}_{*}={\bm{Z}}_{*}-\text{diag}({\bm{Z}}_{*});$
$\displaystyle\Lambda_{1*}\in\lambda_{1}\partial_{{\bm{E}}}\|{\bm{E}}_{*}\|_{2,1,G};$
$\displaystyle-\Lambda_{2*;k,l}\in\nu_{k,l}\partial_{Q_{k,l}}\|{\bm{Q}}_{*;k,l}\|_{1,G};$
$\displaystyle\Lambda_{1*}{\bm{U}}^{T}_{*}+\Lambda_{3*}\mathbf{1}^{T}-\Lambda_{4*}=0;$
$\displaystyle\Lambda_{4*}\in\partial_{{\bm{Z}}}\|{\bm{Z}}_{*}\|_{1}.$

  3. 3.

$\\{{\bm{U}}^{(t)}\\},\\{{\bm{Z}}^{(t)}\\},\\{{\bm{E}}^{(t)}\\},\\{{\bm{Q}}^{(t)}\\}$
and $\\{{\bm{Y}}^{(t)}\\}$ each is Cauchy sequence, and thus converges to its
critical point.

We detail the full proof in Appendix B for the sake of consciousness. In
practice, the ADMM typically defines an inner iteration number, such as $10$,
in advance. This pre-defined value provides a viable fast convergence for
ADMM.

###  5.3 Ablation Denoising Models

In addition to the proposed all-in-one DoT denoising optimization, we separate
the objective into several sub-tasks and investigate the effectiveness of
these individual designs. We also compare with the TV regularizer that is
conventionally used for signal denoising as well as some graph node denoising
tasks [42, 41]. The corresponding objective function and the update rules for
the three ablation models are detailed below.

####  5.3.1 Node Feature Denoising

We define node feature denoising model similar to (4), and reformulated the
problem with ${\bm{Q}}=\bm{{\mathcal{W}}}{\bm{U}}$:

$\min_{{\bm{U}}}\|{\bm{\nu}}{\bm{Q}}\|_{1,G}+\frac{1}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G},\quad\text{s.t.}\quad{\bm{Q}}=\bm{{\mathcal{W}}}{\bm{U}}.$
(15)

The function design is similar to that of [42], except that we optimize the
objective with ADMM instead of the Split-Bregman [44] method. The associated
augmented Lagrangian problem reads

$\displaystyle\min_{{\bm{U}}}\sum_{k,l}\nu_{k,l}\|{\bm{Q}}_{k,l}\|_{1,G}+\frac{1}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G}$
(16)
$\displaystyle+\sum_{k,l}\frac{\mu_{2}}{2}\|{\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}\|_{2}^{2}+\sum_{k,l}\text{tr}\left(\Lambda_{2;k,l}^{\top}({\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}})\right).$

The above objective function can solved with the following ADMM scheme:

  1. 1.

${\bm{U}}^{(t+1)}$:

$\displaystyle=\left({\bm{D}}+\mu_{2}^{(t)}{\bm{I}}\right)^{-1}\left({\bm{D}}{\bm{X}}+\sum_{k,l}{\mathcal{W}}_{k,l}^{\top}(\mu_{2}^{(t)}{\bm{Q}}^{(t)}_{k,l}+\Lambda_{2;k,l}^{(t)})\right);$

  2. 2.

${\bm{Q}}_{k,l}^{(t+1)}[i,:]$:

$\displaystyle=\mathcal{T}_{\nu_{k,l}d_{i}/\mu_{2}^{(t)}}\left({\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}[i,:]-\frac{1}{\mu^{(t)}_{2}}\Lambda_{2;k,l}^{(t)}[i,:]\right);$

  3. 3.

Update $\Lambda_{2;k,l}^{(t+1)}$ and $\mu_{2}^{(t+1)}$:

$\displaystyle\Lambda_{2;k,l}^{(t+1)}=$
$\displaystyle\Lambda_{2;k,l}^{(t)}+\mu_{2}^{(t)}({\bm{Q}}_{k,l}^{(t+1)}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)});$
$\displaystyle\mu_{2}^{(t+1)}=$
$\displaystyle\min\left(\rho\mu_{2}^{(t)},\mu_{2,\max}\right).$

In step one of updating ${\bm{U}}^{(t+1)}$, as the matrix
$({\bm{D}}+\mu_{2}^{(t)}{\bm{I}})$ is diagonal, its inverse can be calculated
swiftly without any approximation operation. Also, the update rule in step two
is identical to that of (13), and the same batch operation can be applied
here.

####  5.3.2 Edge Connection Denoising

We next consider a hybrid ${\mathbb{L}}_{2,1}$ regularization in the case of
noisy edge connections, and it is again a subset of the full objective
function (5.1.1):

$\displaystyle\min_{{\bm{U}},{\bm{Z}},{\bm{E}}}$
$\displaystyle\|{\bm{Z}}\|_{1}+\lambda_{1}\|{\bm{E}}\|_{2,1,G}+\frac{\lambda_{2}}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G},$
(17) $\displaystyle s.t.$
$\displaystyle~{}~{}\text{diag}({\bm{Z}})=0,\;{\bm{U}}={\bm{Z}}{\bm{U}}+{\bm{E}},\;{\bm{Z}}{\bm{1}}={\bm{1}}.$

With ${\widetilde{\bm{L}_{0}}}={\bm{I}}-{\bm{Y}}$, its augmented Lagrangian
reads

$\displaystyle{\mathcal{L}}({\bm{U}})=$
$\displaystyle\frac{1}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G}+\frac{\mu_{1}}{2}\Bigl{\|}{\widetilde{\bm{L}_{0}}}{\bm{U}}-{\bm{E}}\Bigr{\|}_{2}^{2}$
$\displaystyle+\text{tr}\left(\Lambda_{1}^{\top}({\widetilde{\bm{L}_{0}}}{\bm{U}}-{\bm{E}})\right).$
(18)

By ${\partial{\mathcal{L}}}/{\partial{\bm{U}}}=0$, the update rule for
${\bm{U}}^{(t+1)}$ reads

$\displaystyle{\bm{U}}^{(t+1)}=$
$\displaystyle\left(\lambda_{2}{\bm{D}}+\mu_{1}^{(t)}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\widetilde{\bm{L}_{0}}}^{(t)}\right)^{-1}$
$\displaystyle\left(\lambda_{2}{\bm{D}}{\bm{F}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\bm{E}}^{(t)}-\Lambda_{1}^{(t)}{\widetilde{\bm{L}_{0}}}^{(t)}\right)$
$\displaystyle\approx$
$\displaystyle\left(\frac{1}{\lambda_{2}}{\bm{D}}^{-1}-\frac{1}{\lambda^{2}_{2}}\mu_{1}^{(t)}{\bm{D}}^{-1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\widetilde{\bm{L}_{0}}}^{(t)}{\bm{D}}^{-1}\right)$
$\displaystyle\left(\lambda_{2}{\bm{D}}{\bm{F}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\bm{E}}^{(t)}-{\widetilde{\bm{L}_{0}}}^{(t)\top}\Lambda_{1}^{(t)}\right).$

Similar to before, the matrix inverse can be approximated by a Linear or
Cholesky solver for computational simplicity.

####  5.3.3 TV Regularizer for Node Feature

Aside from acting the ${\mathbb{L}}_{1}$-norm on Framelet coefficients in
(15), the regularizer can be applied directly on the noisy raw signal with the
TV regularizer for node denoising. The TV regularizer is conventionally used
for signal denoising as well as some graph node denoising tasks [42, 41].
Formally, the objective function is defined as

$\displaystyle\min_{{\bm{U}}}\frac{\alpha}{2}{\rm
tr}({\bm{U}}^{\top}{\bm{L}}{\bm{U}})+\frac{1}{2}\|{\bm{U}}-{\bm{F}}\|^{2}_{2,G}.$
(19)

where $\alpha$ is a tunable hyper-parameter, and ${\bm{L}}$ denotes the
(unnormalized) graph Laplacian.

The update rule of ${\bm{U}}^{(t+1)}$ of the above function is rather
straight-forward. By ${\partial{\mathbb{L}}}/{\partial{\bm{U}}}=0$, we have

$\displaystyle{\bm{U}}^{(t+1)}=({\bm{D}}+\alpha{\bm{L}})^{-1}{\bm{D}}{\bm{X}}\approx({\bm{I}}-\alpha{\bm{D}}^{-1}{\bm{L}}){\bm{X}}.$

##  6 Numerical Experiments

This section reports the performance of our proposed method in comparison to
vanilla graph convolutional layers as well as three ablation models. The main
experiment tests on three node classification tasks with different types of
noise (feature, structure, hybrid) testified. All experiments run with PyTorch
on NVIDIA® Tesla V100 GPU with 5,120 CUDA cores and 16GB HBM2 mounted on an
HPC cluster. Experimental code in PyTorch can be found at
https://github.com/bzho3923/GNN_DoT.

###  6.1 Experimental Protocol

We validate in this experiment the robustness of proposed methods against
feature and/or structure perturbations. We poison the graph before the
training process with a black-box scheme, which requires no extra training
details.

####  6.1.1 Dataset and Baseline

We conduct experiments on three benchmark networks with undirected edges:

  * •

Cora [57, 58]: a text classification dataset with $1,433$ bag-of-words
representation of $2,708$ documents, and they are connected with $5,429$
citation links (edges). The target is to classify each document into one of
$7$ topics.

  * •

Citeseer [57, 58]: another citation network with similar constructions as
Cora, but contains $3,327$ nodes with $3,703$-dimensional input features and
$4,732$ edges. The number of class is $6$.

  * •

Wiki-CS [59]: a snapshot of Wikipedia database from August 2019 with $216,123$
among $11,701$ articles of $10$ topics corresponding to branches of computer
science. Each article is expressed by $300$ features from GloVe word
embeddings.

For all three datasets, the training set constitutes $20$ random labeled
samples from each class. The validation and test data are $500$ and $1,000$,
respectively.

The architecture for all models is fixed to two convolution layers filling a
denoising module that overlays multiple denoising layers. We adopt one of the
three graph convolutions:

  * •

GCN [6]: graph convolutional network, a ${\mathbb{L}}_{2}$-based global
smoothing model that is widely used for semi-supervised node classification
tasks.

  * •

GAT [7]: graph attention network, a message-passing neural network with self-
attention aggregator.

  * •

UFG [9]: a spectral graph neural network enhanced by fast undecimated Framelet
transform. The framework uses conventional ReLU activation (UFG-R) for graph
signal processing but can be extra robust to graph perturbation with shrinkage
activation (UFG-S).

The denoising module falls into one of the following choices: no-denoising
layer, DoT layers, and one of the three ablation denoising layers (node-ADMM,
edge-ADMM, and node-TV) that we introduced in Section 5.

TABLE I: Denoising performance against hybrid error. The accuracy (in
percentage) is reported with standard deviation from 10 repetitive runs.

Model | Dataset | Noise-free | Noisy | DoT | node-ADMM | edge-ADMM | node-TV
---|---|---|---|---|---|---|---
GCN | Cora |  $82.23$$\pm 0.42$ |  $54.31$$\pm 0.94$ |  $64.98$$\pm 1.10$ |  $63.84$$\pm 1.61$ |  $53.24$$\pm 0.75$ |  $54.58$$\pm 1.35$
Citeseer |  $72.30$$\pm 0.93$ |  $34.24$$\pm 2.13$ |  $48.16$$\pm 0.99$ |  $45.06$$\pm 1.22$ |  $35.22$$\pm 1.91$ |  $34.71$$\pm 1.21$
Wiki-CS |  $78.82$$\pm 0.38$ |  $62.27$$\pm 0.61$ |  $64.98$$\pm 0.31$ |  $64.69$$\pm 0.75$ |  $62.60$$\pm 0.18$ |  $63.22$$\pm 0.08$
GAT | Cora |  $82.28$$\pm 0.53$ |  $53.25$$\pm 2.88$ |  $64.45$$\pm 1.78$ |  $64.31$$\pm 1.47$ |  $55.51$$\pm 1.49$ |  $53.38$$\pm 3.58$
Citeseer |  $71.26$$\pm 0.60$ |  $36.17$$\pm 2.14$ |  $46.22$$\pm 1.23$ |  $45.49$$\pm 1.47$ |  $33.94$$\pm 2.54$ |  $35.41$$\pm 2.64$
Wiki-CS |  $77.58$$\pm 0.12$ |  $59.22$$\pm 1.62$ |  $61.52$$\pm 0.65$ |  $61.25$$\pm 0.91$ |  $61.07$$\pm 0.56$ |  $60.15$$\pm 0.47$
UFG-R | Cora |  $83.59$$\pm 0.63$ |  $51.83$$\pm 2.23$ |  $63.86$$\pm 0.75$ |  $61.34$$\pm 0.95$ |  $49.28$$\pm 1.05$ |  $56.87$$\pm 2.05$
Citeseer |  $72.52$$\pm 0.44$ |  $34.83$$\pm 1.89$ |  $45.57$$\pm 0.99$ |  $45.44$$\pm 0.46$ |  $34.83$$\pm 0.70$ |  $35.86$$\pm 1.40$
Wiki-CS |  $75.09$$\pm 0.43$ |  $57.06$$\pm 0.81$ |  $58.13$$\pm 0.56$ |  $57.81$$\pm 0.54$ |  $57.67$$\pm 0.57$ |  $57.64$$\pm 0.57$
UFG-S | Cora |  $83.59$$\pm 0.63$ |  $55.23$$\pm 2.87$ |  $61.86$$\pm 1.17$ |  $58.48$$\pm 1.05$ |  $55.84$$\pm 2.39$ |  $56.54$$\pm 3.81$
Citeseer |  $72.52$$\pm 0.44$ |  $32.18$$\pm 3.46$ |  $45.06$$\pm 1.51$ |  $43.06$$\pm 1.47$ |  $34.58$$\pm 0.84$ |  $34.59$$\pm 0.88$
Wiki-CS |  $75.09$$\pm 0.43$ |  $56.57$$\pm 0.51$ |  $56.71$$\pm 0.64$ |  $56.65$$\pm 0.84$ |  $56.92$$\pm 0.68$ |  $57.31$$\pm 0.16$

####  6.1.2 Training Setup

We take data perturbation in pre-processing. As values of embedded features
are binary for the two citation networks and decimal for Wiki-CS, we spread
$25\%$ poison with random binary noise for the former datasets and Gaussian
white noise with 0.25 standard deviation for the latter one. The structure
perturbation on the three binary adjacency matrices is unified to a $25\%$
noise ratio with an unchanged number of edge connections. That is, we firstly
randomly delete $12.5\%$ and then create $14.3\%$ of total links to keep the
number of perturbed edge connection same as original graph.

For model training, we fix the number of the hidden neuron at $16$ and the
dropout ratio at $0.5$ for all the models, with Adam optimizer. The network
architecture for baseline models is set to two graph convolutional layers with
an activation function and a dropout layer in between. The choice of the
activation function follows the suggestions provided by their authors, where
$ReLU$ is applied for GCN and UFGConv-R, $eLU$ is used for GAT, and the
shrinkage activation is selected for UFGConv-S. The denoising layers are
addressed before the second convolutional layer following another dropout
layer. We removed the activation operation after the first convolutional
layer, considering the denoising layers play the role of filtering latent
features, and it is similar to an non-linear activation operation. All methods
activate the output with a softmax for multi-class classification tasks. The
number of epochs equals $200$. The key hyperparameters are fine-tuned with
grid search. In particular, we search the learning rate and ($L_{2}$) weight
decay from $\\{0.0005,0.001,0.005,0.01,0.05\\}$ and the dilation of UFG from
$\\{2,3,4\\}$ for the base convolution. When training with the denoising
layers, we fix the optimal learning rate, weight decay, dilation and
$\mu_{3},\mu_{4}$ ${1}$, and tune $\mu_{1},\mu_{2}$ from $\\{1,3,5,7,9\\}$ ),
and $\nu_{0}$ from $\\{5,10,50,100,500\\}$. Other unspecified hyper-
parameters, such as the maximum number of level $L$ in UFG and the number of
head in GAT, follows the default values suggested by PyTorch Geometric.

(a)

(b)

(c)

(d)

(e)

Figure 2: Learning curve with UFG-S convolution on Citeseer, hybrid noise. We
report both loss (top row) and accuracy (bottom row) curves on validation
(left column) and test (right column) sets.

###  6.2 Results Analysis

####  6.2.1 Case 1: Hybrid Perturbation

We report in Table I the average accuracy of all methods over $10$ repetitive
runs. Our DoT outperforms ablation methods in most of the scenarios and
improves up to $41\%$ against the baseline scores (in gray). The only
exception is when denoising Wiki-CS with UFG-S, where DoT is beaten by the TV
regularizer with a slightly higher accuracy.

The superior performance of DoT is mainly driven by the effective node
denoising module. In most cases as reported in the table, the node-ADMM method
achieves the second-highest scores while the performance does not distinct
from DoT’s. In opposite, the effect of edge denoising is less significant.
Most scores in the edge-ADMM column have minor improvements over the
baselines, not mention some of them actually have a hurtful of accuracy
because of the absence of node-ADMM module. One possible explanation of this
phenomenon is that the edge connections of graphs in a graph convolution
merely directly influence the expressivity of graph embedding, especially when
no targeted hostile attack is experienced. Our experiment designs non-targeted
data pollution in advance, so that after a graph convolutional layer, the
structure noise is implicitly reflected by feature noises. Consequently, edge
noise no longer has the way to dominate the model performance, and it is
sufficient to denoising on hidden embeddings from the last convolutional
layer, which collects both type of noises.

While the quantitative measurements proves the promising denoising power of
our proposed DoT, we test its robust learning process in Figure 2. In
particular, we present the learning curve of DoT as well as its competitors’
on Citeseer. Compared to the baseline and its ablations, DoT (in orange)
converges rapidly and reaches stationary points at a lower level on the two
loss curves. The variance are also thinner, which demonstrates a stable and
reliable learning process of the trained model.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 3: Box plots of the model performance over node (row $1$) or edge (row
$2$) perturbations. The plots in the three columns reports the average
performance over the three datasets, Cora (column $1$), Citeseer (column $2$)
and Wiki-CS (column $3$), respectively.

####  6.2.2 Case 2: Individual Perturbation

we supplement model performance under extreme feature OR structure
perturbation and investigate the main contribution of each individual
denoising module. The total noise ratio is increased to a high level of $50\%$
so as to magnify the discrepancy of models. For both node and edge noises, the
same distribution that was introduced previously is applied here again. For
the comparison baseline models, we remove UFG-S with shrinkage activation, as
its smoothing power over large noise has been proven in [9]. Appending a
denoising module on top of it thus becomes less necessary.

We present the model comparison in Figure 3. The results suggest exact
contributors of different noise types. In the first settings of $50\%$ node
noises (the first row), node-ADMM outperforms other methods significantly,
following DoT that mainly obliterating node noises. The disparity comes from
the edge-denoising module of DoT that is designed to minimize the norm of
variance of graph topology information. When there exists little structure
noise in the raw input, the extra penalty term could make a negative
contribution to the overall performance in the slightest way. The edge-ADMM
method barely maintains a similar performance to the baseline methods, which
indicates neither of them identifies the real pains of input. The node-TV
method performs surprisingly badly, though. It could be caused by the trace
regularizer in the objective function design, which does not eliminate or
relieve the long-existing problem of over-smoothing, like many other graph
convolutional layers.

In the case of edge-denoising tasks, DoT and edge-ADMM perform as well as
expectation. The node-ADMM and node-TV methods, although the same strategy of
purifying hidden node representations has proven its effectiveness under a
small level of perturbation, their performance is worse than the two methods
that employ explicit edge denoising modules. This observation illustrates the
importance of designing a specific edge denoising module, especially in
extreme scenarios. Whenever it is difficult to justify whether the structure
noise is too large to be implicitly removed, it is always suggested to include
an edge module.

##  7 Conclusion

This paper introduces a robust optimization scheme for denoising graph
features and structure. The graph noise level is measured by a double-term
regularizer, which pushes Framelet coefficients of the node representation to
be sparse and the adjacency matrix to concentrate its uncertainty on small
communities with limited influence. The optimal graph representation is
reached by a sequence of feed-forward propagation, where each layer functions
as an optimization iteration. The update converges to the exact solution of
robust graph representation, which removes irregular graph volatility but
keeps expressivity. We provide comprehensive theoretical support and extensive
empirical validation to prove the effectiveness of our denoising method in
terms of node and/or edge contamination.

We shall emphasize that in experiments we fix the ADMM iteration up to $10$ to
save the computational cost, which is lower than the typical number of $30$ to
$40$. In other words, our model preserves the potential to achieve higher
scores.

## References

  * [1] J. Westerlund, “Testing for error correction in panel data,” _Oxford Bulletin of Economics and Statistics_ , vol. 69, no. 6, pp. 709–748, 2007.
  * [2] A. C. To, J. R. Moore, and S. D. Glaser, “Wavelet denoising techniques with applications to experimental geophysical data,” _Signal Processing_ , vol. 89, no. 2, pp. 144–160, 2009.
  * [3] B. Goyal, A. Dogra, S. Agrawal, B. S. Sohi, and A. Sharma, “Image denoising review: from classical to state-of-the-art approaches,” _Information Fusion_ , vol. 55, pp. 220–244, 2020.
  * [4] W. Jin, Y. Ma, X. Liu, X. Tang, S. Wang, and J. Tang, “Graph structure learning for robust graph neural networks,” in _ACM SIGKDD_ , 2020, pp. 66–74.
  * [5] M. Zhu, X. Wang, C. Shi, H. Ji, and P. Cui, “Interpreting and unifying graph neural networks with an optimization framework,” in _Proceedings of the Web Conference 2021_ , 2021, pp. 1215–1226.
  * [6] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _ICLR_ , 2017.
  * [7] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio, “Graph attention networks,” in _ICLR_ , 2018.
  * [8] B. Xu, H. Shen, Q. Cao, Y. Qiu, and X. Cheng, “Graph wavelet neural network,” in _ICLR_ , 2018.
  * [9] X. Zheng, B. Zhou, J. Gao, Y. G. Wang, P. Liò, M. Li, and G. Montúfar, “How framelets enhance graph neural networks,” in _ICML_ , 2021.
  * [10] J. Klicpera, A. Bojchevski, and S. Günnemann, “Predict then propagate: graph neural networks meet personalized pagerank,” in _ICLR_ , 2018.
  * [11] K. Xu, C. Li, Y. Tian, T. Sonobe, K.-i. Kawarabayashi, and S. Jegelka, “Representation learning on graphs with jumping knowledge networks,” in _ICML_ , 2018.
  * [12] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, “Geometric deep learning: going beyond Euclidean data,” _IEEE Signal Processing Magazine_ , vol. 34, no. 4, pp. 18–42, July 2017.
  * [13] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, “Graph neural networks: a review of methods and applications,” _AI Open_ , vol. 1, pp. 57–81, 2020.
  * [14] Z. Zhang, P. Cui, and W. Zhu, “Deep learning on graphs: a survey,” _IEEE Transactions on Knowledge and Data Engineering_ , 2020.
  * [15] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip, “A comprehensive survey on graph neural networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 1, pp. 4–24, 2020.
  * [16] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful are graph neural networks?” in _ICLR_ , 2018.
  * [17] X. Jiang, R. Zhu, S. Li, and P. Ji, “Co-embedding of nodes and edges with graph neural networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , pp. 1–1, 2020.
  * [18] S. Chen, Y. C. Eldar, and L. Zhao, “Graph unrolling networks: interpretable neural networks for graph signal denoising,” _IEEE Transactions on Signal Processing_ , 2021.
  * [19] H. Nt and T. Maehara, “Revisiting graph neural networks: all we have is low-pass filters,” _arXiv:1905.09550_ , 2019.
  * [20] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” in _ICLR_ , 2014.
  * [21] M. Henaff, J. Bruna, and Y. LeCun, “Deep convolutional networks on graph-structured data,” _arXiv:1506.05163_ , 2015.
  * [22] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _NIPS_ , 2016\.
  * [23] R. Levie, F. Monti, X. Bresson, and M. M. Bronstein, “Cayleynets: graph convolutional neural networks with complex rational spectral filters,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 1, pp. 97–109, 2018\.
  * [24] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Processing Magazine_ , vol. 30, no. 3, pp. 83–98, 2013.
  * [25] X. Zheng, B. Zhou, M. Li, Y. G. Wang, and J. Gao, “MathNet: Haar-like wavelet multiresolution-analysis for graph representation and learning,” _arXiv:2007.11202_ , 2020.
  * [26] N. Tremblay, P. Gonçalves, and P. Borgnat, “Design of graph filters and filterbanks,” in _Cooperative and Graph Signal Processing_. Elsevier, 2018, pp. 299–324.
  * [27] R. Liao, Z. Zhao, R. Urtasun, and R. Zemel, “Lanczosnet: multi-scale deep graph convolutional networks,” in _ICLR_ , 2019.
  * [28] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Autoregressive moving average graph filtering,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 2, pp. 274–288, 2016.
  * [29] T. Schnake, O. Eberle, J. Lederer, S. Nakajima, K. T. Schutt, K.-R. Mueller, and G. Montavon, “Higher-order explanations of graph neural networks via relevant walks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , pp. 1–1, 2021.
  * [30] D. Zhou, O. Bousquet, T. N. Lal, J. Weston, and B. Schölkopf, “Learning with local and global consistency,” in _NIPS_ , 2004.
  * [31] G. Fu, Y. Hou, J. Zhang, K. Ma, B. F. Kamhoua, and J. Cheng, “Understanding graph neural networks from graph signal denoising perspectives,” _arXiv:2006.04386_ , 2020.
  * [32] S. V. Vaseghi, _Advanced Digital Signal Processing and Noise Reduction_. John Wiley & Sons, 2008.
  * [33] S. Parrilli, M. Poderico, C. V. Angelino, and L. Verdoliva, “A nonlocal sar image denoising algorithm based on llmmse wavelet shrinkage,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 50, no. 2, pp. 606–616, 2011\.
  * [34] M. Bahoura and J. Rouat, “Wavelet speech enhancement based on time–scale adaptation,” _Speech Communication_ , vol. 48, no. 12, pp. 1620–1637, 2006\.
  * [35] P. C. Loizou, _Speech Enhancement: Theory and Practice_. CRC press, 2007.
  * [36] S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” _IEEE Transactions on Image Processing_ , vol. 9, no. 9, pp. 1532–1546, 2000.
  * [37] M. Aminghafari, N. Cheze, and J.-M. Poggi, “Multivariate denoising using wavelets and principal component analysis,” _Computational Statistics & Data Analysis_, vol. 50, no. 9, pp. 2381–2398, 2006.
  * [38] S. Ramani, T. Blu, and M. Unser, “Monte-carlo sure: a black-box optimization of regularization parameters for general denoising algorithms,” _IEEE Transactions on Image Processing_ , vol. 17, no. 9, pp. 1540–1554, 2008.
  * [39] Y. Ding and I. W. Selesnick, “Artifact-free wavelet denoising: non-convex sparse regularization, convex optimization,” _IEEE Signal Processing Letters_ , vol. 22, no. 9, pp. 1364–1368, 2015.
  * [40] X. Liu, W. Jin, Y. Ma, Y. Li, H. Liu, Y. Wang, M. Yan, and J. Tang, “Elastic graph neural networks,” in _ICML_ , 2021.
  * [41] Z. Chen, T. Ma, Z. Jin, Y. Song, and Y. Wang, “BiGCN: a bi-directional low-pass filtering graph neural network,” _arXiv:2101.05519_ , 2021.
  * [42] B. Dong, “Sparse representation on graphs by tight wavelet frames and applications,” _Applied and Computational Harmonic Analysis_ , vol. 42, no. 3, pp. 452–479, 2017.
  * [43] D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear variational problems via finite element approximation,” _Computers & Mathematics with Applications_, vol. 2, no. 1, pp. 17–40, 1976.
  * [44] T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” _SIAM journal on imaging sciences_ , vol. 2, no. 2, pp. 323–343, 2009.
  * [45] X. Xie, J. Wu, G. Liu, Z. Zhong, and Z. Lin, “Differentiable linearized ADMM,” in _ICML_. PMLR, 2019, pp. 6902–6911.
  * [46] J. Li, M. Xiao, C. Fang, Y. Dai, C. Xu, and Z. Lin, “Training neural networks by lifted proximal operator machines,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2020.
  * [47] Y.-X. Wang, J. Sharpnack, A. Smola, and R. Tibshirani, “Trend filtering on graphs,” in _Artificial Intelligence and Statistics_. PMLR, 2015, pp. 1042–1050.
  * [48] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, “Neural message passing for quantum chemistry,” in _ICML_ , 2017.
  * [49] F. Wu, A. Souza, T. Zhang, C. Fifty, T. Yu, and K. Weinberger, “Simplifying graph convolutional networks,” in _ICML_ , 2019.
  * [50] B. Xu, H. Shen, Q. Cao, Y. Qiu, and X. Cheng, “Graph wavelet neural network,” in _ICLR_ , 2019.
  * [51] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, “Graph signal processing: overview, challenges, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [52] J.-F. Cai, S. Osher, and Z. Shen, “Split bregman methods and frame based image restoration,” _Multiscale Modeling & Simulation_, vol. 8, no. 2, pp. 337–369, 2010.
  * [53] J. Nocedal and S. Wright, _Numerical Optimization_. Springer Science & Business Media, 2006.
  * [54] W. Deng, M.-J. Lai, Z. Peng, and W. Yin, “Parallel multi-block admm with o (1/k) convergence,” _Journal of Scientific Computing_ , vol. 71, no. 2, pp. 712–736, 2017.
  * [55] G. H. Golub, _Matrix Computations_ , 3rd ed., ser. Johns Hopkins studies in the mathematical sciences. Baltimore: Johns Hopkins University Press, 1996.
  * [56] C. N. Haddad, _Encyclopedia of Optimization_. Springer US, 2009, ch. Cholesky Factorization.
  * [57] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad, “Collective classification in network data,” _AI Magazine_ , vol. 29, no. 3, pp. 93–93, 2008.
  * [58] Z. Yang, W. Cohen, and R. Salakhudinov, “Revisiting semi-supervised learning with graph embeddings,” in _ICML_ , 2016.
  * [59] P. Mernyei and C. Cangea, “Wiki-cs: a wikipedia-based benchmark for graph neural networks,” _arXiv:2007.02901_ , 2020.
  * [60] K. S. Miller, “On the inverse of the sum of matrices,” _Mathematics Magazine_ , vol. 54, no. 2, pp. 67–72, 1981.
  * [61] M. A. Woodbury, “Inverting modified matrices,” _Memorandum Report_ , vol. 42, no. 106, p. 336, 1950.
  * [62] Z. Lin, R. Liu, and Z. Su, “Linearized alternating direction method with adaptive penalty for low-rank representation,” in _NIPS_ , 2011.

##  Appendix A DoT Update Scheme in Detail

This section derives the calculation details of the update rules for DoT in
Section 5.

###  A.1 Update ${\bm{U}}^{(t+1)}$

We start from updating ${\bm{U}}^{(t+1)}$ with
${\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}_{k,l}^{(t)}$, $\Lambda_{k,l}^{(t)}$
and $\mu^{(t)}_{k,l}$s. We omit the constant term in (5.1.1) and rewrite the
minimization function to

$\displaystyle{\mathcal{L}}({\bm{U}})=\frac{\lambda_{2}}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G}+\frac{\mu^{(t)}_{1}}{2}\Bigl{\|}{\widetilde{\bm{L}_{0}}}^{(t)}{\bm{U}}-{\bm{E}}^{(t)}\Bigr{\|}_{2}^{2}$
$\displaystyle\qquad\qquad+\text{tr}\left(\Lambda_{1}^{(t)\top}({\widetilde{\bm{L}_{0}}}^{(t)}{\bm{U}}-{\bm{E}}^{(t)})\right)$
$\displaystyle\sum_{k,l}\left(\frac{\mu^{(t)}_{2}}{2}\bigl{\|}{\bm{Q}}^{(t)}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}\bigr{\|}_{2}^{2}+\text{tr}\left(\Lambda_{2;k,l}^{(t)\top}({\bm{Q}}^{(t)}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}})\right)\right),$

where ${\widetilde{\bm{L}_{0}}}^{(t)}={\bm{I}}-{\bm{Y}}^{(t)}$. Note that
$\frac{1}{2}\|{\bm{U}}-{\bm{X}}\|^{2}_{2,G}=\frac{1}{2}\text{tr}(({\bm{U}}-{\bm{X}})^{\top}{\bm{D}}({\bm{U}}-{\bm{X}}))$
where ${\bm{D}}=\text{diag}(d_{1},...,d_{N})^{\top}$. We then take
${\partial{\mathcal{L}}}/{\partial{\bm{U}}}=0$, which gives

$\displaystyle{\bm{U}}^{(t+1)}=\left(\lambda_{2}{\bm{D}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\widetilde{\bm{L}_{0}}}^{(t)}+\mu_{2}^{(t)}{\bm{I}}\right)^{-1}$
$\displaystyle\qquad\qquad\left(\lambda_{2}{\bm{D}}{\bm{X}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\bm{E}}^{(t)}-{\widetilde{\bm{L}_{0}}}^{(t)\top}\Lambda_{1}^{(t)}\right.$
$\displaystyle\qquad\qquad\left.+\sum_{k,l}(\mu_{2}^{(t)}{\mathcal{W}}_{k,l}^{\top}{\bm{Q}}^{(t)}_{k,l}+{\mathcal{W}}_{k,l}^{\top}\Lambda_{2;k,l}^{(t)})\right)$
$\displaystyle\approx\left({\bm{I}}-\mu^{(t)}_{2}(\lambda_{2}{\bm{D}}+\mu^{(t)}_{2}{\bm{I}})^{-1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\widetilde{\bm{L}_{0}}}^{(t)}\right)\left(\lambda_{2}{\bm{D}}+\mu^{(t)}_{2}{\bm{I}}\right)^{-1}$
$\displaystyle\quad\left(\lambda_{2}{\bm{D}}{\bm{X}}+\mu^{(t)}_{1}{\widetilde{\bm{L}_{0}}}^{(t)\top}{\bm{E}}^{(t)}-{\widetilde{\bm{L}_{0}}}^{(t)\top}\Lambda_{1}^{(t)}\right.$
$\displaystyle\quad\left.+\sum_{k,l}(\mu_{2}^{(t)}{\mathcal{W}}_{k,l}^{\top}{\bm{Q}}^{(t)}_{k,l}+{\mathcal{W}}_{k,l}^{\top}\Lambda_{2;k,l}^{(t)})\right).$

The above equation contains matrix inversion, where the Maclaurin Series of
first order Taylor expansion [60] is applied for approximation. Alternatively,
one could consider the linear solver [55] or the Cholesky Factorization [56].
in the case when the spectral radius of the inverse matrix is larger than one
and the matrix approximation fails to converge.

###  A.2 Update ${\bm{Z}}^{(t+1)}$

For the third step, we fix ${\bm{U}}^{(t+1)}$, ${\bm{Y}}^{(t+1)}$ and all
other variables from iteration $t$ to update ${\bm{Z}}^{(t+1)}$. The sub-
problem is to minimize the following objective

$\displaystyle{\mathcal{L}}({\bm{Z}})=$
$\displaystyle\|{\bm{Z}}\|_{1}+\frac{\mu^{(t)}_{4}}{2}\|{\bm{Y}}^{(t)}-{\bm{Z}}+\text{diag}({\bm{Z}})\|^{2}$
$\displaystyle+\text{tr}\left(\Lambda^{(t)\top}_{4}\left({\bm{Y}}^{(t)}-{\bm{Z}}+\text{diag}({\bm{Z}})\right)\right).$

That can be solved by the following closed-form solution

$\displaystyle{\bm{Z}}^{(t+1)}=$
$\displaystyle{\bm{R}}-\text{diag}({\bm{R}}),$ $\displaystyle\text{ where
}{\bm{R}}:=$
$\displaystyle\mathcal{T}_{1/\mu^{(t)}_{4}}\left({\bm{Y}}^{(t)}+\frac{1}{\mu^{(t)}_{4}}\Lambda^{(t)}_{4}\right).$

Here the $\mathcal{T}_{\eta}(\cdot)$ is the _soft threshold operator_ defined
as follows:

$\displaystyle\mathcal{T}_{\eta}(x)=\text{sign}(x)\max\\{|x|-\eta,0\\}$

that $\text{ReLU}(x-\eta)-\text{ReLU}(-x-\eta)$.

###  A.3 Update ${\bm{E}}^{(t+1)}$

The associated objective function with respect to ${\bm{E}}$ reads

$\displaystyle{\mathcal{L}}({\bm{E}})=$
$\displaystyle\|{\bm{E}}\|_{2,1,G}+\frac{\mu_{1}^{(t)}}{2}\bigl{\|}{\widetilde{\bm{L}_{0}}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}\bigr{\|}_{2}^{2}$
$\displaystyle+\text{tr}\left(\Lambda_{1}^{(t)\top}\bigl{(}{\widetilde{\bm{L}_{0}}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}\bigr{)}\right)$
$\displaystyle=$
$\displaystyle\|{\bm{E}}\|_{2,1,G}+\frac{\mu_{1}^{(t)}}{2}\left\|{\bm{E}}-\Bigl{(}{\widetilde{\bm{L}_{0}}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{\Lambda_{1}^{(t)}}{\mu_{1}^{(t)}}\Bigr{)}\right\|_{2}^{2},$

where ${\widetilde{\bm{L}_{0}}}^{(t+1)}={\bm{I}}-{\bm{Y}}^{(t+1)}$. The
solution to the $i$th row of ${\bm{E}}^{(t+1)}$ is then

$\displaystyle{\bm{E}}_{i}^{(t+1)}=\mathcal{T}_{1/\mu_{1}^{(t)}}^{i}\left({\widetilde{\bm{L}_{0}}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{\Lambda_{1}^{(t)}}{\mu_{1}^{(t)}}\right).$

Here $\mathcal{T}^{i}$ is a row-wise soft-thresholding for group
${\mathbb{L}}_{2}$-regularization. For the $i$th row of $x$,

$\displaystyle\mathcal{T}_{\eta}^{i}(x)=\frac{x_{i}}{\|x_{i}\|_{2}}\max\\{\|x_{i}\|_{2}-\eta,0\\}.$

###  A.4 Update ${\bm{Y}}^{(t+1)}$

We next use fixed ${\bm{U}}^{(t+1)}$ and all other variables from iteration
$t$ to update ${\bm{Y}}^{(t+1)}$. Similar to the last step, we rewrite the
optimization problem by omitting the constant terms, which gives

$\displaystyle{\mathcal{L}}({\bm{Y}})=$
$\displaystyle\frac{\mu^{(t)}_{1}}{2}\|{\bm{U}}^{(t+1)}-{\bm{Y}}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)}\|_{2}^{2}+\frac{\mu^{(t)}_{3}}{2}\|{\bm{Y}}{\bm{1}}-{\bm{1}}\|_{2}^{2}$
$\displaystyle+\frac{\mu^{(t)}_{4}}{2}\|{\bm{Y}}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)})\|^{2}_{2}$
$\displaystyle+\text{tr}\left(\Lambda_{1}^{(t)\top}({\bm{U}}^{(t+1)}-{\bm{Y}}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t)})\right)$
$\displaystyle+\Lambda_{3}^{(t)\top}({\bm{Y}}{\bm{1}}-{\bm{1}})+\text{tr}\left(\Lambda_{4}^{(t)\top}({\bm{Y}}-{\bm{Z}}^{(t)}+\text{diag}({\bm{Z}}^{(t)}))\right).$

With ${\partial{\mathcal{L}}}/{\partial{\bm{Y}}}=0$, we have

$\displaystyle{\bm{Y}}^{(t+1)}=$
$\displaystyle\left(\mu^{(t)}_{1}({\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)}){\bm{U}}^{(t+1)\top}+\mu^{(t)}_{3}{\bm{1}}{\bm{1}}^{\top}\right.$
$\displaystyle\left.+\mu^{(t)}_{4}{\bm{Z}}^{(t+1)}+\Lambda^{(t)}_{1}{\bm{U}}^{(t+1)\top}-\Lambda^{(t)}_{3}{\bm{1}}^{\top}-\Lambda^{(t)}_{4}\right)$
$\displaystyle\left(\mu^{(t)}_{1}{\bm{U}}^{(t+1)}{\bm{U}}^{(t+1)\top}+\mu^{(t)}_{3}{\bm{1}}{\bm{1}}^{\top}+\mu^{(t)}_{4}{\bm{I}}\right)^{-1}.$

Calculating the above equation involves the inversion of a considerably large
matrix. When $N>>d$, we suggest apply the _Woodbury identity_ [61] to reduce
the $\mathbb{R}^{N\times N}$ matrix to a smaller size of
$\mathbb{R}^{(d+1)\times(d+1)}$. We rewrite

$\displaystyle\mu^{(t)}_{1}{\bm{U}}^{(t+1)}{\bm{U}}^{(t+1)\top}+\mu^{(t)}_{3}{\bm{1}}{\bm{1}}^{\top}+\mu^{(t)}_{4}{\bm{I}}$
$\displaystyle=$
$\displaystyle\widetilde{{\bm{U}}}^{(t+1)}\widetilde{{\bm{U}}}^{(t+1)\top}+\mu^{(t)}_{4}{\bm{I}},$

where
$\widetilde{{\bm{U}}}^{(t+1)}=[\sqrt{\mu^{(t)}_{1}}{\bm{U}}^{(t+1)},\sqrt{\mu^{(t)}_{3}}{\bm{1}}]\in\mathbb{R}^{N\times(d+1)}$.
By the Woodbury identity, we reduce

$\displaystyle(\widetilde{{\bm{U}}}^{(t+1)}\widetilde{{\bm{U}}}^{(t+1)\top}+\mu^{(t)}_{4}{\bm{I}})^{-1}$
$\displaystyle\approx$
$\displaystyle\frac{1}{\mu^{(t)}_{4}}{\bm{I}}-\frac{1}{\mu^{(t)}_{4}}\widetilde{{\bm{U}}}^{(t+1)}\left[\mu^{(t)}_{4}{\bm{I}}+\widetilde{{\bm{U}}}^{(t+1)\top}\widetilde{{\bm{U}}}^{(t+1)}\right]^{-1}\widetilde{{\bm{U}}}^{(t+1)\top},$

where
$\mu^{(t)}_{4}{\bm{I}}+\widetilde{{\bm{U}}}^{(t+1)\top}\widetilde{{\bm{U}}}^{(t+1)}$
is in size $(d+1)\times(d+1)$.

###  A.5 Update ${\bm{Q}}^{(t+1)}$

We first rewrite the objective function (5.1.1) as

$\displaystyle{\mathcal{L}}({\bm{Q}})=\sum_{k,l}\nu_{k,l}\|{\bm{Q}}_{k,l}\|_{1,G}+\frac{\mu_{2}^{(t)}}{2}\left\|{\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}\right\|_{2}^{2}$
$\displaystyle\qquad\qquad+\text{tr}\left(\Lambda_{2;k,l}^{(t)\top}\left({\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}\right)\right)$
$\displaystyle=\sum_{k,l}\nu_{k,l}\|{\bm{Q}}_{k,l}\|_{1,G}+\frac{\mu_{2}^{(t)}}{2}\left\|{\bm{Q}}_{k,l}-\left({\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}-\frac{\Lambda_{2;k,l}^{(t)}}{\mu_{2}^{(t)}}\right)\right\|_{2}^{2}.$

The above formulation suggests a row-wise update of ${\bm{Q}}_{k,l}^{(t+1)}$.
For the $i$th row,

$\displaystyle{\bm{Q}}_{k,l}^{(t+1)}[i,:]=\mathcal{T}_{\nu_{k,l}d_{i}/\mu_{2}^{(t)}}\left({\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}[i,:]-\frac{1}{\mu^{(t)}_{2}}\Lambda_{2;k,l}^{(t)}[i,:]\right).$

###### Remark 5

A batch operation can be considered here in implementation. Denote $d$ the
feature dimension of ${\bm{X}}$, and ${\bm{D}}$ the diagonal of graph degrees.
We define $\Delta$ as a matrix consisting of $d$ repeated columns from
${\bm{D}}$ so that $\Delta_{k,l}=\frac{\nu_{k,l}}{\mu^{(t)}_{2}}\Delta$. With
matrix operation we have
${\bm{Q}}_{k,l}^{(t+1)}=\mathcal{T}_{\Delta_{k,l}}\left({\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}-\frac{1}{\mu^{(t)}_{2}}\Lambda_{2;k,l}^{(t)}\right)$.

###  A.6 Update $\Lambda^{(t+1)}$s

We now update the Lagrangian multipliers with respect to the three
constraints. Fix variables ${\bm{U}}^{(t+1)}$,
${\bm{Q}}_{k,l}^{(t+1)},{\bm{Z}}^{(t+1)},{\bm{E}}^{(t+1)}$ and the parameters
$\mu^{(t)}$s, we update the multiplier $\Lambda_{k,l}^{(t+1)}$s by

$\displaystyle\Lambda_{1}^{(t+1)}$
$\displaystyle=\Lambda_{1}^{(t)}+\mu_{1}^{(t)}({\bm{U}}^{(t+1)}-{\bm{Y}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)})$
$\displaystyle\Lambda_{2;k,l}^{(t+1)}$
$\displaystyle=\Lambda_{2;k,l}^{(t)}+\mu_{2}^{(t)}({\bm{Q}}_{k,l}^{(t+1)}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)})$
$\displaystyle\Lambda_{3}^{(t+1)}$
$\displaystyle=\Lambda_{3}^{(t)}+\mu_{3}^{(t)}({\bm{Y}}^{(t+1)}{\bm{1}}-{\bm{1}})$
$\displaystyle\Lambda_{4}^{(t+1)}$
$\displaystyle=\Lambda_{4}^{(t)}+\mu_{4}^{(t)}({\bm{Y}}^{(t+1)}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)}))$

As $\Lambda$s preserve the integral of all residuals with respect to the
constrain along the update progress, ADMM feeds back the difference stepwisely
to drive a zero error on the constrain.

###  A.7 Update ${\bm{\mu}}^{(t+1)}$s

Finally, we consider adaptive penalty parameters $\mu^{(t+1)}$s to lift the
burden of parameter tuning. Such updating rule has been proven to have a fast
convergence speed, see [62]. For $i=1,2,3,4$, we define

$\displaystyle\begin{aligned}
\mu_{i}^{(t+1)}&=\min\left(\rho\mu_{i}^{(t)},\mu_{i,\max}\right).\end{aligned}$
(20)

Here $\mu_{i,\max}$, for $i=1,2,3,4$, is an upper bound over all the adpative
penalty parameters, and the value of $\rho\,(\geq 1)$ is fixed to an
appropriate constant.

###### Remark 6

The $\rho$ can also be defined as a piece-wise function where the split
condition comes from stopping criteria analysis, for instance, see [45].

##  Appendix B Proof of Theorem 3

Step 1: The proof of claim 1)

Consider the iteration scheme for ${\bm{Z}}^{(t+1)}$. On the $(t+1)$th
iteration, according to the updating rule of ${\bm{E}}^{(t+1)}$, its first-
order optimality condition holds, i.e.,

$\mathbf{0}\in\partial\|{\bm{Z}}^{(t+1)}\|_{1}-\mu^{(t)}_{4}({\bm{Y}}^{(t+1)}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)}))-\Lambda^{(t)}_{4}$

According to the second last rule in (14), we immediately see that

$\displaystyle\Lambda^{(t+1)}_{4}\in\partial\|{\bm{Z}}^{(t+1)}\|_{1}.$ (21)

According to the fact, when $y\in\partial\|x\|$ where $\|\cdot\|$ is a norm,
we have $\|y\|_{\text{dual}}\leq 1$. Specifically $\|y\|_{\text{dual}}=1$ if
$x\not=0$, otherwise $\|y\|_{\text{dual}}\leq 1$. As the dual norm of the
$\ell_{1}$ norm is $\ell_{\infty}$, hence $\Lambda^{(t+1)}$ is bounded.

Now consider each ${\bm{Q}}_{k,l}$. The first order condition of the update
rule (13) is given by

$\mathbf{0}\in\nu_{k,l}\partial\|{\bm{Q}}^{(t+1)}_{k,l}\|_{1,G}+\Lambda_{2;k,l}+\mu^{(t+1)}_{2}({\bm{Q}}_{k,l}-\mathcal{W}_{k,l}{\bm{U}}^{(t+1)}).$

This surely gives, as our argument above,

$\displaystyle\Lambda^{(t+1)}_{2;k,l}\in\nu_{k,l}\partial\|{\bm{Q}}^{(t+1)}_{k,l}\|_{1,G},$
(22)

hence all $\Lambda^{(t+1)}_{2;k,l}$ are bounded.

By the similar argument as the first step, we can claim that

$\displaystyle\Lambda^{(t+1)}_{1}\in\partial_{{\bm{E}}}\|{\bm{E}}^{(t+1)}\|_{2,1,G}.$
(23)

Similarly we can conclude that $\Lambda^{(t+1)}_{1}$ is bounded too.

Applying the similar strategy for $\ell_{2,1}$ norm, we can conclude that
$\Lambda^{(t+1)}_{1}$ is also bounded.

Regarding the boundedness of $\Lambda^{(t)}_{3}$, we note that the first order
condition for the update step 1 for ${\bm{Y}}^{(t+1)}$ gives

$\displaystyle-\mu^{(t)}_{1}({\bm{U}}^{(t+1)}-{\bm{Y}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)}){\bm{U}}^{(t+1)\top}$
$\displaystyle+\mu^{(t)}_{3}({\bm{Y}}^{(t+1)}\mathbf{1}-\mathbf{1})\mathbf{1}^{\top}+\mu^{(t)}_{4}({\bm{Y}}^{(t+1)}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)})$
$\displaystyle-\Lambda^{(t)}_{1}{\bm{U}}^{(t+1)\top}+\Lambda^{(t)}_{3}\mathbf{1}^{\top}+\Lambda^{(t)}_{4}=0,$

which is, based on the Lagrangian updating rules (14),

$\displaystyle-(\Lambda^{(t+1)}_{1}-\Lambda^{(t)}_{1}){\bm{U}}^{(t+1)\top}+(\Lambda^{(t+1)}_{3}-\Lambda^{(t)}_{3}){\bm{1}}^{\top}$
$\displaystyle+(\Lambda^{(t+1)}_{4}-\Lambda^{(t)}_{4})-\Lambda^{(t+1)}_{4}{\bm{U}}^{(t+1)}+\Lambda^{(t)}_{3}{\bm{1}}^{\top}+\Lambda^{(t)}_{4}=0.$

Hence

$\displaystyle\Lambda^{(t+1)}_{3}{\bm{1}}^{\top}=-\Lambda^{(t+1)}_{1}{\bm{U}}^{(t+1)}+\Lambda^{(t+1)}_{4}$
(24)

which shows $\Lambda^{(t+1)}_{3}$ is bounded.

Finally the monotonic property of ADMM gives

$\displaystyle{\mathcal{L}}({\bm{U}}^{(t+1)},{\bm{Z}}^{(t+1)},{\bm{E}}^{(t+1)},{\bm{Q}}^{(t+1)},{\bm{Y}}^{(t+1)};\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},$
$\displaystyle\;\;\;\;\;\;\Lambda^{(t)}_{3},\Lambda^{(t)}_{4},\mu^{(t)}_{1},\mu^{(t)}_{2},\mu^{(t)}_{3},\mu^{(t)}_{4})$
$\displaystyle\leq$
$\displaystyle{\mathcal{L}}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)};\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4},$
$\displaystyle\;\;\;\;\;\;\mu^{(t)}_{1},\mu^{(t)}_{2},\mu^{(t)}_{3},\mu^{(t)}_{4}).$

It is easy to re-write that

$\displaystyle{\mathcal{L}}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)};\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4},$
$\displaystyle\;\;\;\;\;\;\;\;\mu^{(t)}_{1},\mu^{(t)}_{2},\mu^{(t)}_{3},\mu^{(t)}_{4})$
$\displaystyle=$
$\displaystyle{\mathcal{L}}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)};\Lambda^{(t-1)}_{1},\Lambda^{(t-1)}_{2},\Lambda^{(t-1)}_{3},$
$\displaystyle\;\;\;\;\Lambda^{(t-1)}_{4},\mu^{(t-1)}_{1},\mu^{(t-1)}_{2},\mu^{(t-1)}_{3},\mu^{(t-1)}_{4})$
$\displaystyle+\frac{\mu^{(t)}_{1}-\mu^{(t-1)}_{1}}{2}\|{\bm{U}}^{(t)}-{\bm{Y}}^{(t)}{\bm{U}}^{(t)}-{\bm{E}}^{(t)}\|_{2}^{2}$
$\displaystyle+\sum_{k,l}\frac{\mu^{(t)}_{2}-\mu^{(t-1)}_{2}}{2}\|{\bm{Q}}^{(t)}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t)}\|_{2}^{2}$
$\displaystyle+\frac{\mu^{(t)}_{3}-\mu^{(t-1)}_{3}}{2}\|{\bm{Y}}^{(t)}{\bm{1}}-{\bm{1}}\|_{2}^{2}$
$\displaystyle+\frac{\mu^{(t)}_{4}-\mu^{(t-1)}_{4}}{2}\|{\bm{Y}}^{(t)}-{\bm{Z}}^{(t)}+\text{diag}({\bm{Z}}^{(t)})\|^{2}_{2}$
$\displaystyle+\text{tr}\left((\Lambda_{1}^{(t)\top}-\Lambda_{1}^{(t-1)\top})({\bm{U}}^{(t)}-{\bm{Y}}^{(t)}{\bm{U}}^{(t)}-{\bm{E}}^{(t)})\right)$
$\displaystyle+\sum_{k,l}\text{tr}\left((\Lambda_{2;k,l}^{(t)\top}-\Lambda_{2;k,l}^{(t-1)\top})({\bm{Q}}_{k,l}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t)})\right)$
$\displaystyle+(\Lambda_{3}^{(t)\top}-\Lambda_{3}^{(t-1)\top})({\bm{Y}}^{(t)}{\bm{1}}-{\bm{1}})$
$\displaystyle+\text{tr}((\Lambda^{(t)\top}_{4}-\Lambda^{(t-1)\top}_{4})({\bm{Y}}^{(t)}-{\bm{Z}}^{(t)}+\text{diag}({\bm{Z}}^{(t)})))$
$\displaystyle=$ RHS

With the updating rules in Step 6 of the algorithm, we have

RHS
$\displaystyle={\mathcal{L}}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)};\Lambda^{(t-1)}_{1},\Lambda^{(t-1)}_{2},$
$\displaystyle\;\;\;\;\;\Lambda^{(t-1)}_{3},\Lambda^{(t-1)}_{4},\mu^{(t-1)}_{1},\mu^{(t-1)}_{2},\mu^{(t-1)}_{3},\mu^{(t-1)}_{4})$
$\displaystyle+\frac{\mu^{(t)}_{1}+\mu^{(t-1)}_{1}}{2\mu^{2(t-1)}_{1}}\|\Lambda_{1}^{(t)}-\Lambda_{1}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\sum_{k,l}\frac{\mu^{(t)}_{2}+\mu^{(t-1)}_{2}}{2\mu^{2(t-1)}_{2}}\|\Lambda_{2;k,l}^{(t)}-\Lambda_{2;k,l}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\frac{\mu^{(t)}_{3}+\mu^{(t-1)}_{3}}{2\mu^{2(t-1)}_{3}}\|\Lambda_{3}^{(t)}-\Lambda_{3}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\frac{\mu^{(t)}_{4}+\mu^{(t-1)}_{4}}{2\mu^{2(t-1)}_{4}}\|\Lambda^{(t)}_{4}-\Lambda^{(t-1)}_{4}\|^{2}_{F}.$

Summing the above result over $t=0,1,...,T(\geq 1)$ gives

$\displaystyle{\mathcal{L}}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)};\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4},$
$\displaystyle\;\;\;\;\;\;\;\mu^{(t)}_{1},\mu^{(t)}_{2},\mu^{(t)}_{3},\mu^{(t)}_{4})$
$\displaystyle\leq$
$\displaystyle{\mathcal{L}}({\bm{U}}^{(1)},{\bm{Z}}^{(1)},{\bm{E}}^{(1)},{\bm{Q}}^{(1)},{\bm{Y}}^{(1)};\Lambda^{(0)}_{1},\Lambda^{(0)}_{2},\Lambda^{(0)}_{3},\Lambda^{(0)}_{4},$
$\displaystyle\;\;\;\;\mu^{(0)}_{1},\mu^{(0)}_{2},\mu^{(0)}_{3},\mu^{(0)}_{4})$
$\displaystyle+\sum^{T}_{t=1}\frac{\mu^{(t)}_{1}+\mu^{(t-1)}_{1}}{2\mu^{2(t-1)}_{1}}\|\Lambda_{1}^{(t)}-\Lambda_{1}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\sum^{T}_{t=1}\sum_{k,l}\frac{\mu^{(t)}_{2}+\mu^{(t-1)}_{2}}{2\mu^{2(t-1)}_{2}}\|\Lambda_{2;k,l}^{(t)}-\Lambda_{2;k,l}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\sum^{T}_{t=1}\frac{\mu^{(t)}_{3}+\mu^{2(t-1)}_{3}}{2\mu^{2(t-1)}_{3}}\|\Lambda_{3}^{(t)}-\Lambda_{3}^{(t-1)}\|^{2}_{F}$
$\displaystyle+\sum^{T}_{t=1}\frac{\mu^{(t)}_{4}+\mu^{2(t-1)}_{4}}{2\mu^{2(t-1)}_{4}}\|\Lambda^{(t)}_{4}-\Lambda^{(t-1)}_{4}\|^{2}_{F}.$

Note that, for $i=1,2,3,4$

$\displaystyle\sum^{T}_{t=1}\frac{\mu^{(t)}_{i}+\mu^{(t-1)}_{i}}{2\mu^{2(t-1)}_{i}}=\sum^{T}_{t=1}\frac{\rho\mu^{(t-1)}_{i}+\mu^{(t-1)}_{i}}{2\mu^{2(t-1)}_{i}}$
$\displaystyle=$
$\displaystyle\frac{\rho+1}{2}\sum^{T}_{t=1}\frac{1}{\mu^{(t-1)}_{i}}=\frac{\rho+1}{2\mu^{(0)}_{i}}\sum^{T}_{t=1}\frac{1}{\rho^{t-1}}\leq\frac{\rho(\rho+1)}{2\mu^{(0)}_{i}(\rho-1)}.$

Note that
${\mathcal{L}}({\bm{U}}^{(1)},{\bm{Z}}^{(1)},{\bm{E}}^{(1)},{\bm{Q}}^{(1)},{\bm{Y}}^{(1)};\Lambda^{(0)}_{1},\Lambda^{(0)}_{2},\Lambda^{(0)}_{3},$
$\Lambda^{(0)}_{4},\mu^{(0)}_{1},\mu^{(0)}_{2},\mu^{(0)}_{3},\mu^{(0)}_{4})$
is finite, and squences
$\\{\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4}\\}$
and $\sum^{T}_{t=1}\frac{\mu^{(t)}_{i}+\mu^{(t-1)}_{i}}{2\mu^{2(t-1)}_{i}}$,
so the objective is bounded. Now Therefore, we can conclude that all the
relevant sequences $\\{\Gamma_{t}\\}$ are unbounded given the finite objective
function values. This completes the proof of claim 1) of Theorem 1.

Step 2: The proof of claim 2)

According to Bolzano-Weierstrass theorem, it follows from the boundedness of
the sequence $\\{\Gamma_{t}\\}_{t=1}^{\infty}$ that
$\\{\Gamma_{t}\\}_{t=1}^{\infty}$ must have a convergent subsequence. Without
loss of generality, suppose that the subsequence of
$\\{\Gamma_{t}\\}_{t=1}^{\infty}$ is represented by itself, and it converges
to an accumulation point, denoted as $\Gamma_{*}$. That is, we have
established

$\displaystyle\lim\limits_{t\rightarrow\infty}({\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Q}}^{(t)},{\bm{Y}}^{(t)},\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4})$
$\displaystyle=$
$\displaystyle({\bm{U}}_{*},{\bm{Z}}_{*},{\bm{E}}_{*},{\bm{Q}}_{*},{\bm{Y}}_{*},\Lambda_{1*},\Lambda_{2*},\Lambda_{3*},\Lambda_{4*}).$
(25)

The updating rules of Lagrange multipliers (17) imply that

$\displaystyle{\bm{U}}^{(t+1)}-{\bm{Z}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t+1)}=(\Lambda_{1}^{(t+1)}-\Lambda_{1}^{(t)})/\mu_{1}^{(t)};$
$\displaystyle{\bm{Q}}_{k,l}^{(t+1)}-{\mathcal{W}}_{k,l}{\bm{U}}^{(t+1)}=(\Lambda_{2;k,l}^{(t+1)}-\Lambda_{2;k,l}^{(t)})/\mu_{2}^{(t)};$
$\displaystyle{\bm{Y}}^{(t+1)}{\bm{1}}-{\bm{1}}=(\Lambda_{3}^{(t+1)}-\Lambda_{3}^{(t)})/\mu_{3}^{(t)};$
$\displaystyle{\bm{Y}}^{(t+1)}-{\bm{Z}}^{(t+1)}+\text{diag}({\bm{Z}}^{(t+1)})=(\Lambda_{4}^{(t+1)}-\Lambda_{4}^{(t)})/\mu_{4}^{(t)};$

By the boundedness of the sequences
$\\{\Lambda^{(t)}_{1},\Lambda^{(t)}_{2},\Lambda^{(t)}_{3},\Lambda^{(t)}_{4}\\}$,
and the fact $\lim\limits_{t\rightarrow\infty}\mu^{(t)}_{i}=\infty$, we
further have, by taking limit,

$\displaystyle{\bm{U}}_{*}$
$\displaystyle={\bm{Y}}_{*}{\bm{U}}_{*}+{\bm{E}}_{*};\;\mathcal{Q}_{*}=\mathcal{W}_{*}{\bm{U}}_{*}$
$\displaystyle{\bm{Y}}_{*}\mathbf{1}$
$\displaystyle=\mathbf{1};\;{\bm{Y}}_{*}={\bm{Z}}_{*}-\text{diag}({\bm{Z}}_{*})$

The other four KKT conditions can be obtained by taking limits over (21) -
(24). The proof of claim 2) is completed.

Step 3: The proof of claim 3)

We conduct the proof in the following sub-steps.

(1) $\\{{\bm{Y}}^{(t)}\\}$ is a Cauchy:

According to the updating rules (14), we have

$\displaystyle\Lambda_{3}^{(t+1)}=\Lambda_{3}^{(t)}+\mu_{3}^{(t)}({\bm{Y}}^{(t+1)}{\bm{1}}-{\bm{1}})$
$\displaystyle\Lambda_{3}^{(t)}=\Lambda_{3}^{(t-1)}+\mu_{3}^{(t-1)}({\bm{Y}}^{(t)}{\bm{1}}-{\bm{1}})$

This gives

$\displaystyle{\bm{Y}}^{(t+1)}{\bm{1}}-{\bm{Y}}^{(t)}{\bm{1}}=\frac{\Lambda^{(t+1)}_{3}-\Lambda^{(t)}_{3}}{\mu^{(t)}_{3}}+\frac{\Lambda^{(t)}_{3}-\Lambda^{(t-1)}_{3}}{\mu^{(t-1)}_{3}}$

Hence

$\displaystyle\|{\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}\|^{2}_{F}\leq\|({\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}){\bm{1}}\|^{2}_{F}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\mu^{2(t)}_{3}}\|\Lambda^{(t+1)}_{3}-\Lambda^{(t)}_{3}\|^{2}_{F}+\frac{1}{\mu^{2(t-1)}_{3}}\|\Lambda^{(t)}_{3}-\Lambda^{(t-1)}_{3}\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\frac{1}{\mu^{2(t)}_{3}}(\|\Lambda^{(t+1)}_{3}-\Lambda^{(t)}_{3}\|^{2}_{F}+\rho^{2}\|\Lambda^{(t)}_{3}-\Lambda^{(t-1)}_{3}\|^{2}_{F})$
$\displaystyle=:$ $\displaystyle\frac{1}{\mu^{2(t)}_{3}}C_{t}$

where $C_{t}$ is bounded as the sequence $\\{\Gamma_{t}\\}$ is bounded. Then
for any $m>n$, we can establish

$\displaystyle\|{\bm{Y}}^{(m)}-{\bm{Y}}^{(n)}\|^{2}_{F}$
$\displaystyle\leq\sum^{m-1}_{t=n}\|{\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}\|^{2}_{F}$
$\displaystyle\leq
C\sum^{m-1}_{t=n}\frac{1}{\mu^{2(t)}_{1}}=C\sum^{m-1}_{t=n}\frac{1}{\mu^{(0)}_{i}\rho^{2t}}$

where $\rho>1$. This completes the proof of Cauchy property of the sequence
$\\{{\bm{Y}}^{(t)}\\}$.

(2) $\\{{\bm{U}}^{(t)}\\}$ is a Cauchy:

The first order condition for the optimal ${\bm{U}}^{(t+1)}$ reads as,
denoting by ${\bm{L}}^{(t+1)}={\bm{I}}-{\bm{Y}}^{(t+1)}$,

$\displaystyle 0=$
$\displaystyle\lambda_{2}{\bm{D}}({\bm{U}}^{(t+1)}-{\bm{X}})+\mu^{(t)}_{1}{\bm{L}}^{(t+1)\top}({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}-{\bm{E}}^{(t)})$
$\displaystyle+\sum_{k,l}\mu^{(t)}_{2}\mathcal{W}_{k,l}^{\top}({\bm{Q}}^{(t)}_{k,l}-\mathcal{W}_{k,l}{\bm{U}}^{(t+1)})$
$\displaystyle+{\bm{L}}^{(t+1)\top}\Lambda^{(t)}_{1}+\sum_{k,l}\mathcal{W}_{k,l}\Lambda^{(t)}_{2;k,l}$

Similarly by using the updating rules (14) and replacing ${\bm{E}}^{(t)}$ and
${\bm{Q}}^{(t)}_{k,l}$, we will have

$\displaystyle 0=$
$\displaystyle\lambda_{2}{\bm{D}}({\bm{U}}^{(t+1)}-{\bm{X}})+\mu^{(t)}_{1}{\bm{L}}^{(t+1)}{\bm{L}}^{(t+1)\top}({\bm{U}}^{(t+1)}-{\bm{U}}^{(t)})$
$\displaystyle+\sum_{k,l}\mu^{(t)}_{2}\mathcal{W}^{T}_{k,l}\mathcal{W}_{k,l}({\bm{U}}^{(t+1)}-{\bm{U}}^{(t)})$
$\displaystyle-$
$\displaystyle\rho{\bm{L}}^{(t+1)}(\Lambda^{(t)}_{1}-\Lambda^{(t-1)}_{1})-\rho\sum_{k,l}\mathcal{W}^{\top}_{k,l}(\Lambda^{(t)}_{2}-\Lambda^{(t-1)}_{2})$
$\displaystyle+\mu^{(t)}_{1}{\bm{L}}^{(t+1)\top}({\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}){\bm{U}}^{(t)}+{\bm{L}}^{(t+1)T}\Lambda^{(t)}_{1}$
$\displaystyle+\sum_{k,l}\mathcal{W}^{\top}_{k,l}\Lambda^{(t)}_{2;k,l}$

Assuming $\mu^{(t)}_{1}=\mu^{(t)}_{2}=\mu^{(t)}$ will give

$\displaystyle\|{\bm{U}}^{(t+1)}-{\bm{U}}^{(t)}\|^{2}_{F}$ $\displaystyle\leq$
$\displaystyle\frac{1}{\mu^{2(t)}}\|({\bm{L}}^{(t+1)}{\bm{L}}^{(t+1)\top}+\sum_{k,l}\mathcal{W}^{\top}_{k,l}\mathcal{W}_{k,l})^{-1}\|^{2}_{F}\cdot$
$\displaystyle\biggl{\\{}\lambda_{2}\|{\bm{D}}({\bm{U}}^{(t+1)}-{\bm{X}})\|^{2}_{F}+\rho^{2}\|{\bm{L}}^{(t+1)}(\Lambda^{(t)}_{1}-\Lambda^{(t-1)}_{1})\|^{2}_{F}\biggr{.}$
$\displaystyle+\rho^{2}\sum_{k,l}\|\mathcal{W}^{\top}_{k,l}(\Lambda^{(t)}_{2}-\Lambda^{(t-1)}_{2})\|^{2}_{F}+\|{\bm{L}}^{(t+1)\top}\Lambda^{(t)}_{1}\|^{2}_{F}$
$\displaystyle+\mu^{(t)}_{1}\|{\bm{L}}^{(t+1)\top}\|^{2}_{F}\|({\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}){\bm{U}}^{(t)}\|^{2}_{F}$
$\displaystyle\biggl{.}+\sum_{k,l}\|\mathcal{W}^{\top}_{k,l}\Lambda^{(t)}_{2;k,l}\|^{2}_{F}\biggr{\\}}.$

We consider the second last term where we have seen the term $\mu^{(t)}_{1}$.
First we note that
$\|{\bm{L}}^{(t+1)\top}\|^{2}_{F}\leq\|({\bm{I}}-{\bm{Y}}^{(t+1)}){\bm{1}}\|_{F}$
while
$({\bm{I}}-{\bm{Y}}^{(t+1)}){\bm{1}}=(\Lambda^{(t)}_{3}-\Lambda^{(t+1)}_{3})/\mu^{(t)}_{3}$.
Hence

$\mu^{(t)}_{1}\|{\bm{L}}^{(t+1)\top}\|^{2}_{F}\leq\frac{\mu^{(t)}_{1}}{\mu^{(t)}_{3}}\|\Lambda^{(t)}_{3}-\Lambda^{(t+1)}_{3}\|^{2}_{F}=\|\Lambda^{(t)}_{3}-\Lambda^{(t+1)}_{3}\|^{2}_{F}$

if we take $\mu^{(t)}_{1}=\mu^{(t)}_{3}$. This means there exists a constant
$C$ such that

$\|{\bm{U}}^{(t+1)}-{\bm{U}}^{(t)}\|^{2}_{F}\leq\frac{1}{\mu^{2(t)}}C.$

Similar argument as for Cauchy property for $\\{{\bm{Y}}^{(t)}\\}$, we can
claim $\\{{\bm{U}}^{(t)}\\}$ is Cauchy.

(3) $\\{{\bm{E}}^{(t)}\\},\\{{\bm{Z}}^{(t)}\\},\\{{\bm{Q}}^{(t)}_{k,l}\\}$ all
are Cauchy:

As these variables share similar pattern with L1-norm constraint, to save the
space, we only take $\\{{\bm{E}}^{(t)}\\}$ as an example to show its Cauchy
property. We will follow the updating rules (14) again.

$\displaystyle\|{\bm{E}}^{(t+1)}-{\bm{E}}^{(t)}\|^{2}_{F}$ $\displaystyle=$
$\displaystyle\|{\bm{E}}^{(t+1)}-({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\Lambda^{(t)}_{1}/{\mu^{(t)}_{1}})$
$\displaystyle+({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\Lambda^{(t)}_{1}/{\mu^{(t)}_{1}})-{\bm{E}}^{(t)}\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\|{\bm{E}}^{(t+1)}-({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\Lambda^{(t)}_{1}/{\mu^{(t)}_{1}})$
$\displaystyle+({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\Lambda^{(t)}_{1}/{\mu^{(t)}_{1}})$
$\displaystyle+(\Lambda^{(t)}_{1}-\Lambda^{(t-1)}_{1})/{\mu^{(t-1)}}-{\bm{L}}^{(t)}{\bm{U}}^{(t)}\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\|{\bm{E}}^{(t+1)}-({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{1}{\mu^{(t)}_{1}}\Lambda^{(t)}_{1})-{\bm{U}}^{(t)}+{\bm{L}}^{(t)}{\bm{U}}^{(t)}$
$\displaystyle+$
$\displaystyle{\bm{U}}^{(t+1)}-{\bm{Y}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{1}{\mu^{(t)}_{1}}(\Lambda^{(t)}_{1}-\rho(\Lambda^{(t-1)}_{1}-\Lambda^{(t)}_{1}))\|^{2}_{F}$
$\displaystyle\leq$
$\displaystyle\|{\bm{E}}^{(t+1)}-({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{1}{\mu^{(t)}_{1}}\Lambda^{(t)}_{1})\|^{2}_{F}$
$\displaystyle+\|{\bm{Y}}^{(t+1)}\|^{2}_{F}\|{\bm{U}}^{(t+1)}-{\bm{U}}^{(t)}\|^{2}_{F}$
$\displaystyle+\|{\bm{U}}^{(t)}\|^{2}_{F}\|{\bm{Y}}^{(t+1)}-{\bm{Y}}^{(t)}\|^{2}_{F}$
$\displaystyle+\frac{1}{\mu^{2(t)}}\|\Lambda^{(t)}_{1}-\rho(\Lambda^{(t-1)}_{1}-\Lambda^{(t)}_{1})\|^{2}_{F}.$

Now we prove the first term is second order controlled. ${\bm{E}}^{(t+1)}$ is
obtained by $\ell_{1}$-threshold operator, so each row element of the first
term is less than $1/\mu^{(t)}_{1}$, hence

$\|{\bm{E}}^{(t+1)}-({\bm{L}}^{(t+1)}{\bm{U}}^{(t+1)}+\frac{1}{\mu^{(t)}_{1}}\Lambda^{(t)}_{1})\|^{2}_{F}\leq\frac{N}{\mu^{2(t)}_{1}},$

where $N$ is the number of rows of ${\bm{E}}^{(t+1)}$. We have already proved
that $\\{{\bm{Y}}^{(t)}\\}$ and $\\{{\bm{U}}^{(t)}\\}$. With the boundedness
of those variable we can claim that $\\{{\bm{E}}^{(t)}\\}$ is Cauchy.

Final conclusion: Given the Cauchy property, we can claim that the sequence
$\\{{\bm{U}}^{(t)},{\bm{Z}}^{(t)},{\bm{E}}^{(t)},{\bm{Y}}^{(t)},{\bm{Q}}^{(t)}\\}$
are convergent to its critical point.

This completes the proof of Theorem 3.