remarkRemark hypothesisHypothesis claimClaim Inhomogeneous graph trend
filtering via a $\ell_{2,0}$ cardinality penalty X Huang, A. Ang, J Zhang, Y
Wang

# Inhomogeneous graph trend filtering via a $\ell_{2,0}$ cardinality
penalty††thanks:  Submitted to the editors DATE. X. Huang and A. Ang are equal
contributors to this work and designated as co-first authors. This work was
funded by the Fog Research Institute under contract no. FRI-454.

Xiaoqing Huang  Department of Biostatistics and Health Data Science, Indiana
University School of Medicine, Indianapolis, IN 46202, USA.  Andersen Ang
School of Electronics and Computer Science, University of Southampton,
Southampton, SO17 1BJ, UK.  Jie Zhang Department of Medical and Molecular
Genetics, Indiana University, Indianapolis, IN 46202, USA.  Yijie Wang
Department of Computer Science, Indiana University, Bloomington, IN 47408,
USA<EMAIL_ADDRESS>
###### Abstract

We study estimation of piecewise smooth signals over a graph. We propose a
$\ell_{2,0}$-norm penalized Graph Trend Filtering (GTF) model to estimate
piecewise smooth graph signals that exhibits inhomogeneous levels of
smoothness across the nodes. We prove that the proposed GTF model is
simultaneously a k-means clustering on the signal over the nodes and a minimum
graph cut on the edges of the graph, where the clustering and the cut share
the same assignment matrix. We propose two methods to solve the proposed GTF
model: a spectral decomposition method and a method based on simulated
annealing. In the experiment on synthetic and real-world datasets, we show
that the proposed GTF model has a better performances compared with existing
approaches on the tasks of denoising, support recovery and semi-supervised
classification. We also show that the proposed GTF model can be solved more
efficiently than existing models for the dataset with a large edge set.

###### keywords:

graph trend filtering, $\ell_{2,0}$-norm, cardinality, k-means, minimum graph
cut, clustering

65F50, 68U01, 68R01

##  1 Introduction

Estimating signal from noisy observations is a classical problem in signal
processing. This work focuses on Graph Signal Estimation (GSE), which is a
powerful tool for analyzing signals that are defined on irregular and complex
domains (e.g., social networks, brain networks, and sensor networks) and has a
lot of real-world applications [10, 9, 8, 21, 13]. Fig. 1a gives an
illustration of a graph signal. In a graph signal, the signal is residing on
each node in a graph. The signal could be a scalar-valued signal (as shown in
Fig. 1a) or vector-valued signal (as shown Fig. 1c). The edges between nodes
in the graph introduce additional information and assumption for the signals
residing on the nodes. A typical assumption is that the graph signal is
globally smooth with respect to the intrinsic structure of the underlying
graph, i.e., the signals on the neighbor nodes tend to be similar, or in other
words, the local variation of the signals around each node is low. Such global
smoothness of a graph signal is further mathematically defined and
quantitatively measured by the $p$-Dirichlet form [24, 19, 3, 5, 13, 25] and
the graph Laplacian quadratic form ($p$-Dirichlet form with $p=2$) [25]. GSE
based on such global smoothness assumption has been investigated in graph
signal processing [24, 19, 10, 9, 8, 21, 13] and in the context of Laplacian
regularization [3, 32].

Although many existing works [24, 19, 10, 9, 8, 21, 13] assumes the global
smoothness on the graph signal, many real-world graph signals exhibit
inhomogeneous levels of smoothness over the graph [29, 31, 20]. For example,
in gene co-expression networks, gene expressions within a gene community tend
to be homogeneous (i.e., having similar values), but the gene expression
across the two communities can be highly different (i.e., inhomogeneous). In
other words, there is a discontinuity in the first-order behaviour of the
signal value across the graph, see Fig. 1(a) for a toy illustration.

Such an “intra-community homogeneous, inter-community inhomogeneous” property
indicates that the graph signals have large variations between different
regions and have a small variation within regions (as shown in Fig. 1a). In
this paper, we study GSE for piecewise smooth graph signals, and we propose a
model: a $\ell_{2,0}$-norm penalized 1st-order Graph Trend Filtering (GTF)
[31], for extracting patterns of the signal on a graph with localized
discontinuities. Before we state explicitly our contribution, below we first
review about the GTF model.

Figure 1: Toy examples. (a) The ground truth of an inhomogeneous scalar-valued
graph signal. (b) The same graph signal in (a) with added Gaussian noise
$\mathcal{N}(0,0.1)$. (c) Comparison of the recovery of the ground truth in
(a) from 20 different noisy inhomogeneous graph signals as in (b) with
different parameter $\lambda$ in (1) and (2). The proposed model (2) achieves
the lowest mean square error when screening the parameter $\lambda$. (d) An
example of a vector-valued graph signal.

###  1.1 GTF for scalar-valued graph signals

Given a graph $G(V,E)$ of $n$ nodes $V=\\{1,\dots,n\\}$ and edges
$E=\\{(i,j)\,|\,i,j\in V\text{ that are connected}\\}$, let
$\bm{y}\in\mathbb{R}^{n}$ be the observed signal measurement over the nodes,
i.e., the element $y_{i}\in\mathbb{R}$ is the scalar value at node $i$. The
goal of the 1st-order GTF for scalar-valued graph signals is to find the
estimated graph signal $\bm{\beta}\in\mathbb{R}^{n}$, where the $i$th element
in $\bm{\beta}$ denoted by $\beta_{i}$ (the estimated signal at node $i$), by
solving

(1)
$\operatorname*{\textrm{argmin}}_{\bm{\beta}\in\mathbb{R}^{n}}\underbrace{\dfrac{1}{2}\|\bm{y}-\bm{\beta}\|_{2}^{2}}_{\text{fitting
/ denoise}}+\lambda\underbrace{\sum_{(i,j)\in
E}\Omega(\beta_{i},\beta_{j})}_{\text{GTF regularization}}.$

In (1), $\lambda\geq 0$ is a regularization parameter and
$\Omega(\beta_{i},\beta_{j}):\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
is a penalty term that penalizes the difference between the signals on nodes
$\\{i,j\\}$. In other words, $\Omega(\beta_{i},\beta_{j})$ is a penalization
of the neighborhood discrepancy of the graph signal, and it encourages the
graph signal value at node $\\{i,j\\}$ to be similar. Different $\Omega$ have
been proposed in the literature for GTF, see Table 1.

Table 1: Different $\Omega$ proposed in the literature for GTF. Paper | $\Omega(\beta_{i},\beta_{j})$ | Name | Parameter
---|---|---|---
[31] | $|\beta_{i}-\beta_{j}|$ |  $\ell_{1}$ norm | -
[29] | $\displaystyle\int_{0}^{|\beta_{i}-\beta_{j}|}\min\left(1,\frac{\left(\gamma-u/\lambda\right)_{+}}{\gamma-1}\right)du$ | Smoothly Clipped Absolute Deviation | $\gamma\geq 2$
[29] | $\displaystyle\int_{0}^{|\beta_{i}-\beta_{j}|}\left(1-\frac{u}{\gamma\lambda}\right)_{+}du$ | Minimax Concave Penalty | $\gamma$

####  1.1.1 GTF with the counting norm penalty

In this paper, we propose the norm $\Omega$ as
$\Omega(\beta_{i},\beta_{j})\,=\,i_{\neq 0}(\beta_{i}-\beta_{j})$, and propose
the GTF model as

(2)
$\operatorname*{\textrm{argmin}}_{\bm{\beta}\in\mathbb{R}^{n}}\frac{1}{2}\|\bm{y}-\bm{\beta}\|_{2}^{2}+\lambda\sum_{(i,j)\in
E}i_{\neq 0}(\beta_{i}-\beta_{j}),$

where $i_{\neq 0}:\mathbb{R}\rightarrow\mathbb{R}$ is an indicator function
defined as $i_{\neq 0}(x)=0$ if $x=0$ and $i_{\neq 0}(x)=1$ otherwise. In
summary, we propose to use a counting norm for $\Omega$, which counts the
number of edges that connects nodes with different signal value.

####  1.1.2 The advantage of using the proposed norm penalty

The current 1st-order GTF model (1) uses $\ell_{1}$ norm [31], Smoothly
Clipped Absolute Deviation (SCAD) [29], or Minimax Concave Penalty (MCP) [29]
to penalize the difference between the signals over connected nodes (as shown
in Table1). However, the resulting estimates generated by these penalties are
“biased toward zero”: the $\ell_{1}$ penalty enforces the minimization of the
term $|\beta_{i}-\beta_{j}|$ and, therefore, impose $|\beta_{i}-\beta_{j}|$
shrink to zero. As a result, the $\ell_{1}$ penalty forces $\beta_{i}$ and
$\beta_{j}$ to be numerically close even if they are not. While the proposed
counting norm penalty applies the same penalty to $|\beta_{i}-\beta_{j}|$ no
matter how different $\beta_{i}$ and $\beta_{j}$ are. Therefore,
$|\beta_{i}-\beta_{j}|$ is not enforced to shrink to zero when $\beta_{i}$ and
$\beta_{j}$ are different. The same comment also applied to SCAD and MCP
penalties.

We further illustrate such advantage of the proposed GTF model (2) numerically
in Fig. 1c: after screening over the model hyper-parameter ($\lambda$ in (1)
and (2), respectively), the GTF model (2) achieves the lowest mean square
error to recover the ground-truth signal as shown in Fig. 1(a) from the noisy
signal as shown in Fig. 1(b).

###  1.2 GTF for vector-valued graph signals

The previous discussion focused on scalar-valued graph signals. Now we extend
the model (2) to vector-valued graph signals (See Fig. 1d for a pictorial
illustration): let $\bm{y}_{i}\in\mathbb{R}^{d}$ be the vector-valued signal
at node $i$, let
$\bm{Y}=[\bm{y}_{1},\dots,\bm{y}_{n}]^{\top}\in\mathbb{R}^{n\times d}$ be the
data matrix (the rows are the observed vector-valued signals on nodes), we
consider the optimization problem

(3)
$\mathcal{P}_{0}~{}:~{}\operatorname*{\textrm{argmin}}_{\bm{B}\in\mathbb{R}^{n\times
d}}\frac{1}{2}\|\bm{Y}-\bm{B}\|_{F}^{2}+\lambda\underbrace{\sum_{(i,j)\in
E}i_{\neq
0}\Big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\Big{)}}_{\eqqcolon\mathcal{L}(\bm{B})},$

where $\bm{B}=[\bm{b}_{1},...,\bm{b}_{n}]^{\top}\in\mathbb{R}^{n\times d}$ and
$\bm{b}_{i}\in\mathbb{R}^{d}$ is the GTF signal estimate at node $i$. Note
that (2) is equivalent to (3) when the dimensionality $d=1$.

##### On the term $\mathcal{L}(\bm{B})$

We call $\mathcal{L}(\bm{B})$ in (3) the “$\ell_{2,0}$ norm” in this paper. It
can be treated as a “group counting norm”, which is the limiting case of the
matrix $\ell_{p,q}$ norm with $p=2$ and $q\rightarrow 0$ on a matrix $\bm{M}$
whose columns are the pairwise difference $\bm{b}_{i}-\bm{b}_{j}$ of $\bm{B}$
defined by the edge set $E$, i.e.,

$\mathcal{L}(\bm{B})~{}=~{}\|\bm{M}\|_{2,0}~{}=~{}\lim_{q\rightarrow
0}\|\bm{M}\|_{2,q}~{}\coloneqq~{}\lim_{q\rightarrow
0}\Bigg{(}\sum_{j=1}^{|E|}\bigg{(}\sqrt{\sum_{i=1}^{d}M_{ij}^{2}}\bigg{)}^{q}\Bigg{)}^{\frac{1}{q}}.$

If $\mathcal{L}(\bm{B})=0$ in $\mathcal{P}_{0}$, it means that all the node
signal estimates $\bm{b}_{i}$ share the same value, or equivalently each
columns of $\bm{B}$ is a constant vector.

###  1.3 Contributions and paper organization

Grounded on the advantage of using $\mathcal{L}$ (see Section 1.1.2) that it
overcomes the problems of using $\ell_{1}$ norm [31], SCAD [29], and MCP
penalties [29] in GTF, in this paper, we propose the 1st-order GTF model
$\mathcal{P}_{0}$.

  * •

We provide theoretical characterizations of the proposed GTF model. In section
2 we prove that $\mathcal{P}_{0}$ is simultaneously performing a k-means
clustering on signals $\bm{B}$ and a minimum graph cut on $G$ (Theorem 2.1 and
Theorem 2.3). This means that the model can produce a better result by fusing
the information of the node signal with the information of the graph structure
to deduce the signal estimate as well as the clustering membership value.

  * •

In section 3 we propose two algorithms to solve $\mathcal{P}_{0}$. We first
show that the proposed GTF problem is combinatorial in nature, see Theorem 2.1
and Theorem 2.3. We then propose a method to solve the problem efficiently
based on the standard technique of spectral approximation. We finally propose
a Hot-bath-based Simulated Annealing algorithm to find the global optimal
solution of the proposed GTF problem in case exact solution is required.

  * •

In section 4, we apply Theorem 2.1 to a semi-supervised learning problem known
as the Modified Absorption Problem (MAP) [31, 29], to showcase the analysis of
the proposed GTF model is generalizable to other models.

  * •

In section 5 we provide numerical results to support the effectiveness of the
algorithms proposed in this work, and to showcase the proposed GTF model
outperforms existing GTF models. In particular, we show that the proposed GTF
model has a better performances compared with existing approaches on the tasks
of denoising, support recovery and semi-supervised classification. We also
show that the proposed GTF model can be solved more efficiently than existing
models for the dataset with a large edge set.

##### Notation

We use $\\{$small italic, capital italic, small bold, capital bold$\\}$ font
to denote $\\{$scalar, set, vector, matrix$\\}$, respectively. We focus on
simple (without multi-edge) connected unweighted undirected graph $G(V,E)$ of
$n$ nodes with signal
$\bm{Y}=[\bm{y}_{1},\dots,\bm{y}_{n}]^{\top}\in\mathbb{R}^{n\times d}$ on $G$.
We let $\bm{A}\in\\{0,1\\}^{n\times n}$ be an adjacency matrix to present the
connections between nodes in $G$, with $A_{ij}=1$ indicating nodes $i$ and $j$
are connected and $A_{ij}=0$ otherwise. We let $\bm{D}$ be a diagonal matrix
with the degree of each node on its diagonal. Let $\bm{L}$ be the Laplacian of
$G$, we have $\bm{L}=\bm{D}-\bm{A}$. See [11] for more background on graph
theory. We let $[n]$ denotes integers from $1$ to $n$. Let
$\bm{X}\in\\{0,1\\}^{n\times k}$ be the assignment matrix, which can be used
to encode the cluster membership of $\\{\bm{y}_{1},\dots,\bm{y}_{n}\\}$ and
the node partition of $G$.

##  2 Equivalent formulation of $\mathcal{P}_{0}$

In this section we provide theoretical characterizations of the proposed GTF
model. We show that the proposed GTF model is equivalent to a simultaneous
k-means and min-cut (Theorem 2.1 and Theorem 2.3), which means that the model
is capable of fusing the information of the node signal with the information
of the graph structure to deduce the signal estimate as well as the clustering
membership.

###  2.1 Theoretical characterizations of the proposed GTF model

Theorem 2.1 shows that the effect of $\mathcal{L}(\bm{B})>0$ in
$\mathcal{P}_{0}$ is to cut the graph $G$ into $k$ partitions and enforce the
node signal in these $k$ partitions to be homogeneous (i.e., they share the
same signal value).

######  Theorem 2.1 (The effect of $\mathcal{L}(\bm{B})>0$ in
$\mathcal{P}_{0}$).

Given a data matrix $\bm{B}\in\mathbb{R}^{n\times d}$ over the nodes $V$ of a
simple connected unweighted undigraph graph $G(V,E)$, if
$\mathcal{L}(\bm{B})>0$, then it indicates the existence of an edge subset

$S_{=1}\coloneqq\Big{\\{}(i,j)\in E:i_{\neq
0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}=1\Big{\\}}$

that cuts $G$ into $k$ disjoint partitions $C=\\{C_{1},\dots,C_{k}\\}$.
Furthermore, $\mathcal{L}(\bm{B})>0$ enforces the node signals in those $k$
partitions sharing the same value. Encoding such partition by an assignment
matrix $\bm{X}\in\\{0,1\\}^{n\times k}$, then

(4a) $\displaystyle\mathcal{L}(\bm{B})$
$\displaystyle=\frac{1}{2}\operatorname{\textrm{tr}}(\bm{X}^{\top}\bm{L}\bm{X}),$
(4b) $\displaystyle\bm{B}$ $\displaystyle=\bm{X}\bm{\mu},$

where $\bm{\mu}=[\bm{\mu}_{1},...,\bm{\mu}_{k}]^{\top}\in\mathbb{R}^{k\times
d}$ with $\mu_{i}$ is the node signal for all the nodes in the $i$th
partition, and $\bm{L}$ is the graph Laplacian of $G$.

######  Proof 2.2.

First we let $E=S_{=1}\cup S_{=0}$ and $S_{=1}\cap S_{=0}=\varnothing$, where

$S_{=0}\coloneqq\Big{\\{}(i,j)\in E:i_{\neq
0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}=0\Big{\\}}.$

The proof has three parts: (1) show $S_{=1}$ is a graph cut, (2) prove (4a)
and (3) show (4b).

#####  $S_{=1}$ is a graph cut

We show that if we remove all edges in $S_{=1}$ from $G$, then for any edge
$(i,j)\in S_{=1}$, the corresponding nodes $\\{i,j\\}$ are in different
disjointed clusters.

  * •

Assume after removing all edges in $S_{=1}$ from $G$ there exists an edge
$(h,l)\in S_{=1}$ such that $i_{\neq
0}\big{(}\|\bm{b}_{h}-\bm{b}_{l}\|_{2}\big{)}=1$ but the nodes $\\{h,l\\}$ are
in the same partition.

  * •

As $G$ is connected, so nodes $\\{h,l\\}$ being in the same partition implies
there exists a path $p=\\{(h,s),\dots,(t,l)\\}$ between $h$ and $l$.

  * •

All the edges along the path $p$ are in $S_{=0}$ because by assumption
$S_{=0}$ contains all edges after removing $S_{=1}$ from $G$. This implies
$i_{\neq 0}\big{(}\|\bm{b}_{h}-\bm{b}_{s}\|_{2}\big{)}=\dots=i_{\neq
0}\big{(}\|\bm{b}_{t}-\bm{b}_{l}\|_{2}\big{)}=0$ and thus $i_{\neq
0}\big{(}\|\bm{b}_{h}-\bm{b}_{l}\|_{2}\big{)}=0$. A contradiction to $i_{\neq
0}\big{(}\|\bm{b}_{h}-\bm{b}_{l}\|_{2}\big{)}=1$.

##### Prove (4a)

$S_{=1}$ is a cut and it cuts $G$ into $k\in\mathbb{Z}^{+}$ ($k\leq n$ is an
unknown number) partitions. Let $\bm{X}\in\\{0,1\\}^{n\times k}$ be the
assignment matrix (where the $i$th row indicate the cluster membership of the
$i$th row of $\bm{Y}$) that encodes the cut $S_{=1}$, and let $\bm{L}$ be the
Laplacian of $G$, we have the cut size

$\begin{array}[]{rcll}\frac{1}{2}\operatorname{\textrm{tr}}(\bm{X}^{\top}\bm{L}\bm{X})&=&|S_{=1}|&\text{see
\cite[cite]{[\@@bibref{}{Fan2010,books/daglib/0004338}{}{}]}}\vspace{1mm}\\\
&=&\displaystyle\sum_{(i,j)\in S_{=1}}i_{\neq
0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}&\text{by definition of
$S_{=1}$}\\\ &=&\displaystyle\sum_{(i,j)\in S_{=1}}i_{\neq
0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}+\underbrace{\sum_{(i,j)\in
S_{=0}}i_{\neq 0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}}_{=0\text{ by
definition}}&\text{by definition of $S_{=0}$}\\\
&=&\displaystyle\sum_{(i,j)\in E}i_{\neq
0}\big{(}\|\bm{b}_{i}-\bm{b}_{j}\|_{2}\big{)}&E=S_{=1}\cup S_{=0}.\end{array}$

##### Prove (4b)

After removing edges in the cut $S_{=1}$ from $G$, the graph $G$ will be
separated into $k$ disjoint partitions. Denote such partition as
$C=\\{C_{1},...,C_{k}\\}$.

  * •

For each partition,

    * –

the nodes within the same partition are connected (because $G$ is connected);
and

    * –

all the edges connecting these nodes are in $S_{=0}$ (by
$S_{=0}=E\,\backslash\,S_{=1}$).

  * •

By definition of $S_{=0}$, the signal of all the nodes in each partition share
the same value, we denote such value by $\bm{\mu}_{i}$.

  * •

As $\bm{\mu}_{i}$ represents the signal value for all the nodes in partition
$C_{i}$, using the assignment matrix $\bm{X}$ gives (4b).

Now we see that (4b) tells that we assign all the nodes in the same partition
with the same signal value stored in $\bm{\mu}$.

Using Theorem 2.1, we now show that the solution to $\mathcal{P}_{0}$ is
simultaneously a k-means clustering and minimum graph cut.

######  Theorem 2.3 ($\mathcal{P}_{0}$ = k-means and min-cut).

$\mathcal{P}_{0}$ is same as

(5)
$\mathcal{P}_{1}~{}:~{}\operatorname*{\textrm{argmin}}_{\begin{subarray}{c}k\in\mathbb{Z}^{+}\\\
(\bm{X},\bm{\mu})\in\Psi\end{subarray}}\frac{1}{2}\underbrace{\|\bm{Y}-\bm{X}\bm{\mu}\|_{F}^{2}}_{\textnormal{k-means}}+\frac{\lambda}{2}\underbrace{\operatorname{\textrm{tr}}(\bm{X}^{\top}\bm{LX})}_{\textnormal{graph
$k$-cut}},$

where
$\Psi\coloneqq\Big{\\{}\bm{X},\bm{\mu}~{}|~{}\bm{\mu}\in\mathbb{R}^{k\times
d},~{}\bm{X}\in\\{0,1\\}^{n\times k},~{}\bm{X}\bm{1}_{k}=\bm{1}_{n}\Big{\\}}$
and $\bm{1}_{k}\in\mathbb{R}^{k}$ denotes vector of 1s.

######  Proof 2.4.

If $\mathcal{L}(\bm{B})>0$, apply (4a) and (4b) in Theorem 2.1 to substitute
$\mathcal{L}(\bm{B})$ and $\bm{B}$ in $\mathcal{P}_{0}$. If
$\mathcal{L}(\bm{B})=0$, it is trivial $\mathcal{P}_{0}=\mathcal{P}_{1}$ when
$k=1$.

###  2.2 Simultaneous k-means and graph cut

The cost function in $\mathcal{P}_{1}$ has two terms:

  1. 1.

The first term is in the form of k-means clustering on $\bm{Y}$ and $\bm{X}$
is the cluster assignment. This clustering depends solely on the data $\bm{Y}$
and ignores the graph $G$.

  2. 2.

The second term is in the form of minimum graph cut on $G$ and $\bm{X}$ is the
graph cut assignment matrix. This cut depends solely on the graph $G$ and
ignores the data $\bm{Y}$.

  3. 3.

Both k-means and graph cut share the same $\bm{X}$, indicating
$\mathcal{P}_{1}$ aims to find a $k$-partition that trades off between k-means
clustering (using only the signal on $\bm{Y}$) and minimum graph cut (using
only the topology on $G$). This means the information from the signal data and
the information of the graph structure are regularizing each other when we
solve for $\bm{B}$, in which this forms the basis why model $\mathcal{P}_{0}$
is worth studying:

    1. $-$

if the edge structure is consistent with the node data, we are reinforcing the
k-means on the node data using the information of the edge;

    2. $-$

if the edge structure is inconsistent with the signals on nodes, we are
correcting (regularizing) the k-means on the node data using the information
of the edge.

The above arguments provide explanations why the proposed model
$\mathcal{P}_{0}$ gave the best recovery in the toy example in Fig. 1. The
model $\mathcal{P}_{0}$ utilizes the two attributes of a graph signal: the
graph network structure (collective information) and the data on each node
(individual information), and when the model fuses the two type of information
to determine the value of $\bm{B}$.

We now see that model $\mathcal{P}_{0}$ is basically a mixed-integer
programming problem, which is NP-hard to solve:

  * •

We showed that $\mathcal{P}_{1}$ is equivalent to $\mathcal{P}_{0}$.

  * •

By Theorem 2.3, we know that solving $\mathcal{P}_{1}$ is an NP-hard problem
because the k-means clustering in $\mathcal{P}_{1}$ is NP-hard to solve [17].

##  3 The algorithms to solve $\mathcal{P}_{1}$

In the last section we see that problem $\mathcal{P}_{0}$ is basically a
mixed-integer program $\mathcal{P}_{1}$, where the integer $k\geq 1$ (the
number of clusters / partitions) is a unknown in the formulation. This number
can also be treated as the factorization rank in the k-means factorization
$\bm{B}=\bm{X}\bm{\mu}$. It is generally hard to find the optimal $k$ exactly
and thus in principle it is challenging to find optimal exact solution for
$\mathcal{P}_{1}$ (hence $\mathcal{P}_{0}$). In this section, we develop two
algorithms to solve $\mathcal{P}_{1}$ under two end-user scenarios:

  * •

For an approximate solution: a fast approach using spectral approximation.

  * •

For accurate solution: a slower approach based on hot-bath simulated
annealing.

###  3.1 Spectral approximation

To solve problem $\mathcal{P}_{1}$ efficiently, we propose a spectral method
based on the work of Mark Newman [18]. Here is the idea of the spectral
approximation: first we simplify $\mathcal{P}_{1}$ to $\mathcal{P}_{2}$ by
Lemma 3.1 to eliminate $\bm{\mu}$. Then we find an approximate solution
$(k,\bm{X})$ for $\mathcal{P}_{2}$ using brute-force on $k$ and spectral
method as follows:

  * •

We first consider $\mathcal{P}_{2}$ with a fixed $k$, denoted as
$\mathcal{P}_{2}^{k}$.

  * •

We solve $\mathcal{P}_{2}^{k}$ to get $(\bm{X}_{k},\bm{\mu}_{k})$ by the
approximate spectral method (Theorem 3.3 and Algorithm 1).

  * •

We screen through different values of $k$ to find the optimal approximate
solution $(k_{*},\bm{X}_{*})$ for $\mathcal{P}_{2}$ (Algorithm 2).

Note that brute-force search on $k$ is not scalable, but we argue that in
practice $k$ is usually small, and we can always narrow down the search scope
of $k$ using prior knowledge from the application.

##### On the diagonal of the matrix $\bm{X}^{\top}\bm{X}$

Before we proceed, we recall a useful fact. For $\bm{X}$ being a 0-1
assignment matrix, the matrix $\bm{X}^{\top}\bm{X}$ is a diagonal matrix that
its $i$th diagonal element is the size of the $i$th cluster, denoted as
$|C_{i}|$. That is

(6) $\bm{X}^{\top}\bm{X}~{}=~{}\begin{bmatrix}|C_{1}|\\\ &|C_{2}|\\\
&&\ddots\\\
&&&|C_{k}|\end{bmatrix},~{}~{}~{}~{}\operatorname{\textrm{tr}}\bm{X}^{\top}\bm{X}~{}=~{}\sum_{i=1}^{k}|C_{i}|=n.$

####  3.1.1 Eliminate $\bm{\mu}$

By observation, the variable $\bm{\mu}$ in $\mathcal{P}_{1}$ has no
constraint. We simplify $\mathcal{P}_{1}$ by eliminating $\bm{\mu}$ using
$\bm{\mu}_{*}=\bm{X}^{\dagger}\bm{Y}$, where
$\bm{X}^{\dagger}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}$ is the left pseudo-
inverse of $\bm{X}$. This gives the following lemma:

######  Lemma 3.1 (Simplifying $\mathcal{P}_{1}$).

$\mathcal{P}_{1}$ is equivalent to

(7)
$\mathcal{P}_{2}~{}:~{}\operatorname*{\textrm{argmax}}_{\begin{subarray}{c}k\in\mathbb{Z}^{+}\\\
\bm{X}\in\Phi\end{subarray}}~{}q~{}\coloneqq~{}\operatorname{\textrm{tr}}\Big{(}\bm{YY}^{\top}\bm{X}\bm{X}^{\dagger}\Big{)}+\lambda\operatorname{\textrm{tr}}\Big{(}\bm{X}^{\top}(-\bm{L})\bm{X}\Big{)},$

where $\Phi=\Big{\\{}\,\bm{X}~{}|~{}\bm{X}\in\\{0,1\\}^{n\times
k},\bm{X}\bm{1}_{k}=\bm{1}_{n}\,\Big{\\}}$ and
$\bm{X}^{\dagger}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}$.

######  Proof 3.2.

The subproblem on $\bm{\mu}$ in $\mathcal{P}_{1}$ has a closed-form solution
$\bm{\mu}^{*}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{Y}$. Then,

$\begin{array}[]{rl}\|\bm{Y}-\bm{X}\bm{\mu}^{*}\|_{F}^{2}=\|\bm{Y}-\bm{X}\bm{X}^{\dagger}\bm{Y}\|_{F}^{2}&=\|\bm{Y}\|_{F}^{2}-2\operatorname{\textrm{tr}}\Big{(}\bm{Y}^{\top}\bm{X}\bm{X}^{\dagger}\bm{Y}\Big{)}+\operatorname{\textrm{tr}}\Big{(}\big{(}\bm{X}\bm{X}^{\dagger}\bm{Y}\big{)}^{\top}\bm{X}\bm{X}^{\dagger}\bm{Y}\Big{)}\\\
&=\|\bm{Y}\|_{F}^{2}-2\operatorname{\textrm{tr}}\Big{(}\bm{Y}\bm{Y}^{\top}\bm{X}\bm{X}^{\dagger}\Big{)}+\operatorname{\textrm{tr}}\Big{(}\bm{Y}\bm{Y}^{\top}(\bm{X}\bm{X}^{\dagger})^{\top}\bm{X}\bm{X}^{\dagger}\Big{)}\\\
&=\|\bm{Y}\|_{F}^{2}-\operatorname{\textrm{tr}}\Big{(}\bm{Y}\bm{Y}^{\top}\bm{X}\bm{X}^{\dagger}\Big{)}\end{array}$

where the last equality is by expanding
$(\bm{X}\bm{X}^{\dagger})^{\top}\bm{X}\bm{X}^{\dagger}$ using
$\bm{X}^{\dagger}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}$ and note that
$\bm{X}^{\top}\bm{X}$ is symmetric (see (6)). Now $\mathcal{P}_{1}$ is
equivalent to

$\operatorname*{\textrm{argmin}}_{\begin{subarray}{c}k\in\mathbb{Z}^{+}\\\
\bm{X}\in\\{0,1\\}^{n\times k}\\\
\bm{X}\bm{1}_{k}=\bm{1}_{n}\end{subarray}}\frac{1}{2}\|\bm{Y}\|_{F}^{2}-\frac{1}{2}\operatorname{\textrm{tr}}\Big{(}\bm{YY}^{\top}\bm{X}\bm{X}^{\dagger}\Big{)}+\frac{\lambda}{2}\operatorname{\textrm{tr}}\big{(}\bm{X}^{\top}\bm{LX}\big{)},$

Multiply the whole expression by 2, ignore the constant, flip the sign, move
the negative sign to $\bm{L}$ and let
$\Phi=\Big{\\{}\bm{X}~{}|~{}\bm{X}\in\\{0,1\\}^{n\times
k},\bm{X}\bm{1}_{k}=\bm{1}_{n}\Big{\\}}$ gives $\mathcal{P}_{2}$.

####  3.1.2 ${\mathcal{P}_{2}^{k}}$ is approximately a vector partition
problem

Let $\mathcal{P}_{2}$ with a fixed $k$ be $\mathcal{P}_{2}^{k}$. By inspecting
the cost function $q$ of $\mathcal{P}_{2}^{k}$, we see that $\bm{YY}^{T}$ and
$\bm{L}$ are positive semi-definite, which motivates us to use spectral method
to approximate the solution of $\mathcal{P}_{2}^{k}$, based on the following
arguments:

  * •

in practice, the data matrix $\bm{Y}$ is often low-rank or approximately low
rank [28]; and

  * •

in practice, there are only a few clusters.

Therefore, we can approximate $q$ by only keeping the terms that are
associated with the top $k$ largest eigenvalues of $\bm{YY}^{\top}$ (and the
top $k$ smallest eigenvalues of $\bm{L}$). To proceed, below we define two
sets of spectral vectors.

##### Two sets of spectral vectors

Let the eigendecomposition of $\bm{YY}^{\top}$ and $\bm{L}$ be
$\bm{U}\bm{\Sigma}\bm{U}^{\top}$ and $\bm{V}\bm{\Gamma}\bm{V}^{\top}$,
respectively. That is, $(\bm{u}_{i},\sigma_{i})$ denotes the $i$th largest
(eigenvector, eigenvalue) of $\bm{YY}^{\top}$, and similarly for
$(\bm{v}_{i},\gamma_{i})$. Now, given a constant $\alpha\in\mathbb{R}$, define
two sets of $k$-dimensional vectors $\\{\bm{r}_{i}\\}_{i=1}^{n}$ and
$\\{\bm{t}_{i}\\}_{i=1}^{n}$ as

(8) $\displaystyle
r_{i}(j)=\sqrt{\sigma_{j}}U_{ij},~{}~{}t_{i}(j)=\sqrt{\alpha-\gamma_{j}}V_{ij},~{}~{}i\in[n]$

where $r_{i}(j)$ denotes the $j$th element of $\bm{r}_{i}$. Note that we will
select $\alpha$ such that $\alpha>\gamma_{j}$ so that $\bm{t}_{i}$ is a real
vector. We will discuss how to pick $\alpha$ in Section 3.1.3.

Theorem 3.3 shows that maximizing $q$ in $\mathcal{P}_{2}$ is (approximately)
equivalent to jointly assigning vectors $\\{\bm{r}_{i}\\}_{i=1}^{n}$ and
$\\{\bm{t}_{i}\\}_{i=1}^{n}$ into groups so as to maximizing the magnitude of
their weighted sum. Such an assignment problem is known as a vector partition
problem [18, 2] and it can be effectively solved by any k-means package. We
remark that the derivation technique of Theorem 3.3 is based on the work of
Mark Newman [18].

######  Theorem 3.3 (${\mathcal{P}_{2}^{k}}\approx$ Vector partition).

Given an integer $1\leq k\leq n$, a constant $\alpha>\gamma_{k}$, let
$|C_{i}|$ be the size of the partition $C_{i}$, then $\mathcal{P}_{2}$ is
approximately the following vector partition problem

(VPP)
$\displaystyle\operatorname*{\textrm{argmax}}_{\bm{r}_{i},\bm{t}_{i}}~{}\displaystyle\sum_{i=1}^{k}\|\bm{\xi}_{i}\|^{2}_{2}+\|\bm{\zeta}_{i}\|^{2}_{2}-\lambda\alpha
n~{}~{}~{}~{}\textrm{s.t.}~{}~{}~{}~{}\displaystyle\bm{\xi}_{i}=\sum_{i\in
C_{i}}\frac{1}{\sqrt{|C_{i}|}}\bm{r}_{i},~{}~{}\displaystyle\bm{\zeta}_{i}=\sum_{i\in
C_{i}}\sqrt{\lambda}\bm{t}_{i}.$

######  Proof 3.4.

The proof has three parts: (1) rewriting $q$ in $\mathcal{P}_{2}$ using
eigendecomposition, (2) define two sets of spectral vectors and (3)
approximate the problem as a vector partition problem.

##### Rewriting $q$ using eigendecomposition

Let the eigendecomposition on $\bm{YY}^{\top}$ and $\bm{L}$ be
$\bm{U}\bm{\Sigma}\bm{U}^{\top}$ and $\bm{V}\bm{\Gamma}\bm{V}^{\top}$,
respectively. That is, $\bm{U}=[\bm{u}_{1},\bm{u}_{2},\dots]$ stores the
eigenvectors of $\bm{YY}^{\top}$ and $\bm{\Sigma}$ is a diagonal matrix
storing the eigenvalues $\sigma_{i}$ of $\bm{YY}^{\top}$, and likewise for
$\bm{V}$, $\bm{\Gamma}$ for $\bm{L}$. Now in $\mathcal{P}_{2}$ we have
$q=\footnotesize{\operatorname{\textrm{tr}}\big{(}\bm{U}\bm{\Sigma}\bm{U}^{\top}\bm{X}\bm{X}^{\dagger}\big{)}+\lambda\operatorname{\textrm{tr}}\big{(}\bm{X}^{\top}\bm{V}(-\bm{\Gamma})\bm{V}^{\top}\bm{X}\big{)}}$.
Consider $-\bm{\Gamma}=\alpha\bm{I}-\bm{\Gamma}+\alpha\bm{I}$ under a
predefined $\alpha$, we have

$\begin{array}[]{rll}q&=\operatorname{\textrm{tr}}\Big{(}\bm{U}\bm{\Sigma}\bm{U}^{\top}\bm{X}(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\Big{)}+\lambda\operatorname{\textrm{tr}}\Big{(}\bm{X}^{\top}\bm{V}(\alpha\bm{I}-\bm{\Gamma})\bm{V}^{\top}\bm{X}\Big{)}-\lambda\alpha
n\\\
&=\operatorname{\textrm{tr}}\Big{(}(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{U}\bm{\Sigma}\bm{U}^{\top}\bm{X}\Big{)}+\lambda\operatorname{\textrm{tr}}\Big{(}(\alpha\bm{I}-\bm{\Gamma})\bm{V}^{\top}\bm{X}\bm{X}^{\top}\bm{V}\Big{)}-\lambda\alpha
n.\end{array}$

Now, invoke (6), we have

$q=\operatorname{\textrm{tr}}\left(\begin{bmatrix}\frac{1}{|C_{1}|}\\\
&\ddots\\\
&&\frac{1}{|C_{k}|}\end{bmatrix}\bm{X}^{\top}\bm{U}\bm{\Sigma}\bm{U}^{\top}\bm{X}+\lambda\begin{bmatrix}\alpha-\gamma_{1}\\\
&\ddots\\\
&&\alpha-\gamma_{n}\end{bmatrix}\bm{V}^{\top}\bm{X}\bm{X}^{\top}\bm{V}\right)-\lambda\alpha
n.$

Since trace operator focuses on the diagonal, we only need to look at the
diagonal of both $\bm{X}^{\top}\bm{U}\bm{\Sigma}\bm{U}^{\top}\bm{X}$ and
$\bm{V}^{\top}\bm{X}\bm{X}^{\top}\bm{V}$. Let $[\bm{M}]_{ii}$ denotes taking
only the diagonal of $\bm{M}$ and setting the rest to zero, the trace term of
$q$ is equivalent to

$\operatorname{\textrm{tr}}\left(\begin{bmatrix}\frac{1}{|C_{1}|}\\\
&\ddots\\\
&&\frac{1}{|C_{k}|}\end{bmatrix}[\bm{X}^{\top}\bm{U}]_{ii}[\bm{\Sigma}]_{ii}[\bm{U}^{\top}\bm{X}]_{ii}+\lambda\begin{bmatrix}\alpha-\gamma_{1}\\\
&\ddots\\\
&&\alpha-\gamma_{n}\end{bmatrix}[\bm{V}^{\top}\bm{X}]_{ii}[\bm{X}^{\top}\bm{V}]_{ii}\right),$

where we note that this expression is the motivation behind the two sets of
spectral vector defined in (8).

Now let $h=1,\dots,k$ be the index for cluster $C_{h}$, let $j=1,\dots,n$ be
the index of $\sigma_{j}$ and $\gamma_{j}$. For the matrix product
$\bm{C}=\bm{A}^{\top}\bm{B}$, recall that $C_{jh}=A_{ij}B_{ih}$ and note that
$\bm{A}^{\top}\bm{B}$ and $\bm{B}^{\top}\bm{A}$ share the same diagonal, now
$q$ becomes

$q~{}=~{}\sum_{h=1}^{k}\frac{1}{|C_{h}|}\sum_{j=1}^{n}\sigma_{j}\left(\sum_{i=1}^{n}U_{ij}X_{ih}\right)^{2}+\lambda\sum_{h=1}^{k}\sum_{j=1}^{n}(\alpha-\gamma_{j})\left(\sum_{i=1}^{n}V_{ij}X_{ih}\right)^{2}-\lambda\alpha
n,$

which gives

(9)
$q~{}=~{}\sum_{h=1}^{k}\frac{1}{|C_{h}|}\sum_{j=1}^{n}\left(\sum_{i=1}^{n}\sqrt{\sigma_{j}}U_{ij}X_{ih}\right)^{2}+\lambda\sum_{h=1}^{k}\sum_{j=1}^{n}\left(\sum_{i=1}^{n}\sqrt{\alpha-\gamma_{j}}V_{ij}X_{ih}\right)^{2}-\lambda\alpha
n.$

##### Define two sets of spectral vectors

Under a fixed $k$, we now approximately maximize $q$ by focusing on the terms
in $q$ that are related to first $k$ largest eigenvalues of $\bm{YY}^{\top}$
and $-\bm{L}$. First we follows (8) and define two sets of $k$-dimensional
vectors $\\{\bm{r}_{i}\\}_{i=1}^{n}$ and $\\{\bm{t}_{i}\\}_{i=1}^{n}$. Note
that we pick $\alpha>\gamma_{k}$ (the $k$th largest eigenvalues of $\bm{L}$)
so $\bm{t}_{i}$ are real for the first $k$ terms.

##### Approximate the problem as a vector partition problem

Now we apply $\\{\bm{r}_{i}\\}_{i=1}^{n}$ and $\\{\bm{t}_{i}\\}_{i=1}^{n}$ on
(9): in the first term, we drop the terms related to the smallest $n-k$
eigenvalues of $\bm{YY}^{\top}$; in the second term, we drop the terms related
to the smallest $n-k$ of the factors $\alpha-\gamma_{i}$, then

$\displaystyle q~{}$
$\displaystyle\approx~{}\sum_{h=1}^{k}\frac{1}{|C_{h}|}\sum_{j=1}^{k}\Big{(}\sum_{i=1}^{n}\sqrt{\sigma_{j}}U_{ij}X_{ih}\Big{)}^{2}+\lambda\sum_{h=1}^{k}\sum_{j=1}^{k}\Big{(}\sum_{i=1}^{n}\sqrt{\alpha-\gamma_{j}}V_{ij}X_{ih}\Big{)}^{2}-\lambda\alpha
n$
$\displaystyle=~{}\sum_{h=1}^{k}\frac{1}{|C_{h}|}\sum_{j=1}^{k}\left(\sum_{i\in
C_{h}}r_{ij}\right)^{2}+\lambda\sum_{h=1}^{k}\sum_{j=1}^{k}\left(\sum_{i\in
C_{h}}t_{ij}\right)^{2}-\lambda\alpha n$
$\displaystyle=~{}\sum_{h=1}^{k}\underbrace{\sum_{j=1}^{k}\left(\frac{1}{\sqrt{|C_{h}|}}\sum_{i\in
C_{h}}r_{ij}\right)^{2}}_{\|\bm{\xi}_{h}\|_{2}^{2}}+\sum_{h=1}^{k}\underbrace{\sum_{j=1}^{k}\left(\sqrt{\lambda}\sum_{i\in
C_{h}}t_{ij}\right)^{2}}_{\|\bm{\zeta}_{h}\|_{2}^{2}}-\lambda\alpha n$ (10)
$\displaystyle=~{}\sum_{h=1}^{k}\left\|\bm{\xi}_{h}\right\|^{2}+\left\|\bm{\zeta}_{h}\right\|^{2}-\lambda\alpha
n.$

Changing the index notation $h$ to $i$ in (10) completes the proof.

We follow [18] to solve VPP by applying k-means on
$\bm{z}_{i}\coloneqq\big{[}\frac{1}{\sqrt{|C_{i}|}}\bm{r}_{i}^{\top},\sqrt{\lambda}\bm{t}_{i}^{\top}\big{]}^{\top}$
for $i\in[n]$. This gives an approximate solution of $\mathcal{P}_{2}^{k}$.
See Algorithm 1.

Algorithm 1 Spectral method for solving $\mathcal{P}_{2}^{k}$

Input: $\bm{Y}\in\mathbb{R}^{n\times d}$, $\bm{L}\in\mathbb{R}^{n\times n}$,
$\lambda\in\mathbb{R}$, and $k$.

Compute:

Top $k$ largest eigenvalues $\\{\sigma_{i}\\}_{i=1}^{k}$ and the corresponding
eigenvectors $\\{\bm{u}_{i}\\}_{i=1}^{k}$ of $\bm{YY}^{\top}$;

Top $k$ smallest eigenvalues $\\{\gamma_{i}\\}_{i=1}^{k}$ and the
corresponding eigenvectors $\\{\bm{v}_{i}\\}_{i=1}^{k}$ of $\bm{L}$;

$\displaystyle\alpha=\frac{1}{n-k}\sum_{i=k+1}^{n}\gamma_{i}$;

$\bm{z}_{i}\coloneqq\big{[}\frac{1}{\sqrt{|C_{i}|}}\bm{r}_{i}^{\top},\sqrt{\lambda}\bm{t}_{i}^{\top}\big{]}^{\top}$
by (8), $i\in[n]$;

Run k-means on $\\{\bm{z}_{i}\\}_{i=1}^{n}$ and obtain $\bm{X}.$

Return $\bm{X}_{k}=\bm{X}$.

####  3.1.3 Approximation error and the choice of $\alpha$

In the approximation we drop the $n-k$ most negative eigenvalues of $-\bm{L}$
to approximate the matrix $\alpha\bm{I}-\bm{L}$ (i.e., the $n-k$ leading
eigenvalues of $\bm{L}$), and, instead of computing
$\bm{V}(\alpha\bm{I}-\bm{\Gamma})\bm{V}^{\top}$, we compute
$\bm{V}(\alpha\bm{I}^{\prime}-\bm{\Gamma}^{\prime})\bm{V}^{\top}$, where
$\bm{I}^{\prime}$ and $\bm{\Gamma}^{\prime}$ are matrix $\bm{I}$ and
$\bm{\Gamma}$ with the last $n-k$ diagonal elements set to zero. The
approximation error in such a procedure is

$\varepsilon\coloneqq\operatorname{\textrm{tr}}\Big{(}\bm{V}(\alpha\bm{I}-\bm{\Gamma})V^{\top}-\bm{V}(\alpha\bm{I}^{\prime}-\bm{\Gamma}^{\prime})\bm{V}^{\top}\Big{)}=\operatorname{\textrm{tr}}\Big{(}(\alpha\bm{I}-\bm{\Gamma})-(\alpha\bm{I}^{\prime}-\bm{\Gamma}^{\prime})\Big{)}=\sum_{i=k+1}^{n}(\gamma_{i}-\alpha)^{2}.$

Let $\dfrac{\partial\varepsilon}{\partial\alpha}=0$ gives the optimal
$\alpha^{*}=\dfrac{1}{n-k}\displaystyle\sum_{i=k+1}^{n}\gamma_{i}$.

####  3.1.4 The whole approximate spectral method

We now propose Algorithm 2 to solve $\mathcal{P}_{2}$ approximately.

Algorithm 2 Algorithm for solving $\mathcal{P}_{2}$ approximately

Input: $\bm{Y}\in\mathbb{R}^{n\times d}$, $\bm{L}\in\mathbb{R}^{n\times n}$,
$\lambda\in\mathbb{R}$, and $k_{\max}\in\mathbb{Z}^{+}$.

Initialize: $q^{*}=+\infty$.

for all $k=1,2,...,k_{\max}$ do

Using Algorithm 1 to get $\bm{X}_{k}$;

Compute
$q_{k}=\operatorname{\textrm{tr}}(\bm{YY}^{\top}\bm{X}_{k}\bm{X}_{k}^{\dagger})+\lambda\operatorname{\textrm{tr}}\Big{(}\bm{X}_{k}^{\top}(-\bm{L})\bm{X}_{k}\Big{)}$;

Update $q^{*}=q_{k}$ and $\bm{X}_{*}=\bm{X}_{k}$ if $q_{k}\leq q^{*}$.

end for

Return $\bm{X}_{*}$ and
$\bm{\mu}=(\bm{X}_{*}^{\top}\bm{X}_{*})^{-1}\bm{X}_{*}^{\top}\bm{Y}$.

###  3.2 Heat-bath simulated annealing

We remark that Algorithm 2 is based on approximating the eigen-spectrum of
$\bm{Y}\bm{Y}^{\top}$ and $\bm{L}$. Such approach is not suitable for
applications that require high accuracy. Now we propose another approach that
is able (in probability) to find the global optimal solution of
${\mathcal{P}_{1}}$ without approximation. First, we rewrite the mixed-integer
program $\mathcal{P}_{1}$ as

(11)
$\displaystyle\mathcal{P}_{3}~{}:~{}\operatorname*{\textrm{argmin}}_{k,\bm{\delta}}$
$\displaystyle~{}~{}\mathcal{H}(\bm{\delta})\coloneqq\underbrace{\sum_{i=1}^{n}\sum_{j=1}^{k}i_{=j}(\delta_{i})\left\|\bm{y}_{i}-\bm{c}_{j}\right\|_{2}^{2}}_{\textnormal{k-means}}+\lambda\underbrace{\sum_{i,j}\sum_{s\neq
t}A_{ij}i_{=s}(\delta_{i})i_{=t}(\delta_{j})}_{\textnormal{graph $k$-cut}}$
$\displaystyle s.t.$
$\displaystyle~{}~{}\bm{c}_{j}=\dfrac{\displaystyle\sum_{l=1}^{n}i_{=j}(\delta_{l})\bm{y}_{l}}{\displaystyle\sum_{l=1}^{n}i_{=j}(\delta_{l})},\
j\in\\{1,2,...,k\\};$ $\displaystyle~{}~{}\delta_{i}\in\\{1,2,...,k\\},\
i=1,2,...,n,$

where $\mathcal{H}$ denotes Hamiltonian energy,
$\bm{\delta}=[\delta_{1},\delta_{2},...,\delta_{n}]^{\top}\in\\{1,\dots,k\\}^{n}$
is a membership state vector, with the $i$th element
$\delta_{i}\in[k]\coloneqq\\{1,2,,\dots,k\\}$ representing the membership of
signal $\bm{y}_{i}$ and the membership of node $i$, e.g., $\delta_{i}=h$ means
the $\bm{y}_{i}$ and node $i$ belong to cluster $h$. The symbols $s,t\in[k]$
are two dummy cluster labels and the indicator function $i_{=j}(\delta_{i})=1$
if $\delta_{i}=j$, and $=0$ otherwise. The matrix $\bm{A}$ with elements
$A_{ij}\in\\{0,1\\}$ is the adjacency matrix of the graph $G$, and lastly,
$\bm{c}_{j}$ is an intermediate variable in $\mathcal{P}_{3}$.

It is trivial to see the equivalence between ${\mathcal{P}_{1}}$ and
$\mathcal{P}_{3}$: in $\mathcal{P}_{3}$ we uses $\bm{\delta}$ to replace the
assignment matrix $\bm{X}$ in $\mathcal{P}_{1}$. The membership states
$\\{\delta_{i}\\}$ act as a 0-1 on/off switches in both the k-means term and
the graph cut term.

####  3.2.1 Contrasting $\mathcal{P}_{1}$ and $\mathcal{P}_{3}$

Recall that in Algorithm 2 for (approximately) solving ${\mathcal{P}_{1}}$, we
screen through different values of $k$ in a pre-defined search set
$[k_{\max}]$ to find the optimal $k_{*}$ and we assume $k_{*}\in[k_{\max}]$.
For $\mathcal{P}_{3}$, we propose to use heat-bath Simulated Annealing (SA) to
find $k_{*}$. We start with an overestimated value $k$ (which is expected to
be larger than $k_{*}$), then SA will find $k_{*}$ (probabilistically) by a
“cooling process.”

####  3.2.2 Heat-bath simulated annealing

To minimize $\mathcal{H}(\bm{\delta})$, we use a heat-bath Monte Carlo
optimization strategy (abbreviated as heat-bath), which is a strategy similar
to the classical SA (abbreviated as cSA) but is shown to have a better
performance [4, 7]. While cSA and heat-bath are conceptually similar, cSA is
based on the Metropolis approach [15] to compute the transition probabilities
without computing partition functions because the partition functions are hard
to compute for many applications. However, in the setting of GTF presented in
this paper, efficient computation of the partition function is not a
challenge, which allows the use of the alternative heat-bath transition
probabilities.

##### Coordinate-wise update

In the proposed SA approach to solve $\mathcal{P}_{3}$, we perform a
coordinate update of $\bm{\delta}$: we update one element $\delta_{i}$ at a
time while keeping the rest of $\bm{\delta}$ fixed at their most recent value.

##### Details of the heat-bath update

As $\delta_{i}\in[k]$ represents a cluster label, the update on $\delta_{i}$
is a switch process. Let $s,t\in[k]$ be cluster labels, where $s$ is the state
of $\delta_{i}$ before the update and $t$ is the state of $\delta_{i}$ after
the update. Such switching update of $\delta_{i}$ is carried out using single-
spin heat-bath update rule [1]: we update the system energy $\mathcal{H}$ when
making a membership assignment $\delta_{i}$ switch from $s$ to $t$, i.e.,

$\underbrace{\mathcal{H}(\delta_{i}=t)}_{\text{energy after
update}}=\underbrace{\mathcal{H}(\delta_{i}=s)}_{\text{energy before
update}}+\underbrace{\Delta\mathcal{H}(\delta_{i}=s\rightarrow
t)}_{\text{change in energy}}.$

The probability of the membership switch at a temperature $T$ is proportional
to the exponential of the corresponding energy change of the entire system,
i.e

(12)
$\displaystyle\mathbb{P}(\delta_{i}=t;T)=\dfrac{\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
t)\right\\}}{\displaystyle\sum_{l=1}^{k}\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
l)\right\\}}.$

##### A brief summary of how SA finds global minimum

From statistical mechanics, at a high enough temperature, the membership state
$\delta_{i}$ follows a discrete uniform distribution, i.e.,

$\color[rgb]{0,0,0}\lim_{T\rightarrow+\infty}\left(\underbrace{\dfrac{\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
t)\right\\}}{\displaystyle\sum_{l=1}^{k}\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
l)\right\\}}}_{\mathbb{P}(\delta_{i}=t;T)}\right)~{}=~{}\dfrac{\displaystyle\lim_{T\rightarrow+\infty}\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
t)\right\\}}{\displaystyle\sum_{l=1}^{k}\lim_{T\rightarrow+\infty}\exp\left\\{-\frac{1}{T}\Delta\mathcal{H}(\delta_{i}=s\rightarrow
l)\right\\}}~{}=~{}\frac{1}{k},$

hence at such moment, it is capable of jumping out of a local minimum, or in
other words, SA at high temperatures is able to explore all the $k$ possible
membership assignments. When the temperature cools down (say $T\rightarrow
1$), the state $\delta_{i}$ obeys a categorical distribution where
$\delta_{i}$ are more likely to switch to $t$ from the current membership
state $s$ if the energy change $\Delta\mathcal{H}(\delta_{i}=s\rightarrow t)$
is negative and small.

##### The whole SA algorithm

Based on the above discussion, we propose the SA Algorithm 3, where scalars
$T_{start}$ and $T_{end}$ are the start and end temperatures, $c\in[0,1)$ is a
constant to control the cooling rate, and $K$ is the number of sweep times for
each temperature.

Algorithm 3 Hot-bath SA for solving $\mathcal{P}_{3}$.

Input: $T_{start}=100$, $T_{end}=0.001$, $c=0.99$, $K$, $\lambda$.

$T=T_{start}$.

while $T\geq T_{end}$ do

for all $k=1,2,...,K$ do

$R=$ RandomShuffle($\\{1,2,\dots,n\\}$);

for all $i=1,2,...,n$ do

Update $\mathcal{H}(\delta_{R[i]})$ at $T$ based on the probability (12) using
roulette wheel algorithm [6];

Update $\mathcal{H}$ after $\delta_{R[i]}$ is assigned;

end for

end for

$T=T\times c$ % cool down the temperature.

end whileReturn The most recent $\bm{\delta}$.

The major computational cost here comes from computing the term
$\Delta\mathcal{H}(\delta_{i}=s\rightarrow t)$, which can be efficiently
computed in $\mathcal{O}(n)$ [23, 16]. However, we remark that it can take
many iterations for the Hot-bath SA algorithm to converge (coold down).

###  3.3 Section summary

We recall that the main contribution of this work is the new GTF model
$\mathcal{P}_{0}$. In this section we proposed two algorithmic approaches for
solving $\mathcal{P}_{0}$. In general, compared with the spectral method, the
Hot-bath SA algorithm provides a more accurate solution but it also requires
more computational time to find the solution. We refer to the discussion in
Section 5 for more details.

##  4 Extension to graph-based Transductive Learning

In this section, we show that the Theorem 2.1 is generalizable to other
models. Let $K>0$ be a given integer. Consider a $K$-class classification
problem in a semi-supervised learning setting: given a dataset with $n$
samples, we observe a subset of class labels encoded in a class assignment
matrix $\bm{Y}\in\\{0,1\\}^{n\times K}$, and a diagonal indicator matrix
$\bm{M}\in\\{0,1\\}^{n\times n}$ that denotes the samples whose labels have
been observed. Such classification problem, called the Modified Absorption
Problem [26, 27] can be expressed as a 1st-order GTF with $\ell_{2,0}$ norm
penalty as

(13)
$\displaystyle\mathcal{Q}~{}:~{}\operatorname*{\textrm{argmin}}_{\bm{B}\in\mathbb{R}^{n\times
K}}\frac{1}{2}\left\|\bm{M}(\bm{Y}-\bm{B})\right\|_{F}^{2}+\lambda\mathcal{L}(\bm{B})+\epsilon\left\|\bm{R}-\bm{B}\right\|_{F}^{2},$

where $\bm{R}\in\mathbb{R}^{n\times K}$ encodes a fixed prior belief on the
samples, and $\lambda\geq 0,\epsilon\geq 0$ are regularization parameters. The
matrix $\bm{B}\in\mathbb{R}^{n\times K}$ contains the fitted probabilities,
with $B_{ij}$ representing the probability that sample $i$ is of class $j$.
Using $\bm{b}_{i}$ to denote the $i$th row of $\bm{B}$, the middle term of
(13) encourages the probabilities of samples to behave smoothly over a graph
$G(V,E)$ constructed based on the sample features.

Now we give a theorem related to $\mathcal{Q}$. The focus here is to
demonstrate the generality of Theorem 2.1, which can be applied to models with
the 1st-order GTF with $\ell_{2,0}$ norm penalty.

######  Theorem 4.1.

$\mathcal{Q}$ is equivalent to

(14)
$\mathcal{Q}_{1}:{\operatorname*{\textrm{argmin}}_{\begin{subarray}{c}k\in\mathbb{Z}^{+}\\\
\bm{X}\in\\{0,1\\}^{n\times
k},\bm{X}^{\top}\bm{1}_{k}=\bm{1}_{n}\end{subarray}}}\begin{array}[]{l}\frac{1}{2}\Big{\|}\bm{M}\Big{(}\bm{Y}-\bm{X}\bm{B}(\bm{X})\Big{)}\Big{\|}_{F}^{2}+\lambda\operatorname{\textrm{tr}}(\bm{X}^{\top}\bm{LX})+\epsilon\Big{\|}\bm{R}-\bm{X}\bm{B}(\bm{X})\Big{\|}_{F}^{2},\end{array}$

with
$\bm{B}\big{(}\bm{X})=(\bm{X}^{\top}(\bm{M}+\bm{I})\bm{X}\big{)}^{-1}\big{(}\bm{X}^{\top}\bm{MY}+\epsilon\bm{X}^{\top}\bm{R}\big{)}$.

We use the techniques for solving $\mathcal{P}_{2}$ to solve $\mathcal{Q}_{1}$
(hence solving $\mathcal{Q}$). Note that $\mathcal{Q}$ and $\mathcal{P}$ have
a different optimization subproblem, so we cannot use Algorithm 1 directly to
solve the inner minimization of $\mathcal{Q}$. We use the Frank-Wolfe
algorithm to solve such inner minimization of $\mathcal{Q}$.

Figure 2: On identifying boundary edges. Left: the ground truth $\bm{y}_{*}$
on the Minnesota road graph. Middle: the ROC curve of the methods. The curves
for the spectral method and SA method are overlapped. Here the $\ell_{0}$
model achieved perfect ROC. Right: the optimal $k_{*}$ ($=4$ clusters) for the
$\ell_{0}$ model $\mathcal{P}_{2}$.

##  5 Experimental results

We now demonstrate the proposed $\ell_{2,0}$ norm penalized 1st-order GTF
model (i.e., problem (3)) outperforms three state-of-the-art GTF models (see
table 2) over two graph learning tasks: support recovery and signal
estimation. Then we show the proposed MAP-GTF model outperforms other models
in the task of semi-supervised classification. All experiments were run on a
computer with 8 cores 3.7GHz Intel CPU, and 32 GB RAM111The code is available
at: https://github.com/EJIUB/GTF_L04.

Table 2: GTF models with different penalty Penalty | Reference | parameters
---|---|---
$\ell_{2,0}$ | this work | $k,\lambda$
$\ell_{1}$ | [31] | $\lambda$
SCAD | [29] | $\gamma,\tau,\lambda$
MCP | [29] | $\gamma,\tau,\lambda$

###  5.1 Support recovery

####  5.1.1 Dataset

We test all the methods for identifying the boundaries of a piece-wise
constant signal on the Minnesota road graph [22] shown in Fig. 2. The graph
$G$ has $n=2,642$ nodes and $3,304$ edges. Let $\bm{y}_{*}\in\mathbb{R}^{n}$
be the ground truth scalar-valued signal for each node in $G$, we add Gaussian
noise $\bm{\theta}\sim\mathcal{N}(0,\sigma^{2}I)$ with $\sigma^{2}=0.05$ to
$\bm{y}_{*}$ to obtain the simulated observed signal
$\bm{y}=\bm{y}_{*}+\bm{\theta}$. The Signal-to-Noise Ratio (SNR) of
$\bm{y}_{*}$ is $10\log_{10}\frac{\|\bm{y}_{*}\|}{\sigma^{2}nd}\approx 12$.

####  5.1.2 Support recovery: perfect ROC curve for the proposed method

The task here is to identify the boundary edges by recovering the true
supports of $\left\|\bm{b}_{i}-\bm{b}_{j}\right\|_{0},\ \forall(i,j)\in E$. We
evaluate the performance of the aforementioned methods by comparing their
result with the ground truth boundary edges computed from $\bm{y}_{*}$ as
shown in the left panel of Fig. 2.

####  5.1.3 Parameter tuning

In the test, we report the best results for each model by tuning their hyper-
parameters.

  * •

We run Algorithm 2 to use the spectral method to solve the proposed
$\ell_{2,0}$ model (3) for the support recovery problem. The algorithm has two
hyper-parameters: the number of partitions $k$ and the regularization
parameter $\lambda$.

  * •

We run Algorithm 3 to use the SA method to solve the proposed $\ell_{2,0}$
model (11). We set the number of partitions $k=7$. The algorithm has only one
hyper-parameter: the regularization parameter $\lambda$.

  * •

For the $\ell_{1}$, SCAD, and MCP models, we use the source code provided by
[29] to identify supports for this experiment. For SCAD, we set $\gamma=3.7$
as suggested in [29]. For MCP, we set $\gamma=1.4$ as suggested in [29].

Figure 3:  Left: the graph structure of $G_{1}$ (top) and $G_{2}$ (bottom) and
the corresponding ground truth signal input to each node in both graphs.
Right: the corresponding plots of input signal SNR versus reconstructed signal
SNR for the GTF models. Note that in both graphs, the $\ell_{2,0}$ model
consistently outperforms other GFT models.

####  5.1.4 The results

We plot the results in Fig. 2. We compare the methods in terms of ROC curve,
i.e. the true positive rate versus the false positive rate of identification
of a boundary edge correctly. The $\ell_{2,0}$ model solved by both the
proposed spectral method and the SA method obtains a perfect ROC curve, which
significantly outperforms other models. The GTF models with nonconvex
penalties (SCAD and MCP) have similar performance, which is better than the
$\ell_{1}$ model.

It is not surprising that the proposed $\ell_{2,0}$ model achieves a
significantly better result; see the discussion in Remark 2.3.

In Fig. 2 we also show that we can find the $k_{*}$ for $\mathcal{P}_{2}$ by
screening different values of $k$ in Algorithm 2. As the search space of $k$
is small, hence the proposed method runs very fast. In addition, the SA method
can also find the $k_{*}$ as long as we set $k$ larger than $k_{*}$.

###  5.2 Denoising using GTF

We now compare the performance of the GTF models on denoising.

####  5.2.1 Dataset: graph generation

For the dataset, we generate a simple connected undirected unweighted graph
$G(V,E)$ with $200$ nodes. We plant three modules/communities in $G$, which
have 50 nodes, 70 nodes, and 80 nodes, respectively. We introduce two
parameters $(p,q)$: $p$ controls the probability of connections between nodes
within the modules/communities and $q$ controls the probability of connections
between nodes across the modules/communities. Then we generate two graphs
$G_{1},G_{2}$ as in Fig.3. For $G_{1}$, we set $p=0.2$ and $q=0.05$ and for
$G_{2}$, we set $p=0.05$ and $q=0.01$. Clearly, the modular/community
structure of $G_{2}$ is relatively less clear than that of $G_{1}$.

####  5.2.2 Dataset: signal generation

We generate a ground truth signal $\bm{y}_{*}\in\mathbb{R}^{200}$ on the
graphs $G_{1}$ and $G_{2}$ where on the module/community with $50$ nodes has
an input signal value of $1$, on the module/community with $70$ nodes has an
input signal value of $-1$, and on the module/community with $80$ nodes has
input signal value of $0$. We simulate $d=10$ independent noisy realizations
of the graph signal and concatenate them to construct a noisy vector-valued
graph signal
$\bm{Y}=[\bm{y}_{*}+\bm{\theta}_{1},\bm{y}_{*}+\bm{\theta}_{2},...]\in\mathbb{R}^{200\times
10}$, where all $\bm{\theta}_{i}\in\mathbb{R}^{200}$ are independently drew
from $\theta_{i}\sim\mathcal{N}(0,\sigma^{2}I)$.

####  5.2.3 The learning task

The goal here is to use different GTF models to recover the vector-valued
signal under different input signal SNR controlled by $\sigma$. We recall that
the reconstructed signal SNR is computed as
$10\log_{10}(\|\bm{Y}^{*}\|_{F}/\|\bm{B}-\bm{Y}^{*}\|_{F})$, where $\bm{B}$ is
the reconstructed signal.

####  5.2.4 Parameter tuning

We use the same setup of hyper-parameters tuning as in the experiment on
support recovery, and for each GTF model, we screen through the hyper-
parameters and report the best result.

####  5.2.5 The results

In Fig.3, we plot the SNR of the reconstructed signal on the graphs $G_{1}$
(top) and $G_{2}$ (bottom) across different SNR values of the input signal for
all the GTF models. The curve of the $\ell_{2,0}$ model solved by the proposed
method is above all the other GTF models, which demonstrates that our model
can take advantage of the graph structure to recover the true signal. For the
reconstruction of signal on the graph $G_{2}$, which is a harder problem as
the modular/community structure of $G_{2}$ is less clear than $G_{1}$, the
$\ell_{2,0}$ model still has a better reconstruction SNR than all other GTF
models, with a consistent amount of at least 5dB better SNR in the high input
SNR regime. Furthermore, we find that the SA method has better performance
than the spectral method when the noisy signal SNR is smaller for signal
recovery on $G_{2}$.

Table 3: Time comparison (mean(std)) with the number of edges in graphs increasing. |  ($p=0.1,q=0.01$) |  ($p=0.5,q=0.1$) |  ($p=0.9,q=0.2$)
---|---|---|---
| $|E|$ | Time(s) | $|E|$ | Time(s) | $|E|$ | Time(s)
$\ell_{1}$ | 805(23) | 4.71(1.22) | 4707(69) | 56.74(20.43) | 8750(40) | 134.92(46.53)
SCAD | 805(23) | 5.22(1.24) | 4707(69) | 50.32(11.72) | 8750(40) | 106.03(25.34)
MCP | 805(23) | 7.71(2.26) | 4707(69) | 42.71(14.01) | 8750(40) | 105.64(33.92)
Alg. 2 | 805(23) | 0.01(0.01) | 4707(69) | 0.01(0.01) | 8750(40) | 0.01(0.01)
Alg. 3 | 805(23) | 990.54(5.42) | 4707(69) | 993.14(6.41) | 8750(40) | 992.94(8.91)

###  5.3 Computational time

In this section, we compare the computational time of the competing GTF
models.

####  5.3.1 Dataset

We follow the procedures described in sections 5.2.1 and 5.2.2 to generate
graph signals for time comparison. The number of nodes and the community
structure of the graphs are the same as described in section 5.2.1. All
generated graphs have 200 nodes and 3 communities with 50, 70, and 80 nodes,
respectively. We control the number of edges of the graphs by setting
different $(p,q)$ pairs (specified in Table 3). We follow the procedure in
section 5.2.2 to generate the signal on graphs.

####  5.3.2 Parameter

For each GTF model, we use the best parameters obtained in section 5.2.4.

####  5.3.3 The results

In Tabel 3, we show the time comparison for different GTF models on denoising
graph signals. For each $(p,q)$ pair, we generate 10 graph signals with 10
different graphs, which have the same community structure but different
numbers of edges. We show the average and standard deviation of computational
time for each GTF model. We observe that the computational times of previous
GTF models ($\ell_{1}$, SCAD, and MCP) increase with the number of edges in
the graphs. However, for the proposed GTF model, the computational times of
the spectral and SA algorithms are not sensitive to the edge size of the
graph.

Table 4: Misclassification rates averaged over 100 trials, with the variance below the rates. | heart | wine | iris | breast | car | wine-qual. | ads | yeast
---|---|---|---|---|---|---|---|---
# of classes | 2 | 3 | 3 | 2 | 4 | 5 | 2 | 10
# of samples | 303 | 178 | 150 | 569 | 1,728 | 1,599 | 3,279 | 1,484
$\ell_{1}$ |  | 0.158
---
(4e-4)
| 0.040
---
(2e-4)
| 0.051
---
(5e-4)
| 0.058
---
(4e-4)
| 0.170
---
(6e-5)
| 0.574
---
(7e-4)
| 0.152
---
(1e-4)
| 0.391
---
(4e-4)
SCAD |  | 0.150
---
(4e-2)
| 0.041
---
(5e-4)
| 0.049
---
(2e-4)
| 0.045
---
(2e-4)
| 0.149
---
(3e-5)
| 0.574
---
(1e-5)
| 0.153
---
(4e-4)
| 0.673
---
(3e-5)
MCP |  | 0.148
---
(3e-4)
| 0.041
---
(1e-4)
| 0.048
---
(3e-4)
| 0.041
---
(2e-4)
| 0.148
---
(1e-5)
| 0.572
---
(1e-6)
| 0.150
---
(5e-4)
| 0.507
---
(1e-5)
$\ell_{2,0}$ (see Section 4) |  | 0.143
---
(2e-4)
| 0.035
---
(6e-5)
| 0.038
---
(3e-4)
| 0.033
---
(1e-4)
| 0.098
---
(1e-6)
| 0.373
---
(3e-4)
| 0.079
---
(5e-4)
| 0.372
---
(5e-5)

###  5.4 Semi-supervised classification

We now report results on the Modified Absorption Problem (MAP). We compare the
proposed MAP with the $\ell_{2,0}$ norm GTF regularization (13) with other
MAP-GTF models proposed in [31, 29], which include MAP with $\ell_{1}$ norm
GTF regularization [31], MAP with SCAD penalized GTF regularization [29], and
a MAP with MCP penalized GTF regularization [29].

####  5.4.1 Dataset

We test the models on 8 popular classification datasets from [12]. For each
dataset, we construct a $5$-nearest neighbor graph based on the Euclidean
distance between provided features. We randomly pick $20\%$ samples from the
dataset and treat them as labeled samples. For each dataset, we compute the
misclassification rates over 100 repetitions.

####  5.4.2 Parameter tuning

For all the MAP-GTF models, we set $\epsilon=0.01$ (the regularization
parameter in (13)) and $R_{ij}=1/K$. For each model, we screen through the
hyper-parameters and report the best results.

####  5.4.3 The results

Table 4 summarizes the misclassification rate for all the methods over 100
repetitions of randomly selected labeled samples. The proposed MAP-GTF model
(13) achieves the smallest misclassification rate. We further do a t-test
based on the misclassification rates over 100 repetitions and find that the
misclassification rates of the proposed MAP-GTF model (13) for all the
datasets are significantly smaller than other methods ($p$-value is less than
0.001).

####  5.4.4 Concluding remarks

Lastly, we comment on the comparisons of these GTF regularizers on MAP. The
experiment results all show that the $\ell_{1}$ model while being convex and
having the least number of parameters (see Table 2), is consistently having
the worst performance among all the tested models. For the remaining nonconvex
models, the $\ell_{2,0}$ model consistently outperforms SCAD and MCP, and has
fewer parameters to tune. We believe that the superior performance of the MAP
model with the $\ell_{2,0}$ norm GTF regularization in such a broad range of
experiments is convincing evidence for the utility of the $\ell_{2,0}$ norm
GTF regularizer.

##  6 Conclusion

We propose a first-order GTF model with $\ell_{2,0}$ norm penalty. By
converting it into a mixed integer program, we prove that the model is
equivalent to a k-means clustering with graph cut regularization (Theorem
2.3). We propose two methods to solve the problem: approximate spectral method
and simulated annealing. Empirically, we demonstrate the performance
improvements of the proposed GTF model on both synthetic and real datasets for
support recovery, signal estimation, and semi-supervised classification.

We believe that the superior performance of the proposed models in a broad
range of experiments is convincing evidence for the utility of $\ell_{2,0}$
norm penalized GTF. For future works, we think it will be interesting to
investigate fast algorithms for solving the GTF model and study the error
bound of the model. In addition, we believe that higher-order GTF models (see
Fig.3 in [31]) with $\ell_{0}$ norm penalty also have great potential.

##  Appendix A Code

The code is available at: https://github.com/EJIUB/GTF_L0.

## References

  * [1] M. P. Allen and D. J. Tildesley, Advanced Monte Carlo methods, in Computer Simulation of Liquids, Oxford University Press, 06 2017.
  * [2] C. Alpert, A. Kahng, and S.-Z. Yao, Spectral partitioning with multiple eigenvectors, Discrete Applied Mathematics, 90 (1999), pp. 3–26.
  * [3] M. Belkin, I. Matveeva, and P. Niyogi, Regularization and semi-supervised learning on large graphs, in Learning Theory, pp. 624–638.
  * [4] S. Bornholdt, Statistical mechanics of community detection, Physical Review E, 74 (2006), p. 016110.
  * [5] S. Bougleux, A. Elmoataz, and M. Melkemi, Discrete regularization on weighted graphs for image and mesh filtering, in Scale Space and Variational Methods in Computer Vision, pp. 128–139.
  * [6] T. Bäck, Evolutionary Algorithms in Theory and Practice: Evolution Strategies, Evolutionary Programming, Genetic Algorithms, Oxford University Press, 02 1996.
  * [7] E. Carlon, Computational Physics, tech. report, KU Leuven, 2013.
  * [8] S. Chen, A. Sandryhaila, J. M. F. Moura, and J. Kovačević, Signal recovery on graphs: Variation minimization, IEEE Transactions on Signal Processing, 63 (2015), pp. 4609–4624.
  * [9] S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, Discrete signal processing on graphs: Sampling theory, IEEE Transactions on Signal Processing, 63 (2015), pp. 6510–6523.
  * [10] S. Chen, R. Varma, A. Singh, and J. Kovacevic, Signal representations on graphs: Tools and applications, CoRR, abs/1512.05406 (2015).
  * [11] F. R. K. Chung, Spectral Graph Theory, American Mathematical Society, 1997.
  * [12] D. Dua and C. Graff, UCI Machine Learning Repository, 2017.
  * [13] A. Elmoataz, O. Lezoray, and S. Bougleux, Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing, IEEE Transactions on Image Processing, 17 (2008), pp. 1047–1060.
  * [14] N. Fan and P. Pardalos, Multi-way clustering and biclustering by the ratio cut and normalized cut in graphs, Journal of Combinatorial Optimization, 23 (2010), pp. 224–251.
  * [15] W. K. Hastings, Monte carlo sampling methods using markov chains and their applications, Biometrika, 57 (1970), pp. 97–109.
  * [16] D. S. Johnson, C. R. Aragon, L. A. McGeoch, and C. Schevon, Optimization by simulated annealing: An experimental evaluation part i, graph partitioning, Operations Research, 37 (1989), pp. 865–892.
  * [17] M. Mahajan, P. Nimbhorkar, and K. Varadarajan, The planar k-means problem is np-hard, in WALCOM: Algorithms and Computation, pp. 274–285.
  * [18] M. Newman, Finding community structure in networks using the eigenvectors of matrices, Physical Review E, 74 (2006).
  * [19] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, Graph signal processing: Overview, challenges, and applications, Proceedings of the IEEE, 106 (2018), pp. 808–828.
  * [20] G. A. Pavlopoulos, M. Secrier, C. N. Moschopoulos, T. G. Soldatos, S. Kossida, J. Aerts, R. Schneider, and P. G. Bagos, Using graph theory to analyze biological networks, BioData Min, 4 (2011), p. 10.
  * [21] D. Romero, M. Ma, and G. B. Giannakis, Kernel-based reconstruction of graph signals, IEEE Transactions on Signal Processing, 65 (2017), pp. 764–778.
  * [22] R. A. Rossi and N. K. Ahmed, The network data repository with interactive graph analytics and visualization, in Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, AAAI’15, p. 4292–4293.
  * [23] S. Z. Selim and K. Alsultan, A simulated annealing algorithm for the clustering problem, Pattern Recognition, 24 (1991), pp. 1003–1008.
  * [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, 30 (2013), pp. 83–98.
  * [25] D. A. Spielman, Spectral graph theory and its applications, in Proceedings of the 48th Annual IEEE Symposium on Foundations of Computer Science, FOCS ’07, p. 29–38.
  * [26] P. Talukdar and K. Crammer, New regularized algorithms for transductive learning, in Machine Learning and Knowledge Discovery in Databases, 2009, pp. 442–457.
  * [27] P. Talukdar and F. Pereira, Experiments in graph-based semi-supervised learning methods for class-instance acquisition, in Proceedings of the 48th Annual Meeting of the Association for Computational Linguistics, pp. 1473–1481.
  * [28] M. Udell and A. Townsend, Why are big data matrices approximately low rank?, SIAM Journal on Mathematics of Data Science, 1 (2019), pp. 144–160.
  * [29] R. Varma, H. Lee, J. Kovačević, and Y. Chi, Vector-valued graph trend filtering with non-convex penalties, IEEE Transactions on Signal and Information Processing over Networks, 6 (2020), pp. 48–62.
  * [30] V. Vazirani, Approximation algorithms, Springer, 2001.
  * [31] Y. Wang, J. Sharpnack, A. J. Smola, and R. J. Tibshirani, Trend filtering on graphs, Journal of Machine Learning Research, 17 (2016), pp. 1–41.
  * [32] X. Zhu, Z. Ghahramani, and J. Lafferty, Semi-supervised learning using gaussian fields and harmonic functions, in Proceedings of the Twentieth International Conference on International Conference on Machine Learning, ICML’03, p. 912–919.